Merge "main" into "craig_factored"

tractor/approx3.c

@@ -97,23 +97,23 @@ static int c_gauss_2d_approx3(int x0, int x1, int y0, int y1,
             ampsum += amp[k];
 
         // We symmetrize the covariance matrix,
-        // so V,I just have three elements for each K: x**2, xy, y**2.
+        // so V,icov just have three elements for each K: x**2, xy, y**2.
         for (k=0; k<K; k++) {
-            // We also scale the I to make the Gaussian evaluation easier
+            // We also scale the icov to make the Gaussian evaluation easier
             double det;
             double isc;
             double* V = VV + 3*k;
-            double* I = II + 3*k;
+            double* icov = II + 3*k;
             V[0] =  var[k*D*D + 0];
             V[1] = (var[k*D*D + 1] + var[k*D*D + 2])*0.5;
             V[2] =  var[k*D*D + 3];
             det = V[0]*V[2] - V[1]*V[1];
             // we fold the -0.5 in the Gaussian exponent term in here...
             isc = -0.5 / det;
-            I[0] =  V[2] * isc;
+            icov[0] =  V[2] * isc;
             // we also fold in the 2*dx*dy term here
-            I[1] = -V[1] * isc * 2.0;
-            I[2] =  V[0] * isc;
+            icov[1] = -V[1] * isc * 2.0;
+            icov[2] =  V[0] * isc;
             scales[k] = amp[k] / sqrt(tpd * det);
             if (minval > 0.0) {
                 maxD[k] = log(minval * sqrt(tpd*det) / ampsum);

tractor/dense_optimizer.py

@@ -3,7 +3,7 @@
 from tractor.engine import logverb
 from tractor.constrained_optimizer import ConstrainedOptimizer
 
-from numpy.linalg import lstsq
+from numpy.linalg import lstsq, LinAlgError
 
 # from astrometry.util.plotutils import PlotSequence
 # ps = PlotSequence('opt')
@@ -156,7 +156,7 @@ def getUpdateDirection(self, tractor, allderivs, damp=0., priors=True,
             #print('Priors: pb', pb, 'mub', mub)
             for ri,ci,vi,bi in zip(rA, cA, vA, pb):
                 if scale_columns:
-                    colscales2[col] += np.dot(vi, vi)
+                    colscales2[ci] += np.dot(vi, vi)
                 for rij,vij,bij in zip(ri, vi, bi):
                     A[Npixels + rij, ci] = vij
                     B[Npixels + rij] += bij
@@ -186,7 +186,18 @@ def getUpdateDirection(self, tractor, allderivs, damp=0., priors=True,
             del chi
 
         # X, resids, rank, singular_vals
-        X,_,_,_ = lstsq(A, B, rcond=None)
+        try:
+            X,_,_,_ = lstsq(A, B, rcond=None)
+        except LinAlgError as e:
+            print('Exception in lstsq:', e)
+            from collections import Counter
+            import traceback
+            traceback.print_exc()
+            if scale_columns:
+                print('Column scales:', colscales)
+            print('A finite:', Counter(np.isfinite(A.ravel())))
+            print('B finite:', Counter(np.isfinite(B.ravel())))
+            return None
 
         if False:
             Aold = super(ConstrainedDenseOptimizer, self).getUpdateDirection(

tractor/emfit.i

@@ -396,7 +396,7 @@ static int em_fit_2d_reg(PyObject* po_img, int x0, int y0,
         for (k=0; k<K; k++) {
             // ASSUME ordering
             double* V = var + k*D*D;
-            double* I = ivar + k*D*D;
+            double* ivar_k = ivar + k*D*D;
             double det;
             //printf("var[%i]: %.3f,%.3f,%.3f,%.3f\n", k, V[0], V[1], V[2], V[3]);
             det = V[0]*V[3] - V[1]*V[2];
@@ -406,10 +406,10 @@ static int em_fit_2d_reg(PyObject* po_img, int x0, int y0,
                 result = -1;
                 goto cleanup;
             }
-            I[0] =  V[3] / det;
-            I[1] = -V[1] / det;
-            I[2] = -V[2] / det;
-            I[3] =  V[0] / det;
+            ivar_k[0] =  V[3] / det;
+            ivar_k[1] = -V[1] / det;
+            ivar_k[2] = -V[2] / det;
+            ivar_k[3] =  V[0] / det;
             scale[k] = amp[k] / sqrt(tpd * det);
         }
 

tractor/emfit2.c

@@ -146,7 +146,7 @@ static int em_fit_2d_reg2(PyObject* po_img, int x0, int y0,
         for (k=0; k<K; k++) {
             // ASSUME ordering
             double* V = var + k*D*D;
-            double* I = ivar + k*D*D;
+            double* ivar_k = ivar + k*D*D;
             double det;
             // symmetrize
             double V12 = (V[1] + V[2])/2.;
@@ -157,9 +157,9 @@ static int em_fit_2d_reg2(PyObject* po_img, int x0, int y0,
                 result = -1;
                 goto cleanup;
             }
-            I[0] =  V[3] / det;
-            I[1] = -2. * V12  / det;
-            I[3] =  V[0] / det;
+            ivar_k[0] =  V[3] / det;
+            ivar_k[1] = -2. * V12  / det;
+            ivar_k[3] =  V[0] / det;
             scale[k] = amp[k] / sqrt(tpd * det);
 
             // maxD: we want to allow each component (before

tractor/gauss_masked.c

@@ -67,34 +67,34 @@ static int c_gauss_2d_masked(int x0, int y0, int W, int H,
             ampsum += amp[k];
 
         // We symmetrize the covariance matrix,
-        // so V,I just have three elements for each K: x**2, xy, y**2.
+        // so V,icov just have three elements for each K: x**2, xy, y**2.
         for (k=0; k<K; k++) {
-            // We also scale the I to make the Gaussian evaluation easier
+            // We also scale the icov to make the Gaussian evaluation easier
             float det;
             float isc;
             float* V = VV + 3*k;
-            float* I = II + 3*k;
+            float* icov = II + 3*k;
             V[0] =  var[k*D*D + 0];
             V[1] = (var[k*D*D + 1] + var[k*D*D + 2])*0.5;
             V[2] =  var[k*D*D + 3];
             det = V[0]*V[2] - V[1]*V[1];
             // we fold the -0.5 in the Gaussian exponent term in here...
             isc = -0.5 / det;
-            I[0] =  V[2] * isc;
+            icov[0] =  V[2] * isc;
             // we also fold in the 2*dx*dy term here
-            I[1] = -V[1] * isc * 2.0;
-            I[2] =  V[0] * isc;
+            icov[1] = -V[1] * isc * 2.0;
+            icov[2] =  V[0] * isc;
             scales[k] = amp[k] / sqrt(tpd * det);
             maxD[k] = -30.;
 
             if (!(isfinite(V[0]) && isfinite(V[1]) && isfinite(V[2]) &&
-                  isfinite(I[0]) && isfinite(I[1]) && isfinite(I[2]) &&
+                  isfinite(icov[0]) && isfinite(icov[1]) && isfinite(icov[2]) &&
                   isfinite(scales[k]))) {
                 //printf("Warning: infinite variance or scale.  Zeroing.\n");
                 // large variance can cause this... set scale = 0.
                 scales[k] = 0.;
                 V[0] = V[1] = V[2] = 1.;
-                I[0] = I[1] = I[2] = -1.;
+                icov[0] = icov[1] = icov[2] = -1.;
                 det = 1.;
                 isc = 1.;
             }