docs: def fprime(self, w, *args):
          """Return the derivatives of the cost function 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:
              gradients of the cost function for predictions
          """
  
          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
          b = np.ones((n, 1))
          i1 = self.i + 1
          h = self.h
          h1 = h + 1
  
          w2 = w[-h1:].reshape(h1, 1)
          w1 = w[:-h1].reshape(i1, h)
  
          if sparse.issparse(x0):
              x0 = x0.tocsr()[idx0]
              x1 = x1.tocsr()[idx1]
              xb0 = sparse.hstack((x0, b))
              xb1 = sparse.hstack((x1, b))
          else:
              x0 = x0[idx0]
              x1 = x1[idx1]
              xb0 = np.hstack((x0, b))
              xb1 = np.hstack((x1, b))
  
          z0 = np.hstack((sigm(xb0.dot(w1)), b))
          z1 = np.hstack((sigm(xb1.dot(w1)), b))
          y0 = z0.dot(w2)
          y1 = z1.dot(w2)
  
          #e = 1 - sigm(y1 - y0)
          #dy = e * dsigm(y1 - y0)
          e = 1 - (y1 - y0)
          dy = e / n
  
          # Calculate the derivative of the cost function w.r.t. F and w2 where:
          # F -- weights between the input and h layers
          # w2 -- weights between the h and output layers
          dw1 = -(xb1.T.dot(dy.dot(w2[:-1].reshape(1, h)) * dsigm(xb1.dot(w1))) -
                 xb0.T.dot(dy.dot(w2[:-1].reshape(1, h)) * dsigm(xb0.dot(w1)))
                         ).reshape(i1 * h) + self.l1 * w[:-h1] / (i1 * h)
          dw2 = -(z1 - z0).T.dot(dy).reshape(h1) + self.l2 * w[-h1:] / h1
  
          return np.append(dw1, dw2)