diff --git a/.flake8 b/.flake8
new file mode 100644
index 0000000..4195c24
--- /dev/null
+++ b/.flake8
@@ -0,0 +1,5 @@
+[flake8]
+max-line-length = 88
+extend-ignore = E203, E501
+per-file-ignores =
+ ripleyk/__init__.py:F401,F403
diff --git a/.github/workflows/python-package.yml b/.github/workflows/python-package.yml
index 14a4e65..35a4e53 100644
--- a/.github/workflows/python-package.yml
+++ b/.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
\ No newline at end of file
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
new file mode 100644
index 0000000..39124f2
--- /dev/null
+++ b/.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
diff --git a/README.md b/README.md
index d1eb5e7..012d328 100644
--- a/README.md
+++ b/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.
+
+[](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,14 +31,14 @@ The mathematical equations for the calculated Ripley K value and normalised L va
 = D \frac{\sum_{i=1}^{n} \sum_{i\ne j} I[D(i,j)\leq r]}{\omega n^{2}} - 2r)
-### 2D Equations:
+### 2D Equations:
 = A \frac{\sum_{i=1}^{n} \sum_{i\ne j} I[D(i,j)\leq r]}{\omega n^{2}})
 = 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:
 = V \frac{\sum_{i=1}^{n} \sum_{i\ne j} I[D(i,j)\leq r]}{\omega n^{2}})
@@ -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.
diff --git a/requirements-dev.txt b/requirements-dev.txt
new file mode 100644
index 0000000..f8931a2
--- /dev/null
+++ b/requirements-dev.txt
@@ -0,0 +1,4 @@
+pytest>=6.2.0
+pre-commit>=2.15.0
+black>=22.3.0
+flake8>=4.0.0
diff --git a/requirements.txt b/requirements.txt
index 028bbb7..62a0b67 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,2 +1,2 @@
-numpy==1.19.2
-scipy==1.5.2
+numpy>=1.21.0
+scipy>=1.7.0
diff --git a/ripleyk/__init__.py b/ripleyk/__init__.py
index 5d0aed8..e25401a 100644
--- a/ripleyk/__init__.py
+++ b/ripleyk/__init__.py
@@ -1 +1,3 @@
-from .ripleyk import *
\ No newline at end of file
+from .ripleyk import calculate_ripley
+
+__all__ = ["calculate_ripley"]
diff --git a/ripleyk/ripleyk.py b/ripleyk/ripleyk.py
index b6996c6..2c15d58 100644
--- a/ripleyk/ripleyk.py
+++ b/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
-
diff --git a/setup.py b/setup.py
index b0a83fe..35e172f 100644
--- a/setup.py
+++ b/setup.py
@@ -8,7 +8,7 @@
version="0.0.3",
author="Sam Ingram",
author_email="sp_ingram12@yahoo.co.uk",
- description="Calculation of the Ripley K (spatial statistics) value in python",
+ description=("Calculation of the Ripley K (spatial statistics) value in python"),
long_description=long_description,
long_description_content_type="text/markdown",
url="https://github.com/SamPIngram/RipleyK",
@@ -19,8 +19,8 @@
"Operating System :: OS Independent",
],
install_requires=[
- 'numpy>=1.19.2',
- 'scipy>=1.5.2',
+ "numpy>=1.19.2",
+ "scipy>=1.5.2",
],
- python_requires='>=3.5',
+ python_requires=">=3.5",
)
diff --git a/tests/test_ripleyk.py b/tests/test_ripleyk.py
new file mode 100644
index 0000000..68c3c20
--- /dev/null
+++ b/tests/test_ripleyk.py
@@ -0,0 +1,94 @@
+import numpy as np
+import pytest
+from ripleyk import ripleyk
+import random
+
+
+@pytest.fixture(scope="module")
+def random_points_3d():
+ """
+ Generates a set of random 3D points within a sphere of radius 1,
+ based on the data generation code provided in the README.md.
+ `random.seed(0)` is used for reproducibility.
+ """
+ xs = []
+ ys = []
+ zs = []
+ radius = 1
+ random.seed(0)
+ # The loop runs 10,000 times, but the number of points generated
+ # will be lower because some points fall outside the sphere.
+ for _ in range(10000):
+ positioned = False
+ while not positioned:
+ x = random.uniform(-radius, radius)
+ y = random.uniform(-radius, radius)
+ z = random.uniform(-radius, radius)
+ if (x**2) + (y**2) + (z**2) < radius**2:
+ xs.append(x)
+ ys.append(y)
+ zs.append(z)
+ positioned = True
+ return np.array(xs), np.array(ys), np.array(zs)
+
+
+def test_calculate_ripley_2d_single_radius(random_points_3d):
+ xs, ys, _ = random_points_3d
+ radius = 0.5
+ bounding_radius = 1
+ k = ripleyk.calculate_ripley(radius, bounding_radius, d1=xs, d2=ys)
+ assert np.isclose(k, 0.757289, atol=1e-6)
+
+
+def test_calculate_ripley_2d_single_radius_boundary_correct(random_points_3d):
+ xs, ys, _ = random_points_3d
+ radius = 0.5
+ bounding_radius = 1
+ k = ripleyk.calculate_ripley(
+ radius, bounding_radius, d1=xs, d2=ys, boundary_correct=True
+ )
+ assert np.isclose(k, 0.864648, atol=1e-6)
+
+
+def test_calculate_ripley_2d_single_radius_boundary_correct_csr_normalise(
+ random_points_3d,
+):
+ xs, ys, _ = random_points_3d
+ radius = 0.5
+ bounding_radius = 1
+ k = ripleyk.calculate_ripley(
+ radius,
+ bounding_radius,
+ d1=xs,
+ d2=ys,
+ boundary_correct=True,
+ CSR_Normalise=True,
+ )
+ assert np.isclose(k, 0.079249, atol=1e-6)
+
+
+def test_calculate_ripley_3d_multiple_radii(random_points_3d):
+ xs, ys, zs = random_points_3d
+ radii = [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]
+ expected_k = [
+ -9.290137551123609e-06,
+ -0.00010544767161566743,
+ -0.00039142698870800463,
+ -0.0007282757869681578,
+ -0.0013863220751694216,
+ -0.002301632731796177,
+ -0.002895761209242842,
+ -0.004294205083374969,
+ -0.005929855937486295,
+ -0.007915443959695345,
+ ]
+ k = ripleyk.calculate_ripley(
+ radii,
+ 1,
+ d1=xs,
+ d2=ys,
+ d3=zs,
+ boundary_correct=True,
+ CSR_Normalise=True,
+ )
+ assert np.allclose(k, expected_k, atol=1e-6)