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.

Sign up for GitHub

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

Fix Outer for SparseArray #940

Merged
merged 25 commits into from
Nov 27, 2023
Merged

Fix Outer for SparseArray #940

Show file tree
Hide file tree
Changes from 20 commits
Commits
Show all changes
25 commits
Select commit Hold shift + click to select a range
c4b6752
Update tensors.py
Li-Xiang-Ideal Nov 21, 2023
45149f3
Merge branch 'Mathics3:master' into master
Li-Xiang-Ideal Nov 23, 2023
983049d
Update SYMBOLS_MANIFEST.txt
Li-Xiang-Ideal Nov 23, 2023
e171e0c
Fix Outer for SparseArray
Li-Xiang-Ideal Nov 24, 2023
e7b68a3
Update tensors.py
Li-Xiang-Ideal Nov 24, 2023
0789bd6
Update tensors.py
Li-Xiang-Ideal Nov 24, 2023
4159b69
Merge branch 'Mathics3:master' into master
Li-Xiang-Ideal Nov 25, 2023
58dbb8b
Update tensors.py
Li-Xiang-Ideal Nov 25, 2023
9583f38
Update tensors.py
Li-Xiang-Ideal Nov 25, 2023
f799576
Update isort-and-black-checks.yml
Li-Xiang-Ideal Nov 25, 2023
f09b340
Update tensors.py
Li-Xiang-Ideal Nov 25, 2023
6dddce8
Update isort-and-black-checks.yml
Li-Xiang-Ideal Nov 25, 2023
12917a1
Update osx.yml
Li-Xiang-Ideal Nov 25, 2023
bbfc74c
Update osx.yml
Li-Xiang-Ideal Nov 25, 2023
9be0ca7
Update osx.yml
Li-Xiang-Ideal Nov 25, 2023
bb133b1
Update osx.yml
Li-Xiang-Ideal Nov 25, 2023
60611b4
Merge pull request #1 from Li-Xiang-Ideal/OSX-test-patch-1
Li-Xiang-Ideal Nov 25, 2023
8a0057f
Update osx.yml
Li-Xiang-Ideal Nov 25, 2023
e235c04
Update tensors.py
Li-Xiang-Ideal Nov 26, 2023
3af6f9a
Merge branch 'Mathics3:master' into master
Li-Xiang-Ideal Nov 26, 2023
4db5ea0
fix typo
Li-Xiang-Ideal Nov 27, 2023
792d218
move eval from ``mathics.builtin.tensors`` to here
Li-Xiang-Ideal Nov 27, 2023
312eaff
Update clusters.py
Li-Xiang-Ideal Nov 27, 2023
088cbed
Update tensors.py
Li-Xiang-Ideal Nov 27, 2023
daf7d6d
Update clusters.py
Li-Xiang-Ideal Nov 27, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions SYMBOLS_MANIFEST.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -574,6 +574,7 @@ System`LessEqual
System`LetterCharacter
System`LetterNumber
System`LetterQ
System`LeviCivitaTensor
System`Level
System`LevelQ
System`LightBlue
Expand Down
119 changes: 109 additions & 10 deletions mathics/builtin/tensors.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,19 +18,32 @@
of any rank can be handled.
"""

import itertools

from sympy.combinatorics import Permutation
from sympy.utilities.iterables import permutations

from mathics.core.atoms import Integer, String
from mathics.core.atoms import Integer, Integer0, Integer1, String
from mathics.core.attributes import A_FLAT, A_ONE_IDENTITY, A_PROTECTED
from mathics.core.builtin import BinaryOperator, Builtin
from mathics.core.convert.python import from_python
from mathics.core.evaluation import Evaluation
from mathics.core.expression import Expression
from mathics.core.list import ListExpression
from mathics.core.symbols import Atom, Symbol, SymbolFalse, SymbolTrue
from mathics.core.systemsymbols import SymbolRule, SymbolSparseArray
from mathics.core.symbols import (
Atom,
Symbol,
SymbolFalse,
SymbolList,
SymbolTimes,
SymbolTrue,
)
from mathics.core.systemsymbols import (
SymbolAutomatic,
SymbolNormal,
SymbolRule,
SymbolSparseArray,
)
from mathics.eval.parts import get_part


Expand Down Expand Up @@ -278,10 +291,21 @@ class Outer(Builtin):
Outer product of two matrices:
>> Outer[Times, {{a, b}, {c, d}}, {{1, 2}, {3, 4}}]
= {{{{a, 2 a}, {3 a, 4 a}}, {{b, 2 b}, {3 b, 4 b}}}, {{{c, 2 c}, {3 c, 4 c}}, {{d, 2 d}, {3 d, 4 d}}}}

Outer product of two sparse arrays:
>> Outer[Times, SparseArray[{{1, 2} -> a, {2, 1} -> b}], SparseArray[{{1, 2} -> c, {2, 1} -> d}]]
= SparseArray[Automatic, {2, 2, 2, 2}, 0, {{1, 2, 1, 2} -> a c, {1, 2, 2, 1} -> a d, {2, 1, 1, 2} -> b c, {2, 1, 2, 1} -> b d}]

'Outer' of multiple lists:
>> Outer[f, {a, b}, {x, y, z}, {1, 2}]
= {{{f[a, x, 1], f[a, x, 2]}, {f[a, y, 1], f[a, y, 2]}, {f[a, z, 1], f[a, z, 2]}}, {{f[b, x, 1], f[b, x, 2]}, {f[b, y, 1], f[b, y, 2]}, {f[b, z, 1], f[b, z, 2]}}}

'Outer' treats input sparse arrays as lists if f=!=Times, or if the input is a mixture of sparse arrays and lists:
>> Outer[f, SparseArray[{{1, 2} -> a, {2, 1} -> b}], SparseArray[{{1, 2} -> c, {2, 1} -> d}]]
= {{{{f[0, 0], f[0, c]}, {f[0, d], f[0, 0]}}, {{f[a, 0], f[a, c]}, {f[a, d], f[a, 0]}}}, {{{f[b, 0], f[b, c]}, {f[b, d], f[b, 0]}}, {{f[0, 0], f[0, c]}, {f[0, d], f[0, 0]}}}}

>> Outer[Times, SparseArray[{{1, 2} -> a, {2, 1} -> b}], {c, d}]
= {{{0, 0}, {a c, a d}}, {{b c, b d}, {0, 0}}}

Arrays can be ragged:
>> Outer[Times, {{1, 2}}, {{a, b}, {c, d, e}}]
Expand All @@ -304,17 +328,37 @@ class Outer(Builtin):
def eval(self, f, lists, evaluation: Evaluation):
"Outer[f_, lists__]"

# If f=!=Times, or lists contain both SparseArray and List, then convert all SparseArrays to Lists
lists = lists.get_sequence()
head = None
for list in lists:
if isinstance(list, Atom):
sparse_to_list = f != SymbolTimes
contain_sparse = False
comtain_list = False
rocky marked this conversation as resolved.
Show resolved Hide resolved
for _list in lists:
if _list.head.sameQ(SymbolSparseArray):
contain_sparse = True
if _list.head.sameQ(SymbolList):
comtain_list = True
sparse_to_list = sparse_to_list or (contain_sparse and comtain_list)
if sparse_to_list:
break
if sparse_to_list:
new_lists = []
for _list in lists:
if isinstance(_list, Atom):
evaluation.message("Outer", "normal")
return
if sparse_to_list:
if _list.head.sameQ(SymbolSparseArray):
_list = Expression(SymbolNormal, _list).evaluate(evaluation)
new_lists.append(_list)
if head is None:
head = list.head
elif not list.head.sameQ(head):
evaluation.message("Outer", "heads", head, list.head)
head = _list.head
elif not _list.head.sameQ(head):
evaluation.message("Outer", "heads", head, _list.head)
return
if sparse_to_list:
lists = new_lists

def rec(item, rest_lists, current):
evaluation.check_stopped()
Expand All @@ -329,7 +373,62 @@ def rec(item, rest_lists, current):
elements.append(rec(element, rest_lists, current))
return Expression(head, *elements)

return rec(lists[0], lists[1:], [])
def rec_sparse(item, rest_lists, current):
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Maybe this code deserves to be moved to a mathics.eval module.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Could you tell me more about how the codes are organized in Mathics3? So that I can determine which code should be moved to mathics.eval or other locations in future updates.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This aspect is described in https://mathics-development-guide.readthedocs.io/en/latest/extending/code-overview/python-modules.html and mathics.eval specifically in https://mathics-development-guide.readthedocs.io/en/latest/extending/code-overview/python-modules.html

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for the references. But I still have some specific questions left to figure out

  1. Should mathics.eval be as general as possible, or is it just a place for eval() of functions in mathics.builtin?
  2. I don't feel any essential difference between rec() and rec_sparse(). So if rec_sparse() should be moved to mathics.eval.tensor, I'd prefer to treat rec() (in Outer and in Inner) the same way.
  3. If rec() is moved to mathics.eval.tensor, it will require as many as 6 parameters (item, rest_lists, current, f, head, evaluation), which seems a bit weird to me.

Copy link
Member

@rocky rocky Nov 26, 2023 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The question as is leads me to believe that I might not have have conveyed the motivation. So let me start with that and wind up with some guidelines which I hope you can then decide on how to organize.

Currently Mathics3 evaluates expressions by walking a tree (which is specifically a kind of M-expression)

If Mathics3 is to be sped up, we need to turn the tree evaluation into a sequence of instructions. The functions that start eval_ you can think of as the instructions. They are top-level functions like Add, Subtract, Factor, etc. but with parameter checking done and objects like Integers and Strings turned into possibly lower, level objects. Basically there would be at least one for every Wolfram Alpha builtin function.

Of course, to support these high-level built in functions, we need other functions. Those will not start with eval_.

With respect to rec() and rec_sparse() specifically, if you do not feel any essential difference between them, then, sure, they both can be moved somewhere under the the mathics.eval module. Please, can we we come up with a more descriptive name than rec() and rec_sparse() and add docstrings to these functions to as well?

So finally, I have to describe or apologize for the crappiness of this code, since we are trying to encourage people not to follow the bad habits of what came before.

Originally, this was one huge monolith and we have not come close to digging ourselves out of the mess that is there. In the hierarchical OO monolith, there were a number of nested private functions such as rec() (or is it "wreck") that had to duplicated in that organization and because they were or are, well, private functions which probably should not have been private but instead put in a common module where the function can be used both places.

Finally as for the number of parameters. It probably means that there is a structure, object, concept or organization that we are missing that exists that encapsulates those parts into a more comprehensive whole.

However since I don't understand what this rec thing does, clearly I can't offer a suggestion. If you have an idea here and want to go with that, please suggest something or make that change.

Otherwise for the first cut, I would not worry about this initially. We can come back to that. What I have found is that there is so much of a mess, we just can't fix everything all at once. And whatever the ultimate better solution is, having splitting off the common rec thing into a separate module will only help things.

Copy link
Member

@rocky rocky Nov 27, 2023 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I looked at the situation with the rec() and rec_sparse() in more detail. Some of the things I said were a bit misleading, or just wrong (because I did not understand the true situation).

So let me clarify the situation here and correct some things I wrote.

First, yes indeed, this code goes back to the earliest times and that means the coding style is not good; names are misleading and there are neither docstrings nor documentation describing the code.

I suppose "rec" relates to the recursive implementation. The only way one could make this more vague would to use f() or fn(). But thankfully, at least the wrong term has not be used which was sometimes done.

In good programming practice, function names reflect what a function does more than how it is implemented. And here the name is vague too.

As a result, you find a number of internal functions called rec() in different Builtin classes that are really not related to one another. But the two functions we are interested in here which are similar are rec() and rec_sparse() from Builtin Outer.

A better name for these two functions would be create_outer_piece() and create_outer_piece_sparse(). If it is useful to mention these work recursively, that would be done in the docstring or description of the function. For example:

"""Recursively computes the outer product of list `item`  with `rest_lists`; 
`current` contains the results computed so far.
"""

It might also be nice to add type annotations. (This code was written in Python 2.7 when type annotations did not exist. But there comments, docstrings and good programming practice did exist in Python 2.7.)

Here is a suggestion for reorganization...

Create a function called eval_Outer() in mathics.eval.tensors , and move the Python code in current method in Outer called eval() to that renaming and improving the rec() and rec_sparse() methods. They can still be an internal function as it is here.

eval() then just calls eval_Outer .

Later, if we want reorganized things further we are in a better position to do that.

I suspect that you like most people were expecting you'd just make this little change and to get the code to understand sparse arrays and be done with it. I know if I were contributing here, that's what I expect.

The reality is there is a big mess that we picked up that we are still trying to dig ourselves out of. In the past, a great deal of our effort has been just clean up work of the kind mentioned here.

It has been a big mistake to just keep coding in the style and organization that came before. That's what led to this being a big mess rather than just a mess.

To the extent you can help us dig out of this, this is greatly appreciated.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I have to say, my main concern was about having a more or less complicated implementation of the eval_ methods inside the mathics.builtin.tensor module, not about the private, nested functions. Of course, it is better if these functions are properly named and documented, but at this stage, it would be enough for me to move the implementation to mathics.eval.tensor, something like


def eval_outer(f, list, evaluation):
    """Evaluates recursively the outer product"""
    ...

and

In mathics.builtin.tensor


    def eval(self, f, lists, evaluation):
        "Outerf[f_, lists__]"
        return eval_outer(f, lists, evaluation)

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sounds like the C/C++ style header-source separation?

Copy link
Member

@rocky rocky Nov 27, 2023 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, I guess so. I suppose this would be a good way to describe it in the developer documentation.

And like C/C++ header-source separation it is not strictly mandatory to make code work, just accepted practice.

As @mmatera says, this is good. We can merge this as is, or do you want to do the other stuff in addition. That can be done in another PR if you prefer.

Let us know.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I believe that (at least in tensor-related evaluations) there are only two cases where recursion is needed, namely tensor contraction (i.e. Inner, TensorContract, etc.) and direct product (i.e. Outer, TensorProduct, etc.). So when I have time I'll try to combine rec() and rec_sparse() into a single public function called create_outer_piece() or something else, and do the same thing for Inner (obviously Inner shares the same bug as Outer).

However, that is a completely different thing than fixing Outer. So I would prefer this PR to be merged as is, then I will do the above things in another PR.

evaluation.check_stopped()
if isinstance(item, tuple): # (rules)
elements = []
for element in item:
elements.extend(rec_sparse(element, rest_lists, current))
return tuple(elements)
else: # rule
_pos, _val = item.elements
if rest_lists:
return rec_sparse(
rest_lists[0],
rest_lists[1:],
(current[0] + _pos.elements, current[1] * _val),
)
else:
return (
Expression(
SymbolRule,
ListExpression(*(current[0] + _pos.elements)),
current[1] * _val,
),
)

# head != SparseArray
if not head.sameQ(SymbolSparseArray):
return rec(lists[0], lists[1:], [])

# head == SparseArray
dims = []
val = Integer1
data = [] # data = [(rules), ...]
for _list in lists:
_dims, _val, _rules = _list.elements[1:]
dims.extend(_dims)
val *= _val
if _val == Integer0: # _val==0, append (_rules)
data.append(_rules.elements)
else: # _val!=0, append (_rules, other pos->_val)
other_pos = []
for pos in itertools.product(*(range(1, d.value + 1) for d in _dims)):
other_pos.append(ListExpression(*(Integer(i) for i in pos)))
rules_pos = set(rule.elements[0] for rule in _rules.elements)
other_pos = set(other_pos) - rules_pos
other_rules = []
for pos in other_pos:
other_rules.append(Expression(SymbolRule, pos, _val))
data.append(_list.elements[3].elements + tuple(other_rules))
dims = ListExpression(*dims)
return Expression(
SymbolSparseArray,
SymbolAutomatic,
dims,
val,
ListExpression(*rec_sparse(data[0], data[1:], ((), Integer1))),
)


class RotationTransform(Builtin):
Expand Down Expand Up @@ -550,7 +649,7 @@ def eval(self, d, type, evaluation: Evaluation):

if isinstance(d, Integer) and type == SymbolSparseArray:
d = d.get_int_value()
perms = list(permutations([i for i in range(1, d + 1)]))
perms = list(permutations(list(range(1, d + 1))))
rocky marked this conversation as resolved.
Show resolved Hide resolved
rules = [
Expression(
SymbolRule,
Expand Down
Loading