Repository: tskit-dev/msprime (language: python)

Issue:
TPM microsats model transition matrix error
When specifying the following `MicrosatMutationModel` I get an "Each row of the transition matrix must be non-negative and sum to one" error.
```{python}
lo = 50
hi = 500
m = 0.25
p = 0.5
model = msprime.TPM(p=p, m=m, lo=lo, hi=hi)
```
I just followed the instructions from the docs to generate a random model. If this is not a valid parameter combination then this should have been caught I think. @andrewkern any idea what is going on here?
TPM microsats model transition matrix error
When specifying the following `MicrosatMutationModel` I get an "Each row of the transition matrix must be non-negative and sum to one" error.
```{python}
lo = 50
hi = 500
m = 0.25
p = 0.5
model = msprime.TPM(p=p, m=m, lo=lo, hi=hi)
```
I just followed the instructions from the docs to generate a random model. If this is not a valid parameter combination then this should have been caught I think. @andrewkern any idea what is going on here?


Discussion / hints:
looks like that value of `lo` is revealing a bug. if I change that value everything works. it's weird because I can't confirm that the rows don't sum to 1, e.g.

```
>>> x = msprime.mutations.general_microsat_rate_matrix(s=0, u=0.5, v=0, p=p, m=m, lo=50, hi=500)
>>> np.sum(x, axis=1)
array([1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1.])
   ```
my first blush reaction is that something wonky is happening with matrix_mutation_model class on the C side? not sure yet

a little further digging -- this looks like rounding error in the transition matrix construction

```
x = msprime.mutations.general_microsat_rate_matrix(s=0, u=0.5, v=0, p=p, m=m, lo=50, hi=500)
print(x[np.where(x<0)])

[-2.22044605e-16]
```

the negative element in the matrix is on the diagonal of the rate matrix, so the numerical issue is creeping in here most likely

(from [`msprime.mutations.general_microsat_rate_matrix`](https://github.com/tskit-dev/msprime/blob/2d00f51fc65fce036accc0102ffeae6f0f108ea7/msprime/mutations.py#L302))
```python
    row_sums = rate_matrix.sum(axis=1, dtype="float64")
    np.fill_diagonal(rate_matrix, 1.0 - row_sums)
```
@petrelharp any recs on how best to beat this? set small diagonal values to zero? 
Well, here's the problem, I think:
```
>>> 1 == (1 / 1.0000000000000001)
True
```
In the problematic code, the max row sum is `alpha=1.0000000000000002` *before* normalizing; then dividing by this value *does not change* the max row sum (i.e., it's still 1.0000000000000002).

Doing as you suggest and changing [that line](https://github.com/tskit-dev/msprime/blob/2d00f51fc65fce036accc0102ffeae6f0f108ea7/msprime/mutations.py#L302) to
```
    np.fill_diagonal(rate_matrix, np.fmax(0.0, 1.0 - row_sums))
```
fixes the error at least with this set of parameters. I'm not sure if there's a better way other than changing the code that throws the error to be more tolerant of a little slop?

I've got a PR, will throw it up.
looks like that value of `lo` is revealing a bug. if I change that value everything works. it's weird because I can't confirm that the rows don't sum to 1, e.g.

```
>>> x = msprime.mutations.general_microsat_rate_matrix(s=0, u=0.5, v=0, p=p, m=m, lo=50, hi=500)
>>> np.sum(x, axis=1)
array([1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1.])
   ```
my first blush reaction is that something wonky is happening with matrix_mutation_model class on the C side? not sure yet

a little further digging -- this looks like rounding error in the transition matrix construction

```
x = msprime.mutations.general_microsat_rate_matrix(s=0, u=0.5, v=0, p=p, m=m, lo=50, hi=500)
print(x[np.where(x<0)])

[-2.22044605e-16]
```

the negative element in the matrix is on the diagonal of the rate matrix, so the numerical issue is creeping in here most likely

(from [`msprime.mutations.general_microsat_rate_matrix`](https://github.com/tskit-dev/msprime/blob/2d00f51fc65fce036accc0102ffeae6f0f108ea7/msprime/mutations.py#L302))
```python
    row_sums = rate_matrix.sum(axis=1, dtype="float64")
    np.fill_diagonal(rate_matrix, 1.0 - row_sums)
```
@petrelharp any recs on how best to beat this? set small diagonal values to zero? 
Well, here's the problem, I think:
```
>>> 1 == (1 / 1.0000000000000001)
True
```
In the problematic code, the max row sum is `alpha=1.0000000000000002` *before* normalizing; then dividing by this value *does not change* the max row sum (i.e., it's still 1.0000000000000002).

Doing as you suggest and changing [that line](https://github.com/tskit-dev/msprime/blob/2d00f51fc65fce036accc0102ffeae6f0f108ea7/msprime/mutations.py#L302) to
```
    np.fill_diagonal(rate_matrix, np.fmax(0.0, 1.0 - row_sums))
```
fixes the error at least with this set of parameters. I'm not sure if there's a better way other than changing the code that throws the error to be more tolerant of a little slop?

I've got a PR, will throw it up.

Fix patch:
diff --git a/msprime/mutations.py b/msprime/mutations.py
index 01ee60dff..1db11890f 100644
--- a/msprime/mutations.py
+++ b/msprime/mutations.py
@@ -296,10 +296,11 @@ def gamma_microsat(m, i, j, lo, hi):
                 rate_matrix[ii, jj] = 0
     # scale by max row sum; then normalize to 1
     row_sums = rate_matrix.sum(axis=1, dtype="float64")
-    alpha = max(row_sums)
+    alpha = np.max(row_sums)
     rate_matrix /= alpha
     row_sums = rate_matrix.sum(axis=1, dtype="float64")
-    np.fill_diagonal(rate_matrix, 1.0 - row_sums)
+    # the max(0, *) is to avoid floating point error
+    np.fill_diagonal(rate_matrix, np.fmax(0.0, 1.0 - row_sums))
     return rate_matrix
 
 


Test patch:
diff --git a/tests/test_mutations.py b/tests/test_mutations.py
index 882646600..cd708624a 100644
--- a/tests/test_mutations.py
+++ b/tests/test_mutations.py
@@ -358,11 +358,14 @@ def test_SMM(self):
         self.validate_model(model)
         self.validate_stationary(model)
 
-    def test_TPM(self):
-        lo = 2
-        hi = 50
-        m = 0.9
-        p = 0.9
+    @pytest.mark.parametrize(
+        "p, m, lo, hi",
+        [
+            (0.9, 0.9, 2, 50),
+            (0.5, 0.25, 50, 500),
+        ],
+    )
+    def test_TPM(self, p, m, lo, hi):
         model = msprime.TPM(p=p, m=m, lo=lo, hi=hi)
         self.validate_model(model)
         self.validate_stationary(model)