# Copyright 2025 Polyquantique
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
"""
Permutation matrices Python interface
References
----------
[1] Collins, B., & Nagatsu, M. (2025). Weingarten calculus for centered random permutation
matrices. arXiv preprint arXiv:2503.18453.
"""
from math import factorial, prod
from functools import lru_cache
from itertools import product
from collections.abc import Sequence
from fractions import Fraction
from sympy import Expr, Symbol, binomial
from haarpy import set_partitions, meet_operation, join_operation, partial_order
from ._utils import _simplify
[docs]
@lru_cache
def mobius_function(
partition_1: tuple[tuple[int, ...], ...], partition_2: tuple[tuple[int, ...], ...]
) -> int:
"""Returns the Möbius function of two given partitions
Parameters
----------
partition_1 : tuple[tuple[int]
The intersected partition
partition_2 : tuple[tuple[int]]
The partition summed over
Returns
-------
int
The value of the Möbius function
Examples
--------
>>> from haarpy import mobius_function
>>> partition_1 = ((0, 3), (1, 2), (4,), (5,))
>>> partition_2 = ((0, 1, 2), (3, 4, 5))
>>> mobius_function(partition_1, partition_2)
-2
"""
partition_set_1 = tuple(set(block) for block in partition_1)
partition_set_2 = tuple(set(block) for block in partition_2)
partition_intersection = tuple(
sum(1 for block_1 in partition_set_1 if block_1 & block_2) for block_2 in partition_set_2
)
return prod(
(-1) ** (block_count - 1) * factorial(block_count - 1)
for block_count in partition_intersection
)
[docs]
@lru_cache
def weingarten_permutation(
first_partition: tuple[tuple[int, ...], ...],
second_partition: tuple[tuple[int, ...], ...],
dimension: Symbol,
) -> Expr:
"""Returns the Weingarten function for random permutation matrices
Parameters
----------
first_partition : tuple[tuple[int, ...], ...]
A partition of the set of :math:`k` integers :math:`[k] = \{1,2,\dots, k\}`
second_partition : tuple[tuple[int, ...], ...]
A second partition of the same set
dimension : Symbol
The dimension of the random permutation matrices
Returns
-------
Expr
The Weingarten function
Example
-------
>>> from sympy import Symbol
>>> from haarpy import weingarten_permutation
>>> d = Symbol('d')
>>> partition_1 = ((0, 1, 2), (3,))
>>> partition_2 = ((0,), (1, 2, 3))
>>> weingarten_permutation(partition_1, partition_2, 4)
Fraction(5, 24)
>>> weingarten_permutation(partition_1, partition_2, d)
(d + 1)/(d*(d - 3)*(d - 2)*(d - 1))
See Also
--------
:func:`haarpy.partition.meet_operation`
Returns the greatest lower bound of the two input partitions
:func:`haarpy.permutation.mobius_function`
Returns the Möbius function of two given partitions
"""
disjoint_partition_tuple = tuple(
(partition for partition in set_partitions(block))
for block in meet_operation(first_partition, second_partition)
)
inferieur_partition_tuple = (
tuple(block for partition in partition_tuple for block in partition)
for partition_tuple in product(*disjoint_partition_tuple)
)
if isinstance(dimension, int):
weingarten = sum(
Fraction(
mobius_function(partition, first_partition)
* mobius_function(partition, second_partition),
prod(dimension - i for i, _ in enumerate(partition)),
)
for partition in inferieur_partition_tuple
)
else:
weingarten_gen = (
mobius_function(partition, first_partition)
* mobius_function(partition, second_partition)
/ prod(dimension - i for i, _ in enumerate(partition))
for partition in inferieur_partition_tuple
)
weingarten = _simplify(weingarten_gen)
return weingarten
[docs]
@lru_cache
def weingarten_centered_permutation(
first_partition: tuple[tuple[int, ...], ...],
second_partition: tuple[tuple[int, ...], ...],
dimension: Symbol,
) -> Expr:
"""Returns the Weingarten function for centered random permutation matrices
Parameters
----------
first_partition : tuple[tuple[int, ...], ...]
A partition of the set of :math:`k` integers :math:`[k] = \{1,2,\dots, k\}`
second_partition : tuple[tuple[int, ...], ...]
A second partition of the same set
dimension : Symbol
The dimension of the random permutation matrices
Returns
-------
Expr
The Weingarten function
Example
-------
>>> from sympy import Symbol
>>> from haarpy import weingarten_centered_permutation
>>> d = Symbol('d')
>>> partition_1 = ((0,), (1,), (2,))
>>> partition_2 = ((0, 2), (1,))
>>> weingarten_centered_permutation(partition_1, partition_2, 3)
Fraction(-1, 9)
>>> weingarten_centered_permutation(partition_1, partition_2, d)
-2/(d**2*(d - 2)*(d - 1))
See Also
--------
:func:`haarpy.partition.meet_operation`
Returns the greatest lower bound of the two input partitions
:func:`haarpy.partition.join_operation`
Returns the least upper bound of the two input partitions
:func:`haarpy.permutation.mobius_function`
Returns the Möbius function of two given partitions
"""
singleton_set_size = sum(
1 for block in join_operation(first_partition, second_partition) if len(block) == 1
)
disjoint_partition_tuple = tuple(
(partition for partition in set_partitions(block))
for block in meet_operation(first_partition, second_partition)
)
inferieur_partition_tuple = tuple(
tuple(block for partition in partition_tuple for block in partition)
for partition_tuple in product(*disjoint_partition_tuple)
)
if isinstance(dimension, int):
weingarten = sum(
(-1) ** i
* Fraction(binomial(singleton_set_size, i))
* sum(
Fraction(
mobius_function(partition, first_partition)
* mobius_function(partition, second_partition),
dimension**i * prod(dimension - j for j in range(len(partition) - i)),
)
for partition in inferieur_partition_tuple
)
for i in range(singleton_set_size + 1)
)
else:
weingarten_gen = (
(-1) ** i
* binomial(singleton_set_size, i)
* sum(
mobius_function(partition, first_partition)
* mobius_function(partition, second_partition)
/ dimension**i
/ prod(dimension - j for j in range(len(partition) - i))
for partition in inferieur_partition_tuple
)
for i in range(singleton_set_size + 1)
)
weingarten = _simplify(weingarten_gen)
return weingarten
[docs]
@lru_cache
def haar_integral_permutation(
row_indices: tuple[int, ...],
column_indices: tuple[int, ...],
dimension: Symbol,
) -> Expr:
"""Returns the integral over Haar random permutation matrices
Parameters
----------
row_indices : tuple[int, ...]
A sequence of row indices
column_indices : tuple[int, ...]
A sequence of column indices
Returns
-------
Expr
The integral under the Haar measure
Raises
------
TypeError
If the input parameters are not of type ``Sequence``
ValueError
If the input sequences are of different length
Examples
--------
>>> from sympy import Symbol
>>> from haarpy import haar_integral_permutation
>>> d = Symbol('d')
>>> row_indices = (0, 1, 2, 2)
>>> column_indices = (1, 0, 2, 2)
>>> haar_integral_permutation(row_indices, column_indices, 3)
Fraction(1, 6)
>>> haar_integral_permutation(row_indices, column_indices, d)
1/(d*(d - 2)*(d - 1))
See Also
--------
:func:`haarpy.permutation.weingarten_permutation`
Returns the Weingarten function for random permutation matrices
"""
if not (isinstance(row_indices, Sequence) and isinstance(column_indices, Sequence)):
raise TypeError
if len(row_indices) != len(column_indices):
raise ValueError("Wrong tuple format")
def sequence_to_partition(sequence: tuple) -> tuple[tuple[int]]:
return sorted(
sorted(index for index, value in enumerate(sequence) if value == unique)
for unique in set(sequence)
)
row_partition = sequence_to_partition(row_indices)
column_partition = sequence_to_partition(column_indices)
return (
Fraction(1, prod(dimension - i for i, _ in enumerate(row_partition)))
if row_partition == column_partition and isinstance(dimension, int)
else (
1 / prod(dimension - i for i, _ in enumerate(row_partition))
if row_partition == column_partition
else 0
)
)
[docs]
@lru_cache
def haar_integral_centered_permutation(
row_indices: tuple[int, ...],
column_indices: tuple[int, ...],
dimension: Symbol,
) -> Expr:
"""Returns the integral over Haar random centered permutation matrices
Parameters
----------
row_indices : tuple[int, ...]
A sequence of row indices
column_indices : tuple[int, ...]
A sequence of column indices
Returns
-------
Expr
The integral under the Haar measure
Raises
------
TypeError
If the input parameters are not of type ``Sequence``
ValueError
If the input sequences are of different length
Examples
--------
>>> from sympy import Symbol
>>> from haarpy import haar_integral_centered_permutation
>>> d = Symbol('d')
>>> row_indices = (0, 1, 2)
>>> column_indices = (0, 1, 1)
>>> haar_integral_centered_permutation(row_indices, column_indices, 3)
Fraction(-1, 27)
>>> haar_integral_centered_permutation(row_indices, column_indices, d)
-2/(d**3*(d - 1))
See Also
--------
:func:`haarpy.partition.partial_order`
Compares two partitions in terms of partial order
:func:`haarpy.permutation.weingarten_centered_permutation`
Returns the Weingarten function for centered random permutation matrices
"""
if not (isinstance(row_indices, Sequence) and isinstance(column_indices, Sequence)):
raise TypeError
if len(row_indices) != len(column_indices):
raise ValueError("Wrong tuple format")
row_partition = tuple(
tuple(index for index, value in enumerate(row_indices) if value == unique)
for unique in set(row_indices)
)
column_partition = tuple(
tuple(index for index, value in enumerate(column_indices) if value == unique)
for unique in set(column_indices)
)
integral_gen = (
weingarten_centered_permutation(
partition_sigma,
partition_tau,
dimension,
)
for partition_sigma, partition_tau in product(
set_partitions(tuple(i for i, _ in enumerate(row_indices))),
set_partitions(tuple(i for i, _ in enumerate(column_indices))),
)
if partial_order(partition_sigma, row_partition)
and partial_order(partition_tau, column_partition)
)
return sum(integral_gen) if isinstance(dimension, int) else _simplify(integral_gen)