COMPASS  5.0.0
End-to-end AO simulation tool using GPU acceleration
cmats.py
1 
37 
38 import numpy as np
39 import time
40 
41 import shesha.config as conf
42 import shesha.constants as scons
43 
44 from shesha.sutra_wrap import Rtc_FFF as Rtc
45 
46 from shesha.ao.wfs import noise_cov
47 
48 import typing
49 from typing import List
50 
51 
53  M2V: np.ndarray,
54  modalIMat: np.ndarray,
55  modeSelect: np.ndarray = None,
56  modeGains: np.ndarray = None,
57 ) -> np.ndarray:
58  """ Generic numpy modal interaction matrix inversion function
59 
60  :parameters:
61 
62  M2V: (nActu x nModes) : modal basis matrix
63 
64  modalIMat: (nSlopes x nModes) : modal interaction matrix
65 
66  modeSelect: (nModes, dtype=bool): (Optional):
67  mode selection, mode at False is filtered
68 
69  modeGains: (nModes, dtype=bool): (Optional):
70  modal gains to apply. These are gain in the reconstruction sens, ie
71  they are applied multiplicatively on the command matrix
72  """
73  if modeSelect is None:
74  modeSelect = np.ones(modalIMat.shape[1], dtype=bool)
75  if modeGains is None:
76  modeGains = np.ones(modalIMat.shape[1], dtype=np.float32)
77 
78  return M2V.dot(modeGains[:, None] * np.linalg.inv(modalIMat[:, modeSelect].T.dot(
79  modalIMat[:, modeSelect])).dot(modalIMat[:, modeSelect].T))
80 
81 
82 def cmat_init(ncontrol: int, rtc: Rtc, p_controller: conf.Param_controller,
83  p_wfss: List[conf.Param_wfs], p_atmos: conf.Param_atmos,
84  p_tel: conf.Param_tel, p_dms: List[conf.Param_dm],
85  nmodes: int = 0) -> None:
86  """ Compute the command matrix on the GPU
87 
88  :parameters:
89 
90  ncontrol: (int) :
91 
92  rtc: (Rtc) :
93 
94  p_controller: (Param_controller) : controller settings
95 
96  p_wfss: (list of Param_wfs) : wfs settings
97 
98  p_atmos: (Param_atmos) : atmos settings
99 
100  p_tel : (Param_tel) : telescope settings
101 
102  p_dms: (list of Param_dm) : dms settings
103 
104  M2V : (np.ndarray[ndim=2, dtype=np.float32]): (optional) KL to volts matrix (for KL cmat)
105 
106  nmodes: (int) : (optional) number of kl modes
107  """
108  if (p_controller.type == scons.ControllerType.LS):
109  print("Doing imat svd...")
110  t0 = time.time()
111  rtc.d_control[ncontrol].svdec_imat()
112  print("svd done in %f s" % (time.time() - t0))
113  eigenv = np.array(rtc.d_control[ncontrol].d_eigenvals)
114  imat = np.array(rtc.d_control[ncontrol].d_imat)
115  maxcond = p_controller.maxcond
116  if (eigenv[0] < eigenv[eigenv.shape[0] - 1]):
117  mfilt = np.where((eigenv / eigenv[eigenv.shape[0] - 3]) < 1. / maxcond)[0]
118  else:
119  mfilt = np.where((1. / (eigenv / eigenv[2])) > maxcond)[0]
120  nfilt = mfilt.shape[0]
121 
122  print("Building cmat...")
123  t0 = time.time()
124  if not p_controller.do_kl_imat:
125  print("Filtering ", nfilt, " modes")
126  rtc.d_control[ncontrol].build_cmat(nfilt)
127  else:
128  # filter imat
129  D_filt = imat.copy()
130  # Direct inversion
131  Dp_filt = np.linalg.inv(D_filt.T.dot(D_filt)).dot(D_filt.T)
132  if (p_controller.klgain is not None):
133  Dp_filt *= p_controller.klgain[None, :]
134  cmat_filt = p_controller._M2V.dot(Dp_filt)
135  rtc.d_control[ncontrol].set_cmat(cmat_filt)
136 
137  print("cmat done in %f s" % (time.time() - t0))
138 
139  if (p_controller.type == scons.ControllerType.MV):
140  Cn = np.zeros(p_controller._imat.shape[0], dtype=np.float32)
141  ind = 0
142  for k in p_controller.nwfs:
143  Cn[ind:ind + 2 * p_wfss[k]._nvalid] = noise_cov(k, p_wfss[k], p_atmos, p_tel)
144  ind += 2 * p_wfss[k]._nvalid
145 
146  rtc.d_control[ncontrol].load_noisemat(Cn)
147  print("Building cmat...")
148  rtc.d_control[ncontrol].build_cmat(p_controller.maxcond)
149 
150  if (p_controller.TTcond == None):
151  p_controller.set_TTcond(p_controller.maxcond)
152 
153  if ("tt" in [dm.type for dm in p_dms]):
154  rtc.d_control[ncontrol].filter_cmat(p_controller.TTcond)
155  print("Done")
156  p_controller.set_cmat(np.array(rtc.d_control[ncontrol].d_cmat))
157 
158 
160  DtD = D.T.dot(D)
161  return np.linalg.svd(DtD)
162 
163 
164 def Btt_for_cmat(rtc, dms, p_dms, p_geom):
165  """ Compute a command matrix in Btt modal basis (see error breakdown) and set
166  it on the sutra_rtc. It computes by itself the volts to Btt matrix.
167 
168  :parameters:
169 
170  rtc: (Rtc) : rtc object
171 
172  dms: (Dms): dms object
173 
174  p_dms: (list of Param_dm): dms settings
175 
176  p_geom: (Param_geom): geometry settings
177 
178  """
179 
180  IFs = basis.compute_IFsparse(dms, p_dms, p_geom).T
181  n = IFs.shape[1]
182  IFtt = IFs[:, -2:].toarray()
183  IFpzt = IFs[:, :n - 2]
184 
185  Btt, P = basis.compute_btt(IFpzt, IFtt)
186  return Btt, P
187 
188 
189 def get_cmat(D, nfilt, Btt=None, rtc=None, svd=None):
190  """Compute a command matrix from an interaction matrix 'D'
191 
192  usage:
193  get_cmat(D,nfilt)
194  get_cmat(D,nfilt,Btt=BTT,rtc=RTC)
195  get_cmat(D,nfilt,svd=SVD)
196 
197  :parameters:
198  D: (np.ndarray[ndim=2, dtype=np.float32]): interaction matrix
199 
200  nfilt: (int): number of element to filter
201 
202  Btt: (np.ndarray[ndim=2, dtype=np.float32]): Btt modal basis
203 
204  rtc: (Rtc) :
205 
206  svd: (tuple of np.ndarray[ndim=1, dtype=np.float32): svd of D.T*D (obtained from np.linalg.svd)
207  """
208  nfilt = max(nfilt, 0) #nfilt is positive
209  if (Btt is not None):
210  if (svd is not None):
211  raise ValueError("Btt and SVD cannt be used together")
212  if (rtc is None):
213  raise ValueError("Btt cannot be used without rtc")
214  n = Btt.shape[1]
215  index = np.concatenate((np.arange(n - nfilt - 2), np.array([n - 2,
216  n - 1]))).astype(int)
217  Btt_filt = Btt[:, index]
218  # Modal interaction basis
219  Dm = D.dot(Btt_filt)
220  # Direct inversion
221  Dmp = np.linalg.inv(Dm.T.dot(Dm)).dot(Dm.T)
222  # Command matrix
223  cmat = Btt_filt.dot(Dmp)
224  else:
225  if (svd is not None):
226  u = svd[0]
227  s = svd[1]
228  v = svd[2]
229  else:
230  u, s, v = svd_for_cmat(D)
231  s_filt = 1 / s
232  if (nfilt > 0):
233  s_filt[-nfilt:] = 0
234  DtDx = v.T.dot(np.diag(s_filt)).dot(u.T)
235  cmat = DtDx.dot(D.T)
236 
237  return cmat.astype(np.float32)
shesha.ao.cmats.svd_for_cmat
def svd_for_cmat(D)
Definition: cmats.py:159
shesha.sutra_wrap
Definition: sutra_wrap.py:1
shesha.ao.cmats.generic_imat_inversion
np.ndarray generic_imat_inversion(np.ndarray M2V, np.ndarray modalIMat, np.ndarray modeSelect=None, np.ndarray modeGains=None)
Generic numpy modal interaction matrix inversion function.
Definition: cmats.py:67
shesha.constants
Numerical constants for shesha and config enumerations for safe-typing.
Definition: constants.py:1
shesha.config
Parameter classes for COMPASS.
Definition: shesha/shesha/config/__init__.py:1
shesha.ao.wfs
On the fly modification of the WFS.
Definition: wfs.py:1
shesha.ao.cmats.cmat_init
None cmat_init(int ncontrol, Rtc rtc, conf.Param_controller p_controller, List[conf.Param_wfs] p_wfss, conf.Param_atmos p_atmos, conf.Param_tel p_tel, List[conf.Param_dm] p_dms, int nmodes=0)
Compute the command matrix on the GPU.
Definition: cmats.py:104
shesha.ao.cmats.get_cmat
def get_cmat(D, nfilt, Btt=None, rtc=None, svd=None)
Compute a command matrix from an interaction matrix 'D'.
Definition: cmats.py:207
shesha.ao.cmats.Btt_for_cmat
def Btt_for_cmat(rtc, dms, p_dms, p_geom)
Compute a command matrix in Btt modal basis (see error breakdown) and set it on the sutra_rtc.
Definition: cmats.py:178