class VSA:
    """High-level factory interface for creating VSA models.

    This class provides a simple, unified API for creating and using
    different VSA models. It's the recommended entry point for most users.

    Examples:
        >>> # Create a MAP model with default settings
        >>> model = VSA.create('MAP')
        >>>
        >>> # Create FHRR with specific dimension and backend
        >>> model = VSA.create('FHRR', dim=512, backend='torch', device='cuda')
        >>>
        >>> # Use the model
        >>> a, b = model.random(), model.random()
        >>> c = model.bind(a, b)
        >>> similarity = model.similarity(a, model.unbind(c, b))
    """

    # Model registry
    _MODELS = {
        "map": MAPModel,
        "fhrr": FHRRModel,
        "hrr": HRRModel,
        "bsc": BSCModel,
        "bsdc": BSDCModel,
        "bsdc_seg": BSDCSEGModel,
        "bsdc-seg": BSDCSEGModel,  # alias with hyphen to match model_name
        "ghrr": GHRRModel,
        "vtb": VTBModel,
    }

    # Default vector spaces for each model
    _DEFAULT_SPACES = {
        "map": "bipolar",
        "fhrr": "complex",
        "hrr": "real",
        "bsc": "binary",
        "bsdc": "sparse",
        "bsdc_seg": "sparse_segment",
        "bsdc-seg": "sparse_segment",  # alias with hyphen
        "ghrr": "matrix",
        "vtb": "real",
    }

    @classmethod
    def create(
        cls,
        model_type: str,
        dim: int = 10000,
        backend: str | Backend | None = None,
        space: str | None = None,
        seed: int | None = None,
        **kwargs,
    ) -> VSAModel:
        """Create a VSA model with the specified configuration.

        Args:
            model_type: Model name ('MAP', 'FHRR', 'HRR', 'BSC', etc.)
            dim: Dimensionality of hypervectors
            backend: Backend name ('numpy', 'torch', 'jax'), a :class:`Backend`
                instance, or None for the default backend
            space: Vector space name or None for model's default
            seed: Random seed for reproducibility
            **kwargs: Additional arguments passed to backend (e.g., device='cuda')

        Returns:
            VSA model instance

        Raises:
            ValueError: If model_type is not recognized

        Examples:
            >>> model = VSA.create('MAP', dim=10000)
            >>> model = VSA.create('FHRR', dim=512, backend='torch', device='cuda')
            >>> model = VSA.create('MAP', space='real')  # Use real-valued MAP
        """
        # Normalize model type
        model_type_lower = model_type.lower()

        if model_type_lower not in cls._MODELS:
            available = list(cls._MODELS.keys())
            raise ValueError(f"Unknown model type '{model_type}'. Available models: {available}")

        # Get model class
        model_class = cls._MODELS[model_type_lower]

        # Create backend
        backend_kwargs = {k: v for k, v in kwargs.items() if k in ["device"]}
        if isinstance(backend, Backend):
            backend_instance = backend
        elif isinstance(backend, str):
            backend_instance = get_backend(backend, **backend_kwargs)
        elif backend is None:
            backend_instance = None
        else:
            raise TypeError("backend must be a backend name (str) or a Backend instance")

        # Determine space type
        if space is None:
            space = cls._DEFAULT_SPACES.get(model_type_lower)

        # Collect space-specific kwargs
        space_kwargs = {}
        if space == "sparse_segment":
            # For BSDC-SEG: default to dim/10 segments (segment_length=10)
            space_kwargs["segments"] = kwargs.get("segments", max(1, dim // 10))

        # Create space if string provided
        if isinstance(space, str):
            space_instance = create_space(
                space, dimension=dim, backend=backend_instance, seed=seed, **space_kwargs
            )
        else:
            space_instance = space

        # Collect model-specific kwargs
        model_kwargs = {}
        if model_type_lower == 'bsdc':
            # BSDC supports binding_mode parameter
            if 'binding_mode' in kwargs:
                model_kwargs['binding_mode'] = kwargs['binding_mode']

        # Create model
        model = model_class(
            dimension=dim, space=space_instance, backend=backend_instance, seed=seed,
            **model_kwargs
        )

        return model

    @classmethod
    def available_models(cls) -> list[str]:
        """Return list of available model names.

        Returns:
            List of model names that can be used with create()
        """
        return list(cls._MODELS.keys())

    @classmethod
    def model_info(cls, model_type: str) -> dict:
        """Get information about a specific model.

        Args:
            model_type: Model name

        Returns:
            Dictionary with model properties

        Example:
            >>> info = VSA.model_info('FHRR')
            >>> print(info['is_exact_inverse'])  # True
        """
        model_type_lower = model_type.lower()
        if model_type_lower not in cls._MODELS:
            raise ValueError(f"Unknown model type '{model_type}'")

        # Create a temporary instance to query properties
        model = cls.create(model_type, dim=100)

        return {
            "name": model.model_name,
            "is_self_inverse": model.is_self_inverse,
            "is_commutative": model.is_commutative,
            "is_exact_inverse": model.is_exact_inverse,
            "default_space": cls._DEFAULT_SPACES.get(model_type_lower),
            "class": model.__class__.__name__,
        }

class VSAModel(ABC):
    """Abstract base class for VSA models.

    A VSA model defines the core operations:
    - bind: Associate two vectors (creates dissimilar result)
    - unbind: Recover one vector given the other and their binding
    - bundle: Combine multiple vectors (preserves similarity)
    - permute: Reorder vector to represent position/sequence

    Different models have different algebraic properties:
    - Self-inverse binding: bind(a, b) = unbind(a, b)
    - Exact vs approximate inverse
    - Commutativity of binding
    """

    def __init__(
        self,
        space: VectorSpace,
        backend: Backend | None = None
    ):
        """Initialize VSA model.

        Args:
            space: Vector space defining the representation
            backend: Computational backend (defaults to space's backend)
        """
        self.space = space
        self.backend = backend if backend is not None else space.backend
        self.dimension = space.dimension

    # ===== Core VSA Operations =====

    @abstractmethod
    def bind(self, a: Array, b: Array) -> Array:
        """Bind two vectors to create an association.

        Binding creates a new vector that is dissimilar to both inputs
        but preserves structured similarity (similar inputs → similar bindings).

        Args:
            a: First vector
            b: Second vector

        Returns:
            Bound vector representing the association of a and b
        """
        pass

    @abstractmethod
    def unbind(self, a: Array, b: Array) -> Array:
        """Unbind to recover one vector given the other.

        For self-inverse models: unbind(a, b) = bind(a, b)
        For others: approximately recovers a from bind(a, b) and b

        Args:
            a: Bound vector or first operand
            b: Second operand

        Returns:
            Recovered vector (exact or approximate depending on model)
        """
        pass

    @abstractmethod
    def bundle(self, vectors: Sequence[Array]) -> Array:
        """Bundle (superpose) multiple vectors.

        Bundling combines vectors while preserving similarity to all inputs.
        The result is similar to each input vector.

        Args:
            vectors: Sequence of vectors to bundle

        Returns:
            Bundled vector representing the superposition

        Raises:
            ValueError: If vectors is empty
        """
        pass

    @abstractmethod
    def permute(self, vec: Array, k: int = 1) -> Array:
        """Permute vector to represent position or sequence.

        Permutation reorders coordinates and is used to encode position
        or create sequences. It's invertible and preserves similarity.

        Args:
            vec: Vector to permute
            k: Number of positions to shift (default: 1)

        Returns:
            Permuted vector
        """
        pass

    def unpermute(self, vec: Array, k: int = 1) -> Array:
        """Inverse permutation.

        Args:
            vec: Vector to unpermute
            k: Number of positions to shift back (default: 1)

        Returns:
            Unpermuted vector
        """
        return self.permute(vec, -k)

    # ===== Similarity and Cleanup =====

    def similarity(self, a: Array, b: Array) -> float:
        """Compute similarity between two vectors.

        Delegates to the vector space's similarity metric.

        Args:
            a: First vector
            b: Second vector

        Returns:
            Similarity score (space-dependent metric)
        """
        return self.space.similarity(a, b)

    def normalize(self, vec: Array) -> Array:
        """Normalize vector according to space conventions.

        Args:
            vec: Vector to normalize

        Returns:
            Normalized vector
        """
        return self.space.normalize(vec)

    # ===== Vector Generation =====

    def random(self, seed: int | None = None) -> Array:
        """Generate a random vector from the space.

        Args:
            seed: Optional random seed

        Returns:
            Random vector
        """
        return self.space.random(seed=seed)

    def random_sequence(self, n: int, seed: int | None = None) -> list[Array]:
        """Generate n random vectors.

        Args:
            n: Number of vectors to generate
            seed: Optional base seed (each vector gets seed + i)

        Returns:
            List of random vectors
        """
        if seed is not None:
            return [self.random(seed=seed + i) for i in range(n)]
        return [self.random() for _ in range(n)]

    # ===== Compositional Operations =====

    def bind_multiple(self, vectors: Sequence[Array]) -> Array:
        """Bind multiple vectors sequentially.

        For n vectors: bind(bind(bind(v1, v2), v3), ...)

        Args:
            vectors: Sequence of vectors to bind

        Returns:
            Result of sequential binding

        Raises:
            ValueError: If fewer than 2 vectors provided
        """
        if len(vectors) < 2:
            raise ValueError("Need at least 2 vectors to bind")

        result = vectors[0]
        for vec in vectors[1:]:
            result = self.bind(result, vec)
        return result

    # ===== Model Properties =====

    @property
    @abstractmethod
    def is_self_inverse(self) -> bool:
        """Whether binding is self-inverse (bind = unbind)."""
        pass

    @property
    @abstractmethod
    def is_commutative(self) -> bool:
        """Whether binding is commutative (bind(a, b) = bind(b, a))."""
        pass

    @property
    @abstractmethod
    def is_exact_inverse(self) -> bool:
        """Whether unbinding gives exact recovery (no approximation error)."""
        pass

    @property
    @abstractmethod
    def model_name(self) -> str:
        """Return the model name (e.g., 'MAP', 'FHRR', 'HRR')."""
        pass

    def __repr__(self) -> str:
        return (f"{self.__class__.__name__}(dimension={self.dimension}, "
                f"space={self.space.space_name}, backend={self.backend.name})")

class FHRRModel(VSAModel):
    """FHRR (Fourier HRR) model using complex phasors.

    Binding: element-wise complex multiplication (phase addition)
    Unbinding: element-wise multiplication with conjugate (phase subtraction)
    Bundling: element-wise addition + normalization to unit magnitude
    Permutation: circular shift (can also use phase rotation)

    Uses ComplexSpace with unit-magnitude phasors.
    """

    def __init__(
        self,
        dimension: int = 512,
        space: VectorSpace | None = None,
        backend: Backend | None = None,
        seed: int | None = None
    ):
        """Initialize FHRR model.

        Args:
            dimension: Dimensionality of hypervectors
                      (can be smaller than MAP due to better capacity)
            space: Vector space (defaults to ComplexSpace)
            backend: Computational backend
            seed: Random seed for space
        """
        if space is None:
            from ..backends import get_backend
            backend = backend if backend is not None else get_backend()
            space = ComplexSpace(dimension, backend=backend, seed=seed)

        super().__init__(space, backend)

    @property
    def model_name(self) -> str:
        return "FHRR"

    @property
    def is_self_inverse(self) -> bool:
        return False  # Requires conjugate, not same operation

    @property
    def is_commutative(self) -> bool:
        return True  # Complex multiplication is commutative

    @property
    def is_exact_inverse(self) -> bool:
        return True  # Conjugate provides exact inverse

    def bind(self, a: Array, b: Array) -> Array:
        """Bind using element-wise complex multiplication.

        For unit phasors: (a * b)[i] = a[i] * b[i]
        This adds phase angles: ∠(a*b) = ∠a + ∠b

        Args:
            a: First vector (unit phasors)
            b: Second vector (unit phasors)

        Returns:
            Bound vector c = a ⊙ b (element-wise product)
        """
        result = self.backend.multiply(a, b)
        # Normalize to unit magnitude
        return self.normalize(result)

    def unbind(self, a: Array, b: Array) -> Array:
        """Unbind using element-wise multiplication with conjugate.

        To recover original from c = a ⊙ b:
        unbind(c, b) = c ⊙ b* = (a ⊙ b) ⊙ b* = a ⊙ (b ⊙ b*) = a ⊙ 1 = a

        Args:
            a: Bound vector (or first operand)
            b: Second operand

        Returns:
            Unbound vector (exact recovery)
        """
        b_conj = self.backend.conjugate(b)
        result = self.backend.multiply(a, b_conj)
        return self.normalize(result)

    def bundle(self, vectors: Sequence[Array]) -> Array:
        """Bundle using element-wise addition.

        Sum phasors and normalize back to unit magnitude.
        The result points in the "average" direction of inputs.

        Args:
            vectors: Sequence of vectors to bundle

        Returns:
            Bundled vector (normalized to unit magnitude)

        Raises:
            ValueError: If vectors is empty
        """
        if not vectors:
            raise ValueError("Cannot bundle empty sequence")

        vectors = list(vectors)

        # Sum all vectors (phasors add vectorially)
        result = self.backend.sum(self.backend.stack(vectors, axis=0), axis=0)

        # Normalize to unit magnitude
        return self.normalize(result)

    def permute(self, vec: Array, k: int = 1) -> Array:
        """Permute using circular shift.

        For FHRR, permutation can be done as:
        1. Circular shift (coordinate permutation)
        2. Phase rotation (multiply by exp(i*2πk/D))

        We use circular shift for consistency with other models.

        Args:
            vec: Vector to permute
            k: Number of positions to shift

        Returns:
            Permuted vector
        """
        return self.backend.roll(vec, shift=k)

    def fractional_power(self, vec: Array, exponent: float) -> Array:
        """Raise phasor to a fractional power.

        For unit phasor z = exp(iθ): z^α = exp(iαθ)
        This is useful for encoding continuous values.

        Args:
            vec: Vector of unit phasors
            exponent: Power to raise to

        Returns:
            Vector with phases scaled by exponent

        Example:
            >>> base = model.random()
            >>> # Encode value 2.5 using fractional power
            >>> encoded = model.fractional_power(base, 2.5)
        """
        # For unit phasors z = exp(iθ), we want: z^α = exp(iαθ)
        # This is exact and avoids branch cuts from complex logarithms.
        #
        # Implementation:
        # 1. Extract phase θ = arg(z)
        # 2. Scale by exponent: αθ
        # 3. Create new phasor: exp(iαθ)

        # Get phase angles using backend operation
        angles = self.backend.angle(vec)

        # Scale angles by exponent
        scaled_angles = self.backend.multiply_scalar(angles, exponent)

        # Create new phasors: exp(i * scaled_angles)
        # exp(iθ) = cos(θ) + i*sin(θ)
        result = self.backend.exp(1j * scaled_angles)

        # Renormalize to unit magnitude (handle numerical errors)
        return self.normalize(result)

    def __repr__(self) -> str:
        return (f"FHRRModel(dimension={self.dimension}, "
                f"space={self.space.space_name}, "
                f"backend={self.backend.name})")

class GHRRModel(VSAModel):
    """GHRR (Generalized Holographic Reduced Representations) model.

    Binding: element-wise matrix multiplication (phase addition per matrix)
    Unbinding: element-wise multiplication with conjugate transpose
    Bundling: element-wise addition + normalization
    Permutation: circular shift (or use non-commutativity instead)

    Uses MatrixSpace with m×m unitary matrices.
    """

    def __init__(
        self,
        dimension: int = 100,
        matrix_size: int = 3,
        space: VectorSpace | None = None,
        backend: Backend | None = None,
        seed: int | None = None,
        diagonality: float | None = None
    ):
        """Initialize GHRR model.

        Args:
            dimension: Number of matrices in hypervector (can be smaller than
                      scalar models due to better capacity)
            matrix_size: Size m of each m×m matrix (default: 3)
                        Larger m → more non-commutative, better for complex structures
                        m=1 recovers FHRR
            space: Vector space (defaults to MatrixSpace)
            backend: Computational backend
            seed: Random seed for space
            diagonality: Control commutativity in [0, 1]
                        None: Random (default)
                        0.0: Maximally non-commutative
                        1.0: Fully commutative (FHRR-like)
        """
        if space is None:
            from ..backends import get_backend
            backend = backend if backend is not None else get_backend()
            space = MatrixSpace(
                dimension,
                matrix_size=matrix_size,
                backend=backend,
                seed=seed,
                diagonality=diagonality,
            )

        super().__init__(space, backend)

        # Store matrix size for easy access
        self.matrix_size = matrix_size if isinstance(space, MatrixSpace) else 1
        self._diagonality = diagonality

    @property
    def model_name(self) -> str:
        return f"GHRR_m{self.matrix_size}"

    @property
    def is_self_inverse(self) -> bool:
        return False  # Requires conjugate transpose

    @property
    def is_commutative(self) -> bool:
        return False  # Matrix multiplication is non-commutative

    @property
    def is_exact_inverse(self) -> bool:
        return True  # Conjugate transpose provides exact inverse

    @property
    def commutativity_degree(self) -> float:
        """Degree of commutativity in [0, 1].

        For GHRR, this depends on the diagonality of Q matrices.
        More diagonal → more commutative.

        Returns:
            0.0 if maximally non-commutative, 1.0 if fully commutative
        """
        if self._diagonality is not None:
            return self._diagonality

        # For random GHRR, larger m tends toward lower diagonality
        # This is approximate based on Yeung et al. Figure 6
        if self.matrix_size == 1:
            return 1.0  # FHRR is commutative
        elif self.matrix_size == 2:
            return 0.7  # Mostly commutative
        elif self.matrix_size == 3:
            return 0.5  # Balanced
        else:
            return 0.3  # Mostly non-commutative

    def bind(self, a: Array, b: Array) -> Array:
        """Bind using element-wise matrix multiplication.

        For matrices at position j: (a ⊗ b)_j = a_j @ b_j

        This is non-commutative: a ⊗ b ≠ b ⊗ a in general.

        Args:
            a: First hypervector (D, m, m)
            b: Second hypervector (D, m, m)

        Returns:
            Bound hypervector c where c_j = a_j @ b_j for all j
        """
        # Element-wise matrix multiplication using matmul broadcast
        # For (D, m, m) @ (D, m, m), this does D separate m×m multiplications
        result = self.backend.matmul(a, b)

        # Normalization not strictly needed for unitary matrices
        # but helps with numerical stability
        return result

    def unbind(self, a: Array, b: Array) -> Array:
        """Unbind using element-wise multiplication with conjugate transpose.

        To recover original from c = a ⊗ b:
        unbind(c, b) = c_j @ b_j† for all j

        This provides exact recovery: unbind(bind(a, b), b) = a

        Args:
            a: Bound hypervector (or first operand)
            b: Second operand

        Returns:
            Unbound hypervector (exact recovery)
        """
        # Compute b^† (conjugate transpose of each matrix)
        b_conj_t = self.backend.conjugate(self.backend.matrix_transpose(b))

        # Element-wise matrix multiply
        result = self.backend.matmul(a, b_conj_t)

        return result

    def bundle(self, vectors: Sequence[Array]) -> Array:
        """Bundle using element-wise addition.

        Sum all hypervectors element-wise. Each element is an m×m matrix.

        For GHRR: (a + b)_j = a_j + b_j (matrix addition)

        Args:
            vectors: Sequence of hypervectors to bundle

        Returns:
            Bundled hypervector

        Raises:
            ValueError: If vectors is empty
        """
        if not vectors:
            raise ValueError("Cannot bundle empty sequence")

        vectors = list(vectors)

        # Sum all vectors (element-wise matrix addition)
        result = self.backend.sum(self.backend.stack(vectors, axis=0), axis=0)

        # Normalize to project back to unitary matrices
        # This is critical for maintaining quasi-orthogonality (Yeung et al. 2024)
        # Uses polar decomposition via SVD
        result = self.space.normalize(result)

        return result

    def permute(self, vec: Array, k: int = 1) -> Array:
        """Permute using circular shift.

        For GHRR, permutation is less critical since non-commutativity
        can encode order. But still useful for some applications.

        Args:
            vec: Hypervector to permute (D, m, m)
            k: Number of positions to shift

        Returns:
            Permuted hypervector
        """
        # Roll along first dimension (shift which matrix is at which position)
        return self.backend.roll(vec, shift=k, axis=0)

    def test_non_commutativity(self, a: Array, b: Array) -> float:
        """Test degree of non-commutativity for two hypervectors.

        Computes: δ(a ⊗ b, b ⊗ a)

        A similarity of 1.0 means commutative, close to 0 means non-commutative.

        Args:
            a: First hypervector
            b: Second hypervector

        Returns:
            Similarity between a⊗b and b⊗a
        """
        ab = self.bind(a, b)
        ba = self.bind(b, a)
        return self.similarity(ab, ba)

    def compute_diagonality(self, vec: Array) -> float:
        """Compute average diagonality of matrices in hypervector.

        Diagonality metric: Σ|Q_jj| / ΣΣ|Q_jk|

        Args:
            vec: Hypervector (D, m, m)

        Returns:
            Diagonality in [0, 1]
        """
        vec_np = self.backend.to_numpy(vec)

        D = vec_np.shape[0]
        m = vec_np.shape[1]

        total_diag = 0.0
        total_all = 0.0

        for i in range(D):
            matrix = vec_np[i]
            # Diagonal sum
            diag_sum = np.sum(np.abs(np.diag(matrix)))
            # Total sum
            all_sum = np.sum(np.abs(matrix))

            total_diag += diag_sum
            total_all += all_sum

        return total_diag / total_all if total_all > 0 else 0.0

    def __repr__(self) -> str:
        return (f"GHRRModel(dimension={self.dimension}, "
                f"matrix_size={self.matrix_size}, "
                f"space={self.space.space_name}, "
                f"backend={self.backend.name})")

class MAPModel(VSAModel):
    """MAP (Multiply-Add-Permute) model.

    Binding: element-wise multiplication
    Unbinding: element-wise multiplication (self-inverse)
    Bundling: element-wise addition + normalization
    Permutation: circular shift

    Best used with BipolarSpace or RealSpace.
    """

    def __init__(
        self,
        dimension: int = 10000,
        space: VectorSpace | None = None,
        backend: Backend | None = None,
        seed: int | None = None,
    ):
        """Initialize MAP model.

        Args:
            dimension: Dimensionality of hypervectors
            space: Vector space (defaults to BipolarSpace)
            backend: Computational backend
            seed: Random seed for space
        """
        if space is None:
            from ..backends import get_backend

            backend = backend if backend is not None else get_backend()
            space = BipolarSpace(dimension, backend=backend, seed=seed)

        super().__init__(space, backend)

        # Pre-compute permutation indices for efficiency
        self._permutation_indices = list(range(self.dimension))

    @property
    def model_name(self) -> str:
        return "MAP"

    @property
    def is_self_inverse(self) -> bool:
        return True

    @property
    def is_commutative(self) -> bool:
        return True

    @property
    def is_exact_inverse(self) -> bool:
        # Exact for bipolar, approximate for continuous
        return self.space.space_name == "bipolar"

    def bind(self, a: Array, b: Array) -> Array:
        """Bind using element-wise multiplication.

        For bipolar: XOR when represented as {0,1}
        For real: Hadamard product

        Args:
            a: First vector
            b: Second vector

        Returns:
            Bound vector c = a ⊙ b
        """
        result = self.backend.multiply(a, b)
        # Normalize to maintain unit norm for continuous spaces
        if self.space.space_name != "bipolar":
            result = self.normalize(result)
        return result

    def unbind(self, a: Array, b: Array) -> Array:
        """Unbind using element-wise multiplication (self-inverse).

        Since binding is self-inverse: unbind(c, b) = c ⊙ b

        Args:
            a: Bound vector (or first operand)
            b: Second operand

        Returns:
            Unbound vector (exact for bipolar, approximate for continuous)
        """
        # For MAP, binding = unbinding
        return self.bind(a, b)

    def bundle(self, vectors: Sequence[Array]) -> Array:
        """Bundle using element-wise addition.

        For bipolar: majority vote after summing
        For real: sum and normalize

        Args:
            vectors: Sequence of vectors to bundle

        Returns:
            Bundled vector

        Raises:
            ValueError: If vectors is empty
        """
        if not vectors:
            raise ValueError("Cannot bundle empty sequence")

        vectors = list(vectors)

        # Sum all vectors
        result = self.backend.sum(self.backend.stack(vectors, axis=0), axis=0)

        # Normalize according to space
        if self.space.space_name == "bipolar":
            # Majority vote: sign of sum
            result = self.backend.sign(result)
            # Handle zeros (shouldn't happen in practice, but be safe)
            # If sum is 0, randomly choose ±1
            zeros_mask = result == 0
            if self.backend.to_numpy(zeros_mask).any():
                # For any zeros, use the first vector's value
                first_vec = vectors[0]
                result = self.backend.where(zeros_mask, first_vec, result)
        else:
            # For continuous spaces, L2 normalize
            result = self.normalize(result)

        return result

    def permute(self, vec: Array, k: int = 1) -> Array:
        """Permute using circular shift.

        Shifts vector elements by k positions to the right.
        Negative k shifts left.

        Args:
            vec: Vector to permute
            k: Number of positions to shift

        Returns:
            Permuted vector
        """
        return self.backend.roll(vec, shift=k)

    def __repr__(self) -> str:
        return (
            f"MAPModel(dimension={self.dimension}, "
            f"space={self.space.space_name}, "
            f"backend={self.backend.name})"
        )

class HRRModel(VSAModel):
    """HRR (Holographic Reduced Representations) model.

    Binding: circular convolution (via FFT)
    Unbinding: circular correlation (via FFT)
    Bundling: element-wise addition + normalization
    Permutation: circular shift

    Uses RealSpace with Gaussian distribution N(0, 1/D).
    """

    def __init__(
        self,
        dimension: int = 10000,
        space: VectorSpace | None = None,
        backend: Backend | None = None,
        seed: int | None = None,
    ):
        """Initialize HRR model.

        Args:
            dimension: Dimensionality of hypervectors (recommend 1000-10000)
            space: Vector space (defaults to RealSpace)
            backend: Computational backend
            seed: Random seed for space
        """
        if space is None:
            from ..backends import get_backend

            backend = backend if backend is not None else get_backend()
            space = RealSpace(dimension, backend=backend, seed=seed)

        super().__init__(space, backend)

    @property
    def model_name(self) -> str:
        return "HRR"

    @property
    def is_self_inverse(self) -> bool:
        return False  # Requires correlation, not same operation

    @property
    def is_commutative(self) -> bool:
        return True  # Convolution is commutative

    @property
    def is_exact_inverse(self) -> bool:
        return False  # Correlation gives approximate inverse

    def bind(self, a: Array, b: Array) -> Array:
        """Bind using circular convolution.

        Implemented via FFT: conv(a, b) = IFFT(FFT(a) * FFT(b))

        Args:
            a: First vector
            b: Second vector

        Returns:
            Bound vector c = a ⊛ b (circular convolution)
        """
        # Circular convolution in frequency domain
        result = self.backend.circular_convolve(a, b)

        # Do NOT normalize - preserves magnitude for proper unbinding via
        # circular correlation. Normalization would interfere with the
        # mathematical relationship required for unbind recovery.
        return result

    def unbind(self, a: Array, b: Array) -> Array:
        """Unbind using circular correlation (approximate inverse of convolution).

        This is the classic HRR unbinding operation that uses circular correlation
        to approximately recover the original vector from a bound pair.

        Args:
            a: Bound vector c = x ⊛ b (result of circular convolution)
            b: Key vector (second operand in binding)

        Returns:
            Approximate recovery of x (original vector), normalized to unit length

        Notes
        -----
        **Mathematical Foundation:**

        HRR binding via circular convolution:
            c = x ⊛ b

        In frequency domain (Fourier):
            C(ω) = X(ω) · B(ω)

        Unbinding via circular correlation:
            x̂ = c ⋆ b = IFFT(C(ω) · B*(ω))

        Where B*(ω) is the complex conjugate of B(ω).

        Substituting C(ω) = X(ω) · B(ω):
            x̂ = IFFT(X(ω) · B(ω) · B*(ω))
              = IFFT(X(ω) · |B(ω)|²)

        For random vectors with approximately uniform power spectrum (|B(ω)|² ≈ 1),
        this gives x̂ ≈ x.

        **Approximation Quality:**

        Recovery similarity depends on:
        - Dimension D: Higher D → better recovery
        - Noise level: Clean binding → better unbind
        - Bundle size: More items → more interference

        Empirical performance (D=10000):
        - Clean unbind: similarity ≈ 0.70-0.72 (approximate inverse)
        - After bundling 2 items: similarity ≈ 0.57
        - After bundling 10 items: similarity ≈ 0.30
        - After bundling 100 items: similarity decreases further

        Note: Unlike FHRR which achieves exact (1.0) recovery, HRR's circular
        correlation provides only approximate recovery. The ~0.71 similarity
        is sufficient for retrieval tasks but requires cleanup/thresholding.

        References
        ----------
        - Plate (1995): "Holographic Reduced Representations"
        - Plate (2003): "Holographic Reduced Representations" (full book)

        Examples
        --------
        >>> model = VSA.create('HRR', dim=10000)
        >>> x = model.random(seed=1)
        >>> b = model.random(seed=2)
        >>> c = model.bind(x, b)
        >>> x_recovered = model.unbind(c, b)
        >>> similarity = model.similarity(x, x_recovered)
        >>> print(f"Recovery similarity: {similarity:.3f}")  # ~0.71
        """
        # Transform to frequency domain
        fa = self.backend.fft(a)
        fb = self.backend.fft(b)

        # Circular correlation: C(ω) * conj(B(ω))
        # This is the classic HRR unbinding operation (Plate, 1995)
        fr = self.backend.multiply(fa, self.backend.conjugate(fb))

        # Transform back to time domain
        time = self.backend.ifft(fr)

        # Take real part (imaginary part should be near zero due to real inputs)
        result = self.backend.real(time)

        # Normalize to unit length for consistent comparison with other vectors
        return self.normalize(result)

    def bundle(self, vectors: Sequence[Array]) -> Array:
        """Bundle using element-wise addition (superposition).

        For HRR, bundling is simple vector addition without normalization.
        This preserves the magnitude relationships needed for proper unbinding.

        Args:
            vectors: Sequence of vectors to bundle

        Returns:
            Bundled vector (unnormalized sum)

        Raises:
            ValueError: If vectors is empty

        Notes:
            Unlike some VSA models, HRR does NOT normalize after bundling.
            Normalization would interfere with the circular correlation unbinding
            operation. The unbind() method handles normalization of its output.
        """
        if not vectors:
            raise ValueError("Cannot bundle empty sequence")

        vectors = list(vectors)

        # Sum all vectors (simple superposition, no normalization)
        result = self.backend.sum(self.backend.stack(vectors, axis=0), axis=0)

        return result

    def permute(self, vec: Array, k: int = 1) -> Array:
        """Permute using circular shift.

        Shifts vector elements by k positions to the right.
        Negative k shifts left.

        Args:
            vec: Vector to permute
            k: Number of positions to shift

        Returns:
            Permuted vector
        """
        return self.backend.roll(vec, shift=k)

    def __repr__(self) -> str:
        return (
            f"HRRModel(dimension={self.dimension}, "
            f"space={self.space.space_name}, "
            f"backend={self.backend.name})"
        )