def rod_nD(X, parallel, gm=None, median=None, data_scaler=None,
           angles_scalers1=None, angles_scalers2=None):
    """
    Find ROD overall scores when Data is higher than 3D:
      # scale dataset using Robust Scaler
      # decompose the full space into a combinations of 3D subspaces,
      # Apply ROD on each combination,
      # squish scores per subspace, so we compare apples to apples,
      # calculate average of ROD scores of all subspaces per observation.
    Note that if gm, data_scaler, angles_scalers1, angles_scalers2 are None,
    that means it is a `fit()` process and they will be calculated and returned
    to the class to be saved for future prediction. Otherwise, if they are not None,
    then it is a prediction process.

    Parameters
    ----------
    X : array-like, data points
    parallel: bool, True runs the algorithm in parallel
    gm: list (default=None), the geometric median
    median: list (default=None), MAD medians
    data_scaler: obj (default=None), RobustScaler of data
    angles_scalers1: list (default=None), MinMaxScalers of Angles group 1
    angles_scalers2: list (default=None), MinMaxScalers of Angles group 2

    Returns
    -------
    ROD decision scores, gm, median, data_scaler, angles_scalers1, angles_scalers2
    """
    if data_scaler is None:  # for fitting
        data_scaler = RobustScaler()
        X = data_scaler.fit_transform(X)
    else:  # for prediction
        X = data_scaler.transform(X)
    dim = X.shape[1]
    all_subspaces = [X[:, _com] for _com in com(range(dim), 3)]
    all_gms = [None] * len(all_subspaces) if gm is None else gm
    all_meds = [None] * len(all_subspaces) if median is None else median
    all_angles_scalers1 = [None] * len(
        all_subspaces) if angles_scalers1 is None else angles_scalers1
    all_angles_scalers2 = [None] * len(
        all_subspaces) if angles_scalers2 is None else angles_scalers2
    if parallel:
        p = Pool(multiprocessing.cpu_count())
        args = [[a, b, c, d, e] for a, b, c, d, e in
                zip(all_subspaces, all_gms, all_meds,
                    all_angles_scalers1, all_angles_scalers2)]
        results = p.starmap(process_sub, args)
        subspaces_scores, gm, median, angles_scalers1, angles_scalers2 = [], [], [], [], []
        for res in results:
            subspaces_scores.append(list(res[0]))
            gm.append(res[1])
            median.append(res[2])
            angles_scalers1.append(res[3])
            angles_scalers2.append(res[4])
        scores = np.average(np.array(subspaces_scores).T, axis=1).reshape(-1)
        p.close()
        p.join()
        return scores, gm, median, data_scaler, angles_scalers1, angles_scalers2
    subspaces_scores, gm, median, angles_scalers1, angles_scalers2 = [], [], [], [], []
    for subspace, _gm, med, ang_s1, ang_s2 in zip(all_subspaces, all_gms,
                                                  all_meds,
                                                  all_angles_scalers1,
                                                  all_angles_scalers2):
        scores_, gm_, med_, ang_s1_, ang_s2_ = process_sub(subspace=subspace,
                                                           gm=_gm, median=med,
                                                           scaler1=ang_s1,
                                                           scaler2=ang_s2)
        subspaces_scores.append(scores_)
        gm.append(gm_)
        median.append(med_)
        angles_scalers1.append(ang_s1_)
        angles_scalers2.append(ang_s2_)
    scores = np.average(np.array(subspaces_scores).T, axis=1).reshape(-1)
    return scores, gm, median, data_scaler, angles_scalers1, angles_scalers2

    Find ROD overall scores when Data is higher than 3D:
      # scale dataset using Robust Scaler
      # decompose the full space into a combinations of 3D subspaces,
      # Apply ROD on each combination,
      # squish scores per subspace, so we compare apples to apples,
      # calculate average of ROD scores of all subspaces per observation.
    Note that if gm, data_scaler, angles_scalers1, angles_scalers2 are None,
    that means it is a `fit()` process and they will be calculated and returned
    to the class to be saved for future prediction. Otherwise, if they are not None,
    then it is a prediction process.

    Parameters
    ----------
    X : array-like, data points
    parallel: bool, True runs the algorithm in parallel
    gm: list (default=None), the geometric median
    median: list (default=None), MAD medians
    data_scaler: obj (default=None), RobustScaler of data
    angles_scalers1: list (default=None), MinMaxScalers of Angles group 1
    angles_scalers2: list (default=None), MinMaxScalers of Angles group 2

    Returns
    -------
    ROD decision scores, gm, median, data_scaler, angles_scalers1, angles_scalers2
    
pyod/models/rod.py