code: 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