I have a question, and the question is that can we implement this code in a way to be run on gpu?

"""
import numpy as np
import pandas as pd
import statsmodels.api as sm
from scipy.linalg import eig, inv

def mr_presso_single(BetaOutcome, BetaExposure, SdOutcome, SdExposure, data,
                     OUTLIERtest=False, DISTORTIONtest=False,
                     SignifThreshold=0.05, NbDistribution=1000, seed=None):
    if seed is not None:
        np.random.seed(seed)
    
    if SignifThreshold > 1:
        raise ValueError("The significance threshold cannot be greater than 1")

    if not isinstance(data, pd.DataFrame):
        raise TypeError("data must be a pandas DataFrame")
    
    # Ensure that inputs are in the correct format
    if isinstance(BetaExposure, str):
        BetaExposure = [BetaExposure]
    if isinstance(SdExposure, str):
        SdExposure = [SdExposure]
    
    if len(BetaExposure) != len(SdExposure):
        raise ValueError("BetaExposure and SdExposure must have the same number of elements")
    
    # Data preparation and validation
    required_columns = [BetaOutcome] + BetaExposure + [SdOutcome] + SdExposure
    data = data[required_columns].dropna()
    data = data.reset_index(drop=True)  # Reset the index to ensure alignment
    
    # Adjust signs based on the first exposure variable
    sign_exposure = np.sign(data[BetaExposure[0]])
    data[BetaOutcome] = data[BetaOutcome] * sign_exposure
    data[BetaExposure[0]] = data[BetaExposure[0]] * sign_exposure
    
    # Calculate weights
    data['Weights'] = 1 / data[SdOutcome] ** 2

    # Validate number of observations
    if len(data) <= len(BetaExposure) + 2:
        raise ValueError("Not enough instrumental variables")

    if len(data) >= NbDistribution:
        raise ValueError("Not enough elements to compute empirical P-values, increase NbDistribution")

    # Step 1: Compute observed RSS
    RSSobs = getRSS_LOO(BetaOutcome, BetaExposure, data, returnIV=OUTLIERtest)

    # Step 2: Compute expected RSS distribution
    randomData = [getRandomData(BetaOutcome, BetaExposure, SdOutcome, SdExposure, data)
                  for _ in range(NbDistribution)]

    RSSexp = []
    for rd in randomData:
        rss = getRSS_LOO(BetaOutcome, BetaExposure, rd, returnIV=OUTLIERtest)
        if OUTLIERtest:
            RSSexp.append(rss[0])
        else:
            RSSexp.append(rss)

    RSSexp = np.array(RSSexp)

    if OUTLIERtest:
        GlobalTest = {
            'RSSobs': RSSobs[0],
            'Pvalue': np.sum(RSSexp > RSSobs[0]) / NbDistribution
        }
    else:
        GlobalTest = {
            'RSSobs': RSSobs,
            'Pvalue': np.sum(RSSexp > RSSobs) / NbDistribution
        }

    # Step 3: Perform single IV outlier test
    if GlobalTest['Pvalue'] < SignifThreshold and OUTLIERtest:
        OutlierTest = []
        n_snps = len(data)
        for snv in range(n_snps):
            randomSNP = pd.DataFrame([rd.iloc[snv] for rd in randomData])
            if len(BetaExposure) == 1:
                Dif = data.iloc[snv][BetaOutcome] - data.iloc[snv][BetaExposure] * RSSobs[1][snv]
                Exp = randomSNP[BetaOutcome] - randomSNP[BetaExposure].values.flatten() * RSSobs[1][snv]
            else:
                Dif = data.iloc[snv][BetaOutcome] - np.sum(data.iloc[snv][BetaExposure].values * RSSobs[1][:, snv])
                Exp = randomSNP[BetaOutcome] - np.sum(randomSNP[BetaExposure].values * RSSobs[1][:, snv], axis=1)
            pval = np.sum(Exp ** 2 > Dif ** 2) / NbDistribution
            OutlierTest.append({
                'RSSobs': Dif ** 2,
                'Pvalue': pval
            })

        OutlierTest = pd.DataFrame(OutlierTest)
        OutlierTest.index = data.index

        # Bonferroni correction
        OutlierTest['Pvalue'] = np.minimum(OutlierTest['Pvalue'] * n_snps, 1)
    else:
        OUTLIERtest = False

    # Step 4: Perform distortion test
    X_all = data[BetaExposure]
    Y_all = data[BetaOutcome]
    W_all = data['Weights']
    model_all = sm.WLS(Y_all, X_all, weights=W_all)
    results_all = model_all.fit()

    if DISTORTIONtest and OUTLIERtest:
        refOutlier = OutlierTest[OutlierTest['Pvalue'] <= SignifThreshold].index.tolist()

        if len(refOutlier) > 0:
            if len(refOutlier) < len(data):
                def getRandomBias(BetaOutcome, BetaExposure, data, refOutlier):
                    n = len(data)
                    indices = refOutlier + list(
                        np.random.choice(
                            data.index.difference(refOutlier),
                            size=n - len(refOutlier),
                            replace=False
                        )
                    )
                    data_subset = data.loc[indices[:n - len(refOutlier)]]
                    X_subset = data_subset[BetaExposure]
                    Y_subset = data_subset[BetaOutcome]
                    W_subset = data_subset['Weights']
                    model_random = sm.WLS(Y_subset, X_subset, weights=W_subset)
                    results_random = model_random.fit()
                    return results_random.params

                BiasExp = []
                for _ in range(NbDistribution):
                    bias = getRandomBias(BetaOutcome, BetaExposure, data, refOutlier)
                    BiasExp.append(bias)
                BiasExp = pd.DataFrame(BiasExp)

                data_no_outliers = data.drop(refOutlier)
                X_no_outliers = data_no_outliers[BetaExposure]
                Y_no_outliers = data_no_outliers[BetaOutcome]
                W_no_outliers = data_no_outliers['Weights']
                model_no_outliers = sm.WLS(Y_no_outliers, X_no_outliers, weights=W_no_outliers)
                results_no_outliers = model_no_outliers.fit()

                BiasObs = (results_all.params - results_no_outliers.params) / np.abs(results_no_outliers.params)
                BiasExp = (results_all.params.values.reshape(1, -1) - BiasExp) / np.abs(BiasExp)

                pvalues = []
                for idx in range(len(BetaExposure)):
                    pvalue = np.sum(np.abs(BiasExp.iloc[:, idx]) > np.abs(BiasObs[idx])) / NbDistribution
                    pvalues.append(pvalue)
                
                BiasTest = {
                    'Outliers Indices': refOutlier,
                    'Distortion Coefficient (%)': 100 * BiasObs,
                    'Pvalue': pvalues
                }
            else:
                BiasTest = {
                    'Outliers Indices': "All SNPs considered as outliers",
                    'Distortion Coefficient (%)': None,
                    'Pvalue': None
                }
        else:
            BiasTest = {
                'Outliers Indices': "No significant outliers",
                'Distortion Coefficient (%)': None,
                'Pvalue': None
            }
    else:
        DISTORTIONtest = False

    # Step 5: Formatting the results
    if GlobalTest['Pvalue'] == 0:
        GlobalTest['Pvalue'] = f"<{1 / NbDistribution}"

    results = {
        'Global Test': GlobalTest
    }

    if OUTLIERtest:
        OutlierTest['Pvalue'] = OutlierTest['Pvalue'].apply(lambda x: f"<{1 / NbDistribution}" if x == 0 else x)
        results['Outlier Test'] = OutlierTest

        if DISTORTIONtest:
            BiasTest['Pvalue'] = [f"<{1 / NbDistribution}" if p == 0 else p for p in BiasTest['Pvalue']]
            results['Distortion Test'] = BiasTest

    # Main MR results
    MR_results = pd.DataFrame({
        'Exposure': BetaExposure,
        'MR Analysis': 'Raw',
        'Causal Estimate': results_all.params.values,
        'Std Error': results_all.bse.values,
        'T-stat': results_all.tvalues,
        'P-value': results_all.pvalues
    })

    if DISTORTIONtest and 'Outliers Indices' in BiasTest and isinstance(BiasTest['Outliers Indices'], list):
        data_no_outliers = data.drop(BiasTest['Outliers Indices'])
        X_no_outliers = data_no_outliers[BetaExposure]
        Y_no_outliers = data_no_outliers[BetaOutcome]
        W_no_outliers = data_no_outliers['Weights']
        model_no_outliers = sm.WLS(Y_no_outliers, X_no_outliers, weights=W_no_outliers)
        results_no_outliers = model_no_outliers.fit()

        MR_results_no_outliers = pd.DataFrame({
            'Exposure': BetaExposure,
            'MR Analysis': 'Outlier-corrected',
            'Causal Estimate': results_no_outliers.params.values,
            'Std Error': results_no_outliers.bse.values,
            'T-stat': results_no_outliers.tvalues,
            'P-value': results_no_outliers.pvalues
        })
    else:
        MR_results_no_outliers = pd.DataFrame({
            'Exposure': BetaExposure,
            'MR Analysis': 'Outlier-corrected',
            'Causal Estimate': [np.nan] * len(BetaExposure),
            'Std Error': [np.nan] * len(BetaExposure),
            'T-stat': [np.nan] * len(BetaExposure),
            'P-value': [np.nan] * len(BetaExposure)
        })

    MR_main_results = pd.concat([MR_results, MR_results_no_outliers], ignore_index=True)
    results['Main MR results'] = MR_main_results

    # Warning if necessary
    if OUTLIERtest and (len(data) / NbDistribution > SignifThreshold):
        print(f"Warning: Outlier test unstable. The significance threshold of {SignifThreshold} for the outlier test "
              f"is not achievable with only {NbDistribution} simulations to compute the null distribution. "
              f"The current precision is <{len(data) / NbDistribution}. Increase NbDistribution.")

    return results
"""