def METHOD_NAME(X_input: np.ndarray, F: Callable, D: float, dim: int, N: int, lamada_: float = 1) -> float:
    """Calculate action by path integral by Wang's method.
    Quantifying the Waddington landscape and biological paths for development and differentiation. Jin Wang, Kun Zhang,
    Li Xu, and Erkang Wang, PNAS, 2011
    Args:
        X_input: The initial guess of the least action path. Default is a straight line connecting the starting and end path.
        F: The reconstructed vector field function. This is assumed to be time-independent.
        D: The diffusion constant. Note that this can be a space-dependent matrix.
        dim: The feature numbers of the input data.
        N: Number of waypoints along the least action path.
        lamada_: Regularization parameter
    Returns:
        The action function calculated by the Hamilton-Jacobian method.
    """
    X_input = X_input.reshape((int(dim), -1)) if len(X_input.shape) == 1 else X_input
    delta, delta_l = delta_delta_l(X_input)
    V_m = np.zeros((N, 1))
    F_l = np.zeros((N, 1))
    E_eff = np.zeros((N, 1))
    for i in range(N - 1):
        F_m = F(X_input[:, i]).reshape((1, -1))
        V_m[i] = V(F, D, X_input[:, i])
        E_eff[i] = np.sum(F(X_input[:, i]) ** 2) / (4 * D) - V_m[i]
        F_l[i] = F_m.dot(delta[:, i]) / delta_l[i]
    P = np.sum((delta_l - np.linalg.norm(X_input[:, N - 1] - X_input[:, 0]) / N) ** 2)
    S_HJ = np.sum((np.sqrt((E_eff + V_m[:N]) / D) - 1 / (2 * D) * F_l) * delta_l) + lamada_ * P
    print(S_HJ)
    return S_HJ