Bug - wrong activations and delta length in BaseMultilayerPerceptron ?

In multilayer_perceptron.py at line 623 we have

activations = [X]
activations.extend(np.empty((batch_size, n_fan_out)) for n_fan_out in layer_units[1:])
deltas = [np.empty_like(a_layer) for a_layer in activations]
coef_grads = [ np.empty((n_fan_in_, n_fan_out_)) 
                        for n_fan_in_, n_fan_out_ in zip(layer_units[:-1], layer_units[1:]) ]

and at line 268 we have defined

self.n_layers_ = len(layer_units) 

as such activations and deltas should has the same the layer_units.

However, in the doc strings from other methods of this class, it is claimed at both activations and deltas has the length of self.n_layers - 1, and it seems to me the logic within those methods also follows this assumption. For exampe at line 236, we have

last = self.n_layers_ - 2
deltas[last] = activations[-1] - y

This does not look right to me. If the network has only three layer (input, hidden, and output) the last would be refer to the the hidden layer. But if my understanding is right the deltas formula is meant for the output layer.

Start in sklearn/neural_network/multilayer_perceptron.py around the cited activation, delta, and n_layers_ definitions. Trace the training methods and their docstrings to determine whether the output-layer delta indexing and the lengths of activations and deltas are consistent. Done means the implementation and documentation agree and the output-layer calculation is correctly addressed.