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
Computes activation, cost function, and derivative.