Simple NN in python not working, problem with backpropagation algorithm perhaps?

Question

I tried to code classic XOR problem using NN with 2 inputs, 2 hidden neurons, 1 output neuron (using stochastic gradient descent). But whatever i do my NN doesn't work correctly, still getting wrong output and i really don't know where is the problem, i think in backpropagation, maybe i multiplied it in wrong order, i really don't know. The formula i used for updating the outputs weights is LEARING_RATE * ERROR * (OUTPUT * (1 - OUTPUT) * HIDDEN_OUTPUT.T and the formula for hidden weights is LEARING_RATE * HIDDEN_ERROR * (HIDDEN_OUTPUT * (1 - HIDDEN_OUTPUT) * INPUT. I would like to know where is the problem in my code, I will be grateful for any help you can provide :)

def sigmoid(x):
    return 1 / (1 + np.exp(-x))


def sigmoid_derivative(y):
    return np.exp(y) * (1 - np.exp(y))


inp = np.array([ [[0, 0]],
                 [[0, 1]],
                 [[1, 0]],
                 [[1, 1]] ])

w_ih = np.random.uniform(size=(2, len(inp[0][0])))  # weights from input to hidden

bias_hidden = np.random.uniform(size=(1, len(w_ih)))

w_ho = np.random.uniform(size=(1, len(w_ih)))       # weights from hidden to output

bias_output = np.random.uniform(size=(1, len(w_ho)))

des_output = [[0, 1, 1, 0]]

r = 0.1 # learning rate

for i in range(10000):
    index_sample = random.randint(0, len(inp) - 1)

    # Forward propagation
    hidden = sigmoid(np.dot(w_ih, inp[index_sample].T) + bias_hidden.T)
    output = sigmoid(np.dot(w_ho, hidden) + bias_output)
    # Calculating error
    output_error = output - np.array(des_output).T[index_sample]
    hidden_error = np.dot(w_ho.T, output_error)

    # Back propagation
    output_gradient = output_error * sigmoid_derivative(output)
    w_ho += np.dot(output_gradient, hidden.T) * r

    bias_output += np.sum(output_gradient, axis=0, keepdims=True) * r


    hidden_gradient = hidden_error * sigmoid_derivative(hidden)
    w_ih += np.dot(hidden_gradient, inp[index_sample]) * r

    bias_hidden += np.sum(np.array(hidden_gradient), axis=0, keepdims=True) * r


test = np.array([ [[0, 0]],
                  [[0, 1]],
                  [[1, 0]],
                  [[1, 1]] ])
for i in test:
    h1_test = sigmoid(np.dot(w_ih, np.array(i).T) + bias_hidden.T)      # testing neural newtwork
    output_test = sigmoid(np.dot(w_ho, h1_test) + bias_output)

    print(output_test, i)