pcalc.presence_matrix

  1# -*- coding: utf-8 -*-
  2# Copyright (c) 2025 Krishna Kumar
  3# SPDX-License-Identifier: MIT
  4
  5from __future__ import annotations
  6
  7from typing import List, Optional
  8
  9import numpy as np
 10from numpy import typing as npt
 11
 12from .presence import PresenceAssertion
 13from .presence_map import PresenceMap
 14from .time_scale import Timescale
 15
 16
 17class PresenceMatrix:
 18    """
 19    A scale-invariant, optionally lazy matrix representation of element presences over discrete time bins.
 20
 21    Each row corresponds to a single `Presence`, and each column represents a discrete time bin defined by a `Timescale`.
 22    The matrix records the degree of presence of each element within a given boundary over time, where each value is a real number in [0.0, 1.0]:
 23    - 1.0 = full presence for the duration of the time bin
 24    - 0.0 = no presence
 25    - values in between represent partial presence due to overlap
 26
 27    The matrix can be computed eagerly via `materialize()` or accessed lazily through NumPy-style slicing.
 28
 29    ### Structure
 30
 31    - **Rows:** Each row corresponds to one `Presence` instance.
 32    - **Columns:** Each column represents a discrete bin of time defined by the `Timescale` (from `t0` to `t1` using `bin_width`).
 33    - The matrix shape is therefore `(len(presences), timescale.num_bins)`.
 34
 35    ### Usage
 36
 37    ```python
 38    from pcalc import PresenceMatrix
 39    from metamodel.timescale import Timescale
 40
 41    # Define a time window with bin width of 1.0
 42    ts = Timescale(t0=0.0, t1=5.0, bin_width=1.0)
 43
 44    # Construct lazily
 45    matrix = PresenceMatrix(presences, time_scale=ts)
 46
 47    # Access entire row
 48    row = matrix[0]  # row 0 as 1D array
 49
 50    # Access single value
 51    value = matrix[2, 3]  # presence at row 2, column 3
 52
 53    # Access slice of row
 54    window = matrix[1, 1:4]  # row 1, columns 1–3
 55
 56    # Access slice of rows
 57    rows = matrix[0:2]  # rows 0 and 1
 58
 59    # Access single column
 60    col = matrix[:, 2]  # all rows at column 2
 61
 62    # Access column block
 63    col_block = matrix[:, 1:4]  # all rows, columns 1–3
 64
 65    # Materialize full matrix (optional)
 66    dense = matrix.materialize()
 67    ```
 68
 69    Features
 70    --------
 71    - Lazy evaluation by default; full matrix computed only if needed
 72    - NumPy-style indexing: supports full rows, row/column slicing, and scalar access
 73    - Precision-preserving: partial bin overlaps are retained in the matrix
 74    - Compatible with matrix operations and flow metric calculations
 75    - Backed by a sparse `PresenceMap`, which stores slice indices and fractional overlaps
 76
 77    Notes
 78    -----
 79    - Use `materialize()` if you want to extract or operate on the full dense matrix
 80    - Use `matrix[i, j]`, `matrix[i, j:k]`, or `matrix[:, j]` for lightweight slicing without allocating the full matrix
 81    - Stepped slicing (e.g., `matrix[::2]`) is not currently supported
 82
 83"""
 84
 85    def __init__(self, presences: List[PresenceAssertion], time_scale: Timescale, materialize=False):
 86        """
 87        Construct a presence matrix from a list of Presences and time window configuration.
 88
 89        Args:
 90            presences: A list of `Presence` instances, one per element of interest.
 91            time_scale: The discrete `Timescale` that the matrix will be normalized to.
 92            materialize: If True, immediately constructs the full backing matrix.
 93                     Otherwise, matrix values will be computed on demand using the sparse presence map.
 94        """
 95
 96        self.presence_matrix: Optional[npt.NDArray[np.float64]] = None
 97        """
 98        A real valued matrix with shape (num_Presences, num_bins).
 99        When every presence.start and presence.end are whole numbers every value
100        in the matrix is either 0.0 or 1.0. If they are not whole numbers, the presence
101        is mapped to a number between 0.0 and 1.0 at the start and end (or both), with any
102        intermediate value being 1.0. 
103        """
104
105        self.time_scale = time_scale
106        """
107        The time scale of the presence matrix.
108        """
109
110        self.presence_map: List[PresenceMap] = []
111        self.shape = None
112        self.init_presence_map(presences)
113
114        if materialize:
115            self.materialize()
116
117    def init_presence_map(self, presences: List[PresenceAssertion]) -> None:
118        """
119        Initialize the internal presence matrix based on the Presence intervals and binning scheme.
120        Only presences that overlap the timescale endpoints [t0, t1) are mapped.
121        """
122        ts = self.time_scale  # Timescale object: includes t0, t1, bin_width
123
124        for row, presence in enumerate(presences):
125            presence_map = PresenceMap(presence, ts)
126            if presence_map.is_mapped:
127                self.presence_map.append(presence_map)
128
129        num_bins = ts.num_bins
130        num_rows = len(self.presence_map)
131        self.shape = (num_rows, num_bins)
132
133    def is_materialized(self) -> bool:
134        """Return True if the backing matrix has been materialized."""
135        return self.presence_matrix is not None
136
137    def materialize(self) -> npt.NDArray[np.float64]:
138        """
139        Materialize and return the full presence matrix from presence maps.
140        If already materialized, this is a no-op and returns the cached matrix.
141
142        Use this only if want the full matrix to operate on. For most cases, you should
143        be able to work with the sparse representation with the presence map and using
144        array slicing on the matrix object. So think twice about why you are materializing
145        a matrix.
146
147        Returns:
148            The dense presence matrix of shape (num_presences, num_bins).
149        """
150        if self.is_materialized():
151            return self.presence_matrix
152
153        num_rows, num_cols = self.shape
154        matrix = np.zeros((num_rows, num_cols), dtype=float)
155
156        for row, pm in enumerate(self.presence_map):
157            if not pm.is_mapped:
158                continue
159            start = pm.start_bin
160            end = pm.end_bin
161            matrix[row, start] = pm.start_value
162            if end - 1 > start:
163                matrix[row, end - 1] = pm.end_value
164            if end - start > 2:
165                matrix[row, start + 1: end - 1] = 1.0
166
167        self.presence_matrix = matrix
168        return self.presence_matrix
169
170    def drop_materialization(self) -> None:
171        """Discard the cached matrix, returning to sparse-only mode."""
172        self.presence_matrix = None
173
174    def _compute_row_slice(self, row: int, start: int, stop: int) -> np.ndarray:
175        """
176        On-demand computation of a row slice from presence map.
177        If the presence is not mapped or out of bounds, returns zeroes.
178        """
179        pm = self.presence_map[row]
180        width = stop - start
181        output = np.zeros(width, dtype=float)
182
183        if not pm.is_mapped:
184            return output
185
186        # Slice and overlap window
187        for col in range(start, stop):
188            if col < pm.start_bin or col >= pm.end_bin:
189                continue
190            elif col == pm.start_bin:
191                output[col - start] = pm.start_value
192            elif col == pm.end_bin - 1:
193                output[col - start] = pm.end_value
194            else:
195                output[col - start] = 1.0
196
197        return output
198
199    @property
200    def presences(self) -> List[PresenceAssertion]:
201        return [pm.presence for pm in self.presence_map]
202
203    def __getitem__(self, index):
204        """
205        Support NumPy-style indexing:
206        - matrix[i]           → row i
207        - matrix[i:j]         → row slice (returns 2D array)
208        - matrix[i, j]        → scalar
209        - matrix[i, j:k]      → row[i], column slice
210        - matrix[:, j]        → column j (all rows)
211        - matrix[:, j:k]      → column block j:k (all rows)
212        """
213        # matrix[i:j] → top-level row slicing
214        if isinstance(index, slice):
215            if index.step is not None:
216                raise ValueError("PresenceMatrix does not support slicing with a step.")
217            start = index.start if index.start is not None else 0
218            stop = index.stop if index.stop is not None else self.shape[0]
219            if self.presence_matrix is not None:
220                return self.presence_matrix[start:stop]
221            return np.array([self[row] for row in range(start, stop)])
222
223        # matrix[i, j], matrix[i, j:k], matrix[:, j]
224        if isinstance(index, tuple):
225            row_idx, col_idx = index
226
227            # matrix[:, j] or matrix[:, j:k]
228            if isinstance(row_idx, slice):
229                if row_idx.step is not None:
230                    raise ValueError("PresenceMatrix does not support stepped row slicing.")
231                row_start = row_idx.start if row_idx.start is not None else 0
232                row_stop = row_idx.stop if row_idx.stop is not None else self.shape[0]
233
234                if isinstance(col_idx, int):
235                    return np.array([self[row, col_idx] for row in range(row_start, row_stop)])
236
237                elif isinstance(col_idx, slice):
238                    if col_idx.step is not None:
239                        raise ValueError("PresenceMatrix does not support stepped column slicing.")
240                    col_start = col_idx.start if col_idx.start is not None else 0
241                    col_stop = col_idx.stop if col_idx.stop is not None else self.shape[1]
242                    return np.array([self[row, col_start:col_stop] for row in range(row_start, row_stop)])
243
244            # matrix[i, j] or matrix[i, j:k]
245            elif isinstance(row_idx, int):
246                if isinstance(col_idx, int):
247                    return self.presence_matrix[row_idx, col_idx] if self.presence_matrix is not None \
248                        else self._compute_row_slice(row_idx, col_idx, col_idx + 1)[0]
249
250                elif isinstance(col_idx, slice):
251                    if col_idx.step is not None:
252                        raise ValueError("PresenceMatrix does not support stepped column slicing.")
253                    start = col_idx.start if col_idx.start is not None else 0
254                    stop = col_idx.stop if col_idx.stop is not None else self.shape[1]
255                    return self.presence_matrix[row_idx, start:stop] if self.presence_matrix is not None \
256                        else self._compute_row_slice(row_idx, start, stop)
257
258        # matrix[i] → single row
259        if isinstance(index, int):
260            return self.presence_matrix[index] if self.presence_matrix is not None \
261                else self._compute_row_slice(index, 0, self.shape[1])
262
263        raise TypeError(f"Invalid index for PresenceMatrix: {index}")