-
Notifications
You must be signed in to change notification settings - Fork 3
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
DM-42285: Update to support numpy 1.25+ and various optimizations. #39
Changes from all commits
deafd5c
8e07d65
602d56a
8d088c2
367497e
30f501f
0934113
3f4cbd2
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -1,3 +1,3 @@ | ||
__version__ = '3.10.3' | ||
__version__ = '3.11.0' | ||
|
||
__version_info__ = __version__.split('.') |
Large diffs are not rendered by default.
Original file line number | Diff line number | Diff line change |
---|---|---|
|
@@ -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 | ||
|
||
|
@@ -226,10 +227,10 @@ def __call__(self, fitPars, returnExpGray=False): | |
|
||
np.add.at(expGrayColorSplit[:, k], | ||
self.obsExpIndexGO[sel], | ||
EGrayGOS / self.EGrayErr2GO[sel]) | ||
(EGrayGOS / self.EGrayErr2GO[sel]).astype(expGrayColorSplit.dtype)) | ||
np.add.at(expGrayColorSplitWt[:, k], | ||
self.obsExpIndexGO[sel], | ||
1. / self.EGrayErr2GO[sel]) | ||
(1. / self.EGrayErr2GO[sel]).astype(expGrayColorSplitWt.dtype)) | ||
|
||
ok, = np.where(expGrayColorSplitWt[:, k] > 0.0) | ||
expGrayColorSplit[ok, k] /= expGrayColorSplitWt[ok, k] | ||
|
@@ -240,14 +241,14 @@ def __call__(self, fitPars, returnExpGray=False): | |
|
||
np.add.at(dExpGraydC0ColorSplit[:, :, k], | ||
(self.obsExpIndexGO[sel], self.fgcmPars.expCoatingIndex[self.obsExpIndexGO[sel]]), | ||
dDeltadC0GOS / self.EGrayErr2GO[sel]) | ||
(dDeltadC0GOS / self.EGrayErr2GO[sel]).astype(dExpGraydC0ColorSplit.dtype)) | ||
|
||
for i in range(c0s.size): | ||
dExpGraydC0ColorSplit[ok, i, k] /= expGrayColorSplitWt[ok, k] | ||
|
||
np.add.at(dExpGraydC1ColorSplit[:, k], | ||
self.obsExpIndexGO[sel], | ||
dDeltadC0GOS * self.deltaTGO[sel] / self.EGrayErr2GO[sel]) | ||
(dDeltadC0GOS * self.deltaTGO[sel] / self.EGrayErr2GO[sel]).astype(dExpGraydC1ColorSplit.dtype)) | ||
|
||
dExpGraydC1ColorSplit[ok, k] /= expGrayColorSplitWt[ok, k] | ||
|
||
|
@@ -266,7 +267,7 @@ def __call__(self, fitPars, returnExpGray=False): | |
(2.0 * (expGrayColorSplitWt[ok, 0] + expGrayColorSplitWt[ok, 1]) * | ||
(expGrayColorSplit[ok, 0] - expGrayColorSplit[ok, 1]) * | ||
(dExpGraydC0ColorSplit[ok, self.fgcmPars.expCoatingIndex[ok], 0] - | ||
dExpGraydC0ColorSplit[ok, self.fgcmPars.expCoatingIndex[ok], 1]))) | ||
dExpGraydC0ColorSplit[ok, self.fgcmPars.expCoatingIndex[ok], 1])).astype(deriv.dtype)) | ||
|
||
deriv /= (ok.size - len(fitPars)) | ||
|
||
|
@@ -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) | ||
|
@@ -382,7 +383,7 @@ def computeCCDChromaticity(self): | |
self.fgcmLog.warning("Found out-of-bounds value for chromaticity for filter %s, detector %d." % (self.fgcmPars.lutFilterNames[self.fInd], self.cInd + self.ccdStartIndex)) | ||
continue | ||
|
||
self.fgcmPars.compCCDChromaticity[self.cInd, self.fInd] = res.x | ||
self.fgcmPars.compCCDChromaticity[self.cInd, self.fInd] = res.x[0] | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Do you also need to add the There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. No, that's fine. The problem is that numpy doesn't like you setting one value with an array of length 1 anymore. Not sure why, but it doesn't like it. Comparisons to an array of length 1 are fine. |
||
cWasFit[self.cInd, self.fInd] = True | ||
|
||
# And make plots if necessary. | ||
|
Original file line number | Diff line number | Diff line change |
---|---|---|
|
@@ -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: | ||
|
@@ -517,10 +522,10 @@ def _starWorker(self, goodStarsAndObs): | |
|
||
np.add.at(objDeltaAperMeanTemp, | ||
(obsObjIDIndex[goodObs], obsBandIndex[goodObs]), | ||
(obsDeltaAper[goodObs] - self.fgcmPars.compMedDeltaAper[obsExpIndex[goodObs]])/obsMagErr2GO) | ||
((obsDeltaAper[goodObs] - self.fgcmPars.compMedDeltaAper[obsExpIndex[goodObs]])/obsMagErr2GO).astype(objDeltaAperMeanTemp.dtype)) | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Formatting? There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. 🤷♂️ There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Oh, it may have just looked funny on my screen. Ignore! |
||
np.add.at(wtSum, | ||
(obsObjIDIndex[goodObs], obsBandIndex[goodObs]), | ||
1./obsMagErr2GO) | ||
(1./obsMagErr2GO).astype(wtSum.dtype)) | ||
|
||
gd = np.where(wtSum > 0.0) | ||
|
||
|
Original file line number | Diff line number | Diff line change |
---|---|---|
|
@@ -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 | ||
|
||
|
@@ -210,10 +211,10 @@ def computeExpGrayForInitialSelection(self): | |
|
||
np.add.at(expGrayForInitialSelection, | ||
obsExpIndex[goodObs], | ||
EGray[goodObs]) | ||
(EGray[goodObs]).astype(expGrayForInitialSelection.dtype)) | ||
np.add.at(expGrayRMSForInitialSelection, | ||
obsExpIndex[goodObs], | ||
EGray[goodObs]**2.) | ||
(EGray[goodObs]).astype(expGrayRMSForInitialSelection.dtype)**2.) | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Should the There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. The opposite, I think. It should be promoted before squaring. There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Ok...you may want to look for other cases (I think I saw others where the promotion was after squaring). |
||
np.add.at(expNGoodStarForInitialSelection, | ||
obsExpIndex[goodObs], | ||
1) | ||
|
@@ -375,7 +376,7 @@ def computeCCDAndExpGray(self, onlyObsErr=False): | |
|
||
np.add.at(ccdGrayWt, | ||
(obsExpIndex[goodObs],obsCCDIndex[goodObs]), | ||
1./EGrayErr2GO) | ||
(1./EGrayErr2GO).astype(ccdGrayWt.dtype)) | ||
np.add.at(ccdNGoodStars, | ||
(obsExpIndex[goodObs],obsCCDIndex[goodObs]), | ||
1) | ||
|
@@ -385,7 +386,7 @@ def computeCCDAndExpGray(self, onlyObsErr=False): | |
obsBandIndex[goodObs]]) | ||
np.add.at(ccdDeltaStd, | ||
(obsExpIndex[goodObs], obsCCDIndex[goodObs]), | ||
obsDeltaStdGO/EGrayErr2GO) | ||
(obsDeltaStdGO/EGrayErr2GO).astype(ccdDeltaStd.dtype)) | ||
|
||
gd = np.where((ccdNGoodStars >= 3) & (ccdGrayWt > 0.0)) | ||
ccdDeltaStd[gd] /= ccdGrayWt[gd] | ||
|
@@ -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: | ||
|
@@ -462,10 +463,10 @@ def computeCCDAndExpGray(self, onlyObsErr=False): | |
|
||
np.add.at(ccdGrayTemp, | ||
obsCCDIndex[goodObs[i1a]], | ||
EGrayGO[i1a]/EGrayErr2GO[i1a]) | ||
(EGrayGO[i1a]/EGrayErr2GO[i1a]).astype(ccdGrayTemp.dtype)) | ||
np.add.at(ccdGrayRMSTemp, | ||
obsCCDIndex[goodObs[i1a]], | ||
EGrayGO[i1a]**2./EGrayErr2GO[i1a]) | ||
(EGrayGO[i1a]**2./EGrayErr2GO[i1a]).astype(ccdGrayRMSTemp.dtype)) | ||
|
||
gdTemp, = np.where((ccdNGoodStars[eInd, :] >= 3) & | ||
(ccdGrayWt[eInd, :] > 0.0) & | ||
|
@@ -494,7 +495,7 @@ def computeCCDAndExpGray(self, onlyObsErr=False): | |
ccdGrayEvalNStars = np.zeros(self.fgcmPars.nCCD, dtype=np.int32) | ||
np.add.at(ccdGrayEval, | ||
obsCCDIndex[goodObs[i1a]], | ||
ccdGrayEvalStars) | ||
ccdGrayEvalStars.astype(ccdGrayEval.dtype)) | ||
np.add.at(ccdGrayEvalNStars, | ||
obsCCDIndex[goodObs[i1a]], | ||
1) | ||
|
@@ -535,10 +536,10 @@ def computeCCDAndExpGray(self, onlyObsErr=False): | |
# we can do all of this at once. | ||
np.add.at(ccdGray, | ||
(obsExpIndex[goodObs],obsCCDIndex[goodObs]), | ||
EGrayGO/EGrayErr2GO) | ||
(EGrayGO/EGrayErr2GO).astype(ccdGray.dtype)) | ||
np.add.at(ccdGrayRMS, | ||
(obsExpIndex[goodObs],obsCCDIndex[goodObs]), | ||
EGrayGO**2./EGrayErr2GO) | ||
(EGrayGO**2./EGrayErr2GO).astype(ccdGrayRMS.dtype)) | ||
|
||
# need at least 3 or else computation can blow up | ||
gd = np.where((ccdNGoodStars >= 3) & (ccdGrayWt > 0.0) & (ccdGrayRMS > 0.0)) | ||
|
@@ -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) | ||
|
@@ -678,25 +679,25 @@ def computeCCDAndExpGray(self, onlyObsErr=False): | |
|
||
np.add.at(expGrayWt, | ||
goodCCD[0], | ||
1./ccdGrayErr[goodCCD]**2.) | ||
(1./ccdGrayErr[goodCCD]**2.).astype(ccdGrayWt.dtype)) | ||
np.add.at(expGray, | ||
goodCCD[0], | ||
ccdGray[goodCCD]/ccdGrayErr[goodCCD]**2.) | ||
(ccdGray[goodCCD]/ccdGrayErr[goodCCD]**2.).astype(expGray.dtype)) | ||
np.add.at(expGrayRMS, | ||
goodCCD[0], | ||
ccdGray[goodCCD]**2./ccdGrayErr[goodCCD]**2.) | ||
(ccdGray[goodCCD]**2./ccdGrayErr[goodCCD]**2.).astype(expGrayRMS.dtype)) | ||
np.add.at(expNGoodCCDs, | ||
goodCCD[0], | ||
1) | ||
np.add.at(expNGoodTilings, | ||
goodCCD[0], | ||
ccdNGoodTilings[goodCCD]) | ||
(ccdNGoodTilings[goodCCD]).astype(expNGoodTilings.dtype)) | ||
np.add.at(expNGoodStars, | ||
goodCCD[0], | ||
ccdNGoodStars[goodCCD]) | ||
(ccdNGoodStars[goodCCD]).astype(expNGoodStars.dtype)) | ||
np.add.at(expDeltaStd, | ||
goodCCD[0], | ||
ccdDeltaStd[goodCCD]/ccdGrayErr[goodCCD]**2.) | ||
(ccdDeltaStd[goodCCD]/ccdGrayErr[goodCCD]**2.).astype(expDeltaStd.dtype)) | ||
|
||
if self.fgcmPars.nCCD >= 3: | ||
# Regular mode, when we have a multi-detector camera. | ||
|
@@ -811,13 +812,13 @@ def computeExpGrayColorSplit(self): | |
|
||
np.add.at(expGrayColorSplit[:, c], | ||
obsExpIndex[goodObs[use]], | ||
EGrayGO[use] / EGrayErr2GO[use]) | ||
(EGrayGO[use] / EGrayErr2GO[use]).astype(expGrayColorSplit.dtype)) | ||
np.add.at(expGrayWtColorSplit[:, c], | ||
obsExpIndex[goodObs[use]], | ||
1. / EGrayErr2GO[use]) | ||
(1. / EGrayErr2GO[use]).astype(expGrayWtColorSplit.dtype)) | ||
np.add.at(expGrayRMSColorSplit[:, c], | ||
obsExpIndex[goodObs[use]], | ||
EGrayGO[use]**2. / EGrayErr2GO[use]) | ||
(EGrayGO[use]**2. / EGrayErr2GO[use]).astype(expGrayRMSColorSplit.dtype)) | ||
np.add.at(expGrayNGoodStarsColorSplit[:, c], | ||
obsExpIndex[goodObs[use]], | ||
1) | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Is |
||
|
@@ -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) | ||
|
||
|
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Just checking that this doesn't need to be
astype(deriv[1: ].dtype)
?There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Nope, it's just the dtype of the full array that matters.