Modernize ripleyk package with tests, CI, and updated docs

.flake8

@@ -0,0 +1,5 @@
+[flake8]
+max-line-length = 88
+extend-ignore = E203, E501
+per-file-ignores =
+    ripleyk/__init__.py:F401,F403

.github/workflows/python-package.yml

@@ -16,7 +16,7 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        python-version: ["3.9", "3.10", "3.11"]
+        python-version: ["3.9", "3.10", "3.11", "3.12", "3.13"]
 
     steps:
     - uses: actions/checkout@v3
@@ -27,14 +27,8 @@ jobs:
     - name: Install dependencies
       run: |
         python -m pip install --upgrade pip
-        python -m pip install flake8 pytest
-        if [ -f requirements.txt ]; then pip install -r requirements.txt; fi
-    - name: Lint with flake8
-      run: |
-        # stop the build if there are Python syntax errors or undefined names
-        flake8 . --count --select=E9,F63,F7,F82 --show-source --statistics
-        # exit-zero treats all errors as warnings. The GitHub editor is 127 chars wide
-        flake8 . --count --exit-zero --max-complexity=10 --max-line-length=127 --statistics
+        pip install -r requirements-dev.txt
+        pip install -e .
     - name: Test with pytest
       run: |
-        pytest
+        pytest

.pre-commit-config.yaml

@@ -0,0 +1,15 @@
+repos:
+-   repo: https://github.com/pre-commit/pre-commit-hooks
+    rev: v6.0.0
+    hooks:
+    -   id: check-yaml
+    -   id: end-of-file-fixer
+    -   id: trailing-whitespace
+-   repo: https://github.com/psf/black
+    rev: 25.9.0
+    hooks:
+    -   id: black
+-   repo: https://github.com/pycqa/flake8
+    rev: 7.3.0
+    hooks:
+    -   id: flake8

README.md

@@ -1,5 +1,8 @@
 # RipleyK
-Calculation of the Ripley K ([spatial statistics](https://en.wikipedia.org/wiki/Spatial_descriptive_statistics)) value in python. This project is still being developed and currently only supports 'circle' based bounding regions for automated boundary correction. Can support 'rectangle' based bounding regions if you do not require boundary corrections. This package allows quick calculation (using [kd-trees](https://en.wikipedia.org/wiki/K-d_tree)) the RipleyK values for 1D-3D systems. 
+
+[![CI](https://github.com/SamPIngram/RipleyK/actions/workflows/ci.yml/badge.svg)](https://github.com/SamPIngram/RipleyK/actions/workflows/ci.yml)
+
+Calculation of the Ripley K ([spatial statistics](https://en.wikipedia.org/wiki/Spatial_descriptive_statistics)) value in python. This project is still being developed and currently only supports 'circle' based bounding regions for automated boundary correction. Can support 'rectangle' based bounding regions if you do not require boundary corrections. This package allows quick calculation (using [kd-trees](https://en.wikipedia.org/wiki/K-d_tree)) the RipleyK values for 1D-3D systems.
 
 ## Installation
 You can install the RipleyK package using the following pip command:
@@ -9,11 +12,11 @@ pip install ripleyk
 
 To get started from source quickly follow these steps:
 
-1. Clone or download this repository and launch python within the folder.
+1. Clone or download this repository.
 
-2. Make a python environment to run in. Always recommend the use of virtual environments for safety. Must be Python3 and currently only tested with Python 3.8.
+2. Make a python environment to run in. We recommend the use of virtual environments. This package is tested on Python 3.8, 3.9, 3.10, and 3.11.
 
-3. Install requirement.txt file into your new python environment
+3. Install requirements using the `requirements.txt` file into your new python environment:
 ```
 pip install -r requirements.txt
 ```
@@ -28,19 +31,19 @@ The mathematical equations for the calculated Ripley K value and normalised L va
 <br/>
 <img src="https://render.githubusercontent.com/render/math?math=L(r) = D \frac{\sum_{i=1}^{n} \sum_{i\ne j} I[D(i,j)\leq r]}{\omega n^{2}} - 2r">
 
-### 2D Equations: 
+### 2D Equations:
 
 <img src="https://render.githubusercontent.com/render/math?math=K(r) = A \frac{\sum_{i=1}^{n} \sum_{i\ne j} I[D(i,j)\leq r]}{\omega n^{2}}">
 
 <br/>
 <img src="https://render.githubusercontent.com/render/math?math=L(r) = A \frac{\sum_{i=1}^{n} \sum_{i\ne j} I[D(i,j)\leq r]}{\omega n^{2}} - \pi r^{2}">
 
-### 3D Equations: 
+### 3D Equations:
 
 <img src="https://render.githubusercontent.com/render/math?math=K(r) = V \frac{\sum_{i=1}^{n} \sum_{i\ne j} I[D(i,j)\leq r]}{\omega n^{2}}">
 <br/>
 
-<img src="https://render.githubusercontent.com/render/math?math=L(r) = V \frac{\sum_{i=1}^{n} \sum_{i\ne j} I[D(i,j)\leq r]}{\omega n^{2}} - \frac{4}{3} \pi r^{3}">
+<img src="httpshttps://render.githubusercontent.com/render/math?math=L(r) = V \frac{\sum_{i=1}^{n} \sum_{i\ne j} I[D(i,j)\leq r]}{\omega n^{2}} - \frac{4}{3} \pi r^{3}">
 
 Note the term "region" is being used interchangeably for distance, area and volume for the 1D, 2D and 3D descriptions respectively. r is the line (1D) or radius (2D/3D) of the search region. ω is the ratio of overlap of the search region with the whole sample region, this is 1 if the search region is entirely within the sample region and <1 if some of the search region is outside of the sample region. If all of the search region is outside of the sample region this is 0 and the Ripley value will not be calculated. I is the indicator function which will be either 1 if the condition D(i,j)≤r is true or will be 0, where D(i,j) is the euclidean distance between points i and j. The size of the sample region is defined as a distance (D), area (A) or volume (V) within the 1D, 2D and 3D equations respectively.
 
@@ -85,7 +88,7 @@ k = ripleyk.calculate_ripley(radius, bounding_radius, d1=xs, d2=ys)
 print(k)
 ```
 
-You should get a value of ```0.7572``` as k.
+You should get a value of ```0.7573``` as k.
 
 By default this will not include any boundary correction. To add this simply add ```boundary_correct=True``` to the above function. This should always increase the value. Now we get a value of ```0.8646``` as k.
 
@@ -114,9 +117,9 @@ k = ripleyk.calculate_ripley(radii, 1, d1=xs, d2=ys, d3=zs,boundary_correct=True
 print(k)
 ```
 
-You should get the following results: 
+You should get the following results:
 ```
-[-9.290137551123609e-06, -0.00010544767161566743, -0.00039142698870800463, -0.0007282757869681578, -0.0013863220751694216, -0.002301632731796177, -0.002895761209242842, -0.004294205083374969, -0.005929855937486295, -0.007915443959695345]
+[-9.29e-06, -0.000105, -0.000391, -0.000728, -0.001386, -0.002302, -0.002896, -0.004294, -0.005930, -0.007915]
 ```
 
 This can be simply plotted against the radii inputted as seen below (would require installation of [matplotlib](https://pypi.org/project/matplotlib/)):
@@ -134,5 +137,17 @@ plt.show()
 
 
 ## Dependencies
-- [Numpy](https://numpy.org/)
-- [Scipy](https://www.scipy.org/)
+- [Numpy](https://numpy.org/) (>=1.21.0)
+- [Scipy](https://www.scipy.org/) (>=1.7.0)
+
+## Contributing
+Contributions are welcome! To contribute, please follow these steps:
+
+1. Fork the repository.
+2. Create a new branch with a descriptive name.
+3. Make your changes and add tests for them.
+4. Run the tests to ensure everything is working correctly:
+   ```
+   pytest
+   ```
+5. Create a pull request.

requirements-dev.txt

@@ -0,0 +1,4 @@
+pytest>=6.2.0
+pre-commit>=2.15.0
+black>=22.3.0
+flake8>=4.0.0

requirements.txt

@@ -1,2 +1,2 @@
-numpy==1.19.2
-scipy==1.5.2
+numpy>=1.21.0
+scipy>=1.7.0

ripleyk/__init__.py

@@ -1 +1,3 @@
-from .ripleyk import *
+from .ripleyk import calculate_ripley
+
+__all__ = ["calculate_ripley"]

ripleyk/ripleyk.py

@@ -1,113 +1,164 @@
 from scipy import spatial
 import numpy as np
-import math
 
 
 def make_tree(d1=None, d2=None, d3=None):
-    active_dimensions = [dimension for dimension in [d1,d2,d3] if dimension is not None]
+    active_dimensions = [
+        dimension for dimension in [d1, d2, d3] if dimension is not None
+    ]
     assert len(active_dimensions) > 0, "Must have at least 1-dimension to make tree"
-    if len(active_dimensions)==1:
+    if len(active_dimensions) == 1:
         points = np.c_[active_dimensions[0].ravel()]
-    elif len(active_dimensions)==2:
+    elif len(active_dimensions) == 2:
         points = np.c_[active_dimensions[0].ravel(), active_dimensions[1].ravel()]
     else:
-        points = np.c_[active_dimensions[0].ravel(), active_dimensions[1].ravel(), active_dimensions[2].ravel()]
+        points = np.c_[
+            active_dimensions[0].ravel(),
+            active_dimensions[1].ravel(),
+            active_dimensions[2].ravel(),
+        ]
     return spatial.cKDTree(points), len(active_dimensions)
 
-def calculate_overlap(point_of_interest, bounding_size, score_radius, sample_shape, dimensions):
-    if dimensions==1:
-            d = point_of_interest[0][0]
-            if d <= abs(score_radius-bounding_size):
-                vol = score_radius
-            elif d >= score_radius+bounding_size:
-                vol = 0
-            else:
-                vol = max(0, min(bounding_size, d+score_radius)-max(-bounding_size, d-score_radius)) 
-    elif sample_shape=='circle':
-        if dimensions==2:
+
+def calculate_overlap(
+    point_of_interest, bounding_size, score_radius, sample_shape, dimensions
+):
+    if dimensions == 1:
+        d = point_of_interest[0][0]
+        if d <= abs(score_radius - bounding_size):
+            vol = score_radius
+        elif d >= score_radius + bounding_size:
+            vol = 0
+        else:
+            vol = max(
+                0,
+                min(bounding_size, d + score_radius)
+                - max(-bounding_size, d - score_radius),
+            )
+    elif sample_shape == "circle":
+        if dimensions == 2:
             d = np.sum((np.array(point_of_interest)) ** 2) ** 0.5
-            if d <= abs(score_radius-bounding_size):
-                vol =  np.pi * score_radius**2
-            elif d >= score_radius+bounding_size:
+            if d <= abs(score_radius - bounding_size):
+                vol = np.pi * score_radius**2
+            elif d >= score_radius + bounding_size:
                 vol = 0
             else:
                 r2, R2, d2 = score_radius**2, bounding_size**2, d**2
                 r, R = score_radius, bounding_size
-                alpha = np.arccos((d2 + r2 - R2) / (2*d*r))
-                beta = np.arccos((d2 + R2 - r2) / (2*d*R))
-                vol = ( r2 * alpha + R2 * beta -
-                        0.5 * (r2 * np.sin(2*alpha) + R2 * np.sin(2*beta))
-                    ) 
-        elif dimensions==3:
+                alpha = np.arccos((d2 + r2 - R2) / (2 * d * r))
+                beta = np.arccos((d2 + R2 - r2) / (2 * d * R))
+                vol = (
+                    r2 * alpha
+                    + R2 * beta
+                    - 0.5 * (r2 * np.sin(2 * alpha) + R2 * np.sin(2 * beta))
+                )
+        elif dimensions == 3:
             d = np.sum((np.array(point_of_interest)) ** 2) ** 0.5
             if d >= score_radius + bounding_size:
-                vol= 0
-            elif bounding_size > (d + score_radius):  # scoring radius is entirely contained in nucleus
-                vol= (4/3) * np.pi * score_radius**3
+                vol = 0
+            elif bounding_size > (
+                d + score_radius
+            ):  # scoring radius is entirely contained in nucleus
+                vol = (4 / 3) * np.pi * score_radius**3
             else:
-                vol = (np.pi * (score_radius + bounding_size - d) ** 2 * (d ** 2 + (2 * d * score_radius - 3 * score_radius ** 2 +
-                                                                2 * d * bounding_size - 3 * bounding_size ** 2)
-                                                        + 6 * score_radius * bounding_size)) / (12 * d)
-    elif sample_shape=='rectangle':
-        if dimensions==2:
+                vol = (
+                    np.pi
+                    * (score_radius + bounding_size - d) ** 2
+                    * (
+                        d**2
+                        + (
+                            2 * d * score_radius
+                            - 3 * score_radius**2
+                            + 2 * d * bounding_size
+                            - 3 * bounding_size**2
+                        )
+                        + 6 * score_radius * bounding_size
+                    )
+                ) / (12 * d)
+    elif sample_shape == "rectangle":
+        if dimensions == 2:
             pass
-        elif dimensions==3:
+        elif dimensions == 3:
             pass
-
-    assert vol > 0, "Attempted to boundary correct a point not within the sample. Check sample_size and points."
+
+    assert vol > 0, (
+        "Attempted to boundary correct a point not within the sample. "
+        "Check sample_size and points."
+    )
     return vol
 
-def calculate_ripley(radii, sample_size, d1=None, d2=None, d3=None, sample_shape='circle', boundary_correct=False, CSR_Normalise=False):
+
+def calculate_ripley(
+    radii,
+    sample_size,
+    d1=None,
+    d2=None,
+    d3=None,
+    sample_shape="circle",
+    boundary_correct=False,
+    CSR_Normalise=False,
+):
     results = []
     tree, dimensions = make_tree(d1=d1, d2=d2, d3=d3)
     if type(radii) is not list:
         radii = [radii]
     for radius in radii:
         if dimensions == 1:
-            score_vol = radius*2
+            score_vol = radius * 2
             bound_size = sample_size
             counts = 0
             for x in zip(d1):
                 if boundary_correct:
-                    vol = calculate_overlap([x], sample_size, radius, sample_shape, dimensions)
-                    boundary_correction = vol/score_vol
-                    counts += (len(tree.query_ball_point([x], radius))-1)/boundary_correction
+                    vol = calculate_overlap(
+                        [x], sample_size, radius, sample_shape, dimensions
+                    )
+                    boundary_correction = vol / score_vol
+                    counts += (
+                        len(tree.query_ball_point([x], radius)) - 1
+                    ) / boundary_correction
                 else:
-                    counts += len(tree.query_ball_point([x], radius))-1
+                    counts += len(tree.query_ball_point([x], radius)) - 1
         elif dimensions == 2:
             score_vol = np.pi * radius**2
-            if sample_shape=='circle':
+            if sample_shape == "circle":
                 bound_size = np.pi * sample_size**2
-            elif sample_shape=='rectangle':
-                bound_size = sample_size[0]*sample_size[1]
+            elif sample_shape == "rectangle":
+                bound_size = sample_size[0] * sample_size[1]
             counts = 0
             for x, y in zip(d1, d2):
                 if boundary_correct:
-                    vol = calculate_overlap([x,y], sample_size, radius, sample_shape, dimensions)
-                    boundary_correction = vol/score_vol
-                    counts += (len(tree.query_ball_point([x,y], radius))-1)/boundary_correction
+                    vol = calculate_overlap(
+                        [x, y], sample_size, radius, sample_shape, dimensions
+                    )
+                    boundary_correction = vol / score_vol
+                    counts += (
+                        len(tree.query_ball_point([x, y], radius)) - 1
+                    ) / boundary_correction
                 else:
-                    counts += len(tree.query_ball_point([x,y], radius))-1
+                    counts += len(tree.query_ball_point([x, y], radius)) - 1
         else:
-            score_vol = (4/3) * np.pi * radius**3
-            if sample_shape=='circle':
-                bound_size = (4/3) * np.pi * sample_size**3
-            elif sample_shape=='rectangle':
-                bound_size = sample_size[0]*sample_size[1]*sample_size[2]
+            score_vol = (4 / 3) * np.pi * radius**3
+            if sample_shape == "circle":
+                bound_size = (4 / 3) * np.pi * sample_size**3
+            elif sample_shape == "rectangle":
+                bound_size = sample_size[0] * sample_size[1] * sample_size[2]
             counts = 0
             for x, y, z in zip(d1, d2, d3):
                 if boundary_correct:
-                    vol = calculate_overlap([x,y,z], sample_size, radius, sample_shape, dimensions)
-                    boundary_correction = vol/score_vol
-                    counts += (len(tree.query_ball_point([x,y,z], radius))-1)/boundary_correction
+                    vol = calculate_overlap(
+                        [x, y, z], sample_size, radius, sample_shape, dimensions
+                    )
+                    boundary_correction = vol / score_vol
+                    counts += (
+                        len(tree.query_ball_point([x, y, z], radius)) - 1
+                    ) / boundary_correction
                 else:
-                    counts += len(tree.query_ball_point([x,y,z], radius))-1
+                    counts += len(tree.query_ball_point([x, y, z], radius)) - 1
+        k_value = bound_size * counts / len(d1) ** 2
         if CSR_Normalise:
-            results.append((bound_size*counts/len(d1)**2) - score_vol)
-        else:
-            results.append(bound_size*counts/len(d1)**2)
-    if len(results)==1:
+            k_value -= score_vol
+        results.append(k_value)
+    if len(results) == 1:
         return results[0]
     else:
         return results
-