Python rpy2.robjects.packages.importr() Examples

def install_matchit(self):
        from rpy2.robjects.packages import importr
        import rpy2.robjects.packages as rpackages
        from rpy2.robjects.vectors import StrVector
        import rpy2.robjects as robjects

        package_names = ["MatchIt"]

        names_to_install = [x for x in package_names if not rpackages.isinstalled(x)]
        if len(names_to_install) > 0:
            robjects.r.options(download_file_method='curl')
            utils = rpackages.importr('utils')
            utils.chooseCRANmirror(ind=0)
            utils.chooseCRANmirror(ind=0)
            utils.install_packages(StrVector(names_to_install))

        return importr("MatchIt") 

def pi0est(p, lambda_qvalue=None):
    """Wrapper for qvalue::pi0est"""
    qvalue = importr("qvalue")
    rp = rpy2.robjects.vectors.FloatVector(p)
    # with suppress_stdout():
    if lambda_qvalue is None:
        pi0res = qvalue.pi0est(rp)
    else:
        if not isinstance(lambda_qvalue, Iterable):
            lambda_qvalue = [lambda_qvalue]
        rlambda = rpy2.robjects.vectors.FloatVector(lambda_qvalue)
        pi0res = qvalue.pi0est(rp, rlambda)
    pi0 = np.array(pi0res.rx2('pi0'))[0]
    pi0_lambda = np.array(pi0res.rx2('pi0.lambda'))
    lambda_vec = np.array(pi0res.rx2('lambda'))
    pi0_smooth = np.array(pi0res.rx2('pi0.smooth'))
    return pi0, pi0_lambda, lambda_vec, pi0_smooth 

def install_bart(self):
        import rpy2.robjects.packages as rpackages
        from rpy2.robjects.packages import importr
        from rpy2.robjects.vectors import StrVector
        import rpy2.robjects as robjects

        robjects.r.options(download_file_method='curl')

        # install.packages("rJava")
        rj = importr("rJava", robject_translations={'.env': 'rj_env'})
        rj._jinit(parameters="-Xmx16g", force_init=True)
        print("rJava heap size is", np.array(rj._jcall(rj._jnew("java/lang/Runtime"), "J", "maxMemory"))[0] / 1e9,
              "GB.", file=sys.stderr)

        package_names = ["bartMachine"]
        utils = rpackages.importr('utils')
        utils.chooseCRANmirror(ind=0)
        utils.chooseCRANmirror(ind=0)

        names_to_install = [x for x in package_names if not rpackages.isinstalled(x)]
        if len(names_to_install) > 0:
            utils.install_packages(StrVector(names_to_install))

        return importr("bartMachine") 

def install_grf(self):
        from rpy2.robjects.packages import importr
        import rpy2.robjects.packages as rpackages
        from rpy2.robjects.vectors import StrVector
        import rpy2.robjects as robjects

        robjects.r.options(download_file_method='curl')

        package_names = ["grf"]
        utils = rpackages.importr('utils')
        utils.chooseCRANmirror(ind=0)
        utils.chooseCRANmirror(ind=0)

        names_to_install = [x for x in package_names if not rpackages.isinstalled(x)]
        if len(names_to_install) > 0:
            utils.install_packages(StrVector(names_to_install))

        return importr("grf") 

def grdevices():
    return importr('grDevices') 

def test_sumr_compile(queue, func, n, setup_func):
    array, module = setup_func("R")
    r_loopsum = """
#-- purer_sum-begin
function(x)
{
  total = 0;
  for (elt in x) {
    total <- total + elt
  } 
}
#-- purer_sum-end
"""
    
    rdo_loopsum = """
f <- %s
for (i in 1:%i) {
  res <- f(x)
}
"""
    compiler = ro.packages.importr("compiler")
    rcode = rdo_loopsum %("compiler::cmpfun("+r_loopsum+")", n)
    time_beg = time.time()
    #print(rcode)
    module.r(rcode)
    time_end = time.time()
    queue.put(time_end - time_beg) 

def importr(packname, newname = None, verbose = False):
    """ Wrapper around rpy2.robjects.packages.importr, 
    adding the following feature(s):
    
    - package instance added to the pseudo-module 'packages'

    """

    assert isinstance(packname, str)
    packinstance = _importr(packname, on_conflict = 'warn')

    # fix the package name (dots possible in R package names)
    if newname is None:
        newname = packname.replace('.', '_')

    Packages().__dict__[newname] = packinstance

    ## Currently too slow for a serious usage: R's introspection 
    ## of S4 classes is not fast enough
    # d = {}
    # for cn in methods.get_classnames(packname):
    #     class AutoS4(RS4):
    #         __metaclass__ = methods.RS4Auto_Type
    #         __rpackagename__ = packname
    #         __rname__ = cn
    #     newcn = cn.replace('.', '_')
    #     d[newcn] = AutoS4
    # S4Classes().__dict__[newname] = d 
    
    return packinstance 

def set_R_plotting_device(self, device):
        """
        Set which device R should use to produce plots.
        If device == 'svg' then the package 'Cairo'
        must be installed. Because Cairo forces "onefile=TRUE",
        it is not posible to include multiple plots per cell.

        Parameters
        ----------

        device : ['png', 'X11', 'svg']
            Device to be used for plotting.
            Currently only "png" and "X11" are supported,
            with 'png' and 'svg' being most useful in the notebook,
            and 'X11' allowing interactive plots in the terminal.

        """
        device = device.strip()
        if device not in ['png', 'X11', 'svg']:
            raise ValueError(
                "device must be one of ['png', 'X11' 'svg'], got '%s'", device)
        if device == 'svg':
            try:
                self.cairo = rpacks.importr('Cairo')
            except ri.embedded.RRuntimeError as rre:
                if rpacks.isinstalled('Cairo'):
                    msg = ('An error occurred when trying to load the ' +
                           'R package Cairo\'\n%s' % str(rre))
                else:
                    msg = textwrap.dedent("""
                    The R package 'Cairo' is required but it does not appear
                    to be installed/available. Try:

                    import rpy2.robjects.packages as rpacks
                    utils = rpacks.importr('utils')
                    utils.chooseCRANmirror(ind=1)
                    utils.install_packages('Cairo')
                    """)
                raise RInterpreterError(msg)
        self.device = device 

def dtwWrapper(data, rows, columns, k):
    '''
    wrapper function for dynamic time warping.
    includes use of exponential adaptive tuning function
    with temporal correlation if k > 0
    '''

    # not explicitly called, but needs to be in R environment
    DTW = importr("dtw")

    # create a data frame of zeros of size number of ids x number of ids
    # fill it with the calculated distance metric for each pair wise comparison

    df_ = pd.DataFrame(index=rows,
                       columns=columns)
    df_ = df_.fillna(0.0).astype(np.float64)

    # fill the array with dtw-distance values
    pandas2ri.activate()

    for i in rows:
        E.info("DTW %s" % i)
        for j in columns:
            series1 = data.loc[i].values.tolist()
            series2 = data.loc[j].values.tolist()
            DTW_value = (R.dtw(series1,
                               series2)).rx('distance')[0][0]
            cort_value = temporalCorrelate(series1, series2)
            tuned_value = adaptiveTune(cort_value, k)
            time_dist = DTW_value * tuned_value
            df_.loc[i][j] = float(time_dist)
            df_[j][i] = float(time_dist)

    return df_ 

def fit_and_evaluate(x,z,t,y,df):
    '''
    Fit and evaluate non-parametric regression using  Darolles, Fan, Florens and Renault (2011)

    Implemented in the `np` package in R.

    See [the np package documation](https://cran.r-project.org/web/packages/np/np.pdf) for details.
    '''
    npr=importr('np')
    y_R = robjects.FloatVector(list(y.flatten()))
    (x_eval, t_eval), y_true = test_points(df, 10000)
    mod = npr.npregiv(y_R, t, z, x=x, zeval=t_eval, xeval=x_eval,
                    method="Tikhonov", p=0, optim_method ="BFGS")
    return ((y_true - to_array(mod.rx2('phi.eval')))**2).mean() 

def qvalue(p, lambda_qvalue=None):
    """Wrapper for qvalue::qvalue"""
    qvalue = importr("qvalue")
    rp = rpy2.robjects.vectors.FloatVector(p)
    if lambda_qvalue is None:
        q = qvalue.qvalue(rp)
    else:
        if not isinstance(lambda_qvalue, Iterable):
            lambda_qvalue = [lambda_qvalue]
        rlambda = rpy2.robjects.vectors.FloatVector(lambda_qvalue)
        q = qvalue.qvalue(rp, **{'lambda':rlambda})
    qval = np.array(q.rx2('qvalues'))
    pi0 = np.array(q.rx2('pi0'))[0]
    return qval, pi0 

def importable(self) -> bool:
        try:
            import rpy2.robjects.packages as rpackages

            rpackages.importr(self.package_name, robject_translations={".env": "env"})

            return True
        except Exception:
            return False 

def graphics():
    return importr('graphics') 

def stats():
    return importr('stats') 

def wgcna():
    return importr('WGCNA') 

def base():
    return importr('base') 

def installRCausalImpact():
    from rpy2.robjects.packages import importr
    utils = importr('utils')
    utils.install_packages('CausalImpact') 

def installRSynth():
    from rpy2.robjects.packages import importr
    utils = importr('utils')
    utils.install_packages('Synth') 

def format_code(self, code: str, notebook: bool, **options) -> str:
        import rpy2.robjects.packages as rpackages

        styler_r = rpackages.importr(self.package_name)
        formatted_code = styler_r.style_text(
            code, **self._transform_options(styler_r, options)
        )
        return "\n".join(formatted_code) 

def importable(self) -> bool:
        try:
            import rpy2.robjects.packages as rpackages

            rpackages.importr(self.package_name)

            return True
        except Exception:
            return False 

def format_code(self, code: str, notebook: bool, **options) -> str:
        import rpy2.robjects.packages as rpackages

        format_r = rpackages.importr(
            self.package_name, robject_translations={".env": "env"}
        )
        formatted_code = format_r.tidy_source(text=code, output=False, **options)
        return "\n".join(formatted_code[0]) 

def principal_curve(data, basis="pca", n_comps=4, clusters_list=None, copy=False):
    """Computes the principal curve
    Arguments
    ---------
    data: :class:`~anndata.AnnData`
        Annotated data matrix.
    basis: `str` (default: `'pca'`)
        Basis to use for computing the principal curve.
    n_comps: `int` (default: 4)
        Number of pricipal components to be used.
    copy: `bool`, (default: `False`)
        Return a copy instead of writing to adata.
    Returns
    -------
    Returns or updates `adata` with the attributes
    principal_curve: `.uns`
        dictionary containing `projections`, `ixsort` and `arclength`
    """
    adata = data.copy() if copy else data
    import rpy2.robjects as robjects
    from rpy2.robjects.packages import importr

    if clusters_list is not None:
        cell_subset = np.array(
            [label in clusters_list for label in adata.obs["clusters"]]
        )
        X_emb = adata[cell_subset].obsm[f"X_{basis}"][:, :n_comps]
    else:
        cell_subset = None
        X_emb = adata.obsm[f"X_{basis}"][:, :n_comps]

    n_obs, n_dim = X_emb.shape

    # convert array to R matrix
    xvec = robjects.FloatVector(X_emb.T.reshape((X_emb.size)))
    X_R = robjects.r.matrix(xvec, nrow=n_obs, ncol=n_dim)

    fit = importr("princurve").principal_curve(X_R)

    adata.uns["principal_curve"] = dict()
    adata.uns["principal_curve"]["ixsort"] = ixsort = np.array(fit[1]) - 1
    adata.uns["principal_curve"]["projections"] = np.array(fit[0])[ixsort]
    adata.uns["principal_curve"]["arclength"] = np.array(fit[2])
    adata.uns["principal_curve"]["cell_subset"] = cell_subset

    return adata if copy else None 

def fitRCausalImpact(Y_pre, Y_post, treated_units):
        import rpy2.robjects as ro #ro.r is the R instace
        #Automatically convert numpy arrays to R vector/matrix
        from rpy2.robjects import numpy2ri
        numpy2ri.activate()

        control_units = [u for u in range(Y_pre.shape[0]) if u not in treated_units]

        #np.matrix is not automatically converted so use ndarray
        if type(Y_pre).__name__=="matrix":
            Y_pre = Y_pre.A
        if type(Y_post).__name__=="matrix":
            Y_post = Y_post.A

        Y = np.hstack((Y_pre, Y_post))
        Y_sc = np.full(Y.shape, np.nan)
        Y_c = Y[control_units,:]
        T0 = Y_pre.shape[1]
        N,T = Y.shape

        try:
            CausalImpact = importr('CausalImpact')
        except:
            raise RuntimeError("Need the 'CausalImpact' package loaded in the rpy2 R environment. Use test.test_simulation.installRCausalImpact")
        
        for unit in range(N): #
            if unit in treated_units:
                data = np.hstack((Y[unit,:].T, Y_c.T))
            else:
                data = np.hstack((Y[unit,:].T, np.delete(Y_c, unit, 0).T))
            r_data = ro.r.matrix(data, nrow=df.shape[0], ncol=df.shape[1])

            r_casaulimpact_out = CausalImpact.CausalImpact(data=r_data, pre_period=ro.IntVector([1, T0]), post_period=ro.IntVector([T0+1, T]))
            #can't seem to get weights from impact$model$bsts.model
            #r_summary = r_casaulimpact_out[r_casaulimpact_out.names.index('summary')]
            #te = r_summary[r_summary.names.index('AbsEffect')]
            #p = r_summary[r_summary.names.index('p')]
            #te_ci_l = r_summary[r_summary.names.index('AbsEffect.lower')]
            #te_ci_u = r_summary[r_summary.names.index('AbsEffect.upper')]
            point_pred = np.array(r_casaulimpact_out[r_casaulimpact_out.names.index('series')])[:,2] #'point.pred' from 'zoo' object
            Y_sc[unit,:] = point_pred

        class RCausalImpact:
            def __self__(self, Y_sc):
                self.Y_sc = Y_sc

        return RCausalImpact(Y_sc) 

def __init__(
        self,
        freq: str,
        prediction_length: int,
        method_name: str = "ets",
        period: int = None,
        trunc_length: Optional[int] = None,
        params: Optional[Dict] = None,
    ) -> None:
        super().__init__(freq=freq, prediction_length=prediction_length)

        try:
            from rpy2 import robjects, rinterface
            import rpy2.robjects.packages as rpackages
            from rpy2.rinterface import RRuntimeError
        except ImportError as e:
            raise ImportError(str(e) + USAGE_MESSAGE) from e

        self._robjects = robjects
        self._rinterface = rinterface
        self._rinterface.initr()
        self._rpackages = rpackages

        this_dir = os.path.dirname(os.path.realpath(__file__))
        this_dir = this_dir.replace("\\", "/")  # for windows
        r_files = [
            n[:-2] for n in os.listdir(f"{this_dir}/R/") if n[-2:] == ".R"
        ]

        for n in r_files:
            try:
                path = Path(this_dir, "R", f"{n}.R")
                robjects.r(f'source("{path}")'.replace("\\", "\\\\"))
            except RRuntimeError as er:
                raise RRuntimeError(str(er) + USAGE_MESSAGE) from er

        supported_methods = ["ets", "arima", "tbats", "croston", "mlp"]
        assert (
            method_name in supported_methods
        ), f"method {method_name} is not supported please use one of {supported_methods}"

        self.method_name = method_name

        self._stats_pkg = rpackages.importr("stats")
        self._r_method = robjects.r[method_name]

        self.prediction_length = prediction_length
        self.freq = freq
        self.period = period if period is not None else get_seasonality(freq)
        self.trunc_length = trunc_length

        self.params = {
            "prediction_length": self.prediction_length,
            "output_types": ["samples"],
            "frequency": self.period,
        }
        if params is not None:
            self.params.update(params) 

def M3Dstat2pvalue(df, columns, pair):

    stats = importr('stats')

    def meltAndPivot(df, columns):
        melt = pd.melt(df, id_vars=columns)
        melt['value'] = melt['value'].astype(float)
        pivot = pd.DataFrame(melt.pivot_table(
            index=columns, values="value", aggfunc=np.mean))
        return pivot

    def sampler(M3Dstat, vector, n):
        rand = np.random.choice(vector, size=n, replace=True, )
        return len(rand[rand > M3Dstat]) / n

    def addPvalues(df, within_df, repeat):
        df['p_value'] = df['value'].apply(
            sampler, args=(within_df['value'], repeat))
        df['p_value_adj'] = stats.p_adjust(FloatVector(
            df['p_value']), method='BH')

    E.info("dimensions:")
    E.info(df.shape)

    M3D_within_col = []
    M3D_between_col = []

    # create regexes from pair in both directions to find columns
    for p in [pair, pair[::-1]]:
        regex = "%s-\d__vs__%s-\d" % (p[0], p[0])
        M3D_within_col.extend([col for col in df.columns
                               if re.search(regex, col)])
        regex = "%s-\d__vs__%s-\d" % (p[0], p[1])
        M3D_between_col.extend([col for col in df.columns
                                if re.search(regex, col)])

    within = copy.copy(columns)
    within.extend(M3D_within_col)
    between = copy.copy(columns)
    between.extend(M3D_between_col)

    between_pivot = meltAndPivot(df.ix[:, between], columns)
    within_pivot = meltAndPivot(df.ix[:, within], columns)

    addPvalues(between_pivot, within_pivot, 100000.0)
    addPvalues(within_pivot, within_pivot, 100000.0)

    between_pivot.reset_index(inplace=True)
    within_pivot.reset_index(inplace=True)

    return between_pivot, within_pivot