Index ndarrays with tuples, not lists

This used to be a warning,

```
FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.
```

but turned into a hard error some time after numpy 1.21

toise/diffuse.py

@@ -293,7 +293,7 @@ def point_source_background(
         if n_sources is not None:
             expand = [None] * bin_areas.ndim
             expand[0] = slice(None)
-            bin_areas = bin_areas * n_sources[expand]
+            bin_areas = bin_areas * n_sources[tuple(expand)]
 
         # dimensions of the keys in expectations are now energy, radial bin
         if is_zenith_weight(zenith_index, self._aeff):
@@ -475,7 +475,7 @@ def point_source_background(
         psi_bins = self._aeff.bin_edges[-1][:-1]
         expand = [None] * 5
         expand[-1] = slice(None)
-        bin_areas = (np.pi * np.diff(psi_bins**2))[expand]
+        bin_areas = (np.pi * np.diff(psi_bins**2))[tuple(expand)]
         # observation time shorter for triggered transient searches
         if livetime is not None:
             background._livetime = livetime / constants.annum
@@ -494,7 +494,7 @@ def point_source_background(
         if n_sources is not None:
             expand = [None] * bin_areas.ndim
             expand[2] = slice(None)
-            bin_areas = bin_areas * n_sources[expand]
+            bin_areas = bin_areas * n_sources[tuple(expand)]
 
         # dimensions of flux are now 1/m^2 sr
         if isinstance(sel, slice):

toise/fictive.py

@@ -86,7 +86,7 @@ def get_aeff(angular_resolution=0.5, energy_threshold=1e3):
     edges = (e_nu, cos_theta, e_mu[idx:], psi_bins)
 
     return effective_areas.effective_area(
-        edges, base.values[..., idx:, None] * smear[expand], "cos_theta"
+        edges, base.values[..., idx:, None] * smear[tuple(expand)], "cos_theta"
     )
 
 

toise/pointsource.py

@@ -41,7 +41,7 @@ def __init__(self, effective_area, fluence, zenith_selection, with_energy=True):
             expand = [None] * effective_area.values.ndim
             expand[zenith_dim] = slice(None)
             effective_area = (
-                effective_area.values[..., :-1] * zenith_selection[expand]
+                effective_area.values[..., :-1] * zenith_selection[tuple(expand)]
             ).sum(axis=zenith_dim)
         else:
             effective_area = effective_area.values[..., zenith_selection, :, :-1]
@@ -82,7 +82,7 @@ def expectations(self, **kwargs):
         expand[1] = slice(None)
 
         # FIXME: this still neglects the opening angle between neutrino and muon
-        total = (self._rate * (specweight[expand])).sum(axis=(0, 1))
+        total = (self._rate * (specweight[tuple(expand)])).sum(axis=(0, 1))
         # assert total.ndim == 2
 
         if not self._use_energies:

toise/radio_aeff_generation.py

@@ -786,7 +786,7 @@ def point_source_background(
         if n_sources is not None:
             expand = [None] * bin_areas.ndim
             expand[0] = slice(None)
-            bin_areas = bin_areas * n_sources[expand]
+            bin_areas = bin_areas * n_sources[tuple(expand)]
 
         # dimensions of the keys in expectations are now energy, radial bin
         if is_zenith_weight(zenith_index, self._aeff):

toise/surface_veto.py

@@ -498,7 +498,7 @@ def point_source_background(
         if n_sources is not None:
             expand = [None] * bin_areas.ndim
             expand[0] = slice(None)
-            bin_areas = bin_areas * n_sources[expand]
+            bin_areas = bin_areas * n_sources[tuple(expand)]
 
         # dimensions of the keys in expectations are now energy, radial bin
         if is_zenith_weight(zenith_index, self._aeff):