def __cost(self, params, phase, X):
        params = self.__roll(params)
        a = np.concatenate((np.ones((X.shape[0], 1)), X), axis=1)
        calculated_a = [a]
        calculated_z = [0]
        for i, theta in enumerate(params):
            z = calculated_a[-1] * theta.transpose()
            calculated_z.append(z)
            a = self.sigmoid(z)
            if i != len(params) - 1:
                a = np.concatenate((np.ones((a.shape[0], 1)), a), axis=1)
            calculated_a.append(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;
        if self.__lambda != 0:
            J += np.sum([np.sum(np.power(theta[:,1:], 2)) for theta in params])*self.__lambda/(2.0*self.__m)
        if phase == 1:
            return J
        reversed_d = []
        reversed_theta_grad = []
        for i in range(len(params)):
            if i == 0:
                d = calculated_a[-1] - self.__y
            else:
                d = np.multiply(reversed_d[-1]*params[-i][:,1:], self.sigmoid_grad(calculated_z[-1-i]))
            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
            reversed_theta_grad.append(theta_grad)
        theta_grad = self.__unroll(reversed(reversed_theta_grad))
        return theta_grad