Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[WIP] Added ArraySlice class so that ArraySymbol accepts slices #22265

Closed
wants to merge 45 commits into from
Closed

[WIP] Added ArraySlice class so that ArraySymbol accepts slices #22265

wants to merge 45 commits into from

Conversation

redeboer
Copy link
Contributor

@redeboer redeboer commented Oct 13, 2021

References to other Issues or PRs

Brief description of what is fixed or changed

  • Added an ArraySlice class that allows taking slices from an ArraySymbol:

    >>> from sympy import ArraySymbol
    >>> A = ArraySymbol("A")
    >>> A[:, 3:-5:2, 0]
    A[:, 3:-5:2, 0]
  • It's now also possible to take a single index from an array:

    >>> A[0]
    A[0]
  • Negative indices are also allowed (normalized in the same way as MatrixSymbol does):

    >>> B = ArraySymbol("B", 5, 3)
    >>> B[-2, -1]
    B[3, 2]
  • Other changes are related to typing and making array expressions more like matrix expressions. See commit history.

Other comments

Motivation is to have ArraySymbol behave more like numpy.array or tensorflow.Tensor, so that you can use SymPy expressions as a template for these different computational back-ends (through lambdify).

Release Notes

  • tensor
    • Defined the ArraySymbol class as an equivalent to MatrixSlice, so that it's possible to take slices of an ArraySymbol.
    • ArraySymbol now accepts a single index, as well as negative indices.

Pins down existing behaviour so we can safely refactor it
This helps static code analysers like mypy and code navigators like
pylance
This gives a better overview of the way in which self._args should be
interpreted. Several other modules also adopted this style.
All tests that involve ArrayElement use an ArraySymbol for the "name"
argument. For the implementation, it's also safer to have a reference
to the parent ArraySymbol, like when comparing ArrayElements.

Note that MatrixElement also treats the name argument as a parent
(although it doesn't perform a type check):
https://github.com/sympy/sympy/blob/6de9857/sympy/matrices/expressions/matexpr.py#L719
Same as with Python list and numpy.array, something like the following
should be allowed:

>>> some_list = [0, 1, 2]
>>> some_list[:5]
[0, 1, 2]
@sympy-bot
Copy link

sympy-bot commented Oct 13, 2021

Hi, I am the SymPy bot (v162). I'm here to help you write a release notes entry. Please read the guide on how to write release notes.

Your release notes are in good order.

Here is what the release notes will look like:

  • tensor
    • Defined the ArraySymbol class as an equivalent to MatrixSlice, so that it's possible to take slices of an ArraySymbol. (#22265 by @ComPWA and @redeboer)

    • ArraySymbol now accepts a single index, as well as negative indices. (#22265 by @ComPWA and @redeboer)

This will be added to https://github.com/sympy/sympy/wiki/Release-Notes-for-1.10.

Click here to see the pull request description that was parsed.
#### References to other Issues or PRs

- Related to #20570
- Follow-up to #20943

#### Brief description of what is fixed or changed

- Added an `ArraySlice` class that allows taking slices from an `ArraySymbol`:
  ```python
  >>> from sympy import ArraySymbol
  >>> A = ArraySymbol("A")
  >>> A[:, 3:-5:2, 0]
  A[:, 3:-5:2, 0]
  ```

- It's now also possible to take a single index from an array:
  ```python
  >>> A[0]
  A[0]
  ```

- Negative indices are also allowed (normalized in the same way as `MatrixSymbol` does):
  ```python
  >>> B = ArraySymbol("B", 5, 3)
  >>> B[-2, -1]
  B[3, 2]
  ```

- Other changes are related to typing and making array expressions more like matrix expressions. See [commit history](https://github.com/sympy/sympy/pull/22265/commits).

#### Other comments

Motivation is to have `ArraySymbol` behave more like [`numpy.array`](https://numpy.org/doc/stable/reference/generated/numpy.array.html) or [`tensorflow.Tensor`](https://www.tensorflow.org/api_docs/python/tf/Tensor), so that you can use SymPy expressions as a template for these different computational back-ends (through `lambdify`).

#### Release Notes

<!-- BEGIN RELEASE NOTES -->
* tensor
  * Defined the ArraySymbol class as an equivalent to MatrixSlice, so that it's possible to take slices of an ArraySymbol.
  * ArraySymbol now accepts a single index, as well as negative indices.
<!-- END RELEASE NOTES -->

@sympy-bot
Copy link

sympy-bot commented Oct 13, 2021

🟠

Hi, I am the SymPy bot (v162). I've noticed that some of your commits add or delete files. Since this is sometimes done unintentionally, I wanted to alert you about it.

This is an experimental feature of SymPy Bot. If you have any feedback on it, please comment at sympy/sympy-bot#75.

The following commits add new files:

  • 837d238:
    • sympy/tensor/array/expressions/tests/test_array_symbol.py

The following commits delete files:

  • 7878f19:
    • sympy/tensor/array/expressions/tests/test_array_symbol.py

If these files were added/deleted on purpose, you can ignore this message.

@redeboer
Copy link
Contributor Author

redeboer commented Oct 13, 2021

@Upabjojr hope this PR is in line with your ongoing work on array expressions? For instance, not sure about potential code duplication with matrix expressions and I'm also not sure how to correctly work with kind (8c0d65f).

@redeboer
Copy link
Contributor Author

redeboer commented Oct 13, 2021

Is there a 'sympy way' to parametrize tests? Apparently can't do this through pytest:
https://github.com/sympy/sympy/runs/3883692406

@ThePauliPrinciple
Copy link
Contributor

I am interested in writing the code printers for array like objects. In particular, I am currently working on trying to print Indexed objects with numpy and tensorflow. While some pattern matching is needed for Indexed compared to Arrays (e.g. x_i-x_j would need x[:,newaxis]-x[newaxis,:] as the array expression), I'll gladly write the code printers for the Array objects and have the printer for Indexed objects use the printing routines for Arrays when appropriate to avoid code duplication.

Please also see my ideas in: #22219

For this I would however, like to see a newaxis object and if possible the option to select specific indices through an array of integers (numpy supports this as far as I understand and tensorflow allows this through tf.gather and tf.scatter). Would you consider adding those options to this merge request?

I also have a question about this merge request (for my future reference). Is it prefered to have many small commits in a merge request or only a single large commit? I could argue it either way.

Also, is it appropriate to ask for additional features in a merge request? Or should I have written it as a separate issue?

@redeboer
Copy link
Contributor Author

#22265 (comment)

Sounds good, yes! While working on this PR, I was already trying to implement some more arithmetic (only addition etc), but I got the impression that ArraySymbol and its related classes have its limitations. For instance, I couldn't get to write something that could be lambdified to numpy.sum(array, axis=1), which sounds like something your ideas might address.

I also have a question about this merge request (for my future reference). Is it prefered to have many small commits in a merge request or only a single large commit? I could argue it either way.

Not sure either; I was trying some sort of TDD here, and perhaps it helps when reviewing the PR. In the end, the commits can be squashed. (Personally, I'm in favour of squash-merging, so that the commit history within the PR doesn't matter anymore).

Also, is it appropriate to ask for additional features in a merge request? Or should I have written it as a separate issue?

Idk, not a sympy developer ;) It sounds like a bigger idea, so indeed worth to post an issue on it (perhaps just linking to the discussion suffices).

@ThePauliPrinciple
Copy link
Contributor

In the case of trying to sum over an axis, what would the sympy expression be?

@redeboer
Copy link
Contributor Author

redeboer commented Oct 13, 2021

In the case of trying to sum over an axis, what would the sympy expression be?

Beats me 🤷‍♂️ Afaics there are no sympy classes that lambdify to numpy.sum.

Perhaps something can be defined through PermuteDims, ArrayContraction, ArrayTensorProduct, etc., but I didn't get that to work. They lambdify to numpy.transpose and einsum though:

def _print_ArrayTensorProduct(self, expr):
array_list = [j for i, arg in enumerate(expr.args) for j in
(self._print(arg), "[%i, %i]" % (2*i, 2*i+1))]
return "%s(%s)" % (self._module_format(self._module + '.einsum'), ", ".join(array_list))
def _print_ArrayContraction(self, expr):
from ..tensor.array.expressions.array_expressions import ArrayTensorProduct
base = expr.expr
contraction_indices = expr.contraction_indices
if not contraction_indices:
return self._print(base)
if isinstance(base, ArrayTensorProduct):
counter = 0
d = {j: min(i) for i in contraction_indices for j in i}
indices = []
for rank_arg in base.subranks:
lindices = []
for i in range(rank_arg):
if counter in d:
lindices.append(d[counter])
else:
lindices.append(counter)
counter += 1
indices.append(lindices)
elems = ["%s, %s" % (self._print(arg), ind) for arg, ind in zip(base.args, indices)]
return "%s(%s)" % (
self._module_format(self._module + '.einsum'),
", ".join(elems)
)
raise NotImplementedError()
def _print_ArrayDiagonal(self, expr):
diagonal_indices = list(expr.diagonal_indices)
if len(diagonal_indices) > 1:
# TODO: this should be handled in sympy.codegen.array_utils,
# possibly by creating the possibility of unfolding the
# ArrayDiagonal object into nested ones. Same reasoning for
# the array contraction.
raise NotImplementedError
if len(diagonal_indices[0]) != 2:
raise NotImplementedError
return "%s(%s, 0, axis1=%s, axis2=%s)" % (
self._module_format("numpy.diagonal"),
self._print(expr.expr),
diagonal_indices[0][0],
diagonal_indices[0][1],
)
def _print_PermuteDims(self, expr):
return "%s(%s, %s)" % (
self._module_format("numpy.transpose"),
self._print(expr.expr),
self._print(expr.permutation.array_form),
)

@ThePauliPrinciple
Copy link
Contributor

ThePauliPrinciple commented Oct 13, 2021

I think einsum can be correctly used for summing over the axis of a single Array. E.g.

import numpy as np
A=np.ones((3,))
np.einsum('i->',A)

Would yield a scalar=3. I would imagine that it would shortcut to a special case and just run np.sum if it detects this is possible.

One could also imagine something like i=Array.get_axis(1) and then Sum(Array,i) could lead to printing np.sum(Array,axis=1).

If you need any help with the printing routines, or want me to give them an attempt that is no problem. I'm not very experienced with the Array part of sympy though, so I'm looking for your help on that side e.g giving an example of what sympy expression should lead to what print statement, or even simply the sympy expression (if it is clear what it should mean) that should be supported for printing.

@redeboer
Copy link
Contributor Author

If you need any help with the printing routines, or want me to give them an attempt that is no problem. I'm not very experienced with the Array part of sympy though, so I'm looking for your help on that side e.g giving an example of what sympy expression should lead to what print statement, or even simply the sympy expression (if it is clear what it should mean) that should be supported for printing.

Thanks, let's move this discussion to #22269

@ThePauliPrinciple
Copy link
Contributor

ThePauliPrinciple commented Oct 13, 2021

Would you consider ArraySlice to allow None (both numpy and tensorflow allow this as an alternative to newaxis), ArraySymbols and Arrays as valid slices? to replicate this behaviour:

afbeelding

Considering #22269 (reply in thread) allowing symbols as axes would be nice too

@redeboer
Copy link
Contributor Author

So guys, is there anything about the PR that needs to be changed? Sounds like the other points raised here are larger issues.

@Upabjojr
Copy link
Contributor

So guys, is there anything about the PR that needs to be changed? Sounds like the other points raised here are larger issues.

This PR makes the shape of ArraySymbol optional, right? Has there been an in-depth discussion about potential consequences of this change?

@redeboer
Copy link
Contributor Author

redeboer commented Oct 18, 2021

I added a small unit test (837d238) that pinned down existing behaviour. At the time, an ArraySymbol without shape was apparently allowed (see below). So unless there have been policy changes from the master in the meantime, this PR did not introduce this optional behaviour.

assert ArraySymbol("A") == ArraySymbol("A")

@pytest.mark.parametrize(
"shape",
[
(),
(3, 2, 4),
symbols("k m n"),
],
)
def test_constructor(self, shape: tuple):
A = ArraySymbol("A", *shape)
assert A.name == Symbol("A")
assert A.shape == shape


Update: following behaviour is indeed defined in the latest master (ef37f78)

>>> from sympy.tensor.array.expressions.array_expressions import ArraySymbol
>>> A = ArraySymbol("A")
>>> A[0, 0]
A[0, 0]

@Upabjojr
Copy link
Contributor

At the time, an ArraySymbol without shape was apparently allowed (see below).

Probably this is a bug. If you don't pass a shape to ArraySymbol, then *shape will be an empty tuple. Empty tuples are the shape of scalars, so it should not be possible to access elements of the ArraySymbol.

@redeboer
Copy link
Contributor Author

Probably this is a bug. If you don't pass a shape to ArraySymbol, then *shape will be an empty tuple. Empty tuples are the shape of scalars, so it should not be possible to access elements of the ArraySymbol.

Then a test should be written to prove it. I think this should be addressed elsewhere, not in this PR.

@redeboer
Copy link
Contributor Author

I think generally this points in the direction of a larger issue: too much logics going on in the constructors/new methods. For instance, ArrayContraction.__new__ should return an ArrayContraction (this is also mypy's opinion once one starts adding type hints). The logics for other types of input can then be implement in some function like array_contraction. But then again, this PR is not the place to fix that.

@Upabjojr
Copy link
Contributor

I think generally this points in the direction of a larger issue: too much logics going on in the constructors/new methods. For instance, ArrayContraction.__new__ should return an ArrayContraction (this is also mypy's opinion once one starts adding type hints). The logics for other types of input can then be implement in some function like array_contraction. But then again, this PR is not the place to fix that.

Yes, it's related to the problem of automatic normalization. Actually, all of SymPy is affected by this problem. If you call Add(x, x) you get Mul(2, x).

I would put the logic in a .normalize( ) method.

There is already a function meant to perform contractions, it's tensorcontraction (name is a bit confusing, because the word tensor is used instead of array, but it's logically equivalent to ArrayContraction).

@redeboer
Copy link
Contributor Author

Then a test should be written to prove it. I think this should be addressed elsewhere, not in this PR.

I do admit though that I indulged in some other clean-ups. So what I could do is clean up this PR (rebase) and extract some of those commits that were not immediately related to ArraySlice to another PR. Then another PR can be made as well to enforce the use of a shape (and, possibly, the option to have no shape, which I would prefer to have).

@Upabjojr
Copy link
Contributor

Then another PR can be made as well to enforce the use of a shape (and, possibly, the option to have no shape, which I would prefer to have).

Remember that the array expressions module has not been tested that much with undefined shapes. It could require some work to make array symbols with undefined shape work with all expression-tree objects.

@redeboer
Copy link
Contributor Author

Okay so sounds like those issues have to be addressed first. I could write some tests and so to enforce the use of a non-empty shape and PR that separately first. I'll also clean up this PR and have it focus only on defining an ArraySlice class.

@redeboer
Copy link
Contributor Author

Remember that the array expressions module has not been tested that much with undefined shapes. It could require some work to make array symbols with undefined shape work with all expression-tree objects.

Side note: I mainly have an interest in code generation, so I can't help much with the math :S

@redeboer redeboer marked this pull request as draft October 18, 2021 17:03
@asmeurer
Copy link
Member

Probably this is a bug. If you don't pass a shape to ArraySymbol, then *shape will be an empty tuple. Empty tuples are the shape of scalars, so it should not be possible to access elements of the ArraySymbol.

IMO using ArraySymbol(name, *shape) instead of ArraySymbol(name, shape) was a mistake. While it's true that *shape is consistent with MatrixSymbol, it's easier to think of shape as a distinct tuple when working with arrays.

@asmeurer
Copy link
Member

shape as a distinct tuple would also open the possibility of undefined shapes in the future. With *shape you always have to have the dimensionality be explicit.

@Upabjojr
Copy link
Contributor

IMO using ArraySymbol(name, *shape) instead of ArraySymbol(name, shape) was a mistake. While it's true that *shape is consistent with MatrixSymbol, it's easier to think of shape as a distinct tuple when working with arrays.

It can still be fixed. sympy.tensor.array.expressions is still a hidden module, no documentation is available and no links to that module are provided. We can still change the API.

@redeboer
Copy link
Contributor Author

Many points raised here, so I tried to collect them in #22321.

@Upabjojr
Copy link
Contributor

#22324

@redeboer redeboer changed the title Added ArraySlice class so that ArraySymbol accepts slices [WIP] Added ArraySlice class so that ArraySymbol accepts slices Oct 21, 2021
@ComPWA ComPWA closed this by deleting the head repository Sep 28, 2023
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

7 participants