geometry_tools.representation

   1"""Work with group representations into finite-dimensional vector
   2spaces, using numerical matrices.
   3
   4
   5To make a representation, instantiate the `Representation` class, and
   6assign numpy arrays to group generators (which are lowercase letters
   7a-z). Use square brackets to access the image of a word in the
   8generators.
   9
  10
  11```python
  12import numpy as np
  13
  14from geometry_tools import representation
  15from geometry_tools.automata import fsa
  16
  17rep = representation.Representation()
  18rep["a"] = np.array([
  19    [3.0, 0.0],
  20    [0.0, 1/3.0]
  21])
  22
  23rep["aaa"]
  24```
  25
  26
  27
  28
  29    array([[27.        ,  0.        ],
  30           [ 0.        ,  0.03703704]])
  31
  32
  33
  34Generator inverses are automatically assigned to capital letters:
  35
  36
  37```python
  38rep = representation.Representation()
  39rep["a"] = np.array([
  40    [3.0, 0.0],
  41    [0.0, 1/3.0]
  42])
  43
  44rep["aA"]
  45```
  46
  47
  48
  49
  50    array([[1., 0.],
  51           [0., 1.]])
  52
  53
  54
  55A common use-case for this class is to get a list of matrices
  56representing all elements in the group, up to a bounded word
  57length. The fastest way to do this is to use the built-in
  58`Representation.freely_reduced_elements` method, which returns a numpy
  59array containing one matrix for each freely reduced word in the group
  60(up to a specified word length).
  61
  62The array of matrices is *not* typically ordered lexicographically. To
  63get a list of words corresponding to the matrices returned, pass the
  64`with_words` flag when calling
  65`Representation.freely_reduced_elements` (see the documentation for
  66that function for details).
  67
  68
  69```python
  70rep = representation.Representation()
  71rep["a"] = np.array([
  72    [3.0, 0.0],
  73    [0.0, 1/3.0]
  74])
  75rep["b"] = np.array([
  76    [1.0, -1.0],
  77    [1.0, 1.0]
  78]) / np.sqrt(2)
  79
  80
  81rep.freely_reduced_elements(6)
  82```
  83
  84
  85
  86
  87    array([[[ 1.00000000e+00,  0.00000000e+00],
  88            [ 0.00000000e+00,  1.00000000e+00]],
  89
  90           [[ 7.07106781e-01,  7.07106781e-01],
  91            [-7.07106781e-01,  7.07106781e-01]],
  92
  93           [[ 2.12132034e+00,  2.12132034e+00],
  94            [-2.35702260e-01,  2.35702260e-01]],
  95
  96           ...,
  97
  98           [[-2.12132034e+00, -2.12132034e+00],
  99            [ 2.35702260e-01, -2.35702260e-01]],
 100
 101           [[-2.35702260e-01, -2.35702260e-01],
 102            [ 2.12132034e+00, -2.12132034e+00]],
 103
 104           [[-1.07699575e-16, -1.00000000e+00],
 105            [ 1.00000000e+00,  7.51818954e-17]]])
 106
 107
 108
 109You can speed up this process even more if you have access to a finite
 110state automaton which provides a unique word for each element in your
 111group.
 112
 113For instance, to find the image of a Cayley ball of radius 10 under
 114the canonical representation of a (3,3,4) triangle group, you can use
 115the `Representation.automaton_accepted` method as follows:
 116
 117
 118```python
 119from geometry_tools import coxeter
 120from geometry_tools.automata import fsa
 121
 122# create the representation and load the built-in automaton
 123rep = coxeter.TriangleGroup((3,3,4)).canonical_representation()
 124automaton = fsa.load_builtin('cox334.wa')
 125
 126rep.automaton_accepted(automaton, 10)
 127```
 128
 129
 130
 131
 132    array([[[ 1.00000000e+00,  0.00000000e+00,  0.00000000e+00],
 133            [ 0.00000000e+00,  1.00000000e+00,  0.00000000e+00],
 134            [ 0.00000000e+00,  0.00000000e+00,  1.00000000e+00]],
 135
 136           [[ 4.44089210e-16, -1.00000000e+00,  2.41421356e+00],
 137            [ 1.00000000e+00, -1.00000000e+00,  1.00000000e+00],
 138            [ 0.00000000e+00,  0.00000000e+00,  1.00000000e+00]],
 139
 140           [[-1.00000000e+00,  1.00000000e+00,  1.41421356e+00],
 141            [ 0.00000000e+00,  1.00000000e+00,  0.00000000e+00],
 142            [ 0.00000000e+00,  0.00000000e+00,  1.00000000e+00]],
 143
 144           ...,
 145
 146           [[-2.47279221e+01,  2.57279221e+01,  2.23137085e+01],
 147            [-2.91421356e+01,  2.91421356e+01,  2.71421356e+01],
 148            [-1.50710678e+01,  1.50710678e+01,  1.40710678e+01]],
 149
 150           [[-7.24264069e+00,  7.24264069e+00,  6.82842712e+00],
 151            [-1.06568542e+01,  1.16568542e+01,  9.24264069e+00],
 152            [-5.82842712e+00,  6.82842712e+00,  4.82842712e+00]],
 153
 154           [[-2.47279221e+01,  2.57279221e+01,  2.23137085e+01],
 155            [-1.06568542e+01,  1.16568542e+01,  9.24264069e+00],
 156            [-3.05563492e+01,  3.29705627e+01,  2.67279221e+01]]])
 157
 158    """
 159
 160import re
 161import itertools
 162
 163import numpy as np
 164import scipy.special
 165
 166from .automata import fsa
 167from . import lie
 168from . import utils
 169from .utils import words
 170
 171if utils.SAGE_AVAILABLE:
 172    from .utils import sagewrap
 173
 174class Representation:
 175    """Model a representation for a finitely generated group
 176    representation into GL(n).
 177    """
 178
 179    @property
 180    def dim(self):
 181        return self._dim
 182
 183    @property
 184    def dtype(self):
 185        return self._dtype
 186
 187    @property
 188    def base_ring(self):
 189        return self._base_ring
 190
 191    @staticmethod
 192    def wrap_func(numpy_matrix):
 193        return numpy_matrix
 194
 195    @staticmethod
 196    def unwrap_func(wrapped_matrix):
 197        return wrapped_matrix
 198
 199    @staticmethod
 200    def array_wrap_func(numpy_array):
 201        return numpy_array
 202
 203    def parse_word(self, word, simple=None):
 204        if simple is None:
 205            simple = self.parse_simple
 206
 207        if simple:
 208            return word
 209
 210        return re.split("[()*]", word)
 211
 212    def __getitem__(self, word):
 213        return self.element(word)
 214
 215    def element(self, word, parse_simple=True):
 216        matrix = self._word_value(word, parse_simple)
 217        return self.__class__.wrap_func(matrix)
 218
 219    def __init__(self, representation=None,
 220                 generator_names=None,
 221                 relations=[],
 222                 base_ring=None,
 223                 parse_simple=True,
 224                 dtype='float64'):
 225        """
 226        Parameters
 227        ----------
 228        representation : Representation
 229            Representation to copy elements from
 230        generator_names : iterable of strings
 231            Names to use for the generators. These must be initialized
 232            as arrays later to use the representation properly.
 233        dtype : numpy dtype
 234            Default data type for matrix entries. Defaults to 'float64'.
 235
 236        """
 237
 238        self._dim = None
 239        self._dtype = dtype
 240        self._base_ring = base_ring
 241        self.parse_simple = parse_simple
 242        self.relations = [r for r in relations]
 243
 244        if representation is not None:
 245            if base_ring is None:
 246                self._base_ring = representation.base_ring
 247
 248            if generator_names is None:
 249                generator_names = list(representation.generators)
 250
 251            self.generators = {}
 252
 253            # don't recompute inverses for an existing representation
 254            for gen in generator_names:
 255                self._set_generator(gen, representation.generators[gen],
 256                                   compute_inverse=False)
 257
 258            self._dim = representation._dim
 259
 260            self.relations += [r for r in representation.relations]
 261
 262        else:
 263            if generator_names is None:
 264                generator_names = []
 265
 266            self.generators = {name[0].lower():None
 267                               for name in words.asym_gens(generator_names)}
 268
 269            for gen in list(self.generators):
 270                self.generators[gen.upper()] = None
 271
 272    def freely_reduced_elements(self, length, maxlen=True,
 273                                with_words=False):
 274        """Return group elements, one for each freely reduced word
 275           in the generators.
 276
 277        Parameters
 278        ----------
 279        length : int
 280            maximum length of a word to find a representation of
 281        maxlen : bool
 282            If `True` (the default), the returned array has one
 283            element for each word up to length `length`
 284            (inclusive). If `False`, only compute the image of words
 285            of length exactly equal to `length`.
 286        with_words : bool
 287            If `True`, also return the list of words corresponding to
 288            the computed matrices.
 289
 290        Returns
 291        -------
 292        result : ndarray or tuple
 293            If `with_words` is `True`, `result` is a tuple `(elements, words)`,
 294            where `elements` is an ndarray of shape `(k, n, n)`,
 295            containing one matrix for each of the `k` freely reduced
 296            words in the generators, and `words` is a python list of
 297            strings, containing those words. If `with_words` is
 298            `False`, just return `elements`.
 299
 300        """
 301        automaton = fsa.free_automaton(list(self.asym_gens()))
 302        return self.automaton_accepted(automaton, length,
 303                                       maxlen=maxlen,
 304                                       with_words=with_words)
 305
 306    def free_words_of_length(self, length):
 307        """Yield freely reduced words in the generators, of a specified length.
 308
 309        Parameters
 310        ----------
 311        length : int
 312            the length of words to return
 313
 314        Yields
 315        ------
 316        word : string
 317            All freely reduced words in the generators of length
 318            `length`.
 319
 320        """
 321
 322        if length == 0:
 323            yield ""
 324        else:
 325            for word in self.free_words_of_length(length - 1):
 326                for generator in self.generators:
 327                    if len(word) == 0 or generator != words.invert_gen(word[-1]):
 328                        yield word + generator
 329
 330    def free_words_less_than(self, length):
 331        """Yield freely reduced words in the generators, up to a specified
 332        length.
 333
 334        Parameters
 335        ----------
 336        length : int
 337            the maximum length of words to return
 338
 339        Yields
 340        ------
 341        word : string
 342            All freely reduced words in the generators, up to length
 343            `length` (inclusive)
 344
 345        """
 346        for i in range(length):
 347            for word in self.free_words_of_length(i):
 348                yield word
 349
 350    def automaton_accepted(self, automaton, length,
 351                           maxlen=True, with_words=False,
 352                           start_state=None, end_state=None,
 353                           precomputed=None, edge_words=False):
 354        """Return group elements representing words accepted by a
 355           finite-state automaton.
 356
 357        Parameters
 358        ----------
 359        automaton : automata.fsa.FSA
 360            Finite-state automaton which determines the accepted words
 361        length : int
 362            Maximum length of a word to compute
 363        maxlen : bool
 364            if `True` (the default), compute representations of all
 365            accepted words of length up to `length`. If `False`, only
 366            compute representations of words whose length is exactly
 367            `length`.
 368        with_words : bool
 369            Whether to return a list of accepted words along with the
 370            computed array of images.
 371        start_state: object
 372            vertex of the automaton to use as the starting state for
 373            all accepted words. If `None`, then the default start
 374            vertex of the automaton is used. *Note*: at most one of
 375            `start_state` and `end_state` can be specified.
 376        end_state: object
 377            vertex of the automaton where all accepted paths must
 378            end. If `None`, then any ending state is allowed. *Note*:
 379            at most one of `start_state` and `end_state` can be
 380            specified.
 381        precomputed: dict
 382            dictionary giving precomputed values of this
 383            function. Keys are tuples of then form `(length, state)`,
 384            where `length` is an integer and `state` is a vertex of
 385            the automaton. If `None`, use an empty dictionary. In
 386            either case, the dictionary will be populated when the
 387            function is called.
 388        edge_words : bool
 389            If True, view each label of the given automaton as a word
 390            in the generators for this representation. Otherwise,
 391            interpret the label as a single generator.
 392
 393        Returns
 394        -------
 395        result : ndarray or tuple
 396            If `with_words` is `True`, `result` is a tuple `(elements, words)`,
 397            where `elements` is an ndarray containing one
 398            matrix for each accepted word, and `words` is a list of
 399            strings containing the corresponding words. If
 400            `with_words` is `False`, just return `elements`.
 401
 402    """
 403
 404        if start_state is not None and end_state is not None:
 405            raise ValueError("At most one of start_state and end_state "
 406                             "can be specified")
 407
 408        as_start = True
 409        state = None
 410
 411        if start_state is not None:
 412            state = start_state
 413            as_start = True
 414
 415        if end_state is not None:
 416            state = end_state
 417            as_start = False
 418
 419        result = self._automaton_accepted(automaton, length,
 420                                          maxlen=maxlen,
 421                                          with_words=with_words,
 422                                          state=state,
 423                                          as_start=as_start,
 424                                          precomputed=precomputed)
 425
 426        if with_words:
 427            matrix_array, words = result
 428        else:
 429            matrix_array = result
 430
 431        wrapped_matrices = self.__class__.array_wrap_func(matrix_array)
 432
 433        if with_words:
 434            return wrapped_matrices, words
 435
 436        return wrapped_matrices
 437
 438    def _automaton_accepted(self, automaton, length,
 439                            state=None, as_start=True, maxlen=True,
 440                            precomputed=None, with_words=False,
 441                            edge_words=False):
 442
 443        if precomputed is None:
 444            precomputed = {}
 445
 446        if (length, state) in precomputed:
 447            return precomputed[(length, state)]
 448
 449        empty_arr = np.array([]).reshape((0, self.dim, self.dim))
 450
 451        if length == 0:
 452            if state is None or as_start or state in automaton.start_vertices:
 453                id_array = np.array([utils.identity(self.dim,
 454                                                    dtype=self.dtype,
 455                                                    base_ring=self.base_ring)])
 456                if with_words:
 457                    return (id_array, [""])
 458                return id_array
 459
 460            if with_words:
 461                return (empty_arr, [])
 462            return empty_arr
 463
 464        if state is None:
 465            as_start = True
 466            state = automaton.start_vertices[0]
 467
 468        if as_start:
 469            adj_states = automaton.out_dict[state]
 470        else:
 471            adj_states = automaton.in_dict[state]
 472
 473        if len(adj_states) == 0 and not as_start:
 474            if with_words:
 475                return (empty_arr, [])
 476            return empty_arr
 477
 478
 479        matrix_list = []
 480        accepted_words = []
 481        for adj_state, labels in adj_states.items():
 482            for label in labels:
 483                result = self._automaton_accepted(
 484                    automaton, length - 1,
 485                    state=adj_state,
 486                    as_start=as_start,
 487                    maxlen=maxlen,
 488                    precomputed=precomputed,
 489                    with_words=with_words
 490                )
 491                if with_words:
 492                    matrices, words = result
 493                    if as_start:
 494                        words = [label + word for word in words]
 495                    else:
 496                        words = [word + label for word in words]
 497                    accepted_words += words
 498                else:
 499                    matrices = result
 500
 501                if edge_words:
 502                    edge_elt = self._word_value(label)
 503                else:
 504                    edge_elt = self.generators[label]
 505                if as_start:
 506                    matrices = edge_elt @ matrices
 507                else:
 508                    matrices = matrices @ edge_elt
 509                matrix_list.append(matrices)
 510
 511        accepted_matrices = empty_arr
 512
 513        if len(matrix_list) > 0:
 514            accepted_matrices = np.concatenate(matrix_list)
 515
 516        if maxlen:
 517            additional_result = self._automaton_accepted(
 518                automaton, 0,
 519                state=state,
 520                as_start=as_start,
 521                with_words=with_words
 522            )
 523            if with_words:
 524                additional_mats, additional_words = additional_result
 525                accepted_words = additional_words + accepted_words
 526            else:
 527                additional_mats = additional_result
 528
 529            accepted_matrices = np.concatenate(
 530                [additional_mats, accepted_matrices]
 531            )
 532
 533        if with_words:
 534            accepted = (accepted_matrices, accepted_words)
 535        else:
 536            accepted = accepted_matrices
 537
 538        precomputed[(length, state)] = accepted
 539        return accepted
 540
 541    def elements(self, words):
 542        """Get images of an iterable of words.
 543
 544        Parameters
 545        ----------
 546        words : iterable of strings
 547            words to find the image of under the representation
 548
 549        Returns
 550        -------
 551        ndarray
 552            numpy array of shape `(l, n, n)`, where `l` is the length
 553            of `words` and `n` is the dimension of the representation.
 554
 555        """
 556        return self.__class__.array_wrap_func(np.array(
 557            [self._word_value(word) for word in words]
 558        ))
 559
 560    def asym_gens(self):
 561        """Iterate over asymmetric group generator names for this
 562           representation.
 563
 564        Only iterate over group generators (lowercase letters), not
 565        semigroup generators (upper and lowercase letters).
 566
 567        Yields
 568        ------
 569        gen : string
 570            group generator names
 571
 572    """
 573        return words.asym_gens(self.generators.keys())
 574
 575    def dual(self):
 576        return self._compose(
 577            lambda M: utils.invert(M).T
 578        )
 579
 580    def _word_value(self, word, parse_simple=None):
 581        gen_list = self.parse_word(word, simple=parse_simple)
 582
 583        matrix = utils.identity(self._dim, dtype=self.dtype)
 584        for gen in gen_list:
 585            matrix = matrix @ self.generators[gen]
 586        return matrix
 587
 588    def _set_generator(self, generator, matrix, compute_inverse=True,
 589                       base_ring=None):
 590
 591        shape = matrix.shape
 592
 593        if self._dim is None:
 594            self._dim = shape[0]
 595        if shape[0] != shape[1]:
 596            raise ValueError("Matrices representing group elements must be square")
 597        if shape[0] != self._dim:
 598            raise ValueError(
 599                "Every matrix in the representation must have the same shape"
 600            )
 601
 602        if re.search("[*()]", generator):
 603            raise ValueError(
 604                "The characters *, (, and ) are reserved and "
 605                " cannot be used in generator names."
 606            )
 607
 608        if not re.search("[a-zA-Z]", generator):
 609            raise ValueError(
 610                "Generator names must contain at least one letter a-z or A-Z"
 611            )
 612
 613        if generator != generator.lower() and generator != generator.upper():
 614            raise ValueError(
 615                "Generator names cannot be mixed-case"
 616            )
 617
 618        if base_ring is not None:
 619            matrix = utils.change_base_ring(matrix, base_ring)
 620            self._base_ring = base_ring
 621
 622        self.generators[generator] = matrix
 623        if compute_inverse:
 624            self.generators[words.invert_gen(generator)] = utils.invert(matrix)
 625
 626        # always update the dtype (we don't have a hierarchy for this)
 627        self._dtype = matrix.dtype
 628
 629    def set_generator(self, generator, matrix, **kwargs):
 630        self._set_generator(generator,
 631            self.__class__.unwrap_func(matrix),
 632            **kwargs)
 633
 634    def __setitem__(self, generator, matrix):
 635        self.set_generator(generator, matrix, compute_inverse=True)
 636
 637    def change_base_ring(self, base_ring=None):
 638        return self._compose(
 639            lambda M: utils.change_base_ring(M, base_ring),
 640            base_ring=base_ring
 641        )
 642
 643    def astype(self, dtype):
 644        base_ring = self.base_ring
 645        if dtype != np.dtype('O'):
 646            base_ring = None
 647
 648        new_rep = self._compose(
 649            lambda M: M.astype(dtype),
 650            dtype=dtype
 651        )
 652        new_rep._base_ring = base_ring
 653
 654        return new_rep
 655
 656    def conjugate(self, mat, inv_mat=None, unwrap=True, **kwargs):
 657        if not unwrap:
 658            return self._conjugate(
 659                mat, inv_mat, **kwargs
 660            )
 661
 662        if inv_mat is not None:
 663            inv_mat = self.__class__.unwrap_func(inv_mat)
 664
 665        return self._conjugate(
 666            self.__class__.unwrap_func(mat),
 667            inv_mat, **kwargs
 668        )
 669
 670    def _conjugate(self, mat, inv_mat=None, **kwargs):
 671        if inv_mat is None:
 672            inv_mat = utils.invert(mat, **kwargs)
 673
 674        return self._compose(lambda M: inv_mat @ M @ mat)
 675
 676    def _compose(self, hom, compute_inverses=False, dtype=None,
 677                 base_ring=None, **kwargs):
 678        """Get a new representation obtained by composing this representation
 679        with hom."""
 680
 681        if dtype is None:
 682            dtype = self.dtype
 683
 684        if base_ring is None:
 685            base_ring = self.base_ring
 686
 687        composed_rep = self.__class__(dtype=dtype, base_ring=base_ring,
 688                                      relations=self.relations, **kwargs)
 689
 690        if compute_inverses:
 691            generator_iterator = self.asym_gens()
 692        else:
 693            generator_iterator = self.generators.keys()
 694
 695        for g in generator_iterator:
 696            image = self.generators[g]
 697            inv_image = self.generators[words.invert_gen(g)]
 698
 699            try:
 700                composed = hom(image, inv=inv_image)
 701            except TypeError:
 702                composed = hom(image)
 703
 704            composed_rep._set_generator(g, composed,
 705                                        compute_inverse=compute_inverses)
 706
 707        return composed_rep
 708
 709    def compose(self, hom, hom_in_wrapped=False, hom_out_wrapped=False,
 710                **kwargs):
 711
 712        hom_1 = hom
 713        if hom_in_wrapped:
 714            hom_1 = (lambda M: hom(self.__class__.unwrap_func(M)))
 715
 716        hom_2 = hom_1
 717        if hom_out_wrapped:
 718            hom_2 = (lambda M: self.__class__.unwrap_func(hom_1(M)))
 719
 720        composed_rep = self._compose(hom_2, **kwargs)
 721
 722        return composed_rep
 723
 724    def _differential(self, word, generator=None, verbose=False):
 725        if generator is not None:
 726            if verbose:
 727                print(
 728                    "computing differential of '{}' with respect to {}".format(
 729                        word, generator)
 730                )
 731            word_diff = words.fox_word_derivative(generator, word)
 732            matrix_diff = [
 733                coeff * self._word_value(word)
 734                for word, coeff in word_diff.items()
 735            ]
 736            if len(matrix_diff) == 0:
 737                return utils.zeros(self.dim, base_ring=self.base_ring)
 738
 739            return np.sum(matrix_diff, axis=0)
 740
 741        blocks = [self._differential(word, g, verbose=verbose)
 742                  for g in self.asym_gens()]
 743
 744        return np.concatenate(blocks, axis=-1)
 745
 746    def subgroup(self, generators, generator_names=None, relations=[],
 747                 compute_inverse=True):
 748        subrep = self.__class__(relations=relations,
 749                                base_ring=self.base_ring,
 750                                dtype=self.dtype)
 751
 752        try:
 753            generator_pairs = generators.items()
 754        except AttributeError:
 755            generator_pairs = zip(GENERATORS[:len(generators)],
 756                                  generators)
 757
 758        if generator_names is not None:
 759            generator_pairs = zip(generator_names, generators)
 760
 761        for g, word in generator_pairs:
 762            subrep._set_generator(g, self._word_value(word),
 763                                  compute_inverse=compute_inverse)
 764            if not compute_inverse:
 765                subrep._set_generator(
 766                    words.invert_gen(g),
 767                    self._word_value(words.formal_inverse(word))
 768                )
 769
 770        return subrep
 771
 772    def _differentials(self, words, **kwargs):
 773        blocks = [self._differential(word, **kwargs) for word in words]
 774        return np.concatenate(blocks, axis=0)
 775
 776    def differential(self, word, **kwargs):
 777        return self._differential(word, **kwargs)
 778
 779    def differentials(self, words, **kwargs):
 780        return self._differentials(words, **kwargs)
 781
 782    def cocycle_matrix(self, **kwargs):
 783        return self.differentials(self.relations, **kwargs)
 784
 785    def coboundary_matrix(self):
 786        blocks = [utils.identity(self._dim, dtype=self.dtype,
 787                                 base_ring=self.base_ring) - self.generators[gen]
 788                  for gen in self.asym_gens() ]
 789
 790        return np.concatenate(blocks, axis=0)
 791
 792    def gln_adjoint(self, base_ring=None, dtype=None, **kwargs):
 793        if base_ring is None:
 794            base_ring = self.base_ring
 795
 796        if dtype is None:
 797            dtype = self.dtype
 798
 799        gln_adjoint = lie.hom.gln_adjoint(
 800            base_ring=base_ring, dtype=dtype
 801        )
 802
 803        return self._compose(gln_adjoint, base_ring=base_ring,
 804                             dtype=dtype, **kwargs)
 805
 806    def sln_adjoint(self, base_ring=None, dtype=None, **kwargs):
 807        if base_ring is None:
 808            base_ring = self.base_ring
 809
 810        if dtype is None:
 811            dtype = self.dtype
 812
 813        sln_adjoint = lie.hom.sln_adjoint(
 814            base_ring=base_ring, dtype=dtype
 815        )
 816
 817        return self._compose(sln_adjoint, base_ring=base_ring,
 818                             dtype=dtype, **kwargs)
 819
 820    def tensor_product(self, rep):
 821        """Return a tensor product of this representation with `rep`.
 822
 823        Parameters
 824        ----------
 825        rep : Representation
 826            Representation to tensor with.
 827
 828        Raises
 829        ------
 830        ValueError
 831            Raised if `self` and `rep` have differing generating sets.
 832
 833        Returns
 834        -------
 835        tensor: Representation
 836            Representation giving the tensor product of self with `rep`.
 837
 838        """
 839        if set(rep.generators) != set(self.generators):
 840            raise ValueError(
 841                "Cannot take a tensor product of a representation of groups with "
 842                "different presentations"
 843            )
 844        else:
 845            product_rep = Representation()
 846            for gen in self.asym_gens():
 847                tens = np.tensordot(self[gen], rep[gen], axes=0)
 848                elt = np.concatenate(np.concatenate(tens, axis=1), axis=1)
 849                product_rep[gen] = np.array(elt)
 850            return product_rep
 851
 852    def symmetric_square(self):
 853        """Return the symmetric square of this representation.
 854
 855        Returns
 856        -------
 857        square : Representation
 858            Symmetric square of `self`.
 859
 860        """
 861
 862        tensor_rep = self.tensor_product(self)
 863        incl = symmetric_inclusion(self._dim)
 864        proj = symmetric_projection(self._dim)
 865        square_rep = Representation()
 866        for g in self.asym_gens():
 867            square_rep[g] = proj * tensor_rep[g] * incl
 868
 869        return square_rep
 870
 871class SageMatrixRepresentation(Representation):
 872    @staticmethod
 873    def wrap_func(numpy_matrix):
 874        return sagewrap.sage_matrix(numpy_matrix)
 875
 876    @staticmethod
 877    def unwrap_func(sage_matrix):
 878        return np.array(sage_matrix)
 879
 880    @staticmethod
 881    def array_wrap_func(numpy_array):
 882        return sagewrap.sage_matrix_list(numpy_array)
 883
 884    def differential(self, word, **kwargs):
 885        return sagewrap.sage_matrix(self._differential(word, **kwargs))
 886
 887    def differentials(self, words, **kwargs):
 888        return sagewrap.sage_matrix(self._differentials(words, **kwargs))
 889
 890    def coboundary_matrix(self):
 891        return sagewrap.sage_matrix(
 892            Representation.coboundary_matrix(self)
 893        )
 894
 895    def compose(self, hom, hom_in_wrapped=True,
 896                hom_out_wrapped=False, **kwargs):
 897        return Representation.compose(self, hom,
 898                                      hom_in_wrapped=hom_in_wrapped,
 899                                      hom_out_wrapped=hom_out_wrapped,
 900                                      **kwargs)
 901
 902def sym_index(i, j, n):
 903    r"""Return coordinate indices for an isomorphism
 904    $\mathrm{Sym}^2(\mathbb{R}^n) \to \mathbb{R}^{\binom{n}{2} + n}$.
 905
 906    If $\{e_1, \ldots, e_n\}$ is the standard basis for $\mathbb{R}^n$,
 907    the isomorphism is realized by giving $\mathrm{Sym}^2(\mathbb{R}^n)$
 908    the ordered basis
 909    $$
 910        \{e_ne_n, e_{n-1}e_{n-1}, e_{n-1}e_n,
 911        e_{n-1}e_{n-1}, e_{n-1}e_{n-2}, e_{n-1}e_{n-3}, \ldots \}.
 912    $$
 913    Schematically this is given by the symmetric matrix
 914
 915    $$\begin{pmatrix} \ddots \\
 916    & 3 & 4 & 5 \\
 917    & & 1 & 2 \\
 918    & & & 0
 919    \end{pmatrix},
 920    $$
 921    where the (i,j) entry of the matrix gives the index of basis
 922    element $e_ie_j$.
 923
 924    Parameters
 925    ----------
 926    i : int
 927        index of one of the terms in the basis monomial $e_ie_j$ for the
 928        symmetric square
 929    j : int
 930        index of the other term in the basis monomial $e_ie_j$ for the
 931        symmetric square
 932    n : int
 933        dimension of the underlying vector space $\mathbb{R}^n$.
 934
 935    Returns
 936    -------
 937    int
 938        index of the corresponding basis vector in
 939        $\mathbb{R}^{\binom{n}{2} + n}$.
 940
 941    """
 942    if i > j:
 943        i, j = j, i
 944    return int((n - i) * (n - i  - 1) / 2 + (j - i))
 945
 946def tensor_pos(i, n):
 947    r"""Return coordinate indices for an isomorphism
 948    $\mathbb{R}^{n^2} \to \mathbb{R}^n \otimes \mathbb{R}^n$.
 949
 950    If $\{e_1, \ldots, e_n\}$ is the standard basis for
 951    $\mathbb{R}^n$, the isomorphism is realized by giving
 952    $\mathbb{R}^n \otimes \mathbb{R}^n$ the ordered basis
 953    $$
 954    \{e_1 \otimes e_1, e_1 \otimes e_2, \ldots, e_1 \otimes e_n, e_2 \otimes e_1, \ldots, \}
 955    $$
 956    represented schematically by the matrix
 957    $$
 958        \begin{pmatrix}
 959            0 & 1 & \ldots \\
 960            n & n + 1 & \ldots\\
 961            \vdots
 962        \end{pmatrix}.
 963    $$
 964    Here the (i, j) entry of the matrix gives the index of the basis
 965    element $e_i \otimes e_j$.
 966
 967    The inverse of this isomorphism is given by `tensor_index`.
 968
 969    Parameters
 970    ----------
 971    i : int
 972        index of a basis vector in $\mathbb{R}^{n^2}$
 973    n : int
 974        dimension of the underlying vector space $\mathbb{R}^n$
 975
 976    Returns
 977    -------
 978    tuple
 979        tuple `(j, k)` determining the monomial $e_j \otimes e_k$
 980        mapped to by given the basis vector in $\mathbb{R}^{n^2}$.
 981
 982    """
 983    return int(i / n), i % n
 984
 985def tensor_index(i,j,n):
 986    r"""Return coordinate indices for an isomorphism
 987    $\mathbb{R}^n \otimes \mathbb{R}^n \to \mathbb{R}^{n^2}$.
 988
 989    If $\{e_1, \ldots, e_n\}$ is the standard basis for
 990    $\mathbb{R}^n$, the isomorphism is realized by giving
 991    $\mathbb{R}^n \otimes \mathbb{R}^n$ the ordered basis
 992    $$
 993    \{e_1 \otimes e_1, e_1 \otimes e_2, \ldots, e_1 \otimes e_n, e_2 \otimes e_1, \ldots, \}
 994    $$
 995    represented schematically by the matrix
 996    $$
 997        \begin{pmatrix}
 998            0 & 1 & \ldots \\
 999            n & n + 1 & \ldots\\
1000            \vdots
1001        \end{pmatrix}.
1002    $$
1003    Here the (i, j) entry of the matrix gives the index of the basis
1004    element $e_i \otimes e_j$.
1005
1006    The inverse of this isomorphism is given by `tensor_pos`.
1007
1008    Parameters
1009    ----------
1010    i : int
1011        index of one of the terms in a basis vector $e_i \otimes e_j$.
1012    j : int
1013        index of the other term in a basis vector $e_i \times e_j$.
1014    n : int
1015        dimension of the underlying vector space $\mathbb{R}^n$
1016
1017    Returns
1018    -------
1019    int
1020        index of a basis vector in $\mathbb{R}^{n^2}$ mapped to by
1021        $e_i \otimes e_j$.
1022    """
1023    return i * n + j
1024
1025def symmetric_inclusion(n):
1026    r"""Return a matrix representing the linear inclusion
1027    $\mathrm{Sym}^2(\mathbb{R}^n) \to \mathbb{R}^n \otimes
1028    \mathbb{R}^n$.
1029
1030    $\mathrm{Sym}^2(\mathbb{R}^n)$ and
1031    $\mathbb{R}^n \otimes \mathbb{R}^n$
1032    are respectively identified with
1033    $\mathbb{R}^{\binom{n}{2} + n}$ and $\mathbb{R}^{n^2}$ via the
1034    isomorphisms described in `sym_index`, `tensor_index`, and
1035    `tensor_pos`.
1036
1037    If $\{e_1, \ldots, e_n\}$ is the standard basis for
1038    $\mathbb{R}^n$, the returned matrix gives the linear map taking
1039    $e_ie_j$ to $\frac{1}{2}(e_i \otimes e_j + e_j \otimes e_i)$,
1040    with respect to the bases specified above.
1041
1042    Parameters
1043    ----------
1044    n : int
1045        Dimension of the underlying vector space $\mathbb{R}^n$.
1046
1047    Returns
1048    -------
1049    matrix : ndarray
1050        $n^2 \times \binom{n}{2} + n$ array defining this linear map.
1051
1052    """
1053    incl_matrix = np.zeros((n * n, int(n * (n + 1) / 2)))
1054    for i in range(n):
1055        for j in range(n):
1056            si = sym_index(i, j, n)
1057            ti = tensor_index(i, j, n)
1058            incl_matrix[ti][si] = 1/2 + (i == j) * 1/2
1059
1060    return np.array(incl_matrix)
1061
1062def symmetric_projection(n):
1063    r"""Return a matrix representing the linear surjection
1064    $\mathbb{R}^n \otimes \mathbb{R}^n \to \mathrm{Sym}^2(\mathbb{R}^n)$.
1065
1066    If $\mathbb{R}^n$ is given the standard basis $\{e_1, \ldots,
1067    e_n\}$, then this matrix represents the linear map determined by
1068    $e_i \otimes e_j \mapsto e_ie_j$. The spaces
1069    $\mathbb{R}^n \otimes \mathbb{R}^n$ and $\mathrm{Sym}^2(\mathbb{R}^n)$
1070    are given the ordered bases determined by the functions
1071    `sym_index`, `tensor_index`, and `tensor_pos`.
1072
1073    Parameters
1074    ----------
1075    n : int
1076        Dimension of the underlying vector space $\mathbb{R}^n$
1077
1078    Returns
1079    -------
1080    ndarray
1081        $\binom{n}{2} + n \times n$ matrix representing the linear map
1082        in the given bases.
1083
1084    """
1085    proj_matrix = np.zeros((int(n * (n + 1) / 2), n * n))
1086    for i in range(n * n):
1087        u, v = tensor_pos(i,n)
1088        proj_matrix[_sym_index(u, v, n)][i] = 1
1089
1090    return np.array(proj_matrix)