Changes to deal with numpy 1.25 deprecation of treating 1 element arr…

…ays as scalars

galpy/actionAngle/actionAngleAdiabatic.py

@@ -152,9 +152,9 @@ def _evaluate(self, *args, **kwargs):
                 for ii in range(len(R)):
                     targs = (R[ii], vR[ii], vT[ii], z[ii], vz[ii])
                     tjr, tlz, tjz = self(*targs, **copy.copy(kwargs))
-                    ojr[ii] = tjr
-                    ojz[ii] = tjz
-                    olz[ii] = tlz
+                    ojr[ii] = tjr[0]
+                    ojz[ii] = tjz[0]
+                    olz[ii] = tlz[0]
                 return (ojr, olz, ojz)
             else:
                 if kwargs.get("_justjr", False):
@@ -180,7 +180,11 @@ def _evaluate(self, *args, **kwargs):
                     )
                 else:
                     axiJ = self._aAS(R[0], vR[0], vT[0], 0.0, 0.0, _Jz=Jz)
-                    return (axiJ[0], axiJ[1], Jz)
+                    return (
+                        numpy.atleast_1d(axiJ[0]),
+                        numpy.atleast_1d(axiJ[1]),
+                        numpy.atleast_1d(Jz),
+                    )
 
     def _EccZmaxRperiRap(self, *args, **kwargs):
         """
@@ -254,10 +258,10 @@ def _EccZmaxRperiRap(self, *args, **kwargs):
                     tecc, tzmax, trperi, trap = self._EccZmaxRperiRap(
                         *targs, **copy.copy(kwargs)
                     )
-                    oecc[ii] = tecc
-                    ozmax[ii] = tzmax
-                    orperi[ii] = trperi
-                    orap[ii] = trap
+                    oecc[ii] = tecc[0]
+                    ozmax[ii] = tzmax[0]
+                    orperi[ii] = trperi[0]
+                    orap[ii] = trap[0]
                 return (oecc, ozmax, orperi, orap)
             else:
                 if _dim(self._pot) == 3:

galpy/actionAngle/actionAngleAdiabaticGrid.py

@@ -142,7 +142,7 @@ def __init__(
                             numpy.sqrt(2.0 * y[jj] * self._EzZmaxs[ii]),
                             _justjz=True,
                             **kwargs
-                        )[2]
+                        )[2][0]
                         if jj == nEz - 1:
                             jzEzzmax[ii] = jz[ii, jj]
         for ii in range(nR):
@@ -261,7 +261,7 @@ def __init__(
                                 0.0,
                                 _justjr=True,
                                 **kwargs
-                            )[0]
+                            )[0][0]
                         except UnboundError:  # pragma: no cover
                             raise
                         if jj == 0:

galpy/actionAngle/actionAngleStaeckel.py

@@ -200,9 +200,9 @@ def _evaluate(self, *args, **kwargs):
                     except (TypeError, IndexError):
                         tkwargs["delta"] = delta
                     tjr, tlz, tjz = self(*targs, **tkwargs)
-                    ojr[ii] = tjr
-                    ojz[ii] = tjz
-                    olz[ii] = tlz
+                    ojr[ii] = tjr[0]
+                    ojz[ii] = tjz[0]
+                    olz[ii] = tlz[0]
                 return (ojr, olz, ojz)
             else:
                 # Set up the actionAngleStaeckelSingle object
@@ -481,7 +481,7 @@ def _uminumaxvmin(self, *args, **kwargs):
         HISTORY:
            2017-12-12 - Written - Bovy (UofT)
         """
-        delta = kwargs.pop("delta", self._delta)
+        delta = numpy.atleast_1d(kwargs.pop("delta", self._delta))
         if len(args) == 5:  # R,vR.vT, z, vz
             R, vR, vT, z, vz = args
         elif len(args) == 6:  # R,vR.vT, z, vz, phi
@@ -552,19 +552,16 @@ def _uminumaxvmin(self, *args, **kwargs):
                 for ii in range(len(R)):
                     targs = (R[ii], vR[ii], vT[ii], z[ii], vz[ii])
                     tkwargs = copy.copy(kwargs)
-                    try:
-                        tkwargs["delta"] = delta[ii]
-                    except TypeError:
-                        tkwargs["delta"] = delta
+                    tkwargs["delta"] = delta[ii] if len(delta) > 1 else delta[0]
                     tumin, tumax, tvmin = self._uminumaxvmin(*targs, **tkwargs)
-                    oumin[ii] = tumin
-                    oumax[ii] = tumax
-                    ovmin[ii] = tvmin
+                    oumin[ii] = tumin[0]
+                    oumax[ii] = tumax[0]
+                    ovmin[ii] = tvmin[0]
                 return (oumin, oumax, ovmin)
             else:
                 # Set up the actionAngleStaeckelSingle object
                 aASingle = actionAngleStaeckelSingle(
-                    R[0], vR[0], vT[0], z[0], vz[0], pot=self._pot, delta=delta
+                    R[0], vR[0], vT[0], z[0], vz[0], pot=self._pot, delta=delta[0]
                 )
                 umin, umax = aASingle.calcUminUmax()
                 vmin = aASingle.calcVmin()
@@ -975,8 +972,8 @@ def calcUminUmax(self, **kwargs):
                     try:
                         umin = optimize.brentq(
                             _JRStaeckelIntegrandSquared,
-                            rstart,
-                            self._ux - eps,
+                            numpy.atleast_1d(rstart)[0],
+                            numpy.atleast_1d(self._ux)[0] - eps,
                             (
                                 E,
                                 L,
@@ -1011,8 +1008,8 @@ def calcUminUmax(self, **kwargs):
                 )
                 umax = optimize.brentq(
                     _JRStaeckelIntegrandSquared,
-                    self._ux + eps,
-                    rend,
+                    numpy.atleast_1d(self._ux)[0] + eps,
+                    numpy.atleast_1d(rend)[0],
                     (
                         E,
                         L,

galpy/actionAngle/actionAngleVertical.py

@@ -55,13 +55,6 @@ def __init__(self, *args, **kwargs):
             raise OSError("Must specify pot= for actionAngleVertical")
         self._pot = kwargs["pot"]
         return None
-        """
-        self._parse_eval_args(*args,_noOrbUnitsCheck=True,**kwargs)
-        self._z= self._eval_z
-        self._vz= self._eval_vz
-        self._verticalpot= kwargs['pot']
-        return None
-        """
 
     def _evaluate(self, *args, **kwargs):
         """

galpy/actionAngle/actionAngleVerticalInverse.py

@@ -96,14 +96,14 @@ def __init__(
                     ),
                 )
                 js[ii] = 0.0
-                Omegas[ii] = tO
+                Omegas[ii] = tO[0]
                 xmaxs[ii] = 0.0
                 continue
             tJ, tO = self._aAV.actionsFreqs(
                 0.0, numpy.sqrt(2.0 * (E - evaluatelinearPotentials(self._pot, 0.0)))
             )
-            js[ii] = tJ
-            Omegas[ii] = tO
+            js[ii] = tJ[0]
+            Omegas[ii] = tO[0]
             xmaxs[ii] = self._aAV.calcxmax(
                 0.0,
                 numpy.sqrt(2.0 * (E - evaluatelinearPotentials(self._pot, 0.0))),

galpy/df/constantbetaHernquistdf.py

@@ -121,7 +121,7 @@ def fE(self, E):
             )
         if len(Etilde_out) > 0:
             fE[Etilde_out] = 0.0
-        return fE
+        return fE.reshape(E.shape)
 
     def _icmf(self, ms):
         """Analytic expression for the normalized inverse cumulative mass

galpy/df/constantbetadf.py

@@ -343,7 +343,7 @@ def fE(self, E):
         if self._halfint:
             # fE is simply given by the relevant derivative
             out[indx] = self._gradfunc(self._rphi(Eint[indx]))
-            return out / (
+            return out.reshape(E.shape) / (
                 2.0
                 * numpy.pi**1.5
                 * 2 ** (0.5 - self._beta)
@@ -392,7 +392,7 @@ def fE(self, E):
                     for tE in Eint[indx]
                 ]
             )
-            return -out * self._fE_prefactor
+            return -out.reshape(E.shape) * self._fE_prefactor
 
 
 def _fEintegrand_raw(r, pot, E, dmp1nudrmp1, alpha):

galpy/df/isotropicHernquistdf.py

@@ -95,7 +95,7 @@ def fE(self, E):
         # Fix out of bounds values
         if len(Etilde_out) > 0:
             fE[Etilde_out] = 0
-        return fE
+        return fE.reshape(E.shape)
 
     def _dMdE(self, E):
         # E already in internal units here

galpy/df/kingdf.py

@@ -109,7 +109,7 @@ def fE(self, E):
                 * (2.0 * numpy.pi * self._sigma2) ** -1.5
                 * self.rho1
             )
-        return out  # mass density, not /self.M as for number density
+        return out.reshape(E.shape)  # mass density, not /self.M as for number density
 
 
 class _scalefreekingdf:

galpy/df/osipkovmerrittHernquistdf.py

@@ -106,7 +106,7 @@ def fQ(self, Q):
         fQ += 8.0 * self._a2overra2 * sqrtQtilde * (1.0 - 2.0 * Qtilde)
         if len(Qtilde_out) > 0:
             fQ[Qtilde_out] = 0.0
-        return self._fQnorm * fQ
+        return self._fQnorm * fQ.reshape(Q.shape)
 
     def _icmf(self, ms):
         """Analytic expression for the normalized inverse cumulative mass

galpy/df/osipkovmerrittdf.py

@@ -141,7 +141,7 @@ def Lintegrand(t, L2lim, E):
         # Numerical issues can make the integrand's sqrt argument negative, only
         # happens at dMdE ~ 0, so just set to zero
         out[numpy.isnan(out)] = 0.0
-        return out
+        return out.reshape(E.shape)
 
     def _sample_eta(self, r, n=1):
         """Sample the angle eta which defines radial vs tangential velocities"""

galpy/df/sphericaldf.py

@@ -166,7 +166,17 @@ def __call__(self, *args, **kwargs):
             r = numpy.sqrt(R**2.0 + z**2.0)
             vrad = (R * vR + z * vz) / r
             L = numpy.atleast_1d(numpy.sqrt(vtotSq - vrad**2.0) * r)
-        return self._call_internal(E, L, Lz)  # Some function for each sub-class
+        return self._call_internal(E, L, Lz).reshape(
+            args[0].shape
+            if len(args) == 1 and hasattr(args[0], "shape")
+            else (
+                args[0][0].shape
+                if len(args) == 1
+                and hasattr(args[0], "__len__")
+                and hasattr(args[0][0], "shape")
+                else (args[0].shape if hasattr(args[0], "shape") else ())
+            )
+        )
 
     @physical_conversion("energydensity", pop=True)
     def dMdE(self, E):
@@ -192,7 +202,9 @@ def dMdE(self, E):
             2023-05-23 - Written - Bovy (UofT)
 
         """
-        return self._dMdE(numpy.atleast_1d(conversion.parse_energy(E, vo=self._vo)))
+        return self._dMdE(
+            numpy.atleast_1d(conversion.parse_energy(E, vo=self._vo))
+        ).reshape(E.shape if isinstance(E, numpy.ndarray) else ())
 
     def vmomentdensity(self, r, n, m, **kwargs):
         """

galpy/df/streamdf.py

@@ -279,15 +279,6 @@ def _progenitor_setup(self, progenitor, leading, useTMHessian):
                 self._progenitor_jr, self._progenitor_lz, self._progenitor_jz
             )
             self._dOdJp = h
-            print(
-                self._progenitor_jr,
-                self._progenitor_lz,
-                self._progenitor_jz,
-                h,
-                fr,
-                fp,
-                fz,
-            )
             # Replace frequencies with TM frequencies
             self._progenitor_Omegar = fr
             self._progenitor_Omegaphi = fp
@@ -2169,9 +2160,8 @@ def pOparapar(self, Opar, apar, tdisrupt=None):
         apar = conversion.parse_angle(apar)
         if tdisrupt is None:
             tdisrupt = self._tdisrupt
-        if isinstance(Opar, (int, float, numpy.float32, numpy.float64)):
-            Opar = numpy.array([Opar])
-        out = numpy.zeros(len(Opar))
+        Opar = numpy.array(Opar)
+        out = numpy.zeros_like(Opar)
         # Compute ts
         ts = apar / Opar
         # Evaluate
@@ -2535,9 +2525,8 @@ def ptdAngle(self, t, dangle):
            2013-12-05 - Written - Bovy (IAS)
 
         """
-        if isinstance(t, (int, float, numpy.float32, numpy.float64)):
-            t = numpy.array([t])
-        out = numpy.zeros(len(t))
+        t = numpy.array(t)
+        out = numpy.zeros_like(t)
         dO = dangle / t[(t > 0.0) * (t < self._tdisrupt)]
         # p(t|a) = \int dO p(O,t|a) = \int dO p(t|O,a) p(O|a) = \int dO delta (t-a/O)p(O|a) = O*2/a p(O|a); p(O|a) = \int dt p(a|O,t) p(O)p(t) = 1/O p(O)
         out[(t > 0.0) * (t < self._tdisrupt)] = (
@@ -2647,18 +2636,17 @@ def pangledAngle(self, angleperp, dangle, smallest=False):
            2013-12-06 - Written - Bovy (IAS)
 
         """
-        if isinstance(angleperp, (int, float, numpy.float32, numpy.float64)):
-            angleperp = numpy.array([angleperp])
-        out = numpy.zeros(len(angleperp))
+        angleperp = numpy.array(angleperp)
+        out = numpy.zeros_like(angleperp)
         out = numpy.array(
             [
                 integrate.quad(
                     self._pangledAnglet, 0.0, self._tdisrupt, (ap, dangle, smallest)
                 )[0]
-                for ap in angleperp
+                for ap in numpy.atleast_1d(angleperp)
             ]
         )
-        return out
+        return out.reshape(angleperp.shape)
 
     @physical_conversion("angle", pop=True)
     def meanangledAngle(self, dangle, smallest=False):
@@ -2771,9 +2759,8 @@ def _pangledAnglet(self, t, angleperp, dangle, smallest):
             eigIndx = 0
         else:
             eigIndx = 1
-        if isinstance(angleperp, (int, float, numpy.float32, numpy.float64)):
-            angleperp = numpy.array([angleperp])
-            t = numpy.array([t])
+        angleperp = numpy.array(angleperp)
+        t = numpy.array(t)
         out = numpy.zeros_like(angleperp)
         tindx = t < self._tdisrupt
         out[tindx] = (
@@ -3732,7 +3719,7 @@ def _determine_stream_track_single(
     allAcfsTrack[1] = tacfs[1][0]
     allAcfsTrack[2] = tacfs[2][0]
     for jj in range(3, 9):
-        allAcfsTrack[jj] = tacfs[jj]
+        allAcfsTrack[jj] = tacfs[jj][0]
     tjac = calcaAJac(
         progenitorTrack(trackt).vxvv[0],
         aA,
@@ -4010,35 +3997,35 @@ def calcaAJac(
         xv[ii] -= dxv[ii]
         angleIndx = 3
         if actionsFreqsAngles:
-            jac[0, ii] = (tjr - jr) / dxv[ii]
-            jac[1, ii] = (tlz - lz) / dxv[ii]
-            jac[2, ii] = (tjz - jz) / dxv[ii]
-            jac[3, ii] = (tOr - Or) / dxv[ii]
-            jac[4, ii] = (tOphi - Ophi) / dxv[ii]
-            jac[5, ii] = (tOz - Oz) / dxv[ii]
+            jac[0, ii] = (tjr - jr)[0] / dxv[ii]
+            jac[1, ii] = (tlz - lz)[0] / dxv[ii]
+            jac[2, ii] = (tjz - jz)[0] / dxv[ii]
+            jac[3, ii] = (tOr - Or)[0] / dxv[ii]
+            jac[4, ii] = (tOphi - Ophi)[0] / dxv[ii]
+            jac[5, ii] = (tOz - Oz)[0] / dxv[ii]
             angleIndx = 6
         elif freqs:
-            jac[0, ii] = (tOr - Or) / dxv[ii]
-            jac[1, ii] = (tOphi - Ophi) / dxv[ii]
-            jac[2, ii] = (tOz - Oz) / dxv[ii]
+            jac[0, ii] = (tOr - Or)[0] / dxv[ii]
+            jac[1, ii] = (tOphi - Ophi)[0] / dxv[ii]
+            jac[2, ii] = (tOz - Oz)[0] / dxv[ii]
         else:
-            jac[0, ii] = (tjr - jr) / dxv[ii]
-            jac[1, ii] = (tlz - lz) / dxv[ii]
-            jac[2, ii] = (tjz - jz) / dxv[ii]
+            jac[0, ii] = (tjr - jr)[0] / dxv[ii]
+            jac[1, ii] = (tlz - lz)[0] / dxv[ii]
+            jac[2, ii] = (tjz - jz)[0] / dxv[ii]
         if dOdJ:
-            jac2[0, ii] = (tOr - Or) / dxv[ii]
-            jac2[1, ii] = (tOphi - Ophi) / dxv[ii]
-            jac2[2, ii] = (tOz - Oz) / dxv[ii]
+            jac2[0, ii] = (tOr - Or)[0] / dxv[ii]
+            jac2[1, ii] = (tOphi - Ophi)[0] / dxv[ii]
+            jac2[2, ii] = (tOz - Oz)[0] / dxv[ii]
         # For the angles, make sure we do not hit a turning point
-        jac[angleIndx, ii] = (
-            (tar - ar + numpy.pi) % (2.0 * numpy.pi) - numpy.pi
-        ) / dxv[ii]
+        jac[angleIndx, ii] = ((tar - ar + numpy.pi) % (2.0 * numpy.pi) - numpy.pi)[
+            0
+        ] / dxv[ii]
         jac[angleIndx + 1, ii] = (
             (taphi - aphi + numpy.pi) % (2.0 * numpy.pi) - numpy.pi
-        ) / dxv[ii]
-        jac[angleIndx + 2, ii] = (
-            (taz - az + numpy.pi) % (2.0 * numpy.pi) - numpy.pi
-        ) / dxv[ii]
+        )[0] / dxv[ii]
+        jac[angleIndx + 2, ii] = ((taz - az + numpy.pi) % (2.0 * numpy.pi) - numpy.pi)[
+            0
+        ] / dxv[ii]
     if dOdJ:
         jac2[3, :] = jac[3, :]
         jac2[4, :] = jac[4, :]