42 from shesha.ao import imats, cmats, tomo, basis, modopti
46 from typing
import List
49 from shesha.sutra_wrap import (carmaWrap_context, Sensors, Dms, Target, Rtc_brahma,
50 Rtc_cacao_FFF, Atmos, Telescope)
54 def rtc_init(context: carmaWrap_context, tel: Telescope, wfs: Sensors, dms: Dms,
55 atmos: Atmos, p_wfss: list, p_tel: conf.Param_tel, p_geom: conf.Param_geom,
56 p_atmos: conf.Param_atmos, ittime: float, p_centroiders=
None,
57 p_controllers=
None, p_dms=
None, do_refslp=
False, brahma=
False, cacao=
False,
58 tar=
None, dataBase={}, use_DB=
False):
59 """Initialize all the SutraRtc objects : centroiders and controllers
62 context: (carmaWrap_context): context
63 tel: (Telescope) : Telescope object
64 wfs: (Sensors) : Sensors object
65 dms: (Dms) : Dms object
66 atmos: (Atmos) : Atmos object
67 p_wfss: (list of Param_wfs) : wfs settings
68 p_tel: (Param_tel) : telescope settings
69 p_geom: (Param_geom) : geom settings
70 p_atmos: (Param_atmos) : atmos settings
71 ittime: (float) : iteration time [s]
72 p_centroiders : (list of Param_centroider): (optional) centroiders settings
73 p_controllers : (list of Param_controller): (optional) controllers settings
74 p_dms: (list of Param_dms) : (optional) dms settings
75 do_refslp : (bool): (optional) do ref slopes flag, default=False
76 brahma: (bool) : (optional) brahma flag
77 cacao: (bool) : (optional) cacao flag
78 tar: (Target) : (optional)
79 dataBase: (dict): (optional) dict containig paths to files to load
80 use_DB: (bool): use dataBase flag
82 Rtc : (Rtc) : Rtc object
87 rtc = Rtc_brahma(context, wfs, tar,
"rtc_brahma")
89 rtc = Rtc_cacao_FFF(
"compass_calPix",
"compass_loopData")
97 ncentro = len(p_centroiders)
102 ncontrol = len(p_controllers)
106 if p_centroiders
is not None:
107 for i
in range(ncentro):
108 nwfs = p_centroiders[i].nwfs
112 if p_controllers
is not None:
113 if (p_wfss
is not None and p_dms
is not None):
114 for i
in range(ncontrol):
115 if not "dm" in dataBase:
116 imat = imats.imat_geom(wfs, dms, p_wfss, p_dms, p_controllers[i],
121 if p_dms[0].type == scons.DmType.PZT:
122 dm_init.correct_dm(context, dms, p_dms, p_controllers[i], p_geom,
123 imat, dataBase=dataBase, use_DB=use_DB)
126 p_atmos, ittime, p_tel, rtc, dms, wfs, tel, atmos,
127 p_centroiders, do_refslp, dataBase=dataBase,
131 roket_flag =
True in [w.roket
for w
in p_wfss]
133 p_controller = p_controllers[0]
134 Nphi = np.where(p_geom._spupil)[0].size
136 list_dmseen = [p_dms[j].type
for j
in p_controller.ndm]
137 nactu = np.sum([p_dms[j]._ntotact
for j
in p_controller.ndm])
139 rtc.add_controller(context, p_controller.nvalid, p_controller.nslope,
140 p_controller.nactu, p_controller.delay,
141 context.active_device, scons.ControllerType.GEO, dms,
142 p_controller.ndm, p_controller.ndm.size,
143 p_controller.nwfs, p_controller.nwfs.size, Nphi,
True)
157 def rtc_standalone(context: carmaWrap_context, nwfs: int, nvalid: list, nactu: int,
158 centroider_type: list, delay: list, offset: list, scale: list,
159 brahma: bool =
False, fp16: bool =
False, cacao: bool =
False) -> Rtc:
163 print(
"start rtc_standalone")
165 rtc = Rtc_brahma(context,
None,
None,
"rtc_brahma")
169 rtc = Rtc_cacao_FHF(
"compass_calPix",
"compass_loopData")
171 rtc = Rtc_cacao_FFF(
"compass_calPix",
"compass_loopData")
178 for k
in range(nwfs):
181 rtc.add_centroider(context, nvalid[k], offset[k], scale[k],
False,
182 context.active_device, centroider_type[k])
184 nslopes = sum([c.nslopes
for c
in rtc.d_centro])
185 rtc.add_controller(context, sum(nvalid), nslopes, nactu, delay[k],
186 context.active_device,
"generic", idx_centro=np.arange(nwfs),
189 print(
"rtc_standalone set")
194 p_centroider: conf.Param_centroider, p_tel: conf.Param_tel,
195 p_atmos: conf.Param_atmos, wfs: Sensors, rtc: Rtc):
196 """ Initialize a centroider object in Rtc
199 context: (carmaWrap_context): context
200 nwfs : (int) : index of wfs
201 p_wfs : (Param_wfs): wfs settings
202 p_centroider : (Param_centroider) : centroider settings
203 wfs: (Sensors): Sensor object
204 rtc : (Rtc) : Rtc object
206 if (p_wfs.type == scons.WFSType.SH):
207 if (p_centroider.type != scons.CentroiderType.CORR):
208 s_offset = p_wfs.npix // 2. - 0.5
210 if (p_centroider.type_fct == scons.CentroiderFctType.MODEL):
211 if (p_wfs.npix % 2 == 0):
212 s_offset = p_wfs.npix // 2 - 0.5
214 s_offset = p_wfs.npix // 2
216 s_offset = p_wfs.npix // 2 - 0.5
217 s_scale = p_wfs.pixsize
219 elif (p_wfs.type == scons.WFSType.PYRHR
or p_wfs.type == scons.WFSType.PYRLR):
221 s_scale = (p_wfs.Lambda * 1e-6 / p_tel.diam) * \
222 p_wfs.pyr_ampl * CONST.RAD2ARCSEC
224 rtc.add_centroider(context, p_wfs._nvalid, s_offset, s_scale, p_centroider.filter_TT,
225 context.active_device, p_centroider.type, wfs.d_wfs[nwfs])
226 rtc.d_centro[-1].load_validpos(p_wfs._validsubsx, p_wfs._validsubsy,
227 p_wfs._nvalid * p_wfs.nPupils)
229 rtc.d_centro[-1].set_npix(p_wfs.npix)
231 if (p_centroider.type != scons.CentroiderType.MASKEDPIX):
232 p_centroider._nslope = 2 * p_wfs._nvalid
234 p_centroider._nslope = p_wfs._validsubsx.size
236 if (p_centroider.type == scons.CentroiderType.PYR):
238 rtc.d_centro[nwfs].set_pyr_method(p_centroider.method)
239 rtc.d_centro[nwfs].set_pyr_thresh(p_centroider.thresh)
241 elif (p_wfs.type == scons.WFSType.SH):
242 if (p_centroider.type == scons.CentroiderType.TCOG):
243 rtc.d_centro[nwfs].set_threshold(p_centroider.thresh)
244 elif (p_centroider.type == scons.CentroiderType.BPCOG):
245 rtc.d_centro[nwfs].set_nmax(p_centroider.nmax)
246 elif (p_centroider.type == scons.CentroiderType.WCOG
or
247 p_centroider.type == scons.CentroiderType.CORR):
248 r0 = p_atmos.r0 * (p_wfs.Lambda / 0.5)**(6 / 5.)
249 seeing = CONST.RAD2ARCSEC * (p_wfs.Lambda * 1.e-6) / r0
250 npix = seeing // p_wfs.pixsize
252 if p_centroider.type == scons.CentroiderType.WCOG:
253 rtc.d_centro[nwfs].init_weights()
254 rtc.d_centro[nwfs].load_weights(p_centroider.weights,
255 p_centroider.weights.ndim)
257 corrnorm = np.ones((2 * p_wfs.npix, 2 * p_wfs.npix), dtype=np.float32)
258 p_centroider.sizex = 3
259 p_centroider.sizey = 3
260 p_centroider.interpmat = rtc_util.create_interp_mat(
261 p_centroider.sizex, p_centroider.sizey).astype(np.float32)
263 if (p_centroider.weights
is None):
264 raise ValueError(
"p_centroider.weights is None")
265 rtc.d_centro[nwfs].init_corr(p_centroider.sizex, p_centroider.sizey,
266 p_centroider.interpmat)
267 rtc.d_centro[nwfs].load_corr(p_centroider.weights, corrnorm,
268 p_centroider.weights.ndim)
271 def comp_weights(p_centroider: conf.Param_centroider, p_wfs: conf.Param_wfs, npix: int):
273 Compute the weights used by centroider wcog and corr
276 p_centroider : (Param_centroider) : centroider settings
277 p_wfs : (Param_wfs) : wfs settings
280 if (p_centroider.type_fct == scons.CentroiderFctType.MODEL):
282 if (p_wfs.gsalt > 0):
284 tmp2 = utilities.makegaussian(tmp.shape[1],
285 npix * p_wfs._nrebin).astype(np.float32)
286 tmp3 = np.zeros((tmp.shape[1], tmp.shape[1], p_wfs._nvalid),
289 for j
in range(p_wfs._nvalid):
290 tmp3[:, :, j] = np.fft.ifft2(
291 np.fft.fft2(tmp[:, :, j]) * np.fft.fft2(tmp2.T)).real
292 tmp3[:, :, j] *= tmp3.shape[0] * tmp3.shape[1]
293 tmp3[:, :, j] = np.fft.fftshift(tmp3[:, :, j])
295 offset = (p_wfs._Ntot - p_wfs._nrebin * p_wfs.npix) // 2
296 j = offset + p_wfs._nrebin * p_wfs.npix
297 tmp = np.zeros((j - offset + 1, j - offset + 1, tmp3.shape[2]),
299 tmp3 = np.cumsum(tmp3[offset:j, offset:j, :], axis=0)
300 tmp[1:, 1:, :] = np.cumsum(tmp3, axis=1)
301 tmp = np.diff(tmp[::p_wfs._nrebin, ::p_wfs._nrebin, :], axis=0)
302 tmp = np.diff(tmp, axis=1)
304 p_centroider.weights = tmp
306 p_centroider.type_fct = scons.CentroiderFctType.GAUSS
307 print(
"No LGS found, centroider weighting function becomes gaussian")
309 if (p_centroider.type_fct == scons.CentroiderFctType.GAUSS):
310 if p_centroider.width
is None:
311 p_centroider.width = npix
312 if (p_wfs.npix % 2 == 1):
313 p_centroider.weights = utilities.makegaussian(
314 p_wfs.npix, p_centroider.width, p_wfs.npix // 2,
315 p_wfs.npix // 2).astype(np.float32)
316 elif (p_centroider.type == scons.CentroiderType.CORR):
317 p_centroider.weights = utilities.makegaussian(
318 p_wfs.npix, p_centroider.width, p_wfs.npix // 2,
319 p_wfs.npix // 2).astype(np.float32)
321 p_centroider.weights = utilities.makegaussian(
322 p_wfs.npix, p_centroider.width, p_wfs.npix // 2 - 0.5,
323 p_wfs.npix // 2 - 0.5).astype(np.float32)
326 def init_controller(context, i: int, p_controller: conf.Param_controller, p_wfss: list,
327 p_geom: conf.Param_geom, p_dms: list, p_atmos: conf.Param_atmos,
328 ittime: float, p_tel: conf.Param_tel, rtc: Rtc, dms: Dms,
329 wfs: Sensors, tel: Telescope, atmos: Atmos,
330 p_centroiders: List[conf.Param_centroider], do_refslp=
False,
331 dataBase={}, use_DB=
False):
333 Initialize the controller part of rtc
336 context: (carmaWrap_context): context
337 i : (int) : controller index
338 p_controller: (Param_controller) : controller settings
339 p_wfss: (list of Param_wfs) : wfs settings
340 p_geom: (Param_geom) : geom settings
341 p_dms: (list of Param_dms) : dms settings
342 p_atmos: (Param_atmos) : atmos settings
343 ittime: (float) : iteration time [s]
344 p_tel: (Param_tel) : telescope settings
345 rtc: (Rtc) : Rtc objet
346 dms: (Dms) : Dms object
347 wfs: (Sensors) : Sensors object
348 tel: (Telescope) : Telescope object
349 atmos: (Atmos) : Atmos object
350 p_centroiders: (list of Param_centroider): centroiders settings
352 if (p_controller.type != scons.ControllerType.GEO):
353 nwfs = p_controller.nwfs
354 if (len(p_wfss) == 1):
355 nwfs = p_controller.nwfs
357 nvalid = sum([p_wfss[k]._nvalid
for k
in nwfs])
358 p_controller.set_nvalid(int(np.sum([p_wfss[k]._nvalid
for k
in nwfs])))
360 for c
in p_centroiders:
362 tmp = tmp + c._nslope
363 p_controller.set_nslope(int(tmp))
365 nslope = np.sum([c._nslope
for c
in p_centroiders])
366 p_controller.set_nslope(int(nslope))
369 ndms = p_controller.ndm.tolist()
370 nactu = np.sum([p_dms[j]._ntotact
for j
in ndms])
371 p_controller.set_nactu(int(nactu))
373 alt = np.array([p_dms[j].alt
for j
in p_controller.ndm], dtype=np.float32)
375 list_dmseen = [p_dms[j].type
for j
in p_controller.ndm]
376 if (p_controller.type == scons.ControllerType.GEO):
377 Nphi = np.where(p_geom._spupil)[0].size
385 rtc.add_controller(context, p_controller.nvalid, p_controller.nslope,
386 p_controller.nactu, p_controller.delay, context.active_device,
387 p_controller.type, dms, p_controller.ndm, p_controller.ndm.size,
388 p_controller.nwfs, p_controller.nwfs.size, Nphi,
False,
389 p_controller.nstates)
390 print(
"CONTROLLER ADDED")
391 if (p_wfss
is not None and do_refslp):
392 rtc.do_centroids_ref(i)
394 if (p_controller.type == scons.ControllerType.GEO):
397 if (p_controller.type == scons.ControllerType.LS):
399 p_tel, rtc, dms, wfs, tel, atmos, dataBase=dataBase,
402 if (p_controller.type == scons.ControllerType.CURED):
405 if (p_controller.type == scons.ControllerType.MV):
409 elif (p_controller.type == scons.ControllerType.GENERIC):
412 p_controller._imat = imats.imat_geom(wfs, dms, p_wfss, p_dms, p_controller,
415 print(
"p_controller._imat not set")
419 p_controller: conf.Param_controller, p_dms: list, roket=
False):
421 Initialize geometric controller
424 i: (int): controller index
425 rtc: (Rtc): rtc object
426 dms: (Dms): Dms object
427 p_geom: (Param_geom): geometry settings
428 p_controller: (Param_controller): controller settings
429 p_dms: (list of Param_dms): dms settings
430 roket: (bool): Flag to initialize ROKET
432 indx_pup = np.where(p_geom._spupil.flatten(
'F'))[0].astype(np.int32)
433 indx_mpup = np.where(p_geom._mpupil.flatten(
'F'))[0].astype(np.int32)
435 indx_dm = np.zeros((p_controller.ndm.size * indx_pup.size), dtype=np.int32)
436 for dmn
in range(p_controller.ndm.size):
437 tmp_s = (p_geom._ipupil.shape[0] - (p_dms[dmn]._n2 - p_dms[dmn]._n1 + 1)) // 2
438 tmp_e0 = p_geom._ipupil.shape[0] - tmp_s
439 tmp_e1 = p_geom._ipupil.shape[1] - tmp_s
440 pup_dm = p_geom._ipupil[tmp_s:tmp_e0, tmp_s:tmp_e1]
441 indx_dm[cpt:cpt + np.where(pup_dm)[0].size] = np.where(pup_dm.flatten(
'F'))[0]
442 cpt += np.where(pup_dm)[0].size
444 unitpervolt = np.array([p_dms[j].unitpervolt
445 for j
in range(len(p_dms))], dtype=np.float32)
447 rtc.d_control[i].init_proj_sparse(dms, indx_dm, unitpervolt, indx_pup, indx_mpup,
452 p_geom: conf.Param_geom, p_dms: list, p_atmos: conf.Param_atmos,
453 ittime: float, p_tel: conf.Param_tel, rtc: Rtc, dms: Dms,
454 wfs: Sensors, tel: Telescope, atmos: Atmos, dataBase: dict = {},
455 use_DB: bool =
False):
457 Initialize the least square controller
459 i : (int) : controller index
460 p_controller: (Param_controller) : controller settings
461 p_wfss: (list of Param_wfs) : wfs settings
462 p_geom: (Param_geom) : geom settings
463 p_dms: (list of Param_dms) : dms settings
464 p_atmos: (Param_atmos) : atmos settings
465 ittime: (float) : iteration time [s]
466 p_tel: (Param_tel) : telescope settings
467 rtc: (Rtc) : Rtc objet
468 dms: (Dms) : Dms object
469 wfs: (Sensors) : Sensors object
470 tel: (Telescope) : Telescope object
471 atmos: (Atmos) : Atmos object
474 if p_controller.do_kl_imat:
475 IF = basis.compute_IFsparse(dms, p_dms, p_geom).T
476 M2V, _ = basis.compute_btt(IF[:, :-2], IF[:, -2:].toarray())
477 print(
"Filtering ", p_controller.nModesFilt,
" modes based on mode ordering")
478 M2V = M2V[:, list(range(M2V.shape[1] - 2 - p_controller.nModesFilt)) + [-2, -1]]
480 if len(p_controller.klpush) == 1:
481 p_controller.klpush = p_controller.klpush[0] * np.ones(M2V.shape[1])
482 imats.imat_init(i, rtc, dms, p_dms, wfs, p_wfss, p_tel, p_controller, M2V,
483 dataBase=dataBase, use_DB=use_DB)
485 if p_controller.modopti:
486 print(
"Initializing Modal Optimization : ")
487 p_controller.nrec = int(2**np.ceil(np.log2(p_controller.nrec)))
488 if p_controller.nmodes
is None:
489 p_controller.nmodes = sum([p_dms[j]._ntotact
for j
in range(len(p_dms))])
491 IF = basis.compute_IFsparse(dms, p_dms, p_geom).T
492 M2V, _ = basis.compute_btt(IF[:, :-2], IF[:, -2:].toarray())
493 M2V = M2V[:, list(range(p_controller.nmodes - 2)) + [-2, -1]]
495 rtc.d_control[i].init_modalOpti(p_controller.nmodes, p_controller.nrec, M2V,
496 p_controller.gmin, p_controller.gmax,
497 p_controller.ngain, 1. / ittime)
498 ol_slopes = modopti.open_loopSlp(tel, atmos, wfs, rtc, p_controller.nrec, i,
500 rtc.d_control[i].loadopen_loopSlp(ol_slopes)
501 rtc.d_control[i].modalControlOptimization()
503 cmats.cmat_init(i, rtc, p_controller, p_wfss, p_atmos, p_tel, p_dms,
504 nmodes=p_controller.nmodes)
506 rtc.d_control[i].set_gain(p_controller.gain)
508 sum([p_dms[j]._ntotact
for j
in range(len(p_dms))]), dtype=np.float32)
511 mgain[cc:cc + ndm._ntotact] = ndm.gain
513 rtc.d_control[i].set_modal_gains(mgain)
517 p_dms: list, p_wfss: list):
519 Initialize the CURED controller
521 i : (int) : controller index
522 rtc: (Rtc) : Rtc objet
523 p_controller: (Param_controller) : controller settings
524 p_dms: (list of Param_dms) : dms settings
525 p_wfss: (list of Param_wfs) : wfs settings
528 print(
"initializing cured controller")
529 if (scons.DmType.TT
in [p_dms[j].type
for j
in range(len(p_dms))]):
533 rtc.d_control[i].init_cured(p_wfss[0].nxsub, p_wfss[0]._isvalid,
534 p_controller.cured_ndivs, tt_flag)
535 rtc.d_control[i].set_gain(p_controller.gain)
539 p_geom: conf.Param_geom, p_dms: list, p_atmos: conf.Param_atmos,
540 p_tel: conf.Param_tel, rtc: Rtc, dms: Dms, wfs: Sensors,
543 Initialize the MV controller
546 i : (int) : controller index
547 p_controller: (Param_controller) : controller settings
548 p_wfss: (list of Param_wfs) : wfs settings
549 p_geom: (Param_geom) : geom settings
550 p_dms: (list of Param_dms) : dms settings
551 p_atmos: (Param_atmos) : atmos settings
552 p_tel: (Param_tel) : telescope settings
553 rtc: (Rtc) : Rtc objet
554 dms: (Dms) : Dms object
555 wfs: (Sensors) : Sensors object
556 atmos: (Atmos) : Atmos object
558 p_controller._imat = imats.imat_geom(wfs, dms, p_wfss, p_dms, p_controller)
560 rtc.d_control[i].set_imat(p_controller._imat)
561 rtc.d_control[i].set_gain(p_controller.gain)
562 size = sum([p_dms[j]._ntotact
for j
in range(len(p_dms))])
563 mgain = np.ones(size, dtype=np.float32)
564 rtc.d_control[i].set_modal_gains(mgain)
565 tomo.do_tomo_matrices(i, rtc, p_wfss, dms, atmos, wfs, p_controller, p_geom, p_dms,
567 cmats.cmat_init(i, rtc, p_controller, p_wfss, p_atmos, p_tel, p_dms)
573 Initialize the generic controller
576 i: (int): controller index
577 p_controller: (Param_controller): controller settings
578 p_dms: (list of Param_dm): dms settings
579 rtc: (Rtc): Rtc object
581 size = sum([p_dms[j]._ntotact
for j
in range(len(p_dms))])
582 decayFactor = np.ones(size, dtype=np.float32)
583 mgain = np.ones(size, dtype=np.float32) * p_controller.gain
584 matE = np.identity(size, dtype=np.float32)
585 cmat = np.zeros((size, p_controller.nslope), dtype=np.float32)
587 rtc.d_control[i].set_decayFactor(decayFactor)
588 rtc.d_control[i].set_modal_gains(mgain)
589 rtc.d_control[i].set_cmat(cmat)
590 rtc.d_control[i].set_matE(matE)