input: Write a python function to implement feedforward and backpropagation for a simple neural network and update the weights and biases. Use the provided dataset and initial weights. 
  assert np.allclose(output, np.array([[0.06108959]]), atol=1e-2)
  assert np.allclose(d_output.shape, (4, 1))
  assert np.allclose(W1.shape, (3, 4))
output: def simple_neural_network(X, y, W1, b1, W2, b2, learning_rate):
      def sigmoid(x):
          return 1 / (1 + np.exp(-x))
      
      def sigmoid_derivative(x):
          return x * (1 - x)
      
      # Feedforward
      Z1 = np.dot(X, W1) + b1
      A1 = sigmoid(Z1)
      Z2 = np.dot(A1, W2) + b2
      A2 = sigmoid(Z2)
      
      output = A2
      
      # Calculate error
      error = y - output
      d_output = error * sigmoid_derivative(output)
      
      # Backpropagation
      error_hidden_layer = d_output.dot(W2.T)
      d_hidden_layer = error_hidden_layer * sigmoid_derivative(A1)
      
      # Update weights and biases
      W2 += A1.T.dot(d_output) * learning_rate
      b2 += np.sum(d_output, axis=0, keepdims=True) * learning_rate
      W1 += X.T.dot(d_hidden_layer) * learning_rate
      b1 += np.sum(d_hidden_layer, axis=0, keepdims=True) * learning_rate
      
      return output, d_output, W1
  
  import numpy as np
  # Provided dataset and initial weights
  X = np.array([[0, 0], [0, 1], [1, 0], [1, 1]])
  y = np.array([[0], [1], [1], [0]])
  W1 = np.random.uniform(size=(2, 4))
  b1 = np.random.uniform(size=(1, 4))
  W2 = np.random.uniform(size=(4, 1))
  b2 = np.random.uniform(size=(1, 1))
  learning_rate = 0.1
  
  output, d_output, W1 = simple_neural_network(X, y, W1, b1, W2, b2, learning_rate)