COMPASS  5.0.0
End-to-end AO simulation tool using GPU acceleration
modalBasis.py
1 
37 from shesha.ao import basis
38 import shesha.util.utilities as util
39 import shesha.util.make_pupil as mkP
40 import shesha.constants as scons
41 import scipy.ndimage
42 from scipy.sparse.csr import csr_matrix
43 import numpy as np
44 
45 class ModalBasis(object):
46  """ This optimizer class handles all the modal basis and DM Influence functions
47  related operations.
48 
49  Attributes:
50  config : (config) : Configuration parameters module
51 
52  dms : (DmCompass) : DmCompass instance
53 
54  target : (TargetCompass) : TargetCompass instance
55 
56  slaved_actus : TODO : docstring
57 
58  selected_actus : TODO : docstring
59 
60  couples_actus : TODO : docstring
61 
62  index_under_spiders : TODO : docstring
63 
64  modal_basis : (np.ndarray) : Last modal basis computed
65 
66  projection_matrix : (np.ndarray) : Last projection_matrix computed
67  """
68  def __init__(self, config, dms, target):
69  """ Instantiate a ModalBasis object
70 
71  Parameters:
72  config : (config) : Configuration parameters module
73 
74  dms : (DmCompass) : DmCompass instance
75 
76  target : (TargetCompass) : TargetCompass instance
77  """
78  self.config = config
79  self.dms = dms
80  self.target = target
81  self.slaved_actus = None
82  self.selected_actus = None
83  self.couples_actus = None
84  self.index_under_spiders = None
85  self.modal_basis = None
86  self.projection_matrix = None
87 
88  def compute_influ_basis(self, dm_index: int) -> csr_matrix:
89  """ Computes and return the influence function phase basis of the specified DM
90  as a sparse matrix
91 
92  Parameters:
93  dm_index : (int) : Index of the DM
94 
95  Return:
96  influ_sparse : (csr_matrix) : influence function phases
97  """
98  return basis.compute_dm_basis(self.dms.dms.d_dms[dm_index],
99  self.config.p_dms[dm_index],
100  self.config.p_geom)
101 
102  def compute_modes_to_volts_basis(self, modal_basis_type: str, *, merged: bool = False,
103  nbpairs: int = None, return_delta: bool = False) -> np.ndarray:
104  """ Computes a given modal basis ("KL2V", "Btt", "Btt_petal") and return the 2 transfer matrices
105 
106  Parameters:
107  modal_basis_type : (str) : modal basis to compute ("KL2V", "Btt", "Btt_petal")
108 
109  merged : (bool, optional) :
110 
111  nbpairs : (int, optional) :
112 
113  Return:
114  modal_basis : (np.ndarray) : modes to volts matrix
115 
116  projection_matrix : (np.ndarray) : volts to modes matrix (None if "KL")
117  """
118  if (modal_basis_type == "KL2V"):
119  print("Computing KL2V basis...")
120  self.modal_basis = basis.compute_KL2V(
121  self.config.p_controllers[0], self.dms.dms,
122  self.config.p_dms, self.config.p_geom,
123  self.config.p_atmos, self.config.p_tel)
124  fnz = util.first_non_zero(self.modal_basis, axis=0)
125  # Computing the sign of the first non zero element
126  #sig = np.sign(modal_basis[[fnz, np.arange(modal_basis.shape[1])]])
127  sig = np.sign(self.modal_basis[tuple([
128  fnz, np.arange(self.modal_basis.shape[1])
129  ])]) # pour remove le future warning!
130  self.modal_basis *= sig[None, :]
131  projection_matrix = None
132  elif (modal_basis_type == "Btt"):
133  print("Computing Btt basis...")
135  merged=merged, nbpairs=nbpairs,
136  return_delta=return_delta)
137  fnz = util.first_non_zero(self.modal_basis, axis=0)
138  # Computing the sign of the first non zero element
139  #sig = np.sign(modal_basis[[fnz, np.arange(modal_basis.shape[1])]])
140  sig = np.sign(self.modal_basis[tuple([
141  fnz, np.arange(self.modal_basis.shape[1])
142  ])]) # pour remove le future warning!
143  self.modal_basis *= sig[None, :]
144  elif (modal_basis_type == "Btt_petal"):
145  print("Computing Btt with a Petal basis...")
147  else:
148  raise ArgumentError("Unsupported modal basis")
149 
150  return self.modal_basis, self.projection_matrix
151 
152  def compute_btt_basis(self, *, merged: bool = False, nbpairs: int = None,
153  return_delta: bool = False) -> np.ndarray:
154  """ Computes the so-called Btt modal basis. The <merged> flag allows merto merge
155  2x2 the actuators influence functions for actuators on each side of the spider (ELT case)
156 
157  Parameters:
158  merged : (bool, optional) : If True, merge 2x2 the actuators influence functions for
159  actuators on each side of the spider (ELT case). Default
160  is False
161 
162  nbpairs : (int, optional) : Default is None. TODO : description
163 
164  return_delta : (bool, optional) : If False (default), the function returns
165  Btt (modes to volts matrix),
166  and P (volts to mode matrix).
167  If True, returns delta = IF.T.dot(IF) / N
168  instead of P
169 
170  Return:
171  Btt : (np.ndarray) : Btt modes to volts matrix
172 
173  projection_matrix : (np.ndarray) : volts to Btt modes matrix
174  """
175  dms_basis = basis.compute_IFsparse(self.dms.dms, self.config.p_dms, self.config.p_geom)
176  influ_basis = dms_basis[:-2,:]
177  tt_basis = dms_basis[-2:,:].toarray()
178  if (merged):
179  couples_actus, index_under_spiders = self.compute_merged_influ(0,
180  nbpairs=nbpairs)
181  influ_basis2 = influ_basis.copy()
182  index_remove = index_under_spiders.copy()
183  index_remove += list(couples_actus[:, 1])
184  print("Pairing Actuators...")
185  for i in tqdm(range(couples_actus.shape[0])):
186  influ_basis2[couples_actus[i, 0], :] += influ_basis2[
187  couples_actus[i, 1], :]
188  print("Pairing Done")
189  boolarray = np.zeros(influ_basis2.shape[0], dtype=np.bool)
190  boolarray[index_remove] = True
191  self.slaved_actus = boolarray
192  self.selected_actus = ~boolarray
193  self.couples_actus = couples_actus
194  self.index_under_spiders = index_under_spiders
195  influ_basis2 = influ_basis2[~boolarray, :]
196  influ_basis = influ_basis2
197 
198  self.btt, self.projection_matrix = basis.compute_btt(influ_basis.T, tt_basis.T, return_delta=return_delta)
199 
200  if (merged):
201  btt2 = np.zeros((len(boolarray) + 2, self.btt.shape[1]))
202  btt2[np.r_[~boolarray, True, True], :] = self.btt
203  btt2[couples_actus[:, 1], :] = btt2[couples_actus[:, 0], :]
204 
205  P2 = np.zeros((self.btt.shape[1], len(boolarray) + 2))
206  P2[:, np.r_[~boolarray, True, True]] = self.projection_matrix
207  P2[:, couples_actus[:, 1]] = P2[:, couples_actus[:, 0]]
208  self.btt = btt2
209  self.projection_matrix = P2
210 
211  return self.btt, self.projection_matrix
212 
213  def compute_merged_influ(self, dm_index : int, *, nbpairs: int = None) -> np.ndarray:
214  """ Used to compute merged IF from each side of the spider
215  for an ELT case (Petalling Effect)
216 
217  Parameters:
218  dm_index : (int) : DM index
219 
220  nbpairs : (int, optional) : Default is None. TODO : description
221 
222  Return:
223  pairs : (np.ndarray) : TODO description
224 
225  discard : (list) : TODO description
226  """
227  p_geom = self.config.p_geom
228 
229 
230  cent = p_geom.pupdiam / 2. + 0.5
231  p_tel = self.config.p_tel
232  p_tel.t_spiders = 0.51
233  spup = mkP.make_pupil(p_geom.pupdiam, p_geom.pupdiam, p_tel, cent,
234  cent).astype(np.float32).T
235 
236  p_tel.t_spiders = 0.
237  spup2 = mkP.make_pupil(p_geom.pupdiam, p_geom.pupdiam, p_tel, cent,
238  cent).astype(np.float32).T
239 
240  spiders = spup2 - spup
241 
242  (spidersID, k) = scipy.ndimage.label(spiders)
243  spidersi = util.pad_array(spidersID, p_geom.ssize).astype(np.float32)
244  px_list_spider = [np.where(spidersi == i) for i in range(1, k + 1)]
245 
246  # DM positions in iPupil:
247  dm_posx = self.config.p_dms[dm_index]._xpos - 0.5
248  dm_posy = self.config.p_dms[dm_index]._ypos - 0.5
249  dm_pos_mat = np.c_[dm_posx, dm_posy].T # one actu per column
250 
251  pitch = self.config.p_dms[dm_index]._pitch
252  discard = np.zeros(len(dm_posx), dtype=np.bool)
253  pairs = []
254 
255  # For each of the k pieces of the spider
256  for k, px_list in enumerate(px_list_spider):
257  pts = np.c_[px_list[1],
258  px_list[0]] # x,y coord of pixels of the spider piece
259  # line_eq = [a, b]
260  # Which minimizes leqst squares of aa*x + bb*y = 1
261  line_eq = np.linalg.pinv(pts).dot(np.ones(pts.shape[0]))
262  aa, bb = line_eq[0], line_eq[1]
263 
264  # Find any point of the fitted line.
265  # For simplicity, the intercept with one of the axes x = 0 / y = 0
266  if np.abs(bb) < np.abs(aa): # near vertical
267  one_point = np.array([1 / aa, 0.])
268  else: # otherwise
269  one_point = np.array([0., 1 / bb])
270 
271  # Rotation that aligns the spider piece to the horizontal
272  rotation = np.array([[-bb, aa], [-aa, -bb]]) / (aa**2 + bb**2)**.5
273 
274  # Rotated the spider mask
275  rotated_px = rotation.dot(pts.T - one_point[:, None])
276  # Min and max coordinates along the spider length - to filter actuators that are on
277  # 'This' side of the pupil and not the other side
278  min_u, max_u = rotated_px[0].min() - 5. * pitch, rotated_px[0].max(
279  ) + 5. * pitch
280 
281  # Rotate the actuators
282  rotated_actus = rotation.dot(dm_pos_mat - one_point[:, None])
283  sel_good_side = (rotated_actus[0] > min_u) & (rotated_actus[0] < max_u)
284  threshold = 0.05
285  # Actuators below this piece of spider
286  sel_discard = (np.abs(rotated_actus[1]) < threshold * pitch) & sel_good_side
287  discard |= sel_discard
288 
289  # Actuator 'near' this piece of spider
290  sel_pairable = (np.abs(rotated_actus[1]) > threshold * pitch) & \
291  (np.abs(rotated_actus[1]) < 1. * pitch) & \
292  sel_good_side
293 
294  pairable_index = np.where(sel_pairable)[0] # Indices of these actuators
295  u_coord = rotated_actus[
296  0, sel_pairable] # Their linear coord along the spider major axis
297 
298  order = np.sort(u_coord) # Sort by linear coordinate
299  order_index = pairable_index[np.argsort(
300  u_coord)] # And keep track of original indexes
301 
302  # i = 0
303  # while i < len(order) - 1:
304  if (nbpairs is None):
305  i = 0
306  ii = len(order) - 1
307  else:
308  i = len(order) // 2 - nbpairs
309  ii = len(order) // 2 + nbpairs
310  while (i < ii):
311  # Check if next actu in sorted order is very close
312  # Some lonely actuators may be hanging in this list
313  if np.abs(order[i] - order[i + 1]) < .2 * pitch:
314  pairs += [(order_index[i], order_indx[i + 1])]
315  i += 2
316  else:
317  i += 1
318  print('To discard: %u actu' % np.sum(discard))
319  print('%u pairs to slave' % len(pairs))
320  if np.sum(discard) == 0:
321  discard = []
322  else:
323  list(np.where(discard)[0])
324  return np.asarray(pairs), list(np.where(discard)[0])
325 
326  def compute_btt_petal(self) -> np.ndarray:
327  """ Computes a Btt modal basis with Pistons filtered
328 
329  Return:
330  Btt : (np.ndarray) : Btt modes to volts matrix
331 
332  P : (np.ndarray) : volts to Btt modes matrix
333  """
334  pzt_index = np.where([d.type is scons.DmType.PZT for d in self.config.p_dms])[0][0]
335  influ_pzt = self.compute_influ_basis(pzt_index)
336  petal_dm_index = np.where([
337  d.influ_type is scons.InfluType.PETAL for d in self.config.p_dms
338  ])[0][0]
339  influ_petal = self.compute_influ_basis(petal_dm_index)
340  tt_index = np.where([d.type is scons.DmType.TT for d in self.config.p_dms])[0][0]
341  influ_tt = self.compute_influ_basis(tt_index).toarray()
342 
343  self.modal_basis, self.projection_matrix = basis.compute_btt(influ_pzt.T, influ_tt.T, influ_petal=influ_petal)
344  return self.modal_basis, self.projection_matrix
345 
346  def compute_phase_to_modes(self, modal_basis: np.ndarray) -> np.ndarray:
347  """ Return the phase to modes matrix by using the given modal basis
348 
349  Parameters:
350  modal_basis : (np.ndarray) : Modal basis matrix
351 
352  Return:
353  phase_to_modes : (np.ndarray) : phase to modes matrix
354  """
355  nbmode = modal_basis.shape[1]
356  phase = self.target.get_tar_phase(0)
357  phase_to_modes = np.zeros((nbmode, phase.shape[0], phase.shape[1]))
358  S = np.sum(self.config.p_geom._spupil)
359  for i in range(nbmode):
360  self.dms.set_command((modal_basis[:, i]).copy())
361  # self.next(see_atmos=False)
362  self.target.raytrace(0, dms=self.dms, ncpa=False, reset=True)
363  phase = self.target.get_tar_phase(0, pupil=True)
364  # Normalisation pour les unites rms en microns !!!
365  norm = np.sqrt(np.sum((phase)**2) / S)
366  phase_to_modes[i] = phase / norm
367  return phase_to_modes
shesha.supervisor.optimizers.modalBasis.ModalBasis.btt
btt
Definition: modalBasis.py:233
shesha.supervisor.optimizers.modalBasis.ModalBasis.config
config
Definition: modalBasis.py:104
shesha.supervisor.optimizers.modalBasis.ModalBasis.couples_actus
couples_actus
Definition: modalBasis.py:109
shesha.supervisor.optimizers.modalBasis.ModalBasis.compute_merged_influ
np.ndarray compute_merged_influ(self, int dm_index, *int nbpairs=None)
Used to compute merged IF from each side of the spider for an ELT case (Petalling Effect)
Definition: modalBasis.py:252
shesha.ao
Python package for AO operations on COMPASS simulation.
Definition: shesha/shesha/ao/__init__.py:1
shesha.supervisor.optimizers.modalBasis.ModalBasis.modal_basis
modal_basis
Definition: modalBasis.py:111
shesha.supervisor.optimizers.modalBasis.ModalBasis.slaved_actus
slaved_actus
Definition: modalBasis.py:107
shesha.supervisor.optimizers.modalBasis.ModalBasis.index_under_spiders
index_under_spiders
Definition: modalBasis.py:110
shesha.supervisor.optimizers.modalBasis.ModalBasis.__init__
def __init__(self, config, dms, target)
Definition: modalBasis.py:103
shesha.util.utilities
Basic utilities function.
Definition: utilities.py:1
shesha.constants
Numerical constants for shesha and config enumerations for safe-typing.
Definition: constants.py:1
shesha.util.make_pupil
Pupil creation functions.
Definition: make_pupil.py:1
shesha.supervisor.optimizers.modalBasis.ModalBasis.compute_phase_to_modes
np.ndarray compute_phase_to_modes(self, np.ndarray modal_basis)
Return the phase to modes matrix by using the given modal basis.
Definition: modalBasis.py:380
shesha.supervisor.optimizers.modalBasis.ModalBasis.dms
dms
Definition: modalBasis.py:105
shesha.supervisor.optimizers.modalBasis.ModalBasis.compute_modes_to_volts_basis
np.ndarray compute_modes_to_volts_basis(self, str modal_basis_type, *bool merged=False, int nbpairs=None, bool return_delta=False)
Computes a given modal basis ("KL2V", "Btt", "Btt_petal") and return the 2 transfer matrices.
Definition: modalBasis.py:142
shesha.supervisor.optimizers.modalBasis.ModalBasis.target
target
Definition: modalBasis.py:106
shesha.supervisor.optimizers.modalBasis.ModalBasis.projection_matrix
projection_matrix
Instantiate a ModalBasis object.
Definition: modalBasis.py:112
shesha.supervisor.optimizers.modalBasis.ModalBasis.compute_btt_petal
np.ndarray compute_btt_petal(self)
Computes a Btt modal basis with Pistons filtered.
Definition: modalBasis.py:359
shesha.supervisor.optimizers.modalBasis.ModalBasis.selected_actus
selected_actus
Definition: modalBasis.py:108
shesha.supervisor.optimizers.modalBasis.ModalBasis.compute_influ_basis
csr_matrix compute_influ_basis(self, int dm_index)
Computes and return the influence function phase basis of the specified DM as a sparse matrix.
Definition: modalBasis.py:123
shesha.supervisor.optimizers.modalBasis.ModalBasis.compute_btt_basis
np.ndarray compute_btt_basis(self, *bool merged=False, int nbpairs=None, bool return_delta=False)
Computes the so-called Btt modal basis.
Definition: modalBasis.py:199
shesha.supervisor.optimizers.modalBasis.ModalBasis
This optimizer class handles all the modal basis and DM Influence functions related operations.
Definition: modalBasis.py:50