diff --git a/aotools/turbulence/infinitephasescreen.py b/aotools/turbulence/infinitephasescreen.py index 53ae6df..a43a9e3 100644 --- a/aotools/turbulence/infinitephasescreen.py +++ b/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): """ diff --git a/test/test_infinitephasescreen.py b/test/test_infinitephasescreen.py index 21cdde7..d5f272d 100644 --- a/test/test_infinitephasescreen.py +++ b/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