SWY: Monthly quickflow implementation is unclear and unmaintainable

[emlys]

Unsure if there's a bug here or not - but either way, we should refactor and document the _calculate_monthly_quick_flow function. It's not maintainable if we don't understand what it's doing.

• Is the original equation correct? Is it expected to produce these very large and very small numbers?
• If so, can we preemptively threshold the inputs to prevent over/underflow?
• Document and refactor the code so that it will be maintainable going forward.

QF details

According to the user's guide, QF is defined in terms of two cases:

1. For stream pixels: $\text{QF}_{stream,m} = \ P_{stream,m}$

2. For non-stream pixels:
   $\text{QF}_{i,m} = n_{m} \times \left( \left( a_{i,m} - S_{i} \right)\exp\left( - \frac{0.2S_{i}}{a_{i,m}} \right) + \frac{S_{i}^{2}}{a_{i,m}}\exp\left( \frac{0.8S_{i}}{a_{i,m}} \right)E_{1}\left( \frac{S_{i}}{a_{i,m}} \right) \right) \times \left( 25.4\ \left\lbrack \frac{\text{mm}}{\text{in}} \right\rbrack \right)$

This equation can't be directly implemented in code because some of the terms overflow or underflow on valid data.

The edge cases happen in this term of the equation:
$\frac{S_{i}^{2}}{a_{i,m}}\exp\left( \frac{0.8S_{i}}{a_{i,m}} \right)E_{1}\left( \frac{S_{i}}{a_{i,m}} \right)$

a. Where $s_i = 0$, this term should equal 0 (per conversation with Rafa). But, evaluating in numpy you get NaN and a warning because $E_{1}\left( 0 \right) = \inf$.

Solution: Preemptively set the result to 0 where s_i = 0 in order to avoid calculations with infinity.

b. When the ratio $\frac{S_{i}}{a_{i,m}}$ becomes large, this term approaches 0. [NOTE: I don't know how to prove this mathematically, but it holds true when I tested with reasonable values of s_i and a_im]. The $exp$ term becomes very large, while the $E1$ term becomes very small. When $\frac{S_{i}}{a_{i,m}}$ > 880 ish, $\exp\left( \frac{0.8S_{i}}{a_{i,m}} \right)$ will overflow and round up to infinity. Meanwhile, $E_{1}\left( \frac{S_{i}}{a_{i,m}} \right)$ will become so small that it rounds down to 0.

Values of s_i / a_im that make exp overflow aren't common, but are totally possible with valid input data, for example:

CN = 30
P = 6 (millimeters of rain per month)
n = 10 (number of rain events per month)

s_i = 1000/CN - 10 = 23.33
a_im = P/(n * 25.4) = 8/(10 * 25.4) = 0.023

s_i / a_im = 23.33 / 0.023 = 1014

exp(0.8 * 1014) = exp(811) = 1.6 * 10^352, which overflows float64. And E1 underflows:
>>> scipy.special.exp1(811)
0.0

Solution: Catch the overflow warning, and set the term to 0 anywhere that overflow happens.

[dcdenu4]

From current source code:

invest/src/natcap/invest/seasonal_water_yield/seasonal_water_yield.py

Lines 1110 to 1122 in ae2756c

	# Precompute the last two terms in quickflow so we can handle a
	# numerical instability when s_i is large and/or a_im is small
	# on large valid_si/a_im this number will be zero and the latter
	# exponent will also be zero because of a divide by zero. rather than
	# raise that numerical warning, just handle it manually
	E1 = scipy.special.expn(1, valid_si / a_im)
	E1[valid_si == 0] = 0
	nonzero_e1_mask = E1 != 0
	exp_result = numpy.zeros(valid_si.shape)
	exp_result[nonzero_e1_mask] = numpy.exp(
	(0.8 * valid_si[nonzero_e1_mask]) / a_im[nonzero_e1_mask] +
	numpy.log(E1[nonzero_e1_mask]))

[emlys]

From conversation with Doug, Lisa, and Rafa today:

• In the overflow case, we agree that the whole QF equation approaches 0 and we can set it to 0

• Pick an implementation that makes sense - whatever cutoff is easy to implement. Reasoning is that in this case you have little precipitation and it's all being absorbed into the soil, so quickflow is effectively 0.

• If we get any negative QF values, set them to 0 - negative QF doesn't make any sense