COMPASS  5.0.0
End-to-end AO simulation tool using GPU acceleration
roket_cpu.py
1 """
2 Created on Wed Apr 27 09:28:23 2016
3 
4 @author: fferreira
5 """
6 
7 import cProfile
8 import pstats as ps
9 
10 import sys, os
11 import numpy as np
12 import carmaWrap as ch
13 import shesha as ao
14 import time
15 import matplotlib.pyplot as pl
16 pl.ion()
17 import hdf5_util as h5u
18 import pandas
19 from scipy.sparse import csr_matrix
20 
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)
24 
25 #get parameters from file
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)
38 else:
39  raise ValueError("Parameter file extension must be .py or .h5")
40 
41 print("param_file is", param_file)
42 
43 if (len(sys.argv) > 2):
44  savename = sys.argv[2]
45 else:
46  savename = "roket_default.h5"
47 
48 
55 
56 if (hasattr(config, "simul_name")):
57  if (config.simul_name is None):
58  simul_name = ""
59  else:
60  simul_name = config.simul_name
61 else:
62  simul_name = ""
63 print("simul name is", simul_name)
64 
65 matricesToLoad = {}
66 if (simul_name == b""):
67  clean = 1
68 else:
69  clean = 0
70  param_dict = h5u.params_dictionary(config)
71  matricesToLoad = h5u.checkMatricesDataBase(os.environ["SHESHA_ROOT"] + "/data/",
72  config, param_dict)
73 #initialisation:
74 # context
75 #c=ch.carmaWrap_context(7)
76 c = ch.carmaWrap_context(devices=np.array([6], dtype=np.int32))
77 #c.set_active_device(device)
78 
79 # wfs
80 print("->wfs")
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)
83 
84 # atmos
85 print("->atmos")
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,
88  load=matricesToLoad)
89 
90 # dm
91 print("->dm")
92 dms = ao.dm_init(config.p_dms, config.p_wfss, wfs, config.p_geom, config.p_tel)
93 
94 # target
95 print("->target")
96 tar = ao.target_init(c, tel, config.p_target, config.p_atmos, config.p_geom,
97  config.p_tel, config.p_dms)
98 
99 print("->rtc")
100 # rtc
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)
104 
105 if not clean:
106  h5u.validDataBase(os.environ["SHESHA_ROOT"] + "/data/", matricesToLoad)
107 
108 print("====================")
109 print("init done")
110 print("====================")
111 print("objects initialzed on GPU:")
112 print("--------------------------------------------------------")
113 print(atm)
114 print(wfs)
115 print(dms)
116 print(tar)
117 print(rtc)
118 
119 print("----------------------------------------------------")
120 print("iter# | SE SR image | LE SR image | Fitting | LE SR phase var")
121 print("----------------------------------------------------")
122 
123 error_flag = True in [w.roket for w in config.p_wfss]
124 
125 
126 
133 def loop(n):
134  """
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.
138 
139  :param n: (int) : number of iterations
140 
141  :return:
142  com : (np.array((n,nactus))) : full command buffer
143 
144  noise_com : (np.array((n,nactus))) : noise contribution for error breakdown
145 
146  alias_wfs_com : (np.array((n,nactus))) : aliasing estimation in the WFS direction
147 
148  tomo_com : (np.array((n,nactus))) : tomography error estimation
149 
150  H_com : (np.array((n,nactus))) : Filtered modes contribution for error breakdown
151 
152  trunc_com : (np.array((n,nactus))) : Truncature and sampling error of WFS
153 
154  bp_com : (np.array((n,nactus))) : Bandwidth error estimation on target
155 
156  mod_com : (np.array((n,nactus))) : Commanded modes expressed on the actuators
157 
158  fit : (float) : fitting (mean variance of the residual target phase after projection)
159 
160  SR : (float) : final strehl ratio returned by the simulation
161  """
162  if (error_flag):
163  # Initialize buffers for error breakdown
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)
174  fit = np.zeros(n)
175  psf_ortho = tar.get_image(0, 'se') * 0.
176  Ee = np.copy(noise_com)
177  Ff = np.copy(Ee)
178  #gamma = 1.0
179  gRD = np.identity(RD.shape[0]) - config.p_controllers[0].gain * gamma * RD
180  t0 = time.time()
181  for i in range(n):
182  atm.move_atmos()
183 
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)
189  tar.dmtrace(0, dms)
190  else:
191  for t in range(config.p_target.ntargets):
192  tar.atmos_trace(t, atm, tel)
193  tar.dmtrace(t, dms)
194  for w in range(len(config.p_wfss)):
195  wfs.sensors_trace(w, "all", tel, atm, dms)
196  wfs.sensors_compimg(w)
197  rtc.docentroids(0)
198  rtc.docontrol(0)
199  # if( i%500==0 and i>0):
200  # #gamma = centroid_gain(Ff[i-500:i,:],Ee[i-500:i,:])
201  # gRD = np.identity(RD.shape[0])-config.p_controllers[0].gain*gamma*RD
202  if (error_flag and i > -1):
203  #compute the error breakdown for this iteration
204  error_breakdown(com, noise_com, alias_wfs_com, tomo_com, H_com,
205  trunc_com, bp_com, wf_com, mod_com, fit, psf_ortho, i,
206  Ee, Ff, gamma, gRD)
207 
208  rtc.applycontrol(0, dms)
209 
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]))
214  t1 = time.time()
215  print(" loop execution time:", t1 - t0, " (", n, "iterations), ", (t1 - t0) / n,
216  "(mean) ", n / (t1 - t0), "Hz")
217  if (error_flag):
218  #Returns the error breakdown
219  SR2 = np.exp(-tar.get_strehl(0, comp_strehl=False)[3])
220  SR = tar.get_strehl(0, comp_strehl=False)[1]
221  #bp_com[-1,:] = bp_com[-2,:]
222  #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
225 
226 
227 def preloop(n):
228  """
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.
232 
233  :param n: (int) : number of iterations
234 
235  :return:
236  com : (np.array((n,nactus))) : full command buffer
237 
238  noise_com : (np.array((n,nactus))) : noise contribution for error breakdown
239 
240  alias_wfs_com : (np.array((n,nactus))) : aliasing estimation in the WFS direction
241 
242  tomo_com : (np.array((n,nactus))) : tomography error estimation
243 
244  H_com : (np.array((n,nactus))) : Filtered modes contribution for error breakdown
245 
246  trunc_com : (np.array((n,nactus))) : Truncature and sampling error of WFS
247 
248  bp_com : (np.array((n,nactus))) : Bandwidth error estimation on target
249 
250  mod_com : (np.array((n,nactus))) : Commanded modes expressed on the actuators
251 
252  fit : (float) : fitting (mean variance of the residual target phase after projection)
253 
254  SR : (float) : final strehl ratio returned by the simulation
255  """
256  for i in range(0, n):
257  atm.move_atmos()
258 
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)
264  else:
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)
270  rtc.docentroids(0)
271  rtc.docontrol(0)
272 
273  rtc.applycontrol(0, dms)
274 
275 
276 
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):
284  """
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)
289 
290  - Ageom : Aliasing contribution on WFS direction
291  Obtained by computing commands from DM orthogonal phase (projection + slopes_geom)
292 
293  - B : Projection on the target direction
294  Obtained as the commmands output of the geometric controller
295 
296  - C : Wavefront
297  Obtained by computing commands from DM parallel phase (RD*B)
298 
299  - E : Wavefront + aliasing + ech/trunc + tomo
300  Obtained by performing the AO loop iteration without noise on the WFS
301 
302  - F : Wavefront + aliasing + tomo
303  Obtained by performing the AO loop iteration without noise on the WFS and using phase deriving slopes
304 
305  - G : tomo
306 
307  Note : rtc.get_err returns to -CMAT.slopes
308 
309  :parameters:
310  noise_com : np.array((niter,nactu)) : Noise contribution
311  Computed with com-E
312 
313  alias_wfs_com : np.array((niter,nactu)) : Aliasing on WFS direction contribution
314  Computed with Ageom
315 
316  tomo_com : np.array((niter,nactu)) : Tomographic error contribution
317  Computed with C-B
318 
319  H_com : np.array((niter,nactu)) : Filtered modes error
320  Computed with B
321 
322  trunc_com : np.array((niter,nactu)) : sampling/truncature error contribution
323  Computed with E-F
324 
325  bp_com : np.array((niter,nactu)) : Bandwidth error
326 
327  wf_com : np.array((niter,nactu)) : Reconstructed wavefront
328 
329  mod_com : np.array((niter,nactu)) : commanded modes
330 
331  fit : np.array((niter)) : fitting value
332 
333  i : (int) : current iteration number
334 
335  """
336  g = config.p_controllers[0].gain
337  Dcom = rtc.get_command(0)
338  Derr = rtc.get_err(0)
339  com[i, :] = Dcom
340  tarphase = tar.get_phase(0)
341 
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"
348  ): # Select the same pixels with or without noise
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)
356 
357  rtc.docentroids(0)
358  if (config.p_centroiders[0].type == b"tcog"):
359  rtc.setthresh(0, config.p_centroiders[0].thresh)
360 
361  rtc.docontrol(0)
362  E = rtc.get_err(0)
363  Ee[i, :] = E
364  # Apply loop filter to get contribution of noise on commands
365  if (i + 1 < config.p_loop.niter):
366  noise_com[i + 1, :] = gRD.dot(noise_com[i, :]) + g * (Derr - E)
367 
368 
371  rtc.docentroids_geom(0)
372  rtc.docontrol(0)
373  F = rtc.get_err(0)
374  Ff[i, :] = F
375  # Apply loop filter to get contribution of sampling/truncature on commands
376  if (i + 1 < config.p_loop.niter):
377  trunc_com[i + 1, :] = gRD.dot(trunc_com[i, :]) + g * (E - gamma * F)
378 
379 
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)
386  """
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)
399  """
400  rtc.docentroids_geom(0)
401  rtc.docontrol(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) # - (E-F))
406 
407 
410  tar.atmos_trace(0, atm, tel)
411  rtc.docontrol_geo(1, dms, tar, 0)
412  B = rtc.get_command(1)
413 
414 
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]
420  if (i >= N_preloop):
421  psf_ortho += tar.get_image(0, 'se') / niters
422 
423 
426  modes = P.dot(B)
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)
432 
433 
436  C = mod_com[i, :] - mod_com[i - 1, :]
437 
438  bp_com[i, :] = gRD.dot(bp_com[i - 1, :]) - C
439 
440 
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)
448  modes = P.dot(G)
449  modes[-nfiltered - 2:-2] = 0
450  wf_com[i, :] = Btt.dot(modes)
451 
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)
455 
456  # Without anyone noticing...
457  tar.set_phase(0, tarphase)
458  rtc.setCom(0, Dcom)
459 
460 
461 def centroid_gain(E, F):
462 
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]
466 
467  return np.mean(cgains)
468 
469 
470 
477  IF = rtc.get_IFsparse(1).T
478  N = IF.shape[0]
479  n = IF.shape[1]
480  #T = IF[:,-2:].copy()
481  T = rtc.get_IFtt(1)
482  #IF = IF[:,:n-2]
483  n = IF.shape[1]
484 
485  delta = IF.T.dot(IF).toarray() / N
486 
487  # Tip-tilt + piston
488  Tp = np.ones((T.shape[0], T.shape[1] + 1))
489  Tp[:, :2] = T.copy() #.toarray()
490  deltaT = IF.T.dot(Tp) / N
491  # Tip tilt projection on the pzt dm
492  tau = np.linalg.inv(delta).dot(deltaT)
493 
494  # Famille generatrice sans tip tilt
495  G = np.identity(n)
496  tdt = tau.T.dot(delta).dot(tau)
497  subTT = tau.dot(np.linalg.inv(tdt)).dot(tau.T).dot(delta)
498  G -= subTT
499 
500  # Base orthonormee sans TT
501  gdg = G.T.dot(delta).dot(G)
502  U, s, V = np.linalg.svd(gdg)
503  U = U[:, :U.shape[1] - 3]
504  s = s[:s.size - 3]
505  L = np.identity(s.size) / np.sqrt(s)
506  B = G.dot(U).dot(L)
507 
508  # Rajout du TT
509  TT = T.T.dot(T) / N #.toarray()/N
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))
513  mini[0, 1] = 0
514  mini[1, 0] = 0
515  Btt[n:, n - 3:] = mini
516 
517  # Calcul du projecteur actus-->modes
518  delta = np.zeros((n + T.shape[1], n + T.shape[1]))
519  #IF = rtc.get_IFsparse(1).T
520  delta[:-2, :-2] = IF.T.dot(IF).toarray() / N
521  delta[-2:, -2:] = T.T.dot(T) / N
522  P = Btt.T.dot(delta)
523 
524  return Btt.astype(np.float32), P.astype(np.float32)
525 
526 
527 def compute_cmatWithBtt(Btt, nfilt):
528  D = rtc.get_imat(0)
529  #D = ao.imat_geom(wfs,config.p_wfss,config.p_controllers[0],dms,config.p_dms,meth=0)
530  # Filtering on Btt modes
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:]
534 
535  # Modal interaction basis
536  Dm = D.dot(Btt_filt)
537  # Direct inversion
538  Dmp = np.linalg.inv(Dm.T.dot(Dm)).dot(Dm.T)
539  # Command matrix
540  cmat = Btt_filt.dot(Dmp)
541 
542  return Dm.astype(np.float32), cmat.astype(np.float32)
543 
544 
545 def compute_cmatWithBtt2(Btt, nfilt):
546  D = rtc.get_imat(0)
547 
548  # Modal interaction basis
549  Dm = D.dot(Btt)
550  # Filtering on modal imat
551  DmtDm = Dm.T.dot(Dm)
552  U, s, V = np.linalg.svd(DmtDm)
553  s = 1. / s
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)
557  # Command matrix
558  cmat = Btt.dot(Dmp)
559 
560  return Dm.astype(np.float32), cmat.astype(np.float32)
561 
562 
563 
570 
571 
572 def cov_cor(P, noise, trunc, alias, H, bp, tomo):
573  cov = np.zeros((6, 6))
574  cor = np.zeros((6, 6))
575  bufdict = {
576  "0": noise.T,
577  "1": trunc.T,
578  "2": alias.T,
579  "3": H.T,
580  "4": bp.T,
581  "5": tomo.T
582  }
583  for i in range(cov.shape[0]):
584  for j in range(cov.shape[1]):
585  if (j >= i):
586  tmpi = P.dot(bufdict[str(i)])
587  tmpj = P.dot(bufdict[str(j)])
588  cov[i, j] = np.sum(
589  np.mean(tmpi * tmpj, axis=1) -
590  np.mean(tmpi, axis=1) * np.mean(tmpj, axis=1))
591  else:
592  cov[i, j] = cov[j, i]
593 
594  s = np.reshape(np.diag(cov), (cov.shape[0], 1))
595  sst = np.dot(s, s.T)
596  ok = np.where(sst)
597  cor[ok] = cov[ok] / np.sqrt(sst[ok])
598 
599  return cov, cor
600 
601 
602 
608 
609 
610 def save_it(filename):
611  IF = rtc.get_IFsparse(1)
612  TT = rtc.get_IFtt(1)
613 
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)
623 
624  fname = "/home/fferreira/Data/" + filename
625  pdict = {
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,
632  "wf_com": wf_com.T,
633  "P": P,
634  "Btt": Btt,
635  "IF.data": IF.data,
636  "IF.indices": IF.indices,
637  "IF.indptr": IF.indptr,
638  "TT": TT,
639  "dm_dim": dm_dim,
640  "indx_pup": indx_pup,
641  "fitting": fit,
642  "SR": SR,
643  "SR2": SR2,
644  "cov": cov,
645  "cor": cor,
646  "psfortho": np.fft.fftshift(psf_ortho),
647  "E": E,
648  "F": F,
649  "dm.xpos": config.p_dms[0]._xpos,
650  "dm.ypos": config.p_dms[0]._ypos,
651  "R": cmat,
652  "Nact": Nact
653  }
654  h5u.save_h5(fname, "psf", config, psf)
655  #h5u.writeHdf5SingleDataset(fname,com.T,datasetName="com")
656  for k in list(pdict.keys()):
657  h5u.save_hdf5(fname, k, pdict[k])
658 
659 
660 
667 nfiltered = int(config.p_controllers[0].maxcond)
668 niters = config.p_loop.niter
669 N_preloop = 1000
670 config.p_loop.set_niter(niters + N_preloop)
671 Btt, P = compute_btt2()
672 rtc.load_Btt(1, Btt.dot(Btt.T))
673 Dm, cmat = compute_cmatWithBtt(Btt, nfiltered)
674 rtc.set_cmat(0, cmat)
675 R = rtc.get_cmat(0)
676 imat = rtc.get_imat(0)
677 RD = np.dot(R, imat)
678 Nact = ao.create_nact_geom(config.p_dms, 0)
679 gamma = 1. / 0.51495
680 #gamma = centroid_gain(100)
681 #print("gamma = ",gamma)
682 
683 #gRD = np.identity(RD.shape[0])-config.p_controllers[0].gain*gamma*RD
684 #diagRD = np.diag(gRD)
685 #gRD = np.diag(diagRD)
686 #gRD=np.diag(gRD)
687 
688 #imat_geom = ao.imat_geom(wfs,config.p_wfss,config.p_controllers[0],dms,config.p_dms,meth=0)
689 #RDgeom = np.dot(R,imat_geom)
690 #preloop(1000)
691 
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(
693  niters + N_preloop)
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:, :]
700 E = E[N_preloop:, :]
701 F = F[N_preloop:, :]
702 save_it(savename)
roket_cpu.loop
def loop(n)
Definition: roket_cpu.py:161
roket_cpu.save_it
def save_it(filename)
Definition: roket_cpu.py:610
roket_cpu.compute_btt2
def compute_btt2()
Definition: roket_cpu.py:476
roket_cpu.compute_cmatWithBtt
def compute_cmatWithBtt(Btt, nfilt)
Definition: roket_cpu.py:527
roket_cpu.error_breakdown
def error_breakdown(com, noise_com, alias_wfs_com, tomo_com, H_com, trunc_com, bp_com, wf_com, mod_com, fit, psf_ortho, i, Ee, Ff, gamma, gRD)
Definition: roket_cpu.py:334
roket_cpu.preloop
def preloop(n)
Definition: roket_cpu.py:255
roket_cpu.centroid_gain
def centroid_gain(E, F)
Definition: roket_cpu.py:461
roket_cpu.cov_cor
def cov_cor(P, noise, trunc, alias, H, bp, tomo)
Definition: roket_cpu.py:572
roket_cpu.compute_cmatWithBtt2
def compute_cmatWithBtt2(Btt, nfilt)
Definition: roket_cpu.py:545