use non_nan arrays for greater than or less than comparisons

[MichaelMedford]

These changes get the same results for comparisons, but avoid errors like

/opt/anaconda3/envs/astroconda/lib/python3.7/site-packages/astropy/table/column.py:984: RuntimeWarning: invalid value encountered in greater
  result = getattr(super(), op)(other)
/opt/anaconda3/envs/astroconda/lib/python3.7/site-packages/astropy/table/column.py:984: RuntimeWarning: invalid value encountered in greater_equal
  result = getattr(super(), op)(other)

[mwhosek]

Hey Michael, this mostly looks good, but I did have one worry about the change to _make_companions_table(). On this line:

star_systems_phase_non_nan = np.nan_to_num(star_systems['phase'])

You are changing phase values that are nan to 0. However, setting phase = 0 has meaning for the MIST isochrones; it means that the star is on the main sequence, I believe. Could we change the nan to something that doesn't have meaning, like -99? Then, when we do:

bad = np.where( (star_systems_phase_non_nan > 5) & (star_systems_phase_non_nan < 101) & (star_systems_phase_non_nan != 9))

We might add a condition to catch these nan stars too.

Also, do we know why some of the phases are coming up as nans?

[MichaelMedford]

Hey Michael, this mostly looks good, but I did have one worry about the change to _make_companions_table(). On this line:

star_systems_phase_non_nan = np.nan_to_num(star_systems['phase'])

You are changing phase values that are nan to 0. However, setting phase = 0 has meaning for the MIST isochrones; it means that the star is on the main sequence, I believe. Could we change the nan to something that doesn't have meaning, like -99? Then, when we do:

bad = np.where( (star_systems_phase_non_nan > 5) & (star_systems_phase_non_nan < 101) & (star_systems_phase_non_nan != 9))

We might add a condition to catch these nan stars too.

Also, do we know why some of the phases are coming up as nans?

If you would feel more comfortable changing the nan value to -99 and creating a new condition, I'm all for that. My goal was to simply create the exact same behavior without generating warnings. From above:

In [1]: import numpy as np
In [2]: array = np.array([0, 50, np.nan])
In [3]: bad = np.where( (array > 5) & (array < 100) & (array != 9) )
/Users/michael/miniconda3/bin/ipython:1: RuntimeWarning: invalid value encountered in greater
  #!/Users/michael/miniconda3/bin/python
/Users/michael/miniconda3/bin/ipython:1: RuntimeWarning: invalid value encountered in less
  #!/Users/michael/miniconda3/bin/python
In [4]: print(bad)
(array([1]),)
In [5]: array_non_nan = np.nan_to_num(array)
In [6]: bad = np.where( (array_non_nan > 5) & (array_non_nan < 100) & (array_non_nan != 9) )
In [7]: print(bad)
(array([1]),)

However we could change it to the following:

In [9]: array_non_nan = np.nan_to_num(array, nan=-99)
In [10]: bad = np.where( (array_non_nan > 5) & (array_non_nan < 100) & (array_non_nan != 9) & (array_non_nan != -99) )
In [11]: print(bad)
(array([1]),)

I'm not sure why there are nan values floating around, but my guess would be that when you ran star_systems['phase'] = np.round(self.iso_interps['phase'](star_systems['mass'])), there were masses that were outside of the interpolation. So, for example:

In [18]: x = np.linspace(0, 1, 5)
In [19]: y = 2*x + 3
In [20]: interp = interpolate.interp1d(x, y, kind='linear', bounds_error=False, fill_value=np.nan)
In [21]: print(interp(0.5))
4.0
In [22]: print(interp(2))
nan

This currently goes hidden because you have set bounds_error=False. If instead it is set to True, you get:

In [23]: interp = interpolate.interp1d(x, y, kind='linear', bounds_error=True, fill_value=np.nan)
In [24]: print(interp(2))
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-24-78b73f63d070> in <module>()
----> 1 print(interp(2))

~/miniconda3/lib/python3.6/site-packages/scipy/interpolate/polyint.py in __call__(self, x)
     77         """
     78         x, x_shape = self._prepare_x(x)
---> 79         y = self._evaluate(x)
     80         return self._finish_y(y, x_shape)
     81 

~/miniconda3/lib/python3.6/site-packages/scipy/interpolate/interpolate.py in _evaluate(self, x_new)
    661         y_new = self._call(self, x_new)
    662         if not self._extrapolate:
--> 663             below_bounds, above_bounds = self._check_bounds(x_new)
    664             if len(y_new) > 0:
    665                 # Note fill_value must be broadcast up to the proper size

~/miniconda3/lib/python3.6/site-packages/scipy/interpolate/interpolate.py in _check_bounds(self, x_new)
    693                              "range.")
    694         if self.bounds_error and above_bounds.any():
--> 695             raise ValueError("A value in x_new is above the interpolation "
    696                              "range.")
    697 

ValueError: A value in x_new is above the interpolation range.

[mwhosek]

If you would feel more comfortable changing the nan value to -99 and creating a new condition, I'm all for that. My goal was to simply create the exact same behavior without generating warnings.

Yes, let's go with the change you suggest. I want to protect the "phase = 0" assignment.

I'm not sure why there are nan values floating around, but my guess would be that when you ran star_systems['phase'] = np.round(self.iso_interps'phase'), there were masses that were outside of the interpolation.

Right, now I remember this is why we used bounds_error=False to begin with. If a star is outside the mass range of the original isochrone, then we don't have any evolution info for it and what to remove it anyway. So, let's keep bounds_error=False and apply your change.

[MichaelMedford]

All instances of np.non_to_nan have been set to have nan=-99, and following comparisons adjusted to work for this addition.