Pass X through a soft-plus function, , in a numerically
        stable way (using the log-sum-exp trick).

        The softplus transformation is:

        .. math::
            \log(1 + \exp\{X\})

        Parameters
        ----------
            X: ndarray
                shape (N,) array or shape (N, D) array of data.

        Returns
        -------
            spX: ndarray
                array of same shape of X with the result of softmax(X).

def softplus(X):
    """ Pass X through a soft-plus function, , in a numerically
        stable way (using the log-sum-exp trick).

        The softplus transformation is:

        .. math::
            \log(1 + \exp\{X\})

        Parameters
        ----------
            X: ndarray
                shape (N,) array or shape (N, D) array of data.

        Returns
        -------
            spX: ndarray
                array of same shape of X with the result of softmax(X).
    """

    if np.isscalar(X):
        return logsumexp(np.vstack((np.zeros(1), [X])).T, axis=1)[0]

    N = X.shape[0]

    if X.ndim == 1:
        return logsumexp(np.vstack((np.zeros(N), X)).T, axis=1)
    elif X.ndim == 2:
        sftX = np.empty(X.shape, dtype=float)
        for d in range(X.shape[1]):
            sftX[:, d] = logsumexp(np.vstack((np.zeros(N), X[:, d])).T, axis=1)
        return sftX
    else:
        raise ValueError("This only works on up to 2D arrays.")