Merge pull request #31 from andrewpaulreeves/slopecovariance_bugfix

Slopecovariance bugfix

Author: matthewtownson

aotools/turbulence/slopecovariance.py

@@ -139,7 +139,7 @@ def make_covariance_matrix(self):
 
             self._make_covariance_matrix_mp(self.threads)
 
-        self.covariance_matrix = mirror_covariance_matrix(self.covariance_matrix, self.n_subaps)
+        self.covariance_matrix = mirror_covariance_matrix(self.covariance_matrix)
 
         return self.covariance_matrix
 
@@ -184,7 +184,7 @@ def _make_covariance_matrix(self):
                     self.covariance_matrix[
                             cov_mat_coord_x1: cov_mat_coord_x2,
                             cov_mat_coord_y1 + self.n_subaps[wfs_j]: cov_mat_coord_y2 + self.n_subaps[wfs_j]
-                            ] += cov_xy.T * r0_scale
+                            ] += numpy.fliplr(numpy.flipud(cov_xy)) * r0_scale
                     self.covariance_matrix[
                             cov_mat_coord_x1 + self.n_subaps[wfs_i]: cov_mat_coord_x2 + self.n_subaps[wfs_i],
                             cov_mat_coord_y1 + self.n_subaps[wfs_j]: cov_mat_coord_y2 + self.n_subaps[wfs_j]
@@ -239,7 +239,7 @@ def _make_covariance_matrix_mp(self, threads):
                     self.covariance_matrix[
                             cov_mat_coord_x1: cov_mat_coord_x2,
                             cov_mat_coord_y1 + self.n_subaps[wfs_j]: cov_mat_coord_y2 + self.n_subaps[wfs_j]
-                            ] += cov_xy.T * r0_scale
+                            ] += numpy.fliplr(numpy.flipud(cov_xy)) * r0_scale
                     self.covariance_matrix[
                             cov_mat_coord_x1 + self.n_subaps[wfs_i]: cov_mat_coord_x2 + self.n_subaps[wfs_i],
                             cov_mat_coord_y1 + self.n_subaps[wfs_j]: cov_mat_coord_y2 + self.n_subaps[wfs_j]
@@ -462,38 +462,16 @@ def calculate_structure_function(phase, nbOfPoint=None, step=None):
     return sf_x
 
 
-def mirror_covariance_matrix(cov_mat, n_subaps):
+def mirror_covariance_matrix(cov_mat):
     """
     Mirrors a covariance matrix around the axis of the diagonal.
 
     Parameters:
         cov_mat (ndarray): The covariance matrix to mirror
         n_subaps (ndarray): Number of sub-aperture in each WFS
     """
-    total_slopes = cov_mat.shape[0]
-    n_wfs = n_subaps.shape[0]
-
-    n1 = 0
-    for n in range(n_wfs):
-        m1 = 0
-        for m in range(n + 1):
-            if n != m:
-                n2 = n1 + 2 * n_subaps[n]
-                m2 = m1 + 2 * n_subaps[m]
-
-                nn1 = total_slopes - 2 * n_subaps[n] - n1
-                nn2 = nn1 + 2 * n_subaps[n]
-
-                mm1 = total_slopes - 2 * n_subaps[m] - m1
-                mm2 = mm1 + 2 * n_subaps[m]
-
-                cov_mat[nn1: nn2, mm1: mm2] = (
-                    numpy.swapaxes(cov_mat[n1: n2, m1: m2], 1, 0)
-                )
-
-                m1 += 2 * n_subaps[m]
-        n1 += 2 * n_subaps[n]
-    return cov_mat
+
+    return numpy.bitwise_or(cov_mat.view("int32"), cov_mat.T.view("int32")).view("float32")
 
 def create_tomographic_covariance_reconstructor(covariance_matrix, n_onaxis_subaps, svd_conditioning=0):
     """

test/test_slopecovariance.py

@@ -52,6 +52,7 @@ def test_slopecovmat_makecovmat():
     asterism_radius = 10
 
     subap_diameters = [telescope_diameter / nx_subaps] * n_wfs
+
     pupil_masks = [aotools.circle(nx_subaps / 2., nx_subaps)] * n_wfs
     gs_altitudes = [90000] * n_wfs
     gs_positions = [
@@ -73,6 +74,41 @@ def test_slopecovmat_makecovmat():
     covariance_matrix = cov_mat.make_covariance_matrix()
 
 
+def test_slopecovmat_makecovmat_uneven():
+    threads = 1
+
+    n_wfs = 3
+    telescope_diameter = 8.
+    nx_subaps = 10
+
+    n_layers = 3
+    layer_altitudes = numpy.linspace(0, 20000, n_layers)
+    layer_r0s = [1] * n_layers
+    layer_L0s = [25.] * n_layers
+
+    asterism_radius = 10
+    # What if all WFSs don't have the same number of subaps?
+    pupil_masks = [aotools.circle(4.5, nx_subaps), aotools.circle(nx_subaps / 2., nx_subaps), aotools.circle(nx_subaps / 2., nx_subaps)]
+    subap_diameters = [telescope_diameter / nx_subaps] * n_wfs
+    gs_altitudes = [90000] * n_wfs
+    gs_positions = [
+        [asterism_radius, 0],
+        [numpy.sin(numpy.pi / 3.) * asterism_radius, numpy.cos(numpy.pi / 3.) * asterism_radius],
+        [numpy.sin(numpy.pi / 3.) * asterism_radius, -numpy.cos(numpy.pi / 3.) * asterism_radius],
+        [-numpy.sin(numpy.pi / 3.) * asterism_radius, numpy.cos(numpy.pi / 3.) * asterism_radius],
+        [-numpy.sin(numpy.pi / 3.) * asterism_radius, -numpy.cos(numpy.pi / 3.) * asterism_radius],
+        [-asterism_radius, 0]]
+    wfs_magnifications = [1.] * n_wfs
+    pupil_offsets = [[0, 0]] * n_wfs
+    wfs_rotations = [0] * n_wfs
+    wfs_wavelengths = [550e-9] * n_wfs
+
+    cov_mat = aotools.CovarianceMatrix(n_wfs, pupil_masks, telescope_diameter, subap_diameters, gs_altitudes, gs_positions,
+                                    wfs_wavelengths,
+                                    n_layers, layer_altitudes, layer_r0s, layer_L0s, threads)
+
+    covariance_matrix = cov_mat.make_covariance_matrix()
+
 def test_covtomorecon():
     threads = 1