Replacement for `BQM.to_numpy_matrix()`

`BQM.to_numpy_matrix()` was removed in dimod 0.10.0, however creating a dense matrix from a BQM is still of interest to some users. It can currently be done with
```python
import dimod
import numpy as np

Q = np.triu(np.random.default_rng(42).random((5, 5)))

bqm = dimod.BinaryQuadraticModel(Q, 'BINARY')

lin, (row, col, quad), offset = bqm.to_numpy_vectors(sort_indices=True)

newQ = np.zeros((5, 5))
np.fill_diagonal(newQ, lin)
newQ[row, col] = quad

np.testing.assert_array_equal(Q, newQ)
```
**Proposals**
The main barrier is handing the linear terms, because the expectations should be different for binary and spin bqms. I think this leads to a few options

1) Add a `BQM.to_dense()` method (as proposed in https://github.com/dwavesystems/dimod/issues/617) that returns a single matrix.
```python
Q = bqm.to_dense()
```
and raises an exception if the binary quadratic model is SPIN-valued. This would not have an equivalent implementation on quadratic models.

2) Add a `BQM.to_dense()` method that always returns the linear terms as an array
```python
lin, quad, offset = bqm.to_dense()
```
This has the advantage of being valid for both BINARY and SPIN, as well as for quadratic models. Folks who just want the QUBO matrix can then do
```python
quad.fill_diagonal(lin)
```

3) A slightly more sugary version would be implementing `__array__` on both `BQM.linear` and `BQM.quadratic`
```python
lin = np.asarray(bqm.linear)
quad = np.asarray(bqm.quadratic)
```
I like this because it's very explicit about which part of the model is being turned into an array. It also allows for some nice construction like `np.array(bqm.linear).sum()`