Repository: PlasmaPy/PlasmaPy
Pull request #480 (2018-05-30 01:04:31)

Handling nparray Quantities part deux: Coulomb_logarithm
Chipping away at issue #469 , this time working on Coulomb_logarithm and the functions it relies on, namely impact_parameter and _boilerplate. 

As before, the goal this time is to make it so that Coulomb_logarithm works with len>1 array Quantities, such as `T_e = np.linspace(1,10,10) * u.eV`.  Let's see how this works.

One thought I had while working on this is the following.  So far I've found two ways to handle both scalar and vector Quantity inputs that need incompatible syntax:
- way 1: 
```
if np.isscalar(T_e.value):
    # handle scalar case
else: 
    # handle array case
```

vs. 
- way 2:
```
try:
    # handle scalar case
except ValueError: 
    # handle array case
```

I've heard a proverb that the "ask forgiveness not permission" is more pythonic, so maybe way 2 is better?  It makes me a little uneasy... maybe just takes some getting used to.  On way 1, it does feels a little clunky to have to get the `.value` for `np.isscalar`, but it doesn't seem to work correctly if you give it the full Quantity object.  Anybody have a preference or know a better / more preferred way?

Let's see how the tests are doing.

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



 - In the file [`plasmapy/physics/transport/tests/test_collisions.py`](https://github.com/PlasmaPy/PlasmaPy/blob/314427e26c8dc80feb95b97ccf0ce7fab28313a6/plasmapy/physics/transport/tests/test_collisions.py), following are the PEP8 issues :

> [Line 43:1](https://github.com/PlasmaPy/PlasmaPy/blob/314427e26c8dc80feb95b97ccf0ce7fab28313a6/plasmapy/physics/transport/tests/test_collisions.py#L43): [W293](https://duckduckgo.com/?q=pep8%20W293) blank line contains whitespace
> [Line 62:1](https://github.com/PlasmaPy/PlasmaPy/blob/314427e26c8dc80feb95b97ccf0ce7fab28313a6/plasmapy/physics/transport/tests/test_collisions.py#L62): [W293](https://duckduckgo.com/?q=pep8%20W293) blank line contains whitespace
> [Line 89:1](https://github.com/PlasmaPy/PlasmaPy/blob/314427e26c8dc80feb95b97ccf0ce7fab28313a6/plasmapy/physics/transport/tests/test_collisions.py#L89): [W293](https://duckduckgo.com/?q=pep8%20W293) blank line contains whitespace


StanczakDominik: I remember reading something about how one of these ways can be slower than the other... we might be able to figure something out here via `%%timeit`.
codecov[bot]: # [Codecov](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/480?src=pr&el=h1) Report
> :exclamation: No coverage uploaded for pull request base (`master@f944ae2`). [Click here to learn what that means](https://docs.codecov.io/docs/error-reference#section-missing-base-commit).
> The diff coverage is `98.14%`.

[![Impacted file tree graph](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/480/graphs/tree.svg?width=650&height=150&src=pr&token=LzZ5mnXVad)](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/480?src=pr&el=tree)

```diff
@@           Coverage Diff            @@
##             master    #480   +/-   ##
========================================
  Coverage          ?   92.9%           
========================================
  Files             ?      67           
  Lines             ?    6480           
  Branches          ?       0           
========================================
  Hits              ?    6020           
  Misses            ?     460           
  Partials          ?       0
```


| [Impacted Files](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/480?src=pr&el=tree) | Coverage Δ | |
|---|---|---|
| [plasmapy/physics/transport/braginskii.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/480/diff?src=pr&el=tree#diff-cGxhc21hcHkvcGh5c2ljcy90cmFuc3BvcnQvYnJhZ2luc2tpaS5weQ==) | `99.71% <ø> (ø)` | |
| [...lasmapy/physics/transport/tests/test\_collisions.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/480/diff?src=pr&el=tree#diff-cGxhc21hcHkvcGh5c2ljcy90cmFuc3BvcnQvdGVzdHMvdGVzdF9jb2xsaXNpb25zLnB5) | `100% <100%> (ø)` | |
| [plasmapy/physics/transport/collisions.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/480/diff?src=pr&el=tree#diff-cGxhc21hcHkvcGh5c2ljcy90cmFuc3BvcnQvY29sbGlzaW9ucy5weQ==) | `98.61% <96.15%> (ø)` | |

------

[Continue to review full report at Codecov](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/480?src=pr&el=continue).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/480?src=pr&el=footer). Last update [f944ae2...e536a00](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/480?src=pr&el=lastupdated). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments).

samurai688: Alright, I think this addresses comments so far.  It didn't seem to make any difference in Coulomb_logarithm performance to extend bmin to be always the same length as bmax when needed, so I put that in, I agree as @lemmatum said that's good behavior from the function.

@StanczakDominik  .. btw our rapidly shifting network of test coverage plugins seems to be behaving oddly with this PR.  Coveralls is suddenly mad at me about langmuir.py, which I haven't touched. Is it a rebase-needed thing? Though I guess it's better than Codecov, which gives me a cubist representation of unhappiness, but shuts down attempts to learn more due to lack of a "base report".  :)
StanczakDominik: Oh yeah, that's definitely a rebase-needed issue! I'll rebase later today and I think it should clear itself up. Don't worry about it! Thanks for making my day with your vivid description of that plot, though :D

Patch:
diff --git a/plasmapy/physics/transport/collisions.py b/plasmapy/physics/transport/collisions.py
index 2f33222d..5b9a03be 100644
--- a/plasmapy/physics/transport/collisions.py
+++ b/plasmapy/physics/transport/collisions.py
@@ -268,18 +268,21 @@ def Coulomb_logarithm(T,
         ln_Lambda = np.log(bmax / bmin)
     elif method == "GMS-3":
         ln_Lambda = np.log(bmax / bmin)
-        if ln_Lambda < 2:
-            ln_Lambda = 2 * u.dimensionless_unscaled
+        if np.any(ln_Lambda < 2):
+            if np.isscalar(ln_Lambda.value):
+                ln_Lambda = 2 * u.dimensionless_unscaled
+            else:
+                ln_Lambda[ln_Lambda < 2] = 2 * u.dimensionless_unscaled
     elif method in ["GMS-4", "GMS-5", "GMS-6"]:
         ln_Lambda = 0.5 * np.log(1 + bmax ** 2 / bmin ** 2)
     else:
         raise ValueError("Unknown method! Choose from 'classical' and 'GMS-N', N from 1 to 6.")
     # applying dimensionless units
     ln_Lambda = ln_Lambda.to(u.dimensionless_unscaled).value
-    if ln_Lambda < 2 and method in ["classical", "GMS-1", "GMS-2"]:
+    if np.any(ln_Lambda < 2) and method in ["classical", "GMS-1", "GMS-2"]:
         warnings.warn(f"Coulomb logarithm is {ln_Lambda} and {method} relies on weak coupling.",
                       utils.CouplingWarning)
-    elif ln_Lambda < 4:
+    elif np.any(ln_Lambda < 4):
         warnings.warn(f"Coulomb logarithm is {ln_Lambda}, you might have strong coupling effects",
                       utils.CouplingWarning)
     return ln_Lambda
@@ -306,11 +309,22 @@ def _boilerPlate(T, particles, V):
     # obtaining reduced mass of 2 particle collision system
     reduced_mass = atomic.reduced_mass(*particles)
 
-    if V == 0:
+    if np.any(V == 0):
         raise utils.exceptions.PhysicsError("You cannot have a collision for zero velocity!")
     # getting thermal velocity of system if no velocity is given
-    if np.isnan(V):
+    if V is None:
         V = parameters.thermal_speed(T, mass=reduced_mass)
+    elif np.any(np.isnan(V)):
+        if np.isscalar(V.value) and np.isscalar(T.value):
+            V = parameters.thermal_speed(T, mass=reduced_mass)
+        elif np.isscalar(V.value):
+            V = parameters.thermal_speed(T, mass=reduced_mass)
+        elif np.isscalar(T.value):
+            V = V.copy()
+            V[np.isnan(V)] = parameters.thermal_speed(T, mass=reduced_mass)
+        else:
+            V = V.copy()
+            V[np.isnan(V)] = parameters.thermal_speed(T[np.isnan(V)], mass=reduced_mass)
     _check_relativistic(V, 'V')
     return T, masses, charges, reduced_mass, V
 
@@ -541,10 +555,19 @@ def impact_parameter(T,
         # shorter than either of these two impact parameters, so we choose
         # the larger of these two possibilities. That is, between the
         # deBroglie wavelength and the distance of closest approach.
-        if bPerp > lambdaBroglie:
-            bmin = bPerp
-        else:
+        # ARRAY NOTES
+        # T and V should be guaranteed to be same size inputs from _boilerplate
+        # therefore, lambdaBroglie and bPerp are either both scalar or both array
+        # if np.isscalar(bPerp.value) and np.isscalar(lambdaBroglie.value):  # both scalar
+        try:  # assume both scalar
+            if bPerp > lambdaBroglie:
+                bmin = bPerp
+            else:
+                bmin = lambdaBroglie
+        # else:  # both lambdaBroglie and bPerp are arrays
+        except ValueError:  # both lambdaBroglie and bPerp are arrays
             bmin = lambdaBroglie
+            bmin[bPerp > lambdaBroglie] = bPerp[bPerp > lambdaBroglie]
     elif method == "GMS-1":
         # 1st method listed in Table 1 of reference [1]
         # This is just another form of the classical Landau-Spitzer
@@ -592,6 +615,13 @@ def impact_parameter(T,
         bmin = (lambdaBroglie ** 2 + bPerp ** 2) ** (1 / 2)
     else:
         raise ValueError(f"Method {method} not found!")
+    # ARRAY NOTES
+    # it could be that bmin and bmax have different sizes. If Te is a scalar,
+    # T and V will be scalar from _boilerplate, so bmin will scalar.  However
+    # if n_e is an array, than bmax will be an array. if this is the case,
+    # do we want to extend the scalar bmin to equal the length of bmax? Sure.
+    if np.isscalar(bmin.value) and not np.isscalar(bmax.value):
+        bmin = np.repeat(bmin, len(bmax))
     return bmin.to(u.m), bmax.to(u.m)
 
 
@@ -933,7 +963,7 @@ def fundamental_electron_collision_freq(T_e,
                                  n_e,
                                  particles,
                                  z_mean=Z_i,
-                                 V=np.nan, # probably needs to be enabled!
+                                 V=np.nan * u.m / u.s,
                                  method=coulomb_log_method)
         # dividing out by typical Coulomb logarithm value implicit in
         # the collision frequency calculation and replacing with


Test patch:
diff --git a/plasmapy/physics/transport/tests/test_collisions.py b/plasmapy/physics/transport/tests/test_collisions.py
index ba7031fb..2fa7fc9e 100644
--- a/plasmapy/physics/transport/tests/test_collisions.py
+++ b/plasmapy/physics/transport/tests/test_collisions.py
@@ -1,6 +1,7 @@
 import numpy as np
 import pytest
 from astropy import units as u
+from astropy.tests.helper import assert_quantity_allclose
 
 from plasmapy.physics.transport import (Coulomb_logarithm,
                                         impact_parameter_perp,
@@ -20,7 +21,9 @@ class Test_Coulomb_logarithm:
     def setup_class(self):
         """initializing parameters for tests """
         self.temperature1 = 10 * 11604 * u.K
+        self.T_arr = np.array([1, 2]) * u.eV
         self.density1 = 1e20 * u.cm ** -3
+        self.n_arr = np.array([1e20, 2e20]) * u.cm ** -3
         self.temperature2 = 1 * 11604 * u.K
         self.density2 = 1e23 * u.cm ** -3
         self.z_mean = 2.5
@@ -38,12 +41,77 @@ def setup_class(self):
         self.gms6 = 3.635342040477818
         self.gms6_negative = 0.030720859361047514
 
+    def test_handle_invalid_V(self):
+        """Test that V default, V = None, and V = np.nan all give the same result"""
+        methodVal_0 = Coulomb_logarithm(self.T_arr[0],
+                                        self.n_arr[0],
+                                        self.particles,
+                                        z_mean=1 * u.dimensionless_unscaled,
+                                        V=np.nan * u.m / u.s)
+        methodVal_1 = Coulomb_logarithm(self.T_arr[0],
+                                        self.n_arr[0],
+                                        self.particles,
+                                        z_mean=1 * u.dimensionless_unscaled,
+                                        V=None)
+        methodVal_2 = Coulomb_logarithm(self.T_arr[0],
+                                        self.n_arr[0],
+                                        self.particles,
+                                        z_mean=1 * u.dimensionless_unscaled)
+        assert_quantity_allclose(methodVal_0, methodVal_1)
+        assert_quantity_allclose(methodVal_0, methodVal_2)
+
+    def test_handle_zero_V(self):
+        """Test that V == 0 returns a PhysicsError"""
+        with pytest.raises(exceptions.PhysicsError):
+            Coulomb_logarithm(self.T_arr[0],
+                              self.n_arr[0],
+                              self.particles,
+                              z_mean=1 * u.dimensionless_unscaled,
+                              V=0 * u.m / u.s)
+
+    @pytest.mark.parametrize("method", ["classical", "GMS-1", "GMS-2", "GMS-3", "GMS-4", "GMS-5",
+                                        "GMS-6"])
+    def test_handle_numpy_array(self, method):
+        """Tests to verify that can handle Quantities with numpy array as the value"""
+        methodVal = Coulomb_logarithm(self.T_arr,
+                                      self.n_arr,
+                                      self.particles,
+                                      z_mean=1 * u.dimensionless_unscaled,
+                                      V=np.nan * u.m / u.s,
+                                      method=method)
+        methodVal_0 = Coulomb_logarithm(self.T_arr[0],
+                                        self.n_arr[0],
+                                        self.particles,
+                                        z_mean=1 * u.dimensionless_unscaled,
+                                        V=np.nan * u.m / u.s,
+                                        method=method)
+        assert_quantity_allclose(methodVal[0], methodVal_0)
+
+    def test_handle_V_arraysizes(self):
+        """Test that different sized V input array gets handled by _boilerplate"""
+        methodVal_0 = Coulomb_logarithm(self.T_arr[0],
+                                        self.n_arr[0],
+                                        self.particles,
+                                        z_mean=1 * u.dimensionless_unscaled,
+                                        V=np.array([np.nan, 3e7]) * u.m / u.s)
+        methodVal_1 = Coulomb_logarithm(self.T_arr[1],
+                                        self.n_arr[0],
+                                        self.particles,
+                                        z_mean=1 * u.dimensionless_unscaled,
+                                        V=np.array([1e7, np.nan]) * u.m / u.s)
+        methodVal_2 = Coulomb_logarithm(self.T_arr,
+                                        self.n_arr[0],
+                                        self.particles,
+                                        z_mean=1 * u.dimensionless_unscaled,
+                                        V=np.array([np.nan, np.nan]) * u.m / u.s)
+        assert_quantity_allclose(methodVal_0[0], methodVal_2[0])
+        assert_quantity_allclose(methodVal_1[1], methodVal_2[1])
+
     def test_symmetry(self):
         lnLambda = Coulomb_logarithm(self.temperature1, self.density2, self.particles)
         lnLambdaRev = Coulomb_logarithm(self.temperature1, self.density2, self.particles[::-1])
         assert lnLambda == lnLambdaRev
 
-
     def test_Chen_Q_machine(self):
         """
         Tests whether Coulomb logarithm gives value consistent with
@@ -66,7 +134,6 @@ def test_Chen_Q_machine(self):
                   f"{lnLambdaChen} and not {lnLambda}.")
         assert testTrue, errStr
 
-
     def test_Chen_lab(self):
         """
         Tests whether Coulomb logarithm gives value consistent with
@@ -524,7 +591,9 @@ class Test_impact_parameter:
     def setup_class(self):
         """initializing parameters for tests """
         self.T = 11604 * u.K
+        self.T_arr = np.array([1, 2]) * u.eV
         self.n_e = 1e17 * u.cm ** -3
+        self.n_e_arr = np.array([1e17, 2e17]) * u.cm ** -3
         self.particles = ('e', 'p')
         self.z_mean = 2.5
         self.V = 1e4 * u.km / u.s
@@ -587,6 +656,24 @@ def test_bad_method(self):
                              V=np.nan * u.m / u.s,
                              method="meow")
 
+    @pytest.mark.parametrize("method", ["classical", "GMS-1", "GMS-2", "GMS-3", "GMS-4", "GMS-5",
+                                        "GMS-6"])
+    def test_handle_numpy_array(self, method):
+        """Tests to verify that can handle Quantities with numpy array as the value"""
+        methodVal = impact_parameter(self.T_arr,
+                                     self.n_e_arr,
+                                     self.particles,
+                                     z_mean=1 * u.dimensionless_unscaled,
+                                     V=np.nan * u.m / u.s,
+                                     method=method)
+        methodVal_0 = impact_parameter(self.T_arr[0],
+                                       self.n_e_arr[0],
+                                       self.particles,
+                                       z_mean=1 * u.dimensionless_unscaled,
+                                       V=np.nan * u.m / u.s,
+                                       method=method)
+        assert_quantity_allclose((methodVal[0][0], methodVal[1][0]), methodVal_0)
+
 
 class Test_collision_frequency:
     @classmethod