<commit_before>from __future__ import print_function
import time
import random
import numpy as np
from selection.tests.instance import gaussian_instance
from selection.bayesian.initial_soln import selection, instance
from selection.randomized.api import randomization
from selection.bayesian.cisEQTLS.Simes_selection import simes_selection
from selection.bayesian.cisEQTLS.inference_per_gene import selection_probability_variants, \
    sel_prob_gradient_map_lasso, selective_inf_lasso
from scipy.stats import norm as normal
from selection.bayesian.cisEQTLS.Simes_selection import BH_q


def one_trial(n=350, p= 5000, s= 10, snr = 5., seed_n = 19, method="theoretical"):

    random.seed(seed_n)

    sample = instance(n=n, p=p, s=s, sigma=1, rho=0, snr=snr)

    X, y, true_beta, nonzero, noise_variance = sample.generate_response()

    n, p = X.shape

    random_Z = np.random.standard_normal(p)
    sel = selection(X, y, random_Z, method ="theoretical")
    lam, epsilon, active, betaE, cube, initial_soln = sel

    if sel is not None:

        lagrange = lam * np.ones(p)
        active_sign = np.sign(betaE)
        nactive = active.sum()

        feasible_point = np.fabs(betaE)

        noise_variance = 1.

        randomizer = randomization.isotropic_gaussian((p,), 1.)

        generative_X = X[:, active]
        prior_variance = 1000.

        grad_map = sel_prob_gradient_map_lasso(X,
                                               feasible_point,
                                               active,
                                               active_sign,
                                               lagrange,
                                               generative_X,
                                               noise_variance,
                                               randomizer,
                                               epsilon)

        inf = selective_inf_lasso(y, grad_map, prior_variance)

        samples = inf.posterior_samples()

        adjusted_intervals = np.vstack([np.percentile(samples, 5, axis=0), np.percentile(samples, 95, axis=0)])

        selective_mean = np.mean(samples, axis=0)

        Q = np.linalg.inv(prior_variance * (generative_X.dot(generative_X.T)) + noise_variance * np.identity(n))
        post_mean = prior_variance * ((generative_X.T.dot(Q)).dot(y))
        post_var = prior_variance * np.identity(nactive) - ((prior_variance ** 2) * (generative_X.T.dot(Q).dot(generative_X)))
        unadjusted_intervals = np.vstack([post_mean - 1.65 * (post_var.diagonal()), post_mean + 1.65 * (post_var.diagonal())])

        coverage_ad = np.zeros(p)
        coverage_unad = np.zeros(p)
        nerr = 0.
        true_val = true_beta[active]
        active_set = [i for i in range(p) if active[i]]
        active_ind = np.zeros(p)
        active_ind[active_set] = 1

        list_results = []
        list_results.append(true_beta)
        list_results.append(active_ind)

        ad_lower_credible = np.zeros(p)
        ad_upper_credible = np.zeros(p)

        unad_lower_credible = np.zeros(p)
        unad_upper_credible = np.zeros(p)

        ad_mean = np.zeros(p)
        unad_mean = np.zeros(p)

        if nactive > 1:
            try:
                for l in range(nactive):
                    if (adjusted_intervals[0, l] <= true_val[l]) and (true_val[l] <= adjusted_intervals[1, l]):
                        coverage_ad[active_set[l]] += 1
                    if (unadjusted_intervals[0, l] <= true_val[l]) and (true_val[l] <= unadjusted_intervals[1, l]):
                        coverage_unad[active_set[l]] += 1
                    ad_lower_credible[active_set[l]] = adjusted_intervals[0, l]
                    ad_upper_credible[active_set[l]] = adjusted_intervals[1, l]
                    unad_lower_credible[active_set[l]] = unadjusted_intervals[0, l]
                    unad_upper_credible[active_set[l]] = unadjusted_intervals[1, l]
                    ad_mean[active_set[l]] = selective_mean[l]
                    unad_mean[active_set[l]] = post_mean[l]

            except ValueError:
                nerr += 1
                print('ignore iteration raising ValueError')

            list_results.append(coverage_ad)
            list_results.append(coverage_unad)

            ngrid = 1000
            quantiles = np.zeros((ngrid, nactive))
            for i in range(ngrid):
                quantiles[i, :] = np.percentile(samples, (i * 100.) / ngrid, axis=0)

            index_grid = np.argmin(np.abs(quantiles - np.zeros((ngrid, nactive))), axis=0)
            p_value = 2 * np.minimum(np.true_divide(index_grid, ngrid), 1. - np.true_divide(index_grid, ngrid))
            p_BH = BH_q(p_value, 0.10)

            D_BH = np.zeros(p)

            if p_BH is not None:
                for indx in p_BH[1]:
                    D_BH[active_set[indx]] = 1

            list_results.append(D_BH)

            list_results.append(ad_lower_credible)
            list_results.append(ad_upper_credible)
            list_results.append(unad_lower_credible)
            list_results.append(unad_upper_credible)

            list_results.append(ad_mean)<commit_msg>added another argument in test for BH level<commit_after>from __future__ import print_function
import time
import random
import numpy as np
from selection.tests.instance import gaussian_instance
from selection.bayesian.initial_soln import selection, instance
from selection.randomized.api import randomization
from selection.bayesian.cisEQTLS.Simes_selection import simes_selection
from selection.bayesian.cisEQTLS.inference_per_gene import selection_probability_variants, \
    sel_prob_gradient_map_lasso, selective_inf_lasso
from scipy.stats import norm as normal
from selection.bayesian.cisEQTLS.Simes_selection import BH_q


def one_trial(n=350, p= 5000, s= 10, snr = 5., seed_n = 19, bh_level=0.1, method="theoretical"):

    random.seed(seed_n)

    sample = instance(n=n, p=p, s=s, sigma=1, rho=0, snr=snr)

    X, y, true_beta, nonzero, noise_variance = sample.generate_response()

    n, p = X.shape

    random_Z = np.random.standard_normal(p)
    sel = selection(X, y, random_Z, method ="theoretical")
    lam, epsilon, active, betaE, cube, initial_soln = sel

    if sel is not None:

        lagrange = lam * np.ones(p)
        active_sign = np.sign(betaE)
        nactive = active.sum()

        feasible_point = np.fabs(betaE)

        noise_variance = 1.

        randomizer = randomization.isotropic_gaussian((p,), 1.)

        generative_X = X[:, active]
        prior_variance = 1000.

        grad_map = sel_prob_gradient_map_lasso(X,
                                               feasible_point,
                                               active,
                                               active_sign,
                                               lagrange,
                                               generative_X,
                                               noise_variance,
                                               randomizer,
                                               epsilon)

        inf = selective_inf_lasso(y, grad_map, prior_variance)

        samples = inf.posterior_samples()

        adjusted_intervals = np.vstack([np.percentile(samples, 5, axis=0), np.percentile(samples, 95, axis=0)])

        selective_mean = np.mean(samples, axis=0)

        Q = np.linalg.inv(prior_variance * (generative_X.dot(generative_X.T)) + noise_variance * np.identity(n))
        post_mean = prior_variance * ((generative_X.T.dot(Q)).dot(y))
        post_var = prior_variance * np.identity(nactive) - ((prior_variance ** 2) * (generative_X.T.dot(Q).dot(generative_X)))
        unadjusted_intervals = np.vstack([post_mean - 1.65 * (post_var.diagonal()), post_mean + 1.65 * (post_var.diagonal())])

        coverage_ad = np.zeros(p)
        coverage_unad = np.zeros(p)
        nerr = 0.
        true_val = true_beta[active]
        active_set = [i for i in range(p) if active[i]]
        active_ind = np.zeros(p)
        active_ind[active_set] = 1

        list_results = []
        list_results.append(true_beta)
        list_results.append(active_ind)

        ad_lower_credible = np.zeros(p)
        ad_upper_credible = np.zeros(p)

        unad_lower_credible = np.zeros(p)
        unad_upper_credible = np.zeros(p)

        ad_mean = np.zeros(p)
        unad_mean = np.zeros(p)

        if nactive > 1:
            try:
                for l in range(nactive):
                    if (adjusted_intervals[0, l] <= true_val[l]) and (true_val[l] <= adjusted_intervals[1, l]):
                        coverage_ad[active_set[l]] += 1
                    if (unadjusted_intervals[0, l] <= true_val[l]) and (true_val[l] <= unadjusted_intervals[1, l]):
                        coverage_unad[active_set[l]] += 1
                    ad_lower_credible[active_set[l]] = adjusted_intervals[0, l]
                    ad_upper_credible[active_set[l]] = adjusted_intervals[1, l]
                    unad_lower_credible[active_set[l]] = unadjusted_intervals[0, l]
                    unad_upper_credible[active_set[l]] = unadjusted_intervals[1, l]
                    ad_mean[active_set[l]] = selective_mean[l]
                    unad_mean[active_set[l]] = post_mean[l]

            except ValueError:
                nerr += 1
                print('ignore iteration raising ValueError')

            list_results.append(coverage_ad)
            list_results.append(coverage_unad)

            ngrid = 1000
            quantiles = np.zeros((ngrid, nactive))
            for i in range(ngrid):
                quantiles[i, :] = np.percentile(samples, (i * 100.) / ngrid, axis=0)

            index_grid = np.argmin(np.abs(quantiles - np.zeros((ngrid, nactive))), axis=0)
            p_value = 2 * np.minimum(np.true_divide(index_grid, ngrid), 1. - np.true_divide(index_grid, ngrid))
            p_BH = BH_q(p_value, bh_level)

            D_BH = np.zeros(p)

            if p_BH is not None:
                for indx in p_BH[1]:
                    D_BH[active_set[indx]] = 1

            list_results.append(D_BH)

            list_results.append(ad_lower_credible)
            list_results.append(ad_upper_credible)
            list_results.append(unad_lower_credible)
            list_results.append(unad_upper_credible)

            list_results.append(ad_mean)
{'fraction_non_alphanumeric': 0.09959388899632567, 'fraction_numerical': 0.012086637014117193, 'mean_word_length': 2.0163312919218432, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 32, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '1364626', 'n_tokens_mistral': 3425, 'n_tokens_neox': 3148, 'n_words': 601}