Ускоряем и оптимизируем numpy, pandas, scipy и sklearn

• Генерация координат: Программа создает два вектора случайных координат: rx и ry, которые находятся в диапазоне от 0 до 1. Эти координаты представляют собой точки на плоскости.

• Проверка попадания: Точка считается находящейся внутри круга, если расстояние dist от начала координат (0, 0) меньше радиуса, равного 1. Это условие можно проверить с помощью следующей формулы: rx² + ry² < 1.

• Вычисление: Итоговая оценка числа π (pi_est) вычисляется как отношение количества точек, попавших внутрь круга, к общему количеству сгенерированных точек.

Что такое полная бета-функция?

Чтобы понять неполную версию, нужно вспомнить полную. Полная бета-функция $\(B(a, b)\)$ — это определенный интеграл от нуля до единицы, который зависит от двух параметров $\(a\)$ и $\(b\)$:

$$\(B(a, b) = \int_{0}^{1} t^{a-1} (1-t)^{b-1} \, dt\)$$

Что такое неполная бета-функция?

В неполной бета-функции верхний предел интеграла заменяется на переменную $\(x\)$ (где $\(0 \le x \le 1\))$. Это значит, что мы интегрируем функцию не до конца, а только на отрезке от $\(0\)$ до $\(x\)$.

Обозначается она как $\(B_x(a, b)\)$ и определяется следующим образом:

$$\(B_x(a, b) = \int_{0}^{x} t^{a-1} (1-t)^{b-1} \, dt\)$$

from tabulate import tabulateimport numpy as npdef generate_data(rank, num_points, noise_level):    np.random.seed(42)    x = np.linspace(-10, 10, num_points)    y = np.linspace(-10, 10, num_points)    if rank == 1:        z = 3 * x + 5 + noise_level * np.random.randn(num_points)        data = np.column_stack((x, y, z))    elif rank == 2:        z = 2 * x + 3 * y + 5 + noise_level * np.random.randn(num_points)        data = np.column_stack((x, y, z))    elif rank == 3:        z = 2 * x**2 + 3 * y**2 + 5 + noise_level * np.random.randn(num_points)        data = np.column_stack((x, y, z))    return datadef GMT_trend2d(data, rank):    import numpy as np    from sklearn.linear_model import LinearRegression    # scale factor for normally distributed data is 1.4826    # https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.median_abs_deviation.html    MAD_NORMALIZE = 1.4826    # significance value    sig_threshold = 0.51    if rank not in [1,2,3]:        raise Exception('Number of model parameters "rank" should be 1, 2, or 3')    #see gmt_stat.c    def gmtstat_f_q (chisq1, nu1, chisq2, nu2):        import scipy.special as sc        if chisq1 == 0.0:            return 1        if chisq2 == 0.0:            return 0        return sc.betainc(0.5*nu2, 0.5*nu1, chisq2/(chisq2+chisq1))    if rank in [2,3]:        x = data[:,0]        x = np.interp(x, (x.min(), x.max()), (-1, +1))    if rank == 3:        y = data[:,1]        y = np.interp(y, (y.min(), y.max()), (-1, +1))    z = data[:,2]    w = np.ones(z.shape)    if rank == 1:        xy = np.expand_dims(np.zeros(z.shape),1)    elif rank == 2:        xy = np.expand_dims(x,1)    elif rank == 3:        xy = np.stack([x,y]).transpose()    # create linear regression object    mlr = LinearRegression()    chisqs = []    coeffs = []    while True:        # fit linear regression        mlr.fit(xy, z, sample_weight=w)        r = np.abs(z - mlr.predict(xy))        chisq = np.sum((r**2*w))/(z.size-3)        chisqs.append(chisq)        k = 1.5 * MAD_NORMALIZE * np.median(r)        w = np.where(r <= k, 1, (2*k/r) - (k * k/(r**2)))        sig = 1 if len(chisqs)==1 else gmtstat_f_q(chisqs[-1], z.size-3, chisqs[-2], z.size-3)        # Go back to previous model only if previous chisq < current chisq        if len(chisqs)==1 or chisqs[-2] > chisqs[-1]:            coeffs = [mlr.intercept_, *mlr.coef_]        #print ('chisq', chisq, 'significant', sig)        if sig < sig_threshold:            break    # get the slope and intercept of the line best fit    return (coeffs[:rank])def calculate_mse(data, coeffs, rank):    z_actual = data[:, 2]    if rank == 1:        z_predicted = coeffs[0]    elif rank == 2:        # Interpolate x the same way as in GMT_trend2d        x = data[:, 0]        x_interp = np.interp(x, (x.min(), x.max()), (-1, +1))        z_predicted = coeffs[0] + coeffs[1] * x_interp    elif rank == 3:        # Interpolate x and y the same way as in GMT_trend2d        x = data[:, 0]        x_interp = np.interp(x, (x.min(), x.max()), (-1, +1))        y = data[:, 1]        y_interp = np.interp(y, (y.min(), y.max()), (-1, +1))        z_predicted = coeffs[0] + coeffs[1] * x_interp + coeffs[2] * y_interp    mse = np.mean((z_actual - z_predicted) ** 2)    return msedef test_mse(num_points = 100*1000, ranks = [1, 2, 3], noise_levels = [0, 1, 10, 50]):    import warnings    results = []    # Suppress the specific warning    with warnings.catch_warnings():        warnings.simplefilter("ignore", category=RuntimeWarning)        for rank in ranks:            for noise_level in noise_levels:                data = generate_data(rank, num_points, noise_level)                # round the output                coeffs_gmt = [v.round(8) for v in GMT_trend2d(data, rank)]                mse_gmt = np.round(calculate_mse(data, coeffs_gmt, rank), 0)                results.append([rank, noise_level, mse_gmt])    headers = ["Rank", "Noise Level", "GMT_trend2d, MSE"]    print(tabulate(results, headers=headers))test_mse()

using namespace np;using namespace scipy;using namespace sklearn;auto generate_data(auto rank, auto num_points, auto noise_level) {    random::seed(42);    auto x = linspace(-10.0, 10.0, num_points);    auto y = linspace(-10.0, 10.0, num_points);    if (rank == 1) {        auto z = 3 * x + 5 + noise_level * random::randn(num_points);        return column_stack(x, y, z);    }    if (rank == 2) {        auto z = 2 * x + 3 * y + 5 + noise_level * random::randn(num_points);        return column_stack(x, y, z);    }    auto z = 2 * x * x + 3 * y * y + 5 + noise_level * random::randn(num_points);    return column_stack(x, y, z);}auto GMT_trend2d(const Array<float_> &data, int rank) {    float_ MAD_NORMALIZE = 1.4826;    float_ sig_threshold = 0.51;    if (rank != 1 && rank != 2 && rank != 3) {        throw sklearn::RuntimeError("Number of model parameters \"rank\" should be 1, 2, or 3");    }    auto gmtstat_f_q = [](float_ chisq1, float_ nu1, float_ chisq2, float_ nu2) {        if (chisq1 == 0.0) return 1.0;        if (chisq2 == 0.0) return 0.0;        return scipy::special::betainc(0.5 * nu2, 0.5 * nu1, chisq2 / (chisq2 + chisq1));    };    Array<float_> x;    if (rank == 2 || rank == 3) {        auto x_ = data[":,0"];        x = interp(x_, Array<float_>{x_.min(), x_.max()}, Array<float_>{-1, +1});    }    Array<float_> y;    if (rank == 3) {        auto y_ = data[":,1"];        y = interp(y_, Array<float_>{y_.min(), y_.max()}, Array<float_>{-1, +1});    }    auto z = data[":, 2"].copy();    Array<float_> w = ones(z.shape()).copy();    Array<float_> xy;    if (rank == 1) {        xy = expand_dims(zeros(z.shape()), 1);    } else if (rank == 2) {        xy = expand_dims(x, 1);    } else if (rank == 3) {        xy = stack(x, y).transpose();    }    auto mlr = linear_model::LinearRegression{};    std::vector<float_> chisqs;    Array<float_> coeffs;    while (true) {        mlr.fit(xy, z, w);        auto r = abs_sub(z, mlr.predict(xy));        auto chisq = sum_sq_weighted(r, w) / static_cast<float_>(z.size() - 3);        chisqs.push_back(chisq);        auto k = 1.5 * MAD_NORMALIZE * median(r);        w = where_tukey(r, k);        auto sig = (chisqs.size() == 1 ? 1 : gmtstat_f_q(chisqs[chisqs.size() - 1], static_cast<float_>(z.size() - 3), chisqs[chisqs.size() - 2], static_cast<float_>(z.size() - 3)));        if (chisqs.size() == 1 or chisqs[chisqs.size() - 2] > chisqs[chisqs.size() - 1]) {            coeffs = mlr.coeffs_();        }        if (sig < sig_threshold) {            break;        }    }    auto result = Array<float_>{};    for (int i = 0; i < rank; ++i) {        result = append(result, Array<float_>{coeffs.get(i)});    }    return result;auto calculate_mse(const Array<float_> &data, const Array<float_> &coeffs, int rank) {    auto z_actual = data[":,2"].copy();    Array<float_> z_predicted;    if (rank == 1) {        z_predicted = coeffs.get(0) * ones(z_actual.shape());    } else if (rank == 2) {        // Interpolate x the same way as in GMT_trend2d        auto x_ = data[":,0"];        auto x = interp(x_, Array<float_>{x_.min(), x_.max()}, Array<float_>{-1, +1});        z_predicted = coeffs.get(0) + coeffs.get(1) * x;    } else if (rank == 3) {        // Interpolate x and y the same way as in GMT_trend2d        auto x_ = data[":,0"];        auto x = interp(x_, Array<float_>{x_.min(), x_.max()}, Array<float_>{-1, +1});        auto y_ = data[":,1"];        auto y = interp(y_, Array<float_>{y_.min(), y_.max()}, Array<float_>{-1, +1});        z_predicted = coeffs.get(0) + coeffs.get(1) * x + coeffs.get(2) * y;    }    using sklearn::metrics::mean_squared_error;    using sklearn::metrics::MeanSquaredErrorParameters;    MeanSquaredErrorParameters<Array<float_>> params{.y_true = z_actual, .y_pred = z_predicted};    return mean_squared_error(params);}void test_mse(int num_points = 100 * 1000, const std::vector<int> &ranks = {1, 2, 3}, const std::vector<int> &noise_levels = {0, 1, 10, 50}) {    std::vector<std::tuple<int, int, float_>> results;    for (auto rank: ranks) {        for (auto noise_level: noise_levels) {            auto data = generate_data(rank, num_points, noise_level);            auto coeffs_gmt = GMT_trend2d(data, rank);            // round coefficients to 8 decimal places            auto coeffs_rounded = coeffs_gmt.copy();            for (std::size_t i = 0; i < coeffs_rounded.size(); ++i) {                coeffs_rounded.set(i, std::round(coeffs_rounded.get(i) * 1e8) / 1e8);            }            auto mse_gmt = calculate_mse(data, coeffs_rounded, rank);            // round MSE to zero decimal places            mse_gmt = std::round(mse_gmt);            results.emplace_back(rank, noise_level, mse_gmt);        }    }    // print table    std::cout << "Rank\tNoise Level\tGMT_trend2d, MSE\n";    for (const auto &[rank, noise_level, mse]: results) {        std::cout << rank << "\t" << noise_level << "\t" << mse << "\n";    }}