Arithmetic on tensor module elements, manifold objects: Always Return a Copy

[mjungmath]

This question arose from ticket #30239, comment:36.

Should FiniteRankFreeModule and manifold objects always return a mutable copy, even for trivial operations? At least, this would be a consistent behavior.

As pointed out by Matthias, this already holds true for FreeModule:

sage: M = FreeModule(QQ, 3)
sage: v = M([1,2,3])
sage: w = v + 0
sage: w == v
True
sage: w is v
False

I feel quite torn about this, but slightly tend to the copy-version.

Addendum:

For FreeModule, we also have the following behavior:

sage: M = FreeModule(QQ, 3)
sage: M(0)
(0, 0, 0)
sage: M.zero()
(0, 0, 0)
sage: M.zero() is M(0)
False

I don't think that a parent should do that, especially when it already has a zero method. Should that be changed?

CC: @egourgoulhon @tscrim @mkoeppe

Component: misc

Branch/Commit: u/gh-mjungmath/return_copy_always_scalarfields @ b7b0b3f

Issue created by migration from https://trac.sagemath.org/ticket/30302

[mjungmath]

Description changed:

--- 
+++ 
@@ -1,6 +1,6 @@
 This question arose from ticket #30239, [comment:36](https://github.com/sagemath/sage/issues/30239#comment:36).
 
-Should `FiniteRankFreeModule` and manifold object always return a mutable copy, even for trivial operations? At least, this would be a consistent behavior.
+Should `FiniteRankFreeModule` and manifold objects always return a mutable copy, even for trivial operations? At least, this would be a consistent behavior.
 
 As pointed out by Matthias, this already holds true for `FreeModule`:
 

[mjungmath]

Description changed:

--- 
+++ 
@@ -14,7 +14,7 @@
 False
 ```
 
-I feel quite torn about this, but tend to the copy-version.
+I feel quite torn about this, but slightly tend to the copy-version.
 
 **Addendum:**
 

[mkoeppe]

comment:3

Both M() and M(0) return a new mutable element. That's how one creates a new vector, whose components can then be modified. If you want the immutable 0 element, you can use M.zero().

The arithmetic operations should either always create an immutable element or always create a new mutable element. It would be very inconvenient if the result of an operation was sometimes immutable, sometimes a new mutable element, depending on the input. (And it would be highly problematic if sometimes it would return an existing mutable element.)

[egourgoulhon]

comment:4

My concern here would be performance. For complicated symbolic expressions, creating a copy can have a significant CPU cost. The main concern is about scalar fields, because tensor fields arithmetics always end up to scalar fields arithmetics, the scalar fields being the tensor components in a given frame. For the time being, ScalarField._add_() starts with

        if self._is_zero:
            return other
        if other._is_zero:
            return self

If we change this to return a copy (I understand the arguments put forward by Matthias), then I would advocate for extensive benchmarks, with complicated tensor fields (those in the doctests are too simple), such as in those mentionned here.

[mkoeppe]

comment:5

Symbolic expressions are immutable and therefore do not need to be copied.

The cost of making a copy should be dominated by copying the dictionaries in the Components objects.

[mjungmath]

comment:6

If making a copy is so costly, we should try to optimize it.

If I remember correctly, until Sage 9 or so, all elements emerged from arithmetics were created from scratch, even the trivial ones. Even if copying would slow the current code down a bit, it should be still faster than versions before Sage 9.

[egourgoulhon]

comment:7

Replying to @mjungmath:

If making a copy is so costly, we should try to optimize it.

Copying a symbolic expression with hundreds of terms will always remain slower than returning self or other.

If I remember correctly, until Sage 9 or so, all elements emerged from arithmetics were created from scratch, even the trivial ones.

That's not true: already in Sage 7.4, we have the code snippet shown in comment:4.

[mjungmath]

comment:8

Replying to @mjungmath:

If I remember correctly, until Sage 9 or so, all elements emerged from arithmetics were created from scratch, even the trivial ones. Even if copying would slow the current code down a bit, it should be still faster than versions before Sage 9.

My memory was partially wrong. For scalar fields, the _is_zero variable was already implemented, for tensor fields it was not.

Furtermore, I've compared some computation times:

Current state:

sage: M = Manifold(2, 'M', structure='topological') # the 2-dimensional sphere S^2
sage: U = M.open_subset('U') # complement of the North pole
sage: c_xy.<x,y> = U.chart() # stereographic coordinates from the North pole
sage: V = M.open_subset('V') # complement of the South pole
sage: c_uv.<u,v> = V.chart() # stereographic coordinates from the South pole
sage: M.declare_union(U,V)   # S^2 is the union of U and V
sage: xy_to_uv = c_xy.transition_map(c_uv, (x/(x^2+y^2), y/(x^2+y^2)),
....:                                intersection_name='W',
....:                                restrictions1= x^2+y^2!=0,
....:                                restrictions2= u^2+v^2!=0)
sage: uv_to_xy = xy_to_uv.inverse()
sage: f = M.scalar_field({c_xy: 1/(1+x^2+y^2), c_uv: (u^2+v^2)/(1+u^2+v^2)},
....:                    name='f')
sage: %timeit f+0
The slowest run took 177.32 times longer than the fastest. This could mean that an intermediate result is being cached.
100000 loops, best of 5: 4.22 µs per loop
sage: %timeit f*1
The slowest run took 50.53 times longer than the fastest. This could mean that an intermediate result is being cached.
100000 loops, best of 5: 11.8 µs per loop

Returning a copy:

sage: M = Manifold(2, 'M', structure='topological') # the 2-dimensional sphere S^2
sage: U = M.open_subset('U') # complement of the North pole
sage: c_xy.<x,y> = U.chart() # stereographic coordinates from the North pole
sage: V = M.open_subset('V') # complement of the South pole
sage: c_uv.<u,v> = V.chart() # stereographic coordinates from the South pole
sage: M.declare_union(U,V)   # S^2 is the union of U and V
sage: xy_to_uv = c_xy.transition_map(c_uv, (x/(x^2+y^2), y/(x^2+y^2)),
....:                                intersection_name='W',
....:                                restrictions1= x^2+y^2!=0,
....:                                restrictions2= u^2+v^2!=0)
sage: uv_to_xy = xy_to_uv.inverse()
sage: f = M.scalar_field({c_xy: 1/(1+x^2+y^2), c_uv: (u^2+v^2)/(1+u^2+v^2)},
....:                    name='f')
sage: %timeit f+0
The slowest run took 23.78 times longer than the fastest. This could mean that an intermediate result is being cached.
10000 loops, best of 5: 60.6 µs per loop
sage: %timeit f*1
The slowest run took 24.43 times longer than the fastest. This could mean that an intermediate result is being cached.
10000 loops, best of 5: 33.4 µs per loop

That's a heavy loss of comutational time.

[mjungmath]

comment:9

Replying to @egourgoulhon:

Replying to @mjungmath:

If making a copy is so costly, we should try to optimize it.

Copying a symbolic expression with hundreds of terms will always remain slower than returning self or other.

If I remember correctly, until Sage 9 or so, all elements emerged from arithmetics were created from scratch, even the trivial ones.

That's not true: already in Sage 7.4, we have the code snippet shown in comment:4.

Yes sorry, I was mistaken. See above.

[mkoeppe]

comment:10

Replying to @mjungmath:

Returning a copy:

... could you share the changes that you made for benchmarking this?

[egourgoulhon]

comment:11

Replying to @mkoeppe:

Symbolic expressions are immutable and therefore do not need to be copied.

You are right. Actually they are not copied when copying chart functions (and hence would not be copied when copying scalar fields): the code of ChartFunction.copy() is indeed:

        resu = type(self)(self.parent())
        for kk, vv in self._express.items():
            resu._express[kk] = vv
        resu._expansion_symbol = self._expansion_symbol
        resu._order = self._order
        return resu

Here if kk is 'SR', vv is a symbolic expression.

[mjungmath]

Branch: u/gh-mjungmath/return_copy_always_scalarfields

[mkoeppe]

comment:13

Another implementation strategy if you consider Components objects an implementation detail that is not exposed to the user: Change them to "copy-on-write" semantics

---

New commits:

4f3b69f	Trac 30291: simple checks for trivial cases in `_mul_`, `_add_` and _div_
b7b0b3f	return copies