45 from scipy.sparse
import csr_matrix
47 from typing
import List
48 from tqdm
import trange
51 def compute_KL2V(p_controller: conf.Param_controller, dms: Dms, p_dms: list,
52 p_geom: conf.Param_geom, p_atmos: conf.Param_atmos,
53 p_tel: conf.Param_tel):
54 """ Compute the Karhunen-Loeve to Volt matrix
55 (transfer matrix between the KL space and volt space for a pzt dm)
59 p_controller: (Param_controller) : p_controller settings
61 dms : (shesha_dms) : Dms object
63 p_dms: (list of Param_dm) : dms settings
65 p_geom : (Param_geom) : geometry parameters
67 p_atmos : (Param_atmos) : atmos parameters
69 p_tel : (Param_tel) : telescope parameters
73 KL2V : (np.array(np.float32,dim=2)) : KL to Volt matrix
75 ntotact = np.array([p_dms[i]._ntotact
for i
in range(len(p_dms))], dtype=np.int64)
76 KL2V = np.zeros((np.sum(ntotact), np.sum(ntotact)), dtype=np.float32)
81 for i
in range(p_controller.ndm.size):
82 ndm = p_controller.ndm[i]
83 if (p_dms[ndm].type == scons.DmType.PZT):
84 tmp = (p_geom._ipupil.shape[0] - (p_dms[ndm]._n2 - p_dms[ndm]._n1 + 1)) // 2
85 tmp_e0 = p_geom._ipupil.shape[0] - tmp
86 tmp_e1 = p_geom._ipupil.shape[1] - tmp
87 pup = p_geom._ipupil[tmp:tmp_e0, tmp:tmp_e1]
88 indx_valid = np.where(pup.flatten(
"F") > 0)[0].astype(np.int32)
89 p2m = p_tel.diam / p_geom.pupdiam
90 norm = -(p2m * p_tel.diam / (2 * p_atmos.r0))**(5. / 3)
92 dms.d_dms[ndm].compute_KLbasis(p_dms[ndm]._xpos, p_dms[ndm]._ypos,
93 indx_valid, indx_valid.size, norm, 1.0)
95 KL2V[indx_act:indx_act + ntotact[ndm], indx_act:indx_act + ntotact[ndm]] = \
96 np.fliplr(dms.d_dms[ndm].d_KLbasis)
97 indx_act += ntotact[ndm]
98 elif (p_dms[ndm].type == scons.DmType.TT):
100 if (p_controller.nmodes
is not None and
101 p_controller.nmodes < KL2V.shape[1] - 2 * nTT):
102 KL2V = KL2V[:, :p_controller.nmodes]
104 KL2V = KL2V[:, :KL2V.shape[1] - 2 * nTT]
106 raise ValueError(
"More than 1 TipTilt found! Stupid")
108 KL2V[:, :KL2V.shape[1] - 2] = KL2V[:, 2:]
109 KL2V[:, KL2V.shape[1] - 2:] = np.zeros((np.sum(ntotact), 2), dtype=np.float32)
110 KL2V[np.sum(ntotact) - 2:, KL2V.shape[1] - 2:] = np.identity(2, dtype=np.float32)
116 """ Compute a the DM basis as a sparse matrix :
117 - push on each actuator
118 - get the corresponding dm shape
119 - apply pupil mask and store in a column
122 g_dm: (Dm) : Dm object
124 p_dm: (Param_dm) : dm settings
126 p_geom: (Param_geom) : geom settings
130 IFbasis = (csr_matrix) : DM IF basis
132 tmp = (p_geom._ipupil.shape[0] - (p_dm._n2 - p_dm._n1 + 1)) // 2
133 tmp_e0 = p_geom._ipupil.shape[0] - tmp
134 tmp_e1 = p_geom._ipupil.shape[1] - tmp
135 pup = p_geom._ipupil[tmp:tmp_e0, tmp:tmp_e1]
136 indx_valid = np.where(pup.flatten(
"F") > 0)[0].astype(np.int32)
139 for i
in trange(p_dm._ntotact):
141 g_dm.comp_oneactu(i, 1.0)
142 shape = np.array(g_dm.d_shape)
143 IFvec = csr_matrix(shape.flatten(
"F")[indx_valid])
147 row = np.append(0, IFvec.getnnz())
149 val = np.append(val, IFvec.data)
150 col = np.append(col, IFvec.indices)
151 row = np.append(row, row[-1] + IFvec.getnnz())
153 IFbasis = csr_matrix((val, col, row))
158 """ Compute the influence functions of all DMs as a sparse matrix :
159 - push on each actuator
160 - get the corresponding dm shape
161 - apply pupil mask and store in a column
165 g_dm: (Dms) : Dms object
167 p_dms: (Param_dms) : dms settings
169 p_geom: (Param_geom) : geom settings
173 IFbasis = (csr_matrix) : DM IF basis
183 val = np.append(val, IFi.data)
184 col = np.append(col, IFi.indices)
185 row = np.append(row, row[-1] + IFi.indptr[1:])
186 IFsparse = csr_matrix((val, col, row))
190 def command_on_Btt(rtc: Rtc, dms: Dms, p_dms: list, p_geom: conf.Param_geom, nfilt: int):
191 """ Compute a command matrix in Btt modal basis (see error breakdown) and set
192 it on the sutra_rtc. It computes by itself the volts to Btt matrix.
196 rtc: (Rtc) : rtc object
198 dms: (Dms): dms object
200 p_dms: (list of Param_dm): dms settings
202 p_geom: (Param_geom): geometry settings
204 nfilt: (int): number of modes to filter
209 IFtt = IFs[:, -2:].copy().toarray()
210 IFpzt = IFs[:, :n - 2]
217 """ Compute a command matrix on the Btt basis and load it in the GPU
221 rtc: (Rtc): rtc object
223 Btt: (np.ndarray[ndim=2, dtype=np.float32]) : volts to Btt matrix
225 nfilt: (int): number of modes to filter
227 D = np.array(rtc.d_control[0].d_imat)
230 Btt_filt = np.zeros((Btt.shape[0], Btt.shape[1] - nfilt))
231 Btt_filt[:, :Btt_filt.shape[1] - 2] = Btt[:, :Btt.shape[1] - (nfilt + 2)]
232 Btt_filt[:, Btt_filt.shape[1] - 2:] = Btt[:, Btt.shape[1] - 2:]
237 Dmp = np.linalg.inv(Dm.T.dot(Dm)).dot(Dm.T)
239 cmat = Btt_filt.dot(Dmp)
240 rtc.d_control[0].set_cmat(cmat.astype(np.float32))
242 return cmat.astype(np.float32)
245 def command_on_KL(rtc: Rtc, dms: Dms, p_controller: conf.Param_controller,
246 p_dms: List[conf.Param_dm], p_geom: conf.Param_geom,
247 p_atmos: conf.Param_atmos, p_tel: conf.Param_tel, nfilt: int):
248 """ Compute a command matrix in KL modal basis and set
249 it on the sutra_rtc. It computes by itself the volts to KL matrix.
253 rtc: (Rtc) : rtc object
255 dms: (Dms): dms object
257 p_dms: (list of Param_dm): dms settings
259 p_geom: (Param_geom): geometry settings
261 p_atmos : (Param_atmos) : atmos parameters
263 p_tel : (Param_tel) : telescope parameters
265 nfilt: (int): number of modes to filter
267 KL2V =
compute_KL2V(p_controller, dms, p_dms, p_geom, p_atmos, p_tel)
272 """ Compute a command matrix on the KL basis and load it in the GPU
276 rtc: (Rtc): rtc object
278 KL2V: (np.ndarray[ndim=2, dtype=np.float32]) : volts to KL matrix
280 nfilt: (int): number of modes to filter
282 D = np.array(rtc.d_control[0].d_imat)
283 KL2V_filt = np.zeros((KL2V.shape[0], KL2V.shape[1] - nfilt))
284 KL2V_filt[:, :KL2V_filt.shape[1] - 2] = KL2V[:, :KL2V.shape[1] - (nfilt + 2)]
285 KL2V_filt[:, KL2V_filt.shape[1] - 2:] = KL2V[:, KL2V.shape[1] - 2:]
288 Dm = D.dot(KL2V_filt)
290 Dmp = np.linalg.inv(Dm.T.dot(Dm)).dot(Dm.T)
292 cmat = KL2V_filt.dot(Dmp)
293 rtc.d_control[0].set_cmat(cmat.astype(np.float32))
295 return cmat.astype(np.float32)
299 actu_y_pos: np.ndarray, periodic=
'n'):
301 Values you are looking for are:
315 elif periodic ==
'n-1':
318 raise ValueError(
'periodic can only be "n" or "n-1" to set boundary condition.')
319 xnorm = (np.round((actu_x_pos - np.min(actu_x_pos)) / pitch).astype(np.int32)) % n
320 ynorm = (np.round((actu_y_pos - np.min(actu_y_pos)) / pitch).astype(np.int32)) % n
323 data = np.zeros((n, n, n, n), dtype=np.float32)
326 data[i, j, i, j] = 1.
328 data = np.fft.fftn(data, axes=(2, 3))
330 takeSine = np.zeros((n, n), dtype=bool)
331 takeSine[0, n // 2 + 1:] = 1
333 takeSine[n // 2, n // 2 + 1:] = 1
335 takeSine[n // 2 + 1:, :] = 1
341 actuPush = data[:, :, xnorm, ynorm]
348 def compute_btt(IFpzt, IFtt, influ_petal=None, return_delta=False):
349 """ Returns Btt to Volts and Volts to Btt matrices
353 IFpzt : (csr_matrix) : influence function matrix of pzt DM, sparse and arrange as (Npts in pup x nactus)
355 IFtt : (np.ndarray(ndim=2,dtype=np.float32)) : Influence function matrix of the TT mirror arrange as (Npts in pup x 2)
357 influ_petal : (np.ndarray) : Influence function matrix of M4 petals.
358 Default is None, if set, the Btt produced is also orthogonal
359 to petal modes, then only driven by petal DM
361 return_delta : (bool, optional) : If True, returns delta instead of P. Default is False
365 Btt : (np.ndarray(ndim=2,dtype=np.float32)) : Btt to Volts matrix
367 P : (np.ndarray(ndim=2,dtype=np.float32)) : Volts to Btt matrix
372 raise ValueError(
"Influence functions must be arrange as (Npts_pup x nactus)")
374 delta = IFpzt.T.dot(IFpzt).toarray() / N
377 Tp = np.ones((IFtt.shape[0], IFtt.shape[1] + 1))
378 Tp[:, :2] = IFtt.copy(
380 deltaT = IFpzt.T.dot(Tp) / N
382 tau = np.linalg.inv(delta).dot(deltaT)
385 if influ_petal
is not None:
387 nseg = influ_petal.toarray().shape[0]
388 petal_modes = -1 / (nseg - 1) * np.ones((nseg, (nseg - 1)))
389 petal_modes += nseg / (nseg - 1) * np.eye(nseg)[:, 0:(
391 tau_petal = IFpzt.T.dot(influ_petal).toarray().dot(petal_modes.T)
392 tau = np.concatenate((tau_petal, np.linalg.inv(delta).dot(deltaT)), axis=1)
396 tdt = tau.T.dot(delta).dot(tau)
397 subTT = tau.dot(np.linalg.inv(tdt)).dot(tau.T).dot(delta)
401 gdg = G.T.dot(delta).dot(G)
402 U, s, V = np.linalg.svd(gdg)
403 U = U[:, :U.shape[1] - nfilt]
404 s = s[:s.size - nfilt]
405 L = np.identity(s.size) / np.sqrt(s)
409 TT = IFtt.T.dot(IFtt) / N
410 Btt = np.zeros((n + 2, n - 1))
411 Btt[:B.shape[0], :B.shape[1]] = B
412 mini = 1. / np.sqrt(np.abs(TT))
416 if influ_petal
is not None:
417 Btt[:n, -7:-2] = tau_petal
420 Delta = np.zeros((n + IFtt.shape[1], n + IFtt.shape[1]))
422 Delta[:-2, :-2] = delta
429 return Btt.astype(np.float32), P.astype(np.float32)