Repository: PlasmaPy/PlasmaPy
Pull request #477 (2018-05-26 16:32:28)

Modify gyroradius function to work with numpy array Quantities
First baby step of progress towards issue #469, this one is focused solely on the gyroradius function.

Gyroradius seemed to have trouble due to its unique feature, in which it will accept a T_i input, or a Vperp input, but not both.  I put in some conditional statements to hopefully extend this behavior to numpy array Quantities.  

The main new behavior that I have tried to test for and achieve is:
- If only Ti is input and is a no-nan-Qarray, function should work (result is a Qarray calculated with Ti)
- If only Vperp is input and is a no-nan-Qarray, function should work (result is a Qarray calculated with Vperp)
- If both Vperp and Ti are input and are both Qarray, but with nans in each array such that for each element there is exactly 1 nan, function should work (result is a Qarray calculated with the non-nans)
- if, for some reason, Vperp is a Qarray of all nan and Ti is valid, function should work. The result I went with was to return a Qarray of len(Vperp) all calculated with the valid scalar Ti, but I'm open to other returns here. Hopefully this case won't get activated much.
- same as above but vice versa with Vperp, Ti

And hopefully raise ValueError on all the various invalid combinations. Of note is
- if Ti is a valid scalar, and Vperp is a Qarray that has some nan but not all, raise ValueError
- same as above but vice versa with Vperp, Ti

This one is debateable, I think doing it this way of raising the ValueError keeps to the spirit of the original "you must provide either Ti or Vperp but not both" and is hopefully more consistent and least confusing, but I could also imagine a way in which the scalar would be used only in the cases where there are the nans in the array.

---

Actually, it looks like gyroradius and hall parameter are the only two plasma parameter functions that currently have issues when given numpy array Quantities.  So, we'd be 50% done!...  Except that hall parameter calls collision_frequency, which opens up the collisions module, which looks to be in a world of pain.



pep8speaks: Hello @samurai688! Thanks for submitting your pull request.

 - In the file [`plasmapy/physics/parameters.py`](https://github.com/PlasmaPy/PlasmaPy/blob/36348688cc883fcfe11ac7c94bbe65637153a548/plasmapy/physics/parameters.py), following are the PEP8 issues :

> [Line 876:1](https://github.com/PlasmaPy/PlasmaPy/blob/36348688cc883fcfe11ac7c94bbe65637153a548/plasmapy/physics/parameters.py#L876): [W293](https://duckduckgo.com/?q=pep8%20W293) blank line contains whitespace
> [Line 884:1](https://github.com/PlasmaPy/PlasmaPy/blob/36348688cc883fcfe11ac7c94bbe65637153a548/plasmapy/physics/parameters.py#L884): [W293](https://duckduckgo.com/?q=pep8%20W293) blank line contains whitespace
> [Line 886:5](https://github.com/PlasmaPy/PlasmaPy/blob/36348688cc883fcfe11ac7c94bbe65637153a548/plasmapy/physics/parameters.py#L886): [E303](https://duckduckgo.com/?q=pep8%20E303) too many blank lines (2)

 - In the file [`plasmapy/physics/tests/test_parameters.py`](https://github.com/PlasmaPy/PlasmaPy/blob/36348688cc883fcfe11ac7c94bbe65637153a548/plasmapy/physics/tests/test_parameters.py), following are the PEP8 issues :

> [Line 601:1](https://github.com/PlasmaPy/PlasmaPy/blob/36348688cc883fcfe11ac7c94bbe65637153a548/plasmapy/physics/tests/test_parameters.py#L601): [W293](https://duckduckgo.com/?q=pep8%20W293) blank line contains whitespace
> [Line 605:9](https://github.com/PlasmaPy/PlasmaPy/blob/36348688cc883fcfe11ac7c94bbe65637153a548/plasmapy/physics/tests/test_parameters.py#L605): [E128](https://duckduckgo.com/?q=pep8%20E128) continuation line under-indented for visual indent
> [Line 606:1](https://github.com/PlasmaPy/PlasmaPy/blob/36348688cc883fcfe11ac7c94bbe65637153a548/plasmapy/physics/tests/test_parameters.py#L606): [W293](https://duckduckgo.com/?q=pep8%20W293) blank line contains whitespace
> [Line 612:1](https://github.com/PlasmaPy/PlasmaPy/blob/36348688cc883fcfe11ac7c94bbe65637153a548/plasmapy/physics/tests/test_parameters.py#L612): [W293](https://duckduckgo.com/?q=pep8%20W293) blank line contains whitespace
> [Line 618:1](https://github.com/PlasmaPy/PlasmaPy/blob/36348688cc883fcfe11ac7c94bbe65637153a548/plasmapy/physics/tests/test_parameters.py#L618): [W293](https://duckduckgo.com/?q=pep8%20W293) blank line contains whitespace

PlasmaPy follows the [PEP8 style guide](https://www.python.org/dev/peps/pep-0008/), which may be checked locally using [pycodestyle](http://pycodestyle.pycqa.org/en/latest/intro.html)

Patch:
diff --git a/plasmapy/physics/parameters.py b/plasmapy/physics/parameters.py
index 7914740e..f01bea4e 100644
--- a/plasmapy/physics/parameters.py
+++ b/plasmapy/physics/parameters.py
@@ -874,12 +874,45 @@ def gyroradius(B, particle='e-', *, Vperp=np.nan * u.m / u.s, T_i=np.nan * u.K):
 
     """
 
-    if np.isfinite(Vperp) and np.isfinite(T_i):
-        raise ValueError("Cannot have both Vperp and T_i as arguments to "
-                         "gyroradius")
-
-    elif not np.isfinite(Vperp) and np.isfinite(T_i):
-        Vperp = thermal_speed(T_i, particle=particle)
+    isfinite_Ti = np.isfinite(T_i)
+    isfinite_Vperp = np.isfinite(Vperp)
+
+    # check 1: ensure either Vperp or T_i invalid, keeping in mind that
+    # the underlying values of the astropy quantity may be numpy arrays
+    if np.any(np.logical_not(np.logical_xor(isfinite_Vperp, isfinite_Ti))):
+        raise ValueError("Must give Vperp or T_i, but not both, as arguments to gyroradius")
+
+    # check 2: get Vperp as the thermal speed if is not already a valid input
+    if np.isscalar(Vperp.value) and np.isscalar(T_i.value):  # both T_i and Vperp are scalars
+        # we know exactly one of them is nan from check 1
+        if isfinite_Ti:
+            # T_i is valid, so use it to determine Vperp
+            Vperp = thermal_speed(T_i, particle=particle)
+        # else: Vperp is alread valid, do nothing
+    elif np.isscalar(Vperp.value):  # only T_i is an array
+        # this means either Vperp must be nan, or T_i must be array of all nan,
+        # or else we couldn't have gotten through check 1
+        if isfinite_Vperp:
+            # Vperp is valid, T_i is a vector that is all nan
+            # uh...
+            Vperp = np.repeat(Vperp, len(T_i))
+        else:
+            # normal case where Vperp is scalar nan and T_i is valid array
+            Vperp = thermal_speed(T_i, particle=particle)
+    elif np.isscalar(T_i.value):  # only Vperp is an array
+        # this means either T_i must be nan, or V_perp must be array of all nan,
+        # or else we couldn't have gotten through check 1
+        if isfinite_Ti:
+            # T_i is valid, V_perp is an array of all nan
+            # uh...
+            Vperp = thermal_speed(np.repeat(T_i, len(Vperp)), particle=particle)
+        # else: normal case where T_i is scalar nan and Vperp is already a valid array
+        # so, do nothing
+    else:  # both T_i and Vperp are arrays
+        # we know all the elementwise combinations have one nan and one finite, due to check 1
+        # use the valid Vperps, and replace the others with those calculated from T_i
+        Vperp = Vperp.copy()  # avoid changing Vperp's value outside function
+        Vperp[isfinite_Ti] = thermal_speed(T_i[isfinite_Ti], particle=particle)
 
     omega_ci = gyrofrequency(B, particle)
 


Test patch:
diff --git a/plasmapy/physics/tests/test_parameters.py b/plasmapy/physics/tests/test_parameters.py
index ff7e5cd6..b9974119 100644
--- a/plasmapy/physics/tests/test_parameters.py
+++ b/plasmapy/physics/tests/test_parameters.py
@@ -38,17 +38,25 @@
 T_i = 1e6 * u.K
 
 B_arr = np.array([0.001, 0.002]) * u.T
-rho_arr = np.array([5e-10, 2e-10]) * u.kg / u.m ** 3
-
 B_nanarr = np.array([0.001, np.nan]) * u.T
+B_allnanarr = np.array([np.nan, np.nan]) * u.T
+
+rho_arr = np.array([5e-10, 2e-10]) * u.kg / u.m ** 3
 rho_infarr = np.array([np.inf, 5e19]) * u.m ** -3
 rho_negarr = np.array([-5e19, 6e19]) * u.m ** -3
+
+T_arr = np.array([1e6, 2e6]) * u.K
 T_nanarr = np.array([1e6, np.nan]) * u.K
+T_nanarr2 = np.array([np.nan, 2e6]) * u.K
+T_allnanarr = np.array([np.nan, np.nan]) * u.K
 T_negarr = np.array([1e6, -5151.]) * u.K
 
-mu = m_p.to(u.u).value
-
 V = 25.2 * u.m / u.s
+V_arr = np.array([25, 50]) * u.m / u.s
+V_nanarr = np.array([25, np.nan]) * u.m / u.s
+V_allnanarr = np.array([np.nan, np.nan]) * u.m / u.s
+
+mu = m_p.to(u.u).value
 
 
 class Test_mass_density:
@@ -588,6 +596,48 @@ def test_gyroradius():
         gyroradius(1.1 * u.T, particle="p", Vperp=1.2 * u.m, T_i=1.1 * u.K)
 
 
+class Test_gyroradius:
+
+    def test_handle_numpy_array(self):
+        # Tests to verify that can handle Quantities with numpy array as the value:
+        assert gyroradius(B_arr, Vperp=V_arr)[0] == gyroradius(B_arr[0], Vperp=V_arr[0])
+        assert gyroradius(B_arr, T_i=T_arr)[0] == gyroradius(B_arr[0], T_i=T_arr[0])
+
+    def test_handle_mixed_Qarrays(self):
+        # If both Vperp or Ti are input as Qarrays, but only one of the two is valid
+        # at each element, then that's fine, the function should work:
+        assert gyroradius(B_arr, Vperp=V_nanarr, T_i=T_nanarr2)[0] == gyroradius(B_arr[0],
+                                                                                 Vperp=V_nanarr[0],
+                                                                                 T_i=T_nanarr2[0])
+
+    def test_raise_two_valid_inputs(self):
+        # If both Vperp or Ti are nan-less, Qarrays or not, should raise ValueError:
+        with pytest.raises(ValueError):
+            gyroradius(B_arr, Vperp=V, T_i=T_arr)
+        with pytest.raises(ValueError):
+            gyroradius(B_arr, Vperp=V_arr, T_i=T_i)
+
+    def test_all_valid_and_one_valid(self):
+        # If one of (Vperp, Ti) is a valid and one is Qarray with at least one valid, ValueError:
+        with pytest.raises(ValueError):
+            gyroradius(B_arr, Vperp=V, T_i=T_nanarr)
+        with pytest.raises(ValueError):
+            gyroradius(B_arr, Vperp=V_nanarr, T_i=T_i)
+
+    def test_scalar_and_nan_qarray(self):
+        # If either Vperp or Ti is a valid scalar and the other is a Qarray of all nans,
+        # should do something valid and not raise a ValueError
+        assert np.all(np.isfinite(gyroradius(B_arr, Vperp=V, T_i=T_allnanarr)))
+        assert np.all(np.isfinite(gyroradius(B_arr, Vperp=V_allnanarr, T_i=T_i)))
+
+    def test_keeps_arguments_unchanged(self):
+        Vperp1 = u.Quantity([np.nan, 1], unit=u.m/u.s)
+        Vperp2 = u.Quantity([np.nan, 1], unit=u.m/u.s)  # an exact copy
+        T_i = u.Quantity([1, np.nan], unit=u.K)
+
+        gyroradius(B_arr, Vperp=Vperp1, T_i = T_i)
+        assert_quantity_allclose(Vperp1, Vperp2)
+
 def test_plasma_frequency():
     r"""Test the plasma_frequency function in parameters.py."""