docs: def func(self, w, *args):
          """Return the costs of the neural network for predictions.
  
          Args:
              w (array of float): weight vectors such that:
                  w[:-h1] -- weights between the input and h layers
                  w[-h1:] -- weights between the h and output layers
              args: features (args[0]) and target (args[1])
  
          Returns:
              combined cost of RMSE, L1, and L2 regularization
          """
          x0 = args[0]
          x1 = args[1]
  
          n0 = x0.shape[0]
          n1 = x1.shape[0]
  
          # n -- number of pairs to evaluate
          n = max(n0, n1) * 10
          idx0 = np.random.choice(range(n0), size=n)
          idx1 = np.random.choice(range(n1), size=n)
  
          # b -- bias for the input and h layers
          b0 = np.ones((n0, 1))
          b1 = np.ones((n1, 1))
          i1 = self.i + 1
          h = self.h
          h1 = h + 1
  
          # Predict for features -- cannot use predict_raw() because here
          # different weights can be used.
          if sparse.issparse(x0):
              p0 = np.hstack((sigm(sparse.hstack((x0, b0)).dot(w[:-h1].reshape(
                                 i1, h))), b0)).dot(w[-h1:].reshape(h1, 1))
              p1 = np.hstack((sigm(sparse.hstack((x1, b1)).dot(w[:-h1].reshape(
                                 i1, h))), b1)).dot(w[-h1:].reshape(h1, 1))
          else:
              p0 = np.hstack((sigm(np.hstack((x0, b0)).dot(w[:-h1].reshape(
                                 i1, h))), b0)).dot(w[-h1:].reshape(h1, 1))
              p1 = np.hstack((sigm(np.hstack((x1, b1)).dot(w[:-h1].reshape(
                                 i1, h))), b1)).dot(w[-h1:].reshape(h1, 1))
  
          p0 = p0[idx0]
          p1 = p1[idx1]
  
          # Return the cost that consists of the sum of squared error +
          # L2-regularization for weights between the input and h layers +
          # L2-regularization for weights between the h and output layers.
          #return .5 * (sum((1 - sigm(p1 - p0)) ** 2) + self.l1 * sum(w[:-h1] ** 2) +
          return .5 * (sum((1 - p1 + p0) ** 2) / n +
                       self.l1 * sum(w[:-h1] ** 2) / (i1 * h) +
                       self.l2 * sum(w[-h1:] ** 2) / h1)