docs: def __cost(self, params, phase, X):
          """Computes activation, cost function, and derivative."""
          params = self.__roll(params)
          a = np.concatenate((np.ones((X.shape[0], 1)), X), axis=1) # This is a1
          calculated_a = [a] # a1 is at index 0, a_n is at index n-1
          calculated_z = [0] # There is no z1, z_n is at index n-1
          for i, theta in enumerate(params): # calculated_a now contains a1, a2, a3 if there was only one hidden layer (two theta matrices)
              z = calculated_a[-1] * theta.transpose() # z_n = a_n-1 * Theta_n-1'
              calculated_z.append(z) # Save the new z_n
              a = self.sigmoid(z) # a_n = sigmoid(z_n)
              if i != len(params) - 1: # Don't append extra ones for the output layer
                  a = np.concatenate((np.ones((a.shape[0], 1)), a), axis=1) # Append the extra column of ones for all other layers
              calculated_a.append(a) # Save the new a
  
          if phase == 0:
              if self.__num_labels > 1:
                  return np.argmax(calculated_a[-1], axis=1)
              return np.round(calculated_a[-1])
  
          J = np.sum(-np.multiply(self.__y, np.log(calculated_a[-1]))-np.multiply(1-self.__y, np.log(1-calculated_a[-1])))/self.__m; # Calculate cost
          if self.__lambda != 0: # If we're using regularization...
              J += np.sum([np.sum(np.power(theta[:,1:], 2)) for theta in params])*self.__lambda/(2.0*self.__m) # ...add it from all theta matrices
  
          if phase == 1:
              return J
  
          reversed_d = []
          reversed_theta_grad = []
          for i in range(len(params)): # For once per theta matrix...
              if i == 0: # ...if it's the first one...
                  d = calculated_a[-1] - self.__y # ...initialize the error...
              else: # ...otherwise d_n-1 = d_n * Theta_n-1[missing ones] .* sigmoid(z_n-1)
                  d = np.multiply(reversed_d[-1]*params[-i][:,1:], self.sigmoid_grad(calculated_z[-1-i])) # With i=1/1 hidden layer we're getting Theta2 at index -1, and z2 at index -2
              reversed_d.append(d)
              theta_grad = reversed_d[-1].transpose() * calculated_a[-i-2] / self.__m
              if self.__lambda != 0:
                  theta_grad += np.concatenate((np.zeros((params[-1-i].shape[0], 1)), params[-1-i][:,1:]), axis=1) * self.__lambda / self.__m# regularization
              reversed_theta_grad.append(theta_grad)
          theta_grad = self.__unroll(reversed(reversed_theta_grad))
          return theta_grad