fall back to least squares solver when Cholesky fails for infinite phasescreen

Author: bernadett92

aotools/turbulence/infinitephasescreen.py

@@ -143,16 +143,14 @@ def makeAMatrix(self):
         Calculates the "A" matrix, that uses the existing data to find a new 
         component of the new phase vector.
         """
-        # Cholsky solve can fail - if so do brute force inversion
         try:
             cf = linalg.cho_factor(self.cov_mat_zz)
             inv_cov_zz = linalg.cho_solve(cf, numpy.identity(self.cov_mat_zz.shape[0]))
+            self.A_mat = self.cov_mat_xz.dot(inv_cov_zz)
         except linalg.LinAlgError:
-            # print("Cholesky solve failed. Performing SVD inversion...")
-            # inv_cov_zz = numpy.linalg.pinv(self.cov_mat_zz)
-            raise linalg.LinAlgError("Could not invert Covariance Matrix to for A and B Matrices. Try with a larger pixel scale or smaller L0")
-
-        self.A_mat = self.cov_mat_xz.dot(inv_cov_zz)
+            print("Cholesky solve failed. Performing least squares inversion...")
+            inv_cov_zz = numpy.linalg.lstsq(self.cov_mat_zz, numpy.identity(self.cov_mat_zz.shape[0]), rcond=1e-8)
+            self.A_mat = self.cov_mat_xz.dot(inv_cov_zz[0])
 
     def makeBMatrix(self):
         """

test/test_infinitephasescreen.py

@@ -1,4 +1,5 @@
 from aotools.turbulence import infinitephasescreen
+import numpy as np
 
 def testVKInitScreen():
 
@@ -9,6 +10,15 @@ def testVKAddRow():
     scrn = infinitephasescreen.PhaseScreenVonKarman(128, 4./64, 0.2, 50, n_columns=4)
     scrn.add_row()
 
+def testLeastSquaresSolver():
+    airmass = 1.0 / np.cos(30.0 / 180. * np.pi)
+    r0 = 0.9759 * 0.5 / (0.7 * 4.848) * airmass ** (-3. / 5.)
+    r0wavelength = r0 * (500 / 500.0) ** (6. / 5.)
+    screen = infinitephasescreen.PhaseScreenVonKarman(120, 1./120, r0wavelength, 50, 1, n_columns=2)
+    mean_scrn = np.sqrt(np.mean(screen.scrn**2))
+    for i in range(500):
+        screen.add_row()
+    assert(0.5 <= mean_scrn/np.sqrt(np.mean(screen.scrn**2)) <= 1.5)
 
 
 # Test of Kolmogoroc screen