[Another vectorised implementation that works with an intermediate boolean array of n_clusters x n_samples rather than n_samples x n_samples:
```python
cluster_assignments = np.repeat([np.arange(np.max(labels) + 1)], labels.size, axis=0).T == labels
C += cluster_assignments[labels]
```
This has the benefit of allowing you to first convert `cluster_assignments` to a sparse matrix which may be very worthwhile for `C` updating. But this is all premature optimisation on my part.]