DM-42285: Update to support numpy 1.25+ and various optimizations.

fgcm/fgcmChisq.py

@@ -8,6 +8,7 @@
 from .fgcmUtilities import MaxFitIterations
 from .fgcmUtilities import Cheb2dField
 from .fgcmUtilities import objFlagDict
+from .fgcmUtilities import histogram_rev_sorted
 
 from .fgcmNumbaUtilities import numba_test, add_at_1d, add_at_2d, add_at_3d
 
@@ -606,7 +607,7 @@ def _magWorker(self, goodStarsAndObs):
                 obsYGO = snmm.getArray(self.fgcmStars.obsYHandle)[goodObs]
                 expCcdHash = (obsExpIndexGO[ok] * (self.fgcmPars.nCCD + 1) +
                               obsCCDIndexGO[ok])
-                h, rev = esutil.stat.histogram(expCcdHash, rev=True)
+                h, rev = histogram_rev_sorted(expCcdHash)
                 use, = np.where(h > 0)
                 for i in use:
                     i1a = rev[rev[i]: rev[i + 1]]

fgcm/fgcmChromaticity.py

@@ -5,6 +5,7 @@
 import time
 import matplotlib.pyplot as plt
 import scipy.optimize as optimize
+from .fgcmUtilities import histogram_rev_sorted
 
 from .sharedNumpyMemManager import SharedNumpyMemManager as snmm
 
@@ -354,7 +355,7 @@ def computeCCDChromaticity(self):
         ccdFilterHash = (obsLUTFilterIndexGO.astype(np.int64)*(self.fgcmPars.nCCD + 1) +
                          obsCCDIndexGO.astype(np.int64))
 
-        h, rev = esutil.stat.histogram(ccdFilterHash, rev=True)
+        h, rev = histogram_rev_sorted(ccdFilterHash)
 
         # Make a simple cut here.
         use, = np.where(h >= 10)

fgcm/fgcmConnectivity.py

@@ -4,6 +4,7 @@
 import matplotlib.pyplot as plt
 
 from .fgcmUtilities import objFlagDict
+from .fgcmUtilities import histogram_rev_sorted
 from .sharedNumpyMemManager import SharedNumpyMemManager as snmm
 
 class FgcmConnectivity(object):
@@ -76,8 +77,11 @@ def plotConnectivity(self):
                                 ((objFlag[obsObjIDIndex] & mask) == 0))
 
             # Split into nights with a histogram
-            h, rev = esutil.stat.histogram(self.fgcmPars.expNightIndex[obsExpIndex[goodObs]],
-                                           min=0, max=self.fgcmPars.nCampaignNights, rev=True)
+            h, rev = histogram_rev_sorted(
+                self.fgcmPars.expNightIndex[obsExpIndex[goodObs]],
+                min=0,
+                max=self.fgcmPars.nCampaignNights,
+            )
 
             done = False
             while not done:

fgcm/fgcmDeltaAper.py

@@ -8,6 +8,7 @@
 import matplotlib.pyplot as plt
 
 from .fgcmUtilities import dataBinner
+from .fgcmUtilities import histogram_rev_sorted
 
 import multiprocessing
 
@@ -89,7 +90,7 @@ def computeDeltaAperExposures(self):
         self.fgcmPars.compMedDeltaAper[:] = self.illegalValue
         self.fgcmPars.compEpsilon[:] = self.illegalValue
 
-        h, rev = esutil.stat.histogram(obsExpIndex[goodObs], rev=True, min=0)
+        h, rev = histogram_rev_sorted(obsExpIndex[goodObs], min=0)
         expIndices, = np.where(h >= self.minStarPerExp)
 
         for expIndex in expIndices:
@@ -248,7 +249,11 @@ def computeEpsilonMap(self):
         # Do the mapping
         self.fgcmLog.info("Computing delta-aper epsilon spatial map at nside %d" % self.deltaAperFitSpatialNside)
         ipring = hpg.angle_to_pixel(self.deltaAperFitSpatialNside, objRA, objDec, nest=False)
-        h, rev = esutil.stat.histogram(ipring, min=0, max=hpg.nside_to_npixel(self.deltaAperFitSpatialNside) - 1, rev=True)
+        h, rev = histogram_rev_sorted(
+            ipring,
+            min=0,
+            max=hpg.nside_to_npixel(self.deltaAperFitSpatialNside) - 1,
+        )
 
         offsetMap = np.zeros(h.size, dtype=[('nstar_fit', 'i4', (self.fgcmStars.nBands, )),
                                             ('epsilon', 'f4', (self.fgcmStars.nBands, ))])
@@ -365,7 +370,7 @@ def computeEpsilonPerCcd(self):
 
         filterCcdHash = ccdIndexGO*(self.fgcmPars.nLUTFilter + 1) + lutFilterIndexGO
 
-        h, rev = esutil.stat.histogram(filterCcdHash, rev=True)
+        h, rev = histogram_rev_sorted(filterCcdHash)
 
         # Arbitrary minimum number here
         gdHash, = np.where(h > 10)
@@ -392,7 +397,7 @@ def computeEpsilonPerCcd(self):
 
             xyBinHash = xBin[i1a]*(self.deltaAperFitPerCcdNy + 1) + yBin[i1a]
 
-            h2, rev2 = esutil.stat.histogram(xyBinHash, rev=True)
+            h2, rev2 = histogram_rev_sorted(xyBinHash)
 
             gdHash2, = np.where(h2 > 10)
             for j in gdHash2:

fgcm/fgcmExposureSelector.py

@@ -4,6 +4,7 @@
 import esutil
 
 from .fgcmUtilities import expFlagDict, logFlaggedExposuresPerBand, checkFlaggedExposuresPerBand
+from .fgcmUtilities import histogram_rev_sorted
 
 from .sharedNumpyMemManager import SharedNumpyMemManager as snmm
 
@@ -143,8 +144,11 @@ def selectCalibratableNights(self):
         self.fgcmLog.info('Flagging exposures on %d bad nights with too few photometric exposures.' % (badNights.size))
 
         # and we need to use *all* the exposures to flag bad nights
-        h,rev = esutil.stat.histogram(self.fgcmPars.expNightIndex,min=0,
-                                      max=self.fgcmPars.nCampaignNights-1,rev=True)
+        h, rev = histogram_rev_sorted(
+            self.fgcmPars.expNightIndex,
+            min=0,
+            max=self.fgcmPars.nCampaignNights - 1,
+        )
 
         for badNight in badNights:
             i1a=rev[rev[badNight]:rev[badNight+1]]

fgcm/fgcmGray.py

@@ -12,6 +12,7 @@
 from .fgcmUtilities import Cheb2dField
 from .fgcmUtilities import computeDeltaRA
 from .fgcmUtilities import expFlagDict
+from .fgcmUtilities import histogram_rev_sorted
 
 from .sharedNumpyMemManager import SharedNumpyMemManager as snmm
 
@@ -403,7 +404,7 @@ def computeCCDAndExpGray(self, onlyObsErr=False):
             FGrayGO = 10.**(EGrayGO/(-2.5))
             FGrayErrGO = (np.log(10.)/2.5)*np.sqrt(EGrayErr2GO)*FGrayGO
 
-            h, rev = esutil.stat.histogram(obsExpIndex[goodObs], rev=True)
+            h, rev = histogram_rev_sorted(obsExpIndex[goodObs])
 
             use, = np.where(h >= 3)
             for i in use:
@@ -574,7 +575,7 @@ def computeCCDAndExpGray(self, onlyObsErr=False):
             expCcdHash = (obsExpIndex[goodObs]*(self.fgcmPars.nCCD + 1) +
                           obsCCDIndex[goodObs])
 
-            h, rev = esutil.stat.histogram(expCcdHash, rev=True)
+            h, rev = histogram_rev_sorted(expCcdHash)
 
             # Anything with 2 or fewer stars will be marked bad
             use, = np.where(h >= 3)
@@ -1206,7 +1207,7 @@ def computeCCDAndExpDeltaMagBkg(self):
         goodObs, = np.where(obsFlag == 0)
 
         # Do the exposures first
-        h, rev = esutil.stat.histogram(obsExpIndex[goodObs], rev=True)
+        h, rev = histogram_rev_sorted(obsExpIndex[goodObs])
 
         expDeltaMagBkg[:] = self.illegalValue
 
@@ -1228,7 +1229,7 @@ def computeCCDAndExpDeltaMagBkg(self):
         # Do the exp/ccd second
         expCcdHash = obsExpIndex[goodObs]*(self.fgcmPars.nCCD + 1) + obsCCDIndex[goodObs]
 
-        h, rev = esutil.stat.histogram(expCcdHash, rev=True)
+        h, rev = histogram_rev_sorted(expCcdHash)
 
         # We need at least 3 for a median, and of those from the percentile...
         use, = np.where(h > int(3./self.deltaMagBkgOffsetPercentile))
@@ -1289,7 +1290,7 @@ def computeExposureReferenceOffsets(self):
             obsYGO = snmm.getArray(self.fgcmStars.obsYHandle)[goodObs]
             expCcdHash = (obsExpIndexGO[ok] * (self.fgcmPars.nCCD + 1) +
                           obsCCDIndexGO[ok])
-            h, rev = esutil.stat.histogram(expCcdHash, rev=True)
+            h, rev = histogram_rev_sorted(expCcdHash)
             use, = np.where(h > 0)
             for i in use:
                 i1a = rev[rev[i]: rev[i + 1]]
@@ -1318,7 +1319,7 @@ def computeExposureReferenceOffsets(self):
 
         # And then this can be split per exposure.
 
-        h, rev = esutil.stat.histogram(obsExpIndexGO[goodRefObsGO], rev=True)
+        h, rev = histogram_rev_sorted(obsExpIndexGO[goodRefObsGO])
 
         use, = np.where(h >= self.minStarPerExp)
 

fgcm/fgcmMakeStars.py

@@ -6,6 +6,7 @@
 import hpgeom as hpg
 
 from .fgcmLogger import FgcmLogger
+from .fgcmUtilities import histogram_rev_sorted
 
 
 class FgcmMakeStars(object):
@@ -336,14 +337,14 @@ def makePrimaryStars(self, raArray, decArray, filterNameArray,
 
         # Split into pixels
         ipring = hpg.angle_to_pixel(self.starConfig['coarseNSide'], raArray, decArray, nest=False)
-        hpix, revpix = esutil.stat.histogram(ipring, rev=True)
+        hpix, revpix = histogram_rev_sorted(ipring)
 
         gdpix, = np.where(hpix > 0)
         self.fgcmLog.info("Matching primary stars in %d pixels" % (gdpix.size))
 
         for ii, gpix in enumerate(gdpix):
             # This is the array of all the observations in the coarse pixel
-            p1a=revpix[revpix[gpix]: revpix[gpix + 1]]
+            p1a = revpix[revpix[gpix]: revpix[gpix + 1]]
 
             if p1a.size == 0:
                 continue
@@ -445,7 +446,7 @@ def makePrimaryStars(self, raArray, decArray, filterNameArray,
 
                 # This is the official working version, but slower
                 fakeId = np.arange(p1a.size)
-                hist, rev = esutil.stat.histogram(fakeId[i1], rev=True)
+                hist, rev = histogram_rev_sorted(fakeId[i1])
 
                 if (hist.max() == 1):
                     self.fgcmLog.warning("  No matches found for pixel %d, %s band" %
@@ -681,7 +682,7 @@ def makeMatchedStars(self, raArray, decArray, filterNameArray):
             i1 = matches[1]
 
         self.fgcmLog.info("Collating observations")
-        nObsPerObj, obsInd = esutil.stat.histogram(i1, rev=True)
+        nObsPerObj, obsInd = histogram_rev_sorted(i1)
 
         if (nObsPerObj.size != self.objCat.size):
             raise ValueError("Number of primary stars (%d) does not match observations (%d)." %
@@ -729,13 +730,13 @@ def makeMatchedStars(self, raArray, decArray, filterNameArray):
             self.objCat['dec'][gd],
             nest=False
         )
-        hist, rev = esutil.stat.histogram(ipring, rev=True)
+        hist, rev = histogram_rev_sorted(ipring)
 
         high,=np.where(hist > self.starConfig['densMaxPerPixel'])
         ok,=np.where(hist > 0)
         self.fgcmLog.info("There are %d/%d pixels with high stellar density" % (high.size, ok.size))
         for i in range(high.size):
-            i1a=rev[rev[high[i]]:rev[high[i]+1]]
+            i1a = rev[rev[high[i]]: rev[high[i] + 1]]
             cut=np.random.choice(i1a,size=i1a.size-self.starConfig['densMaxPerPixel'],replace=False)
             objClass[gd[cut]] = 0
 
@@ -810,7 +811,7 @@ def makeReferenceMatches(self, refLoader):
             self.objIndexCat['dec'],
             nest=False
         )
-        hpix, revpix = esutil.stat.histogram(ipring, rev=True)
+        hpix, revpix = histogram_rev_sorted(ipring)
 
         pixelCats = []
         nBands = len(self.starConfig['referenceFilterNames'])

fgcm/fgcmQeSysSlope.py

@@ -7,6 +7,8 @@
 import scipy.optimize
 from astropy.time import Time
 
+from .fgcmUtilities import histogram_rev_sorted
+
 from .sharedNumpyMemManager import SharedNumpyMemManager as snmm
 
 class FgcmQeSysSlope(object):
@@ -80,15 +82,15 @@ def computeQeSysSlope(self, name):
         obsMagStdGO -= deltaQESlopeGO
 
         # split per wash interval
-        washH, washRev = esutil.stat.histogram(self.fgcmPars.expWashIndex[obsExpIndexGO], rev=True, min=0)
+        washH, washRev = histogram_rev_sorted(self.fgcmPars.expWashIndex[obsExpIndexGO], min=0)
         washIndices, = np.where(washH > 0)
 
         for washIndex in washIndices:
             i1a = washRev[washRev[washIndex]: washRev[washIndex + 1]]
 
             # Split per band, and compute the delta-T and delta-Mag
 
-            bandH, bandRev = esutil.stat.histogram(obsBandIndex[goodObs[i1a]], rev=True, min=0)
+            bandH, bandRev = histogram_rev_sorted(obsBandIndex[goodObs[i1a]], min=0)
             bandIndices, = np.where(bandH > 0)
 
             deltaTAll = None

fgcm/fgcmRetrieveAtmosphere.py

@@ -10,6 +10,8 @@
 from .sharedNumpyMemManager import SharedNumpyMemManager as snmm
 
 from .fgcmUtilities import retrievalFlagDict
+from .fgcmUtilities import histogram_rev_sorted
+
 
 class FgcmRetrieveAtmosphere(object):
     """
@@ -121,7 +123,7 @@ def r1ToPwv(self, fgcmRetrieval):
 
         # next, we median together each exposure...
         minExpIndex = np.min(expIndexArray[zUse])
-        h, rev = esutil.stat.histogram(expIndexArray[zUse], min=minExpIndex, rev=True)
+        h, rev = histogram_rev_sorted(expIndexArray[zUse], min=minExpIndex)
 
         gd, = np.where(h >= self.minCCDPerExp)
 
@@ -142,8 +144,7 @@ def r1ToPwv(self, fgcmRetrieval):
         # next, we do the median smoothing using pwvRetrievalSmoothBlock
         # self.pwvRetrievalSmoothBlock = fgcmConfig.pwvRetrievalSmoothBlock
 
-        h, rev = esutil.stat.histogram(self.fgcmPars.expNightIndex[rLnPwvStruct['EXPINDEX']],
-                                       rev=True)
+        h, rev = histogram_rev_sorted(self.fgcmPars.expNightIndex[rLnPwvStruct['EXPINDEX']])
 
         # we do this on any night that we have at least 1
         gd, = np.where(h > 0)
@@ -177,7 +178,7 @@ def r1ToPwv(self, fgcmRetrieval):
         nightIndexWithLnPwv = np.unique(self.fgcmPars.expNightIndex[rLnPwvStruct['EXPINDEX']])
 
         a, b = esutil.numpy_util.match(nightIndexWithLnPwv, self.fgcmPars.expNightIndex)
-        h, rev = esutil.stat.histogram(self.fgcmPars.expNightIndex[b], rev=True)
+        h, rev = histogram_rev_sorted(self.fgcmPars.expNightIndex[b])
 
         gd, = np.where(h > 0)
 
@@ -385,8 +386,7 @@ def r0ToNightlyTau(self, fgcmRetrieval):
             extDelta = (-2.5*np.log10(r0[expIndexArray[use], ccdIndexArray[use]]) +
                          2.5*np.log10(I0Ref))
 
-            h,rev = esutil.stat.histogram(self.fgcmPars.expNightIndex[expIndexArray[use]],
-                                          rev=True, min=0)
+            h, rev = histogram_rev_sorted(self.fgcmPars.expNightIndex[expIndexArray[use]], min=0)
 
             gd, = np.where(h > self.tauRetrievalMinCCDPerNight)
             if not self.quietMode:
@@ -534,8 +534,7 @@ def slopeFunc(X, slope):
             extDelta = (-2.5*np.log10(r0Gray[use]) +
                          2.5*np.log10(I0Ref))
 
-            h,rev = esutil.stat.histogram(self.fgcmPars.expNightIndex[expIndexArray[use]],
-                                          rev=True, min=0)
+            h, rev = histogram_rev_sorted(self.fgcmPars.expNightIndex[expIndexArray[use]], min=0)
 
             gd, = np.where(h > self.tauRetrievalMinCCDPerNight)
             if not self.quietMode:

fgcm/fgcmSigmaCal.py

@@ -5,6 +5,7 @@
 import time
 
 from .fgcmUtilities import retrievalFlagDict
+from .fgcmUtilities import histogram_rev_sorted
 
 import multiprocessing
 
@@ -209,8 +210,11 @@ def run(self, applyGray=True):
                         continue
 
                     # These have already been limited to the plot percentile range
-                    h, rev = esutil.stat.histogram(objMagStdMean[goodStars[plotIndices[band][ok]], bandIndex],
-                                                   nbin=nPlotBin, rev=True)
+                    h, rev = histogram_rev_sorted(
+                        objMagStdMean[goodStars[plotIndices[band][ok]], bandIndex],
+                        nbin=nPlotBin,
+                    )
+
                     for j, nInBin in enumerate(h):
                         if nInBin < 100:
                             continue

fgcm/fgcmStars.py

@@ -8,6 +8,7 @@
 from .fgcmUtilities import objFlagDict
 from .fgcmUtilities import obsFlagDict
 from .fgcmUtilities import getMemoryString
+from .fgcmUtilities import histogram_rev_sorted
 
 from .sharedNumpyMemManager import SharedNumpyMemManager as snmm
 
@@ -1639,7 +1640,7 @@ def performSuperStarOutlierCuts(self, fgcmPars, reset=False):
                            (fgcmPars.nCCD+1) +
                            obsCCDIndex[goodObs])
 
-        h, rev = esutil.stat.histogram(epochFilterHash, rev=True)
+        h, rev = histogram_rev_sorted(epochFilterHash)
 
         nbad = 0
 
@@ -1711,7 +1712,7 @@ def performFocalPlaneOutlierCuts(self, fgcmPars, reset=False, ignoreRef=False):
         # compute EGray, GO for Good Obs
         EGrayGO, EGrayErr2GO = self.computeEGray(goodObs, onlyObsErr=True, ignoreRef=ignoreRef)
 
-        h, rev = esutil.stat.histogram(obsExpIndex[goodObs], rev=True)
+        h, rev = histogram_rev_sorted(obsExpIndex[goodObs])
 
         nbad = 0
 
@@ -1770,7 +1771,7 @@ def applySuperStarFlat(self,fgcmPars):
                                (fgcmPars.nCCD+1) +
                                obsCCDIndex)
 
-            h, rev = esutil.stat.histogram(epochFilterHash, rev=True)
+            h, rev = histogram_rev_sorted(epochFilterHash)
 
             for i in range(h.size):
                 if h[i] == 0: continue