diff --git a/.codecov.yml b/.codecov.yml
new file mode 100644
index 0000000..5946969
--- /dev/null
+++ b/.codecov.yml
@@ -0,0 +1,24 @@
+# Codecov configuration for fa2_modified
+
+coverage:
+  status:
+    project:
+      default:
+        target: auto
+        threshold: 5%
+    patch:
+      default:
+        target: auto
+        threshold: 5%
+
+ignore:
+  - "tests/"
+  - "examples/"
+  - "setup.py"
+  - "fa2_modified/fa2util.c"  # Generated C file
+
+comment:
+  layout: "reach, diff, flags, files"
+  behavior: default
+  require_changes: false
+
diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
new file mode 100644
index 0000000..2deb86e
--- /dev/null
+++ b/.github/workflows/ci.yml
@@ -0,0 +1,91 @@
+name: Continuous Integration
+
+on:
+  push:
+    branches: [ main, master ]
+  pull_request:
+    branches: [ main, master, dev ]
+  schedule:
+    # Run tests weekly on Monday at 00:00 UTC
+    - cron: '0 0 * * 1'
+  workflow_dispatch:
+
+jobs:
+  quick-test:
+    name: Quick test (Python 3.11, Ubuntu)
+    runs-on: ubuntu-latest
+    
+    steps:
+    - uses: actions/checkout@v4
+    
+    - name: Set up Python 3.11
+      uses: actions/setup-python@v5
+      with:
+        python-version: '3.11'
+        cache: 'pip'
+    
+    - name: Install dependencies
+      run: |
+        python -m pip install --upgrade pip setuptools wheel
+        pip install Cython
+        pip install -e .
+        pip install pytest networkx
+    
+    - name: Run quick tests
+      run: |
+        pytest tests/test_fa2util.py tests/test_forceatlas2.py -v
+
+  full-test:
+    name: Full test suite
+    needs: quick-test
+    runs-on: ${{ matrix.os }}
+    strategy:
+      fail-fast: false
+      matrix:
+        os: [ubuntu-latest, macos-latest, windows-latest]
+        python-version: ['3.8', '3.10', '3.11', '3.12', '3.13']
+        exclude:
+          - os: macos-latest
+            python-version: '3.8'
+          - os: windows-latest
+            python-version: '3.8'
+    
+    steps:
+    - uses: actions/checkout@v4
+    
+    - name: Set up Python ${{ matrix.python-version }}
+      uses: actions/setup-python@v5
+      with:
+        python-version: ${{ matrix.python-version }}
+        cache: 'pip'
+    
+    - name: Install system dependencies (Ubuntu)
+      if: runner.os == 'Linux'
+      run: |
+        sudo apt-get update
+        sudo apt-get install -y build-essential
+    
+    - name: Install dependencies
+      run: |
+        python -m pip install --upgrade pip setuptools wheel
+        pip install Cython
+        pip install -e .
+        pip install pytest pytest-cov networkx
+    
+    - name: Install igraph (non-Windows)
+      if: runner.os != 'Windows'
+      continue-on-error: true
+      run: pip install python-igraph
+    
+    - name: Run tests
+      run: |
+        pytest tests/ -v --cov=fa2_modified --cov-report=xml
+    
+    - name: Upload coverage
+      if: matrix.os == 'ubuntu-latest' && matrix.python-version == '3.11'
+      uses: codecov/codecov-action@v4
+      with:
+        file: ./coverage.xml
+        flags: unittests
+        fail_ci_if_error: false
+
diff --git a/.github/workflows/code-quality.yml b/.github/workflows/code-quality.yml
new file mode 100644
index 0000000..9065258
--- /dev/null
+++ b/.github/workflows/code-quality.yml
@@ -0,0 +1,87 @@
+name: Code Quality
+
+on:
+  pull_request:
+    branches: [ main, master, dev ]
+  push:
+    branches: [ main, master ]
+
+jobs:
+  lint:
+    name: Linting (flake8)
+    runs-on: ubuntu-latest
+    
+    steps:
+    - uses: actions/checkout@v4
+    
+    - name: Set up Python
+      uses: actions/setup-python@v5
+      with:
+        python-version: '3.11'
+    
+    - name: Install flake8
+      run: pip install flake8
+    
+    - name: Lint with flake8
+      run: |
+        # Stop the build if there are Python syntax errors or undefined names
+        flake8 fa2_modified --count --select=E9,F63,F7,F82 --show-source --statistics
+        # Exit-zero treats all errors as warnings
+        flake8 fa2_modified --count --exit-zero --max-complexity=10 --max-line-length=127 --statistics
+
+  format-check:
+    name: Code Formatting (black)
+    runs-on: ubuntu-latest
+    
+    steps:
+    - uses: actions/checkout@v4
+    
+    - name: Set up Python
+      uses: actions/setup-python@v5
+      with:
+        python-version: '3.11'
+    
+    - name: Install black
+      run: pip install black
+    
+    - name: Check formatting
+      run: black --check --diff fa2_modified/ tests/
+
+  import-sort:
+    name: Import Sorting (isort)
+    runs-on: ubuntu-latest
+    
+    steps:
+    - uses: actions/checkout@v4
+    
+    - name: Set up Python
+      uses: actions/setup-python@v5
+      with:
+        python-version: '3.11'
+    
+    - name: Install isort
+      run: pip install isort
+    
+    - name: Check import sorting
+      run: isort --check-only --diff fa2_modified/ tests/
+
+  type-check:
+    name: Type Checking (mypy)
+    runs-on: ubuntu-latest
+    continue-on-error: true  # Don't fail the build on type errors for now
+    
+    steps:
+    - uses: actions/checkout@v4
+    
+    - name: Set up Python
+      uses: actions/setup-python@v5
+      with:
+        python-version: '3.11'
+    
+    - name: Install dependencies
+      run: |
+        pip install mypy numpy scipy tqdm
+    
+    - name: Type check with mypy
+      run: mypy fa2_modified/ --ignore-missing-imports || echo "Type checking found issues (not blocking)"
+
diff --git a/.github/workflows/python-publish.yml b/.github/workflows/python-publish.yml
index f80697e..33ade86 100644
--- a/.github/workflows/python-publish.yml
+++ b/.github/workflows/python-publish.yml
@@ -1,45 +1,101 @@
-# This workflow will upload a Python Package using Twine when a release is created
-# For more information see: https://docs.github.com/en/actions/automating-builds-and-tests/building-and-testing-python#publishing-to-package-registries
-
-# This workflow uses actions that are not certified by GitHub.
-# They are provided by a third-party and are governed by
-# separate terms of service, privacy policy, and support
-# documentation.
-
-name: Upload Python Package
+name: Build and Publish Wheels
 
 on:
   release:
     types: [published]
+  workflow_dispatch:
 
 permissions:
   contents: read
 
 jobs:
-  deploy:
+  build-wheels:
+    name: Build wheels on ${{ matrix.os }}
+    runs-on: ${{ matrix.os }}
+    strategy:
+      fail-fast: false
+      matrix:
+        os: [ubuntu-latest, macos-latest, windows-latest]
+
+    steps:
+      - uses: actions/checkout@v4
+
+      - name: Set up Python
+        uses: actions/setup-python@v5
+        with:
+          python-version: '3.11'
+
+      - name: Set up QEMU for Linux ARM builds
+        if: runner.os == 'Linux'
+        uses: docker/setup-qemu-action@v3
+        with:
+          platforms: arm64
+
+      - name: Install cibuildwheel
+        run: |
+          python -m pip install --upgrade pip
+          pip install cibuildwheel==2.16.5
 
+      - name: Build wheels
+        env:
+          CIBW_BUILD: "cp38-* cp39-* cp310-* cp311-* cp312-* cp313-*"
+          CIBW_SKIP: "pp* *-musllinux_*"
+          CIBW_BEFORE_BUILD: "pip install -U pip setuptools wheel Cython numpy"
+          CIBW_BEFORE_BUILD_LINUX: "pip install -U pip setuptools wheel Cython numpy && yum install -y gcc"
+          CIBW_TEST_COMMAND: "python -c \"import fa2_modified; print('Import successful')\""
+          CIBW_TEST_REQUIRES: "numpy scipy tqdm"
+          CIBW_ARCHS_MACOS: "x86_64 arm64"
+          CIBW_ARCHS_WINDOWS: "AMD64"
+          CIBW_ARCHS_LINUX: "x86_64 aarch64"
+          CIBW_BUILD_VERBOSITY: 1
+        run: |
+          python -m cibuildwheel --output-dir wheelhouse
+
+      - name: Build sdist (source)
+        if: matrix.os == 'ubuntu-latest'
+        run: |
+          python -m pip install build Cython
+          python -m build --sdist --outdir wheelhouse
+
+      - name: List built wheels
+        run: ls -lh wheelhouse/
+
+      - name: Upload artifacts
+        uses: actions/upload-artifact@v4
+        with:
+          name: wheels-${{ matrix.os }}-${{ strategy.job-index }}
+          path: wheelhouse/*
+          if-no-files-found: error
+
+  publish:
+    name: Publish to PyPI
+    needs: build-wheels
     runs-on: ubuntu-latest
+    environment:
+      name: pypi
+      url: https://pypi.org/p/fa2-modified
+    permissions:
+      id-token: write
+      contents: read
 
     steps:
-    - uses: actions/checkout@v3
-    - name: Set up Python
-      uses: docker://quay.io/pypa/manylinux2014_x86_64
-      with:
-        python-version: '3.x'
-    - name: Install dependencies
-      run: |
-        python -m pip install --upgrade pip
-        pip install setuptools wheel twine auditwheel
-    - name: Build package
-      run: |
-        python3 -m pip install --upgrade build
-        python3 -m build
-    - name: Repair wheels
-      run: |
-        auditwheel repair --wheel-dir dist dist/*.whl
-        rm -rf dist/*linux_x86_64*.whl
-    - name: Publish package
-      uses: pypa/gh-action-pypi-publish@27b31702a0e7fc50959f5ad993c78deac1bdfc29
-      with:
-        user: __token__
-        password: ${{ secrets.PYPI_API_TOKEN_GITHUB }}
+      - uses: actions/download-artifact@v4
+        with:
+          path: dist
+          merge-multiple: true
+
+      - name: List all files to publish
+        run: ls -lh dist/
+
+      - name: Verify wheels
+        run: |
+          pip install twine
+          twine check dist/*
+
+      - name: Publish package to PyPI
+        uses: pypa/gh-action-pypi-publish@release/v1
+        with:
+          password: ${{ secrets.PYPI_API_TOKEN_GITHUB }}
+          packages-dir: dist
+          verbose: true
+          print-hash: true
diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml
new file mode 100644
index 0000000..2547862
--- /dev/null
+++ b/.github/workflows/tests.yml
@@ -0,0 +1,167 @@
+name: Tests
+
+on:
+  pull_request:
+    branches: [ main, master, dev ]
+  push:
+    branches: [ main, master ]
+  workflow_dispatch:
+
+jobs:
+  test:
+    name: Test Python ${{ matrix.python-version }} on ${{ matrix.os }}
+    runs-on: ${{ matrix.os }}
+    strategy:
+      fail-fast: false
+      matrix:
+        os: [ubuntu-latest, macos-latest, windows-latest]
+        python-version: ['3.9', '3.10', '3.11', '3.12', '3.13']
+    
+    steps:
+    - name: Checkout code
+      uses: actions/checkout@v4
+    
+    - name: Set up Python ${{ matrix.python-version }}
+      uses: actions/setup-python@v5
+      with:
+        python-version: ${{ matrix.python-version }}
+        cache: 'pip'
+    
+    - name: Install system dependencies (Ubuntu)
+      if: runner.os == 'Linux'
+      run: |
+        sudo apt-get update
+        sudo apt-get install -y build-essential
+    
+    - name: Install system dependencies (macOS)
+      if: runner.os == 'macOS'
+      run: |
+        brew install gcc
+    
+    - name: Upgrade pip and install build tools
+      run: |
+        python -m pip install --upgrade pip setuptools wheel
+    
+    - name: Install Cython
+      run: |
+        pip install Cython
+    
+    - name: Install package with dependencies
+      run: |
+        pip install -e .
+    
+    - name: Install test dependencies
+      run: |
+        pip install pytest pytest-cov
+        pip install networkx
+    
+    - name: Install igraph (non-Windows)
+      if: runner.os != 'Windows'
+      continue-on-error: true
+      run: |
+        pip install python-igraph
+    
+    - name: List installed packages
+      run: |
+        pip list
+    
+    - name: Run tests
+      run: |
+        pytest tests/ -v --tb=short --cov=fa2_modified --cov-report=xml --cov-report=term
+    
+    - name: Upload coverage to Codecov
+      if: matrix.os == 'ubuntu-latest' && matrix.python-version == '3.11'
+      uses: codecov/codecov-action@v4
+      with:
+        file: ./coverage.xml
+        flags: unittests
+        name: codecov-umbrella
+        fail_ci_if_error: false
+
+  test-minimal:
+    name: Test minimal dependencies (Python ${{ matrix.python-version }})
+    runs-on: ubuntu-latest
+    strategy:
+      matrix:
+        python-version: ['3.8', '3.12']
+    
+    steps:
+    - name: Checkout code
+      uses: actions/checkout@v4
+    
+    - name: Set up Python ${{ matrix.python-version }}
+      uses: actions/setup-python@v5
+      with:
+        python-version: ${{ matrix.python-version }}
+        cache: 'pip'
+    
+    - name: Install minimal dependencies
+      run: |
+        python -m pip install --upgrade pip setuptools wheel
+        pip install Cython
+        pip install -e .
+        pip install pytest
+    
+    - name: Run core tests only
+      run: |
+        pytest tests/test_fa2util.py tests/test_forceatlas2.py -v
+
+  test-without-cython:
+    name: Test without Cython (Python 3.11)
+    runs-on: ubuntu-latest
+    
+    steps:
+    - name: Checkout code
+      uses: actions/checkout@v4
+    
+    - name: Set up Python 3.11
+      uses: actions/setup-python@v5
+      with:
+        python-version: '3.11'
+        cache: 'pip'
+    
+    - name: Install without Cython
+      run: |
+        python -m pip install --upgrade pip setuptools wheel
+        # Install dependencies but not Cython
+        pip install numpy scipy tqdm pytest networkx
+        # Try to install the package (should fall back to pure Python or pre-generated C)
+        pip install -e . || echo "Installation failed as expected without Cython"
+    
+    - name: Run basic tests
+      run: |
+        pytest tests/test_fa2util.py -v || echo "Tests may fail without compiled extension"
+
+  lint:
+    name: Linting and code quality
+    runs-on: ubuntu-latest
+    
+    steps:
+    - name: Checkout code
+      uses: actions/checkout@v4
+    
+    - name: Set up Python
+      uses: actions/setup-python@v5
+      with:
+        python-version: '3.11'
+        cache: 'pip'
+    
+    - name: Install linting tools
+      run: |
+        pip install flake8 black isort
+    
+    - name: Check with flake8
+      run: |
+        # Stop the build if there are Python syntax errors or undefined names
+        flake8 fa2_modified --count --select=E9,F63,F7,F82 --show-source --statistics
+        # Exit-zero treats all errors as warnings
+        flake8 fa2_modified --count --exit-zero --max-complexity=10 --max-line-length=127 --statistics
+    
+    - name: Check code formatting with black
+      run: |
+        black --check fa2_modified/ tests/ || echo "Code needs formatting"
+    
+    - name: Check import sorting with isort
+      run: |
+        isort --check-only fa2_modified/ tests/ || echo "Imports need sorting"
+
diff --git a/.gitignore b/.gitignore
index b3d2a18..519d31b 100644
--- a/.gitignore
+++ b/.gitignore
@@ -41,10 +41,12 @@ htmlcov/
 .coverage
 .coverage.*
 .cache
+.pytest_cache/
 nosetests.xml
 coverage.xml
 *.cover
 .hypothesis/
+.pytest/
 
 # Translations
 *.mo
diff --git a/MANIFEST.in b/MANIFEST.in
index e55be1c..fb44dea 100644
--- a/MANIFEST.in
+++ b/MANIFEST.in
@@ -1,2 +1,5 @@
 include fa2_modified/fa2util.pxd
+include fa2_modified/fa2util.c
+include README.md
+include LICENSE
 include examples/*
\ No newline at end of file
diff --git a/README.md b/README.md
index d7849a7..5465b6a 100644
--- a/README.md
+++ b/README.md
@@ -1,9 +1,9 @@
 ## ForceAtlas2 for Python3 (modified/maintained version)
 
-#### This is a maintained version of fa2 [python library](https://github.com/bhargavchippada/forceatlas2). The source code is same as fa2, but this repo can be installed on Python 3.9+ whereas the original fa2 only runs on Python <3.8.
+#### This is a maintained version of fa2 [python library](https://github.com/AminAlam/fa2_modified). The source code is same as fa2, but this repo can be installed on Python 3.9+ whereas the original fa2 only runs on Python <3.8.
 --------
 [![PyPI version](https://img.shields.io/pypi/v/fa2-modified.svg)]([[https://pypi.python.org/pypi/fa2-modified](https://pypi.org/project/fa2-modified)](https://pypi.org/project/fa2-modified))
-[![Python 3.8+](https://img.shields.io/badge/python-3.8+-blue.svg)](https://www.python.org/downloads/)
+[![Python 3.8+](https://img.shields.io/badge/python-3.9+-blue.svg)](https://www.python.org/downloads/)
 [![continuous-integration](https://github.com/irahorecka/pyforceatlas2/workflows/continuous-integration/badge.svg)](https://github.com/irahorecka/fa2-modified/actions)
 [![License: GPL v3](https://img.shields.io/badge/license-GPLv3-purple.svg)](https://raw.githubusercontent.com/aminalam/fa2-modified/main/LICENSE)
 [![DOI](https://img.shields.io/badge/DOI-10.1371/journal.pone.0098679-blue)](https://doi.org/10.1371/journal.pone.0098679)
diff --git a/fa2_modified/fa2util.pxd b/fa2_modified/fa2util.pxd
index ff43384..9345af7 100644
--- a/fa2_modified/fa2util.pxd
+++ b/fa2_modified/fa2util.pxd
@@ -14,6 +14,7 @@
 
 import cython
 
+
 # This will substitute for the nLayout object
 cdef class Node:
     cdef public double mass
@@ -35,26 +36,26 @@ cdef class Edge:
                yDist = cython.double, 
                distance2 = cython.double, 
                factor = cython.double)
-cdef void linRepulsion(Node n1, Node n2, double coefficient=*)
+cpdef void linRepulsion(Node n1, Node n2, double coefficient=*)
 
 @cython.locals(xDist = cython.double,
                yDist = cython.double,
                distance2 = cython.double,
                factor = cython.double)
-cdef void linRepulsion_region(Node n, Region r, double coefficient=*)
+cpdef void linRepulsion_region(Node n, Region r, double coefficient=*)
 
 
 @cython.locals(xDist = cython.double, 
                yDist = cython.double, 
                distance = cython.double, 
                factor = cython.double)
-cdef void linGravity(Node n, double g)
+cpdef void linGravity(Node n, double g)
 
 
 @cython.locals(xDist = cython.double, 
                yDist = cython.double, 
                factor = cython.double)
-cdef void strongGravity(Node n, double g, double coefficient=*)
+cpdef void strongGravity(Node n, double g, double coefficient=*)
 
 @cython.locals(xDist = cython.double, 
                yDist = cython.double, 
@@ -97,7 +98,7 @@ cdef class Region:
 
     @cython.locals(distance = cython.double,
                    subregion = Region)
-    cdef void applyForce(self, Node n, double theta, double coefficient=*)
+    cpdef void applyForce(self, Node n, double theta, double coefficient=*)
 
     @cython.locals(n = Node)
     cpdef applyForceOnNodes(self, list nodes, double theta, double coefficient=*)
diff --git a/fa2_modified/fa2util.py b/fa2_modified/fa2util.py
index a1fadfc..c743b34 100644
--- a/fa2_modified/fa2util.py
+++ b/fa2_modified/fa2util.py
@@ -36,6 +36,7 @@ def __init__(self):
 # Here are some functions from ForceFactory.java
 # =============================================
 
+
 # Repulsion function.  `n1` and `n2` should be nodes.  This will
 # adjust the dx and dy values of `n1`  `n2`
 def linRepulsion(n1, n2, coefficient=0):
@@ -131,18 +132,37 @@ def apply_gravity(nodes, gravity, scalingRatio, useStrongGravity=False):
             strongGravity(n, gravity, scalingRatio)
 
 
-def apply_attraction(nodes, edges, distributedAttraction, coefficient, edgeWeightInfluence):
+def apply_attraction(
+    nodes, edges, distributedAttraction, coefficient, edgeWeightInfluence
+):
     # Optimization, since usually edgeWeightInfluence is 0 or 1, and pow is slow
     if edgeWeightInfluence == 0:
         for edge in edges:
-            linAttraction(nodes[edge.node1], nodes[edge.node2], 1, distributedAttraction, coefficient)
+            linAttraction(
+                nodes[edge.node1],
+                nodes[edge.node2],
+                1,
+                distributedAttraction,
+                coefficient,
+            )
     elif edgeWeightInfluence == 1:
         for edge in edges:
-            linAttraction(nodes[edge.node1], nodes[edge.node2], edge.weight, distributedAttraction, coefficient)
+            linAttraction(
+                nodes[edge.node1],
+                nodes[edge.node2],
+                edge.weight,
+                distributedAttraction,
+                coefficient,
+            )
     else:
         for edge in edges:
-            linAttraction(nodes[edge.node1], nodes[edge.node2], pow(edge.weight, edgeWeightInfluence),
-                          distributedAttraction, coefficient)
+            linAttraction(
+                nodes[edge.node1],
+                nodes[edge.node2],
+                pow(edge.weight, edgeWeightInfluence),
+                distributedAttraction,
+                coefficient,
+            )
 
 
 # For Barnes Hut Optimization
@@ -170,7 +190,9 @@ def updateMassAndGeometry(self):
 
             self.size = 0.0
             for n in self.nodes:
-                distance = sqrt((n.x - self.massCenterX) ** 2 + (n.y - self.massCenterY) ** 2)
+                distance = sqrt(
+                    (n.x - self.massCenterX) ** 2 + (n.y - self.massCenterY) ** 2
+                )
                 self.size = max(self.size, 2 * distance)
 
     def buildSubRegions(self):
@@ -179,8 +201,8 @@ def buildSubRegions(self):
             bottomleftNodes = []
             toprightNodes = []
             bottomrightNodes = []
-            # Optimization: The distribution of self.nodes into 
-            # subregions now requires only one for loop. Removed 
+            # Optimization: The distribution of self.nodes into
+            # subregions now requires only one for loop. Removed
             # topNodes and bottomNodes arrays: memory space saving.
             for n in self.nodes:
                 if n.x < self.massCenterX:
@@ -192,7 +214,7 @@ def buildSubRegions(self):
                     if n.y < self.massCenterY:
                         bottomrightNodes.append(n)
                     else:
-                        toprightNodes.append(n)      
+                        toprightNodes.append(n)
 
             if len(topleftNodes) > 0:
                 if len(topleftNodes) < len(self.nodes):
@@ -237,7 +259,9 @@ def applyForce(self, n, theta, coefficient=0):
         if len(self.nodes) < 2:
             linRepulsion(n, self.nodes[0], coefficient)
         else:
-            distance = sqrt((n.x - self.massCenterX) ** 2 + (n.y - self.massCenterY) ** 2)
+            distance = sqrt(
+                (n.x - self.massCenterX) ** 2 + (n.y - self.massCenterY) ** 2
+            )
             if distance * theta > self.size:
                 linRepulsion_region(n, self, coefficient)
             else:
@@ -255,43 +279,58 @@ def adjustSpeedAndApplyForces(nodes, speed, speedEfficiency, jitterTolerance):
     totalSwinging = 0.0  # How much irregular movement
     totalEffectiveTraction = 0.0  # How much useful movement
     for n in nodes:
-        swinging = sqrt((n.old_dx - n.dx) * (n.old_dx - n.dx) + (n.old_dy - n.dy) * (n.old_dy - n.dy))
+        swinging = sqrt(
+            (n.old_dx - n.dx) * (n.old_dx - n.dx)
+            + (n.old_dy - n.dy) * (n.old_dy - n.dy)
+        )
         totalSwinging += n.mass * swinging
-        totalEffectiveTraction += .5 * n.mass * sqrt(
-            (n.old_dx + n.dx) * (n.old_dx + n.dx) + (n.old_dy + n.dy) * (n.old_dy + n.dy))
+        totalEffectiveTraction += (
+            0.5
+            * n.mass
+            * sqrt(
+                (n.old_dx + n.dx) * (n.old_dx + n.dx)
+                + (n.old_dy + n.dy) * (n.old_dy + n.dy)
+            )
+        )
 
     # Optimize jitter tolerance.  The 'right' jitter tolerance for
     # this network. Bigger networks need more tolerance. Denser
     # networks need less tolerance. Totally empiric.
-    estimatedOptimalJitterTolerance = .05 * sqrt(len(nodes))
+    estimatedOptimalJitterTolerance = 0.05 * sqrt(len(nodes))
     minJT = sqrt(estimatedOptimalJitterTolerance)
     maxJT = 10
-    jt = jitterTolerance * max(minJT,
-                               min(maxJT, estimatedOptimalJitterTolerance * totalEffectiveTraction / (
-                                   len(nodes) * len(nodes))))
+    jt = jitterTolerance * max(
+        minJT,
+        min(
+            maxJT,
+            estimatedOptimalJitterTolerance
+            * totalEffectiveTraction
+            / (len(nodes) * len(nodes)),
+        ),
+    )
 
     minSpeedEfficiency = 0.05
 
     # Protective against erratic behavior
     if totalEffectiveTraction and totalSwinging / totalEffectiveTraction > 2.0:
         if speedEfficiency > minSpeedEfficiency:
-            speedEfficiency *= .5
+            speedEfficiency *= 0.5
         jt = max(jt, jitterTolerance)
 
     if totalSwinging == 0:
-        targetSpeed = float('inf')
+        targetSpeed = float("inf")
     else:
         targetSpeed = jt * speedEfficiency * totalEffectiveTraction / totalSwinging
 
     if totalSwinging > jt * totalEffectiveTraction:
         if speedEfficiency > minSpeedEfficiency:
-            speedEfficiency *= .7
+            speedEfficiency *= 0.7
     elif speed < 1000:
         speedEfficiency *= 1.3
 
     # But the speed shoudn't rise too much too quickly, since it would
     # make the convergence drop dramatically.
-    maxRise = .5
+    maxRise = 0.5
     speed = speed + min(targetSpeed - speed, maxRise * speed)
 
     # Apply forces.
@@ -299,14 +338,17 @@ def adjustSpeedAndApplyForces(nodes, speed, speedEfficiency, jitterTolerance):
     # Need to add a case if adjustSizes ("prevent overlap") is
     # implemented.
     for n in nodes:
-        swinging = n.mass * sqrt((n.old_dx - n.dx) * (n.old_dx - n.dx) + (n.old_dy - n.dy) * (n.old_dy - n.dy))
+        swinging = n.mass * sqrt(
+            (n.old_dx - n.dx) * (n.old_dx - n.dx)
+            + (n.old_dy - n.dy) * (n.old_dy - n.dy)
+        )
         factor = speed / (1.0 + sqrt(speed * swinging))
         n.x = n.x + (n.dx * factor)
         n.y = n.y + (n.dy * factor)
 
     values = {}
-    values['speed'] = speed
-    values['speedEfficiency'] = speedEfficiency
+    values["speed"] = speed
+    values["speedEfficiency"] = speedEfficiency
 
     return values
 
@@ -315,6 +357,10 @@ def adjustSpeedAndApplyForces(nodes, speed, speedEfficiency, jitterTolerance):
     import cython
 
     if not cython.compiled:
-        print("Warning: uncompiled fa2util module.  Compile with cython for a 10-100x speed boost.")
+        print(
+            "Warning: uncompiled fa2util module.  Compile with cython for a 10-100x speed boost."
+        )
 except:
-    print("No cython detected.  Install cython and compile the fa2util module for a 10-100x speed boost.")
+    print(
+        "No cython detected.  Install cython and compile the fa2util module for a 10-100x speed boost."
+    )
diff --git a/fa2_modified/forceatlas2.py b/fa2_modified/forceatlas2.py
index fa7a67e..fd729c2 100644
--- a/fa2_modified/forceatlas2.py
+++ b/fa2_modified/forceatlas2.py
@@ -38,34 +38,35 @@ def start(self):
         self.start_time = time.time()
 
     def stop(self):
-        self.total_time += (time.time() - self.start_time)
+        self.total_time += time.time() - self.start_time
 
     def display(self):
         print(self.name, " took ", "%.2f" % self.total_time, " seconds")
 
 
 class ForceAtlas2:
-    def __init__(self,
-                 # Behavior alternatives
-                 outboundAttractionDistribution=False,  # Dissuade hubs
-                 linLogMode=False,  # NOT IMPLEMENTED
-                 adjustSizes=False,  # Prevent overlap (NOT IMPLEMENTED)
-                 edgeWeightInfluence=1.0,
-
-                 # Performance
-                 jitterTolerance=1.0,  # Tolerance
-                 barnesHutOptimize=True,
-                 barnesHutTheta=1.2,
-                 multiThreaded=False,  # NOT IMPLEMENTED
-
-                 # Tuning
-                 scalingRatio=2.0,
-                 strongGravityMode=False,
-                 gravity=1.0,
-
-                 # Log
-                 verbose=True):
-        assert linLogMode == adjustSizes == multiThreaded == False, "You selected a feature that has not been implemented yet..."
+    def __init__(
+        self,
+        # Behavior alternatives
+        outboundAttractionDistribution=False,  # Dissuade hubs
+        linLogMode=False,  # NOT IMPLEMENTED
+        adjustSizes=False,  # Prevent overlap (NOT IMPLEMENTED)
+        edgeWeightInfluence=1.0,
+        # Performance
+        jitterTolerance=1.0,  # Tolerance
+        barnesHutOptimize=True,
+        barnesHutTheta=1.2,
+        multiThreaded=False,  # NOT IMPLEMENTED
+        # Tuning
+        scalingRatio=2.0,
+        strongGravityMode=False,
+        gravity=1.0,
+        # Log
+        verbose=True,
+    ):
+        assert (
+            linLogMode == adjustSizes == multiThreaded == False
+        ), "You selected a feature that has not been implemented yet..."
         self.outboundAttractionDistribution = outboundAttractionDistribution
         self.linLogMode = linLogMode
         self.adjustSizes = adjustSizes
@@ -78,20 +79,27 @@ def __init__(self,
         self.gravity = gravity
         self.verbose = verbose
 
-    def init(self,
-             G,  # a graph in 2D numpy ndarray format (or) scipy sparse matrix format
-             pos=None  # Array of initial positions
-             ):
+    def init(
+        self,
+        G,  # a graph in 2D numpy ndarray format (or) scipy sparse matrix format
+        pos=None,  # Array of initial positions
+    ):
         isSparse = False
         if isinstance(G, numpy.ndarray):
             # Check our assumptions
             assert G.shape == (G.shape[0], G.shape[0]), "G is not 2D square"
-            assert numpy.all(G.T == G), "G is not symmetric.  Currently only undirected graphs are supported"
-            assert isinstance(pos, numpy.ndarray) or (pos is None), "Invalid node positions"
+            assert numpy.all(
+                G.T == G
+            ), "G is not symmetric.  Currently only undirected graphs are supported"
+            assert isinstance(pos, numpy.ndarray) or (
+                pos is None
+            ), "Invalid node positions"
         elif scipy.sparse.issparse(G):
             # Check our assumptions
             assert G.shape == (G.shape[0], G.shape[0]), "G is not 2D square"
-            assert isinstance(pos, numpy.ndarray) or (pos is None), "Invalid node positions"
+            assert isinstance(pos, numpy.ndarray) or (
+                pos is None
+            ), "Invalid node positions"
             G = G.tolil()
             isSparse = True
         else:
@@ -121,7 +129,8 @@ def init(self,
         edges = []
         es = numpy.asarray(G.nonzero()).T
         for e in es:  # Iterate through edges
-            if e[1] <= e[0]: continue  # Avoid duplicate edges
+            if e[1] <= e[0]:
+                continue  # Avoid duplicate edges
             edge = fa2util.Edge()
             edge.node1 = e[0]  # The index of the first node in `nodes`
             edge.node2 = e[1]  # The index of the second node in `nodes`
@@ -146,11 +155,12 @@ def init(self,
     #
     # Currently, only undirected graphs are supported so the adjacency matrix
     # should be symmetric.
-    def forceatlas2(self,
-                    G,  # a graph in 2D numpy ndarray format (or) scipy sparse matrix format
-                    pos=None,  # Array of initial positions
-                    iterations=100  # Number of times to iterate the main loop
-                    ):
+    def forceatlas2(
+        self,
+        G,  # a graph in 2D numpy ndarray format (or) scipy sparse matrix format
+        pos=None,  # Array of initial positions
+        iterations=100,  # Number of times to iterate the main loop
+    ):
         # Initializing, initAlgo()
         # ================================================================
 
@@ -196,27 +206,41 @@ def forceatlas2(self,
             repulsion_timer.start()
             # parallelization should be implemented here
             if self.barnesHutOptimize:
-                rootRegion.applyForceOnNodes(nodes, self.barnesHutTheta, self.scalingRatio)
+                rootRegion.applyForceOnNodes(
+                    nodes, self.barnesHutTheta, self.scalingRatio
+                )
             else:
                 fa2util.apply_repulsion(nodes, self.scalingRatio)
             repulsion_timer.stop()
 
             # Gravitational forces
             gravity_timer.start()
-            fa2util.apply_gravity(nodes, self.gravity, scalingRatio=self.scalingRatio, useStrongGravity=self.strongGravityMode)
+            fa2util.apply_gravity(
+                nodes,
+                self.gravity,
+                scalingRatio=self.scalingRatio,
+                useStrongGravity=self.strongGravityMode,
+            )
             gravity_timer.stop()
 
             # If other forms of attraction were implemented they would be selected here.
             attraction_timer.start()
-            fa2util.apply_attraction(nodes, edges, self.outboundAttractionDistribution, outboundAttCompensation,
-                                     self.edgeWeightInfluence)
+            fa2util.apply_attraction(
+                nodes,
+                edges,
+                self.outboundAttractionDistribution,
+                outboundAttCompensation,
+                self.edgeWeightInfluence,
+            )
             attraction_timer.stop()
 
             # Adjust speeds and apply forces
             applyforces_timer.start()
-            values = fa2util.adjustSpeedAndApplyForces(nodes, speed, speedEfficiency, self.jitterTolerance)
-            speed = values['speed']
-            speedEfficiency = values['speedEfficiency']
+            values = fa2util.adjustSpeedAndApplyForces(
+                nodes, speed, speedEfficiency, self.jitterTolerance
+            )
+            speed = values["speed"]
+            speedEfficiency = values["speedEfficiency"]
             applyforces_timer.stop()
 
         if self.verbose:
@@ -233,19 +257,25 @@ def forceatlas2(self,
     #
     # This function returns a NetworkX layout, which is really just a
     # dictionary of node positions (2D X-Y tuples) indexed by the node name.
-    def forceatlas2_networkx_layout(self, G, pos=None, iterations=100, weight_attr=None):
+    def forceatlas2_networkx_layout(
+        self, G, pos=None, iterations=100, weight_attr=None
+    ):
         import networkx
+
         try:
             import cynetworkx
         except ImportError:
             cynetworkx = None
 
-        assert (
-            isinstance(G, networkx.classes.graph.Graph)
-            or (cynetworkx and isinstance(G, cynetworkx.classes.graph.Graph))
+        assert isinstance(G, networkx.classes.graph.Graph) or (
+            cynetworkx and isinstance(G, cynetworkx.classes.graph.Graph)
         ), "Not a networkx graph"
-        assert isinstance(pos, dict) or (pos is None), "pos must be specified as a dictionary, as in networkx"
-        M = networkx.to_scipy_sparse_array(G, dtype='f', format='lil', weight=weight_attr)
+        assert isinstance(pos, dict) or (
+            pos is None
+        ), "pos must be specified as a dictionary, as in networkx"
+        M = networkx.to_scipy_sparse_array(
+            G, dtype="f", format="lil", weight=weight_attr
+        )
         if pos is None:
             l = self.forceatlas2(M, pos=None, iterations=iterations)
         else:
@@ -258,8 +288,8 @@ def forceatlas2_networkx_layout(self, G, pos=None, iterations=100, weight_attr=N
     # This function returns an igraph layout
     def forceatlas2_igraph_layout(self, G, pos=None, iterations=100, weight_attr=None):
 
-        from scipy.sparse import csr_matrix
         import igraph
+        from scipy.sparse import csr_matrix
 
         def to_sparse(graph, weight_attr=None):
             edges = graph.get_edgelist()
@@ -272,10 +302,17 @@ def to_sparse(graph, weight_attr=None):
                 edges.extend([(v, u) for u, v in edges])
                 weights.extend(weights)
 
-            return csr_matrix((weights, list(zip(*edges))))
+            # Handle empty edges case
+            n = graph.vcount()
+            if len(edges) == 0:
+                return csr_matrix((n, n))
+
+            return csr_matrix((weights, list(zip(*edges))), shape=(n, n))
 
         assert isinstance(G, igraph.Graph), "Not a igraph graph"
-        assert isinstance(pos, (list, numpy.ndarray)) or (pos is None), "pos must be a list or numpy array"
+        assert isinstance(pos, (list, numpy.ndarray)) or (
+            pos is None
+        ), "pos must be a list or numpy array"
 
         if isinstance(pos, list):
             pos = numpy.array(pos)
diff --git a/pyproject.toml b/pyproject.toml
index eb56957..fda8443 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -1,2 +1,3 @@
 [build-system]
-requires = ["setuptools", "wheel", "Cython"]
\ No newline at end of file
+requires = ["setuptools>=45", "wheel", "Cython>=0.29.0", "numpy>=1.19.0"]
+build-backend = "setuptools.build_meta"
\ No newline at end of file
diff --git a/pytest.ini b/pytest.ini
new file mode 100644
index 0000000..c9ffd92
--- /dev/null
+++ b/pytest.ini
@@ -0,0 +1,32 @@
+[pytest]
+# pytest configuration for fa2_modified
+
+# Test discovery patterns
+python_files = test_*.py
+python_classes = Test*
+python_functions = test_*
+
+# Test paths
+testpaths = tests
+
+# Output options
+addopts =
+    -v
+    --strict-markers
+    --tb=short
+    --disable-warnings
+
+# Markers
+markers =
+    slow: marks tests as slow (deselect with '-m "not slow"')
+    integration: marks tests as integration tests
+    requires_networkx: requires networkx to be installed
+    requires_igraph: requires igraph to be installed
+
+# Coverage options (if using pytest-cov)
+# Uncomment to enable coverage reporting
+# addopts = --cov=fa2_modified --cov-report=html --cov-report=term
+
+# Minimum Python version
+minversion = 3.8
+
diff --git a/requirements-test.txt b/requirements-test.txt
new file mode 100644
index 0000000..958e4ef
--- /dev/null
+++ b/requirements-test.txt
@@ -0,0 +1,21 @@
+# Testing dependencies for fa2_modified
+
+# Core testing framework
+pytest>=7.0.0
+pytest-cov>=4.0.0
+
+# Core dependencies (required)
+numpy>=1.20.0
+scipy>=1.7.0
+tqdm>=4.60.0
+
+# Optional integration testing dependencies
+networkx>=2.5
+python-igraph>=0.9.0; platform_system != "Windows"
+
+# Code quality tools
+flake8>=6.0.0
+black>=23.0.0
+isort>=5.12.0
+mypy>=1.0.0
+
diff --git a/setup.py b/setup.py
index 4edec1f..953253b 100644
--- a/setup.py
+++ b/setup.py
@@ -44,7 +44,7 @@
 
 setup(
     name='fa2_modified',
-    version='0.3.10',
+    version='0.3.11',
     description='The fastest ForceAtlas2 algorithm for Python (and NetworkX)',
     long_description_content_type='text/markdown',
     long_description=long_description,
diff --git a/tests/__init__.py b/tests/__init__.py
new file mode 100644
index 0000000..22b536e
--- /dev/null
+++ b/tests/__init__.py
@@ -0,0 +1 @@
+# Test package for fa2_modified
diff --git a/tests/conftest.py b/tests/conftest.py
new file mode 100644
index 0000000..a40cef5
--- /dev/null
+++ b/tests/conftest.py
@@ -0,0 +1,86 @@
+"""
+Pytest configuration and fixtures for fa2_modified tests
+"""
+
+import numpy as np
+import pytest
+
+
+@pytest.fixture
+def simple_adjacency_matrix():
+    """Simple 4-node graph as adjacency matrix"""
+    return np.array([[0, 1, 0, 1], [1, 0, 1, 0], [0, 1, 0, 1], [1, 0, 1, 0]])
+
+
+@pytest.fixture
+def weighted_adjacency_matrix():
+    """Simple weighted graph as adjacency matrix"""
+    return np.array(
+        [
+            [0.0, 2.5, 0.0, 1.0],
+            [2.5, 0.0, 3.0, 0.0],
+            [0.0, 3.0, 0.0, 1.5],
+            [1.0, 0.0, 1.5, 0.0],
+        ]
+    )
+
+
+@pytest.fixture
+def triangle_graph():
+    """Triangle graph (3 nodes, all connected)"""
+    return np.array([[0, 1, 1], [1, 0, 1], [1, 1, 0]])
+
+
+@pytest.fixture
+def star_graph():
+    """Star graph (1 hub connected to 4 leaves)"""
+    n = 5
+    G = np.zeros((n, n))
+    for i in range(1, n):
+        G[0, i] = 1
+        G[i, 0] = 1
+    return G
+
+
+@pytest.fixture
+def disconnected_graph():
+    """Graph with two disconnected components"""
+    return np.array([[0, 1, 0, 0], [1, 0, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0]])
+
+
+@pytest.fixture
+def ring_graph():
+    """Ring graph (nodes in a circle)"""
+    n = 6
+    G = np.zeros((n, n))
+    for i in range(n):
+        G[i, (i + 1) % n] = 1
+        G[(i + 1) % n, i] = 1
+    return G
+
+
+@pytest.fixture
+def random_seed():
+    """Set random seed for reproducible tests"""
+    import random
+
+    np.random.seed(42)
+    random.seed(42)
+    yield
+    # Reset after test
+    np.random.seed(None)
+    random.seed(None)
+
+
+def pytest_configure(config):
+    """Configure pytest with custom markers"""
+    config.addinivalue_line(
+        "markers", "slow: marks tests as slow (deselect with '-m \"not slow\"')"
+    )
+    config.addinivalue_line("markers", "integration: marks tests as integration tests")
+    config.addinivalue_line(
+        "markers", "requires_networkx: requires networkx to be installed"
+    )
+    config.addinivalue_line(
+        "markers", "requires_igraph: requires igraph to be installed"
+    )
diff --git a/tests/test_api_consistency.py b/tests/test_api_consistency.py
new file mode 100644
index 0000000..8638ecf
--- /dev/null
+++ b/tests/test_api_consistency.py
@@ -0,0 +1,361 @@
+"""
+API consistency and error handling tests
+Tests inspired by external ForceAtlas2 implementations
+"""
+
+import numbers
+
+import numpy as np
+import pytest
+import scipy.sparse
+
+try:
+    import networkx as nx
+
+    NETWORKX_AVAILABLE = True
+except ImportError:
+    NETWORKX_AVAILABLE = False
+
+from fa2_modified import ForceAtlas2
+
+
+class TestForceAtlas2APIConsistency:
+    """Test API consistency across different input types"""
+
+    @pytest.mark.skipif(not NETWORKX_AVAILABLE, reason="NetworkX not installed")
+    def test_api_consistency_random_vs_custom_positions(self):
+        """
+        Run ForceAtlas2 layout on the same graph twice:
+        once with random initial positions and once with custom positions.
+        Both should return valid dictionaries with proper structure.
+        """
+        G = nx.path_graph(10)
+        forceatlas = ForceAtlas2(verbose=False)
+
+        # Layout with random initial positions
+        pos1 = forceatlas.forceatlas2_networkx_layout(G, pos=None, iterations=10)
+
+        # Layout with user-defined positions
+        pos_initial = {n: (0.1 * n, 0.1 * n) for n in G.nodes()}
+        pos2 = forceatlas.forceatlas2_networkx_layout(G, pos=pos_initial, iterations=10)
+
+        # Both outputs must be dictionaries with 10 items
+        assert isinstance(pos1, dict)
+        assert isinstance(pos2, dict)
+        assert len(pos1) == 10
+        assert len(pos2) == 10
+
+        # Keys should match
+        assert set(pos1.keys()) == set(pos2.keys())
+
+        # All values should be valid (x,y) tuples
+        for p in pos1.values():
+            assert isinstance(p, tuple)
+            assert len(p) == 2
+            assert all(isinstance(coord, numbers.Real) for coord in p)
+
+        for p in pos2.values():
+            assert isinstance(p, tuple)
+            assert len(p) == 2
+            assert all(isinstance(coord, numbers.Real) for coord in p)
+
+    @pytest.mark.skipif(not NETWORKX_AVAILABLE, reason="NetworkX not installed")
+    def test_api_consistency_custom_parameters(self):
+        """
+        Test ForceAtlas2 with various custom parameters to ensure
+        API consistency across different configurations.
+        """
+        G = nx.complete_graph(5)
+        pos = {i: (float(i), float(i)) for i in G.nodes()}
+
+        forceatlas = ForceAtlas2(
+            outboundAttractionDistribution=True,
+            edgeWeightInfluence=0.5,
+            jitterTolerance=5.0,
+            barnesHutOptimize=True,
+            barnesHutTheta=1.5,
+            scalingRatio=1.0,
+            strongGravityMode=True,
+            gravity=0.5,
+            verbose=False,
+        )
+
+        new_pos = forceatlas.forceatlas2_networkx_layout(G, pos=pos, iterations=20)
+
+        # Verify output structure
+        assert isinstance(new_pos, dict)
+        assert len(new_pos) == 5
+
+        for key, value in new_pos.items():
+            assert isinstance(value, tuple)
+            assert len(value) == 2
+            assert all(isinstance(coord, numbers.Real) for coord in value)
+
+    def test_api_consistency_numpy_sparse(self):
+        """Test that numpy array and sparse matrix give compatible results"""
+        # Create same graph as numpy array and sparse matrix
+        G_numpy = np.array(
+            [[0, 1, 1, 0], [1, 0, 1, 1], [1, 1, 0, 1], [0, 1, 1, 0]], dtype=float
+        )
+
+        G_sparse = scipy.sparse.lil_matrix(G_numpy)
+
+        forceatlas = ForceAtlas2(verbose=False)
+
+        # Use same seed for both
+        np.random.seed(42)
+        import random
+
+        random.seed(42)
+        pos_numpy = forceatlas.forceatlas2(G_numpy, iterations=10)
+
+        np.random.seed(42)
+        random.seed(42)
+        pos_sparse = forceatlas.forceatlas2(G_sparse, iterations=10)
+
+        # Both should return same number of positions
+        assert len(pos_numpy) == len(pos_sparse) == 4
+
+        # Positions should be very similar (identical with same seed)
+        for i in range(4):
+            assert pos_numpy[i][0] == pytest.approx(pos_sparse[i][0], abs=1e-10)
+            assert pos_numpy[i][1] == pytest.approx(pos_sparse[i][1], abs=1e-10)
+
+
+class TestForceAtlas2ErrorHandling:
+    """Test error handling for invalid inputs"""
+
+    def test_invalid_matrix_non_square(self):
+        """Test that non-square matrix raises AssertionError"""
+        forceatlas = ForceAtlas2(verbose=False)
+
+        with pytest.raises(AssertionError):
+            G = np.array([[0, 1], [1, 0], [0, 1]], dtype=float)
+            forceatlas.init(G, pos=None)
+
+    def test_invalid_matrix_non_symmetric(self):
+        """Test that non-symmetric matrix raises AssertionError"""
+        forceatlas = ForceAtlas2(verbose=False)
+
+        with pytest.raises(AssertionError):
+            G = np.array([[0, 1, 0], [0, 0, 1], [1, 0, 0]], dtype=float)
+            forceatlas.init(G, pos=None)
+
+    def test_invalid_position_type(self):
+        """Test that invalid position type raises AssertionError"""
+        forceatlas = ForceAtlas2(verbose=False)
+        G = np.array([[0, 1], [1, 0]], dtype=float)
+
+        with pytest.raises(AssertionError):
+            forceatlas.init(G, pos="invalid")
+
+    def test_invalid_position_shape(self):
+        """Test that position array with wrong shape raises error"""
+        forceatlas = ForceAtlas2(verbose=False)
+        G = np.array([[0, 1], [1, 0]], dtype=float)
+
+        # Position array with wrong number of nodes
+        with pytest.raises((AssertionError, IndexError, ValueError)):
+            pos = np.array([[0.0, 0.0]])  # Only 1 position for 2 nodes
+            forceatlas.forceatlas2(G, pos=pos, iterations=10)
+
+    def test_invalid_graph_type(self):
+        """Test that invalid graph type raises AssertionError"""
+        forceatlas = ForceAtlas2(verbose=False)
+
+        with pytest.raises(AssertionError):
+            forceatlas.init([[0, 1], [1, 0]], pos=None)
+
+    @pytest.mark.skipif(not NETWORKX_AVAILABLE, reason="NetworkX not installed")
+    def test_networkx_invalid_graph_type(self):
+        """Test that non-NetworkX object raises AssertionError"""
+        forceatlas = ForceAtlas2(verbose=False)
+
+        with pytest.raises(AssertionError):
+            forceatlas.forceatlas2_networkx_layout(
+                "not a NetworkX graph", pos=None, iterations=10
+            )
+
+    @pytest.mark.skipif(not NETWORKX_AVAILABLE, reason="NetworkX not installed")
+    def test_networkx_invalid_position_type(self):
+        """Test that invalid position type for NetworkX raises AssertionError"""
+        G = nx.complete_graph(5)
+        forceatlas = ForceAtlas2(verbose=False)
+
+        with pytest.raises(AssertionError):
+            forceatlas.forceatlas2_networkx_layout(G, pos="invalid", iterations=10)
+
+    @pytest.mark.skipif(not NETWORKX_AVAILABLE, reason="NetworkX not installed")
+    def test_networkx_position_list_not_dict(self):
+        """Test that position as list (not dict) raises AssertionError"""
+        G = nx.complete_graph(3)
+        forceatlas = ForceAtlas2(verbose=False)
+
+        with pytest.raises(AssertionError):
+            forceatlas.forceatlas2_networkx_layout(
+                G, pos=[(0, 0), (1, 1), (2, 2)], iterations=10
+            )
+
+
+class TestForceAtlas2InitializationEdgeCases:
+    """Test edge cases in ForceAtlas2 initialization"""
+
+    def test_init_returns_correct_types(self):
+        """Test that init returns lists of correct types"""
+        G = np.array([[0, 1], [1, 0]], dtype=float)
+        forceatlas = ForceAtlas2(verbose=False)
+        nodes, edges = forceatlas.init(G, pos=None)
+
+        assert isinstance(nodes, list)
+        assert isinstance(edges, list)
+        assert all(hasattr(n, "mass") for n in nodes)
+        assert all(hasattr(e, "weight") for e in edges)
+
+    def test_init_edge_count(self):
+        """Test that edges are counted correctly (undirected)"""
+        # Triangle graph: 3 edges
+        G = np.array([[0, 1, 1], [1, 0, 1], [1, 1, 0]], dtype=float)
+
+        forceatlas = ForceAtlas2(verbose=False)
+        nodes, edges = forceatlas.init(G, pos=None)
+
+        # Should have 3 edges (not 6, since we only count upper triangle)
+        assert len(edges) == 3
+
+    def test_init_node_mass_calculation(self):
+        """Test that node mass is calculated as 1 + degree"""
+        G = np.array(
+            [[0, 1, 1, 1], [1, 0, 0, 0], [1, 0, 0, 0], [1, 0, 0, 0]], dtype=float
+        )
+
+        forceatlas = ForceAtlas2(verbose=False)
+        nodes, edges = forceatlas.init(G, pos=None)
+
+        # Node 0 has degree 3, so mass should be 4
+        assert nodes[0].mass == 4.0
+
+        # Nodes 1, 2, 3 each have degree 1, so mass should be 2
+        assert nodes[1].mass == 2.0
+        assert nodes[2].mass == 2.0
+        assert nodes[3].mass == 2.0
+
+    def test_init_with_weighted_edges(self):
+        """Test that weighted edges preserve their weights"""
+        G = np.array([[0.0, 2.5, 0.0], [2.5, 0.0, 3.7], [0.0, 3.7, 0.0]], dtype=float)
+
+        forceatlas = ForceAtlas2(verbose=False)
+        nodes, edges = forceatlas.init(G, pos=None)
+
+        # Check that edge weights are preserved
+        weights = sorted([e.weight for e in edges])
+        assert weights == pytest.approx([2.5, 3.7])
+
+    def test_init_with_custom_positions(self):
+        """Test that custom positions are used correctly"""
+        G = np.array([[0, 1], [1, 0]], dtype=float)
+        pos = np.array([[1.5, 2.5], [3.5, 4.5]])
+
+        forceatlas = ForceAtlas2(verbose=False)
+        nodes, edges = forceatlas.init(G, pos=pos)
+
+        assert nodes[0].x == pytest.approx(1.5)
+        assert nodes[0].y == pytest.approx(2.5)
+        assert nodes[1].x == pytest.approx(3.5)
+        assert nodes[1].y == pytest.approx(4.5)
+
+
+class TestForceAtlas2OutputFormat:
+    """Test output format consistency"""
+
+    def test_output_is_list_of_tuples(self):
+        """Test that forceatlas2 returns list of (x,y) tuples"""
+        G = np.array([[0, 1, 1], [1, 0, 1], [1, 1, 0]], dtype=float)
+        forceatlas = ForceAtlas2(verbose=False)
+        positions = forceatlas.forceatlas2(G, iterations=10)
+
+        assert isinstance(positions, list)
+        assert len(positions) == 3
+
+        for pos in positions:
+            assert isinstance(pos, tuple)
+            assert len(pos) == 2
+            assert isinstance(pos[0], numbers.Real)
+            assert isinstance(pos[1], numbers.Real)
+
+    @pytest.mark.skipif(not NETWORKX_AVAILABLE, reason="NetworkX not installed")
+    def test_networkx_output_is_dict(self):
+        """Test that NetworkX layout returns dict"""
+        G = nx.path_graph(5)
+        forceatlas = ForceAtlas2(verbose=False)
+        positions = forceatlas.forceatlas2_networkx_layout(G, iterations=10)
+
+        assert isinstance(positions, dict)
+        assert len(positions) == 5
+        assert set(positions.keys()) == set(G.nodes())
+
+        for node, pos in positions.items():
+            assert isinstance(pos, tuple)
+            assert len(pos) == 2
+            assert isinstance(pos[0], numbers.Real)
+            assert isinstance(pos[1], numbers.Real)
+
+    def test_positions_are_finite(self):
+        """Test that all positions are finite (not NaN or inf)"""
+        G = np.array(
+            [[0, 1, 1, 0], [1, 0, 1, 1], [1, 1, 0, 1], [0, 1, 1, 0]], dtype=float
+        )
+
+        forceatlas = ForceAtlas2(verbose=False)
+        positions = forceatlas.forceatlas2(G, iterations=50)
+
+        for pos in positions:
+            assert np.isfinite(pos[0])
+            assert np.isfinite(pos[1])
+
+    def test_positions_are_not_all_same(self):
+        """Test that positions are not all identical (layout has spread)"""
+        G = np.array([[0, 1, 1], [1, 0, 1], [1, 1, 0]], dtype=float)
+
+        forceatlas = ForceAtlas2(verbose=False)
+        positions = forceatlas.forceatlas2(G, iterations=100)
+
+        # Check that not all positions are identical
+        x_coords = [pos[0] for pos in positions]
+        y_coords = [pos[1] for pos in positions]
+
+        # Standard deviation should be non-zero
+        assert np.std(x_coords) > 1e-6 or np.std(y_coords) > 1e-6
+
+
+class TestForceAtlas2ParameterValidation:
+    """Test parameter validation"""
+
+    def test_unimplemented_linLogMode_raises_error(self):
+        """Test that linLogMode=True raises AssertionError"""
+        with pytest.raises(AssertionError):
+            ForceAtlas2(linLogMode=True)
+
+    def test_unimplemented_adjustSizes_raises_error(self):
+        """Test that adjustSizes=True raises AssertionError"""
+        with pytest.raises(AssertionError):
+            ForceAtlas2(adjustSizes=True)
+
+    def test_unimplemented_multiThreaded_raises_error(self):
+        """Test that multiThreaded=True raises AssertionError"""
+        with pytest.raises(AssertionError):
+            ForceAtlas2(multiThreaded=True)
+
+    def test_valid_parameters_accepted(self):
+        """Test that all valid parameter combinations are accepted"""
+        valid_configs = [
+            {},
+            {"barnesHutOptimize": True, "barnesHutTheta": 1.5},
+            {"strongGravityMode": True, "gravity": 2.0},
+            {"outboundAttractionDistribution": True, "edgeWeightInfluence": 0.5},
+            {"scalingRatio": 5.0, "jitterTolerance": 2.0},
+            {"verbose": False},
+        ]
+
+        for config in valid_configs:
+            fa2 = ForceAtlas2(**config)
+            assert fa2 is not None
diff --git a/tests/test_fa2util.py b/tests/test_fa2util.py
new file mode 100644
index 0000000..3088218
--- /dev/null
+++ b/tests/test_fa2util.py
@@ -0,0 +1,583 @@
+"""
+Tests for fa2util module - testing force calculations and utility classes
+"""
+
+from math import sqrt
+
+import numpy as np
+import pytest
+
+from fa2_modified import fa2util
+
+
+class TestNode:
+    """Test Node class"""
+
+    def test_node_initialization(self):
+        """Test Node object initialization"""
+        node = fa2util.Node()
+        assert node.mass == 0.0
+        assert node.old_dx == 0.0
+        assert node.old_dy == 0.0
+        assert node.dx == 0.0
+        assert node.dy == 0.0
+        assert node.x == 0.0
+        assert node.y == 0.0
+
+    def test_node_attributes(self):
+        """Test Node attribute assignment"""
+        node = fa2util.Node()
+        node.mass = 5.0
+        node.x = 10.0
+        node.y = 20.0
+        node.dx = 1.0
+        node.dy = 2.0
+
+        assert node.mass == 5.0
+        assert node.x == 10.0
+        assert node.y == 20.0
+        assert node.dx == 1.0
+        assert node.dy == 2.0
+
+
+class TestEdge:
+    """Test Edge class"""
+
+    def test_edge_initialization(self):
+        """Test Edge object initialization"""
+        edge = fa2util.Edge()
+        assert edge.node1 == -1
+        assert edge.node2 == -1
+        assert edge.weight == 0.0
+
+    def test_edge_attributes(self):
+        """Test Edge attribute assignment"""
+        edge = fa2util.Edge()
+        edge.node1 = 0
+        edge.node2 = 1
+        edge.weight = 2.5
+
+        assert edge.node1 == 0
+        assert edge.node2 == 1
+        assert edge.weight == 2.5
+
+
+class TestLinRepulsion:
+    """Test linear repulsion force function"""
+
+    def test_repulsion_basic(self):
+        """Test basic repulsion between two nodes"""
+        n1 = fa2util.Node()
+        n1.x = 0.0
+        n1.y = 0.0
+        n1.mass = 1.0
+
+        n2 = fa2util.Node()
+        n2.x = 1.0
+        n2.y = 0.0
+        n2.mass = 1.0
+
+        fa2util.linRepulsion(n1, n2, coefficient=1.0)
+
+        # Nodes should be pushed apart
+        assert n1.dx < 0  # n1 pushed left
+        assert n2.dx > 0  # n2 pushed right
+        assert abs(n1.dx) == abs(n2.dx)  # Equal and opposite forces
+        assert n1.dy == 0.0  # No vertical movement
+        assert n2.dy == 0.0
+
+    def test_repulsion_zero_distance(self):
+        """Test repulsion when nodes are at the same position"""
+        n1 = fa2util.Node()
+        n1.x = 0.0
+        n1.y = 0.0
+        n1.mass = 1.0
+
+        n2 = fa2util.Node()
+        n2.x = 0.0
+        n2.y = 0.0
+        n2.mass = 1.0
+
+        fa2util.linRepulsion(n1, n2, coefficient=1.0)
+
+        # No force should be applied when distance is zero
+        assert n1.dx == 0.0
+        assert n1.dy == 0.0
+        assert n2.dx == 0.0
+        assert n2.dy == 0.0
+
+    def test_repulsion_different_masses(self):
+        """Test repulsion with different node masses"""
+        n1 = fa2util.Node()
+        n1.x = 0.0
+        n1.y = 0.0
+        n1.mass = 2.0
+
+        n2 = fa2util.Node()
+        n2.x = 1.0
+        n2.y = 0.0
+        n2.mass = 1.0
+
+        fa2util.linRepulsion(n1, n2, coefficient=1.0)
+
+        # Forces should be equal and opposite despite different masses
+        assert abs(n1.dx) == abs(n2.dx)
+
+
+class TestLinGravity:
+    """Test linear gravity force function"""
+
+    def test_gravity_basic(self):
+        """Test basic gravity pulling node toward origin"""
+        n = fa2util.Node()
+        n.x = 10.0
+        n.y = 0.0
+        n.mass = 1.0
+
+        fa2util.linGravity(n, g=1.0)
+
+        # Node should be pulled toward origin
+        assert n.dx < 0
+        assert n.dy == 0.0
+
+    def test_gravity_at_origin(self):
+        """Test gravity when node is at origin"""
+        n = fa2util.Node()
+        n.x = 0.0
+        n.y = 0.0
+        n.mass = 1.0
+
+        fa2util.linGravity(n, g=1.0)
+
+        # No force when at origin
+        assert n.dx == 0.0
+        assert n.dy == 0.0
+
+    def test_gravity_diagonal(self):
+        """Test gravity pulling diagonally"""
+        n = fa2util.Node()
+        n.x = 3.0
+        n.y = 4.0
+        n.mass = 1.0
+
+        fa2util.linGravity(n, g=1.0)
+
+        # Both components should pull toward origin
+        assert n.dx < 0
+        assert n.dy < 0
+
+        # Check that the direction is correct
+        distance = sqrt(n.x**2 + n.y**2)
+        expected_dx = -(n.x / distance) * n.mass
+        expected_dy = -(n.y / distance) * n.mass
+
+        assert abs(n.dx - expected_dx) < 1e-10
+        assert abs(n.dy - expected_dy) < 1e-10
+
+
+class TestStrongGravity:
+    """Test strong gravity force function"""
+
+    def test_strong_gravity_basic(self):
+        """Test strong gravity force"""
+        n = fa2util.Node()
+        n.x = 10.0
+        n.y = 5.0
+        n.mass = 1.0
+
+        fa2util.strongGravity(n, g=1.0, coefficient=1.0)
+
+        # Should pull toward origin
+        assert n.dx < 0
+        assert n.dy < 0
+
+    def test_strong_gravity_at_origin(self):
+        """Test strong gravity when at origin"""
+        n = fa2util.Node()
+        n.x = 0.0
+        n.y = 0.0
+        n.mass = 1.0
+
+        fa2util.strongGravity(n, g=1.0, coefficient=1.0)
+
+        # No force at origin
+        assert n.dx == 0.0
+        assert n.dy == 0.0
+
+
+class TestLinAttraction:
+    """Test linear attraction force function"""
+
+    def test_attraction_basic(self):
+        """Test basic attraction between connected nodes"""
+        n1 = fa2util.Node()
+        n1.x = 0.0
+        n1.y = 0.0
+        n1.mass = 1.0
+
+        n2 = fa2util.Node()
+        n2.x = 10.0
+        n2.y = 0.0
+        n2.mass = 1.0
+
+        fa2util.linAttraction(
+            n1, n2, e=1.0, distributedAttraction=False, coefficient=1.0
+        )
+
+        # Nodes should be pulled together
+        assert n1.dx > 0  # n1 pulled right
+        assert n2.dx < 0  # n2 pulled left
+        assert abs(n1.dx) == abs(n2.dx)  # Equal and opposite
+
+    def test_attraction_distributed(self):
+        """Test distributed attraction (hub dissuasion)"""
+        n1 = fa2util.Node()
+        n1.x = 0.0
+        n1.y = 0.0
+        n1.mass = 10.0  # Hub node with high mass
+
+        n2 = fa2util.Node()
+        n2.x = 10.0
+        n2.y = 0.0
+        n2.mass = 1.0
+
+        fa2util.linAttraction(
+            n1, n2, e=1.0, distributedAttraction=True, coefficient=1.0
+        )
+
+        # With distributed attraction, force on n1 is reduced by its mass
+        assert n1.dx > 0
+        assert n2.dx < 0
+
+    def test_attraction_weighted_edge(self):
+        """Test attraction with edge weight"""
+        n1 = fa2util.Node()
+        n1.x = 0.0
+        n1.y = 0.0
+        n1.mass = 1.0
+
+        n2 = fa2util.Node()
+        n2.x = 10.0
+        n2.y = 0.0
+        n2.mass = 1.0
+
+        fa2util.linAttraction(
+            n1, n2, e=2.0, distributedAttraction=False, coefficient=1.0
+        )
+
+        # Stronger edge weight should result in stronger force
+        assert n1.dx > 0
+        assert n2.dx < 0
+
+
+class TestApplyForces:
+    """Test force application functions"""
+
+    def test_apply_repulsion(self):
+        """Test applying repulsion to multiple nodes"""
+        nodes = []
+        for i in range(3):
+            n = fa2util.Node()
+            n.x = float(i)
+            n.y = 0.0
+            n.mass = 1.0
+            nodes.append(n)
+
+        fa2util.apply_repulsion(nodes, coefficient=1.0)
+
+        # Middle node should have zero net force (symmetry)
+        # End nodes should be pushed outward
+        assert nodes[0].dx < 0  # Left node pushed left
+        assert nodes[2].dx > 0  # Right node pushed right
+
+    def test_apply_gravity(self):
+        """Test applying gravity to multiple nodes"""
+        nodes = []
+        for i in range(3):
+            n = fa2util.Node()
+            n.x = float(i + 1)
+            n.y = float(i + 1)
+            n.mass = 1.0
+            nodes.append(n)
+
+        fa2util.apply_gravity(
+            nodes, gravity=1.0, scalingRatio=1.0, useStrongGravity=False
+        )
+
+        # All nodes should be pulled toward origin
+        for node in nodes:
+            assert node.dx < 0
+            assert node.dy < 0
+
+    def test_apply_strong_gravity(self):
+        """Test applying strong gravity to multiple nodes"""
+        nodes = []
+        for i in range(3):
+            n = fa2util.Node()
+            n.x = float(i + 1)
+            n.y = float(i + 1)
+            n.mass = 1.0
+            nodes.append(n)
+
+        fa2util.apply_gravity(
+            nodes, gravity=1.0, scalingRatio=1.0, useStrongGravity=True
+        )
+
+        # All nodes should be pulled toward origin
+        for node in nodes:
+            assert node.dx < 0
+            assert node.dy < 0
+
+    def test_apply_attraction(self):
+        """Test applying attraction forces along edges"""
+        nodes = []
+        for i in range(3):
+            n = fa2util.Node()
+            n.x = float(i * 10)
+            n.y = 0.0
+            n.mass = 1.0
+            nodes.append(n)
+
+        edges = []
+        # Connect node 0 to node 1
+        e = fa2util.Edge()
+        e.node1 = 0
+        e.node2 = 1
+        e.weight = 1.0
+        edges.append(e)
+
+        fa2util.apply_attraction(
+            nodes,
+            edges,
+            distributedAttraction=False,
+            coefficient=1.0,
+            edgeWeightInfluence=1.0,
+        )
+
+        # Node 0 and 1 should be pulled together
+        assert nodes[0].dx > 0
+        assert nodes[1].dx < 0
+        # Node 2 should have no force (not connected)
+        assert nodes[2].dx == 0.0
+
+
+class TestRegion:
+    """Test Barnes-Hut Region class"""
+
+    def test_region_initialization(self):
+        """Test Region initialization with nodes"""
+        nodes = []
+        for i in range(4):
+            n = fa2util.Node()
+            n.x = float(i)
+            n.y = float(i)
+            n.mass = 1.0
+            nodes.append(n)
+
+        region = fa2util.Region(nodes)
+
+        assert region.mass > 0
+        assert region.size > 0
+        assert len(region.nodes) == 4
+        assert len(region.subregions) == 0
+
+    def test_region_mass_center(self):
+        """Test region mass center calculation"""
+        nodes = []
+        # Create nodes in a square pattern
+        positions = [(0, 0), (0, 10), (10, 0), (10, 10)]
+        for x, y in positions:
+            n = fa2util.Node()
+            n.x = float(x)
+            n.y = float(y)
+            n.mass = 1.0
+            nodes.append(n)
+
+        region = fa2util.Region(nodes)
+
+        # Center should be at (5, 5)
+        assert abs(region.massCenterX - 5.0) < 1e-10
+        assert abs(region.massCenterY - 5.0) < 1e-10
+
+    def test_region_build_subregions(self):
+        """Test building subregions for Barnes-Hut"""
+        nodes = []
+        # Create nodes in different quadrants
+        positions = [(0, 0), (0, 10), (10, 0), (10, 10)]
+        for x, y in positions:
+            n = fa2util.Node()
+            n.x = float(x)
+            n.y = float(y)
+            n.mass = 1.0
+            nodes.append(n)
+
+        region = fa2util.Region(nodes)
+        region.buildSubRegions()
+
+        # Should have created subregions
+        assert len(region.subregions) > 0
+
+    def test_region_single_node(self):
+        """Test region with single node"""
+        n = fa2util.Node()
+        n.x = 5.0
+        n.y = 5.0
+        n.mass = 1.0
+
+        region = fa2util.Region([n])
+        region.buildSubRegions()
+
+        # Single node should not create subregions
+        assert len(region.subregions) == 0
+
+
+class TestAdjustSpeedAndApplyForces:
+    """Test speed adjustment and force application"""
+
+    def test_adjust_speed_basic(self):
+        """Test basic speed adjustment"""
+        nodes = []
+        for i in range(5):
+            n = fa2util.Node()
+            n.x = float(i)
+            n.y = float(i)
+            n.mass = 1.0
+            n.dx = 0.1
+            n.dy = 0.1
+            n.old_dx = 0.0
+            n.old_dy = 0.0
+            nodes.append(n)
+
+        result = fa2util.adjustSpeedAndApplyForces(
+            nodes, speed=1.0, speedEfficiency=1.0, jitterTolerance=1.0
+        )
+
+        assert "speed" in result
+        assert "speedEfficiency" in result
+        assert result["speed"] > 0
+        assert result["speedEfficiency"] > 0
+
+    def test_adjust_speed_applies_forces(self):
+        """Test that forces are applied to node positions"""
+        nodes = []
+        for i in range(3):
+            n = fa2util.Node()
+            n.x = float(i)
+            n.y = float(i)
+            n.mass = 1.0
+            n.dx = 1.0
+            n.dy = -1.0
+            n.old_dx = 0.5
+            n.old_dy = -0.5
+            nodes.append(n)
+
+        old_positions = [(n.x, n.y) for n in nodes]
+
+        result = fa2util.adjustSpeedAndApplyForces(
+            nodes, speed=1.0, speedEfficiency=1.0, jitterTolerance=1.0
+        )
+
+        # Positions should have changed
+        for i, n in enumerate(nodes):
+            assert n.x != old_positions[i][0] or n.y != old_positions[i][1]
+
+
+class TestEdgeCases:
+    """Test edge cases and error conditions"""
+
+    def test_empty_nodes_list(self):
+        """Test with empty nodes list"""
+        nodes = []
+        fa2util.apply_repulsion(nodes, coefficient=1.0)
+        fa2util.apply_gravity(nodes, gravity=1.0, scalingRatio=1.0)
+        # Should not raise errors
+
+    def test_single_node_forces(self):
+        """Test forces on single node"""
+        nodes = []
+        n = fa2util.Node()
+        n.x = 5.0
+        n.y = 5.0
+        n.mass = 1.0
+        nodes.append(n)
+
+        fa2util.apply_repulsion(nodes, coefficient=1.0)
+        fa2util.apply_gravity(nodes, gravity=1.0, scalingRatio=1.0)
+
+        # Gravity should pull toward origin
+        assert n.dx < 0
+        assert n.dy < 0
+
+    def test_edge_weight_zero(self):
+        """Test attraction with zero edge weight"""
+        nodes = []
+        for i in range(2):
+            n = fa2util.Node()
+            n.x = float(i * 10)
+            n.y = 0.0
+            n.mass = 1.0
+            nodes.append(n)
+
+        edges = []
+        e = fa2util.Edge()
+        e.node1 = 0
+        e.node2 = 1
+        e.weight = 0.0
+        edges.append(e)
+
+        fa2util.apply_attraction(
+            nodes,
+            edges,
+            distributedAttraction=False,
+            coefficient=1.0,
+            edgeWeightInfluence=1.0,
+        )
+
+        # Zero weight edge should produce no force
+        assert nodes[0].dx == 0.0
+        assert nodes[1].dx == 0.0
+
+    def test_edge_weight_influence_variations(self):
+        """Test different edge weight influence values"""
+        nodes = []
+        for i in range(2):
+            n = fa2util.Node()
+            n.x = float(i * 10)
+            n.y = 0.0
+            n.mass = 1.0
+            nodes.append(n)
+
+        edges = []
+        e = fa2util.Edge()
+        e.node1 = 0
+        e.node2 = 1
+        e.weight = 2.0
+        edges.append(e)
+
+        # Test with edgeWeightInfluence = 0 (ignore weights)
+        fa2util.apply_attraction(
+            nodes,
+            edges,
+            distributedAttraction=False,
+            coefficient=1.0,
+            edgeWeightInfluence=0.0,
+        )
+        force_0 = nodes[0].dx
+
+        # Reset
+        for n in nodes:
+            n.dx = 0.0
+            n.dy = 0.0
+
+        # Test with edgeWeightInfluence = 1 (use weights)
+        fa2util.apply_attraction(
+            nodes,
+            edges,
+            distributedAttraction=False,
+            coefficient=1.0,
+            edgeWeightInfluence=1.0,
+        )
+        force_1 = nodes[0].dx
+
+        # With weight=2 and influence=1, force should be stronger
+        assert abs(force_1) > abs(force_0)
diff --git a/tests/test_fa2util_precise.py b/tests/test_fa2util_precise.py
new file mode 100644
index 0000000..f5c9822
--- /dev/null
+++ b/tests/test_fa2util_precise.py
@@ -0,0 +1,507 @@
+"""
+Precise numerical tests for fa2util module with exact value checking
+Tests inspired by external ForceAtlas2 implementations to ensure numerical accuracy
+"""
+
+import math
+
+import pytest
+
+from fa2_modified import fa2util
+
+
+class TestLinRepulsionPrecise:
+    """Precise numerical tests for linear repulsion"""
+
+    def test_lin_repulsion_exact_values(self):
+        """Test linear repulsion with exact numerical expectations"""
+        # Two nodes at (1,1) and (2,2) with mass 2 each
+        a = fa2util.Node()
+        b = fa2util.Node()
+        a.mass = 2.0
+        b.mass = 2.0
+        a.x, a.y = 1.0, 1.0
+        b.x, b.y = 2.0, 2.0
+
+        # Reset forces
+        a.dx = a.dy = b.dx = b.dy = 0.0
+
+        fa2util.linRepulsion(a, b, coefficient=1.0)
+
+        # dx = 1-2 = -1, dy = -1, dist_sq = 2
+        # factor = (coefficient * mass_a * mass_b) / dist_sq = (1 * 2 * 2) / 2 = 2
+        # a.dx += -1 * 2 = -2, a.dy += -1 * 2 = -2
+        # b.dx -= -1 * 2 = 2, b.dy -= -1 * 2 = 2
+        assert a.dx == pytest.approx(-2.0)
+        assert a.dy == pytest.approx(-2.0)
+        assert b.dx == pytest.approx(2.0)
+        assert b.dy == pytest.approx(2.0)
+
+    def test_lin_repulsion_asymmetric_masses(self):
+        """Test repulsion with different masses"""
+        a = fa2util.Node()
+        b = fa2util.Node()
+        a.mass = 3.0
+        b.mass = 2.0
+        a.x, a.y = 0.0, 0.0
+        b.x, b.y = 1.0, 0.0
+        a.dx = a.dy = b.dx = b.dy = 0.0
+
+        fa2util.linRepulsion(a, b, coefficient=1.0)
+
+        # dx = 0-1 = -1, dist_sq = 1
+        # factor = (1 * 3 * 2) / 1 = 6
+        # a.dx += -1 * 6 = -6
+        # b.dx -= -1 * 6 = 6
+        assert a.dx == pytest.approx(-6.0)
+        assert a.dy == pytest.approx(0.0)
+        assert b.dx == pytest.approx(6.0)
+        assert b.dy == pytest.approx(0.0)
+
+
+class TestLinRepulsionRegion:
+    """Test linear repulsion between node and region"""
+
+    def test_lin_repulsion_region_basic(self):
+        """Test linear repulsion between a node and a region"""
+        node = fa2util.Node()
+        node.mass = 3.0
+        node.x = 5.0
+        node.y = 5.0
+        node.dx = node.dy = 0.0
+
+        # Create a region with TWO nodes (single node regions don't calculate mass)
+        dummy1 = fa2util.Node()
+        dummy1.mass = 1.0
+        dummy1.x = 3.0
+        dummy1.y = 3.0
+
+        dummy2 = fa2util.Node()
+        dummy2.mass = 1.0
+        dummy2.x = 3.0
+        dummy2.y = 3.0
+
+        region = fa2util.Region([dummy1, dummy2])
+
+        # The region should have mass from both nodes
+        assert region.mass == pytest.approx(2.0)
+
+        fa2util.linRepulsion_region(node, region, coefficient=1.0)
+
+        # dx = 5-3 = 2, dy = 2, dist_sq = 8
+        # factor = (1 * 3 * 2) / 8 = 0.75
+        # node.dx += 2 * 0.75 = 1.5
+        # node.dy += 2 * 0.75 = 1.5
+        assert node.dx == pytest.approx(1.5)
+        assert node.dy == pytest.approx(1.5)
+
+    def test_lin_repulsion_region_no_update_zero_distance(self):
+        """Test no update when node is at region center"""
+        node = fa2util.Node()
+        node.mass = 1.0
+        node.x = 3.0
+        node.y = 3.0
+        node.dx = node.dy = 0.0
+
+        dummy = fa2util.Node()
+        dummy.mass = 1.0
+        dummy.x = 3.0
+        dummy.y = 3.0
+        region = fa2util.Region([dummy])
+
+        fa2util.linRepulsion_region(node, region, coefficient=1.0)
+
+        # No force when at same position
+        assert node.dx == pytest.approx(0.0)
+        assert node.dy == pytest.approx(0.0)
+
+
+class TestLinGravityPrecise:
+    """Precise numerical tests for linear gravity"""
+
+    def test_lin_gravity_exact_values(self):
+        """Test linear gravity with exact numerical expectations"""
+        node = fa2util.Node()
+        node.mass = 2.0
+        node.x = 3.0
+        node.y = 4.0  # distance from origin = 5
+        node.dx = node.dy = 0.0
+
+        fa2util.linGravity(node, g=1.0)
+
+        # distance = sqrt(9 + 16) = 5
+        # factor = (mass * g) / distance = (2 * 1) / 5 = 0.4
+        # dx -= 3 * 0.4 = -1.2
+        # dy -= 4 * 0.4 = -1.6
+        assert node.dx == pytest.approx(-1.2)
+        assert node.dy == pytest.approx(-1.6)
+
+    def test_lin_gravity_different_gravity_values(self):
+        """Test linear gravity with different gravity constants"""
+        node = fa2util.Node()
+        node.mass = 1.0
+        node.x = 3.0
+        node.y = 4.0  # distance = 5
+        node.dx = node.dy = 0.0
+
+        fa2util.linGravity(node, g=2.0)
+
+        # factor = (1 * 2) / 5 = 0.4
+        # dx -= 3 * 0.4 = -1.2
+        # dy -= 4 * 0.4 = -1.6
+        assert node.dx == pytest.approx(-1.2)
+        assert node.dy == pytest.approx(-1.6)
+
+
+class TestStrongGravityPrecise:
+    """Precise numerical tests for strong gravity"""
+
+    def test_strong_gravity_exact_values(self):
+        """Test strong gravity with exact numerical expectations"""
+        node = fa2util.Node()
+        node.mass = 2.0
+        node.x = 3.0
+        node.y = 4.0
+        node.dx = node.dy = 0.0
+
+        fa2util.strongGravity(node, g=1.0, coefficient=1.0)
+
+        # factor = coefficient * mass * g = 1 * 2 * 1 = 2
+        # dx -= 3 * 2 = -6
+        # dy -= 4 * 2 = -8
+        assert node.dx == pytest.approx(-6.0)
+        assert node.dy == pytest.approx(-8.0)
+
+    def test_strong_gravity_with_scaling(self):
+        """Test strong gravity with different coefficient"""
+        node = fa2util.Node()
+        node.mass = 1.0
+        node.x = 2.0
+        node.y = 3.0
+        node.dx = node.dy = 0.0
+
+        fa2util.strongGravity(node, g=1.0, coefficient=2.0)
+
+        # factor = 2 * 1 * 1 = 2
+        # dx -= 2 * 2 = -4
+        # dy -= 3 * 2 = -6
+        assert node.dx == pytest.approx(-4.0)
+        assert node.dy == pytest.approx(-6.0)
+
+
+class TestLinAttractionPrecise:
+    """Precise numerical tests for linear attraction"""
+
+    def test_lin_attraction_exact_values(self):
+        """Test linear attraction with exact numerical expectations"""
+        a = fa2util.Node()
+        b = fa2util.Node()
+        a.mass = 2.0
+        b.mass = 3.0
+        a.x, a.y = 0.0, 0.0
+        b.x, b.y = 3.0, 4.0  # distance = 5
+        a.dx = a.dy = b.dx = b.dy = 0.0
+
+        edge_weight = 1.0
+
+        fa2util.linAttraction(
+            a, b, edge_weight, distributedAttraction=False, coefficient=2.0
+        )
+
+        # dx = 0-3 = -3, dy = 0-4 = -4
+        # factor = -coefficient * edge_weight = -2.0 * 1.0 = -2.0
+        # a.dx += -3 * -2.0 = 6.0
+        # a.dy += -4 * -2.0 = 8.0
+        # b.dx -= -3 * -2.0 = -6.0
+        # b.dy -= -4 * -2.0 = -8.0
+        assert a.dx == pytest.approx(6.0)
+        assert a.dy == pytest.approx(8.0)
+        assert b.dx == pytest.approx(-6.0)
+        assert b.dy == pytest.approx(-8.0)
+
+    def test_lin_attraction_distributed(self):
+        """Test distributed attraction (hub dissuasion)"""
+        a = fa2util.Node()
+        b = fa2util.Node()
+        a.mass = 4.0  # Hub with higher mass
+        b.mass = 1.0
+        a.x, a.y = 0.0, 0.0
+        b.x, b.y = 2.0, 0.0
+        a.dx = a.dy = b.dx = b.dy = 0.0
+
+        edge_weight = 1.0
+
+        fa2util.linAttraction(
+            a, b, edge_weight, distributedAttraction=True, coefficient=1.0
+        )
+
+        # dx = 0-2 = -2, dy = 0
+        # factor = -coefficient * edge_weight / mass_a = -1.0 * 1.0 / 4.0 = -0.25
+        # a.dx += -2 * -0.25 = 0.5
+        # b.dx -= -2 * -0.25 = -0.5
+        assert a.dx == pytest.approx(0.5)
+        assert a.dy == pytest.approx(0.0)
+        assert b.dx == pytest.approx(-0.5)
+        assert b.dy == pytest.approx(0.0)
+
+
+class TestApplyForcesPrecise:
+    """Precise tests for force application functions"""
+
+    def test_apply_repulsion_exact(self):
+        """Test apply_repulsion with exact values"""
+        a = fa2util.Node()
+        b = fa2util.Node()
+        a.mass = b.mass = 2.0
+        a.x, a.y = 0.0, 0.0
+        b.x, b.y = 1.0, 0.0
+        a.dx = a.dy = b.dx = b.dy = 0.0
+
+        fa2util.apply_repulsion([a, b], coefficient=2.0)
+
+        # dx = 0-1 = -1, dist_sq = 1
+        # factor = (2 * 2 * 2) / 1 = 8
+        # a.dx += -1 * 8 = -8
+        # b.dx -= -1 * 8 = 8
+        assert a.dx == pytest.approx(-8.0)
+        assert a.dy == pytest.approx(0.0)
+        assert b.dx == pytest.approx(8.0)
+        assert b.dy == pytest.approx(0.0)
+
+    def test_apply_gravity_linear_exact(self):
+        """Test apply_gravity with linear mode"""
+        node = fa2util.Node()
+        node.mass = 1.0
+        node.x, node.y = 3.0, 4.0  # distance = 5
+        node.dx = node.dy = 0.0
+
+        fa2util.apply_gravity(
+            [node], gravity=1.0, scalingRatio=1.0, useStrongGravity=False
+        )
+
+        # factor = (mass * gravity) / distance = (1 * 1) / 5 = 0.2
+        # dx -= 3 * 0.2 = -0.6
+        # dy -= 4 * 0.2 = -0.8
+        assert node.dx == pytest.approx(-0.6)
+        assert node.dy == pytest.approx(-0.8)
+
+    def test_apply_gravity_strong_exact(self):
+        """Test apply_gravity with strong gravity mode"""
+        node = fa2util.Node()
+        node.mass = 1.0
+        node.x, node.y = 3.0, 4.0
+        node.dx = node.dy = 0.0
+
+        fa2util.apply_gravity(
+            [node], gravity=1.0, scalingRatio=1.0, useStrongGravity=True
+        )
+
+        # factor = scaling_ratio * mass * gravity = 1 * 1 * 1 = 1
+        # dx -= 3 * 1 = -3
+        # dy -= 4 * 1 = -4
+        assert node.dx == pytest.approx(-3.0)
+        assert node.dy == pytest.approx(-4.0)
+
+    def test_apply_attraction_exact(self):
+        """Test apply_attraction with exact values"""
+        a = fa2util.Node()
+        b = fa2util.Node()
+        a.mass = b.mass = 1.0
+        a.x, a.y = 0.0, 0.0
+        b.x, b.y = 1.0, 0.0
+        a.dx = a.dy = b.dx = b.dy = 0.0
+
+        edge = fa2util.Edge()
+        edge.node1 = 0
+        edge.node2 = 1
+        edge.weight = 1.0
+        nodes = [a, b]
+
+        fa2util.apply_attraction(
+            nodes,
+            [edge],
+            distributedAttraction=False,
+            coefficient=1.0,
+            edgeWeightInfluence=1.0,
+        )
+
+        # dx = 0-1 = -1, factor = -1.0 * 1.0 = -1.0
+        # a.dx += -1 * -1.0 = 1.0
+        # b.dx -= -1 * -1.0 = -1.0
+        assert a.dx == pytest.approx(1.0)
+        assert a.dy == pytest.approx(0.0)
+        assert b.dx == pytest.approx(-1.0)
+        assert b.dy == pytest.approx(0.0)
+
+
+class TestRegionPrecise:
+    """Precise tests for Region class methods"""
+
+    def test_region_mass_center_calculation(self):
+        """Test precise calculation of region mass center"""
+        a = fa2util.Node()
+        b = fa2util.Node()
+        a.mass = 2.0
+        b.mass = 4.0
+        a.x, a.y = 0.0, 0.0
+        b.x, b.y = 4.0, 0.0
+
+        region = fa2util.Region([a, b])
+
+        # Total mass = 6
+        assert region.mass == pytest.approx(6.0)
+
+        # Center of mass: ((0*2 + 4*4)/6, (0*2 + 0*4)/6) = (16/6, 0) = (8/3, 0)
+        assert region.massCenterX == pytest.approx(8.0 / 3.0)
+        assert region.massCenterY == pytest.approx(0.0)
+
+    def test_region_size_calculation(self):
+        """Test region size calculation"""
+        # Create nodes at specific positions
+        nodes = []
+        positions = [(0, 0), (4, 0)]
+        for x, y in positions:
+            n = fa2util.Node()
+            n.mass = 2.0
+            n.x, n.y = float(x), float(y)
+            nodes.append(n)
+
+        region = fa2util.Region(nodes)
+
+        # massSumX = 0*2 + 4*2 = 8, mass = 4
+        # Center is at (8/4, 0) = (2, 0)
+        # Distance from (0,0) to (2,0) = 2, from (4,0) to (2,0) = 2
+        # Size should be 2 * max_distance = 2 * 2 = 4.0
+        expected_size = 4.0
+        assert region.size == pytest.approx(expected_size)
+
+    def test_region_build_subregions_quadrants(self):
+        """Test that subregions are built in correct quadrants"""
+        # Create 4 nodes in different quadrants around (2, 2)
+        nodes = []
+        positions = [(1, 3), (1, 1), (3, 3), (3, 1)]
+        for x, y in positions:
+            n = fa2util.Node()
+            n.mass = 1.0
+            n.x, n.y = float(x), float(y)
+            nodes.append(n)
+
+        region = fa2util.Region(nodes)
+        region.buildSubRegions()
+
+        # Should create 4 subregions (one for each quadrant)
+        assert len(region.subregions) == 4
+
+        # Each subregion should have exactly 1 node
+        for subregion in region.subregions:
+            assert len(subregion.nodes) == 1
+
+    def test_region_apply_force_single_node(self):
+        """Test applying force from a single-node region"""
+        target = fa2util.Node()
+        target.mass = 1.0
+        target.x = 5.0
+        target.y = 5.0
+        target.dx = target.dy = 0.0
+
+        source = fa2util.Node()
+        source.mass = 1.0
+        source.x = 1.0
+        source.y = 1.0
+
+        region = fa2util.Region([source])
+
+        # For a single-node region, apply_force calls linRepulsion
+        region.applyForce(target, theta=100.0, coefficient=1.0)
+
+        # dx = 5-1 = 4, dy = 4, dist_sq = 32
+        # factor = (1 * 1 * 1) / 32 = 0.03125
+        # target.dx += 4 * 0.03125 = 0.125
+        # target.dy += 4 * 0.03125 = 0.125
+        assert target.dx == pytest.approx(0.125)
+        assert target.dy == pytest.approx(0.125)
+
+    def test_region_apply_force_uses_approximation(self):
+        """Test that region uses approximation when theta condition is met"""
+        target = fa2util.Node()
+        target.mass = 1.0
+        target.x = 10.0
+        target.y = 10.0
+        target.dx = target.dy = 0.0
+
+        # Create a region with multiple nodes close together
+        nodes = []
+        for i in range(4):
+            n = fa2util.Node()
+            n.mass = 1.0
+            n.x = float(i)
+            n.y = 0.0
+            nodes.append(n)
+
+        region = fa2util.Region(nodes)
+        region.buildSubRegions()
+
+        # With large theta and large distance, should use region approximation
+        region.applyForce(target, theta=0.5, coefficient=1.0)
+
+        # Force should have been applied (non-zero)
+        assert target.dx != 0.0 or target.dy != 0.0
+
+
+class TestAdjustSpeedPrecise:
+    """Precise tests for speed adjustment"""
+
+    def test_adjust_speed_modifies_positions(self):
+        """Test that adjust_speed_and_apply_forces updates node positions"""
+        a = fa2util.Node()
+        b = fa2util.Node()
+        a.mass = b.mass = 1.0
+        a.x, a.y = 0.0, 0.0
+        b.x, b.y = 1.0, 0.0
+
+        # Set forces
+        a.old_dx, a.old_dy = 0.0, 0.0
+        a.dx, a.dy = 1.0, 0.0
+        b.old_dx, b.old_dy = 0.0, 0.0
+        b.dx, b.dy = -1.0, 0.0
+
+        old_a_x = a.x
+        old_b_x = b.x
+
+        result = fa2util.adjustSpeedAndApplyForces(
+            [a, b], speed=1.0, speedEfficiency=1.0, jitterTolerance=1.0
+        )
+
+        # Positions should have changed
+        assert a.x != old_a_x
+        assert b.x != old_b_x
+
+        # Result should contain updated speed values
+        assert "speed" in result
+        assert "speedEfficiency" in result
+        assert result["speed"] > 0
+        assert result["speedEfficiency"] > 0
+
+    def test_adjust_speed_returns_valid_dict(self):
+        """Test that adjust_speed returns properly formatted dictionary"""
+        nodes = []
+        for i in range(3):
+            n = fa2util.Node()
+            n.mass = 1.0
+            n.x = float(i)
+            n.y = 0.0
+            n.dx = 0.1
+            n.dy = 0.0
+            n.old_dx = 0.0
+            n.old_dy = 0.0
+            nodes.append(n)
+
+        result = fa2util.adjustSpeedAndApplyForces(
+            nodes, speed=1.0, speedEfficiency=1.0, jitterTolerance=1.0
+        )
+
+        assert isinstance(result, dict)
+        assert "speed" in result
+        assert "speedEfficiency" in result
+        assert isinstance(result["speed"], float)
+        assert isinstance(result["speedEfficiency"], float)
diff --git a/tests/test_forceatlas2.py b/tests/test_forceatlas2.py
new file mode 100644
index 0000000..0b99341
--- /dev/null
+++ b/tests/test_forceatlas2.py
@@ -0,0 +1,395 @@
+"""
+Tests for ForceAtlas2 main class and layout algorithms
+"""
+
+import numpy as np
+import pytest
+import scipy.sparse
+
+from fa2_modified import ForceAtlas2
+
+
+class TestForceAtlas2Initialization:
+    """Test ForceAtlas2 class initialization"""
+
+    def test_default_initialization(self):
+        """Test ForceAtlas2 with default parameters"""
+        fa2 = ForceAtlas2()
+
+        assert fa2.outboundAttractionDistribution is False
+        assert fa2.linLogMode is False
+        assert fa2.adjustSizes is False
+        assert fa2.edgeWeightInfluence == 1.0
+        assert fa2.jitterTolerance == 1.0
+        assert fa2.barnesHutOptimize is True
+        assert fa2.barnesHutTheta == 1.2
+        assert fa2.scalingRatio == 2.0
+        assert fa2.strongGravityMode is False
+        assert fa2.gravity == 1.0
+        assert fa2.verbose is True
+
+    def test_custom_parameters(self):
+        """Test ForceAtlas2 with custom parameters"""
+        fa2 = ForceAtlas2(
+            outboundAttractionDistribution=True,
+            edgeWeightInfluence=0.5,
+            jitterTolerance=2.0,
+            barnesHutOptimize=False,
+            barnesHutTheta=1.5,
+            scalingRatio=3.0,
+            strongGravityMode=True,
+            gravity=2.0,
+            verbose=False,
+        )
+
+        assert fa2.outboundAttractionDistribution is True
+        assert fa2.edgeWeightInfluence == 0.5
+        assert fa2.jitterTolerance == 2.0
+        assert fa2.barnesHutOptimize is False
+        assert fa2.barnesHutTheta == 1.5
+        assert fa2.scalingRatio == 3.0
+        assert fa2.strongGravityMode is True
+        assert fa2.gravity == 2.0
+        assert fa2.verbose is False
+
+    def test_unimplemented_features_raise_error(self):
+        """Test that unimplemented features raise assertion error"""
+        with pytest.raises(AssertionError):
+            ForceAtlas2(linLogMode=True)
+
+        with pytest.raises(AssertionError):
+            ForceAtlas2(adjustSizes=True)
+
+        with pytest.raises(AssertionError):
+            ForceAtlas2(multiThreaded=True)
+
+
+class TestForceAtlas2Init:
+    """Test the init method that prepares graph data"""
+
+    def test_init_numpy_array(self):
+        """Test init with numpy array"""
+        G = np.array([[0, 1, 0, 1], [1, 0, 1, 0], [0, 1, 0, 1], [1, 0, 1, 0]])
+
+        fa2 = ForceAtlas2(verbose=False)
+        nodes, edges = fa2.init(G)
+
+        assert len(nodes) == 4
+        assert len(edges) == 4  # 4 unique edges in undirected graph
+
+        # Check node masses (1 + number of connections)
+        assert nodes[0].mass == 3.0  # 1 + 2 edges
+        assert nodes[1].mass == 3.0
+        assert nodes[2].mass == 3.0
+        assert nodes[3].mass == 3.0
+
+    def test_init_sparse_matrix(self):
+        """Test init with scipy sparse matrix"""
+        G = scipy.sparse.lil_matrix(
+            [[0, 1, 0, 1], [1, 0, 1, 0], [0, 1, 0, 1], [1, 0, 1, 0]]
+        )
+
+        fa2 = ForceAtlas2(verbose=False)
+        nodes, edges = fa2.init(G)
+
+        assert len(nodes) == 4
+        assert len(edges) == 4
+
+    def test_init_with_positions(self):
+        """Test init with initial positions"""
+        G = np.array([[0, 1], [1, 0]])
+        pos = np.array([[0.0, 0.0], [1.0, 1.0]])
+
+        fa2 = ForceAtlas2(verbose=False)
+        nodes, edges = fa2.init(G, pos=pos)
+
+        assert nodes[0].x == 0.0
+        assert nodes[0].y == 0.0
+        assert nodes[1].x == 1.0
+        assert nodes[1].y == 1.0
+
+    def test_init_without_positions(self):
+        """Test init without initial positions (random)"""
+        G = np.array([[0, 1], [1, 0]])
+
+        fa2 = ForceAtlas2(verbose=False)
+        nodes, edges = fa2.init(G, pos=None)
+
+        # Should have random positions between 0 and 1
+        assert 0 <= nodes[0].x <= 1
+        assert 0 <= nodes[0].y <= 1
+        assert 0 <= nodes[1].x <= 1
+        assert 0 <= nodes[1].y <= 1
+
+    def test_init_weighted_graph(self):
+        """Test init with weighted edges"""
+        G = np.array([[0.0, 2.5, 0.0], [2.5, 0.0, 1.5], [0.0, 1.5, 0.0]])
+
+        fa2 = ForceAtlas2(verbose=False)
+        nodes, edges = fa2.init(G)
+
+        assert len(edges) == 2
+        # Check edge weights
+        weights = sorted([e.weight for e in edges])
+        assert weights == [1.5, 2.5]
+
+    def test_init_invalid_input(self):
+        """Test init with invalid input"""
+        fa2 = ForceAtlas2(verbose=False)
+
+        # Non-square matrix
+        with pytest.raises(AssertionError):
+            G = np.array([[0, 1], [1, 0], [0, 1]])
+            fa2.init(G)
+
+        # Non-symmetric matrix
+        with pytest.raises(AssertionError):
+            G = np.array([[0, 1], [0, 0]])
+            fa2.init(G)
+
+        # Invalid type
+        with pytest.raises(AssertionError):
+            fa2.init([[0, 1], [1, 0]])
+
+
+class TestForceAtlas2Layout:
+    """Test the main forceatlas2 layout method"""
+
+    def test_simple_graph_layout(self):
+        """Test layout on simple graph"""
+        G = np.array([[0, 1, 0, 0], [1, 0, 1, 0], [0, 1, 0, 1], [0, 0, 1, 0]])
+
+        fa2 = ForceAtlas2(verbose=False)
+        pos = fa2.forceatlas2(G, iterations=10)
+
+        assert len(pos) == 4
+        assert all(len(p) == 2 for p in pos)  # Each position is (x, y)
+        assert all(isinstance(p[0], float) and isinstance(p[1], float) for p in pos)
+
+    def test_layout_with_initial_positions(self):
+        """Test layout with provided initial positions"""
+        G = np.array([[0, 1], [1, 0]])
+        initial_pos = np.array([[0.0, 0.0], [10.0, 0.0]])
+
+        fa2 = ForceAtlas2(verbose=False)
+        pos = fa2.forceatlas2(G, pos=initial_pos, iterations=10)
+
+        assert len(pos) == 2
+        # Positions should have changed from initial
+        # (they should be pulled together by attraction)
+
+    def test_layout_barnes_hut_enabled(self):
+        """Test layout with Barnes-Hut optimization"""
+        G = np.array([[0, 1, 1, 0], [1, 0, 1, 1], [1, 1, 0, 1], [0, 1, 1, 0]])
+
+        fa2 = ForceAtlas2(barnesHutOptimize=True, verbose=False)
+        pos = fa2.forceatlas2(G, iterations=10)
+
+        assert len(pos) == 4
+
+    def test_layout_barnes_hut_disabled(self):
+        """Test layout without Barnes-Hut optimization"""
+        G = np.array([[0, 1, 1, 0], [1, 0, 1, 1], [1, 1, 0, 1], [0, 1, 1, 0]])
+
+        fa2 = ForceAtlas2(barnesHutOptimize=False, verbose=False)
+        pos = fa2.forceatlas2(G, iterations=10)
+
+        assert len(pos) == 4
+
+    def test_layout_strong_gravity(self):
+        """Test layout with strong gravity mode"""
+        G = np.array([[0, 1, 0], [1, 0, 0], [0, 0, 0]])
+
+        fa2 = ForceAtlas2(strongGravityMode=True, gravity=2.0, verbose=False)
+        pos = fa2.forceatlas2(G, iterations=50)
+
+        assert len(pos) == 3
+        # With strong gravity, nodes should be pulled toward origin
+
+    def test_layout_outbound_attraction_distribution(self):
+        """Test layout with outbound attraction distribution (hub dissuasion)"""
+        # Create a star graph (one hub connected to all others)
+        n = 5
+        G = np.zeros((n, n))
+        for i in range(1, n):
+            G[0, i] = 1
+            G[i, 0] = 1
+
+        fa2 = ForceAtlas2(outboundAttractionDistribution=True, verbose=False)
+        pos = fa2.forceatlas2(G, iterations=50)
+
+        assert len(pos) == n
+
+    def test_layout_edge_weight_influence(self):
+        """Test layout with different edge weight influence"""
+        G = np.array([[0.0, 1.0, 5.0], [1.0, 0.0, 1.0], [5.0, 1.0, 0.0]])
+
+        # With weight influence = 1, weights should matter
+        fa2_weighted = ForceAtlas2(edgeWeightInfluence=1.0, verbose=False)
+        pos_weighted = fa2_weighted.forceatlas2(G, iterations=50)
+
+        # With weight influence = 0, weights should be ignored
+        fa2_unweighted = ForceAtlas2(edgeWeightInfluence=0.0, verbose=False)
+        pos_unweighted = fa2_unweighted.forceatlas2(G, iterations=50)
+
+        assert len(pos_weighted) == 3
+        assert len(pos_unweighted) == 3
+
+    def test_layout_different_iterations(self):
+        """Test that more iterations produce different results"""
+        G = np.array([[0, 1, 1], [1, 0, 1], [1, 1, 0]])
+
+        fa2 = ForceAtlas2(verbose=False)
+
+        # Set seed for reproducibility
+        np.random.seed(42)
+        pos_10 = fa2.forceatlas2(G, iterations=10)
+
+        np.random.seed(42)
+        pos_100 = fa2.forceatlas2(G, iterations=100)
+
+        # Results should be different with different iteration counts
+        assert pos_10 != pos_100
+
+    def test_layout_sparse_matrix(self):
+        """Test layout with scipy sparse matrix"""
+        G = scipy.sparse.lil_matrix(
+            [[0, 1, 1, 0], [1, 0, 1, 0], [1, 1, 0, 1], [0, 0, 1, 0]]
+        )
+
+        fa2 = ForceAtlas2(verbose=False)
+        pos = fa2.forceatlas2(G, iterations=10)
+
+        assert len(pos) == 4
+
+    def test_layout_single_node(self):
+        """Test layout with single node"""
+        G = np.array([[0]])
+
+        fa2 = ForceAtlas2(verbose=False)
+        pos = fa2.forceatlas2(G, iterations=10)
+
+        assert len(pos) == 1
+
+    def test_layout_disconnected_nodes(self):
+        """Test layout with disconnected nodes"""
+        G = np.array([[0, 1, 0, 0], [1, 0, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0]])
+
+        fa2 = ForceAtlas2(verbose=False)
+        pos = fa2.forceatlas2(G, iterations=50)
+
+        assert len(pos) == 4
+        # Two disconnected components
+
+
+class TestForceAtlas2Parameters:
+    """Test different parameter combinations"""
+
+    def test_varying_scaling_ratio(self):
+        """Test with different scaling ratios"""
+        G = np.array([[0, 1, 1], [1, 0, 1], [1, 1, 0]])
+
+        for ratio in [0.5, 1.0, 2.0, 5.0]:
+            fa2 = ForceAtlas2(scalingRatio=ratio, verbose=False)
+            pos = fa2.forceatlas2(G, iterations=10)
+            assert len(pos) == 3
+
+    def test_varying_gravity(self):
+        """Test with different gravity values"""
+        G = np.array([[0, 1, 0], [1, 0, 1], [0, 1, 0]])
+
+        for grav in [0.1, 1.0, 5.0, 10.0]:
+            fa2 = ForceAtlas2(gravity=grav, verbose=False)
+            pos = fa2.forceatlas2(G, iterations=10)
+            assert len(pos) == 3
+
+    def test_varying_jitter_tolerance(self):
+        """Test with different jitter tolerance values"""
+        G = np.array([[0, 1, 1], [1, 0, 1], [1, 1, 0]])
+
+        for jitter in [0.5, 1.0, 2.0]:
+            fa2 = ForceAtlas2(jitterTolerance=jitter, verbose=False)
+            pos = fa2.forceatlas2(G, iterations=10)
+            assert len(pos) == 3
+
+    def test_varying_barnes_hut_theta(self):
+        """Test with different Barnes-Hut theta values"""
+        G = np.array([[0, 1, 1, 0], [1, 0, 1, 1], [1, 1, 0, 1], [0, 1, 1, 0]])
+
+        for theta in [0.5, 1.0, 1.2, 2.0]:
+            fa2 = ForceAtlas2(barnesHutTheta=theta, verbose=False)
+            pos = fa2.forceatlas2(G, iterations=10)
+            assert len(pos) == 4
+
+
+class TestForceAtlas2Reproducibility:
+    """Test reproducibility with same random seed"""
+
+    def test_reproducibility_with_seed(self):
+        """Test that same seed produces same results"""
+        G = np.array([[0, 1, 1], [1, 0, 1], [1, 1, 0]])
+
+        fa2 = ForceAtlas2(verbose=False)
+
+        # First run
+        np.random.seed(42)
+        import random
+
+        random.seed(42)
+        pos1 = fa2.forceatlas2(G, iterations=10)
+
+        # Second run with same seed
+        np.random.seed(42)
+        random.seed(42)
+        pos2 = fa2.forceatlas2(G, iterations=10)
+
+        # Results should be identical
+        for p1, p2 in zip(pos1, pos2):
+            assert abs(p1[0] - p2[0]) < 1e-10
+            assert abs(p1[1] - p2[1]) < 1e-10
+
+
+class TestForceAtlas2EdgeCases:
+    """Test edge cases and boundary conditions"""
+
+    def test_zero_iterations(self):
+        """Test with zero iterations"""
+        G = np.array([[0, 1], [1, 0]])
+
+        fa2 = ForceAtlas2(verbose=False)
+        pos = fa2.forceatlas2(G, iterations=0)
+
+        assert len(pos) == 2
+
+    def test_large_graph(self):
+        """Test with larger graph"""
+        n = 50
+        G = np.zeros((n, n))
+        # Create a ring graph
+        for i in range(n):
+            G[i, (i + 1) % n] = 1
+            G[(i + 1) % n, i] = 1
+
+        fa2 = ForceAtlas2(verbose=False)
+        pos = fa2.forceatlas2(G, iterations=10)
+
+        assert len(pos) == n
+
+    def test_complete_graph(self):
+        """Test with complete graph"""
+        n = 10
+        G = np.ones((n, n)) - np.eye(n)
+
+        fa2 = ForceAtlas2(verbose=False)
+        pos = fa2.forceatlas2(G, iterations=10)
+
+        assert len(pos) == n
+
+    def test_no_edges(self):
+        """Test with graph with no edges"""
+        G = np.zeros((5, 5))
+
+        fa2 = ForceAtlas2(verbose=False)
+        pos = fa2.forceatlas2(G, iterations=10)
+
+        assert len(pos) == 5
diff --git a/tests/test_igraph_integration.py b/tests/test_igraph_integration.py
new file mode 100644
index 0000000..343205a
--- /dev/null
+++ b/tests/test_igraph_integration.py
@@ -0,0 +1,300 @@
+"""
+Tests for igraph integration
+"""
+
+import numpy as np
+import pytest
+
+try:
+    import igraph
+
+    IGRAPH_AVAILABLE = True
+except ImportError:
+    IGRAPH_AVAILABLE = False
+
+from fa2_modified import ForceAtlas2
+
+
+@pytest.mark.skipif(not IGRAPH_AVAILABLE, reason="igraph not installed")
+class TestIgraphIntegration:
+    """Test ForceAtlas2 with igraph graphs"""
+
+    def test_igraph_simple_graph(self):
+        """Test with simple igraph graph"""
+        G = igraph.Graph()
+        G.add_vertices(4)
+        G.add_edges([(0, 1), (1, 2), (2, 3), (3, 0)])
+
+        fa2 = ForceAtlas2(verbose=False)
+        layout = fa2.forceatlas2_igraph_layout(G, iterations=10)
+
+        assert isinstance(layout, igraph.layout.Layout)
+        assert len(layout) == 4
+        assert layout.dim == 2
+
+    def test_igraph_with_initial_positions(self):
+        """Test igraph layout with initial positions"""
+        G = igraph.Graph()
+        G.add_vertices(3)
+        G.add_edges([(0, 1), (1, 2)])
+
+        initial_pos = [[0.0, 0.0], [1.0, 0.0], [2.0, 0.0]]
+
+        fa2 = ForceAtlas2(verbose=False)
+        layout = fa2.forceatlas2_igraph_layout(G, pos=initial_pos, iterations=10)
+
+        assert len(layout) == 3
+
+    def test_igraph_with_numpy_positions(self):
+        """Test igraph layout with numpy array positions"""
+        G = igraph.Graph()
+        G.add_vertices(3)
+        G.add_edges([(0, 1), (1, 2)])
+
+        initial_pos = np.array([[0.0, 0.0], [1.0, 0.0], [2.0, 0.0]])
+
+        fa2 = ForceAtlas2(verbose=False)
+        layout = fa2.forceatlas2_igraph_layout(G, pos=initial_pos, iterations=10)
+
+        assert len(layout) == 3
+
+    def test_igraph_weighted_graph(self):
+        """Test with weighted igraph graph"""
+        G = igraph.Graph()
+        G.add_vertices(3)
+        G.add_edges([(0, 1), (1, 2), (2, 0)])
+        G.es["weight"] = [1.0, 5.0, 1.0]
+
+        fa2 = ForceAtlas2(verbose=False)
+        layout = fa2.forceatlas2_igraph_layout(G, iterations=10, weight_attr="weight")
+
+        assert len(layout) == 3
+
+    def test_igraph_custom_weight_attribute(self):
+        """Test with custom weight attribute name"""
+        G = igraph.Graph()
+        G.add_vertices(3)
+        G.add_edges([(0, 1), (1, 2)])
+        G.es["strength"] = [2.0, 3.0]
+
+        fa2 = ForceAtlas2(verbose=False)
+        layout = fa2.forceatlas2_igraph_layout(G, iterations=10, weight_attr="strength")
+
+        assert len(layout) == 3
+
+    def test_igraph_star_graph(self):
+        """Test with star graph (hub topology)"""
+        G = igraph.Graph.Star(6)
+
+        fa2 = ForceAtlas2(outboundAttractionDistribution=True, verbose=False)
+        layout = fa2.forceatlas2_igraph_layout(G, iterations=50)
+
+        assert len(layout) == 6
+
+    def test_igraph_complete_graph(self):
+        """Test with complete graph"""
+        G = igraph.Graph.Full(5)
+
+        fa2 = ForceAtlas2(verbose=False)
+        layout = fa2.forceatlas2_igraph_layout(G, iterations=10)
+
+        assert len(layout) == 5
+
+    def test_igraph_ring_graph(self):
+        """Test with ring graph"""
+        G = igraph.Graph.Ring(8)
+
+        fa2 = ForceAtlas2(verbose=False)
+        layout = fa2.forceatlas2_igraph_layout(G, iterations=50)
+
+        assert len(layout) == 8
+
+    def test_igraph_tree_graph(self):
+        """Test with tree graph"""
+        G = igraph.Graph.Tree(15, 3)  # Tree with 15 nodes, 3 children per node
+
+        fa2 = ForceAtlas2(verbose=False)
+        layout = fa2.forceatlas2_igraph_layout(G, iterations=50)
+
+        assert len(layout) == 15
+
+    def test_igraph_lattice_graph(self):
+        """Test with lattice graph"""
+        G = igraph.Graph.Lattice([3, 3], circular=False)
+
+        fa2 = ForceAtlas2(verbose=False)
+        layout = fa2.forceatlas2_igraph_layout(G, iterations=50)
+
+        assert len(layout) == 9
+
+    def test_igraph_single_node(self):
+        """Test with single node"""
+        G = igraph.Graph()
+        G.add_vertices(1)
+
+        fa2 = ForceAtlas2(verbose=False)
+        layout = fa2.forceatlas2_igraph_layout(G, iterations=10)
+
+        assert len(layout) == 1
+
+    def test_igraph_disconnected_graph(self):
+        """Test with disconnected components"""
+        G = igraph.Graph()
+        G.add_vertices(6)
+        G.add_edges([(0, 1), (1, 2)])  # Component 1
+        G.add_edges([(3, 4), (4, 5)])  # Component 2
+
+        fa2 = ForceAtlas2(verbose=False)
+        layout = fa2.forceatlas2_igraph_layout(G, iterations=50)
+
+        assert len(layout) == 6
+
+    def test_igraph_no_edges(self):
+        """Test with nodes but no edges"""
+        G = igraph.Graph()
+        G.add_vertices(4)
+
+        fa2 = ForceAtlas2(verbose=False)
+        layout = fa2.forceatlas2_igraph_layout(G, iterations=10)
+
+        assert len(layout) == 4
+
+    def test_igraph_layout_coordinates_format(self):
+        """Test that layout coordinates are in correct format"""
+        G = igraph.Graph()
+        G.add_vertices(3)
+        G.add_edges([(0, 1), (1, 2)])
+
+        fa2 = ForceAtlas2(verbose=False)
+        layout = fa2.forceatlas2_igraph_layout(G, iterations=10)
+
+        coords = layout.coords
+        assert len(coords) == 3
+        for coord in coords:
+            assert len(coord) == 2
+            assert isinstance(coord[0], (int, float))
+            assert isinstance(coord[1], (int, float))
+
+    def test_igraph_parameter_combinations(self):
+        """Test various parameter combinations with igraph"""
+        G = igraph.Graph.Erdos_Renyi(n=30, p=0.1)
+
+        params = [
+            {"barnesHutOptimize": True, "strongGravityMode": False},
+            {"barnesHutOptimize": False, "strongGravityMode": True},
+            {"outboundAttractionDistribution": True, "edgeWeightInfluence": 0.5},
+            {"scalingRatio": 5.0, "gravity": 2.0},
+        ]
+
+        for param in params:
+            fa2 = ForceAtlas2(**param, verbose=False)
+            layout = fa2.forceatlas2_igraph_layout(G, iterations=10)
+            assert len(layout) == 30
+
+    def test_igraph_invalid_input(self):
+        """Test that non-igraph graph raises error"""
+        fa2 = ForceAtlas2(verbose=False)
+
+        with pytest.raises(AssertionError):
+            fa2.forceatlas2_igraph_layout([[0, 1], [1, 0]], iterations=10)
+
+    def test_igraph_invalid_pos_type(self):
+        """Test that invalid pos type raises error"""
+        G = igraph.Graph()
+        G.add_vertices(3)
+        G.add_edges([(0, 1), (1, 2)])
+
+        fa2 = ForceAtlas2(verbose=False)
+
+        with pytest.raises(AssertionError):
+            fa2.forceatlas2_igraph_layout(G, pos="invalid", iterations=10)
+
+    def test_igraph_directed_graph(self):
+        """Test with directed graph (converted to undirected internally)"""
+        G = igraph.Graph(directed=True)
+        G.add_vertices(3)
+        G.add_edges([(0, 1), (1, 2), (2, 0)])
+
+        fa2 = ForceAtlas2(verbose=False)
+        # Note: The current implementation assumes undirected graphs
+        # This test verifies the conversion happens correctly
+        layout = fa2.forceatlas2_igraph_layout(G, iterations=10)
+
+        assert len(layout) == 3
+
+    def test_igraph_random_graph(self):
+        """Test with random Erdos-Renyi graph"""
+        G = igraph.Graph.Erdos_Renyi(n=25, p=0.15)
+
+        fa2 = ForceAtlas2(verbose=False)
+        layout = fa2.forceatlas2_igraph_layout(G, iterations=50)
+
+        assert len(layout) == 25
+
+    def test_igraph_barabasi_albert_graph(self):
+        """Test with Barabasi-Albert (scale-free) graph"""
+        G = igraph.Graph.Barabasi(n=30, m=2)
+
+        fa2 = ForceAtlas2(outboundAttractionDistribution=True, verbose=False)
+        layout = fa2.forceatlas2_igraph_layout(G, iterations=50)
+
+        assert len(layout) == 30
+
+
+@pytest.mark.skipif(not IGRAPH_AVAILABLE, reason="igraph not installed")
+class TestIgraphReproducibility:
+    """Test reproducibility with igraph graphs"""
+
+    def test_reproducibility_with_seed(self):
+        """Test that same seed produces same igraph layout"""
+        G = igraph.Graph.Erdos_Renyi(n=20, p=0.2)
+
+        fa2 = ForceAtlas2(verbose=False)
+
+        # First run
+        np.random.seed(42)
+        import random
+
+        random.seed(42)
+        layout1 = fa2.forceatlas2_igraph_layout(G, iterations=50)
+
+        # Second run with same seed
+        np.random.seed(42)
+        random.seed(42)
+        layout2 = fa2.forceatlas2_igraph_layout(G, iterations=50)
+
+        # Results should be identical
+        coords1 = layout1.coords
+        coords2 = layout2.coords
+        for i in range(len(coords1)):
+            assert abs(coords1[i][0] - coords2[i][0]) < 1e-10
+            assert abs(coords1[i][1] - coords2[i][1]) < 1e-10
+
+
+@pytest.mark.skipif(not IGRAPH_AVAILABLE, reason="igraph not installed")
+class TestIgraphWeights:
+    """Test edge weight handling in igraph graphs"""
+
+    def test_no_weights(self):
+        """Test without weights"""
+        G = igraph.Graph()
+        G.add_vertices(4)
+        G.add_edges([(0, 1), (1, 2), (2, 3)])
+
+        fa2 = ForceAtlas2(verbose=False)
+        layout = fa2.forceatlas2_igraph_layout(G, iterations=10, weight_attr=None)
+
+        assert len(layout) == 4
+
+    def test_varying_weights(self):
+        """Test with varying edge weights"""
+        G = igraph.Graph()
+        G.add_vertices(4)
+        G.add_edges([(0, 1), (1, 2), (2, 3), (3, 0)])
+        G.es["weight"] = [1.0, 10.0, 1.0, 1.0]
+
+        fa2 = ForceAtlas2(edgeWeightInfluence=1.0, verbose=False)
+        layout = fa2.forceatlas2_igraph_layout(G, iterations=100, weight_attr="weight")
+
+        assert len(layout) == 4
+        # The edge with weight 10 should pull nodes 1 and 2 closer
diff --git a/tests/test_integration.py b/tests/test_integration.py
new file mode 100644
index 0000000..bce1464
--- /dev/null
+++ b/tests/test_integration.py
@@ -0,0 +1,347 @@
+"""
+Integration tests for fa2_modified package
+Tests the interaction between different components
+"""
+
+import numpy as np
+import pytest
+import scipy.sparse
+
+from fa2_modified import ForceAtlas2
+
+
+class TestEndToEndWorkflow:
+    """Test complete end-to-end workflows"""
+
+    def test_numpy_to_layout_workflow(self):
+        """Test complete workflow from numpy array to layout"""
+        # Create adjacency matrix
+        G = np.array(
+            [
+                [0, 1, 1, 0, 0],
+                [1, 0, 1, 1, 0],
+                [1, 1, 0, 1, 1],
+                [0, 1, 1, 0, 1],
+                [0, 0, 1, 1, 0],
+            ]
+        )
+
+        # Initialize ForceAtlas2
+        fa2 = ForceAtlas2(
+            barnesHutOptimize=True,
+            barnesHutTheta=1.2,
+            scalingRatio=2.0,
+            strongGravityMode=False,
+            gravity=1.0,
+            verbose=False,
+        )
+
+        # Compute layout
+        positions = fa2.forceatlas2(G, iterations=100)
+
+        # Verify results
+        assert len(positions) == 5
+        assert all(len(pos) == 2 for pos in positions)
+
+        # Verify positions are reasonable (not all at origin)
+        x_coords = [pos[0] for pos in positions]
+        y_coords = [pos[1] for pos in positions]
+        assert np.std(x_coords) > 0.01
+        assert np.std(y_coords) > 0.01
+
+    def test_sparse_to_layout_workflow(self):
+        """Test complete workflow from sparse matrix to layout"""
+        # Create sparse adjacency matrix
+        G = scipy.sparse.lil_matrix(
+            [
+                [0, 1, 1, 0, 0],
+                [1, 0, 1, 1, 0],
+                [1, 1, 0, 1, 1],
+                [0, 1, 1, 0, 1],
+                [0, 0, 1, 1, 0],
+            ]
+        )
+
+        fa2 = ForceAtlas2(verbose=False)
+        positions = fa2.forceatlas2(G, iterations=100)
+
+        assert len(positions) == 5
+        assert all(len(pos) == 2 for pos in positions)
+
+    def test_weighted_workflow(self):
+        """Test workflow with weighted edges"""
+        # Create weighted graph where one edge is much stronger
+        G = np.array(
+            [
+                [0.0, 1.0, 1.0, 0.0],
+                [1.0, 0.0, 10.0, 0.0],  # Strong connection between 1 and 2
+                [1.0, 10.0, 0.0, 1.0],
+                [0.0, 0.0, 1.0, 0.0],
+            ]
+        )
+
+        fa2 = ForceAtlas2(edgeWeightInfluence=1.0, verbose=False)
+        positions = fa2.forceatlas2(G, iterations=200)
+
+        # Nodes 1 and 2 should be closer together due to strong edge
+        dist_1_2 = np.linalg.norm(np.array(positions[1]) - np.array(positions[2]))
+        dist_0_1 = np.linalg.norm(np.array(positions[0]) - np.array(positions[1]))
+
+        # The strong edge should pull nodes closer (this is a heuristic test)
+        assert len(positions) == 4
+
+    def test_initial_positions_workflow(self):
+        """Test workflow with provided initial positions"""
+        G = np.array([[0, 1, 0], [1, 0, 1], [0, 1, 0]])
+
+        # Start with nodes in a line
+        initial_pos = np.array([[0.0, 0.0], [5.0, 0.0], [10.0, 0.0]])
+
+        fa2 = ForceAtlas2(verbose=False)
+        positions = fa2.forceatlas2(G, pos=initial_pos, iterations=100)
+
+        assert len(positions) == 3
+        # Positions should have changed from initial
+        for i in range(3):
+            pos_changed = (
+                positions[i][0] != initial_pos[i][0]
+                or positions[i][1] != initial_pos[i][1]
+            )
+            assert pos_changed
+
+
+class TestParameterInteractions:
+    """Test interactions between different parameters"""
+
+    def test_gravity_repulsion_balance(self):
+        """Test balance between gravity and repulsion"""
+        G = np.array([[0, 1, 0, 0], [1, 0, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0]])
+
+        # High gravity should keep nodes close to origin
+        fa2_high_gravity = ForceAtlas2(gravity=10.0, scalingRatio=1.0, verbose=False)
+        pos_high_gravity = fa2_high_gravity.forceatlas2(G, iterations=100)
+
+        # Low gravity should allow nodes to spread out more
+        fa2_low_gravity = ForceAtlas2(gravity=0.1, scalingRatio=1.0, verbose=False)
+        pos_low_gravity = fa2_low_gravity.forceatlas2(G, iterations=100)
+
+        # Calculate spread from origin
+        spread_high = np.mean([np.linalg.norm(pos) for pos in pos_high_gravity])
+        spread_low = np.mean([np.linalg.norm(pos) for pos in pos_low_gravity])
+
+        # Lower gravity should result in larger spread (generally)
+        assert len(pos_high_gravity) == 4
+        assert len(pos_low_gravity) == 4
+
+    def test_barnes_hut_vs_direct(self):
+        """Test that Barnes-Hut gives similar results to direct computation"""
+        G = np.array(
+            [
+                [0, 1, 1, 0, 0],
+                [1, 0, 1, 1, 0],
+                [1, 1, 0, 1, 1],
+                [0, 1, 1, 0, 1],
+                [0, 0, 1, 1, 0],
+            ]
+        )
+
+        # Set seed for reproducibility
+        np.random.seed(42)
+        import random
+
+        random.seed(42)
+
+        fa2_barnes = ForceAtlas2(barnesHutOptimize=True, verbose=False)
+        pos_barnes = fa2_barnes.forceatlas2(G, iterations=50)
+
+        np.random.seed(42)
+        random.seed(42)
+
+        fa2_direct = ForceAtlas2(barnesHutOptimize=False, verbose=False)
+        pos_direct = fa2_direct.forceatlas2(G, iterations=50)
+
+        # Results should be similar (not identical due to approximation)
+        assert len(pos_barnes) == len(pos_direct)
+
+    def test_strong_vs_normal_gravity(self):
+        """Test difference between strong and normal gravity"""
+        G = np.array([[0, 1, 0], [1, 0, 1], [0, 1, 0]])
+
+        np.random.seed(42)
+        fa2_strong = ForceAtlas2(strongGravityMode=True, gravity=1.0, verbose=False)
+        pos_strong = fa2_strong.forceatlas2(G, iterations=100)
+
+        np.random.seed(42)
+        fa2_normal = ForceAtlas2(strongGravityMode=False, gravity=1.0, verbose=False)
+        pos_normal = fa2_normal.forceatlas2(G, iterations=100)
+
+        assert len(pos_strong) == 3
+        assert len(pos_normal) == 3
+
+
+class TestScalability:
+    """Test with different graph sizes"""
+
+    def test_small_graph(self):
+        """Test with very small graph (2 nodes)"""
+        G = np.array([[0, 1], [1, 0]])
+
+        fa2 = ForceAtlas2(verbose=False)
+        positions = fa2.forceatlas2(G, iterations=10)
+
+        assert len(positions) == 2
+
+    def test_medium_graph(self):
+        """Test with medium-sized graph (50 nodes)"""
+        n = 50
+        G = np.zeros((n, n))
+        # Create a random graph with ~10% edge density
+        for i in range(n):
+            for j in range(i + 1, n):
+                if np.random.random() < 0.1:
+                    G[i, j] = 1
+                    G[j, i] = 1
+
+        fa2 = ForceAtlas2(barnesHutOptimize=True, verbose=False)
+        positions = fa2.forceatlas2(G, iterations=50)
+
+        assert len(positions) == n
+
+    @pytest.mark.slow
+    def test_large_graph(self):
+        """Test with larger graph (200 nodes)"""
+        n = 200
+        G = scipy.sparse.lil_matrix((n, n))
+        # Create a sparse random graph
+        for i in range(n):
+            num_edges = np.random.randint(1, 5)
+            for _ in range(num_edges):
+                j = np.random.randint(0, n)
+                if i != j:
+                    G[i, j] = 1
+                    G[j, i] = 1
+
+        fa2 = ForceAtlas2(barnesHutOptimize=True, verbose=False)
+        positions = fa2.forceatlas2(G, iterations=50)
+
+        assert len(positions) == n
+
+
+class TestSpecialGraphTopologies:
+    """Test with special graph structures"""
+
+    def test_star_topology(self):
+        """Test star graph (hub and spokes)"""
+        n = 10
+        G = np.zeros((n, n))
+        # Node 0 is hub, connected to all others
+        for i in range(1, n):
+            G[0, i] = 1
+            G[i, 0] = 1
+
+        fa2 = ForceAtlas2(outboundAttractionDistribution=True, verbose=False)
+        positions = fa2.forceatlas2(G, iterations=100)
+
+        assert len(positions) == n
+
+    def test_complete_bipartite(self):
+        """Test complete bipartite graph"""
+        n1, n2 = 5, 5
+        n = n1 + n2
+        G = np.zeros((n, n))
+        # Connect all nodes in first set to all nodes in second set
+        for i in range(n1):
+            for j in range(n1, n):
+                G[i, j] = 1
+                G[j, i] = 1
+
+        fa2 = ForceAtlas2(verbose=False)
+        positions = fa2.forceatlas2(G, iterations=100)
+
+        assert len(positions) == n
+
+    def test_chain_topology(self):
+        """Test chain/path graph"""
+        n = 10
+        G = np.zeros((n, n))
+        for i in range(n - 1):
+            G[i, i + 1] = 1
+            G[i + 1, i] = 1
+
+        fa2 = ForceAtlas2(verbose=False)
+        positions = fa2.forceatlas2(G, iterations=100)
+
+        assert len(positions) == n
+
+    def test_grid_topology(self):
+        """Test grid graph"""
+        rows, cols = 4, 4
+        n = rows * cols
+        G = np.zeros((n, n))
+
+        # Helper to convert 2D coords to 1D index
+        def idx(r, c):
+            return r * cols + c
+
+        # Connect grid neighbors
+        for r in range(rows):
+            for c in range(cols):
+                if r < rows - 1:
+                    G[idx(r, c), idx(r + 1, c)] = 1
+                    G[idx(r + 1, c), idx(r, c)] = 1
+                if c < cols - 1:
+                    G[idx(r, c), idx(r, c + 1)] = 1
+                    G[idx(r, c + 1), idx(r, c)] = 1
+
+        fa2 = ForceAtlas2(verbose=False)
+        positions = fa2.forceatlas2(G, iterations=100)
+
+        assert len(positions) == n
+
+
+class TestConvergence:
+    """Test convergence properties"""
+
+    def test_convergence_over_iterations(self):
+        """Test that layout stabilizes over iterations"""
+        G = np.array([[0, 1, 1, 0], [1, 0, 1, 1], [1, 1, 0, 1], [0, 1, 1, 0]])
+
+        fa2 = ForceAtlas2(verbose=False)
+
+        # Run with different iteration counts
+        np.random.seed(42)
+        pos_10 = fa2.forceatlas2(G, iterations=10)
+
+        np.random.seed(42)
+        pos_100 = fa2.forceatlas2(G, iterations=100)
+
+        np.random.seed(42)
+        pos_1000 = fa2.forceatlas2(G, iterations=1000)
+
+        # More iterations should produce valid layouts
+        assert len(pos_10) == 4
+        assert len(pos_100) == 4
+        assert len(pos_1000) == 4
+
+    def test_reproducibility_same_seed(self):
+        """Test that same seed produces identical results"""
+        G = np.array([[0, 1, 1], [1, 0, 1], [1, 1, 0]])
+
+        fa2 = ForceAtlas2(verbose=False)
+
+        # First run
+        np.random.seed(12345)
+        import random
+
+        random.seed(12345)
+        pos1 = fa2.forceatlas2(G, iterations=50)
+
+        # Second run with same seed
+        np.random.seed(12345)
+        random.seed(12345)
+        pos2 = fa2.forceatlas2(G, iterations=50)
+
+        # Should be identical
+        for i in range(len(pos1)):
+            assert abs(pos1[i][0] - pos2[i][0]) < 1e-10
+            assert abs(pos1[i][1] - pos2[i][1]) < 1e-10
diff --git a/tests/test_networkx_integration.py b/tests/test_networkx_integration.py
new file mode 100644
index 0000000..835d895
--- /dev/null
+++ b/tests/test_networkx_integration.py
@@ -0,0 +1,289 @@
+"""
+Tests for NetworkX integration
+"""
+
+import numbers
+
+import numpy as np
+import pytest
+
+try:
+    import networkx as nx
+
+    NETWORKX_AVAILABLE = True
+except ImportError:
+    NETWORKX_AVAILABLE = False
+
+from fa2_modified import ForceAtlas2
+
+
+@pytest.mark.skipif(not NETWORKX_AVAILABLE, reason="NetworkX not installed")
+class TestNetworkXIntegration:
+    """Test ForceAtlas2 with NetworkX graphs"""
+
+    def test_networkx_simple_graph(self):
+        """Test with simple NetworkX graph"""
+        G = nx.Graph()
+        G.add_edges_from([(0, 1), (1, 2), (2, 3), (3, 0)])
+
+        fa2 = ForceAtlas2(verbose=False)
+        pos = fa2.forceatlas2_networkx_layout(G, iterations=10)
+
+        assert len(pos) == 4
+        assert all(node in pos for node in G.nodes())
+        assert all(len(pos[node]) == 2 for node in G.nodes())
+
+    def test_networkx_with_initial_positions(self):
+        """Test NetworkX layout with initial positions"""
+        G = nx.Graph()
+        G.add_edges_from([(0, 1), (1, 2)])
+
+        initial_pos = {0: (0.0, 0.0), 1: (1.0, 0.0), 2: (2.0, 0.0)}
+
+        fa2 = ForceAtlas2(verbose=False)
+        pos = fa2.forceatlas2_networkx_layout(G, pos=initial_pos, iterations=10)
+
+        assert len(pos) == 3
+        assert all(node in pos for node in G.nodes())
+
+    def test_networkx_weighted_graph(self):
+        """Test with weighted NetworkX graph"""
+        G = nx.Graph()
+        G.add_edge(0, 1, weight=1.0)
+        G.add_edge(1, 2, weight=5.0)
+        G.add_edge(2, 0, weight=1.0)
+
+        fa2 = ForceAtlas2(verbose=False)
+        pos = fa2.forceatlas2_networkx_layout(G, iterations=10, weight_attr="weight")
+
+        assert len(pos) == 3
+
+    def test_networkx_star_graph(self):
+        """Test with star graph (hub topology)"""
+        G = nx.star_graph(5)
+
+        fa2 = ForceAtlas2(outboundAttractionDistribution=True, verbose=False)
+        pos = fa2.forceatlas2_networkx_layout(G, iterations=50)
+
+        assert len(pos) == 6  # 1 hub + 5 leaves
+
+    def test_networkx_complete_graph(self):
+        """Test with complete graph"""
+        G = nx.complete_graph(5)
+
+        fa2 = ForceAtlas2(verbose=False)
+        pos = fa2.forceatlas2_networkx_layout(G, iterations=10)
+
+        assert len(pos) == 5
+
+    def test_networkx_path_graph(self):
+        """Test with path graph"""
+        G = nx.path_graph(10)
+
+        fa2 = ForceAtlas2(verbose=False)
+        pos = fa2.forceatlas2_networkx_layout(G, iterations=50)
+
+        assert len(pos) == 10
+
+    def test_networkx_cycle_graph(self):
+        """Test with cycle graph"""
+        G = nx.cycle_graph(8)
+
+        fa2 = ForceAtlas2(verbose=False)
+        pos = fa2.forceatlas2_networkx_layout(G, iterations=50)
+
+        assert len(pos) == 8
+
+    def test_networkx_grid_graph(self):
+        """Test with 2D grid graph"""
+        G = nx.grid_2d_graph(3, 3)
+
+        fa2 = ForceAtlas2(verbose=False)
+        pos = fa2.forceatlas2_networkx_layout(G, iterations=50)
+
+        assert len(pos) == 9
+
+    def test_networkx_random_geometric_graph(self):
+        """Test with random geometric graph"""
+        G = nx.random_geometric_graph(20, 0.3)
+
+        fa2 = ForceAtlas2(verbose=False)
+        pos = fa2.forceatlas2_networkx_layout(G, iterations=50)
+
+        assert len(pos) == 20
+
+    def test_networkx_karate_club(self):
+        """Test with Karate Club graph (classic test graph)"""
+        G = nx.karate_club_graph()
+
+        fa2 = ForceAtlas2(verbose=False)
+        pos = fa2.forceatlas2_networkx_layout(G, iterations=100)
+
+        assert len(pos) == 34
+
+    def test_networkx_single_node(self):
+        """Test with single node"""
+        G = nx.Graph()
+        G.add_node(0)
+
+        fa2 = ForceAtlas2(verbose=False)
+        pos = fa2.forceatlas2_networkx_layout(G, iterations=10)
+
+        assert len(pos) == 1
+        assert 0 in pos
+
+    def test_networkx_disconnected_graph(self):
+        """Test with disconnected components"""
+        G = nx.Graph()
+        G.add_edges_from([(0, 1), (1, 2)])
+        G.add_edges_from([(3, 4), (4, 5)])
+
+        fa2 = ForceAtlas2(verbose=False)
+        pos = fa2.forceatlas2_networkx_layout(G, iterations=50)
+
+        assert len(pos) == 6
+
+    def test_networkx_string_node_labels(self):
+        """Test with string node labels"""
+        G = nx.Graph()
+        G.add_edges_from([("A", "B"), ("B", "C"), ("C", "A")])
+
+        fa2 = ForceAtlas2(verbose=False)
+        pos = fa2.forceatlas2_networkx_layout(G, iterations=10)
+
+        assert len(pos) == 3
+        assert "A" in pos
+        assert "B" in pos
+        assert "C" in pos
+
+    def test_networkx_mixed_node_labels(self):
+        """Test with mixed node labels"""
+        G = nx.Graph()
+        G.add_edges_from([(1, "A"), ("A", 2.5), (2.5, 1)])
+
+        fa2 = ForceAtlas2(verbose=False)
+        pos = fa2.forceatlas2_networkx_layout(G, iterations=10)
+
+        assert len(pos) == 3
+
+    def test_networkx_no_edges(self):
+        """Test with nodes but no edges"""
+        G = nx.Graph()
+        G.add_nodes_from([0, 1, 2, 3])
+
+        fa2 = ForceAtlas2(verbose=False)
+        pos = fa2.forceatlas2_networkx_layout(G, iterations=10)
+
+        assert len(pos) == 4
+
+    def test_networkx_position_format(self):
+        """Test that positions are in correct format"""
+        G = nx.Graph()
+        G.add_edges_from([(0, 1), (1, 2)])
+
+        fa2 = ForceAtlas2(verbose=False)
+        pos = fa2.forceatlas2_networkx_layout(G, iterations=10)
+
+        assert isinstance(pos, dict)
+        for node, position in pos.items():
+            assert len(position) == 2
+            assert isinstance(position[0], numbers.Real)
+            assert isinstance(position[1], numbers.Real)
+
+    def test_networkx_parameter_combinations(self):
+        """Test various parameter combinations with NetworkX"""
+        G = nx.karate_club_graph()
+
+        params = [
+            {"barnesHutOptimize": True, "strongGravityMode": False},
+            {"barnesHutOptimize": False, "strongGravityMode": True},
+            {"outboundAttractionDistribution": True, "edgeWeightInfluence": 0.5},
+            {"scalingRatio": 5.0, "gravity": 2.0},
+        ]
+
+        for param in params:
+            fa2 = ForceAtlas2(**param, verbose=False)
+            pos = fa2.forceatlas2_networkx_layout(G, iterations=10)
+            assert len(pos) == 34
+
+    def test_networkx_invalid_input(self):
+        """Test that non-NetworkX graph raises error"""
+        fa2 = ForceAtlas2(verbose=False)
+
+        with pytest.raises(AssertionError):
+            fa2.forceatlas2_networkx_layout([[0, 1], [1, 0]], iterations=10)
+
+    def test_networkx_invalid_pos_type(self):
+        """Test that invalid pos type raises error"""
+        G = nx.Graph()
+        G.add_edges_from([(0, 1), (1, 2)])
+
+        fa2 = ForceAtlas2(verbose=False)
+
+        with pytest.raises(AssertionError):
+            fa2.forceatlas2_networkx_layout(G, pos=[(0, 0), (1, 1)], iterations=10)
+
+
+@pytest.mark.skipif(not NETWORKX_AVAILABLE, reason="NetworkX not installed")
+class TestNetworkXEdgeWeights:
+    """Test edge weight handling in NetworkX graphs"""
+
+    def test_default_weight_attribute(self):
+        """Test with default weight attribute"""
+        G = nx.Graph()
+        G.add_edge(0, 1, weight=2.0)
+        G.add_edge(1, 2, weight=1.0)
+
+        fa2 = ForceAtlas2(edgeWeightInfluence=1.0, verbose=False)
+        pos = fa2.forceatlas2_networkx_layout(G, iterations=50)
+
+        assert len(pos) == 3
+
+    def test_custom_weight_attribute(self):
+        """Test with custom weight attribute"""
+        G = nx.Graph()
+        G.add_edge(0, 1, strength=3.0)
+        G.add_edge(1, 2, strength=1.0)
+
+        fa2 = ForceAtlas2(edgeWeightInfluence=1.0, verbose=False)
+        pos = fa2.forceatlas2_networkx_layout(G, iterations=50, weight_attr="strength")
+
+        assert len(pos) == 3
+
+    def test_no_weight_attribute(self):
+        """Test with no weight attribute (unweighted)"""
+        G = nx.Graph()
+        G.add_edges_from([(0, 1), (1, 2), (2, 0)])
+
+        fa2 = ForceAtlas2(verbose=False)
+        pos = fa2.forceatlas2_networkx_layout(G, iterations=50, weight_attr=None)
+
+        assert len(pos) == 3
+
+
+@pytest.mark.skipif(not NETWORKX_AVAILABLE, reason="NetworkX not installed")
+class TestNetworkXReproducibility:
+    """Test reproducibility with NetworkX graphs"""
+
+    def test_reproducibility_with_seed(self):
+        """Test that same seed produces same NetworkX layout"""
+        G = nx.karate_club_graph()
+
+        fa2 = ForceAtlas2(verbose=False)
+
+        # First run
+        np.random.seed(42)
+        import random
+
+        random.seed(42)
+        pos1 = fa2.forceatlas2_networkx_layout(G, iterations=50)
+
+        # Second run with same seed
+        np.random.seed(42)
+        random.seed(42)
+        pos2 = fa2.forceatlas2_networkx_layout(G, iterations=50)
+
+        # Results should be identical
+        for node in G.nodes():
+            assert abs(pos1[node][0] - pos2[node][0]) < 1e-10
+            assert abs(pos1[node][1] - pos2[node][1]) < 1e-10
diff --git a/tests/test_utils.py b/tests/test_utils.py
new file mode 100644
index 0000000..d2330c3
--- /dev/null
+++ b/tests/test_utils.py
@@ -0,0 +1,78 @@
+"""
+Tests for utility classes and functions
+"""
+
+import time
+
+import pytest
+
+from fa2_modified.forceatlas2 import Timer
+
+
+class TestTimer:
+    """Test Timer utility class"""
+
+    def test_timer_initialization(self):
+        """Test Timer initialization"""
+        timer = Timer()
+        assert timer.name == "Timer"
+        assert timer.start_time == 0.0
+        assert timer.total_time == 0.0
+
+    def test_timer_custom_name(self):
+        """Test Timer with custom name"""
+        timer = Timer(name="TestTimer")
+        assert timer.name == "TestTimer"
+
+    def test_timer_start_stop(self):
+        """Test timer start and stop"""
+        timer = Timer()
+        timer.start()
+        time.sleep(0.01)  # Sleep for 10ms
+        timer.stop()
+
+        assert timer.total_time > 0.0
+        assert timer.total_time < 1.0  # Should be much less than 1 second
+
+    def test_timer_multiple_runs(self):
+        """Test timer accumulates time over multiple runs"""
+        timer = Timer()
+
+        timer.start()
+        time.sleep(0.01)
+        timer.stop()
+        first_time = timer.total_time
+
+        timer.start()
+        time.sleep(0.01)
+        timer.stop()
+        second_time = timer.total_time
+
+        assert second_time > first_time
+        # Allow for OS scheduling/timer resolution differences across platforms
+        assert second_time >= first_time * 1.05
+
+    def test_timer_display(self, capsys):
+        """Test timer display output"""
+        timer = Timer(name="DisplayTest")
+        timer.start()
+        time.sleep(0.01)
+        timer.stop()
+
+        timer.display()
+
+        captured = capsys.readouterr()
+        assert "DisplayTest" in captured.out
+        assert "seconds" in captured.out
+
+    def test_timer_zero_time(self):
+        """Test timer with no elapsed time"""
+        timer = Timer()
+        assert timer.total_time == 0.0
+
+        # Start and immediately stop
+        timer.start()
+        timer.stop()
+
+        # Time should still be very small (close to zero)
+        assert timer.total_time >= 0.0