Python package

setup.py

@@ -0,0 +1,25 @@
+from setuptools import setup, find_packages
+
+setup(
+    name="VoA",
+    version="0.0.1",
+    author="Mateusz Przyborowski, Krzysztof Suwada",
+    description="Tools proposed in the paper 'Validation of Association'",
+    long_description=open("README.md").read(),
+    long_description_content_type="text/markdown",
+    url="https://github.com/nexocodecom/VoAPython",
+    packages=find_packages(),
+    install_requires=[
+        "numpy >= 1.23.5",
+        "numba >= 0.58.1",
+        "plotly >= 5.18.0",
+        # Add any other dependencies here
+    ],
+    classifiers=[
+        "Programming Language :: Python :: 3",
+        "License :: OSI Approved :: MIT License",
+        "Operating System :: OS Independent",
+    ],
+    python_requires='>=3.8',
+)
+

tests/test_voa.py

@@ -0,0 +1,37 @@
+import numpy as np
+import voa
+
+def test_copula():
+    true_copula = np.array(
+        [
+            [-0.1505847, -0.2742042, -0.3742406, -0.4714045, -0.5741693, 0.69006556, 0.5741693, 0.4714045, 0.3742406,
+             0.2742042, 0.1505847],
+            [-0.2742042, -0.4993070, -0.6814663, 1.1200012, 0.8286734, 1.25656172, 1.0455226, 0.8583951, 0.6814663,
+             0.4993070, 0.2742042],
+            [-0.3742406, -0.6814663, -0.9300817, 0.4485395, 0.1078143, 0.20579830, 1.4269545, 1.1715584, 0.9300817,
+             0.6814663, 0.3742406],
+            [-0.4714045, -0.8583951, 0.4485395, 1.4395893, 0.9643376, 0.80237742, 1.7974341, 1.4757296, 1.1715584,
+             0.8583951, 0.4714045],
+            [-0.5741693, -1.0455226, 0.1078143, 0.9643376, 0.4270426, 0.05847053, 0.8811133, 0.4165482, -0.1078143,
+             1.0455226, 0.5741693],
+            [0.6900656, -0.5863955, 0.2057983, 0.8023774, 1.3448223, 0.63245553, 1.3448223, 0.8023774, 0.2057983,
+             1.2565617, 0.6900656],
+            [0.5741693, 1.0455226, 1.4269545, 1.7974341, 2.1892691, 1.34482230, 1.7351985, 0.9643376, 0.1078143,
+             0.8286734, -0.5741693],
+            [0.4714045, 0.8583951, 1.1715584, 1.4757296, 1.7974341, 0.80237742, 0.9643376, 1.4395893, 0.4485395,
+             1.1200012, -0.4714045],
+            [0.3742406, 0.6814663, 0.9300817, 1.1715584, 1.4269545, 1.71498585, 1.6425832, 2.0686374, 0.8705564,
+             1.5173982, -0.3742406],
+            [0.2742042, 0.4993070, 0.6814663, 0.8583951, 1.0455226, 1.25656172, 0.8286734, 1.1200012, 1.5173982,
+             2.1858551, -0.2742042],
+            [0.1505847, 0.2742042, 0.3742406, 0.4714045, 0.5741693, 0.69006556, -0.5741693, -0.4714045, -0.3742406,
+             -0.2742042, -0.1505847],
+        ]
+    )
+
+    r = np.array(range(1, 11))
+    s = np.array([3, 6, 2, 9, 4, 1, 7, 5, 8, 10])
+    test_copula = voa.calculate_copula_grid(r, s)
+    assert np.allclose(true_copula, test_copula, atol=1e-6), "Test failed: Copula grid is incorrect"
+
+

voa/__init__.py

@@ -0,0 +1 @@
+from .copula import calculate_copula_grid, create_Q_grid, calculate_copula_mc_grid, create_Q_plot

voap.py → voa/copula.py

@@ -1,10 +1,13 @@
 import numpy as np
 import numba
-# from tqdm.auto import tqdm
+import time
+import plotly.graph_objs as go
+import plotly.io as pio
+from typing import List, Dict, Tuple, Any
 
 
 @numba.njit(parallel=True)
-def Q_function(Q_grid, C_grid):
+def Q_function(Q_grid: np.ndarray, C_grid: np.ndarray) -> Dict[str, np.ndarray]:
     n = len(C_grid) - 1
     z = []
 
@@ -20,31 +23,32 @@ def Q_function(Q_grid, C_grid):
     return {"x": Q_grid[:, 0], "y": Q_grid[:, 1], "z": np.array(z)}
 
 
-def create_Q_grid(n=10):
+def create_Q_grid(n: int = 10) -> np.ndarray:
     u = np.linspace(0.5 / (n + 1), (0.5 + n) / (n + 1), n + 1)
     grid = np.stack(np.meshgrid(u, u)).T.reshape(-1, 2)
 
     return grid
 
+
 @numba.njit(parallel=True)
-def calculate_copula_part(t, c):
+def calculate_copula_part(t: np.ndarray, c: np.ndarray) -> np.ndarray:
     n = len(t)
-    n_inv = 1 / n # n**-1 ## numba doesn't get n**-1 xD
-    c = np.zeros(c.shape)
+    n_inv = 1.0 / float(n)
     for i in numba.prange(1, int(n / 2) + 2):
         for j in numba.prange(int(n / 2) + 2):
-            # c[i, j] = c[i - 1, j] + (j >= t[i - 1]) * n_inv
             c[i, j] = c[0, j] + n_inv * np.sum(j >= t[:i])
     return c
 
+
 @numba.njit
-def wn(n, c, i, j):
+def wn(n: int, c: float, i: int, j: int) -> float:
     u = (i + 0.5) / (n + 1)
     v = (j + 0.5) / (n + 1)
     return n**0.5 * (c - u * v) * (u * v * (1 - u) * (1 - v)) ** -0.5
 
+
 @numba.njit
-def fill_matrix(n, ks, c, x_prim, y_prim):
+def fill_matrix(n: int, ks: np.ndarray, c: np.ndarray, x_prim: bool, y_prim: bool) -> np.ndarray:
     sign = (-1) ** (x_prim + y_prim)
 
     for i in range(int(np.floor((n + 1) / 2)) + 1):
@@ -54,18 +58,16 @@ def fill_matrix(n, ks, c, x_prim, y_prim):
             ks[x, y] = sign * wn(n, c[i, j], i, j)
     return ks
 
+
 @numba.njit
-def arma_copula(rx, ry):
+def arma_copula(rx: np.ndarray, ry: np.ndarray) -> np.ndarray:
     n = len(rx)
 
     ctab = np.zeros((n + 1, n + 1))
     ctabs22 = np.zeros((n + 1, n + 1))
     ctabs12 = np.zeros((n + 1, n + 1))
     ctabs21 = np.zeros((n + 1, n + 1))
 
-    rsx = np.zeros(n)
-    rsy = np.zeros(n)
-
     ks = np.zeros((n + 1, n + 1))
 
     rsx = len(rx) + 1 - rx
@@ -88,14 +90,15 @@ def arma_copula(rx, ry):
 
     return ks
 
+
 @numba.njit
-def calculate_copula_grid(x, y):
+def calculate_copula_grid(x: np.ndarray, y: np.ndarray) -> np.ndarray:
     return arma_copula(
         np.argsort(np.argsort(x)) + 1, np.argsort(np.argsort(y)) + 1
     )
 
-# @numba.njit(parallel=True)
-def calculate_copula_mc_grid(x, y, mc=100, seed=0):
+
+def calculate_copula_mc_grid(x: np.ndarray, y: np.ndarray, mc: int = 100, seed: int = 0) -> np.ndarray:
     rng = np.random.RandomState(seed)
     k = len(x)
     g = np.zeros((k + 1, len(y) + 1))
@@ -107,6 +110,54 @@ def calculate_copula_mc_grid(x, y, mc=100, seed=0):
 
     return g
 
+def create_Q_plot(X: List[float], Y: List[float], k_plot_grid: int = 100, MC: int = 100, display: bool = True) -> Dict:
+    if len(X) != len(Y):
+        raise ValueError("Size of X and Y do not match")
+
+    # Start the timer
+    start_time = time.time()
+
+    # Calculate the copula grid using Monte Carlo estimation
+    C_grid = calculate_copula_mc_grid(np.array(X), np.array(Y), MC)
+
+    # Create the Q grid
+    Q_grid = create_Q_grid(k_plot_grid)
+
+    # Calculate the Q function on the grid
+    plot_points = Q_function(Q_grid, C_grid)
+
+    # Stop the timer
+    time_taken = time.time() - start_time
+    print(f"Time taken for calculations: {time_taken:.2f} seconds")
+
+    # Create a contour plot using Plotly
+    contour_plot = go.Contour(
+        z=plot_points['z'],
+        x=plot_points['x'],
+        y=plot_points['y'],
+        contours=dict(
+            coloring='heatmap'
+        )
+    )
+
+    layout = go.Layout(
+        title='Q Function Contour Plot',
+        xaxis_title='X',
+        yaxis_title='Y',
+        height=600,
+        width=600
+    )
+
+    fig = go.Figure(data=[contour_plot], layout=layout)
+
+    # Display the plot if the display flag is True
+    if display:
+        fig.show()
 
-# def create_Q_plot(x, y):
-#     return Q_function(create_Q_grid(), calculate_copula_mc_grid(x, y))
+    # Return the plot and data as a dictionary
+    return {
+        'Q_plot': fig,
+        'C_grid': C_grid,
+        'Q_grid': Q_grid,
+        'plot_points': plot_points
+    }