here is my code:

"""
import numpy as np
import pandas as pd
from numpy.linalg import inv
from statsmodels.api import OLS, add_constant
from statsmodels.tools import add_constant
from scipy.stats import norm
import warnings

def mr_presso(
    BetaOutcome,
    BetaExposure,
    SdOutcome,
    SdExposure,
    data,
    OUTLIERtest=False,
    DISTORTIONtest=False,
    SignifThreshold=0.05,
    NbDistribution=1000,
    seed=None
):
    """
    Python implementation of the MR-PRESSO method for Mendelian Randomization.

    Parameters:
    - BetaOutcome: str, name of the outcome beta column
    - BetaExposure: list of str, names of the exposure beta columns
    - SdOutcome: str, name of the outcome standard deviation column
    - SdExposure: list of str, names of the exposure standard deviation columns
    - data: pandas DataFrame, input data containing the necessary columns
    - OUTLIERtest: bool, whether to perform the outlier test
    - DISTORTIONtest: bool, whether to perform the distortion test
    - SignifThreshold: float, significance threshold for tests
    - NbDistribution: int, number of distributions for empirical p-values
    - seed: int or None, random seed for reproducibility

    Returns:
    - res: dict, containing MR-PRESSO results
    """
    if seed is not None:
        np.random.seed(seed)
    
    # Input validation
    if SignifThreshold > 1:
        raise ValueError("The significance threshold cannot be greater than 1")
    
    if len(BetaExposure) != len(SdExposure):
        raise ValueError("BetaExposure and SdExposure must have the same number of elements")
    
    if not isinstance(data, pd.DataFrame):
        raise TypeError("data must be a pandas DataFrame")
    
    # Check for required columns
    required_columns = [BetaOutcome] + BetaExposure + [SdOutcome] + SdExposure
    missing_columns = set(required_columns) - set(data.columns)
    if missing_columns:
        raise ValueError(f"Missing columns in data: {missing_columns}")
    
    # Select and clean data
    data = data[required_columns].dropna().copy()
    
    # Scale BetaOutcome and BetaExposure by the sign of the first BetaExposure
    first_exposure_sign = np.sign(data[BetaExposure[0]])
    data[BetaOutcome] *= first_exposure_sign
    for exp in BetaExposure:
        data[exp] *= first_exposure_sign
    
    # Calculate weights
    data['Weights'] = 1 / (data[SdOutcome] ** 2)
    
    # Check number of instrumental variables
    if data.shape[0] <= len(BetaExposure) + 2:
        raise ValueError("Not enough instrumental variables")
    
    if data.shape[0] >= NbDistribution:
        raise ValueError("Not enough elements to compute empirical P-values, increase NbDistribution")
    
    # Define matrix power function equivalent to R's "%^%"
    def matrix_power_eigen(matrix, power):
        eigenvalues, eigenvectors = np.linalg.eigh(matrix)
        # Handle small eigenvalues to prevent numerical issues
        eigenvalues = np.where(np.abs(eigenvalues) < 1e-10, 0, eigenvalues)
        D_power = np.diag(eigenvalues ** power)
        return eigenvectors @ D_power @ inv(eigenvectors)
    
    # Function to get RSS with Leave-One-Out (LOO)
    def getRSS_LOO(BetaOutcome, BetaExposure, data, returnIV):
        dataW = data[[BetaOutcome] + BetaExposure].values * np.sqrt(data['Weights'].values[:, np.newaxis])
        X = dataW[:, 1:]
        Y = dataW[:, 0]
        n = X.shape[0]
        p = X.shape[1]
        
        RSS = 0
        CausalEstimate_LOO = np.zeros(n)
        
        # Precompute (X'X)^-1 X'Y for each LOO sample
        for i in range(n):
            X_loo = X[~np.eye(n, dtype=bool)[i]]
            Y_loo = Y[~np.eye(n, dtype=bool)[i]]
            try:
                XtX_inv = inv(X_loo.T @ X_loo)
            except np.linalg.LinAlgError:
                XtX_inv = matrix_power_eigen(X_loo.T @ X_loo, -1)
            beta_loo = XtX_inv @ X_loo.T @ Y_loo
            CausalEstimate_LOO[i] = beta_loo
            if len(BetaExposure) == 1:
                RSS += (Y[i] - X[i, :] * beta_loo) ** 2
            else:
                RSS += (Y[i] - np.dot(X[i, :], beta_loo)) ** 2
        
        RSS = np.sum(RSS)
        
        if returnIV:
            return [RSS, CausalEstimate_LOO]
        else:
            return RSS
    
    # Function to get random data for permutations
    def getRandomData(BetaOutcome, BetaExposure, SdOutcome, SdExposure, data):
        data_random = data.copy()
        for exp, sd_exp in zip(BetaExposure, SdExposure):
            data_random[exp] = np.random.normal(loc=data_random[exp], scale=data_random[sd_exp])
        data_random[BetaOutcome] = np.random.normal(
            loc=0,  # Under null hypothesis
            scale=data_random[SdOutcome]
        )
        return data_random
    
    # Step 1: Compute observed RSS
    RSSobs = getRSS_LOO(BetaOutcome, BetaExposure, data, returnIV=OUTLIERtest)
    
    # Step 2: Compute distribution of expected RSS
    RSSexp = []
    for _ in range(NbDistribution):
        random_data = getRandomData(BetaOutcome, BetaExposure, SdOutcome, SdExposure, data)
        rss = getRSS_LOO(BetaOutcome, BetaExposure, random_data, returnIV=OUTLIERtest)
        if OUTLIERtest:
            RSSexp.append(rss[0])
        else:
            RSSexp.append(rss)
    RSSexp = np.array(RSSexp)
    
    # Compute Global Test
    if OUTLIERtest:
        Observed_RSS = RSSobs[0]
        Pvalue_global = np.sum(RSSexp > Observed_RSS) / NbDistribution
        GlobalTest = {'RSSobs': Observed_RSS, 'Pvalue': Pvalue_global}
    else:
        Observed_RSS = RSSobs
        Pvalue_global = np.sum(RSSexp > Observed_RSS) / NbDistribution
        GlobalTest = {'RSSobs': Observed_RSS, 'Pvalue': Pvalue_global}
    
    # Step 3: Compute Outlier Test
    if GlobalTest['Pvalue'] < SignifThreshold and OUTLIERtest:
        OutlierTest_list = []
        for snv in range(data.shape[0]):
            # For each SNP, compute the difference between observed and expected
            Dif = data.iloc[snv][BetaOutcome] - data.iloc[snv][BetaExposure] * RSSobs[1][snv]
            # Since RSSexp is RSS, not per SNP, this part is unclear. Adjusting to compute p-value based on Dif
            pval = np.sum((RSSexp > Dif ** 2)) / NbDistribution
            OutlierTest_list.append({'RSSobs': Dif ** 2, 'Pvalue': pval})
        
        OutlierTest = pd.DataFrame(OutlierTest_list, index=data.index)
        # Apply Bonferroni correction
        OutlierTest['Pvalue'] = np.minimum(OutlierTest['Pvalue'] * data.shape[0], 1)
        OutlierTest.index.name = 'SNP'
    else:
        OUTLIERtest = False
    
    # Step 4: Compute Distortion Test
    if DISTORTIONtest and OUTLIERtest:
        # Identify outliers
        refOutlier = OutlierTest[OutlierTest['Pvalue'] <= SignifThreshold].index.tolist()
        
        BiasTest = {'Outliers Indices': None, 'Distortion Coefficient': None, 'Pvalue': None}
        
        if len(refOutlier) > 0:
            if len(refOutlier) < data.shape[0]:
                BiasExp = []
                for _ in range(NbDistribution):
                    # Replace outliers with random SNPs
                    indices = refOutlier + list(np.random.choice(data.index.difference(refOutlier), size=(data.shape[0] - len(refOutlier)), replace=False))
                    mod_random = OLS(data.loc[indices, BetaOutcome], data.loc[indices, BetaExposure]).fit(weights=data.loc[indices, 'Weights'])
                    BiasExp.append(mod_random.params.values)
                BiasExp = np.array(BiasExp)
                
                # Model without outliers
                mod_noOutliers = OLS(data.loc[data.index.difference(refOutlier), BetaOutcome],
                                     data.loc[data.index.difference(refOutlier), BetaExposure]
                                    ).fit(weights=data.loc[data.index.difference(refOutlier), 'Weights'])
                
                # Observed Bias
                BiasObs = (mod_random.params.values - mod_noOutliers.params.values) / np.abs(mod_noOutliers.params.values)
                
                # Expected Bias
                BiasExp = (mod_random.params.values - BiasExp) / np.abs(BiasExp)
                
                # P-value for distortion
                Pvalue_distortion = np.sum(np.abs(BiasExp) > np.abs(BiasObs)) / NbDistribution
                
                BiasTest = {
                    'Outliers Indices': refOutlier,
                    'Distortion Coefficient': 100 * BiasObs,
                    'Pvalue': Pvalue_distortion
                }
            else:
                BiasTest = {
                    'Outliers Indices': "All SNPs considered as outliers",
                    'Distortion Coefficient': np.nan,
                    'Pvalue': np.nan
                }
        else:
            BiasTest = {
                'Outliers Indices': "No significant outliers",
                'Distortion Coefficient': np.nan,
                'Pvalue': np.nan
            }
    # Step 5: Formatting the results
    GlobalTest['Pvalue'] = "<{:.2e}".format(1 / NbDistribution) if GlobalTest['Pvalue'] == 0 else GlobalTest['Pvalue']
    
    res_dict = {'Global Test': GlobalTest}
    
    if OUTLIERtest:
        # Replace 0 p-values with < threshold
        OutlierTest['Pvalue'] = OutlierTest['Pvalue'].apply(
            lambda x: "<{:.2e}".format(n / NbDistribution) if x == 0 else x
        )
        res_dict['Outlier Test'] = OutlierTest
        if DISTORTIONtest:
            if 'BiasTest' in locals():
                if BiasTest['Pvalue'] == 0:
                    BiasTest['Pvalue'] = "<{:.2e}".format(1 / NbDistribution)
                res_dict['Distortion Test'] = BiasTest
            else:
                res_dict['Distortion Test'] = {'Outliers Indices': None, 'Distortion Coefficient': None, 'Pvalue': None}
            # Warning for outlier test stability
            if (data.shape[0] / NbDistribution) > SignifThreshold:
                warnings.warn(
                    f"Outlier test unstable. The significance threshold of {SignifThreshold} "
                    f"for the outlier test is not achievable with only {NbDistribution} to compute the null distribution. "
                    f"The current precision is <{data.shape[0] / NbDistribution}. Increase NbDistribution."
                )
    
    # Original MR Results
    mod_all = OLS(data[BetaOutcome], data[BetaExposure]).fit(weights=data['Weights'])
    original_mr = pd.DataFrame({
        'Exposure': BetaExposure,
        'MR Analysis': 'Raw',
        'Causal Estimate': mod_all.params.values,
        'Sd': mod_all.bse.values,
        'T-stat': mod_all.tvalues.values,
        'P-value': mod_all.pvalues.values
    })
    
    # Outlier Corrected MR Results
    if OUTLIERtest and 'mod_noOutliers' in locals():
        outlier_corrected = pd.DataFrame({
            'Exposure': BetaExposure,
            'MR Analysis': 'Outlier-corrected',
            'Causal Estimate': mod_noOutliers.params.values,
            'Sd': mod_noOutliers.bse.values,
            'T-stat': mod_noOutliers.tvalues.values,
            'P-value': mod_noOutliers.pvalues.values
        })
    else:
        warnings.warn(
            "No outliers were identified, therefore the results for the outlier-corrected MR are set to NA"
        )
        outlier_corrected = pd.DataFrame({
            'Exposure': BetaExposure,
            'MR Analysis': 'Outlier-corrected',
            'Causal Estimate': [np.nan] * len(BetaExposure),
            'Sd': [np.nan] * len(BetaExposure),
            'T-stat': [np.nan] * len(BetaExposure),
            'P-value': [np.nan] * len(BetaExposure)
        })
    
    # Combine Original and Outlier Corrected MR Results
    MR = pd.concat([original_mr, outlier_corrected], ignore_index=True)
    
    res_dict['Main MR results'] = MR
    res_dict['MR-PRESSO results'] = res_dict.pop('Global Test')
    if 'Outlier Test' in res_dict:
        res_dict['MR-PRESSO results']['Outlier Test'] = res_dict.pop('Outlier Test')
    if 'Distortion Test' in res_dict:
        res_dict['MR-PRESSO results']['Distortion Test'] = res_dict.pop('Distortion Test')
    
    return res_dict

# Example Usage
if __name__ == "__main__":
    # Example data
    data = pd.DataFrame({
        'BetaOutcome': [0.2, 0.5, 0.3],
        'BetaExposure1': [0.1, 0.4, 0.2],
        'BetaExposure2': [0.2, 0.3, 0.1],
        'SdOutcome': [0.05, 0.07, 0.06],
        'SdExposure1': [0.02, 0.03, 0.025],
        'SdExposure2': [0.03, 0.04, 0.035]
    })
    

"""

here is the dataset variable names I have used in my work. 
    results = mr_presso(
        BetaOutcome='Y_effect',
        BetaExposure=['E1_effect', 'E2_effect'],
        SdOutcome='Y_se',
        SdExposure=['E1_se', 'E2_se'],
        data=df,
        OUTLIERtest=True,
        DISTORTIONtest=True,
        SignifThreshold=0.05,
        NbDistribution=1000,
        seed=123
    )
    
    # Displaying the results
    for key, value in results.items():
        print(f"\n{key}:\n{value}")

fix the error:
 TypeError                                 Traceback (most recent call last)
TypeError: only length-1 arrays can be converted to Python scalars

The above exception was the direct cause of the following exception:

ValueError                                Traceback (most recent call last)
<ipython-input-5-d4732ac837f5> in <cell line: 1>()
----> 1 results = mr_presso(
      2     BetaOutcome='Y_effect',
      3     BetaExposure=['E1_effect', 'E2_effect'],
      4     SdOutcome='Y_se',
      5     SdExposure=['E1_se', 'E2_se'],

1 frames
<ipython-input-4-8e426ad3fac5> in getRSS_LOO(BetaOutcome, BetaExposure, data, returnIV)
    103                 XtX_inv = matrix_power_eigen(X_loo.T @ X_loo, -1)
    104             beta_loo = XtX_inv @ X_loo.T @ Y_loo
--> 105             CausalEstimate_LOO[i] = beta_loo
    106             if len(BetaExposure) == 1:
    107                 RSS += (Y[i] - X[i, :] * beta_loo) ** 2

ValueError: setting an array element with a sequence.