A few unrelated bugfixes

crowdsource/crowdsource_base.py

@@ -13,6 +13,7 @@
 entirely simultaneously.
 """
 
+import time
 import sys
 import os
 import numpy
@@ -295,10 +296,13 @@ def fit_once(im, x, y, psfs, weight=None,
     repeat = 1 if not psfderiv else 3
     nskypar = nskyx * nskyy
     npixim = im.shape[0]*im.shape[1]
-    xloc = numpy.zeros(repeat*numpy.sum(sz*sz).astype('i4') +
-                       nskypar*npixim, dtype='i4')
-    # yloc = numpy.zeros(len(xloc), dtype='i4')
-    # no longer need yloc; csc entries are built directly.
+    zsz = (repeat*numpy.sum(sz*sz) + nskypar*npixim).astype('i4')
+    if zsz >= 2**32:
+        raise ValueError(
+            'Number of pixels being fit is too large (>2**32); '
+            'failing early.  This usually indicates a problem with '
+            'the choice of PSF size & too many sources.')
+    xloc = numpy.zeros(zsz, dtype='i4')
     values = numpy.zeros(len(xloc), dtype='f4')
     colnorm = numpy.zeros(len(x)*repeat+nskypar, dtype='f4')
     first = 0
@@ -569,33 +573,38 @@ def sky_im(im, weight=None, npix=20, order=1):
     return bg
 
 
-def get_sizes(x, y, imbs, weight=None, blist=None):
+def get_sizes(x, y, imbs, weight=None, blist=None, blistsz=299,
+              cutofflist=None):
+    if cutofflist is None:
+        cutofflist = [
+            (-numpy.inf, 19), (1000, 59), (20000, 149)]
     x = numpy.round(x).astype('i4')
     y = numpy.round(y).astype('i4')
     peakbright = imbs[x, y]
-    sz = numpy.zeros(len(x), dtype='i4')
-    cutoff = 1000
-    sz[peakbright > cutoff] = 59
-    sz[peakbright <= cutoff] = 19  # for the moment...
-    # for very bright things, use a bigger PSF
-    # but if there are too many of these, don't bother.
-    cutoff2 = 20000
-    if ((numpy.sum(peakbright > cutoff2) < numpy.sum(peakbright > cutoff)/2)
-            or (numpy.sum(peakbright > cutoff) < 100)):
-        sz[peakbright > cutoff2] = 149
-    else:
-        print('Too many bright sources, using smaller PSF stamp size...')
 
     if weight is not None:
-        sz[weight[x, y] == 0] = 149  # saturated/off edge sources get big PSF
+        # treat saturated / off edge sources as very bright.
+        peakbright[weight[x, y] == 0] = cutofflist[-1][0] + 1
+
+    sz = numpy.zeros(len(x), dtype='i4')
+    nbright = list()
+    for cutoff, tsz in cutofflist:
+        m = peakbright > cutoff
+        sz[m] = tsz
+        nbright.append(numpy.sum(m))
+
+    if ((len(nbright) > 2) and (nbright[-1] > 100) and
+            (nbright[-1] > nbright[-2] / 2)):
+        print('Too many bright sources, using smaller PSF stamp size...')
+        sz[peakbright > cutofflist[-2][0]] = cutofflist[-2][1]
 
     # sources near listed sources get very big PSF
     if blist is not None and len(x) > 0:
         for xb, yb in zip(blist[0], blist[1]):
             dist2 = (x-xb)**2 + (y-yb)**2
             indclose = numpy.argmin(dist2)
             if dist2[indclose] < 5**2:
-                sz[indclose] = 299
+                sz[indclose] = blistsz
     return sz
 
 
@@ -809,8 +818,16 @@ def fit_im(im, psf, weight=None, dq=None, psfderiv=True,
         psfs = [numpy.zeros((len(xa), minsz, minsz), dtype='f4')
                 for i in range(repeat)]
         flux = numpy.zeros(len(xa)*repeat, dtype='f4')
+        if verbose:
+            subreg_iter = 0
+            t0 = time.time()
+            print("Starting subregion iterations")
         for (bdxf, bdxl, bdxaf, bdxal, bdyf, bdyl, bdyaf, bdyal) in (
                 subregions(im.shape, ntilex, ntiley)):
+            if verbose:
+                print(f"Subregion iteration {subreg_iter} starting; "
+                      f"dt={time.time()-t0}", flush=True)
+                subreg_iter += 1
             mbda = in_bounds(xa, ya, [bdxaf-0.5, bdxal-0.5],
                              [bdyaf-0.5, bdyal-0.5])
             mbd = in_bounds(xa, ya, [bdxf-0.5, bdxl-0.5],
@@ -831,6 +848,8 @@ def fit_im(im, psf, weight=None, dq=None, psfderiv=True,
                 im[sall]-sky[sall], xa[mbda]-bdxaf, ya[mbda]-bdyaf, psfsbda,
                 psfderiv=tpsfderiv, weight=weightbda, guess=guessmbda,
                 nskyx=nskyx, nskyy=nskyy)
+            if numpy.all(numpy.isnan(tmodel)):
+                raise ValueError("Model is all NaNs")
             model[spri] = tmodel[sfit]
             msky[spri] = tmsky[sfit]
             ind = numpy.flatnonzero(mbd)
@@ -1029,7 +1048,7 @@ def compute_iso_fit(impsfstack, psfstack, weightstack, apcor, psfderiv):
 
 
 def sky_model_basis(i, j, nskyx, nskyy, nx, ny):
-    import basisspline
+    from crowdsource import basisspline
     if (nskyx < 3) or (nskyy < 3):
         raise ValueError('Invalid sky model.')
     expandx = (nskyx-1.)/(3-1)

setup.py

@@ -1,14 +1,14 @@
 from setuptools import setup
 
 setup(
-    name = 'crowdsourcephoto',
+    name = 'crowdsource',
     packages = ['crowdsource'],
-    version = '0.5.6',
+    version = '0.5.8',
     description = 'Crowded field photometry pipeline',
     author = 'Andrew Saydjari',
     author_email = 'aksaydjari@gmail.com',
     url = 'https://github.com/schlafly/crowdsource',
-    download_url = 'https://github.com/schlafly/crowdsource/archive/refs/tags/v0.5.6.tar.gz',
+    download_url = 'https://github.com/schlafly/crowdsource/archive/refs/tags/v0.5.8.tar.gz',
     license = 'MIT',
     install_requires=[
         'astropy',