routing: regularize bimodal model

If the success and fail amounts indicate that a channel doesn't obey a
bimodal distribution, we fall back to a uniform/linear success
probability model. This also helps to avoid numerical normalization
issues with the bimodal model.

This is achieved by adding a very small summand to the balance
distribution P(x) ~ exp(-x/s) + exp((x-c)/s), 1/c that helps to
regularize the probability distribution. The distribution becomes finite
for intermediate balances where the exponentials would be evaluated to
an exact zero (float) otherwise. This regularization is effective in
edge cases and leads to falling back to a uniform model should the
bimodal model fail.

This affects the normalization to be s * (-2 * exp(-c/s) + 2 + 1/s) and
the primitive function to receive an extra term x/(cs).

The previously added fuzz seed is expected to be resolved with this.

routing/probability_bimodal.go

@@ -416,34 +416,38 @@ func cannotSend(failAmount, capacity lnwire.MilliSatoshi, now,
 
 // primitive computes the indefinite integral of our assumed (normalized)
 // liquidity probability distribution. The distribution of liquidity x here is
-// the function P(x) ~ exp(-x/s) + exp((x-c)/s), i.e., two exponentials residing
+// the function P(x) ~ exp(-x/s) + exp((x-c)/s) + 1/c, i.e., two exponentials
-// at the ends of channels. This means that we expect liquidity to be at either
+// residing at the ends of channels. This means that we expect liquidity to be
-// side of the channel with capacity c. The s parameter (scale) defines how far
+// at either side of the channel with capacity c. The s parameter (scale)
-// the liquidity leaks into the channel. A very low scale assumes completely
+// defines how far the liquidity leaks into the channel. A very low scale
-// unbalanced channels, a very high scale assumes a random distribution. More
+// assumes completely unbalanced channels, a very high scale assumes a random
-// details can be found in
+// distribution. More details can be found in
 // https://github.com/lightningnetwork/lnd/issues/5988#issuecomment-1131234858.
+// Additionally, we add a constant term 1/c to the distribution to avoid
+// normalization issues and to fall back to a uniform distribution should the
+// previous success and fail amounts contradict a bimodal distribution.
 func (p *BimodalEstimator) primitive(c, x float64) float64 {
 	s := float64(p.BimodalScaleMsat)
 
 	// The indefinite integral of P(x) is given by
-	// Int P(x) dx = H(x) = s * (-e(-x/s) + e((x-c)/s)),
+	// Int P(x) dx = H(x) = s * (-e(-x/s) + e((x-c)/s) + x/(c*s)),
 	// and its norm from 0 to c can be computed from it,
-	// norm = [H(x)]_0^c = s * (-e(-c/s) + 1 -(1 + e(-c/s))).
+	// norm = [H(x)]_0^c = s * (-e(-c/s) + 1 + 1/s -(-1 + e(-c/s))) =
+	// = s * (-2*e(-c/s) + 2 + 1/s).
+	// The prefactors s are left out, as they cancel out in the end.
+	// norm can only become zero, if c is zero, which we sorted out before
+	// calling this method.
 	ecs := math.Exp(-c / s)
-	exs := math.Exp(-x / s)
+	norm := -2*ecs + 2 + 1/s
 
 	// It would be possible to split the next term and reuse the factors
 	// from before, but this can lead to numerical issues with large
 	// numbers.
 	excs := math.Exp((x - c) / s)
-
+	exs := math.Exp(-x / s)
-	// norm can only become zero, if c is zero, which we sorted out before
-	// calling this method.
-	norm := -2*ecs + 2
 
 	// We end up with the primitive function of the normalized P(x).
-	return (-exs + excs) / norm
+	return (-exs + excs + x/(c*s)) / norm
 }
 
 // integral computes the integral of our liquidity distribution from the lower

routing/probability_bimodal_test.go

@@ -269,10 +269,30 @@ func TestSmallScale(t *testing.T) {
 	// An amount that's close to the success amount should have a very high
 	// probability.
 	amtCloseSuccess := successAmount + 1
-	_, err := estimator.probabilityFormula(
+	p, err := estimator.probabilityFormula(
 		capacity, successAmount, failAmount, amtCloseSuccess,
 	)
-	require.ErrorContains(t, err, "normalization factor is zero")
+	require.NoError(t, err)
+	require.InDelta(t, 1.0, p, defaultTolerance)
+
+	// An amount that's close to the fail amount should have a very low
+	// probability.
+	amtCloseFail := failAmount - 1
+	p, err = estimator.probabilityFormula(
+		capacity, successAmount, failAmount, amtCloseFail,
+	)
+	require.NoError(t, err)
+	require.InDelta(t, 0.0, p, defaultTolerance)
+
+	// In the region where the bimodal model doesn't give good forecasts, we
+	// fall back to a uniform model, which interpolates probabilities
+	// linearly.
+	amtLinear := successAmount + (failAmount-successAmount)*1/4
+	p, err = estimator.probabilityFormula(
+		capacity, successAmount, failAmount, amtLinear,
+	)
+	require.NoError(t, err)
+	require.InDelta(t, 0.75, p, defaultTolerance)
 }
 
 // TestIntegral tests certain limits of the probability distribution integral.