Possible bug in the implementation of histogram in histogram based gradient boosting.

### Describe the bug

I was looking at the code in the file `sklearn/ensemble/_hist_gradient_boosting/histogram.pyx` and I believe there might be a bug. The note at the top of the file describes the following:

```
# Notes:
# - IN views are read-only, OUT views are write-only
# - In a lot of functions here, we pass feature_idx and the whole 2d
#   histograms arrays instead of just histograms[feature_idx]. This is because
#   Cython generated C code will have strange Python interactions (likely
#   related to the GIL release and the custom histogram dtype) when using 1d
#   histogram arrays that come from 2d arrays.
# - The for loops are un-wrapped, for example:
#
#   for i in range(n):
#       array[i] = i
#
#   will become
#
#   for i in range(n // 4):
#       array[i] = i
#       array[i + 1] = i + 1
#       array[i + 2] = i + 2
#       array[i + 3] = i + 3
#
#   This is to hint gcc that it can auto-vectorize these 4 operations and
#   perform them all at once.
```

This is used in several places in the file for example here:

```

cpdef void _build_histogram(
        const int feature_idx,
        const unsigned int [::1] sample_indices,  # IN
        const X_BINNED_DTYPE_C [::1] binned_feature,  # IN
        const G_H_DTYPE_C [::1] ordered_gradients,  # IN
        const G_H_DTYPE_C [::1] ordered_hessians,  # IN
        hist_struct [:, ::1] out) nogil:  # OUT
    """Return histogram for a given feature."""
    cdef:
        unsigned int i = 0
        unsigned int n_node_samples = sample_indices.shape[0]
        unsigned int unrolled_upper = (n_node_samples // 4) * 4

        unsigned int bin_0
        unsigned int bin_1
        unsigned int bin_2
        unsigned int bin_3
        unsigned int bin_idx

    for i in range(0, unrolled_upper, 4):
        bin_0 = binned_feature[sample_indices[i]]
        bin_1 = binned_feature[sample_indices[i + 1]]
        bin_2 = binned_feature[sample_indices[i + 2]]
        bin_3 = binned_feature[sample_indices[i + 3]]

        out[feature_idx, bin_0].sum_gradients += ordered_gradients[i]
        out[feature_idx, bin_1].sum_gradients += ordered_gradients[i + 1]
        out[feature_idx, bin_2].sum_gradients += ordered_gradients[i + 2]
        out[feature_idx, bin_3].sum_gradients += ordered_gradients[i + 3]

        out[feature_idx, bin_0].sum_hessians += ordered_hessians[i]
        out[feature_idx, bin_1].sum_hessians += ordered_hessians[i + 1]
        out[feature_idx, bin_2].sum_hessians += ordered_hessians[i + 2]
        out[feature_idx, bin_3].sum_hessians += ordered_hessians[i + 3]

        out[feature_idx, bin_0].count += 1
        out[feature_idx, bin_1].count += 1
        out[feature_idx, bin_2].count += 1
        out[feature_idx, bin_3].count += 1

    for i in range(unrolled_upper, n_node_samples):
        bin_idx = binned_feature[sample_indices[i]]
        out[feature_idx, bin_idx].sum_gradients += ordered_gradients[i]
        out[feature_idx, bin_idx].sum_hessians += ordered_hessians[i]
        out[feature_idx, bin_idx].count += 1
```

However, there is an issue. Here we are indexing over bins for example `bin_0`, `bin_1`, and so on. Nevertheless, it is possible that the bins are the same (this is in fact not too unlikely because there are only 256 possible bins). The issue is that if this happens, and it is the case that the compiler executes the function as the note describes, then we will get race conditions because it is possible that we modify the same memory in `out` simultaneously because of the incorrect parallel execution.  I am not sure if this is a bug or not because it depends on how this is compiled -- but at the very least the code and the comments are incorrect. 

### Steps/Code to Reproduce

There is no code to reproduce the problem. 

### Expected Results

The expected result is that the computation should have the same output as done without vectorization and roop unrolling but I am not sure this is the case. 

### Actual Results

It is unclear what the actual results are but there is an issue at the very least with what is written in the code. 

### Versions

It is present in branch 1.2.x of the repo (and the main branch as of Dec 22 2022)

COMMENT:
Histograms functions are hard to vectorize, see
- https://stackoverflow.com/questions/12985949/methods-to-vectorise-histogram-in-simd
- Shahbahrami, Asadollah et al. "SIMD Vectorization of Histogram Functions." 2007
  IEEE International Conf. on Application-specific Systems, Architectures and
  Processors (ASAP) (2007): 174-179.

Labels: Bug, Needs Triage