Repository: zhiyanov/DCoNA (language: python)

Issue:
Create build.yml



Discussion / hints:

Fix patch:
diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml
new file mode 100644
index 0000000..740b766
--- /dev/null
+++ b/.github/workflows/build.yml
@@ -0,0 +1,98 @@
+name: Build
+
+on: 
+    push:
+        branches: [ shell ]
+    pull_request:
+        branches: [ shell ]
+    workflow_dispatch:
+
+
+jobs:
+  build_wheels:
+    name: Build wheel for ${{ matrix.python }}-${{ matrix.buildplat[1] }}
+    runs-on: ${{ matrix.buildplat[0] }}
+    strategy:
+      # Ensure that a wheel builder finishes even if another fails
+      fail-fast: false
+      matrix:
+        buildplat:
+          - [ubuntu-latest, manylinux_x86_64]
+          - [macos-latest, macosx_x86_64]
+          - [macos-latest, macosx_arm64]
+          - [macos-latest, macosx_universal2]
+          - [windows-latest, win_amd64]
+          - [windows-latest, win32]
+          
+        python: ["cp38", "cp39", "cp310", "pp38"]
+        
+        exclude:
+          # Don't build PyPy 32-bit windows
+          - buildplat: [windows-latest, win32]
+            python: "pp38"
+          - buildplat: [macos-latest, macosx_arm64]
+            python: "pp38"
+          - buildplat: [macos-latest, macosx_universal2]
+            python: "pp38"
+          - buildplat: [windows-latest, win32]
+            python: "cp310"
+
+    steps:
+      - uses: actions/checkout@v2
+
+      - name: Setup Python
+        uses: actions/setup-python@v2
+        with:
+          python-version: '3.9'  # update once build dependencies are available
+
+      - name: Build wheels
+        uses: pypa/cibuildwheel@v2.5.0
+        env:
+          CIBW_BUILD: ${{ matrix.python }}-${{ matrix.buildplat[1] }}
+          CIBW_ARCHS: all
+          CIBW_BEFORE_BUILD_WINDOWS: "pip install delvewheel"
+          CIBW_REPAIR_WHEEL_COMMAND_WINDOWS: "delvewheel repair -w {dest_dir} {wheel}"
+          
+      - uses: actions/upload-artifact@v3
+        with:
+          name: wheels_artifact
+          path: ./wheelhouse/*.whl
+          
+          
+  build_sdist:
+    name: Build source distribution
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v2
+
+      - name: Build sdist
+        run: pipx run build --sdist
+            
+      - uses: actions/upload-artifact@v3
+        with:
+          name: source_artifact
+          path: dist/*.tar.gz
+
+
+  publish:
+    name: Publish on PyPI
+    runs-on: ubuntu-latest
+    needs: [build_wheels, build_sdist]
+    steps:
+      - uses: actions/download-artifact@v3
+        with:
+          name: wheels_artifact
+          path: dist
+      - uses: actions/download-artifact@v3
+        with:
+          name: source_artifact
+          path: dist
+          
+      - name: Publish package
+        #uses: pypa/gh-action-pypi-publish@27b31702a0e7fc50959f5ad993c78deac1bdfc29
+        uses: pypa/gh-action-pypi-publish@master
+        with:
+          user: __token__
+          password: ${{ secrets.TEST_PYPI_API_TOKEN }}
+          repository_url: https://test.pypi.org/legacy/
+          skip_existing: true
\ No newline at end of file
diff --git a/.gitignore b/.gitignore
index 2a05372..80b0c65 100644
--- a/.gitignore
+++ b/.gitignore
@@ -14,6 +14,7 @@ dist/
 downloads/
 eggs/
 .eggs/
+wheelhouse/
 # Not default behaviour
 # lib/
 lib64/
@@ -143,5 +144,5 @@ scripts/
 # Build dir
 build/
 
-# .swp files: they are being created during running vim
-.swp
+# vim
+*.swp
diff --git a/.gitmodules b/.gitmodules
index dd664e9..4f7ea07 100644
--- a/.gitmodules
+++ b/.gitmodules
@@ -1,4 +1,4 @@
 [submodule "native/extern/pybind11"]
-	path = native/extern/pybind11
+    path = dcona/native/extern/pybind11
     url = https://github.com/pybind/pybind11.git
     branch = stable
diff --git a/__init__.py b/__init__.py
deleted file mode 100644
index e69de29..0000000
diff --git a/build.sh b/build.sh
index c86bc88..df63077 100755
--- a/build.sh
+++ b/build.sh
@@ -2,8 +2,8 @@
 
 mkdir build
 cd build
-cmake ../native
+cmake ../dcona/native
 make
-cp *.so ../core/extern/
+cp *.so ../dcona/core/extern/
 cd ..
 rm -rf build
diff --git a/core/__init__.py b/core/__init__.py
deleted file mode 100644
index e69de29..0000000
diff --git a/dcona/__init__.py b/dcona/__init__.py
new file mode 100644
index 0000000..364e563
--- /dev/null
+++ b/dcona/__init__.py
@@ -0,0 +1,3 @@
+from .lib.ztest import ztest
+from .lib.zscore import zscore
+from .lib.hypergeom import hypergeom
diff --git a/dcona/__main__.py b/dcona/__main__.py
new file mode 100644
index 0000000..cffdc10
--- /dev/null
+++ b/dcona/__main__.py
@@ -0,0 +1,145 @@
+def ztest_cli(config_path):
+    import pandas as pd 
+    from . import lib
+
+    data_path, description_path, \
+    reference_group, experimental_group, \
+    correlation, alternative, score, \
+    interaction_path, \
+    repeats_number, \
+    output_dir_path, \
+    process_number = \
+    lib.utils.read_config(config_path)
+    
+    lib.dump.check_directory_existence(output_dir_path)
+    
+    data_df = pd.read_csv(data_path, sep=",", index_col=0)
+    description_df = pd.read_csv(description_path, sep=",")
+    if interaction_path:
+        interaction_df = pd.read_csv(interaction_path, sep=",")
+    else:
+        interaction_df = None
+    
+    result = lib.ztest(
+        data_df, description_df,
+        reference_group, experimental_group,
+        correlation, alternative,
+        interaction=interaction_df,
+        repeats_number=repeats_number,
+        output_dir=output_dir_path,
+        process_number=process_number
+    )
+
+    if not (result is None): 
+        path_to_file = output_dir_path.rstrip("/") + "/{}_ztest.csv".format(correlation)
+        result.to_csv(path_to_file, sep=",", index=None)
+        print(f"File saved at: {path_to_file}") 
+    
+def zscore_cli(config_path):   
+    import pandas as pd
+    from . import lib
+
+    data_path, description_path, \
+    reference_group, experimental_group, \
+    correlation, alternative, score, \
+    interaction_path, \
+    repeats_number, \
+    output_dir_path, \
+    process_number = \
+    lib.utils.read_config(config_path)
+    
+    lib.dump.check_directory_existence(output_dir_path)
+    
+    data_df = pd.read_csv(data_path, sep=",", index_col=0)
+    description_df = pd.read_csv(description_path, sep=",")
+    if interaction_path:
+        interaction_df = pd.read_csv(interaction_path, sep=",")
+    else:
+        interaction_df = None
+    
+    result = lib.zscore(
+        data_df, description_df,
+        reference_group, experimental_group,
+        correlation=correlation,
+        score=score,
+        alternative=alternative,
+        interaction=interaction_df,
+        repeats_number=repeats_number,
+        output_dir=output_dir_path,
+        process_number=process_number
+    )
+                            
+    if not (result is None): 
+        path_to_file = output_dir_path.rstrip("/") + \
+            "/{}_{}_{}_zscore.csv".format(
+                correlation, score, alternative)
+        result.to_csv(path_to_file, sep=",", index=None)
+        print(f"File saved at: {path_to_file}") 
+
+def hypergeom_cli(config_path):   
+    import pandas as pd
+    from . import lib
+
+    data_path, description_path, \
+    reference_group, experimental_group, \
+    correlation, alternative, score, \
+    interaction_path, \
+    repeats_number, \
+    output_dir_path, \
+    process_number = \
+    lib.utils.read_config(config_path)
+    
+    lib.dump.check_directory_existence(output_dir_path)
+
+    path_to_file = output_dir.rstrip("/") + f"/{correlation}_ztest.csv"
+    ztest_df = pd.read_csv(path_to_file, sep=",")
+    
+    if interaction_path:
+        oriented = True
+    else:
+        oriented = False
+
+    result = lib.hypergeom(
+        ztest_df,
+        alternative=alternative,
+        oriented=oriented,
+        output_dir=output_dir_path
+    )
+                            
+    if not (result is None): 
+        path_to_file = output_dir.rstrip("/") + \
+            f"/{alternative}_hypergeom.csv"
+        result.to_csv(
+            path_to_file,
+            sep=",",
+            index=None
+        )
+        print(f"File saved at: {path_to_file}") 
+
+def main():
+    import argparse
+
+    parser = argparse.ArgumentParser(
+        prog='dcona',
+        description='DCoNA (Differential Correlation Network Analysis)',
+        epilog='https://github.com/zhiyanov/DCoNA'
+    )
+    parser.add_argument(
+        "tool", choices=['ztest', 'zscore', 'hypergeom'],
+        help="One of DCoNA tools"
+    )
+    parser.add_argument(
+        "config_path",
+        help="Path to JSON config file"
+    )
+
+    args = parser.parse_args()
+    if args.tool=="ztest":
+        ztest_cli(args.config_path)
+    elif args.tool=="zscore":
+        zscore_cli(args.config_path)
+    else:
+        hypergeom_cli(args.config_path)
+
+if __name__=="__main__":
+    main()
diff --git a/dcona/core/__init__.py b/dcona/core/__init__.py
new file mode 100644
index 0000000..8e918f9
--- /dev/null
+++ b/dcona/core/__init__.py
@@ -0,0 +1,2 @@
+from . import extern
+from . import utils
diff --git a/core/correlations/__init__.py b/dcona/core/correlations/__init__.py
similarity index 100%
rename from core/correlations/__init__.py
rename to dcona/core/correlations/__init__.py
diff --git a/core/correlations/correlation_utils.py b/dcona/core/correlations/correlation_utils.py
similarity index 100%
rename from core/correlations/correlation_utils.py
rename to dcona/core/correlations/correlation_utils.py
diff --git a/core/extern/__init__.py b/dcona/core/extern/__init__.py
similarity index 100%
rename from core/extern/__init__.py
rename to dcona/core/extern/__init__.py
diff --git a/core/extern/pcorrelations.py b/dcona/core/extern/pcorrelations.py
similarity index 100%
rename from core/extern/pcorrelations.py
rename to dcona/core/extern/pcorrelations.py
diff --git a/core/extern/ppipelines.py b/dcona/core/extern/ppipelines.py
similarity index 98%
rename from core/extern/ppipelines.py
rename to dcona/core/extern/ppipelines.py
index 5d1e6e6..8bc68c3 100644
--- a/core/extern/ppipelines.py
+++ b/dcona/core/extern/ppipelines.py
@@ -12,6 +12,7 @@
 from .putils import \
     reorder
 
+
 def get_num_ind(indexes, *args):
     index_hash = {
         ind: num for num, ind in enumerate(indexes)
@@ -32,7 +33,7 @@ def ztest_pipeline(
     source_indexes=None,
     target_indexes=None,
     correlation="spearman",
-    correlation_alternative=None,
+    correlation_alternative=False,
     alternative="two-sided",
     repeats_num=1000,
     process_num=1,
@@ -238,8 +239,4 @@ def score_pipeline(
 
         indexes = np.arange(data.shape[0], dtype="int32")
     
-    return (
-        indexes,
-        scores,
-        pvalues
-    )
+    return indexes, scores, pvalues
diff --git a/core/extern/pscores.py b/dcona/core/extern/pscores.py
similarity index 100%
rename from core/extern/pscores.py
rename to dcona/core/extern/pscores.py
diff --git a/core/extern/putils.py b/dcona/core/extern/putils.py
similarity index 100%
rename from core/extern/putils.py
rename to dcona/core/extern/putils.py
diff --git a/core/utils/__init__.py b/dcona/core/utils/__init__.py
similarity index 51%
rename from core/utils/__init__.py
rename to dcona/core/utils/__init__.py
index 5adb22e..16281fe 100644
--- a/core/utils/__init__.py
+++ b/dcona/core/utils/__init__.py
@@ -1,2 +1,1 @@
 from .utils import *
-from .dump import *
diff --git a/core/utils/utils.py b/dcona/core/utils/utils.py
similarity index 94%
rename from core/utils/utils.py
rename to dcona/core/utils/utils.py
index eac1722..70c7538 100644
--- a/core/utils/utils.py
+++ b/dcona/core/utils/utils.py
@@ -1,10 +1,10 @@
+import pandas as pd
 import numpy as np
 import math
 import tqdm
 import sys
 import os
-
-import core.extern
+import itertools
 
 BOOTSTRAP_REPEATS = 10**3
 
@@ -49,5 +49,4 @@ def bootstrap_sample(
         if (statistic != None):
             yield statistic(*samples)
         else:
-            yield sample
-
+            yield sample 
diff --git a/dcona/lib/__init__.py b/dcona/lib/__init__.py
new file mode 100644
index 0000000..d575914
--- /dev/null
+++ b/dcona/lib/__init__.py
@@ -0,0 +1,3 @@
+from .ztest import ztest
+from .zscore import zscore
+from .hypergeom import hypergeom
diff --git a/core/utils/dump.py b/dcona/lib/dump.py
similarity index 75%
rename from core/utils/dump.py
rename to dcona/lib/dump.py
index 5625065..4c0d4f0 100644
--- a/core/utils/dump.py
+++ b/dcona/lib/dump.py
@@ -2,17 +2,21 @@
 import sys
 import os
 
-import core.extern
+from ..core import extern as cextern
 
-CHUNK_LENGTH = 10**6
+CHUNK_LENGTH = 10**5
 
 
 def check_directory_existence(
     directory_path,
-    message="Output directory does not exist"
+    message_isdir="Output directory does not exist",
+    message_access="Permission denied (output directory)"
 ):
     if os.path.isdir(directory_path) == False:
-        print(message)
+        print(message_isdir)
+        sys.exit()
+    if os.access(directory_path, os.W_OK | os.X_OK) == False:
+        print(message_access)
         sys.exit()
 
 def save_by_chunks(
@@ -22,13 +26,12 @@ def save_by_chunks(
     index_transform=None,
     chunk_length=CHUNK_LENGTH
 ):
-    df_inds = []
-    rows = len(indexes)
     MODE, HEADER = 'w', True
     
     # Splitting rows into chunks
     chunk_number = (len(indexes) + chunk_length - 1) // chunk_length
-    for i in tqdm.tqdm(range(chunk_number), desc="Save progress"):
+
+    for i in tqdm.tqdm(range(chunk_number), desc="Save progress: ", ascii=True):
         start = i * chunk_length
         end = (i + 1) * chunk_length
         if i >= chunk_number - 1:
@@ -42,7 +45,7 @@ def save_by_chunks(
             if not (index_transform is None):
                 ind = index_transform[ind]
 
-            s, t = core.extern.paired_index(ind, len(df_indexes))
+            s, t = cextern.paired_index(ind, len(df_indexes))
             source_indexes.append(df_indexes[s])
             target_indexes.append(df_indexes[t])
         
@@ -62,3 +65,4 @@ def save_by_chunks(
         )
         
         MODE, HEADER = 'a', False
+
diff --git a/dcona/lib/hypergeom.py b/dcona/lib/hypergeom.py
new file mode 100644
index 0000000..089a9b5
--- /dev/null
+++ b/dcona/lib/hypergeom.py
@@ -0,0 +1,92 @@
+import numpy as np
+import pandas as pd
+import scipy.stats
+
+from . import dump
+
+FDR_THRESHOLD = 0.05
+
+
+def hypergeom(
+    ztest_df,
+    alternative="two-sided",
+    oriented=True,
+    output_dir=None
+):
+    if alternative == "less":
+        report_df = ztest_df[
+            (report_df["AdjPvalue"] < FDR_THRESHOLD) &
+            (report_df["Statistic"] < 0)
+        ]
+    elif alternative == "greater":
+        report_df = ztest_df[
+            (report_df["AdjPvalue"] < FDR_THRESHOLD) &
+            (report_df["Statistic"] > 0)
+        ]
+    else:
+        report_df = ztest_df[report_df["AdjPvalue"] < FDR_THRESHOLD]
+    
+    if oriented:
+        report_occurrence = get_occurrence(list(report_df["Source"]))
+        initial_occurrence = get_occurrence(list(data_df["Source"]))
+    else:
+        ll = list(report_df["Source"])
+        ll.extend(list(report_df["Target"]))
+        report_occurrence = get_occurrence(ll)
+
+        ll = list(data_df["Source"])
+        ll.extend(list(data_df["Target"]))
+        initial_occurrence = get_occurrence(ll)
+
+    report_interaction_number = 0
+    for index in report_occurrence:
+        report_interaction_number += report_occurrence[index]
+
+    initial_interaction_number = 0
+    for index in initial_occurrence:
+        initial_interaction_number += initial_occurrence[index]
+
+    output_df = pd.DataFrame()
+    output_df["Molecule"] = [index for index in report_occurrence]
+    output_df["Diff"] = [report_occurrence[molecule] for molecule in output_df["Molecule"]]
+    output_df["Total"] = [initial_occurrence[molecule] for molecule in output_df["Molecule"]]
+    output_df["Proportion"] = output_df["Diff"] / output_df["Total"]
+    output_df["Pvalue"] = 1 - scipy.stats.hypergeom.cdf(
+        output_df["Diff"] - 1,
+        initial_interaction_number,
+        report_interaction_number,
+        output_df["Total"]
+    )
+
+    adjusted_pvalue = np.array(output_df["Pvalue"] * len(output_df["Pvalue"]) / \
+        scipy.stats.rankdata(output_df["Pvalue"]))
+    adjusted_pvalue[adjusted_pvalue > 1] = 1
+    adjusted_pvalue = adjusted_pvalue.flatten()
+    output_df["AdjPvalue"] = adjusted_pvalue
+    
+    if output_dir:
+        dump.check_directory_existence(output_dir)
+        path_to_file = output_dir.rstrip("/") + \
+            f"/{alternative}_hypergeom.csv"
+        output_df.to_csv(
+            path_to_file,
+            sep=",",
+            index=None
+        )
+
+        print(f"File saved at: {path_to_file}")
+        return None
+    else:
+        return output_df
+
+    return output_df
+
+def get_occurrence(array):
+    occurrence = {}
+    for elem in array:
+        if elem in occurrence:
+            occurrence[elem] += 1
+        else:
+            occurrence[elem] = 1
+    
+    return occurrence
diff --git a/dcona/lib/utils.py b/dcona/lib/utils.py
new file mode 100644
index 0000000..e3e2701
--- /dev/null
+++ b/dcona/lib/utils.py
@@ -0,0 +1,60 @@
+import itertools
+import pandas as pd
+import numpy as np
+import json
+
+
+def generate_pairs(objects):    
+    if isinstance(objects, (list, set, pd.core.series.Series, np.ndarray)):
+        pairs = list(itertools.combinations(objects, 2))
+    elif isinstance(objects, str):
+        pairs = objects.replace(",", " ").split(" ")
+        pairs = list(filter(lambda a: a != "", pairs))
+        pairs = list(itertools.combinations(pairs, 2))
+    else:
+        try:
+            pairs = list(itertools.combinations(objects, 2))
+        except TypeError:
+            print("Use appropriate data type: list, np.array, string, etc.")
+            raise
+        except:
+            raise
+    
+    pairs_df = pd.DataFrame(pairs, columns=["Source", "Target"])
+    return pairs_df
+
+# Reads config from a json file
+def read_config(
+    path_to_config
+):
+    config = json.load(open(path_to_config, "r"))
+    data_path = config["data_path"]
+    description_path = config["description_path"]
+    reference_group = config["reference_group"]
+    experimental_group = config["experimental_group"]
+    
+    correlation = config["correlation"]
+    alternative = config["alternative"]
+    score = config["score"]
+    
+    output_dir_path = config["output_dir_path"]
+    repeats_number = config["repeats_number"]
+    process_number = config["process_number"]
+        
+    if ("interaction_path" in config) and (config["interaction_path"] != ""):
+        interaction_path = config["interaction_path"]
+    else:
+        interaction_path = none
+        
+    if ("score" in config) and (config["score"] != ""):
+        score = config["score"]
+    else:
+        score = None
+   
+    return data_path, description_path, \
+        reference_group, experimental_group, \
+        correlation, alternative, score, \
+        interaction_path, \
+        repeats_number, \
+        output_dir_path, \
+        process_number
diff --git a/dcona/lib/zscore.py b/dcona/lib/zscore.py
new file mode 100644
index 0000000..bff8296
--- /dev/null
+++ b/dcona/lib/zscore.py
@@ -0,0 +1,131 @@
+import numpy as np
+import pandas as pd
+import scipy.stats
+
+from ..core import extern
+from . import utils
+from . import dump
+
+
+def zscore(
+    data_df, description_df,
+    reference_group, experimental_group,
+    correlation="spearman", score="mean", alternative="two-sided",
+    interaction=None,
+    repeats_number=None,
+    output_dir=None,
+    process_number=None 
+):
+    if process_number is None:
+        process_number = 1
+
+    # if gene names are in dataframe column, relocate them in df.index
+    # if not pd.api.types.is_number(data_df.iloc[0, 0]):
+    #     data_df = data_df.copy()
+    #     data_df.set_index(data_df.columns[0], inplace=True)
+        
+    if interaction is not None:
+        if isinstance(interaction, pd.core.frame.DataFrame):
+            interaction_df = interaction
+        else:
+            interaction_df = utils.generate_pairs(interaction)
+        
+        if repeats_number is None:
+            repeats_number = int(len(interaction_df) / 0.05)
+    else:
+        interaction_df = None
+        
+        if repeats_number is None:
+            repeats_number = int(len(data_df) / 0.05)
+
+    data_df, sources, scores, \
+    pvalues, adjusted_pvalue = \
+    _zscore(
+        data_df, description_df, interaction_df, \
+        reference_group, experimental_group, \
+        correlation, score, alternative, \
+        repeats_number, process_number
+    )
+
+    output_df = pd.DataFrame(data={
+        "Source": data_df.index.to_numpy()[sources],
+        "Score": scores,
+        "Pvalue": pvalues, 
+        "AdjPvalue": adjusted_pvalue, 
+    })
+    output_df = output_df.sort_values(["AdjPvalue", "Pvalue"])
+                            
+    if output_dir:
+        dump.check_directory_existence(output_dir)
+        path_to_file = output_dir.rstrip("/") + \
+            f"/{correlation}_{score}_{alternative}_zscore.csv"
+        output_df.to_csv(
+            path_to_file,
+            sep=",",
+            index=None
+        )
+
+        print(f"File saved at: {path_to_file}")
+        return None
+    else:
+        return output_df
+
+def _zscore(
+    data_df, description_df, interaction_df,
+    reference_group, experimental_group,
+    correlation, score, alternative,
+    repeats_number, process_number
+):
+    if (correlation != "spearman"):
+        correlation = "pearson"
+
+    if interaction_df is not None:
+        data_molecules = set(data_df.index.to_list())
+        interaction_df = interaction_df[interaction_df["Source"].isin(data_molecules)]
+        interaction_df = interaction_df[interaction_df["Target"].isin(data_molecules)]
+
+        source_indexes = interaction_df["Source"]
+        target_indexes = interaction_df["Target"]
+
+        data_molecules = data_molecules.intersection(
+            set(source_indexes) | set(target_indexes)
+        )
+        data_df = data_df.loc[data_df.index.isin(data_molecules)]
+    else:
+        interaction_df = None
+        source_indexes = None
+        target_indexes = None 
+
+    reference_indexes = description_df.loc[
+        description_df["Group"] == reference_group,
+        "Sample"
+    ].to_list()
+    experimental_indexes = description_df.loc[
+        description_df["Group"] == experimental_group,
+        "Sample"
+    ].to_list()
+
+    print("Computation of z-scores")
+    sources, scores, pvalues = \
+    extern.score_pipeline(
+        data_df,
+        reference_indexes,
+        experimental_indexes,
+        source_indexes,
+        target_indexes,
+        correlation=correlation,
+        score=score,
+        alternative=alternative,
+        repeats_num=repeats_number,
+        process_num=process_number
+    )
+
+    print("Computation of adjusted p-value")
+    adjusted_pvalue = pvalues * len(pvalues) / \
+        scipy.stats.rankdata(pvalues)
+    adjusted_pvalue[adjusted_pvalue > 1] = 1
+    adjusted_pvalue = adjusted_pvalue.flatten()
+    
+    return data_df, sources, scores, \
+        pvalues, adjusted_pvalue
+
diff --git a/native/CMakeLists.txt b/dcona/native/CMakeLists.txt
similarity index 100%
rename from native/CMakeLists.txt
rename to dcona/native/CMakeLists.txt
diff --git a/native/src/correlations/correlations.cpp b/dcona/native/src/correlations/correlations.cpp
similarity index 100%
rename from native/src/correlations/correlations.cpp
rename to dcona/native/src/correlations/correlations.cpp
diff --git a/native/src/correlations/correlations.h b/dcona/native/src/correlations/correlations.h
similarity index 100%
rename from native/src/correlations/correlations.h
rename to dcona/native/src/correlations/correlations.h
diff --git a/native/src/lib/correlations.cpp b/dcona/native/src/lib/correlations.cpp
similarity index 100%
rename from native/src/lib/correlations.cpp
rename to dcona/native/src/lib/correlations.cpp
diff --git a/native/src/lib/pipelines.cpp b/dcona/native/src/lib/pipelines.cpp
similarity index 97%
rename from native/src/lib/pipelines.cpp
rename to dcona/native/src/lib/pipelines.cpp
index e0cdfa7..92fb93e 100644
--- a/native/src/lib/pipelines.cpp
+++ b/dcona/native/src/lib/pipelines.cpp
@@ -15,6 +15,8 @@
 
 #include "../pipelines/pipelines.h"
 #include "../utils/utils.h"
+// #include "../tqdm/tqdm.h"
+#include "../tqdm/progressbar.h"
 
 namespace py = pybind11;
 
@@ -107,6 +109,9 @@ std::tuple<
     // Bootstrap pvalue computations    
     float *dpr, *rcp, *ecp, *sp, *pp;
     int *rip, *eip;
+    
+    std::cout << "Permutation progress: ";
+    progressbar bar(repeats_number + 1);
     for (int r = 0; r < repeats_number + 1; ++r) {
         if (r == 0) {
             rcp = ref_corrs_ptr;
@@ -117,8 +122,6 @@ std::tuple<
             rip = ref_ind_ptr;
             eip = exp_ind_ptr;
         } else {
-            std::cout << "Bootstrap iteration: " << r << "/"
-                << repeats_number << "\n";
             std::shuffle(indexes.begin(), indexes.end(), random_gen);
             for (int i = 0; i < ref_ind_size; ++i) {
                 boot_ref_ind_ptr[i] = indexes[i];
@@ -221,7 +224,10 @@ std::tuple<
                 boot_pvalue_ptr[i] = 0;
             }
         }
+
+        bar.update();
     }
+    std::cout << "\n";
     
     if (repeats_number > 0) {
         for (int i = 0; i < index_size; ++i) {
@@ -369,6 +375,9 @@ std::tuple<
     float *dpr, *rcp, *ecp, *sp, *pp;
     int *rip, *eip;
     float *scp;
+    
+    std::cout << "Permutation progress: ";
+    progressbar bar(repeats_number + 1);
     for (int r = 0; r < repeats_number + 1; ++r) {
         if (r == 0) {
             rcp = boot_ref_corrs_ptr;
@@ -381,8 +390,6 @@ std::tuple<
 
             scp = scores_ptr;
         } else {
-            std::cout << "Bootstrap iteration: " << r << "/"
-                << repeats_number << "\n";
             std::shuffle(indexes.begin(), indexes.end(), random_gen);
             for (int i = 0; i < ref_ind_size; ++i) {
                 boot_ref_ind_ptr[i] = indexes[i];
@@ -492,7 +499,10 @@ std::tuple<
                 pvalues_ptr[i] = 0;
             }
         }
+
+        bar.update();
     }
+    std::cout << "\n";
 
     if (repeats_number > 0) {
         for (int i = 0; i < sources_size; ++i) {
@@ -598,6 +608,9 @@ std::tuple<
     // Bootstrap pvalue computations    
     float *dpr, *rcp, *ecp, *sp, *pp;
     int *rip, *eip;
+    
+    std::cout << "Permutation progress: ";
+    progressbar bar(repeats_number + 1);
     for (int r = 0; r < repeats_number + 1; ++r) {
         if (r == 0) {
             rcp = ref_corrs_ptr;
@@ -608,8 +621,6 @@ std::tuple<
             rip = ref_ind_ptr;
             eip = exp_ind_ptr;
         } else {
-            std::cout << "Bootstrap iteration: " << r << "/"
-                << repeats_number << "\n";
             std::shuffle(indexes.begin(), indexes.end(), random_gen);
             for (int i = 0; i < ref_ind_size; ++i) {
                 boot_ref_ind_ptr[i] = indexes[i];
@@ -621,7 +632,7 @@ std::tuple<
             rcp = boot_ref_corrs_ptr;
             ecp = boot_exp_corrs_ptr;
             sp  = boot_stat_ptr;
-            pp  = (float *) nullptr;
+            pp  = nullptr;
             
             rip = boot_ref_ind_ptr;
             eip = boot_exp_ind_ptr;
@@ -712,7 +723,10 @@ std::tuple<
                 boot_pvalue_ptr[i] = 0;
             }
         }
+
+        bar.update();
     }
+    std::cout << "\n";
     
     if (repeats_number > 0) {
         for (int i = 0; i < pairs_num; ++i) {
@@ -821,6 +835,9 @@ std::tuple<
     float *dpr, *rcp, *ecp, *sp, *pp;
     int *rip, *eip;
     float *scp;
+    
+    std::cout << "Permutation progress: ";
+    progressbar bar(repeats_number + 1);
     for (int r = 0; r < repeats_number + 1; ++r) {
         if (r == 0) {
             rcp = boot_ref_corrs_ptr;
@@ -833,8 +850,6 @@ std::tuple<
 
             scp = scores_ptr;
         } else {
-            std::cout << "Bootstrap iteration: " << r << "/"
-                << repeats_number << "\n";
             std::shuffle(indexes.begin(), indexes.end(), random_gen);
             for (int i = 0; i < ref_ind_size; ++i) {
                 boot_ref_ind_ptr[i] = indexes[i];
@@ -966,7 +981,10 @@ std::tuple<
                 pvalues_ptr[i] = 0;
             }
         }
+
+        bar.update();
     }
+    std::cout << "\n";
 
     if (repeats_number > 0) {
         for (int i = 0; i < sources_size; ++i) {
diff --git a/native/src/lib/scores.cpp b/dcona/native/src/lib/scores.cpp
similarity index 100%
rename from native/src/lib/scores.cpp
rename to dcona/native/src/lib/scores.cpp
diff --git a/native/src/lib/utils.cpp b/dcona/native/src/lib/utils.cpp
similarity index 100%
rename from native/src/lib/utils.cpp
rename to dcona/native/src/lib/utils.cpp
diff --git a/native/src/pipelines/pipelines.cpp b/dcona/native/src/pipelines/pipelines.cpp
similarity index 100%
rename from native/src/pipelines/pipelines.cpp
rename to dcona/native/src/pipelines/pipelines.cpp
diff --git a/native/src/pipelines/pipelines.h b/dcona/native/src/pipelines/pipelines.h
similarity index 100%
rename from native/src/pipelines/pipelines.h
rename to dcona/native/src/pipelines/pipelines.h
diff --git a/native/src/scores/scores.cpp b/dcona/native/src/scores/scores.cpp
similarity index 100%
rename from native/src/scores/scores.cpp
rename to dcona/native/src/scores/scores.cpp
diff --git a/native/src/scores/scores.h b/dcona/native/src/scores/scores.h
similarity index 100%
rename from native/src/scores/scores.h
rename to dcona/native/src/scores/scores.h
diff --git a/dcona/native/src/tqdm/progressbar.h b/dcona/native/src/tqdm/progressbar.h
new file mode 100644
index 0000000..829ff28
--- /dev/null
+++ b/dcona/native/src/tqdm/progressbar.h
@@ -0,0 +1,184 @@
+// The MIT License (MIT)
+//
+// Copyright (c) 2019 Luigi Pertoldi
+//
+// Permission is hereby granted, free of charge, to any person obtaining a copy
+// of this software and associated documentation files (the "Software"), to
+// deal in the Software without restriction, including without limitation the
+// rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
+// sell copies of the Software, and to permit persons to whom the Software is
+// furnished to do so, subject to the following conditions:
+//
+// The above copyright notice and this permission notice shall be included in
+// all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
+// IN THE SOFTWARE.
+//
+// ============================================================================
+//  ___   ___   ___   __    ___   ____  __   __   ___    __    ___
+// | |_) | |_) / / \ / /`_ | |_) | |_  ( (` ( (` | |_)  / /\  | |_)
+// |_|   |_| \ \_\_/ \_\_/ |_| \ |_|__ _)_) _)_) |_|_) /_/--\ |_| \_
+//
+// Very simple progress bar for c++ loops with internal running variable
+//
+// Author: Luigi Pertoldi
+// Created: 3 dic 2016
+//
+// Notes: The bar must be used when there's no other possible source of output
+//        inside the for loop
+//
+
+#ifndef __PROGRESSBAR_HPP
+#define __PROGRESSBAR_HPP
+
+#include <iostream>
+#include <string>
+#include <stdexcept>
+
+class progressbar {
+
+    public:
+      // default destructor
+      ~progressbar()                             = default;
+
+      // delete everything else
+      progressbar           (progressbar const&) = delete;
+      progressbar& operator=(progressbar const&) = delete;
+      progressbar           (progressbar&&)      = delete;
+      progressbar& operator=(progressbar&&)      = delete;
+
+      // default constructor, must call set_niter later
+      inline progressbar();
+      inline progressbar(int n, bool showbar=true);
+
+      // reset bar to use it again
+      inline void reset();
+     // set number of loop iterations
+      inline void set_niter(int iter);
+      // chose your style
+      inline void set_done_char(const std::string& sym) {done_char = sym;}
+      inline void set_todo_char(const std::string& sym) {todo_char = sym;}
+      inline void set_opening_bracket_char(const std::string& sym) {opening_bracket_char = sym;}
+      inline void set_closing_bracket_char(const std::string& sym) {closing_bracket_char = sym;}
+      // to show only the percentage
+      inline void show_bar(bool flag = true) {do_show_bar = flag;}
+      // main function
+      inline void update();
+
+    private:
+      int progress;
+      int n_cycles;
+      int last_perc;
+      bool do_show_bar;
+      bool update_is_called;
+
+      std::string done_char;
+      std::string todo_char;
+      std::string opening_bracket_char;
+      std::string closing_bracket_char;
+};
+
+inline progressbar::progressbar() :
+    progress(0),
+    n_cycles(0),
+    last_perc(0),
+    do_show_bar(true),
+    update_is_called(false),
+    done_char("#"),
+    todo_char(" "),
+    opening_bracket_char("["),
+    closing_bracket_char("]") {}
+
+inline progressbar::progressbar(int n, bool showbar) :
+    progress(0),
+    n_cycles(n),
+    last_perc(0),
+    do_show_bar(showbar),
+    update_is_called(false),
+    done_char("#"),
+    todo_char(" "),
+    opening_bracket_char("["),
+    closing_bracket_char("]") {}
+
+inline void progressbar::reset() {
+    progress = 0,
+    update_is_called = false;
+    last_perc = 0;
+    return;
+}
+
+inline void progressbar::set_niter(int niter) {
+    if (niter <= 0) throw std::invalid_argument(
+        "progressbar::set_niter: number of iterations null or negative");
+    n_cycles = niter;
+    return;
+}
+
+inline void progressbar::update() {
+
+    if (n_cycles == 0) throw std::runtime_error(
+            "progressbar::update: number of cycles not set");
+
+    if (!update_is_called) {
+        if (do_show_bar == true) {
+            std::cout << opening_bracket_char;
+            for (int _ = 0; _ < 50; _++) std::cout << todo_char;
+            std::cout << closing_bracket_char << " 0%";
+        }
+        else std::cout << "0%";
+    }
+    update_is_called = true;
+
+    int perc = 0;
+
+    // compute percentage, if did not change, do nothing and return
+    perc = progress*100./(n_cycles-1);
+    if (perc < last_perc) return;
+
+    // update percentage each unit
+    if (perc == last_perc + 1) {
+        // erase the correct  number of characters
+        if      (perc <= 10)                std::cout << "\b\b"   << perc << '%';
+        else if (perc  > 10 and perc < 100) std::cout << "\b\b\b" << perc << '%';
+        else if (perc == 100)               std::cout << "\b\b\b" << perc << '%';
+    }
+    if (do_show_bar == true) {
+        // update bar every ten units
+        if (perc % 2 == 0) {
+            // erase closing bracket
+            std::cout << std::string(closing_bracket_char.size(), '\b');
+            // erase trailing percentage characters
+            if      (perc  < 10)               std::cout << "\b\b\b";
+            else if (perc >= 10 && perc < 100) std::cout << "\b\b\b\b";
+            else if (perc == 100)              std::cout << "\b\b\b\b\b";
+
+            // erase 'todo_char'
+            for (int j = 0; j < 50-(perc-1)/2; ++j) {
+                std::cout << std::string(todo_char.size(), '\b');
+            }
+
+            // add one additional 'done_char'
+            if (perc == 0) std::cout << todo_char;
+            else           std::cout << done_char;
+
+            // refill with 'todo_char'
+            for (int j = 0; j < 50-(perc-1)/2-1; ++j) std::cout << todo_char;
+
+            // readd trailing percentage characters
+            std::cout << closing_bracket_char << ' ' << perc << '%';
+        }
+    }
+    last_perc = perc;
+    ++progress;
+    std::cout << std::flush;
+
+    return;
+}
+
+#endif
diff --git a/dcona/native/src/tqdm/tqdm.h b/dcona/native/src/tqdm/tqdm.h
new file mode 100644
index 0000000..7455302
--- /dev/null
+++ b/dcona/native/src/tqdm/tqdm.h
@@ -0,0 +1,484 @@
+#pragma once
+
+/**
+Customisable progressbar decorator for iterators.
+Includes a default range iterator printing to stderr.
+
+TODO:
+* handle s(tep) with operator+/-(=)
+* chrono and delay printing
+* class status printer
+* iterator-less: fake range iterable, update() increments value
+
+Usage:
+  # include "tqdm/tqdm.h"
+  for(int i : tqdm::range(4))
+  // same as:
+  //   for (int i : tqdm::tqdm({0, 1, 2, 3}))
+  // or:
+  //   std::vector<int> v{0, 1, 2, 3};
+  //   for (int &i : tqdm::tqdm(v.begin(), v.end())
+    ...
+
+@author Casper dC-L <github.com/casperdcl>
+*/
+
+#include <cassert>      // assert
+#include <cinttypes>    // PRIu64
+#include <cstddef>      // ptrdiff_t, size_t
+#include <cstdint>      // int64_t
+#include <cstdio>       // printf
+#include <iterator>     // iterator
+#include <limits>       // numeric_limits
+#include <stdexcept>    // throw
+#include <string>       // string
+#include <type_traits>  // is_pointer, ...
+#include <utility>      // swap
+#include "utils.h"
+
+#ifndef SIZE_T_MAX
+constexpr size_t SIZE_T_MAX = std::numeric_limits<size_t>::max();
+#endif
+
+namespace tqdm {
+
+struct Params {
+  std::string desc;
+  size_t total = -1;
+  bool leave = true;
+  FILE *f = stderr;
+  int ncols = -1;
+  float mininterval = 0.1f, maxinterval = 10.0f;
+  unsigned miniters = -1;
+  std::string ascii = " 123456789#";
+  bool disable = false;
+  std::string unit = "it";
+  bool unit_scale = false;
+  bool dynamic_ncols = false;
+  float smoothing = 0.3f;
+  std::string bar_format;
+  size_t initial = 0;
+  int position = -1;
+  bool gui = false;
+};
+
+template <typename _Iterator>
+class Tqdm : public MyIteratorWrapper<_Iterator> {
+private:
+  using TQDM_IT = MyIteratorWrapper<_Iterator>;
+  _Iterator e;  // end
+  Params self;  // ha, ha
+
+public:
+  /**
+   containter-like methods
+   */
+  // actually current value
+  // virtual _Iterator begin() { return this->get(); }
+  Tqdm &begin() { return *this; }
+  const Tqdm &begin() const { return *this; }
+  // virtual _Iterator end() { return e; }
+  Tqdm end() const { return Tqdm(e, e); }
+
+  explicit operator _Iterator() { return this->get(); }
+
+  /** constructors
+   */
+  explicit Tqdm(_Iterator begin, _Iterator end)
+      : TQDM_IT(begin), e(end), self() {
+    self.total = size_t(end - begin);
+  }
+
+  explicit Tqdm(_Iterator begin, size_t total)
+      : TQDM_IT(begin), e(begin + total), self() {
+    self.total = total;
+  }
+
+  // Tqdm(const Tqdm& other)
+  //     : TQDM_IT(other.get()),
+  //       e(other.end().get()),
+  //       self(other.self),
+  // {
+  //   // std::memcpy(this, &other, sizeof(Tqdm));
+  // }
+
+  template <typename _Container,
+            typename = typename std::enable_if<
+                !std::is_same<_Container, Tqdm>::value>::type>
+  Tqdm(_Container &v) : TQDM_IT(std::begin(v)), e(std::end(v)), self() {
+    self.total = e - this->get();
+  }
+
+  explicit operator bool() const { return this->get() != e; }
+
+  /** TODO: magic methods */
+  virtual void _incr() const override {
+    if (this->get() == e)
+      throw std::out_of_range(
+          "exhausted");  // TODO: don't throw, just double total
+
+    TQDM_IT::_incr();
+    if (this->get() == e) {
+      printf("\nfinished: %" PRIu64 "/%" PRIu64 "\n",
+        static_cast<std::uint64_t>(self.total),
+        static_cast<std::uint64_t>(self.total));
+    } else
+      printf("\r%" PRIi64 " left", (int64_t)(e - this->get()));
+  }
+  virtual void _incr() override { ((Tqdm const &)*this)._incr(); }
+};
+
+template <typename _Iterator, typename _Tqdm = Tqdm<_Iterator>>
+_Tqdm tqdm(_Iterator begin, _Iterator end) {
+  return _Tqdm(begin, end);
+}
+
+template <typename _Iterator, typename _Tqdm = Tqdm<_Iterator>>
+_Tqdm tqdm(_Iterator begin, size_t total) {
+  return _Tqdm(begin, total);
+}
+
+template <typename _Container,
+          typename _Tqdm = Tqdm<typename _Container::iterator>>
+_Tqdm tqdm(_Container &v) {
+  return _Tqdm(v);
+}
+
+template <size_t N, typename T, typename _Tqdm = Tqdm<T *>>
+_Tqdm tqdm(T (&tab)[N]) {
+  return _Tqdm(tab, N);
+}
+
+template <typename SizeType = int>
+using RangeTqdm = Tqdm<RangeIterator<SizeType>>;
+template <typename SizeType> RangeTqdm<SizeType> range(SizeType n) {
+  return RangeTqdm<SizeType>(RangeIterator<SizeType>(n),
+                             RangeIterator<SizeType>(n));
+}
+template <typename SizeType>
+RangeTqdm<SizeType> range(SizeType start, SizeType end) {
+  return RangeTqdm<SizeType>(RangeIterator<SizeType>(start, end),
+                             RangeIterator<SizeType>(start, end));
+}
+template <typename SizeType>
+RangeTqdm<SizeType> range(SizeType start, SizeType end, SizeType step) {
+  return RangeTqdm<SizeType>(RangeIterator<SizeType>(start, end, step),
+                             RangeIterator<SizeType>(start, end, step));
+}
+
+}  // tqdm
+
+/** Things to port:
+
+class TqdmTypeError(TypeError):
+class TqdmKeyError(KeyError):
+class TqdmDeprecationWarning(Exception):
+ASCII_FMT = " 123456789#"
+UTF_FMT = u" " + u''.join(map(_unich, range(0x258F, 0x2587, -1)))
+class tqdm(object):
+    @staticmethod
+    def format_sizeof(num, suffix=''):
+    @staticmethod
+    def format_interval(t):
+    @staticmethod
+    def status_printer(file):
+    @staticmethod
+    def format_meter(n, total, elapsed, ncols=None, prefix='',
+                     ascii=False, unit='it', unit_scale=False, rate=None,
+                     bar_format=None):
+    def __new__(cls, *args, **kwargs):
+    @classmethod
+    def _get_free_pos(cls, instance=None):
+    @classmethod
+    def _decr_instances(cls, instance):
+    @classmethod
+    def write(cls, s, file=sys.stdout, end="\n"):
+    @classmethod
+    def pandas(tclass, *targs, **tkwargs):
+    def __init__(self, iterable=None, desc=None, total=None, leave=True,
+                 file=sys.stderr, ncols=None, mininterval=0.1,
+                 maxinterval=10.0, miniters=None, ascii=None, disable=False,
+                 unit='it', unit_scale=False, dynamic_ncols=False,
+                 smoothing=0.3, bar_format=None, initial=0, position=None,
+                 gui=False, **kwargs):
+        """
+        Parameters
+        ----------
+        iterable  : iterable, optional
+            Iterable to decorate with a progressbar.
+            Leave blank to manually manage the updates.
+        desc  : str, optional
+            Prefix for the progressbar.
+        total  : int, optional
+            The number of expected iterations. If unspecified,
+            len(iterable) is used if possible. As a last resort, only basic
+            progress statistics are displayed (no ETA, no progressbar).
+            If `gui` is True and this parameter needs subsequent updating,
+            specify an initial arbitrary large positive integer,
+            e.g. int(9e9).
+        leave  : bool, optional
+            If [default: True], keeps all traces of the progressbar
+            upon termination of iteration.
+        file  : `io.TextIOWrapper` or `io.StringIO`, optional
+            Specifies where to output the progress messages
+            [default: sys.stderr]. Uses `file.write(str)` and `file.flush()`
+            methods.
+        ncols  : int, optional
+            The width of the entire output message. If specified,
+            dynamically resizes the progressbar to stay within this bound.
+            If unspecified, attempts to use environment width. The
+            fallback is a meter width of 10 and no limit for the counter and
+            statistics. If 0, will not print any meter (only stats).
+        mininterval  : float, optional
+            Minimum progress update interval, in seconds [default: 0.1].
+        maxinterval  : float, optional
+            Maximum progress update interval, in seconds [default: 10.0].
+        miniters  : int, optional
+            Minimum progress update interval, in iterations.
+            If specified, will set `mininterval` to 0.
+        ascii  : bool or str, optional
+            If unspecified or False, use unicode (smooth blocks) to fill
+            the meter.
+            The fallback is to use ASCII characters " 123456789#".
+        disable  : bool, optional
+            Whether to disable the entire progressbar wrapper
+            [default: False].
+        unit  : str, optional
+            String that will be used to define the unit of each iteration
+            [default: it].
+        unit_scale  : bool, optional
+            If set, the number of iterations will be reduced/scaled
+            automatically and a metric prefix following the
+            International System of Units standard will be added
+            (kilo, mega, etc.) [default: False].
+        dynamic_ncols  : bool, optional
+            If set, constantly alters `ncols` to the environment (allowing
+            for window resizes) [default: False].
+        smoothing  : float, optional
+            Exponential moving average smoothing factor for speed estimates
+            (ignored in GUI mode). Ranges from 0 (average speed) to 1
+            (current/instantaneous speed) [default: 0.3].
+        bar_format  : str, optional
+            Specify a custom bar string formatting. May impact performance.
+            If unspecified, will use '{l_bar}{bar}{r_bar}', where l_bar is
+            '{desc}{percentage:3.0f}%|' and r_bar is
+            '| {n_fmt}/{total_fmt} [{elapsed_str}<{remaining_str},
+{rate_fmt}]'
+            Possible vars: bar, n, n_fmt, total, total_fmt, percentage,
+            rate, rate_fmt, elapsed, remaining, l_bar, r_bar, desc.
+        initial  : int, optional
+            The initial counter value. Useful when restarting a progress
+            bar [default: 0].
+        position  : int, optional
+            Specify the line offset to print this bar (starting from 0)
+            Automatic if unspecified.
+            Useful to manage multiple bars at once (eg, from threads).
+        gui  : bool, optional
+            WARNING: internal parameter - do not use.
+            Use tqdm_gui(...) instead. If set, will attempt to use
+            matplotlib animations for a graphical output [default: False].
+
+        Returns
+        -------
+        out  : decorated iterator.
+        """
+        if disable:
+            self.iterable = iterable
+            self.disable = disable
+            self.pos = self._get_free_pos(self)
+            self._instances.remove(self)
+            return
+
+        if kwargs:
+            self.disable = True
+            self.pos = self._get_free_pos(self)
+            self._instances.remove(self)
+            raise (TqdmDeprecationWarning("""\
+`nested` is deprecated and automated. Use position instead for manual control.
+""", fp_write=getattr(file, 'write', sys.stderr.write))
+                if "nested" in kwargs else
+                TqdmKeyError("Unknown argument(s): " + str(kwargs)))
+
+        # Preprocess the arguments
+        if total is None and iterable is not None:
+            try:
+                total = len(iterable)
+            except (TypeError, AttributeError):
+                total = None
+
+        if ((ncols is None) and (file in (sys.stderr, sys.stdout))) or \
+                dynamic_ncols:  # pragma: no cover
+            if dynamic_ncols:
+                dynamic_ncols = _environ_cols_wrapper()
+                ncols = dynamic_ncols(file)
+            else:
+                ncols = _environ_cols_wrapper()(file)
+
+        if miniters is None:
+            miniters = 0
+            dynamic_miniters = True
+        else:
+            dynamic_miniters = False
+
+        if mininterval is None:
+            mininterval = 0
+
+        if maxinterval is None:
+            maxinterval = 0
+
+        if ascii is None:
+            ascii = not _supports_unicode(file)
+
+        if bar_format and not ((ascii is True) or _is_ascii(ascii)):
+            # Convert bar format into unicode since terminal uses unicode
+            bar_format = _unicode(bar_format)
+
+        if smoothing is None:
+            smoothing = 0
+
+        # Store the arguments
+        self.iterable = iterable
+        self.desc = desc + ': ' if desc else ''
+        self.total = total
+        self.leave = leave
+        self.fp = file
+        self.ncols = ncols
+        self.mininterval = mininterval
+        self.maxinterval = maxinterval
+        self.miniters = miniters
+        self.dynamic_miniters = dynamic_miniters
+        self.ascii = ascii
+        self.disable = disable
+        self.unit = unit
+        self.unit_scale = unit_scale
+        self.gui = gui
+        self.dynamic_ncols = dynamic_ncols
+        self.smoothing = smoothing
+        self.avg_time = None
+        self._time = time
+        self.bar_format = bar_format
+
+        # Init the iterations counters
+        self.last_print_n = initial
+        self.n = initial
+
+        # if nested, at initial sp() call we replace '\r' by '\n' to
+        # not overwrite the outer progress bar
+        self.pos = self._get_free_pos(self) if position is None else position
+
+        if not gui:
+            # Initialize the screen printer
+            self.sp = self.status_printer(self.fp)
+            if self.pos:
+                self.moveto(self.pos)
+            self.sp(self.format_meter(self.n, total, 0,
+                    (dynamic_ncols(file) if dynamic_ncols else ncols),
+                    self.desc, ascii, unit, unit_scale, None, bar_format))
+            if self.pos:
+                self.moveto(-self.pos)
+
+        # Init the time counter
+        self.start_t = self.last_print_t = self._time()
+
+    def __len__(self):
+    def __enter__(self):
+    def __exit__(self, *exc):
+    def __del__(self):
+    def __repr__(self):
+    def __lt__(self, other):
+    def __le__(self, other):
+    def __eq__(self, other):
+    def __ne__(self, other):
+    def __gt__(self, other):
+    def __ge__(self, other):
+    def __hash__(self):
+    def __iter__(self):
+    def update(self, n=1):
+        """
+        Manually update the progress bar, useful for streams
+        such as reading files.
+        E.g.:
+        >>> t = tqdm(total=filesize) # Initialise
+        >>> for current_buffer in stream:
+        ...    ...
+        ...    t.update(len(current_buffer))
+        >>> t.close()
+        The last line is highly recommended, but possibly not necessary if
+        `t.update()` will be called in such a way that `filesize` will be
+        exactly reached and printed.
+
+        Parameters
+        ----------
+        n  : int
+            Increment to add to the internal counter of iterations
+            [default: 1].
+        """
+        if self.disable:
+            return
+
+        if n < 0:
+            raise ValueError("n ({0}) cannot be negative".format(n))
+        self.n += n
+
+        delta_it = self.n - self.last_print_n  # should be n?
+        if delta_it >= self.miniters:
+            # We check the counter first, to reduce the overhead of time()
+            cur_t = self._time()
+            delta_t = cur_t - self.last_print_t
+            if delta_t >= self.mininterval:
+                elapsed = cur_t - self.start_t
+                # EMA (not just overall average)
+                if self.smoothing and delta_t:
+                    self.avg_time = delta_t / delta_it \
+                        if self.avg_time is None \
+                        else self.smoothing * delta_t / delta_it + \
+                        (1 - self.smoothing) * self.avg_time
+
+                if not hasattr(self, "sp"):
+                    raise TqdmDeprecationWarning("""\
+Please use `tqdm_gui(...)` instead of `tqdm(..., gui=True)`
+""", fp_write=getattr(self.fp, 'write', sys.stderr.write))
+
+                if self.pos:
+                    self.moveto(self.pos)
+
+                # Print bar's update
+                self.sp(self.format_meter(
+                    self.n, self.total, elapsed,
+                    (self.dynamic_ncols(self.fp) if self.dynamic_ncols
+                     else self.ncols),
+                    self.desc, self.ascii, self.unit, self.unit_scale,
+                    1 / self.avg_time if self.avg_time else None,
+                    self.bar_format))
+
+                if self.pos:
+                    self.moveto(-self.pos)
+
+                # If no `miniters` was specified, adjust automatically to the
+                # maximum iteration rate seen so far.
+                # e.g.: After running `tqdm.update(5)`, subsequent
+                # calls to `tqdm.update()` will only cause an update after
+                # at least 5 more iterations.
+                if self.dynamic_miniters:
+                    if self.maxinterval and delta_t > self.maxinterval:
+                        self.miniters = self.miniters * self.maxinterval \
+                            / delta_t
+                    elif self.mininterval and delta_t:
+                        self.miniters = self.smoothing * delta_it \
+                            * self.mininterval / delta_t + \
+                            (1 - self.smoothing) * self.miniters
+                    else:
+                        self.miniters = self.smoothing * delta_it + \
+                            (1 - self.smoothing) * self.miniters
+
+                # Store old values for next call
+                self.last_print_n = self.n
+                self.last_print_t = cur_t
+
+    def close(self):
+    def unpause(self):
+    def set_description(self, desc=None):
+    def moveto(self, n):
+    def clear(self, nomove=False):
+    def refresh(self):
+*/
diff --git a/dcona/native/src/tqdm/utils.h b/dcona/native/src/tqdm/utils.h
new file mode 100644
index 0000000..34f7820
--- /dev/null
+++ b/dcona/native/src/tqdm/utils.h
@@ -0,0 +1,397 @@
+#pragma once
+
+#ifndef CUR_OS
+
+#if !defined(IS_WIN) && (defined(_WIN32) || defined(_WIN64))
+#define IS_WIN
+#elif !defined IS_NIX &&                                                     \
+    (defined(unix) || defined(__unix) || defined(__unix__) ||                \
+     defined(__linux__) || defined(__FreeBSD__) || defined(__APPLE__) ||     \
+     defined(__MACH__))
+#define IS_NIX
+#endif
+
+#endif  // CUR_OS
+
+#include <cassert>      // assert
+#include <cstddef>      // ptrdiff_t, size_t
+#include <unistd.h>     // STDERR_FILENO
+#include <iterator>     // iterator
+#include <type_traits>  // is_pointer, ...
+#include <atomic>       // atomic
+#include <cstring>      // strlen
+#include <cerrno>       // EAGAIN
+#include <poll.h>       // poll
+
+/** TODO: port from python
+ * colorama win
+ * weakset
+ * _is_utf(encoding)
+ * _is_ascii(s)
+ * _supports_unicode(file)
+ * _environ_cols()
+ * _sh(const char *cmd[], ...)
+ */
+
+#ifndef constexpr
+#define constexpr static const
+#endif
+
+#ifndef explicit
+#define explicit
+#endif
+
+#ifndef noexcept
+#define noexcept
+#endif
+
+namespace tqdm {
+
+template <typename _Iterator>
+/**
+Wrapper for pointers and std containter iterators.
+@author Casper da Costa-Luis
+*/
+class MyIteratorWrapper
+    : public std::iterator<
+          std::forward_iterator_tag,
+          typename std::iterator_traits<_Iterator>::value_type> {
+  mutable _Iterator p;  // TODO: remove this mutable
+
+public:
+  // already done by std::iterator
+  typedef typename std::iterator_traits<_Iterator>::value_type value_type;
+
+  explicit MyIteratorWrapper(_Iterator x) : p(x) {}
+  // default construct gives end
+  MyIteratorWrapper() : p(nullptr) {}
+  explicit MyIteratorWrapper(const MyIteratorWrapper &mit) : p(mit.p) {}
+
+  // override this in Tqdm class
+  virtual void _incr() { ++p; }
+  // override this in Tqdm class
+  virtual void _incr() const { ++p; }
+
+  MyIteratorWrapper &operator++() {
+    // assert(this->bool() && "Out-of-bounds iterator increment");
+    _incr();
+    return *this;
+  }
+  const MyIteratorWrapper &operator++() const {
+    _incr();
+    return *this;
+  }
+  MyIteratorWrapper operator++(int)const {
+    MyIteratorWrapper tmp(*this);
+    _incr();
+    return tmp;
+  }
+  template <class Other>
+  // two-way comparison: v.begin() == v.cbegin() and vice versa
+  bool operator==(const MyIteratorWrapper<Other> &rhs) const {
+    return p == rhs.p;
+  }
+  template <class Other>
+  bool operator!=(const MyIteratorWrapper<Other> &rhs) const {
+    return p != rhs.p;
+  }
+  template <class Other>
+  size_t operator-(const MyIteratorWrapper<Other> &rhs) {
+    return p - rhs.p;
+  }
+  // template <typename = typename std::enable_if<
+  //               !std::is_const<value_type>::value>::type>
+  value_type &operator*() {
+    // assert(this->bool() && "Invalid iterator dereference!");
+    return *p;
+  }
+  const value_type &operator*() const {
+    // assert(this->bool() && "Invalid iterator dereference!");
+    return *p;
+  }
+  // template <typename = typename std::enable_if<
+  //               !std::is_const<value_type>::value>::type>
+  value_type &operator->() {
+    // assert(this->bool() && "Invalid iterator dereference!");
+    return *p;
+  }
+  const value_type &operator->() const {
+    // assert(this->bool() && "Invalid iterator dereference!");
+    return *p;
+  }
+  // @return the underlying iterator
+  _Iterator &get() { return p; }
+  const _Iterator &get() const { return p; }
+  // TODO: const _Iterator &get() const { return p; }, etc ...
+  //
+  void swap(MyIteratorWrapper &other) noexcept { std::swap(p, other.p); }
+
+  // One way conversion: iterator -> const_iterator
+  // template <typename =
+  //               typename
+  //               std::enable_if<!std::is_const<_Iterator>::value>::type>
+  // operator MyIteratorWrapper<typename std::add_const<_Iterator>::type>()
+  // const {
+  //   return MyIteratorWrapper<const _Iterator>(p);
+  // }
+  template <typename = typename std::is_pointer<_Iterator>>
+  explicit operator bool() const {
+    return p != nullptr;
+  }
+};
+
+template <typename _Iterator,
+          typename _MyIteratorWrapper = MyIteratorWrapper<_Iterator>>
+_MyIteratorWrapper myIteratorWrapper(_Iterator x) {
+  return _MyIteratorWrapper(x);
+}
+
+template <typename IntType = int>
+class RangeIterator
+    : public std::iterator<std::forward_iterator_tag, IntType> {
+private:
+  mutable IntType current;
+  IntType total;
+  IntType step;
+
+public:
+  RangeIterator(IntType total) : current(0), total(total), step(1) {}
+  RangeIterator(IntType start, IntType total)
+      : current(start), total(total), step(1) {}
+  RangeIterator(IntType start, IntType total, IntType step)
+      : current(start), total(total), step(step) {}
+  IntType &operator*() { return current; }
+  const IntType &operator*() const { return current; }
+  RangeIterator &operator++() {
+    current += step;
+    return *this;
+  }
+  const RangeIterator &operator++() const {
+    current += step;
+    return *this;
+  }
+  RangeIterator operator++(int)const {
+    RangeIterator tmp(*this);
+    operator++();
+    return tmp;
+  }
+  explicit operator bool() const { return current < total; }
+  size_t size_remaining() const { return (total - current) / step; }
+
+  /** here be dragons */
+
+  // only use as (it != end), not as (end != it)
+  bool operator!=(const RangeIterator &) const { return current < total; }
+  bool operator==(const RangeIterator &) const { return current >= total; }
+  IntType operator-(const RangeIterator &it) const {
+    // it's used in `end - begin`, but `end` is only a sentinel
+    // so let's use `begin `to be consistent
+    return it.size_remaining();
+  }
+};
+
+const char *_term_move_up() {
+  return
+#if defined(IS_WIN) && !defined(colorama)
+      ""
+#else
+      "\x1b[A"
+#endif
+      ;
+}
+
+static void wait_for_write(int fd) {
+  struct pollfd pfd;
+  pfd.fd = fd;
+  pfd.events = POLLOUT;
+  (void)::poll(&pfd, 1, -1);
+}
+
+// Write a buffer fully or not at all.
+// If false is returned, caller may check errno to see if it's EAGAIN
+// or a real error.
+bool write_harder(int fd, const char *buf, size_t len) {
+  bool did_anything = false;
+
+  while (len) {
+    ssize_t res = ::write(fd, buf, len);
+    if (res == -1) {
+      if (errno == EAGAIN) {
+        if (!did_anything) {
+          return false;
+        }
+        wait_for_write(fd);
+        continue;
+      }
+      return false;
+    }
+    assert(res != 0);
+    did_anything = true;
+    buf += res;
+    len -= res;
+  }
+  return true;
+}
+
+class AbstractLine;
+
+template <class Node> class AtomicList;
+
+// CRTP
+template <class Node> class AtomicNode {
+  friend class AtomicList<Node>;
+
+  std::atomic<Node *> intrusive_link_next;
+  std::atomic<Node *> intrusive_link_prev;
+
+  AtomicNode(Node *next, Node *prev);
+
+public:
+  // Node is initially unattached
+  AtomicNode();
+  ~AtomicNode();
+};
+
+// A non-owning intrusive linked list,
+// using atomics to ensure thread- and signal- safety.
+template <class Node> class AtomicList {
+  AtomicNode<Node> meta;
+
+public:
+  AtomicList();
+  ~AtomicList();
+
+  void append(Node *node);
+};
+
+template <class Node> AtomicNode<Node>::AtomicNode(Node *next, Node *prev) {
+  intrusive_link_next.store(next);
+  intrusive_link_prev.store(prev);
+}
+template <class Node> AtomicNode<Node>::AtomicNode() {
+  intrusive_link_next.store(nullptr);
+  intrusive_link_prev.store(nullptr);
+}
+template <class Node> AtomicNode<Node>::~AtomicNode() {}
+
+class AbstractLine : public AtomicNode<AbstractLine> {
+  friend class Sink;
+
+  struct {
+    bool dirty : 1;
+  } flags;
+
+public:
+  AbstractLine() : flags{} {}
+  // Due to how vtables work, it is cheaper to *not* inline this.
+  virtual ~AbstractLine(){};
+
+  virtual void write(int fd) = 0;
+
+protected:
+  void not_dirty() { this->flags.dirty = false; }
+};
+
+class StaticTextLine : public AbstractLine {
+  const char *text;
+
+public:
+  template <size_t n> StaticTextLine(const char (&lit)[n]) : text(lit) {}
+  void write(int fd) override {
+    bool ok = write_harder(fd, this->text, strlen(this->text));
+    if (ok)
+      this->not_dirty();
+  }
+};
+
+struct SinkOptions {
+  // Only mandatory field. Everything else can just be zeroed.
+  int fd;
+
+  int tty_width;
+  int tty_height;
+
+  // Additional options will be added in future.
+  SinkOptions(int fd) : fd(fd){};
+};
+
+class Sink;
+// We do still need a global list of sinks in order to handle signals.
+// This is still a win over making a single global list of AbstractLine
+// instances, since we can skip entirely any Sink which does not express
+// interest in asynchronous updates.
+static AtomicList<Sink> all_sinks;
+
+class Sink : public AtomicNode<Sink> {
+  SinkOptions opts;
+  AtomicList<AbstractLine> lines;
+
+public:
+  explicit Sink(SinkOptions o) : opts(o) { all_sinks.append(this); }
+  Sink(Sink &&) = delete;
+  Sink &operator=(Sink &&) = delete;
+};
+
+Sink standard_sink(SinkOptions(STDERR_FILENO));
+
+// static void wait_for_write(int fd);
+
+// Write a buffer fully or not at all.
+// If false is returned, caller may check errno to see if it's EAGAIN
+// or a real error.
+bool write_harder(int fd, const char *buf, size_t len);
+
+// To more easily maintain the doubly-linked structure, loop to itself
+// rather than using NULL pointers.
+template <class Node>
+AtomicList<Node>::AtomicList()
+    : meta(static_cast<Node *>(&meta), static_cast<Node *>(&meta)) {}
+
+// Nothing to do - we didn't allocate any objects, merely borrow.
+template <class Node> AtomicList<Node>::~AtomicList() {
+  // TODO: We *really* shouldn't get here with a non-empty node set.
+  // Should we set all nodes to NULL to indicate they have?
+  // Otherwise we're stuck with dangling pointers in the edges ...
+}
+
+template <class Node> void AtomicList<Node>::append(Node *node) {
+  (void)node;
+#if 0
+  Node *singular = &meta;
+  // TODO when deletion is added, need to disable that
+  // for the duration of this function (using a spinlock).
+  // But, if the spinlock is currently held by the main thread,
+  // and we are a signal handler, we need to advance anyway.
+  // I think our spinlocks need to be recursive, which means tracking
+  // the thread ID.
+  //
+  // Actually, we might not need to disable it, just make sure it
+  // always points to a later node, perhaps? But what if there
+  // is no later node? Maybe use the list-itself-is-a-node trick?
+  assert (node->intrusive_link_next.load() == nullptr);
+  assert (node->intrusive_link_prev.load() == nullptr);
+  // Note that tail is a pointer to a (usually) NULL pointer.
+  std::atomic<Node *> *tail = this->approx_tail;
+  assert (tail != nullptr);
+  while (true)
+  {
+      // The tail may have moved, if someone else is also appending.
+      while (Node *tmp = tail->load())
+      {
+          tail = &tmp->intrusive_link_next;
+      }
+      // First argument is a reference, but we can't re-use
+      // the new value, because we want to advance.
+      Node *expected = nullptr;
+      if (tail->compare_exchange_weak(expected, node))
+      {
+          break;
+      }
+  }
+  // If we're wrong, nobody cares until the next append,
+  // which will fix this anyway.
+  this->approx_tail = &node->intrusive_link_next;
+#endif
+}
+
+}  // tqdm
diff --git a/native/src/utils/utils.cpp b/dcona/native/src/utils/utils.cpp
similarity index 100%
rename from native/src/utils/utils.cpp
rename to dcona/native/src/utils/utils.cpp
diff --git a/native/src/utils/utils.h b/dcona/native/src/utils/utils.h
similarity index 100%
rename from native/src/utils/utils.h
rename to dcona/native/src/utils/utils.h
diff --git a/exhaustive/hypergeom.py b/exhaustive/hypergeom.py
deleted file mode 100644
index 734e68d..0000000
--- a/exhaustive/hypergeom.py
+++ /dev/null
@@ -1,132 +0,0 @@
-import numpy as np
-import pandas as pd
-import scipy.stats
-
-import sys
-import time
-import tqdm
-import json
-
-# Import python package
-import core.extern
-import core.utils
-
-# Arg parser
-import argparse
-
-
-# Argument parser
-parser = argparse.ArgumentParser()
-parser.add_argument("config_path")
-args = parser.parse_args()
-
-# Load config file
-CONFIG_PATH = args.config_path
-config = json.load(open(CONFIG_PATH, "r"))
-
-DATA_PATH = config["data_path"]
-DESCRIPTION_PATH = config["description_path"]
-# INTERACTION_PATH = config["interaction_path"]
-OUTPUT_DIR_PATH = config["output_dir_path"]
-
-REFERENCE_GROUP = config["reference_group"]
-EXPERIMENTAL_GROUP = config["experimental_group"]
-
-CORRELATION = config["correlation"]
-ALTERNATIVE = config["alternative"]
-PROCESS_NUMBER = config["process_number"]
-
-FDR_THRESHOLD = config["fdr_treshold"]
-
-core.utils.check_directory_existence(OUTPUT_DIR_PATH)
-
-def get_occurrence(array):
-    occurrence = {}
-    for elem in array:
-        if elem in occurrence:
-            occurrence[elem] += 1
-        else:
-            occurrence[elem] = 1
-    
-    return occurrence
-
-# Main part
-report_stat = np.load(OUTPUT_DIR_PATH.rstrip("/") + \
-    "/{}_ztest_stat.npy".format(CORRELATION)
-)
-report_pvalue = np.load(OUTPUT_DIR_PATH.rstrip("/") + \
-    "/{}_ztest_pvalue.npy".format(CORRELATION)
-)
-report_fdr = np.load(OUTPUT_DIR_PATH.rstrip("/") + \
-    "/{}_ztest_fdr.npy".format(CORRELATION)
-)
-molecules = pd.read_csv(DATA_PATH, sep=",", index_col=0).index.to_numpy()
-
-if ALTERNATIVE == "less":
-    indexes = np.where(
-        (report_fdr < FDR_THRESHOLD) & 
-        (report_stat < 0)
-    )[0]
-elif ALTERNATIVE == "greater":
-    indexes = np.where(
-        (report_fdr < FDR_THRESHOLD) & 
-        (report_stat > 0)
-    )[0]
-else:
-    # ALTERNATIVE == "two-sided"
-    indexes = np.where(report_fdr < FDR_THRESHOLD)[0]
-
-report_stat = report_stat[indexes]
-report_pvalue = report_pvalue[indexes]
-report_fdr = report_fdr[indexes]
-source_indexes = []
-target_indexes = []
-for ind in tqdm.tqdm(indexes):
-    s, t = core.extern.paired_index(ind, len(molecules))
-    source_indexes.append(molecules[s])
-    target_indexes.append(molecules[t])
-source_indexes = np.array(source_indexes, dtype=np.str)
-target_indexes = np.array(source_indexes, dtype=np.str)
-
-# Report indexing
-ll = list(source_indexes)
-ll.extend(list(target_indexes))
-report_occurrence = get_occurrence(ll)
-
-report_interaction_number = 0
-for index in report_occurrence:
-    report_interaction_number += report_occurrence[index]
-
-initial_interaction_number = len(molecules) * (len(molecules) - 1) // 2
-
-# Analisis
-output_df = pd.DataFrame()
-output_df["Molecule"] = [index for index in report_occurrence]
-output_df["Diff"] = [report_occurrence[molecule] for molecule in output_df["Molecule"]]
-output_df["Total"] = [len(molecules) - 1 for molecule in output_df["Molecule"]]
-output_df["Proportion"] = output_df["Diff"] / output_df["Total"]
-
-output_df["Pvalue"] = 1 - scipy.stats.hypergeom.cdf(
-    output_df["Diff"] - 1,
-    initial_interaction_number, # M
-    report_interaction_number, # n
-    output_df["Total"] # N
-)
-
-adjusted_pvalue = np.array(output_df["Pvalue"] * len(output_df["Pvalue"]) / \
-    scipy.stats.rankdata(output_df["Pvalue"]))
-adjusted_pvalue[adjusted_pvalue > 1] = 1
-adjusted_pvalue = adjusted_pvalue.flatten()
-output_df["FDR"] = adjusted_pvalue
-output_df = output_df[output_df["FDR"] < FDR_THRESHOLD]
-
-output_df["NegProportion"] = 1 - output_df["Proportion"]
-output_df = output_df.sort_values(["FDR", "Pvalue", "NegProportion"])
-output_df = output_df.drop(columns=["NegProportion"])
-
-output_df.to_csv(
-    OUTPUT_DIR_PATH.rstrip("/") +
-    "/{}_{}_hypergeom.csv".format(CORRELATION, ALTERNATIVE),
-    sep=",",
-    index=None
-)
diff --git a/exhaustive/warmup.py b/exhaustive/warmup.py
deleted file mode 100644
index 7a68cdd..0000000
--- a/exhaustive/warmup.py
+++ /dev/null
@@ -1,145 +0,0 @@
-import numpy as np
-import pandas as pd
-import scipy.stats
-
-import sys
-import time
-import tqdm
-import json
-
-import time
-
-# Import python package
-import core.extern
-import core.utils
-
-# Arg parser
-import argparse
-
-
-# Argument parser
-parser = argparse.ArgumentParser()
-parser.add_argument("config_path")
-args = parser.parse_args()
-
-# Load config file
-CONFIG_PATH = args.config_path
-config = json.load(open(CONFIG_PATH, "r"))
-
-DATA_PATH = config["data_path"]
-DESCRIPTION_PATH = config["description_path"]
-# INTERACTION_PATH = config["interaction_path"]
-OUTPUT_DIR_PATH = config["output_dir_path"]
-
-REFERENCE_GROUP = config["reference_group"]
-EXPERIMENTAL_GROUP = config["experimental_group"]
-
-CORRELATION = config["correlation"]
-ALTERNATIVE = config["alternative"]
-SCORE = config["score"]
-
-# REPEATS_NUMBER = config["repeats_number"]
-REPEATS_NUMBER = 10
-PROCESS_NUMBER = config["process_number"]
-
-FDR_THRESHOLD = config["fdr_treshold"]
-
-core.utils.check_directory_existence(OUTPUT_DIR_PATH)
-
-# Main part
-data_df = pd.read_csv(DATA_PATH, sep=",", index_col=0)
-description_df = pd.read_csv(DESCRIPTION_PATH, sep=",")
-
-if (CORRELATION == "spearman"):
-    correlation = core.extern.spearmanr
-elif (CORRELATION == "pearson"):
-    correlation = core.extern.pearsonr
-elif (CORRELATION == "spearman_test"):
-    CORRELATION = "spearman"
-    correlation = core.extern.spearmanr_test
-
-interaction_df=None
-source_indexes=None
-target_indexes=None
-
-reference_indexes = description_df.loc[
-    description_df["Group"] == REFERENCE_GROUP,
-    "Sample"
-].to_list()
-experimental_indexes = description_df.loc[
-    description_df["Group"] == EXPERIMENTAL_GROUP,
-    "Sample"
-].to_list()
-
-# Test mode
-# data_df = data_df.iloc[:100]
-
-print("Bootstrap phase")
-start = time.time()
-
-sources, scores, pvalues = \
-core.extern.score_pipeline(
-    data_df,
-    reference_indexes,
-    experimental_indexes,
-    source_indexes,
-    target_indexes,
-    correlation=CORRELATION,
-    score=SCORE,
-    alternative=ALTERNATIVE,
-    repeats_num=REPEATS_NUMBER,
-    process_num=PROCESS_NUMBER
-)
-end = time.time()
-print("Computational time: {:.3f}".format(end - start))
-
-# Expected computation time
-def score_comlexity(size):
-    if SCORE == "mean":
-        return size
-    else:
-        return size * np.log(size)
-
-def correlation_complexity(sample_size):
-    if CORRELATION == "spearman":
-        return sample_size
-    else:
-        return sample_size * np.log(sample_size)
-
-warmup_computation_complexity = \
-    correlation_complexity(len(reference_indexes)) * len(data_df)**2 // 2 + \
-    correlation_complexity(len(experimental_indexes)) * len(data_df)**2 // 2 + \
-    len(data_df) * score_comlexity(len(data_df))
-
-computation_time = []
-for i in tqdm.tqdm(range(100, len(data_df), 100)):
-    pairs_number = i**2 // 2
-    computation_complexity = \
-        correlation_complexity(len(reference_indexes)) * pairs_number + \
-        correlation_complexity(len(experimental_indexes)) * pairs_number + \
-        i * score_comlexity(i)
-    
-    recommended_repeats_number = i / FDR_THRESHOLD
-    
-    time_complexity = (end - start) * \
-        computation_complexity / warmup_computation_complexity * \
-        recommended_repeats_number / REPEATS_NUMBER / 360
-
-    computation_time.append(time_complexity)
-
-import matplotlib.pyplot as plt
-
-computation_time = np.array(computation_time)
-computation_time = computation_time[computation_time < 24]
-plt.plot(
-    list(range(100, len(data_df), 100))[:len(computation_time)],
-    computation_time
-)
-
-plt.ylabel("Time (hours)")
-plt.xlabel("Number of genes")
-
-plt.savefig(
-    OUTPUT_DIR_PATH.rstrip("/") + \
-    "/{}_warmup.png".format(CORRELATION),
-)
diff --git a/exhaustive/zscore.py b/exhaustive/zscore.py
deleted file mode 100644
index 2b134c0..0000000
--- a/exhaustive/zscore.py
+++ /dev/null
@@ -1,113 +0,0 @@
-import numpy as np
-import pandas as pd
-import scipy.stats
-
-import sys
-import time
-import tqdm
-import json
-
-import time
-
-# Import python package
-import core.extern
-import core.utils
-
-# Arg parser
-import argparse
-
-
-# Argument parser
-parser = argparse.ArgumentParser()
-parser.add_argument("config_path")
-args = parser.parse_args()
-
-# Load config file
-CONFIG_PATH = args.config_path
-config = json.load(open(CONFIG_PATH, "r"))
-
-DATA_PATH = config["data_path"]
-DESCRIPTION_PATH = config["description_path"]
-# INTERACTION_PATH = config["interaction_path"]
-OUTPUT_DIR_PATH = config["output_dir_path"]
-
-REFERENCE_GROUP = config["reference_group"]
-EXPERIMENTAL_GROUP = config["experimental_group"]
-
-CORRELATION = config["correlation"]
-ALTERNATIVE = config["alternative"]
-SCORE = config["score"]
-
-REPEATS_NUMBER = config["repeats_number"]
-PROCESS_NUMBER = config["process_number"]
-
-FDR_THRESHOLD = config["fdr_treshold"]
-
-core.utils.check_directory_existence(OUTPUT_DIR_PATH)
-
-# Main part
-data_df = pd.read_csv(DATA_PATH, sep=",", index_col=0)
-description_df = pd.read_csv(DESCRIPTION_PATH, sep=",")
-
-if (CORRELATION == "spearman"):
-    correlation = core.extern.spearmanr
-elif (CORRELATION == "pearson"):
-    correlation = core.extern.pearsonr
-elif (CORRELATION == "spearman_test"):
-    CORRELATION = "spearman"
-    correlation = core.extern.spearmanr_test
-
-interaction_df=None
-source_indexes=None
-target_indexes=None
-
-reference_indexes = description_df.loc[
-    description_df["Group"] == REFERENCE_GROUP,
-    "Sample"
-].to_list()
-experimental_indexes = description_df.loc[
-    description_df["Group"] == EXPERIMENTAL_GROUP,
-    "Sample"
-].to_list()
-
-# Test mode
-# data_df = data_df.iloc[:100]
-
-print("Bootstrap phase")
-start = time.time()
-
-sources, scores, pvalues = \
-core.extern.score_pipeline(
-    data_df,
-    reference_indexes,
-    experimental_indexes,
-    source_indexes,
-    target_indexes,
-    correlation=CORRELATION,
-    score=SCORE,
-    alternative=ALTERNATIVE,
-    repeats_num=REPEATS_NUMBER,
-    process_num=PROCESS_NUMBER
-)
-print("Computational time: {:.3f}".format(time.time() - start))
-
-adjusted_pvalue = pvalues * len(pvalues) / \
-    scipy.stats.rankdata(pvalues)
-adjusted_pvalue[adjusted_pvalue > 1] = 1
-adjusted_pvalue = adjusted_pvalue.flatten()
-
-# Generate report
-print("Report phase")
-output_df = pd.DataFrame()
-output_df["Source"] = data_df.index.to_numpy()[sources]
-output_df["Score"] = scores
-output_df["Pvalue"] = pvalues
-output_df["FDR"] = adjusted_pvalue
-output_df = output_df.sort_values(["FDR", "Pvalue"])
-
-output_df.to_csv(
-    OUTPUT_DIR_PATH.rstrip("/") + \
-    "/{}_{}_{}_zscore.csv".format(CORRELATION, SCORE, ALTERNATIVE),
-    sep=",",
-    index=None
-)
diff --git a/native/extern/pybind11 b/native/extern/pybind11
deleted file mode 160000
index 70a58c5..0000000
--- a/native/extern/pybind11
+++ /dev/null
@@ -1,1 +0,0 @@
-Subproject commit 70a58c577eaf067748c2ec31bfd0b0a614cffba6
diff --git a/network/hypergeom.py b/network/hypergeom.py
deleted file mode 100644
index 513e658..0000000
--- a/network/hypergeom.py
+++ /dev/null
@@ -1,136 +0,0 @@
-import numpy as np
-import pandas as pd
-import scipy.stats
-
-import sys
-import time
-import tqdm
-import json
-
-# Import python package
-import core.extern
-
-# Arg parser
-import argparse
-
-
-# Argument parser
-parser = argparse.ArgumentParser()
-parser.add_argument("config_path")
-parser.add_argument("-o", "--oriented", action="store_true")
-args = parser.parse_args()
-
-# Load config file
-CONFIG_PATH = args.config_path
-config = json.load(open(CONFIG_PATH, "r"))
-
-DATA_PATH = config["data_path"]
-DESCRIPTION_PATH = config["description_path"]
-INTERACTION_PATH = config["interaction_path"]
-OUTPUT_DIR_PATH = config["output_dir_path"]
-
-REFERENCE_GROUP = config["reference_group"]
-EXPERIMENTAL_GROUP = config["experimental_group"]
-
-CORRELATION = config["correlation"]
-ALTERNATIVE = config["alternative"]
-PROCESS_NUMBER = config["process_number"]
-
-FDR_THRESHOLD = config["fdr_treshold"]
-
-ORIENTED = args.oriented
-
-def get_occurrence(array):
-    occurrence = {}
-    for elem in array:
-        if elem in occurrence:
-            occurrence[elem] += 1
-        else:
-            occurrence[elem] = 1
-    
-    return occurrence
-
-
-# Main part
-data_df = pd.read_csv(DATA_PATH, sep=",", index_col=0)
-
-report_df = pd.read_csv(
-    OUTPUT_DIR_PATH.rstrip("/") + 
-    # "/{}_report.csv".format(CORRELATION),
-    "/{}_ztest.csv".format(CORRELATION),
-    sep=","
-)
-if ALTERNATIVE == "less":
-    report_df = report_df[
-        (report_df["FDR"] < FDR_THRESHOLD) &
-        (report_df["Statistic"] < 0)
-    ]
-elif ALTERNATIVE == "greater":
-    report_df = report_df[
-        (report_df["FDR"] < FDR_THRESHOLD) &
-        (report_df["Statistic"] > 0)
-    ]
-else:
-    # ALTERNATIVE == "two-sided"
-    report_df = report_df[report_df["FDR"] < FDR_THRESHOLD]
-
-# Report indexing
-if ORIENTED:
-    report_occurrence = get_occurrence(list(report_df["Source"]))
-else:
-    ll = list(report_df["Source"])
-    ll.extend(list(report_df["Target"]))
-    report_occurrence = get_occurrence(ll)
-
-report_interaction_number = 0
-for index in report_occurrence:
-    report_interaction_number += report_occurrence[index]
-
-# Description indexing
-interaction_df = pd.read_csv(INTERACTION_PATH, sep=",")
-
-data_molecules = set(data_df.index.to_list())
-interaction_df = interaction_df[interaction_df["Source"].isin(data_molecules)]
-interaction_df = interaction_df[interaction_df["Target"].isin(data_molecules)]
-
-if ORIENTED:
-    initial_occurrence = get_occurrence(list(interaction_df["Source"]))
-else:
-    ll = list(interaction_df["Source"])
-    ll.extend(list(interaction_df["Target"]))
-    initial_occurrence = get_occurrence(ll)
-
-initial_interaction_number = 0
-for index in initial_occurrence:
-    initial_interaction_number += initial_occurrence[index]
-
-# Analisis
-output_df = pd.DataFrame()
-output_df["Molecule"] = [index for index in report_occurrence]
-output_df["Diff"] = [report_occurrence[molecule] for molecule in output_df["Molecule"]]
-output_df["Total"] = [initial_occurrence[molecule] for molecule in output_df["Molecule"]]
-output_df["Proportion"] = output_df["Diff"] / output_df["Total"]
-output_df["Pvalue"] = 1 - scipy.stats.hypergeom.cdf(
-    output_df["Diff"] - 1,
-    initial_interaction_number,
-    report_interaction_number,
-    output_df["Total"]
-)
-
-adjusted_pvalue = np.array(output_df["Pvalue"] * len(output_df["Pvalue"]) / \
-    scipy.stats.rankdata(output_df["Pvalue"]))
-adjusted_pvalue[adjusted_pvalue > 1] = 1
-adjusted_pvalue = adjusted_pvalue.flatten()
-output_df["FDR"] = adjusted_pvalue
-output_df = output_df[output_df["FDR"] < FDR_THRESHOLD]
-
-output_df["NegProportion"] = 1 - output_df["Proportion"]
-output_df = output_df.sort_values(["FDR", "Pvalue", "NegProportion"])
-output_df = output_df.drop(columns=["NegProportion"])
-
-output_df.to_csv(
-    OUTPUT_DIR_PATH.rstrip("/") +
-    "/{}_{}_hypergeom.csv".format(CORRELATION, ALTERNATIVE),
-    sep=",",
-    index=None
-)
diff --git a/network/zscore.py b/network/zscore.py
deleted file mode 100644
index e01eb90..0000000
--- a/network/zscore.py
+++ /dev/null
@@ -1,121 +0,0 @@
-import numpy as np
-import pandas as pd
-import scipy.stats
-
-import sys
-import time
-import tqdm
-import json
-
-import time
-
-# Import python package
-import core.extern
-
-# Arg parser
-import argparse
-
-
-# Argument parser
-parser = argparse.ArgumentParser()
-parser.add_argument("config_path")
-args = parser.parse_args()
-
-# Load config file
-CONFIG_PATH = args.config_path
-config = json.load(open(CONFIG_PATH, "r"))
-
-DATA_PATH = config["data_path"]
-DESCRIPTION_PATH = config["description_path"]
-INTERACTION_PATH = config["interaction_path"]
-OUTPUT_DIR_PATH = config["output_dir_path"]
-
-REFERENCE_GROUP = config["reference_group"]
-EXPERIMENTAL_GROUP = config["experimental_group"]
-
-CORRELATION = config["correlation"]
-ALTERNATIVE = config["alternative"]
-SCORE = config["score"]
-
-REPEATS_NUMBER = config["repeats_number"]
-PROCESS_NUMBER = config["process_number"]
-
-FDR_THRESHOLD = config["fdr_treshold"]
-
-
-# Main part
-data_df = pd.read_csv(DATA_PATH, sep=",", index_col=0)
-description_df = pd.read_csv(DESCRIPTION_PATH, sep=",")
-
-if (CORRELATION == "spearman"):
-    correlation = core.extern.spearmanr
-elif (CORRELATION == "pearson"):
-    correlation = core.extern.pearsonr
-elif (CORRELATION == "spearman_test"):
-    CORRELATION = "spearman"
-    correlation = core.extern.spearmanr_test
-
-interaction_df = pd.read_csv(INTERACTION_PATH, sep=",")
-
-data_molecules = set(data_df.index.to_list())
-interaction_df = interaction_df[interaction_df["Source"].isin(data_molecules)]
-interaction_df = interaction_df[interaction_df["Target"].isin(data_molecules)]
-
-# Test mode
-# interaction_df = interaction_df.iloc[:10]
-
-source_indexes = interaction_df["Source"]
-target_indexes = interaction_df["Target"]
-
-data_molecules = data_molecules.intersection(
-    set(source_indexes) | set(target_indexes)
-)
-data_df = data_df.loc[data_df.index.isin(data_molecules)]
-
-reference_indexes = description_df.loc[
-    description_df["Group"] == REFERENCE_GROUP,
-    "Sample"
-].to_list()
-experimental_indexes = description_df.loc[
-    description_df["Group"] == EXPERIMENTAL_GROUP,
-    "Sample"
-].to_list()
-
-print("Bootstrap phase")
-start = time.time()
-
-sources, scores, pvalues = \
-core.extern.score_pipeline(
-    data_df,
-    reference_indexes,
-    experimental_indexes,
-    source_indexes,
-    target_indexes,
-    correlation=CORRELATION,
-    score=SCORE,
-    alternative=ALTERNATIVE,
-    repeats_num=REPEATS_NUMBER,
-    process_num=PROCESS_NUMBER
-)
-print(time.time() - start)
-
-adjusted_pvalue = pvalues * len(pvalues) / \
-    scipy.stats.rankdata(pvalues)
-adjusted_pvalue[adjusted_pvalue > 1] = 1
-adjusted_pvalue = adjusted_pvalue.flatten()
-
-# Generate report
-print("Report phase")
-output_df = pd.DataFrame()
-output_df["Source"] = data_df.index.to_numpy()[sources]
-output_df["Score"] = scores
-output_df["Pvalue"] = pvalues
-output_df["FDR"] = adjusted_pvalue
-output_df = output_df.sort_values(["FDR", "Pvalue"])
-
-output_df.to_csv(
-    OUTPUT_DIR_PATH.rstrip("/") + \
-    "/{}_{}_{}_zscore.csv".format(CORRELATION, SCORE, ALTERNATIVE),
-    sep=",",
-    index=None
-)
diff --git a/pyproject.toml b/pyproject.toml
new file mode 100644
index 0000000..eb462ea
--- /dev/null
+++ b/pyproject.toml
@@ -0,0 +1,5 @@
+[build-system]
+requires = ["pybind11 @ git+https://github.com/pybind/pybind11.git@master", 
+			"setuptools",
+			"wheel"]
+build-backend = "setuptools.build_meta"
\ No newline at end of file
diff --git a/setup.py b/setup.py
new file mode 100644
index 0000000..83db15e
--- /dev/null
+++ b/setup.py
@@ -0,0 +1,89 @@
+try:
+    import pybind11
+    from pybind11.setup_helpers import Pybind11Extension
+except ImportError:
+    from setuptools import Extension as Pybind11Extension
+
+from setuptools import setup, Extension, find_packages
+
+with open("README.md", "r") as readme_file:
+    readme = readme_file.read()
+
+requirements = ['pybind11', 'pandas', 'numpy', 'scipy', 'tqdm']
+
+
+__version__ = "0.0.9"
+
+
+src = 'dcona/native/src/'
+
+ext_modules = [
+    Extension("dcona.core.extern.utils",
+        ["dcona/native/src/utils/utils.cpp",
+         "dcona/native/src/lib/utils.cpp"],
+         include_dirs=[pybind11.get_include()],
+         language='c++',
+         extra_compile_args=['-std=c++11', '-O3', '-Wall']
+    ),
+
+    Extension("dcona.core.extern.correlations",
+        ["dcona/native/src/utils/utils.cpp",
+         "dcona/native/src/lib/correlations.cpp",
+         "dcona/native/src/correlations/correlations.cpp"],
+         include_dirs=[pybind11.get_include()],
+         language='c++',
+         extra_compile_args=['-std=c++11', '-O3', '-Wall']
+    ),
+        
+    Extension("dcona.core.extern.tests",
+        ["dcona/native/src/utils/utils.cpp",
+         "dcona/native/src/correlations/correlations.cpp",
+         "dcona/native/src/tests/tests.cpp",
+         "dcona/native/src/lib/tests.cpp"],
+         include_dirs=[pybind11.get_include()],
+         language='c++',
+         extra_compile_args=['-std=c++11', '-O3', '-Wall']
+    ),
+        
+    Extension("dcona.core.extern.scores",
+        ["dcona/native/src/utils/utils.cpp",
+         "dcona/native/src/scores/scores.cpp",
+         "dcona/native/src/lib/scores.cpp"],
+         include_dirs=[pybind11.get_include()],
+         language='c++',
+         extra_compile_args=['-std=c++11', '-O3', '-Wall']
+    ),
+        
+    Extension("dcona.core.extern.pipelines",
+        ["dcona/native/src/utils/utils.cpp",
+         "dcona/native/src/correlations/correlations.cpp",
+         "dcona/native/src/tests/tests.cpp",
+         "dcona/native/src/scores/scores.cpp",
+         "dcona/native/src/pipelines/pipelines.cpp",
+         "dcona/native/src/lib/pipelines.cpp"],
+         include_dirs=[pybind11.get_include()],
+         language='c++',
+         extra_compile_args=['-std=c++11', '-O3', '-Wall']
+    ),
+]
+
+
+setup(
+    name="dcona",
+    version=__version__,
+    author="Anton Zhiyanov, Narek Engibaryan",
+    author_email="zhiyanovap@gmail.com",
+    url="https://github.com/zhiyanov/DCoNA",
+    description="Differential Correlation Network Analysis",
+    long_description=readme,
+    long_description_content_type="text/markdown",
+    packages=find_packages(),
+    install_requires=requirements,
+    ext_modules=ext_modules,
+    python_requires=">=3.8",
+    entry_points={
+        'console_scripts': [
+            'dcona = dcona.__main__:main',
+        ],
+    },
+)


Test patch:
diff --git a/core/correlations/correlation_tests.py b/dcona/core/correlations/correlation_tests.py
similarity index 100%
rename from core/correlations/correlation_tests.py
rename to dcona/core/correlations/correlation_tests.py
diff --git a/core/extern/ptests.py b/dcona/core/extern/ptests.py
similarity index 100%
rename from core/extern/ptests.py
rename to dcona/core/extern/ptests.py
diff --git a/dcona/lib/ztest.py b/dcona/lib/ztest.py
new file mode 100644
index 0000000..7025793
--- /dev/null
+++ b/dcona/lib/ztest.py
@@ -0,0 +1,227 @@
+import numpy as np
+import pandas as pd
+import scipy.stats
+
+from ..core import extern
+from . import utils
+from . import dump
+
+
+def ztest(
+    data_df, description_df,
+    reference_group, experimental_group,
+    correlation="spearman", alternative="two-sided",
+    interaction=None,
+    repeats_number=None,
+    output_dir=None,
+    process_number=None
+):
+    if process_number is None:
+        process_number = 1
+
+    # if gene names are in dataframe column, relocate them in df.index
+    # if not pd.api.types.is_number(data_df.iloc[0, 0]):
+    #     data_df = data_df.copy()
+    #     data_df.set_index(data_df.columns[0], inplace=True)
+        
+    if interaction is not None:
+        if isinstance(interaction, pd.core.frame.DataFrame):
+            interaction_df = interaction
+        else:
+            interaction_df = utils.generate_pairs(interaction)
+
+        if repeats_number is None:
+            repeats_number = int(len(interaction_df) / 0.05)
+    else:
+        interaction_df = None
+
+        if repeats_number is None:
+            repeats_number = 0
+    
+    sorted_indexes, indexes, df_indexes, \
+    ref_corrs, ref_pvalues, exp_corrs, exp_pvalues, \
+    stat, pvalue, adjusted_pvalue, \
+    boot_pvalue = _ztest(
+        data_df, description_df, interaction_df, \
+        reference_group, experimental_group, \
+        correlation, alternative, \
+        repeats_number, process_number
+    )
+    
+    if output_dir:
+        dump.check_directory_existence(output_dir)
+        
+        if repeats_number > 0: 
+            df_template = pd.DataFrame(columns=[
+                "Source", "Target", "RefCorr", "RefPvalue", 
+                "ExpCorr", "ExpPvalue", "Statistic",
+                "Pvalue", "AdjPvalue", "PermutePvalue"
+            ])
+
+            df_columns = [
+                ref_corrs, ref_pvalues,
+                exp_corrs, exp_pvalues, stat,
+                pvalue, adjusted_pvalue, boot_pvalue
+            ]
+        else:
+            df_template = pd.DataFrame(columns=[
+                "Source", "Target", "RefCorr", "RefPvalue", 
+                "ExpCorr", "ExpPvalue", "Statistic",
+                "Pvalue", "AdjPvalue"
+            ])
+
+            df_columns = [
+                ref_corrs, ref_pvalues,
+                exp_corrs, exp_pvalues, stat,
+                pvalue, adjusted_pvalue
+            ]
+
+
+        path_to_file = output_dir.rstrip("/") + f"/{correlation}_ztest.csv"
+        dump.save_by_chunks(
+            sorted_indexes,
+            df_indexes, df_template, df_columns,
+            path_to_file,
+            index_transform=indexes
+        )
+        
+        print(f"File saved at: {path_to_file}")
+        return None
+    else:
+        source_indexes = []
+        target_indexes = []
+
+        for ind in sorted_indexes:
+            s, t = extern.paired_index(ind, len(df_indexes))
+            source_indexes.append(df_indexes[s])
+            target_indexes.append(df_indexes[t])
+        
+        if repeats_number > 0:
+            output_df = pd.DataFrame(data={
+                "Source": source_indexes,
+                "Target": target_indexes,
+                "RefCorr": ref_corrs[sorted_indexes], 
+                "RefPvalue": ref_pvalues[sorted_indexes], 
+                "ExpCorr": exp_corrs[sorted_indexes], 
+                "ExpPvalue": exp_pvalues[sorted_indexes], 
+                "Statistic": stat[sorted_indexes],
+                "Pvalue": pvalue[sorted_indexes], 
+                "AdjPvalue": adjusted_pvalue[sorted_indexes],
+                "PermutePvalue": boot_pvalue[sorted_indexes] 
+            })
+        else:
+            output_df = pd.DataFrame(data={
+                "Source": source_indexes,
+                "Target": target_indexes,
+                "RefCorr": ref_corrs[sorted_indexes], 
+                "RefPvalue": ref_pvalues[sorted_indexes], 
+                "ExpCorr": exp_corrs[sorted_indexes], 
+                "ExpPvalue": exp_pvalues[sorted_indexes], 
+                "Statistic": stat[sorted_indexes],
+                "Pvalue": pvalue[sorted_indexes], 
+                "AdjPvalue": adjusted_pvalue[sorted_indexes]
+            })
+
+        return output_df
+
+def _ztest(
+    data_df, description_df, interaction_df,
+    reference_group, experimental_group,
+    correlation, alternative,
+    repeats_number, process_number
+):
+    if (correlation != "spearman"):
+        correlation = "pearson"
+
+    if interaction_df is not None:
+        data_molecules = set(data_df.index.to_list())
+        interaction_df = interaction_df[interaction_df["Source"].isin(data_molecules)]
+        interaction_df = interaction_df[interaction_df["Target"].isin(data_molecules)]
+
+        source_indexes = interaction_df["Source"]
+        target_indexes = interaction_df["Target"]
+
+        data_molecules = data_molecules.intersection(
+            set(source_indexes) | set(target_indexes)
+        )
+        data_df = data_df.loc[data_df.index.isin(data_molecules)]
+    else:
+        interaction_df = None
+        source_indexes = None
+        target_indexes = None
+
+    reference_indexes = description_df.loc[
+        description_df["Group"] == reference_group,
+        "Sample"
+    ].to_list()
+    experimental_indexes = description_df.loc[
+        description_df["Group"] == experimental_group,
+        "Sample"
+    ].to_list()
+
+    print("Computation of z-tests")
+    ref_corrs, ref_pvalues, \
+    exp_corrs, exp_pvalues, \
+    stat, pvalue, boot_pvalue = \
+    extern.ztest_pipeline(
+        data_df,
+        reference_indexes,
+        experimental_indexes,
+        source_indexes,
+        target_indexes,
+        correlation=correlation,
+        alternative=alternative,
+        repeats_num=repeats_number,
+        process_num=process_number,
+        correlation_alternative="two-sided"
+    )
+
+    print("Computation of adjusted p-value")
+    adjusted_pvalue = pvalue * len(pvalue) / \
+        scipy.stats.rankdata(pvalue)
+    adjusted_pvalue[adjusted_pvalue > 1] = 1
+    adjusted_pvalue = adjusted_pvalue.flatten()
+
+    indexes = np.arange(len(adjusted_pvalue))
+
+    ref_corrs = ref_corrs[indexes]
+    ref_pvalues = ref_pvalues[indexes]
+
+    exp_corrs = exp_corrs[indexes]
+    exp_pvalues = exp_pvalues[indexes]
+
+    stat = stat[indexes]
+    pvalue = pvalue[indexes]
+    boot_pvalue = boot_pvalue[indexes]
+
+    adjusted_pvalue = adjusted_pvalue[indexes]
+    df_indexes = data_df.index.to_numpy()
+
+    ref_corrs = ref_corrs[indexes]
+    ref_pvalues = ref_pvalues[indexes]
+
+    exp_corrs = exp_corrs[indexes]
+    exp_pvalues = exp_pvalues[indexes]
+
+    stat = stat[indexes]
+    pvalue = pvalue[indexes]
+    boot_pvalue = boot_pvalue[indexes]
+    adjusted_pvalue = adjusted_pvalue[indexes]
+    df_indexes = data_df.index.to_numpy()
+
+    print("Creating (adj-pv, p-value) array")
+    FDR_pvalue = np.core.records.fromarrays(
+        [adjusted_pvalue, pvalue],
+        names='FDR, pvalue'
+    )
+
+    print("Sorting (adj-pv, p-value) array")
+    sorted_indexes = np.argsort(FDR_pvalue, order=('FDR', 'pvalue'))
+    del FDR_pvalue
+    print("(adj-pv, p-value) array is sorted")
+
+    return sorted_indexes, indexes, df_indexes, \
+        ref_corrs, ref_pvalues, exp_corrs, exp_pvalues, \
+        stat, pvalue, adjusted_pvalue, \
+        boot_pvalue
+
diff --git a/native/src/lib/tests.cpp b/dcona/native/src/lib/tests.cpp
similarity index 100%
rename from native/src/lib/tests.cpp
rename to dcona/native/src/lib/tests.cpp
diff --git a/native/src/tests/tests.cpp b/dcona/native/src/tests/tests.cpp
similarity index 100%
rename from native/src/tests/tests.cpp
rename to dcona/native/src/tests/tests.cpp
diff --git a/native/src/tests/tests.h b/dcona/native/src/tests/tests.h
similarity index 100%
rename from native/src/tests/tests.h
rename to dcona/native/src/tests/tests.h
diff --git a/exhaustive/ztest.py b/exhaustive/ztest.py
deleted file mode 100644
index b5a0afc..0000000
--- a/exhaustive/ztest.py
+++ /dev/null
@@ -1,193 +0,0 @@
-import numpy as np
-import pandas as pd
-import scipy.stats
-
-import sys
-import time
-import tqdm
-import json
-import os
-import time
-
-# Import python package
-import core.extern
-import core.utils
-
-# Arg parser
-import argparse
-
-
-# Argument parser
-parser = argparse.ArgumentParser()
-parser.add_argument("config_path")
-args = parser.parse_args()
-
-# Load config file
-CONFIG_PATH = args.config_path
-config = json.load(open(CONFIG_PATH, "r"))
-
-DATA_PATH = config["data_path"]
-DESCRIPTION_PATH = config["description_path"]
-# INTERACTION_PATH = config["interaction_path"]
-OUTPUT_DIR_PATH = config["output_dir_path"]
-
-REFERENCE_GROUP = config["reference_group"]
-EXPERIMENTAL_GROUP = config["experimental_group"]
-
-CORRELATION = config["correlation"]
-ALTERNATIVE = config["alternative"]
-REPEATS_NUMBER = config["repeats_number"]
-PROCESS_NUMBER = config["process_number"]
-
-FDR_THRESHOLD = config["fdr_treshold"]
-
-core.utils.check_directory_existence(OUTPUT_DIR_PATH)
-
-# Main part
-data_df = pd.read_csv(DATA_PATH, sep=",", index_col=0)
-description_df = pd.read_csv(DESCRIPTION_PATH, sep=",")
-
-if (CORRELATION == "spearman"):
-    correlation = core.extern.spearmanr
-elif (CORRELATION == "pearson"):
-    correlation = core.extern.pearsonr
-elif (CORRELATION == "spearman_test"):
-    CORRELATION = "spearman"
-    correlation = core.extern.spearmanr_test
-
-interaction_df=None
-source_indexes=None
-target_indexes=None
-
-reference_indexes = description_df.loc[
-    description_df["Group"] == REFERENCE_GROUP,
-    "Sample"
-].to_list()
-experimental_indexes = description_df.loc[
-    description_df["Group"] == EXPERIMENTAL_GROUP,
-    "Sample"
-].to_list()
-
-# Test mode
-# data_df = data_df.iloc[:1000]
-
-print("Bootstrap phase")
-start = time.time()
-
-ref_corrs, ref_pvalues, \
-exp_corrs, exp_pvalues, \
-stat, pvalue, boot_pvalue = \
-core.extern.ztest_pipeline(
-    data_df,
-    reference_indexes,
-    experimental_indexes,
-    source_indexes,
-    target_indexes,
-    correlation=CORRELATION,
-    alternative=ALTERNATIVE,
-    repeats_num=REPEATS_NUMBER,
-    process_num=PROCESS_NUMBER,
-    correlation_alternative="two-sided"
-)
-print("Computational time: {:.3f}".format(time.time() - start))
-
-adjusted_pvalue = pvalue * len(pvalue) / \
-    scipy.stats.rankdata(pvalue)
-adjusted_pvalue[adjusted_pvalue > 1] = 1
-adjusted_pvalue = adjusted_pvalue.flatten()
-
-indexes = np.arange(len(adjusted_pvalue))
-
-ref_corrs = ref_corrs[indexes]
-ref_pvalues = ref_pvalues[indexes]
-
-exp_corrs = exp_corrs[indexes]
-exp_pvalues = exp_pvalues[indexes]
-
-stat = stat[indexes]
-pvalue = pvalue[indexes]
-boot_pvalue = boot_pvalue[indexes]
-
-adjusted_pvalue = adjusted_pvalue[indexes]
-df_indexes = data_df.index.to_numpy()
-
-# Store output fast in numpy format
-np.save(OUTPUT_DIR_PATH.rstrip("/") + \
-    "/{}_ztest_ref_corr.npy".format(CORRELATION),
-    ref_corrs
-)
-np.save(OUTPUT_DIR_PATH.rstrip("/") + \
-    "/{}_ztest_ref_pvalue.npy".format(CORRELATION),
-    ref_pvalues
-)
-np.save(OUTPUT_DIR_PATH.rstrip("/") + \
-    "/{}_ztest_exp_corr.npy".format(CORRELATION),
-    exp_corrs
-)
-np.save(OUTPUT_DIR_PATH.rstrip("/") + \
-    "/{}_ztest_exp_pvalue.npy".format(CORRELATION),
-    exp_pvalues
-)
-np.save(OUTPUT_DIR_PATH.rstrip("/") + \
-    "/{}_ztest_stat.npy".format(CORRELATION),
-    stat
-)
-np.save(OUTPUT_DIR_PATH.rstrip("/") + \
-    "/{}_ztest_pvalue.npy".format(CORRELATION),
-    pvalue
-)
-np.save(OUTPUT_DIR_PATH.rstrip("/") + \
-    "/{}_ztest_bootpv.npy".format(CORRELATION),
-    boot_pvalue
-)
-np.save(OUTPUT_DIR_PATH.rstrip("/") + \
-    "/{}_ztest_fdr.npy".format(CORRELATION),
-    adjusted_pvalue
-)
-
-# Remove low fdr interactions 
-indexes = np.where(adjusted_pvalue < FDR_THRESHOLD)[0]
-
-ref_corrs = ref_corrs[indexes]
-ref_pvalues = ref_pvalues[indexes]
-
-exp_corrs = exp_corrs[indexes]
-exp_pvalues = exp_pvalues[indexes]
-
-stat = stat[indexes]
-pvalue = pvalue[indexes]
-boot_pvalue = boot_pvalue[indexes]
-adjusted_pvalue = adjusted_pvalue[indexes]
-df_indexes = data_df.index.to_numpy()
-
-# Save mode
-print("Creating FDR/pvalue array")
-FDR_pvalue = np.core.records.fromarrays(
-    [adjusted_pvalue, pvalue],
-    names='FDR, pvalue'
-)
-
-print("Sorting FDR/pvalue array")
-sorted_indexes = np.argsort(FDR_pvalue, order=('FDR','pvalue'))
-del FDR_pvalue
-print("FDR/pvalue array is sorted")
-
-df_template = pd.DataFrame(columns=[
-    "Source", "Target", "RefCorr", "RefPvalue", 
-    "ExpCorr", "ExpPvalue", "Statistic",
-    "Pvalue", "Bootpv", "FDR"
-])
-df_columns = [
-    ref_corrs, ref_pvalues,
-    exp_corrs, exp_pvalues, stat,
-    pvalue, boot_pvalue, adjusted_pvalue
-]
-
-path_to_file = OUTPUT_DIR_PATH.rstrip("/") + "/{}_ztest.csv".format(CORRELATION)
-core.utils.save_by_chunks(
-    sorted_indexes, 
-    df_indexes, df_template, df_columns,
-    path_to_file,
-    # index_transform=None
-    index_transform=indexes
-)
diff --git a/network/ztest.py b/network/ztest.py
deleted file mode 100644
index c7b0124..0000000
--- a/network/ztest.py
+++ /dev/null
@@ -1,128 +0,0 @@
-import numpy as np
-import pandas as pd
-import scipy.stats
-
-import sys
-import time
-import tqdm
-import json
-
-import time
-
-# Import python package
-import core.extern
-
-# Arg parser
-import argparse
-
-
-# Argument parser
-parser = argparse.ArgumentParser()
-parser.add_argument("config_path")
-args = parser.parse_args()
-
-# Load config file
-CONFIG_PATH = args.config_path
-config = json.load(open(CONFIG_PATH, "r"))
-
-DATA_PATH = config["data_path"]
-DESCRIPTION_PATH = config["description_path"]
-INTERACTION_PATH = config["interaction_path"]
-OUTPUT_DIR_PATH = config["output_dir_path"]
-
-REFERENCE_GROUP = config["reference_group"]
-EXPERIMENTAL_GROUP = config["experimental_group"]
-
-CORRELATION = config["correlation"]
-ALTERNATIVE = config["alternative"]
-# REPEATS_NUMBER = config["repeats_number"]
-REPEATS_NUMBER = 0
-PROCESS_NUMBER = config["process_number"]
-
-FDR_THRESHOLD = config["fdr_treshold"]
-
-
-# Main part
-data_df = pd.read_csv(DATA_PATH, sep=",", index_col=0)
-description_df = pd.read_csv(DESCRIPTION_PATH, sep=",")
-
-if (CORRELATION == "spearman"):
-    correlation = core.extern.spearmanr
-elif (CORRELATION == "pearson"):
-    correlation = core.extern.pearsonr
-elif (CORRELATION == "spearman_test"):
-    CORRELATION = "spearman"
-    correlation = core.extern.spearmanr_test
-
-interaction_df = pd.read_csv(INTERACTION_PATH, sep=",")
-
-data_molecules = set(data_df.index.to_list())
-interaction_df = interaction_df[interaction_df["Source"].isin(data_molecules)]
-interaction_df = interaction_df[interaction_df["Target"].isin(data_molecules)]
-
-# Test mode
-# interaction_df = interaction_df.iloc[:10]
-
-source_indexes = interaction_df["Source"]
-target_indexes = interaction_df["Target"]
-
-data_molecules = data_molecules.intersection(
-    set(source_indexes) | set(target_indexes)
-)
-data_df = data_df.loc[data_df.index.isin(data_molecules)]
-
-reference_indexes = description_df.loc[
-    description_df["Group"] == REFERENCE_GROUP,
-    "Sample"
-].to_list()
-experimental_indexes = description_df.loc[
-    description_df["Group"] == EXPERIMENTAL_GROUP,
-    "Sample"
-].to_list()
-
-print("Bootstrap phase")
-start = time.time()
-ref_corrs, ref_pvalues, \
-exp_corrs, exp_pvalues, \
-stat, pvalue, boot_pvalue = \
-core.extern.ztest_pipeline(
-    data_df,
-    reference_indexes,
-    experimental_indexes,
-    source_indexes,
-    target_indexes,
-    correlation=CORRELATION,
-    alternative=ALTERNATIVE,
-    repeats_num=REPEATS_NUMBER,
-    process_num=PROCESS_NUMBER,
-    correlation_alternative="two-sided"
-)
-print(time.time() - start)
-
-adjusted_pvalue = pvalue * len(pvalue) / \
-    scipy.stats.rankdata(pvalue)
-adjusted_pvalue[adjusted_pvalue > 1] = 1
-adjusted_pvalue = adjusted_pvalue.flatten()
-
-# Generate report
-print("Report phase")
-output_df = pd.DataFrame(interaction_df)
-
-output_df["RefCorr"] = ref_corrs 
-output_df["RefPvalue"] = ref_pvalues 
-
-output_df["ExpCorr"] = exp_corrs 
-output_df["ExpPvalues"] = exp_pvalues
-
-output_df["Statistic"] = stat
-output_df["Pvalue"] = pvalue
-output_df["Bootpv"] = boot_pvalue
-output_df["FDR"] = adjusted_pvalue
-output_df = output_df.sort_values(["FDR", "Pvalue"])
-
-output_df.to_csv(
-    OUTPUT_DIR_PATH.rstrip("/") + \
-    "/{}_ztest.csv".format(CORRELATION),
-    sep=",",
-    index=None
-)