2 Created on Tue Jul 12 09:28:23 2016
12 import carmaWrap
as ch
15 import matplotlib.pyplot
as plt
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)
47 if (len(sys.argv) > 2):
48 savename = sys.argv[2]
50 savename =
"roket_default.h5"
52 print(
"save file is ", savename)
61 if (hasattr(config,
"simul_name")):
62 if (config.simul_name
is None):
65 simul_name = config.simul_name
68 print(
"simul name is", simul_name)
71 if (simul_name == b
""):
75 param_dict = h5u.params_dictionary(config)
76 matricesToLoad = h5u.checkMatricesDataBase(os.environ[
"SHESHA_ROOT"] +
"/data/",
81 c = ch.carmaWrap_context(devices=np.array([0, 1], dtype=np.int32))
86 wfs, tel = ao.wfs_init(config.p_wfss, config.p_atmos, config.p_tel, config.p_geom,
87 config.p_target, config.p_loop, config.p_dms)
91 atm = ao.atmos_init(c, config.p_atmos, config.p_tel, config.p_geom, config.p_loop,
92 config.p_wfss, wfs, config.p_target, rank=0, clean=clean,
97 dms = ao.dm_init(config.p_dms, config.p_wfss, wfs, config.p_geom, config.p_tel)
101 tar = ao.target_init(c, tel, config.p_target, config.p_atmos, config.p_geom,
102 config.p_tel, config.p_dms)
106 rtc = ao.rtc_init(tel, wfs, config.p_wfss, dms, config.p_dms, config.p_geom,
107 config.p_rtc, config.p_atmos, atm, config.p_tel, config.p_loop,
108 clean=clean, simul_name=simul_name, load=matricesToLoad)
111 h5u.validDataBase(os.environ[
"SHESHA_ROOT"] +
"/data/", matricesToLoad)
113 print(
"====================")
115 print(
"====================")
116 print(
"objects initialzed on GPU:")
117 print(
"--------------------------------------------------------")
124 print(
"----------------------------------------------------")
125 print(
"iter# | SE SR image | LE SR image | Fitting | LE SR phase var")
126 print(
"----------------------------------------------------")
128 error_flag =
True in [w.roket
for w
in config.p_wfss]
140 Performs the main AO loop for n interations. First, initialize buffers
141 for error breakdown computations. Then, at the end of each iteration, just
142 before applying the new DM shape, calls the error_breakdown function.
144 :param n: (int) : number of iterations
147 com : (np.array((n,nactus))) : full command buffer
149 noise_com : (np.array((n,nactus))) : noise contribution for error breakdown
151 alias_wfs_com : (np.array((n,nactus))) : aliasing estimation in the WFS direction
153 tomo_com : (np.array((n,nactus))) : tomography error estimation
155 H_com : (np.array((n,nactus))) : Filtered modes contribution for error breakdown
157 trunc_com : (np.array((n,nactus))) : Truncature and sampling error of WFS
159 bp_com : (np.array((n,nactus))) : Bandwidth error estimation on target
161 mod_com : (np.array((n,nactus))) : Commanded modes expressed on the actuators
163 fit : (float) : fitting (mean variance of the residual target phase after projection)
165 SR : (float) : final strehl ratio returned by the simulation
169 nactu = rtc.get_command(0).size
170 nslopes = rtc.get_centroids(0).size
171 com = np.zeros((n, nactu), dtype=np.float32)
172 noise_com = np.zeros((n, nactu), dtype=np.float32)
173 alias_wfs_com = np.copy(noise_com)
174 wf_com = np.copy(noise_com)
175 tomo_com = np.copy(noise_com)
176 trunc_com = np.copy(noise_com)
177 H_com = np.copy(noise_com)
178 mod_com = np.copy(noise_com)
179 bp_com = np.copy(noise_com)
185 for i
in range(-10, n):
188 if (config.p_controllers[0].type == b
"geo"):
189 for t
in range(config.p_target.ntargets):
190 tar.atmos_trace(t, atm, tel)
191 rtc.docontrol_geo(0, dms, tar, 0)
192 rtc.applycontrol(0, dms)
195 for t
in range(config.p_target.ntargets):
196 tar.atmos_trace(t, atm, tel)
198 for w
in range(len(config.p_wfss)):
199 wfs.sensors_trace(w,
"all", tel, atm, dms)
200 wfs.sensors_compimg(w)
205 if (error_flag
and i > -1):
209 roket.computeBreakdown()
210 rtc.applycontrol(0, dms)
212 if ((i + 1) % 100 == 0
and i > -1):
213 strehltmp = tar.get_strehl(0)
214 print(i + 1,
"\t", strehltmp[0],
"\t", strehltmp[1],
"\t",
215 np.exp(-strehltmp[2]),
"\t", np.exp(-strehltmp[3]))
217 print(
" loop execution time:", t1 - t0,
" (", n,
"iterations), ", (t1 - t0) / n,
218 "(mean) ", n / (t1 - t0),
"Hz")
221 SR2 = np.exp(-tar.get_strehl(0, comp_strehl=
False)[3])
222 SR = tar.get_strehl(0, comp_strehl=
False)[1]
230 Performs the main AO loop for n interations. First, initialize buffers
231 for error breakdown computations. Then, at the end of each iteration, just
232 before applying the new DM shape, calls the error_breakdown function.
234 :param n: (int) : number of iterations
237 com : (np.array((n,nactus))) : full command buffer
239 noise_com : (np.array((n,nactus))) : noise contribution for error breakdown
241 alias_wfs_com : (np.array((n,nactus))) : aliasing estimation in the WFS direction
243 tomo_com : (np.array((n,nactus))) : tomography error estimation
245 H_com : (np.array((n,nactus))) : Filtered modes contribution for error breakdown
247 trunc_com : (np.array((n,nactus))) : Truncature and sampling error of WFS
249 bp_com : (np.array((n,nactus))) : Bandwidth error estimation on target
251 mod_com : (np.array((n,nactus))) : Commanded modes expressed on the actuators
253 fit : (float) : fitting (mean variance of the residual target phase after projection)
255 SR : (float) : final strehl ratio returned by the simulation
257 for i
in range(0, n):
260 if (config.p_controllers[0].type == b
"geo"):
261 for t
in range(config.p_target.ntargets):
262 tar.atmos_trace(t, atm, tel)
263 rtc.docontrol_geo(0, dms, tar, 0)
264 rtc.applycontrol(0, dms)
266 for t
in range(config.p_target.ntargets):
267 tar.atmos_trace(t, atm, tel)
268 for w
in range(len(config.p_wfss)):
269 wfs.sensors_trace(w,
"all", tel, atm, dms)
270 wfs.sensors_compimg(w)
274 rtc.applycontrol(0, dms)
284 IF = rtc.get_IFsparse(1).T
292 delta = IF.T.dot(IF).toarray() / N
295 Tp = np.ones((T.shape[0], T.shape[1] + 1))
297 deltaT = IF.T.dot(Tp) / N
299 tau = np.linalg.inv(delta).dot(deltaT)
303 tdt = tau.T.dot(delta).dot(tau)
304 subTT = tau.dot(np.linalg.inv(tdt)).dot(tau.T).dot(delta)
308 gdg = G.T.dot(delta).dot(G)
309 U, s, V = np.linalg.svd(gdg)
310 U = U[:, :U.shape[1] - 3]
312 L = np.identity(s.size) / np.sqrt(s)
317 Btt = np.zeros((n + 2, n - 1))
318 Btt[:B.shape[0], :B.shape[1]] = B
319 mini = 1. / np.sqrt(TT)
322 Btt[n:, n - 3:] = mini
325 delta = np.zeros((n + T.shape[1], n + T.shape[1]))
327 delta[:-2, :-2] = IF.T.dot(IF).toarray() / N
328 delta[-2:, -2:] = T.T.dot(T) / N
331 return Btt.astype(np.float32), P.astype(np.float32)
338 Btt_filt = np.zeros((Btt.shape[0], Btt.shape[1] - nfilt))
339 Btt_filt[:, :Btt_filt.shape[1] - 2] = Btt[:, :Btt.shape[1] - (nfilt + 2)]
340 Btt_filt[:, Btt_filt.shape[1] - 2:] = Btt[:, Btt.shape[1] - 2:]
345 Dmp = np.linalg.inv(Dm.T.dot(Dm)).dot(Dm.T)
347 cmat = Btt_filt.dot(Dmp)
349 return Dm.astype(np.float32), cmat.astype(np.float32)
359 U, s, V = np.linalg.svd(DmtDm)
361 s[s.shape[0] - nfilt - 2:s.shape[0] - 2] = 0.
362 DmtDm1 = U.dot(np.diag(s)).dot(U.T)
363 Dmp = DmtDm1.dot(Dm.T)
367 return Dm.astype(np.float32), cmat.astype(np.float32)
379 def cov_cor(P, noise, trunc, alias, H, bp, tomo):
380 cov = np.zeros((6, 6))
389 for i
in range(cov.shape[0]):
390 for j
in range(cov.shape[1]):
392 tmpi = P.dot(bufdict[str(i)])
393 tmpj = P.dot(bufdict[str(j)])
395 np.mean(tmpi * tmpj, axis=1) -
396 np.mean(tmpi, axis=1) * np.mean(tmpj, axis=1))
398 cov[i, j] = cov[j, i]
400 s = np.reshape(np.diag(cov), (cov.shape[0], 1))
402 cor = cov / np.sqrt(sst)
416 IF = rtc.get_IFsparse(1)
418 noise_com = roket.getContributor(
"noise")
419 trunc_com = roket.getContributor(
"nonlinear")
420 alias_wfs_com = roket.getContributor(
"aliasing")
421 H_com = roket.getContributor(
"filtered")
422 bp_com = roket.getContributor(
"bandwidth")
423 tomo_com = roket.getContributor(
"tomo")
424 fit = roket.getContributor(
"fitting")
426 tmp = (config.p_geom._ipupil.shape[0] -
427 (config.p_dms[0]._n2 - config.p_dms[0]._n1 + 1)) / 2
428 tmp_e0 = config.p_geom._ipupil.shape[0] - tmp
429 tmp_e1 = config.p_geom._ipupil.shape[1] - tmp
430 pup = config.p_geom._ipupil[tmp:tmp_e0, tmp:tmp_e1]
431 indx_pup = np.where(pup.flatten() > 0)[0].astype(np.int32)
432 dm_dim = config.p_dms[0]._n2 - config.p_dms[0]._n1 + 1
433 cov, cor =
cov_cor(P, noise_com, trunc_com, alias_wfs_com, H_com, bp_com, tomo_com)
434 psf = tar.get_image(0,
"le", fluxNorm=
False)
435 psfortho = roket.get_tar_imageortho()
436 covv = roket.get_covv()
437 covm = roket.get_covm()
439 fname =
"/home/fferreira/Data/" + filename
441 "noise": noise_com.T,
442 "aliasing": alias_wfs_com.T,
443 "tomography": tomo_com.T,
444 "filtered modes": H_com.T,
445 "non linearity": trunc_com.T,
446 "bandwidth": bp_com.T,
450 "IF.indices": IF.indices,
451 "IF.indptr": IF.indptr,
454 "indx_pup": indx_pup,
460 "psfortho": psfortho,
464 h5u.save_h5(fname,
"psf", config, psf)
466 for k
in list(pdict.keys()):
467 h5u.save_hdf5(fname, k, pdict[k])
477 nfiltered = config.p_controllers[0].maxcond
478 niters = config.p_loop.niter
481 rtc.load_Btt(1, Btt.dot(Btt.T))
483 rtc.set_cmat(0, cmat)
485 imat = rtc.get_imat(0)
486 RD = np.dot(R, imat).astype(np.float32)
488 gRD = (np.identity(RD.shape[0]) - config.p_controllers[0].gain * RD).astype(np.float32)
489 roket = ao.roket_init(rtc, wfs, tar, dms, tel, atm, 0, 1, Btt.shape[0], Btt.shape[1],
490 nfiltered, niters, Btt, P, gRD, RD)