Birch is not finding clusters correctly

### Describe the bug

Birch do not correctly find clusters if the number of samples of one cluster is much larger than another cluster even if the distance between clusters is much larger than threshold.

### Steps/Code to Reproduce

```python
from sklearn.datasets import make_blobs
from sklearn.cluster import Birch
import numpy as np

X1, y1 = make_blobs(n_samples=10000, centers=[[0.6, 0.6]], cluster_std=0.01, random_state=0)
X2, y2 = make_blobs(n_samples=100, centers=[[0.1, 0.1]], cluster_std=0.01, random_state=0)

X = np.append(X1, X2, axis = 0)

brc = Birch(n_clusters=None, threshold = 0.1)
brc.fit(X)

brc.subcluster_centers_
```

### Expected Results

Expected two subclusters  since the distance between clusters is much larger than used threshold.

### Actual Results

Only one subcluster is detected

`brc.subcluster_centers_`:

```
array([[0.59494558, 0.59509375]])
```


### Versions

```shell
System:
    python: 3.8.6 | packaged by conda-forge | (default, Jan 25 2021, 23:21:18)  [GCC 9.3.0]
executable: /home/envs/env1_custom/bin/python
   machine: Linux-4.14.252

Python dependencies:
          pip: 22.0.3
   setuptools: 60.9.3
      sklearn: 1.0.2
        numpy: 1.20.2
        scipy: 1.7.0
       Cython: None
       pandas: 1.2.4
   matplotlib: 3.5.1
       joblib: 1.0.1
threadpoolctl: 2.1.0

Built with OpenMP: True
```

COMMENT:
I digged into the implementation. You know that the way BIRCH works is by adding the point to a node and see if this new node respects the validity condition, ie here if the radius of the new node is inferior to the threshold. As you know, the radius is based on an average of the existing points.

The problem we have in this example is that the great number of point in X1 crushes the values in X2. Therefore the radius never reaches values above the threshold.

Here is a function giving you the radius of the new cluster by integrating n points from X2, given that a subcluster based on X1 already exists. You see that the critical number of samples required is 209 (and you see that if you put 209 point in X2, you end up with 2 clusters):

```python
In [1]: def dist(n):
   ...:     return ((7200 + 0.02 * n) / (10000 + n)) - ((2 * (6000 + n * 0.1) ** 2)/ (10000 + n) ** 2)
   ...: 

In [2]: dist(100)
Out[2]: 0.004901480247034518

In [3]: dist(200)
Out[3]: 0.009611687812379865

In [4]: dist(209)
Out[4]: 0.010026512112923935
```

So why does it work if you swap X2 and X1? Because then the first subcluster created is X2, and X1 grows enough to reach the critical mass and make the radius big enough for the subcluster to split.

Solving this problem therefore requires to change the absorption formula. 

I found one [paper that suggests more stable ways to compute features called BETULA features](https://arxiv.org/pdf/2006.12881.pdf) but I implemented them and they are iso to the original birch formula.

We can imagine other absorption criteria for the algorithm, but, as mentioned in the original paper, BIRCH stays dependent on the order of arrival of the data. The only solution I can see without touching the current formula would be to find a way to precondition the data, to shuffle them to prevent the formation of big subclusters.

Does anybody have an opinion on this?

Labels: Bug, module:cluster