TauFactor

https://taufactor.readthedocs.io/en/latest/
Admin
TauFactor is a Python library for computing tortuosity and transport metrics in porous electrode...
Tokens:177,566
Snippets:2,057
Trust Score:9.6
Update:5 months ago
Show doc for...
Tokens:
Raw
### Install Dependencies and Run Tests

Source: https://taufactor.readthedocs.io/en/latest/contributing.html

Install TauFactor and its development dependencies in editable mode. Then, run linters and tests to ensure code quality and correctness.

```bash
pip install -e .[dev]
ruff check .
pytest
```

--------------------------------

### Example Usage of Tortuosity Calculation

Source: https://taufactor.readthedocs.io/en/latest/_sources/notebooks/01-tortuosity.ipynb.txt

Demonstrates how to use the tortuosity calculation functions with sample data. This example assumes you have already obtained the path length and straight-line distance.

```python
path_len = 150.0  # Example path length
straight_dist = 100.0  # Example straight-line distance

tortuosity = tortuosity_direct(path_len, straight_dist)
print(f"Tortuosity: {tortuosity}")
```

--------------------------------

### Initialize SORSolver with PyTorch

Source: https://taufactor.readthedocs.io/en/latest/_modules/taufactor/taufactor.html

Initializes the SORSolver with image data, over-relaxation factor, and device settings. Requires PyTorch to be installed. Computes initial field, conductive neighbors, and optionally reactive neighbors.

```python
import numpy as np
from abc import ABC, abstractmethod
from IPython.display import clear_output
from timeit import default_timer as timer
import matplotlib.pyplot as plt
try:
    import torch
except Exception:
    torch = None
import warnings
from .metrics import extract_through_feature




class SORSolver(ABC):
    """
    A minimal, clean template for SOR solvers.
    Subclasses override a few well-defined hooks.
    Args:
            img: labelled input image defining (non-)conducting phases.
            oemga: Over-relaxation factor for SOR scheme.
            device: The device to perform computations ('cpu' or 'cuda').
    """



    def __init__(self, img: np.ndarray, omega: float | None = None, precision=None, device='cuda'):
        if torch is None:
            raise ImportError(
                "PyTorch is required to use TauFactor solvers. Install pytorch following "
                "https://taufactor.readthedocs.io/en/latest/installation.html"
            )
        self.cpu_img = self._expand_to_4d(img)
        self.batch_size, self.Nx, self.Ny, self.Nz = self.cpu_img.shape
        self.device = self._init_device(device)
        self.precision = precision or torch.float

        # Overrelaxation factor for SOR
        if omega is None:
            omega = 2 - torch.pi / (1.5 * self.Nx)

        # Initialise pytorch tensors
        torch_img = torch.tensor(self.cpu_img, dtype=self.precision, device=self.device)
        mask = self.return_mask(torch_img)
        vol_x = torch.mean(mask, (2, 3)) # volume fraction
        self.field = self.init_field(mask)
        self.factor = self.init_conductive_neighbours(torch_img)

        # Optional for electrode simulations
        reac_nn = self.init_reactive_neighbours(torch_img)
        if reac_nn is not None:
            a_x = (torch.sum(reac_nn, (2, 3)) / (self.Ny*self.Nz*self.dx)) # surface area
            # Pre-compute reaction prefactor
            k_0 = torch.mean(vol_x, 1) / torch.mean(a_x*self.dx, 1) / self.Nx**2
            reac_nn = reac_nn * k_0[:, None, None, None]
            self.factor += reac_nn
            self.factor[self.factor == 0] = torch.inf
            self.a_x = a_x.cpu().numpy()
            self.k_0 = k_0.cpu().numpy()

        self.vol_x = vol_x.cpu().numpy()
        self.cb = self._init_chequerboard(omega)

        # Init params
        self.converged = False
        self.old_tau = 0
        self.iter = 0
        self.tau = None
        self.tau_x = None
        self.D_eff = None

```

--------------------------------

### Example Usage of Tortuosity Calculation

Source: https://taufactor.readthedocs.io/en/latest/_sources/notebooks/01-tortuosity.ipynb.txt

Demonstrates how to use the `calculate_tortuosity` function with sample values for path length and straight distance.

```python
path_len = 15.0  # Example path length
straight_dist = 10.0  # Example straight distance
tortuosity = calculate_tortuosity(path_len, straight_dist)
print(f"The tortuosity is: {tortuosity}")
```

--------------------------------

### Install TauFactor from Source (Development)

Source: https://taufactor.readthedocs.io/en/latest/_sources/installation.md.txt

Clone the development branch of TauFactor and install it locally. Adjust environment.yml for your CUDA version.

```bash
git clone https://github.com/tldr-group/taufactor.git
cd taufactor
git checkout development
```

```bash
conda env create -f environment.yml
conda activate taufactor
pip install -e .
```

--------------------------------

### Install PyTorch with CUDA

Source: https://taufactor.readthedocs.io/en/latest/_sources/installation.md.txt

Use this command to install PyTorch with CUDA support for Linux machines. Ensure your PyTorch version is compatible.

```bash
conda install pytorch pytorch-cuda=11.7 -c pytorch -c nvidia
```

--------------------------------

### Set Up TauFactor Development Environment with Conda

Source: https://taufactor.readthedocs.io/en/latest/installation.html

Create and activate a Conda environment using the provided environment file, then install TauFactor in editable mode. Adjust environment.yml for your CUDA version.

```bash
conda env create -f environment.yml
conda activate taufactor
pip install -e .

```

--------------------------------

### Install PyTorch with Conda for Linux/CUDA

Source: https://taufactor.readthedocs.io/en/latest/installation.html

Use this command to install PyTorch with CUDA support on a Linux machine. Ensure your PyTorch version meets the requirements.

```bash
conda install pytorch pytorch-cuda=11.7 -c pytorch -c nvidia

```

--------------------------------

### Example Usage of Tortuosity Calculation

Source: https://taufactor.readthedocs.io/en/latest/_sources/notebooks/01-tortuosity.ipynb.txt

Demonstrates how to use the `calculate_tortuosity` function with sample values. The path length must be greater than or equal to the straight-line distance.

```python
path = 10.0  # Example path length
straight = 5.0  # Example straight-line distance
tortuosity = calculate_tortuosity(path, straight)
print(f"The tortuosity is: {tortuosity}")
```

--------------------------------

### Install TauFactor Stable Release

Source: https://taufactor.readthedocs.io/en/latest/_sources/installation.md.txt

Install the most recent stable release of TauFactor using pip. This is the recommended installation method.

```bash
pip install taufactor
```

--------------------------------

### Initialize ElectrodeSolver Parameters

Source: https://taufactor.readthedocs.io/en/latest/_modules/taufactor/electrode.html

Initializes the ElectrodeSolver with image data, boundary conditions, and device settings. Sets up parameters for simulation.

```python
class ElectrodeSolver(SORSolver):
    """
    Default solver for two phase images. Once solve method is
    called, tau, D_eff and D_rel are available as attributes.
    """


    def __init__(self, img, conductive_label=1, reactive_label=0, \
                 omega=None, spacing=None, device='cuda'):
        """
        Initialise parameters, solution field and other tools that can be re-used
        for multiple solves.

        :param img: input img with 1s conductive and 0s non-conductive
        :param bc: Upper and lower boundary conditions. Leave as default.
        :param D_0: reference material diffusivity
        :param device: pytorch device, can be cuda or cpu 
        """
        self.left_bc = 1.0
        self.electrode_bc = 0.0
        self.cond_label=conductive_label
        self.reac_label=reactive_label
        self.dx = spacing or 1
        super().__init__(img, omega=omega, device=device)
        self.c_x = 0
```

--------------------------------

### Initialize and Solve with Periodic Boundary Conditions

Source: https://taufactor.readthedocs.io/en/latest/notebooks/00-overview.html

Initializes and solves using the PeriodicSolver, which applies periodic boundary conditions instead of the default mirror boundaries.

```python
s = tau.PeriodicSolver(img==0)
s.solve()
```

--------------------------------

### Initialize and Run Tortuosity Solver

Source: https://taufactor.readthedocs.io/en/latest/notebooks/01-tortuosity.html

Initializes the PeriodicElectrodeSolver with the batch of structures and solves for tortuosity. Sets a convergence criterion and iteration limit.

```python
s = tau.PeriodicElectrodeSolver(batch, device='cuda', spacing=Lx/nx)
tau_pore = s.solve(iter_limit=50000, conv_crit=1e-3)
```

--------------------------------

### Simulate and Save Impedance (Smooth-Rough)

Source: https://taufactor.readthedocs.io/en/latest/notebooks/10-tau-x-validation.html

Initializes and solves a PeriodicImpedanceSolver for the 'smooth-rough' configuration, then saves the resulting object.

```python
e = tau.PeriodicImpedanceSolver(fields["smooth-rough"], mode='nyquist')
e.solve(iter_limit=500000, conv_crit=-1, verbose='plot', plot_interval=100)
pkl_path = save_solver_object(e, "impedance_solve_smooth-rough")
```

--------------------------------

### Initialize Device

Source: https://taufactor.readthedocs.io/en/latest/_modules/taufactor/taufactor.html

Initializes and returns a PyTorch device. Defaults to CPU if CUDA is requested but not available, issuing a warning.

```python
    @staticmethod
    def _init_device(device) -> torch.device:
        # check device is available
        if torch.device(device).type.startswith('cuda') and not torch.cuda.is_available():
            device = torch.device('cpu')
            warnings.warn("CUDA not available, defaulting device to cpu. "
                          "To avoid this warning, explicitly set the device when "
                          "initialising the solver with device='cpu' ")
        else:
            device = torch.device(device)
        return device
```

--------------------------------

### Simulate and Save Impedance (Wedge)

Source: https://taufactor.readthedocs.io/en/latest/notebooks/10-tau-x-validation.html

Initializes and solves a PeriodicImpedanceSolver for the 'wedge' configuration with an increased iteration limit, then saves the resulting object.

```python
f = tau.PeriodicImpedanceSolver(fields["wedge"], mode='nyquist')
f.solve(iter_limit=1000000, conv_crit=-1, verbose='plot', plot_interval=100)
pkl_path = save_solver_object(f, "impedance_solve_wedge")
```

--------------------------------

### Calculate Tortuosity with `tortuosity` library

Source: https://taufactor.readthedocs.io/en/latest/_sources/notebooks/01-tortuosity.ipynb.txt

Use the `tortuosity` library to calculate the tortuosity of a path. Ensure the library is installed.

```python
import tortuosity
import numpy as np

# Example path (list of 3D coordinates)
path = np.array([
    [0, 0, 0],
    [1, 1, 1],
    [2, 0, 2],
    [3, 1, 3],
    [4, 0, 4]
])

tortuosity_value = tortuosity.tortuosity(path)
print(f"Tortuosity: {tortuosity_value}")
```

--------------------------------

### Initialize Field

Source: https://taufactor.readthedocs.io/en/latest/_modules/taufactor/electrode.html

Placeholder method for initializing the field. Returns None.

```python
def init_field(self, img):
        return None
```

--------------------------------

### Calculate Tortuosity

Source: https://taufactor.readthedocs.io/en/latest/_sources/notebooks/01-tortuosity.ipynb.txt

This snippet calculates the tortuosity of a path. It requires the path length and the straight-line distance between the start and end points.

```python
def calculate_tortuosity(path_length, straight_distance):
    """Calculates tortuosity.

    Args:
        path_length (float): The actual length of the path.
        straight_distance (float): The straight-line distance between the start and end points.

    Returns:
        float: The tortuosity value.
    """
    if straight_distance == 0:
        return float('inf')  # Avoid division by zero if start and end points are the same
    return path_length / straight_distance
```

--------------------------------

### Calculate Interfacial Areas of a Sphere

Source: https://taufactor.readthedocs.io/en/latest/_sources/notebooks/02-surface-areas.ipynb.txt

Computes the interfacial areas of a sphere. Set `normalize=True` to get specific interfacial areas in [1/m].

```python
print(tau.interfacial_areas(sphere, normalize=True))
```

--------------------------------

### Initialize Conductive Neighbors (Placeholder)

Source: https://taufactor.readthedocs.io/en/latest/_modules/taufactor/electrode.html

Placeholder method for initializing conductive neighbors. Returns None.

```python
def init_conductive_neighbours(self, img):
        return None
```

--------------------------------

### Initialize Linear Field

Source: https://taufactor.readthedocs.io/en/latest/_modules/taufactor/taufactor.html

Sets an initial linear field across the volume, scaled by a mask. This is used to establish a starting point for the simulation, typically with boundary conditions applied.

```python
[docs]
    def init_field(self, mask):
        """Sets an initial linear field across the volume"""
        sh = 1 / (2 * self.Nx)
        vec = torch.linspace(self.top_bc + sh, self.bot_bc - sh, self.Nx,
                             dtype=self.precision, device=self.device)
        for i in range(2):
            vec = torch.unsqueeze(vec, -1)
        vec = torch.unsqueeze(vec, 0)
        vec = vec.repeat(self.batch_size, 1, self.Ny, self.Nz, )
        return self._pad(mask * vec, [2*self.top_bc, 2*self.bot_bc])
```

--------------------------------

### Initialize and Solve Periodic Electrode Solver

Source: https://taufactor.readthedocs.io/en/latest/_sources/notebooks/01-tortuosity.ipynb.txt

Initializes the PeriodicElectrodeSolver with specified batch, device, and spacing, then solves for the field with given iteration limits, convergence criteria, and verbose settings. Use 'cuda' for GPU acceleration.

```python
s = tau.PeriodicElectrodeSolver(batch, device='cuda', spacing=Lx/nx)
 tau_pore = s.solve(iter_limit=50000, conv_crit=1e-3, verbose='debug', plot_interval=1)
```

--------------------------------

### Permuting Image Axes for Transport Direction

Source: https://taufactor.readthedocs.io/en/latest/_sources/notebooks/00-overview.ipynb.txt

Before solving, ensure the transport direction aligns with the desired image axis. Use the `flux_direction` utility to visualize and guide this permutation if necessary.

```python
img = tau.util.load_image("../data/phantom.png")
img = tau.util.permute(img, (1, 0, 2))
flux_direction(img)
```

--------------------------------

### Initialize Reactive Neighbors

Source: https://taufactor.readthedocs.io/en/latest/_modules/taufactor/electrode.html

Initializes reactive neighbors by identifying reactive labels in the image. Pads the image, sums by rolling, and filters based on conductive labels. Requires torch.

```python
def init_reactive_neighbours(self, img):
        img2 = torch.zeros_like(img)
        img2[img==self.reac_label] = 1
        img2 = self._pad(img2)[:, :, 1:-1, 1:-1]
        reac_nn = self._sum_by_rolling(img2)
        reac_nn = reac_nn[:, 1:-1]
        reac_nn[img != self.cond_label] = 0.0
        return reac_nn
```

--------------------------------

### Calculate Tortuosity using Euclidean Distance

Source: https://taufactor.readthedocs.io/en/latest/_sources/notebooks/01-tortuosity.ipynb.txt

Calculate tortuosity by dividing the total path length by the straight-line distance between the start and end points. This method assumes a simple, continuous path.

```python
import numpy as np

def calculate_tortuosity(path):
    """Calculates tortuosity based on path length and Euclidean distance."""
    path_length = np.sum(np.linalg.norm(np.diff(path, axis=0), axis=1))
    euclidean_distance = np.linalg.norm(path[-1] - path[0])
    if euclidean_distance == 0:
        return np.inf  # Avoid division by zero if start and end points are the same
    return path_length / euclidean_distance

# Example path (list of 3D coordinates)
path = np.array([
    [0, 0, 0],
    [1, 1, 1],
    [2, 0, 2],
    [3, 1, 3],
    [4, 0, 4]
])

tortuosity_value = calculate_tortuosity(path)
print(f"Tortuosity: {tortuosity_value}")
```

--------------------------------

### Simulate and Save Impedance (Small-Large)

Source: https://taufactor.readthedocs.io/en/latest/notebooks/10-tau-x-validation.html

Initializes and solves a PeriodicImpedanceSolver for the 'small-large' configuration with an increased iteration limit, then saves the resulting object.

```python
d = tau.PeriodicImpedanceSolver(fields["small-large"], mode='nyquist')
d.solve(iter_limit=1000000, conv_crit=-1, verbose='plot', plot_interval=100)
pkl_path = save_solver_object(d, "impedance_solve_small-large")
```

--------------------------------

### Initialize and Solve with Multi-Phase Conductivity

Source: https://taufactor.readthedocs.io/en/latest/notebooks/00-overview.html

Initializes and solves using the MultiPhaseSolver, allowing for analysis of images with more than two conductive phases. Conductivity values for each phase are provided.

```python
# assign conductivity values, where key is segmented label in 'img'
# and value is conductivity
cond = {1:0.32, 2:0.44}

# create a multiphase solver object and set an iteration limit
s = tau.MultiPhaseSolver(img, cond=cond, iter_limit=1000)
# call solve function
s.solve()
```

--------------------------------

### Creating a Tau Factor Solver Instance

Source: https://taufactor.readthedocs.io/en/latest/_sources/notebooks/00-overview.ipynb.txt

Initializes a Tau Factor solver. The `img==0` argument suggests a specific configuration or initial state for the solver, likely related to the input image data.

```python
s = tau.Solver(img==0)
```

--------------------------------

### Calculate Pairwise Interfacial Areas (Unnormalized)

Source: https://taufactor.readthedocs.io/en/latest/notebooks/02-surface-areas.html

Computes the total pairwise interfacial areas for a multi-label structure without normalization. Use this to get the raw count of shared faces between different labels.

```python
tau.interfacial_areas(cubes, normalize = False)
```

--------------------------------

### Run Benchmark Case

Source: https://taufactor.readthedocs.io/en/latest/_modules/taufactor/benchmark.html

Executes a single benchmark case, including structure resolution, solver instantiation, and the solve operation. It handles device-specific setup like CUDA cache clearing and synchronization, and resets PyTorch Dynamo.

```python
def run_benchmark_case(
    N: int,
    device: str,
    conv_crit: float,
    structure: str | Callable = "fcc",
    features: int | None = None,
    iter_limit: int = 10000,
    solver: str | type | None = None,
    solver_kwargs: dict | None = None,
    solve_kwargs: dict | None = None,
) -> dict:
    """Run a single structure benchmark case.

    Args:
        structure: Either one of the predefined keys in ``STRUCTURE_REGISTRY``
            or a custom callable hook returning a structure array.
    """
    cube, structure_name = resolve_structure(structure, N=N, features=features)

    solver_cls = resolve_solver(solver)
    solver_kwargs = dict(solver_kwargs or {})
    solve_kwargs = dict(solve_kwargs or {})

    if device == "cuda":
        torch.cuda.empty_cache()
    torch._dynamo.reset()
    gc.collect()

    if device == "cuda":
        torch.cuda.synchronize()
    start_time = time.perf_counter()

    buf = io.StringIO()
    with redirect_stdout(buf):
        solver = solver_cls(cube, device=device, **solver_kwargs)
        solver.solve(iter_limit=iter_limit, conv_crit=conv_crit, **solve_kwargs)
```

--------------------------------

### __init__

Source: https://taufactor.readthedocs.io/en/latest/_modules/taufactor/taufactor.html

Initializes the MultiPhaseSolver. It validates diffusivities, sets up internal structures for volume fractions and mean diffusivities, and handles missing diffusivity definitions.

```APIDOC
## __init__(img, diffusivities=None, D_scaling=1, omega=None, device='cuda')

### Description
Initializes the MultiPhaseSolver. It validates diffusivities, sets up internal structures for volume fractions and mean diffusivities, and handles missing diffusivity definitions.

### Parameters
- **img** (numpy.ndarray): Labeled image.
- **diffusivities** (dict[int, float], optional): Map `label -> diffusivity`. Defaults to `{1: 1}`.
- **D_scaling** (float, optional): Scaling factor for diffusivity. Defaults to 1.
- **omega** (float, optional): Relaxation parameter for SOR. Defaults to None.
- **device** (str | torch.device, optional): Compute device. Defaults to `'cuda'`.

### Raises
- `ValueError`: If any diffusivity is negative or non-finite.
- `TypeError`: If phase labels are not integers or diffusivities are not floats.

### Attributes Initialized
- **Ds** (dict[int, float]): Internal map of label to diffusivity.
- **conductive_labels** (list[int]): List of labels with positive diffusivity.
- **VF** (dict[int, numpy.ndarray]): Volume fraction per label and batch.
- **D_0** (float): Diffusivity scaling factor.
- **D_mean** (numpy.ndarray): Phase-weighted mean diffusivity per batch.
```

--------------------------------

### Simulate and Save Impedance (Diagonals)

Source: https://taufactor.readthedocs.io/en/latest/notebooks/10-tau-x-validation.html

Initializes and solves a PeriodicImpedanceSolver for the 'diagonals' configuration, then saves the resulting object.

```python
a = tau.PeriodicImpedanceSolver(fields['diagonals'], mode='nyquist')
a.solve(iter_limit=500000, conv_crit=-1, verbose='plot', plot_interval=100)
pkl_path = save_solver_object(a, "impedance_solve_diagonals")
```

--------------------------------

### Simulate and Save Impedance (Diag-Pillars)

Source: https://taufactor.readthedocs.io/en/latest/notebooks/10-tau-x-validation.html

Initializes and solves a PeriodicImpedanceSolver for the 'diag-pillars' configuration, then saves the resulting object.

```python
b = tau.PeriodicImpedanceSolver(fields["diag-pillars"], mode='nyquist')
b.solve(iter_limit=500000, conv_crit=-1, verbose='plot', plot_interval=100)
pkl_path = save_solver_object(b, "impedance_solve_diag-pillars")
```

--------------------------------

### Initialize Impedance Solver

Source: https://taufactor.readthedocs.io/en/latest/_modules/taufactor/electrode.html

Initializes the ImpedanceSolver with image data, conductive/reactive labels, frequency parameters, spacing, mode, and device. It calculates initial conductive and reactive neighbors, volume fraction, and effective area. Raises a TypeError if batch_size > 1.

```python
def __init__(self, img, conductive_label=1, reactive_label=0, 
                 omega=None, spacing=None, mode='tau_e', device='cuda'):
        self.left_bc = 1.0
        self.electrode_bc = 0.0
        self.cond_label=conductive_label
        self.reac_label=reactive_label
        self.dx = spacing or 1
        super().__init__(img, omega=omega, device=device, precision=torch.double)
        if self.batch_size > 1:
            raise TypeError(
                "Error: The ImpedanceSolver is only implemented for batch_size=1! "
                "TODO: vectorize impedance solver for batching."
            )
        
        torch_img = torch.tensor(self.cpu_img, dtype=self.precision, device=self.device)
        self.mask = torch.zeros_like(torch_img)
        self.mask[torch_img==self.cond_label] = 1
        self.cond_nn, self.reac_nn = self.count_neighbours(torch_img, self.mask)

        # Volume fraction (slice-wise)
        self.vol_x = torch.mean(self.mask, (0, 2, 3)).cpu().numpy()
        self.a_x = (torch.sum(self.reac_nn, (0, 2, 3)) / (self.Ny*self.Nz*self.dx)).cpu().numpy()

        # Define frequency, resistance and capacitance
        self.resistance = 1
        self.capacitance = 1
        self.freq_0 = np.mean(self.vol_x) / np.mean(self.a_x*self.dx) / self.Nx**2
        if mode == 'tau_e':
            self.frequency = np.array([self.freq_0/2, self.freq_0/4])
        elif mode == 'nyquist':
            freq_range = self.freq_0 * 2 ** (np.arange(-3, 10, 0.5))
            self.frequency = freq_range[::-1]
        self.Z_mean = compute_impedance(np.full(self.Nx, 1/np.mean(self.vol_x)), 
                                        np.full(self.Nx, np.mean(self.a_x*self.dx)), 
                                        self.frequency)

        # solving params
        self.sum_iter = 0
        self.c_x = 0
```

--------------------------------

### Simulate and Save Impedance (Large-Small)

Source: https://taufactor.readthedocs.io/en/latest/notebooks/10-tau-x-validation.html

Initializes and solves a PeriodicImpedanceSolver for the 'large-small' configuration with an increased iteration limit, then saves the resulting object.

```python
c = tau.PeriodicImpedanceSolver(fields["large-small"], mode='nyquist')
c.solve(iter_limit=1000000, conv_crit=-1, verbose='plot', plot_interval=100)
pkl_path = save_solver_object(c, "impedance_solve_large-small")
```

--------------------------------

### _init_device

Source: https://taufactor.readthedocs.io/en/latest/_modules/taufactor/taufactor.html

Initializes and returns a PyTorch device object. It checks for CUDA availability and defaults to CPU if CUDA is not available.

```APIDOC
## _init_device

### Description
Initializes and returns a PyTorch device object, defaulting to CPU if CUDA is not available.

### Method
`_init_device(device) -> torch.device`

### Parameters
#### Path Parameters
- **device** (str or torch.device) - The desired device (e.g., 'cuda', 'cpu').

### Returns
- **device** (torch.device) - The initialized PyTorch device.
```

--------------------------------

### Loading and Displaying an Image

Source: https://taufactor.readthedocs.io/en/latest/_sources/notebooks/00-overview.ipynb.txt

This snippet shows how to load an image file using `tau.util.load_image` and then display it. This is a common first step in image processing workflows.

```python
img = tau.util.load_image("../data/phantom.png")
plt.imshow(img)
```

--------------------------------

### Basic Surface Area Calculations

Source: https://taufactor.readthedocs.io/en/latest/_sources/notebooks/02-surface-areas.ipynb.txt

Demonstrates fundamental surface area calculations for common geometric shapes. Ensure necessary libraries are imported.

```python
import numpy as np
import pandas as pd
from typing import List, Tuple

import plotly.express as px
import plotly.graph_objects as go

from plotly.subplots import make_subplots

from sklearn.neighbors import NearestNeighbors

from scipy.spatial import Delaunay
from scipy.spatial.qhull import QhullError

from shapely.geometry import Polygon, MultiPolygon, Point, LineString
from shapely.ops import unary_union

from .utils import plot_polygons, plot_points, plot_lines, plot_multipolygons, plot_point_cloud, plot_mesh, plot_surface, plot_triangulation, plot_mesh_and_surface, plot_mesh_and_triangulation, plot_mesh_and_point_cloud, plot_mesh_and_lines, plot_mesh_and_polygons, plot_mesh_and_multipolygons, plot_mesh_and_triangulation_and_point_cloud, plot_mesh_and_triangulation_and_lines, plot_mesh_and_triangulation_and_polygons, plot_mesh_and_triangulation_and_multipolygons, plot_mesh_and_triangulation_and_point_cloud_and_lines, plot_mesh_and_triangulation_and_point_cloud_and_polygons, plot_mesh_and_triangulation_and_point_cloud_and_multipolygons, plot_mesh_and_triangulation_and_lines_and_polygons, plot_mesh_and_triangulation_and_lines_and_multipolygons, plot_mesh_and_triangulation_and_polygons_and_multipolygons, plot_mesh_and_triangulation_and_point_cloud_and_lines_and_polygons, plot_mesh_and_triangulation_and_point_cloud_and_lines_and_multipolygons, plot_mesh_and_triangulation_and_point_cloud_and_polygons_and_multipolygons, plot_mesh_and_triangulation_and_lines_and_polygons_and_multipolygons, plot_mesh_and_triangulation_and_point_cloud_and_lines_and_polygons_and_multipolygons


def triangle_area(p1: np.ndarray, p2: np.ndarray, p3: np.ndarray) -> float:
    """Calculate the area of a triangle given three points."""
    return 0.5 * np.abs(p1[0] * (p2[1] - p3[1]) + p2[0] * (p3[1] - p1[1]) + p3[0] * (p1[1] - p2[1]))

def polygon_area(vertices: np.ndarray) -> float:
    """Calculate the area of a polygon given its vertices."""
    n = len(vertices)
    if n < 3:
        return 0.0
    area = 0.0
    for i in range(n):
        j = (i + 1) % n
        area += vertices[i][0] * vertices[j][1]
        area -= vertices[j][0] * vertices[i][1]
    return np.abs(area) / 2.0

def surface_area_of_mesh(vertices: np.ndarray, faces: np.ndarray) -> float:
    """Calculate the total surface area of a mesh."""
    total_area = 0.0
    for face in faces:
        p1 = vertices[face[0]]
        p2 = vertices[face[1]]
        p3 = vertices[face[2]]
        total_area += triangle_area(p1, p2, p3)
    return total_area

def surface_area_of_point_cloud(points: np.ndarray) -> float:
    """Calculate the surface area of a point cloud using Delaunay triangulation."""
    try:
        tri = Delaunay(points)
        return surface_area_of_mesh(points, tri.simplices)
    except QhullError:
        return 0.0

def surface_area_of_polygon(polygon: Polygon) -> float:
    """Calculate the surface area of a Shapely Polygon."""
    return polygon.area

def surface_area_of_multipolygon(multipolygon: MultiPolygon) -> float:
    """Calculate the surface area of a Shapely MultiPolygon."""
    return multipolygon.area

def surface_area_of_line(line: LineString) -> float:
    """Calculate the surface area of a Shapely LineString (which is 0)."""
    return line.area

def surface_area_of_point(point: Point) -> float:
    """Calculate the surface area of a Shapely Point (which is 0)."""
    return point.area

def calculate_surface_areas(data: dict) -> dict:
    """Calculate surface areas for various geometric types in a dictionary."""
    results = {}
    for key, value in data.items():
        if isinstance(value, np.ndarray):
            if value.ndim == 2 and value.shape[1] == 3:  # Assuming points or vertices
                try:
                    results[f'{key}_area'] = surface_area_of_point_cloud(value)
                except QhullError:
                    results[f'{key}_area'] = 0.0
            elif value.ndim == 2 and value.shape[0] == 3: # Assuming faces for a mesh
                # This case is ambiguous without vertices, assuming it's not directly calculable here
                pass
        elif isinstance(value, Polygon):
            results[f'{key}_area'] = surface_area_of_polygon(value)
        elif isinstance(value, MultiPolygon):
            results[f'{key}_area'] = surface_area_of_multipolygon(value)
        elif isinstance(value, LineString):
            results[f'{key}_area'] = surface_area_of_line(value)
        elif isinstance(value, Point):
            results[f'{key}_area'] = surface_area_of_point(value)
        elif isinstance(value, dict) and 'vertices' in value and 'faces' in value:
            results[f'{key}_area'] = surface_area_of_mesh(value['vertices'], value['faces'])
    return results


# Example Usage:

# 1. Point Cloud (approximated surface area)
points = np.array([[0, 0, 0], [1, 0, 0], [0, 1, 0], [0, 0, 1], [1, 1, 0], [1, 0, 1], [0, 1, 1], [1, 1, 1]])

# 2. Mesh (vertices and faces)
vertices = np.array([[0, 0, 0], [1, 0, 0], [0, 1, 0], [0, 0, 1]])
faces = np.array([[0, 1, 2], [0, 1, 3], [0, 2, 3], [1, 2, 3]]) # Tetrahedron faces
mesh_data = {'vertices': vertices, 'faces': faces}

# 3. Shapely Polygon
polygon = Polygon([(0, 0), (1, 0), (1, 1), (0, 1)])

# 4. Shapely MultiPolygon
poly1 = Polygon([(0, 0), (1, 0), (1, 1), (0, 1)])
poly2 = Polygon([(2, 2), (3, 2), (3, 3), (2, 3)])
multipolygon = MultiPolygon([poly1, poly2])

# 5. Shapely LineString
line = LineString([(0, 0), (1, 1)])

# 6. Shapely Point
point = Point(0, 0)

# Combine all data
all_data = {
    'point_cloud': points,
    'mesh': mesh_data,
    'simple_polygon': polygon,
    'complex_multipolygon': multipolygon,
    'a_line': line,
    'a_point': point
}

# Calculate all surface areas
surface_areas = calculate_surface_areas(all_data)
print(surface_areas)

# Plotting examples (assuming plotting functions are available)
# plot_point_cloud(points)
# plot_mesh(vertices, faces)
# plot_polygons([polygon])
# plot_multipolygons([multipolygon])
# plot_lines([line])
# plot_points([point])

```

--------------------------------

### Initialize and Solve PeriodicElectrodeSolver

Source: https://taufactor.readthedocs.io/en/latest/_sources/notebooks/01-tortuosity.ipynb.txt

Initializes the solver with specified batch, device, and spacing. Solves for tortuosity with a given iteration limit, convergence criterion, and plotting interval.

```python
s = tau.PeriodicElectrodeSolver(batch, device='cuda', spacing=Lx/nx)

```

```python
tau_pore = s.solve(iter_limit=200000, conv_crit=1e-5, verbose='plot', plot_interval=100)
```

--------------------------------

### Create Multi-Feature Test Structures with Labeling

Source: https://taufactor.readthedocs.io/en/latest/notebooks/12-solver-benchmark.html

Creates test structures with multiple features and applies periodic labeling. Useful for benchmarking solvers that handle complex geometries or periodic boundary conditions.

```python
from taufactor.metrics import label_periodic
from scipy.ndimage import generate_binary_structure

fields = {}
features = 2
fields["blocks"] = create_stacked_blocks(Nx, features=features)
fields["2d_diagonals"] = create_2d_diagonals(Nx, features=features)
fields["2d_zigzag"] = create_2d_zigzag(Nx, features=features)
fields["3d_diagonals"] = create_3d_diagonals(Nx, features=features)

fields["blocks"][Nx//2:,:Nx//2,Nx//2:] *= 2
connectivity = generate_binary_structure(3, 1)
fields["2d_diagonals"] = label_periodic(fields["2d_diagonals"], 1, connectivity, periodic=(False, True, True))[0]
fields["2d_zigzag"] = label_periodic(fields["2d_zigzag"], 1, connectivity, periodic=(False, True, True))[0]
fields["3d_diagonals"] = label_periodic(fields["3d_diagonals"], 1, connectivity, periodic=(False, True, True))[0]

colors = ['white', 'gray', 'red', 'blue', 'green']
fig, axes = plt.subplots(nrows=1, ncols=4, figsize=(15, 4), dpi=100)
for i, structure in enumerate(fields.keys()):
    axes[i].remove()
    axes[i] = fig.add_subplot(1, 5, i + 1, projection='3d')
    for j in np.unique(fields[structure]):
        if j>0:
            axes[i].voxels(np.transpose(fields[structure], (2,1,0))==j, facecolors=colors[j], edgecolor='black', alpha=1.0)
    axes[i].set_zlim(0, Nx)
    axes[i].set_aspect('equal')
    axes[i].set_title(structure)
    axes[i].set_xlim(0, Nx)
    axes[i].set_ylim(0, Nx)
plt.tight_layout()
plt.show()
```

--------------------------------

### Create and Checkout New Branch

Source: https://taufactor.readthedocs.io/en/latest/contributing.html

Navigate into the cloned directory and create a new branch for your development work. This isolates your changes.

```bash
cd taufactor
git checkout -b my_dev_branch
```

--------------------------------

### run_benchmark_study

Source: https://taufactor.readthedocs.io/en/latest/_autosummary/taufactor.benchmark.html

Conducts a comprehensive convergence benchmark across various structure sizes, devices, and convergence criteria, saving results to a file.

```APIDOC
## run_benchmark_study(Ns=(100, 128, 200, 256, 300, 384, 400), devices=('cuda',), conv_crit_values=(0.001,), structure='fcc', features=1, outfile='taufactor_benchmark_results.txt', write_file=True, iter_limit=10000, solver=None, solver_kwargs=None, solve_kwargs=None)

### Description
Run a convergence benchmark on synthetic structures.

### Parameters
- **Ns** (list[int] | tuple[int, ...], optional) - A list or tuple of structure sizes to benchmark. Defaults to (100, 128, 200, 256, 300, 384, 400).
- **devices** (list[str] | tuple[str, ...], optional) - A list or tuple of devices to run the benchmarks on. Defaults to ('cuda',).
- **conv_crit_values** (list[float] | tuple[float, ...], optional) - A list or tuple of convergence criteria to test. Defaults to (0.001,).
- **structure** (str | Callable, optional) - The structure to benchmark. Defaults to 'fcc'. Can be a predefined name ('fcc', 'blocks', 'diagonal2d', 'zigzag', 'diagonal3d') or a custom callable hook.
- **features** (int, optional) - Number of features to apply to the structure. Defaults to 1.
- **outfile** (str, optional) - The name of the output file for benchmark results. Defaults to 'taufactor_benchmark_results.txt'.
- **write_file** (bool, optional) - Whether to write the results to a file. Defaults to True.
- **iter_limit** (int, optional) - The maximum number of iterations for each benchmark case. Defaults to 10000.
- **solver** (str | type | None, optional) - The solver to use. Defaults to None.
- **solver_kwargs** (dict, optional) - Keyword arguments for the solver. Defaults to None.
- **solve_kwargs** (dict, optional) - Keyword arguments for the solve method. Defaults to None.

### Returns
- list[dict] - A list of dictionaries, where each dictionary contains the results of a benchmark case.
```
Tokens:
Raw