function met = hydroMetrics(Q_sim, Q_obs)
% HYDROMETRICS 
    mask = isfinite(Q_sim) & isfinite(Q_obs);
    sim = Q_sim(mask);
    obs = Q_obs(mask);

    met = struct('NSE', NaN, 'PBIAS', NaN, 'RMSE', NaN, 'KGE', NaN, 'n', numel(sim));

    if isempty(sim)
        return;
    end

    den = sum((obs - mean(obs)).^2);
    if den > eps
        met.NSE = 1 - sum((sim - obs).^2) / den;
    end

    denBias = sum(obs);
    if abs(denBias) > eps
        met.PBIAS = 100 * sum(sim - obs) / denBias;
    end

    met.RMSE = sqrt(mean((sim - obs).^2));

    if numel(sim) >= 2
        ms = mean(sim);
        mo = mean(obs);
        so = std(sim);
        oo = std(obs);
        if oo > eps && abs(mo) > eps
            if so > eps
                r = sum((sim - ms) .* (obs - mo)) / ((numel(sim) - 1) * so * oo);
            else
                r = 0;
            end
            alpha = so / oo;
            beta = ms / mo;
            met.KGE = 1 - sqrt((r - 1)^2 + (alpha - 1)^2 + (beta - 1)^2);
        end
    end
end
