Source code for mercurial.hierarchy.complexity

"""Complexity measure Λ(x) for LADDER level classification (Definition 3.2)."""

from typing import Tuple

import numpy as np
from scipy.linalg import svd
from scipy.spatial.distance import pdist




class ComplexityMeasure:
    """
    Computes Λ(x) = α log₂ N_dof + β I_int + γ Σ

    where:
    - N_dof: effective degrees of freedom (rank of covariance matrix)
    - I_int: integrated information (Φ, simplified from IIT)
    - Σ: organizational entropy (based on mutual information between components)
    """

    def __init__(self, alpha: float = 1.0, beta: float = 1.0, gamma: float = 0.1):
        self.alpha = alpha
        self.beta = beta
        self.gamma = gamma



    def effective_dof(self, pattern_state: np.ndarray) -> float:
        """
        Compute N_dof ≈ rank of covariance matrix (numerical rank).
        Uses singular value decomposition and count singular values above noise floor.
        """
        # Ensure 2D array: (samples, dimensions)
        if pattern_state.ndim == 1:
            pattern_state = pattern_state.reshape(-1, 1)
        if pattern_state.shape[0] == 1:
            return 1.0

        # Center the data
        centered = pattern_state - np.mean(pattern_state, axis=0)
        # Compute SVD
        U, S, Vt = svd(centered, full_matrices=False)
        # Noise floor: assume smallest singular value is noise
        noise_floor = np.median(S) * 0.01 if len(S) > 0 else 0.0
        # Count singular values above noise floor
        rank = np.sum(S > noise_floor)
        return float(max(rank, 1))




    def integrated_information(
        self, pattern_state: np.ndarray, partition_minsize: int = 2
    ) -> float:
        """
        Approximate Φ (I_int) using a simplified IIT measure.
        Φ = min_{partition} (I(whole) - Σ I(parts))
        We approximate by splitting the dimensions into two random partitions.
        """
        if pattern_state.ndim == 1:
            pattern_state = pattern_state.reshape(-1, 1)
        n_samples, n_dims = pattern_state.shape
        if n_dims < partition_minsize * 2:
            # Too few dimensions, return mutual information between all pairs
            return self._mutual_information_all(pattern_state)

        # Compute total mutual information (whole system)
        I_whole = self._mutual_information_all(pattern_state)

        # Try several random bipartitions and take minimum
        n_partitions = min(10, n_dims)
        min_partition_info = I_whole
        for _ in range(n_partitions):
            perm = np.random.permutation(n_dims)
            part1 = perm[: n_dims // 2]
            part2 = perm[n_dims // 2 :]
            # Information of parts: sum of mutual information within each part
            I_parts = self._mutual_information_all(
                pattern_state[:, part1]
            ) + self._mutual_information_all(pattern_state[:, part2])
            min_partition_info = min(min_partition_info, I_parts)

        return I_whole - min_partition_info


    def _mutual_information_all(self, data: np.ndarray) -> float:
        """
        Estimate total mutual information (sum of pairwise mutual information)
        as a proxy for integrated information.
        """
        if data.shape[1] <= 1:
            return 0.0
        # Compute pairwise mutual information using Gaussian approximation
        # For simplicity, we use correlation as proxy
        corr = np.corrcoef(data.T)
        # Mutual information for Gaussian: -0.5 * log(1 - ρ²)
        # Avoid negative or zero correlations
        rho_sq = np.clip(corr**2, 0.0, 0.9999)
        mi_matrix = -0.5 * np.log(1 - rho_sq)
        # Sum upper triangle (excluding diagonal)
        n = mi_matrix.shape[0]
        total_mi = np.sum(mi_matrix[np.triu_indices(n, k=1)])
        return total_mi



    def organizational_entropy(self, pattern_state: np.ndarray) -> float:
        """
        Compute Σ (organizational entropy) based on the distribution of pairwise distances.
        Lower Σ means more ordered organization.
        """
        if pattern_state.ndim == 1:
            pattern_state = pattern_state.reshape(-1, 1)
        if pattern_state.shape[0] < 2:
            return 0.0

        # Pairwise Euclidean distances between samples
        dists = pdist(pattern_state)
        # Histogram the distances to estimate entropy
        hist, bin_edges = np.histogram(dists, bins="auto", density=True)
        hist = hist[hist > 0]
        # Shannon entropy of distance distribution
        entropy = -np.sum(hist * np.log(hist + 1e-12))
        # Normalize by log(number of bins) to get 0-1 range
        max_entropy = np.log(len(hist)) if len(hist) > 0 else 1.0
        return entropy / max_entropy if max_entropy > 0 else 0.0




    def compute(self, pattern_state: np.ndarray) -> float:
        """
        Compute Λ = α log₂ N_dof + β I_int + γ Σ
        """
        N = self.effective_dof(pattern_state)
        I_int = self.integrated_information(pattern_state)
        Sigma = self.organizational_entropy(pattern_state)

        term1 = self.alpha * np.log2(N) if N > 0 else 0.0
        term2 = self.beta * I_int
        term3 = self.gamma * Sigma

        return term1 + term2 + term3






class LevelClassifier:
    """
    Assigns a LADDER level to a pattern based on its complexity measure.
    Uses predefined thresholds (calibratable).
    """

    # Predefined complexity thresholds for levels 0-18 (from LADDER)
    # These are initial estimates; can be calibrated with data
    LEVEL_THRESHOLDS = [
        0.0,  # Level 0
        0.5,
        1.0,
        1.5,
        2.0,
        2.5,
        3.0,
        3.5,
        4.0,
        4.5,
        5.0,  # Levels 1-10
        6.0,
        7.0,
        8.0,
        9.0,
        10.0,
        11.0,
        12.0,
        13.0,
        14.0,  # Levels 11-18
    ]

    LEVEL_NAMES = [
        "Neutral Substrate",
        "Quantum-Gravitational",
        "Elementary Particles",
        "Composite Particles",
        "Atomic Species",
        "Molecular Aggregates",
        "Mesoscale Structures",
        "Macroscopic Materials",
        "Functional Assemblies",
        "Integrated Biological",
        "Networked Systems",
        "Planetary Systems",
        "Interplanetary",
        "Interstellar",
        "Galactic",
        "Intergalactic",
        "Cosmic Web",
        "Observable Universe",
        "Entire Branch",
    ]



    @classmethod
    def classify(cls, complexity: float) -> Tuple[int, str]:
        """Return (level_index, level_name) for given complexity."""
        for i, thresh in enumerate(cls.LEVEL_THRESHOLDS):
            if complexity < thresh:
                return max(0, i - 1), cls.LEVEL_NAMES[max(0, i - 1)]
        return len(cls.LEVEL_THRESHOLDS) - 1, cls.LEVEL_NAMES[-1]




    @classmethod
    def get_level_range(
        cls, pattern_state: np.ndarray, measure: ComplexityMeasure
    ) -> Tuple[int, int]:
        """
        Compute complexity and return level index.
        Also returns a confidence interval (simplified).
        """
        comp = measure.compute(pattern_state)
        level, name = cls.classify(comp)
        # Rough confidence: ±1 level for uncertainty
        return level, level, comp




    def neural_complexity(self, E: np.ndarray, I: np.ndarray) -> float:
        """
        Compute Λ from neural activity (spectral and temporal features).
        """
        # Use spectral entropy as proxy
        fft_E = np.fft.fft(E - np.mean(E))
        power_E = np.abs(fft_E[: len(fft_E) // 2]) ** 2
        power_E = power_E / (np.sum(power_E) + 1e-12)
        entropy_E = -np.sum(power_E * np.log2(power_E + 1e-12))

        # Also incorporate variance
        var_E = np.var(E)
        # Normalise
        return entropy_E + var_E