def kge_2009(simulated_array, observed_array, s=(1, 1, 1), replace_nan=None,
             replace_inf=None, remove_neg=False, remove_zero=False, return_all=False):
    """Compute the Kling-Gupta efficiency (2009).

    .. image:: /pictures/KGE_2009.png

    **Range:** -inf < KGE (2009) < 1, larger is better.

    **Notes:** Gupta et al. (2009) created this metric to demonstrate the relative importance of
    the three components of the NSE, which are correlation, bias and variability. This was done
    with hydrologic modeling as the context. This metric is meant to address issues with the NSE.

    Parameters
    ----------
    simulated_array: one dimensional ndarray
        An array of simulated data from the time series.

    observed_array: one dimensional ndarray
        An array of observed data from the time series.

    s: tuple of length three
        Represents the scaling factors to be used for re-scaling the Pearson product-moment
        correlation coefficient (r), Alpha, and Beta, respectively.

    replace_nan: float, optional
        If given, indicates which value to replace NaN values with in the two arrays. If None, when
        a NaN value is found at the i-th position in the observed OR simulated array, the i-th value
        of the observed and simulated array are removed before the computation.

    replace_inf: float, optional
        If given, indicates which value to replace Inf values with in the two arrays. If None, when
        an inf value is found at the i-th position in the observed OR simulated array, the i-th
        value of the observed and simulated array are removed before the computation.

    remove_neg: boolean, optional
        If True, when a negative value is found at the i-th position in the observed OR simulated
        array, the i-th value of the observed AND simulated array are removed before the
        computation.

    remove_zero: boolean, optional
        If true, when a zero value is found at the i-th position in the observed OR simulated
        array, the i-th value of the observed AND simulated array are removed before the
        computation.

    return_all: bool
        If True, returns all of the components of the KGE metric, which are r, alpha, and beta, respectively.

    Returns
    -------
    float (tuple of float)
        The Kling-Gupta (2009) efficiency value, unless the return_all parameter is True.

    Examples
    --------

    >>> import HydroErr as he
    >>> import numpy as np

    >>> sim = np.array([5, 7, 9, 2, 4.5, 6.7])
    >>> obs = np.array([4.7, 6, 10, 2.5, 4, 6.8])
    >>> he.kge_2009(sim, obs)
    0.912223072345668

    >>> he.kge_2009(sim, obs, return_all=True)  # Returns (r, alpha, beta, kge)
    (0.9615951377405804, 0.927910707932087, 1.0058823529411764, 0.9181073779138655)

    References
    ----------
    - Gupta, H. V., Kling, H., Yilmaz, K. K., & Martinez, G. F. (2009). Decomposition of the mean
      squared error and NSE performance criteria: Implications for improving hydrological modelling.
      Journal of Hydrology, 377(1-2), 80-91.
    """

    # Treats data
    simulated_array, observed_array = treat_values(
        simulated_array,
        observed_array,
        replace_nan=replace_nan,
        replace_inf=replace_inf,
        remove_neg=remove_neg,
        remove_zero=remove_zero
    )

    # Means
    sim_mean = np.mean(simulated_array)
    obs_mean = np.mean(observed_array)

    # Standard Deviations
    sim_sigma = np.std(simulated_array, ddof=1)
    obs_sigma = np.std(observed_array, ddof=1)

    # Pearson R
    top_pr = np.sum((observed_array - obs_mean) * (simulated_array - sim_mean))
    bot1_pr = np.sqrt(np.sum((observed_array - obs_mean) ** 2))
    bot2_pr = np.sqrt(np.sum((simulated_array - sim_mean) ** 2))
    pr = top_pr / (bot1_pr * bot2_pr)

    # Ratio between mean of simulated and observed data
    if obs_mean != 0:
        beta = sim_mean / obs_mean
    else:
        beta = np.nan

    # Relative variability between simulated and observed values
    if obs_sigma != 0:
        alpha = sim_sigma / obs_sigma
    else:
        alpha = np.nan

    if not np.isnan(beta) and not np.isnan(alpha):
        kge = 1 - np.sqrt(
            (s[0] * (pr - 1)) ** 2 + (s[1] * (alpha - 1)) ** 2 + (s[2] * (beta - 1)) ** 2)
    else:
        if obs_mean == 0:
            warnings.warn(
                'Warning: The observed data mean is 0. Therefore, Beta is infinite and the KGE '
                'value cannot be computed.')
        if obs_sigma == 0:
            warnings.warn(
                'Warning: The observed data standard deviation is 0. Therefore, Alpha is infinite '
                'and the KGE value cannot be computed.')
        kge = np.nan

    assert type(return_all) == bool, "expected <type 'bool'> for parameter return_all, got {}".format(type(return_all))

    if return_all:
        return pr, alpha, beta, kge
    else:
        return kge
Compute the Kling-Gupta efficiency (2009).

    .. image:: /pictures/KGE_2009.png

    **Range:** -inf < KGE (2009) < 1, larger is better.

    **Notes:** Gupta et al. (2009) created this metric to demonstrate the relative importance of
    the three components of the NSE, which are correlation, bias and variability. This was done
    with hydrologic modeling as the context. This metric is meant to address issues with the NSE.

    Parameters
    ----------
    simulated_array: one dimensional ndarray
        An array of simulated data from the time series.

    observed_array: one dimensional ndarray
        An array of observed data from the time series.

    s: tuple of length three
        Represents the scaling factors to be used for re-scaling the Pearson product-moment
        correlation coefficient (r), Alpha, and Beta, respectively.

    replace_nan: float, optional
        If given, indicates which value to replace NaN values with in the two arrays. If None, when
        a NaN value is found at the i-th position in the observed OR simulated array, the i-th value
        of the observed and simulated array are removed before the computation.

    replace_inf: float, optional
        If given, indicates which value to replace Inf values with in the two arrays. If None, when
        an inf value is found at the i-th position in the observed OR simulated array, the i-th
        value of the observed and simulated array are removed before the computation.

    remove_neg: boolean, optional
        If True, when a negative value is found at the i-th position in the observed OR simulated
        array, the i-th value of the observed AND simulated array are removed before the
        computation.

    remove_zero: boolean, optional
        If true, when a zero value is found at the i-th position in the observed OR simulated
        array, the i-th value of the observed AND simulated array are removed before the
        computation.

    return_all: bool
        If True, returns all of the components of the KGE metric, which are r, alpha, and beta, respectively.

    Returns
    -------
    float (tuple of float)
        The Kling-Gupta (2009) efficiency value, unless the return_all parameter is True.

    Examples
    --------

    >>> import HydroErr as he
    >>> import numpy as np

    >>> sim = np.array([5, 7, 9, 2, 4.5, 6.7])
    >>> obs = np.array([4.7, 6, 10, 2.5, 4, 6.8])
    >>> he.kge_2009(sim, obs)
    0.912223072345668

    >>> he.kge_2009(sim, obs, return_all=True)  # Returns (r, alpha, beta, kge)
    (0.9615951377405804, 0.927910707932087, 1.0058823529411764, 0.9181073779138655)

    References
    ----------
    - Gupta, H. V., Kling, H., Yilmaz, K. K., & Martinez, G. F. (2009). Decomposition of the mean
      squared error and NSE performance criteria: Implications for improving hydrological modelling.
      Journal of Hydrology, 377(1-2), 80-91.