Opened 4 months ago

Closed 3 months ago

#30061 closed defect (fixed)

Speed up constructing high-dimensional Euclidean spaces

Reported by: mkoeppe Owned by:
Priority: major Milestone: sage-9.2
Component: geometry Keywords:
Cc: egourgoulhon, tscrim, nbruin, gh-mwageringel, gh-mjungmath Merged in:
Authors: Eric Gourgoulhon Reviewers: Matthias Koeppe
Report Upstream: N/A Work issues:
Branch: 3ce0c15 (Commits) Commit: 3ce0c15fc9b3d71842946ad6a0b0449771b1ada7
Dependencies: #30065, #30074 Stopgaps:

Description (last modified by mkoeppe)

The n-dimensional Euclidean space is available in Sage in many variants.

This ticket brings the speed of the most basic operation (constructing the space) of EuclideanSpace (from sage-manifolds) closer to that of the other variants.

Spaces without scalar product:

sage: VectorSpace(QQ, 5).category()
Category of finite dimensional vector spaces with basis over (number fields and quotient fields and metric spaces)
sage: %time VectorSpace(QQ, 5)
CPU times: user 28 µs, sys: 9 µs, total: 37 µs
Wall time: 40.1 µs
Vector space of dimension 5 over Rational Field
sage: %time VectorSpace(QQ, 80)
CPU times: user 218 µs, sys: 1 µs, total: 219 µs
Wall time: 223 µs
Vector space of dimension 80 over Rational Field
sage: %time VectorSpace(QQ, 1000)
CPU times: user 208 µs, sys: 1 µs, total: 209 µs
Wall time: 213 µs
Vector space of dimension 1000 over Rational Field

sage: for n in 5, 80, 1000, 4000, 10000, 100000: print("n = {}".format(n)); print("Construction: ", timeit("V = VectorSpace(QQ, {})".format(n),number=1,repeat=1)); u = [ x for x i
....: n range(n) ]; v = [ x + 1 for x in range(n) ]; V = VectorSpace(QQ, n); print("Distance:     ", timeit("norm(V(u) - V(v))"))
n = 5
Construction:  1 loop, best of 1: 56.1 ms per loop
Distance:      625 loops, best of 3: 81.7 μs per loop
n = 80
Construction:  1 loop, best of 1: 159 μs per loop
Distance:      625 loops, best of 3: 299 μs per loop
n = 1000
Construction:  1 loop, best of 1: 236 μs per loop
Distance:      125 loops, best of 3: 3.18 ms per loop
n = 4000
Construction:  1 loop, best of 1: 147 μs per loop
Distance:      25 loops, best of 3: 11.3 ms per loop
n = 10000
Construction:  1 loop, best of 1: 149 μs per loop
Distance:      25 loops, best of 3: 28.7 ms per loop
n = 100000
Construction:  1 loop, best of 1: 162 μs per loop
Distance:      5 loops, best of 3: 296 ms per loop
sage: CombinatorialFreeModule(QQ, range(5)).category()
Category of finite dimensional vector spaces with basis over Rational Field
sage: %time CombinatorialFreeModule(QQ, range(5))
CPU times: user 80 µs, sys: 0 ns, total: 80 µs
Wall time: 86.1 µs
Free module generated by {0, 1, 2, 3, 4} over Rational Field
sage: %time CombinatorialFreeModule(QQ, range(80))
CPU times: user 244 µs, sys: 10 µs, total: 254 µs
Wall time: 259 µs
Free module generated by {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79} over Rational Field
sage: %time CombinatorialFreeModule(QQ, range(1000))
CPU times: user 322 µs, sys: 26 µs, total: 348 µs
Wall time: 354 µs
sage: AffineSpace(RR, 5).category()
Category of schemes over Real Field with 53 bits of precision

sage: %time AffineSpace(RR, 5)
CPU times: user 47 µs, sys: 1 µs, total: 48 µs
Wall time: 51 µs
Affine Space of dimension 5 over Real Field with 53 bits of precision
sage: %time AffineSpace(RR, 80)
CPU times: user 2 ms, sys: 39 µs, total: 2.04 ms
Wall time: 2.04 ms
Affine Space of dimension 80 over Real Field with 53 bits of precision
sage: %time AffineSpace(RR, 1000)
CPU times: user 62.8 ms, sys: 1.45 ms, total: 64.3 ms
Wall time: 63.5 ms
Affine Space of dimension 1000 over Real Field with 53 bits of precision

Spaces with scalar product:

sage: VectorSpace(QQ, 5, inner_product_matrix=matrix.identity(5)).category()
Category of finite dimensional vector spaces with basis over (number fields and quotient fields and metric spaces)
sage: %time VectorSpace(QQ, 5, inner_product_matrix=matrix.identity(5))
CPU times: user 6.55 ms, sys: 666 µs, total: 7.22 ms
Wall time: 6.88 ms
Ambient quadratic space of dimension 5 over Rational Field
Inner product matrix:
[1 0 0 0 0]
[0 1 0 0 0]
[0 0 1 0 0]
[0 0 0 1 0]
[0 0 0 0 1]
sage: %time VectorSpace(QQ, 80, inner_product_matrix=matrix.identity(80))
CPU times: user 14.1 ms, sys: 439 µs, total: 14.6 ms
Wall time: 14.4 ms
Ambient quadratic space of dimension 80 over Rational Field
Inner product matrix:
sage: %time VectorSpace(QQ, 1000, inner_product_matrix=matrix.identity(1000))
CPU times: user 1.44 s, sys: 34.4 ms, total: 1.47 s
Wall time: 1.48 s
Ambient quadratic space of dimension 1000 over Rational Field
Inner product matrix: ...

sage: for n in 5, 80, 1000, 4000: print("n = {}".format(n)); print("Construction: ", timeit("V = VectorSpace(QQ, {}, inner_product_matrix=matrix.identity({}))".format(n, n),number=1,repeat=1)); u = [ x for x in range(n) ]; v = [ x + 1 for x in range(n) ]; V = VectorSpace(QQ, n, inner_product_matrix=matrix.identity(n)); print("Distance:     ", timeit("t = V(u)-V(v); sqrt(t.inner_product(t))"))
n = 5
Construction:  1 loop, best of 1: 61.9 ms per loop
Distance:      625 loops, best of 3: 49.4 μs per loop
n = 80
Construction:  1 loop, best of 1: 9.75 ms per loop
Distance:      625 loops, best of 3: 243 μs per loop
n = 1000
Construction:  1 loop, best of 1: 1.51 s per loop
Distance:      25 loops, best of 3: 14.4 ms per loop
n = 4000
Construction:  1 loop, best of 1: 23.8 s per loop
Distance:      5 loops, best of 3: 225 ms per loop

NB: It is not in the category MetricSpaces, and thus the element methods dist and abs (?!) are missing... The methods __abs__, norm, and dot_product are unrelated to the inner product matrix; only inner_product uses inner_product_matrix.

sage: EuclideanSpace(5).category()
Category of smooth manifolds over Real Field with 53 bits of precision

sage: %time EuclideanSpace(5)
CPU times: user 177 ms, sys: 10.4 ms, total: 187 ms
Wall time: 196 ms
5-dimensional Euclidean space E^5
sage: %time EuclideanSpace(80)
CPU times: user 32.8 s, sys: 758 ms, total: 33.6 s
Wall time: 31.5 s
80-dimensional Euclidean space E^80
sage: %time EuclideanSpace(1000)
(timeout)

With this ticket (and its dependencies #30065, #30074):

sage: %time EuclideanSpace(5)
CPU times: user 4.87 ms, sys: 372 µs, total: 5.24 ms
Wall time: 4.93 ms
5-dimensional Euclidean space E^5
sage: %time EuclideanSpace(80)
CPU times: user 106 ms, sys: 4.52 ms, total: 110 ms
Wall time: 92 ms
80-dimensional Euclidean space E^80
sage: %time EuclideanSpace(1000)
CPU times: user 1.04 s, sys: 33.5 ms, total: 1.07 s
Wall time: 986 ms

Some caching is happening too. The second time:

sage: %time EuclideanSpace(1000)
CPU times: user 207 ms, sys: 5.63 ms, total: 213 ms
Wall time: 212 ms
1000-dimensional Euclidean space E^1000

Scalar product without a space:

sage: %time DiagonalQuadraticForm(QQ, [1]*5)
CPU times: user 60 µs, sys: 1e+03 ns, total: 61 µs
Wall time: 62.9 µs
Quadratic form in 5 variables over Rational Field with coefficients: 
[ 1 0 0 0 0 ]
[ * 1 0 0 0 ]
[ * * 1 0 0 ]
[ * * * 1 0 ]
[ * * * * 1 ]
sage: %time DiagonalQuadraticForm(QQ, [1]*80)
CPU times: user 1.35 ms, sys: 31 µs, total: 1.39 ms
Wall time: 1.44 ms
Quadratic form in 80 variables over Rational Field with coefficients: 
sage: %time DiagonalQuadraticForm(QQ, [1]*1000)
CPU times: user 144 ms, sys: 2.85 ms, total: 147 ms
Wall time: 147 ms
Quadratic form in 1000 variables over Rational Field with coefficients: 

Change History (70)

comment:1 Changed 4 months ago by mkoeppe

Most of the time spent seems to come from SR (Maxima). Is it possible to speed things up for this easy special case, avoiding symbolic computation completely?

comment:2 follow-ups: Changed 4 months ago by egourgoulhon

EuclideanSpace(n) creates the following objects:

  1. a real smooth manifold E of dimension n
  2. a chart X on E, representing the Cartesian coordinates
  3. the coordinate vector frame associated with X
  4. a Riemannian metric g on E, with its components w.r.t. X initialized to diag(1,1,...,1)
  5. the Levi-Civita connection associated with g, with its components w.r.t. X initialized to 0

Most of the CPU time is spent in step 2, actually in creating the symbolic variables (elements of SR) representing the Cartesian coordinates. The latter operation is equivalent to

sage: def create_coords(n):
....:     coords = SR.var(["x{}".format(i) for i in range(n)], domain='real')
....:     for x in coords:
....:         assume(x, 'real')

and we have:

sage: %time EuclideanSpace(5)
CPU times: user 345 ms, sys: 28 ms, total: 373 ms
Wall time: 373 ms
5-dimensional Euclidean space E^5
sage: %time create_coords(5)
CPU times: user 456 ms, sys: 338 µs, total: 456 ms
Wall time: 402 ms
sage: %time EuclideanSpace(20)
CPU times: user 3.66 s, sys: 58.9 ms, total: 3.72 s
Wall time: 3.48 s
20-dimensional Euclidean space E^20
sage: %time create_coords(20)
CPU times: user 3.38 s, sys: 75.3 ms, total: 3.45 s
Wall time: 3.21 s
sage: %time EuclideanSpace(40)
CPU times: user 12.6 s, sys: 253 ms, total: 12.8 s
Wall time: 11.9 s
40-dimensional EuclideanSR.var(["x{}".format(i) for i in range(n)], space E^40
sage: %time create_coords(40)
CPU times: user 12.2 s, sys: 226 ms, total: 12.5 s
Wall time: 11.4 s

Actually, it turns out that what takes most of CPU time is demanding that the symbolic variables are real. Indeed, if we skip domain='real' and assume(x, 'real') (each uses roughly half of the CPU time and, although they look redundant, both are actually necessary for efficient simplifications), we get very small times:

sage: %time coords = SR.var(["x{}".format(i) for i in range(40)])
CPU times: user 336 µs, sys: 9 µs, total: 345 µs
Wall time: 350 µs

I don't know why SR.var('x', domain='real') and assume(x, 'real') are so slow...

NB: in the current setting, chart coordinates are stored as SR variables; but since SymPy can be used as the symbolic backend on manifolds, via

sage: E = EuclideanSpace(5)
sage: E.set_calculus_method('sympy')

chart coordinates could be stored as SymPy variables, instead of SR ones.

Last edited 4 months ago by egourgoulhon (previous) (diff)

comment:3 in reply to: ↑ 2 Changed 4 months ago by mkoeppe

Thanks very much for the explanation and analysis!

Replying to egourgoulhon:

I don't know why SR.var('x', domain='real') and assume(x, 'real') are so slow...

I have created #30065 for this

comment:4 in reply to: ↑ 2 ; follow-up: Changed 4 months ago by mkoeppe

Replying to egourgoulhon:

in the current setting, chart coordinates are stored as SR variables; but since SymPy can be used as the symbolic backend on manifolds, via

sage: E = EuclideanSpace(5)
sage: E.set_calculus_method('sympy')

chart coordinates could be stored as SymPy variables, instead of SR ones.

Thanks, I'll try this.

Actually, would it make sense (at least for simple cases) to have a calculus method that only uses Sage's rational functions?

comment:5 follow-up: Changed 4 months ago by mkoeppe

In a related direction, would it be of interest to have a category of "algebraic" differential manifolds, as a differential geometry view on real varieties and semialgebraic sets (and their boundaries)?

comment:6 in reply to: ↑ 4 ; follow-ups: Changed 4 months ago by egourgoulhon

Replying to mkoeppe:

Replying to egourgoulhon:

in the current setting, chart coordinates are stored as SR variables; but since SymPy can be used as the symbolic backend on manifolds, via

sage: E = EuclideanSpace(5)
sage: E.set_calculus_method('sympy')

chart coordinates could be stored as SymPy variables, instead of SR ones.

Thanks, I'll try this.

Actually, would it make sense (at least for simple cases) to have a calculus method that only uses Sage's rational functions?

Yes, as soon as you don't have to take a square root, as for instance in evaluating the volume n-form on a pseudo-Riemannian manifold (aka Levi-Civita tensor). In particular, this would forbid the computation of the triple scalar product within spherical coordinates in the Euclidean 3-space.

Actually, it should be relatively easy to add a new calculus method, in addition to the two ones already implemented (SR and SymPy), via the classes CalculusMethod and ChartFunction.

Last edited 4 months ago by egourgoulhon (previous) (diff)

comment:7 in reply to: ↑ 5 ; follow-up: Changed 4 months ago by egourgoulhon

Replying to mkoeppe:

In a related direction, would it be of interest to have a category of "algebraic" differential manifolds, as a differential geometry view on real varieties and semialgebraic sets (and their boundaries)?

Yes, I think so.

comment:8 in reply to: ↑ 6 Changed 4 months ago by mkoeppe

Replying to egourgoulhon:

Actually, would it make sense (at least for simple cases) to have a calculus method that only uses Sage's rational functions?

Yes, as soon as you don't have to take a square root, as for instance in evaluating the volume n-form on a pseudo-Riemannian manifold (aka Levi-Civita tensor). In particular, this would forbid the computation of the triple scalar product within spherical coordinates in the Euclidean 3-space.

Actually, it should be relatively easy to add a new calculus method, in addition to the two ones already implemented (SR and SymPy), via the classes CalculusMethod and ChartFunction.

Thanks! I have created #30070 for this.

comment:9 in reply to: ↑ 7 ; follow-up: Changed 4 months ago by mkoeppe

Replying to egourgoulhon:

Replying to mkoeppe:

In a related direction, would it be of interest to have a category of "algebraic" differential manifolds, as a differential geometry view on real varieties and semialgebraic sets (and their boundaries)?

Yes, I think so.

OK, I have created #30069 for the case of real algebraic manifolds.

Does sage-manifolds currently implement manifolds with boundary, or is it planned to implement them?

comment:10 in reply to: ↑ 6 Changed 4 months ago by egourgoulhon

Replying to egourgoulhon:

Replying to mkoeppe:

Actually, would it make sense (at least for simple cases) to have a calculus method that only uses Sage's rational functions?

Yes, as soon as you don't have to take a square root, as for instance in evaluating the volume n-form on a pseudo-Riemannian manifold (aka Levi-Civita tensor). In particular, this would forbid the computation of the triple scalar product within spherical coordinates in the Euclidean 3-space.

Actually, for spherical coordinates, there is already some issue with rational functions at the level of the metric tensor itself, since the components of the Euclidean metric involve the sine function.

comment:11 in reply to: ↑ 9 ; follow-up: Changed 4 months ago by egourgoulhon

Replying to mkoeppe:

OK, I have created #30069 for the case of real algebraic manifolds.

Very good.

Does sage-manifolds currently implement manifolds with boundary, or is it planned to implement them?

No, manifolds with boundary are not implemented yet; it would be nice to have them.

comment:12 Changed 4 months ago by mkoeppe

With #30065 and #30074:

sage: %prun EuclideanSpace(1000)

         53513603 function calls (53508842 primitive calls) in 37.037 seconds

   Ordered by: internal time

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
  2002006    4.658    0.000    4.786    0.000 tensorfield.py:684(_del_restrictions)
  3003010    3.263    0.000    3.766    0.000 tensorfield.py:667(_del_derived)
  2002000    3.076    0.000   18.325    0.000 free_module_tensor.py:1191(_set_comp_unsafe)
  2002000    3.047    0.000   25.990    0.000 tensorfield_paral.py:820(set_comp)
  2002000    2.571    0.000   22.942    0.000 free_module_tensor.py:1263(set_comp)
  2002000    2.478    0.000    7.194    0.000 comp.py:864(__setitem__)
  2002000    2.147    0.000    2.147    0.000 free_module_tensor.py:1474(del_other_comp)
  2003000    2.029    0.000    2.398    0.000 comp.py:616(_check_indices)
  2002000    1.825    0.000   20.150    0.000 tensorfield_paral.py:730(_set_comp_unsafe)
  2002006    1.697    0.000    9.256    0.000 tensorfield_paral.py:708(_del_derived)
        1    1.532    1.532   36.171   36.171 free_module_basis.py:566(__init__)
  1002000    1.095    0.000    1.108    0.000 scalarfield.py:1154(is_trivial_zero)
  1001003    0.931    0.000    7.942    0.000 vectorfield.py:1612(_del_derived)
  8071534    0.813    0.000    0.813    0.000 {built-in method builtins.isinstance}
        1    0.722    0.722   15.342   15.342 free_module_basis.py:375(__init__)
  1001003    0.675    0.000    0.675    0.000 vectorfield.py:292(_del_dependencies)
  2004002    0.656    0.000    0.656    0.000 finite_rank_free_module.py:2037(irange)
  1001000    0.603    0.000    5.144    0.000 diff_form.py:1241(_del_derived)
     6800    0.532    0.000    0.598    0.000 maxima_lib.py:412(_eval_line)
  2045310    0.378    0.000    0.378    0.000 {built-in method builtins.hasattr}
  1001003    0.335    0.000    5.137    0.000 multivectorfield.py:966(_del_derived)
  5005023    0.334    0.000    0.334    0.000 {method 'clear' of 'dict' objects}

comment:13 follow-up: Changed 4 months ago by mkoeppe

  • Cc tscrim added
  • Description modified (diff)
  • Summary changed from EuclideanSpace constructor is slow to Constructing Euclidean spaces, many variants

comment:14 Changed 4 months ago by mkoeppe

  • Description modified (diff)

comment:15 Changed 4 months ago by egourgoulhon

  • Description modified (diff)

comment:16 in reply to: ↑ 13 Changed 4 months ago by egourgoulhon

Replying to mkoeppe:

Regarding the new ticket description, I would not say that the first three examples, namely VectorSpace(QQ, 5), CombinatorialFreeModule(QQ, range(5)) and AffineSpace(RR, 5) do construct an Euclidean space, because their ouputs are not endowed with a scalar product, are they?

comment:17 Changed 4 months ago by mkoeppe

  • Description modified (diff)

Of course, I agree. The summary already showed the categories of the spaces, but I've added two headers to emphasize this point.

comment:18 in reply to: ↑ 11 Changed 4 months ago by mkoeppe

Replying to egourgoulhon:

Replying to mkoeppe:

Does sage-manifolds currently implement manifolds with boundary, or is it planned to implement them?

No, manifolds with boundary are not implemented yet; it would be nice to have them.

I have created #30080 for this

comment:19 Changed 4 months ago by mkoeppe

  • Description modified (diff)

comment:20 Changed 4 months ago by mkoeppe

  • Dependencies set to #30065, #30074

comment:21 follow-up: Changed 4 months ago by mkoeppe

For further speedups of EuclideanSpace(1000), it seems one would need to work on sage.tensor.modules.free_module_basis.FreeModuleCoBasis.__init__ and FreeModuleBasis.__init__. They dominate the computation now, using non-sparse operations, calling sage.manifolds.differentiable.tensorfield_paral.TensorFieldParal.set_comp 10002 times. Does it make sense to make this sparse in some way? Would major speedups be expected from cythonizing this module?

comment:22 in reply to: ↑ 21 Changed 4 months ago by egourgoulhon

Replying to mkoeppe:

For further speedups of EuclideanSpace(1000), it seems one would need to work on sage.tensor.modules.free_module_basis.FreeModuleCoBasis.__init__ and FreeModuleBasis.__init__. They dominate the computation now, using non-sparse operations, calling sage.manifolds.differentiable.tensorfield_paral.TensorFieldParal.set_comp 10002 times. Does it make sense to make this sparse in some way?

Yes, actually in the current framework, any uninitialized component of a vector/tensor is considered to be zero. So in FreeModuleBasis.__init__, we could skip the loop initializing the components of the basis vector to zero:

        for i in fmodule.irange():
            v = fmodule.element_class(fmodule)
-           for j in fmodule.irange():
-               v.set_comp(self)[j] = fmodule._ring.zero()
            v.set_comp(self)[i] = fmodule._ring.one()
            vl.append(v)

This would make the number of calls to TensorFieldParal.set_comp to be 1000, instead in of 10002 ! Similarly in FreeModuleCoBasis.__init_:

        for i in self._fmodule.irange():
            v = self._fmodule.linear_form()
-           for j in self._fmodule.irange():
-               v.set_comp(basis)[j] = 0
-           v.set_comp(basis)[i] = 1
+           v.set_comp(basis)[i] = self._fmodule._ring.one()
            vl.append(v)

The above unnecessary initializations to zero have been implemented at the early stage of the manifolds project, when we were not certain to keep the storage convention of having uninitialized components to be zero. Given the heavy use of that convention now, it's pretty safe to remove these lines.

Last edited 4 months ago by egourgoulhon (previous) (diff)

comment:23 follow-up: Changed 4 months ago by egourgoulhon

  • Branch set to public/manifolds/init_module_basis-30061
  • Commit set to fccbf438290e2404a79a4c7a650721120157a6c0

The changes suggested in comment:22 are implemented in the above branch.


New commits:

fccbf43Improve initialization of elements of a free module basis

comment:24 Changed 4 months ago by git

  • Commit changed from fccbf438290e2404a79a4c7a650721120157a6c0 to 3ce0c15fc9b3d71842946ad6a0b0449771b1ada7

Branch pushed to git repo; I updated commit sha1. New commits:

134da39sage.symbolic.assumptions.GenericDeclaration.assume: Validate against cached valid features first
1fbfb68Update instead of overwriting
a99aa3fMerge branch 't/30065/faster_maxima_assumptions' into t/30061/public/manifolds/init_module_basis-30061
8b31261sage.symbolic.assumptions.GenericDeclaration: Make it a UniqueRepresentation
51ea23dsage.symbolic.assumptions.GenericDeclaration.assume: Check first if already in _assumptions
9119e82sage.symbolic.assumptions._assumptions: Change from list to OrderDict
cda81ffsage.symbolic.assumptions.GenericDeclaration: Make it hashable by inheriting comparisons from UniqueRepresentation
642836dsage.symbolic.assumptions: Use dict instead of OrderedDict for _assumptions
3ce0c15Merge branch 't/30074/even_faster_maxima_assumptions' into t/30061/public/manifolds/init_module_basis-30061

comment:25 Changed 4 months ago by mkoeppe

Merged in the tickets listed as dependencies

comment:26 in reply to: ↑ 23 Changed 4 months ago by mkoeppe

Replying to egourgoulhon:

The changes suggested in comment:22 are implemented in the above branch.

Great! Together with the merged tickets, this now gives:

sage: %time EuclideanSpace(5)
CPU times: user 4.87 ms, sys: 372 µs, total: 5.24 ms
Wall time: 4.93 ms
5-dimensional Euclidean space E^5
sage: %time EuclideanSpace(80)
CPU times: user 106 ms, sys: 4.52 ms, total: 110 ms
Wall time: 92 ms
80-dimensional Euclidean space E^80
sage: %time EuclideanSpace(1000)
CPU times: user 1.04 s, sys: 33.5 ms, total: 1.07 s
Wall time: 986 ms

comment:27 Changed 4 months ago by mkoeppe

And some caching is happening too. The second time:

sage: %time EuclideanSpace(1000)
CPU times: user 207 ms, sys: 5.63 ms, total: 213 ms
Wall time: 212 ms
1000-dimensional Euclidean space E^1000

comment:28 Changed 4 months ago by mkoeppe

That's faster than VectorSpace(QQ, 1000, inner_product_matrix=matrix.identity(1000))

comment:29 Changed 4 months ago by mkoeppe

Most of the caching is coming from #30074. Making the real-declared symbolic variables known to Maxima is costly the first time; after that, it's free.

comment:30 Changed 4 months ago by mkoeppe

First time:

sage: %prun EuclideanSpace(10000)
         4614917 function calls (4585885 primitive calls) in 42.685 seconds

   Ordered by: internal time

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
    80096   35.217    0.000   35.779    0.000 maxima_lib.py:412(_eval_line)
       25    0.529    0.021    0.529    0.021 {built-in method sage.libs.ecl.ecl_eval}
    10000    0.425    0.000   39.373    0.004 {method 'var' of 'sage.symbolic.ring.SymbolicRing' objects}
    20000    0.377    0.000   38.550    0.002 assumptions.py:213(assume)
    20048    0.323    0.000    0.382    0.000 maxima_lib.py:274(max_to_string)
    20006    0.275    0.000    0.288    0.000 free_module_tensor.py:258(__init__)
    20000    0.208    0.000    0.460    0.000 assumptions.py:396(preprocess_assumptions)
       56    0.207    0.004    0.207    0.004 {built-in method _imp.create_dynamic}
    50000    0.190    0.000    0.190    0.000 latex.py:463(__add__)
   209119    0.165    0.000    0.165    0.000 {built-in method builtins.hasattr}
    20007    0.164    0.000    0.180    0.000 tensorfield_paral.py:689(_init_derived)
  10032/1    0.139    0.000   42.050   42.050 {sage.misc.classcall_metaclass.typecall}
    10139    0.129    0.000    0.130    0.000 {built-in method builtins.repr}
    20001    0.125    0.000    0.143    0.000 interface.py:620(__getattr__)
      714    0.119    0.000    0.119    0.000 {built-in method marshal.loads}
   546190    0.114    0.000    0.114    0.000 {built-in method builtins.isinstance}
        1    0.102    0.102   39.733   39.733 chart.py:1610(_init_coordinates)
    10003    0.091    0.000    0.107    0.000 chart_func.py:318(__init__)
    10000    0.091    0.000    0.460    0.000 {method '_maxima_init_' of 'sage.structure.sage_object.SageObject' objects}
    40024    0.090    0.000    0.099    0.000 interface.py:491(_next_var_name)
    20000    0.087    0.000   13.550    0.001 interface.py:587(function_call)
    40024    0.086    0.000   15.460    0.000 maxima_lib.py:492(set)
    30024    0.086    0.000   13.245    0.000 interface.py:251(__call__)
    30024    0.082    0.000   13.139    0.000 interface.py:706(__init__)
    30000    0.080    0.000    0.237    0.000 scalarfield.py:1154(is_trivial_zero)
    20000    0.079    0.000    3.032    0.000 interface.py:534(_convert_args_kwds)
    20024    0.078    0.000   16.639    0.001 interface.py:981(__del__)
    20000    0.077    0.000    0.416    0.000 free_module_tensor.py:1191(_set_comp_unsafe)
   160195    0.074    0.000    0.074    0.000 {method 'find' of 'str' objects}
    30000    0.069    0.000    0.126    0.000 chart_func.py:810(is_trivial_zero)
    10000    0.065    0.000    0.118    0.000 diff_form.py:1241(_del_derived)
    10001    0.061    0.000    0.069    0.000 interface.py:419(_relation_symbols)
    20048    0.059    0.000    0.059    0.000 {method 'python' of 'sage.libs.ecl.EclObject' objects}
    90434    0.057    0.000    0.057    0.000 {method 'split' of 'str' objects}
    20006    0.051    0.000    0.053    0.000 tensorfield.py:684(_del_restrictions)
    30024    0.051    0.000   13.056    0.000 maxima_lib.py:561(_create)
    56/50    0.049    0.001    0.054    0.001 {built-in method _imp.exec_dynamic}
    20000    0.045    0.000    0.045    0.000 infinity.py:1020(gen)
    30010    0.044    0.000    0.051    0.000 tensorfield.py:667(_del_derived)
    20000    0.044    0.000   21.876    0.001 interface.py:653(__call__)
    20001    0.044    0.000    0.070    0.000 vectorfield_module.py:1765(dual_exterior_power)
    20000    0.042    0.000    0.251    0.000 comp.py:864(__setitem__)
    30000    0.042    0.000    0.050    0.000 comp.py:616(_check_indices)
    20000    0.041    0.000    0.043    0.000 interface.py:603(_function_call_string)
    10001    0.040    0.000    0.332    0.000 vectorfield.py:1533(__init__)
    10000    0.040    0.000    0.224    0.000 expression_conversions.py:155(__call__)
    20024    0.039    0.000   16.518    0.001 maxima_lib.py:516(clear)

comment:31 Changed 4 months ago by mkoeppe

Second time:

         2221257 function calls (2221199 primitive calls) in 2.982 seconds

   Ordered by: internal time

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
    50000    0.297    0.000    0.297    0.000 latex.py:463(__add__)
    30010    0.243    0.000    0.251    0.000 tensorfield.py:667(_del_derived)
    20006    0.214    0.000    0.226    0.000 free_module_tensor.py:258(__init__)
    10000    0.178    0.000    0.266    0.000 {method 'var' of 'sage.symbolic.ring.SymbolicRing' objects}
    10003    0.091    0.000    0.107    0.000 chart_func.py:318(__init__)
    30000    0.086    0.000    0.254    0.000 scalarfield.py:1154(is_trivial_zero)
    20000    0.085    0.000    0.656    0.000 free_module_tensor.py:1191(_set_comp_unsafe)
    10000    0.077    0.000    0.136    0.000 diff_form.py:1241(_del_derived)
   130046    0.075    0.000    0.075    0.000 {built-in method builtins.hasattr}
    30000    0.074    0.000    0.134    0.000 chart_func.py:810(is_trivial_zero)
    20007    0.073    0.000    0.090    0.000 tensorfield_paral.py:689(_init_derived)
    20000    0.068    0.000    0.085    0.000 assumptions.py:396(preprocess_assumptions)
    20006    0.055    0.000    0.057    0.000 tensorfield.py:684(_del_restrictions)
    20001    0.049    0.000    0.076    0.000 vectorfield_module.py:1765(dual_exterior_power)
   270213    0.048    0.000    0.048    0.000 {built-in method builtins.isinstance}
    30000    0.045    0.000    0.054    0.000 comp.py:616(_check_indices)
        1    0.044    0.044    0.412    0.412 chart.py:1610(_init_coordinates)
    20000    0.044    0.000    0.271    0.000 comp.py:864(__setitem__)
    10001    0.042    0.000    0.276    0.000 vectorfield.py:1533(__init__)
    30001    0.040    0.000    0.040    0.000 chart_func.py:460(expr)
    20000    0.034    0.000    0.714    0.000 free_module_tensor.py:1263(set_comp)
    10001    0.034    0.000    0.243    0.000 diff_form.py:1137(__init__)
    20000    0.033    0.000    0.746    0.000 tensorfield_paral.py:820(set_comp)
    30000    0.033    0.000    0.033    0.000 {method '_latex_' of 'sage.symbolic.expression.Expression' objects}
    10001    0.033    0.000    0.059    0.000 free_module_alt_form.py:298(_new_comp)
    30000    0.033    0.000    0.033    0.000 free_module_tensor.py:923(set_name)
    10001    0.028    0.000    0.054    0.000 free_module_element.py:230(_new_comp)
    10000    0.027    0.000    0.065    0.000 comp.py:3061(_ordered_indices)
    20000    0.026    0.000    0.026    0.000 free_module_tensor.py:1474(del_other_comp)
    10001    0.026    0.000    0.026    0.000 chart.py:311(<genexpr>)
    10000    0.025    0.000    0.178    0.000 comp.py:4784(__setitem__)
   120055    0.025    0.000    0.042    0.000 {built-in method builtins.len}
    10001    0.024    0.000    0.024    0.000 euclidean.py:858(<genexpr>)
    30010    0.024    0.000    0.024    0.000 tensorfield.py:654(_init_derived)
    10001    0.023    0.000    0.023    0.000 vectorframe.py:1764(<genexpr>)
    10001    0.023    0.000    0.038    0.000 vectorfield_module.py:1709(exterior_power)
        3    0.022    0.007    0.339    0.113 free_module_basis.py:208(set_name)
    20005    0.022    0.000    0.052    0.000 comp.py:496(__init__)
        1    0.022    0.022    1.688    1.688 free_module_basis.py:566(__init__)
    20006    0.022    0.000    0.318    0.000 tensorfield_paral.py:708(_del_derived)

comment:32 Changed 4 months ago by mkoeppe

  • Cc nbruin added

comment:33 Changed 4 months ago by mkoeppe

First time profile, filtered to methods in maxima_lib and sage.libs.ecl:

sage: %prun EuclideanSpace(5000)
    40096    7.970    0.000    8.203    0.000 maxima_lib.py:412(_eval_line)
       25    0.511    0.020    0.511    0.020 {built-in method sage.libs.ecl.ecl_eval}
    10048    0.133    0.000    0.158    0.000 maxima_lib.py:274(max_to_string)
    20024    0.037    0.000    3.505    0.000 maxima_lib.py:492(set)
    15024    0.021    0.000    2.934    0.000 maxima_lib.py:561(_create)
    10024    0.017    0.000    3.597    0.000 maxima_lib.py:516(clear)
    40096    0.015    0.000    0.017    0.000 maxima_lib.py:463(<listcomp>)
     5023    0.007    0.000    0.275    0.000 maxima_lib.py:541(get)

comment:34 follow-up: Changed 4 months ago by mkoeppe

The time spent on maxima_lib can probably still be reduced a little by using ECL objects directly for "context", but I will not work on this at the moment.

comment:35 in reply to: ↑ 34 Changed 4 months ago by mkoeppe

Replying to mkoeppe:

The time spent on maxima_lib can probably still be reduced a little by using ECL objects directly for "context", but I will not work on this at the moment.

This is #30086 (just for saving the thought and some preliminary work on it).

comment:36 Changed 4 months ago by mkoeppe

  • Authors set to Eric Gourgoulhon

comment:37 Changed 4 months ago by mkoeppe

Next I think we can benchmark some other basic operations in the spaces

comment:38 Changed 4 months ago by mkoeppe

  • Description modified (diff)

comment:39 Changed 4 months ago by mkoeppe

  • Description modified (diff)

comment:40 follow-up: Changed 4 months ago by mkoeppe

Does EuclideanSpace know how to compute geodesics, and thus distances?

comment:41 in reply to: ↑ 40 ; follow-up: Changed 4 months ago by egourgoulhon

Replying to mkoeppe:

Does EuclideanSpace know how to compute geodesics, and thus distances?

Yes it does, via IntegratedGeodesic. Note however that for Euclidean spaces, geodesics are trivial (straight lines).

comment:42 Changed 4 months ago by egourgoulhon

By the way, an overview of EuclideanSpace capabilities is in this thematic tutorial; see also the associated notebooks and this section of the reference manual.

comment:43 Changed 4 months ago by mkoeppe

Thanks for the pointers. It's on my list of things to learn more systematically about sage-manifolds.

comment:44 in reply to: ↑ 41 ; follow-up: Changed 4 months ago by mkoeppe

Replying to egourgoulhon:

Replying to mkoeppe:

Does EuclideanSpace know how to compute geodesics, and thus distances?

Yes it does, via IntegratedGeodesic.

Should it (and presumably other manifold classes) be added to the category MetricSpaces then (#30062)?

comment:45 in reply to: ↑ 44 ; follow-up: Changed 4 months ago by egourgoulhon

Replying to mkoeppe:

Should it (and presumably other manifold classes) be added to the category MetricSpaces then (#30062)?

Yes. It will be easy to implement the distance function by means of the Cartesian chart.

comment:46 Changed 4 months ago by mkoeppe

  • Cc gh-mwageringel added

comment:47 in reply to: ↑ 45 Changed 4 months ago by mkoeppe

Replying to egourgoulhon:

Replying to mkoeppe:

Should it (and presumably other manifold classes) be added to the category MetricSpaces then (#30062)?

Yes. It will be easy to implement the distance function by means of the Cartesian chart.

Thanks. I have added this task to #30062.

comment:48 Changed 3 months ago by mkoeppe

  • Cc gh-mjungmath added

comment:49 follow-up: Changed 3 months ago by mkoeppe

Should we set this ticket to "needs_review"?

comment:50 follow-up: Changed 3 months ago by gh-mjungmath

Can you please briefly summarize in the description what has actually been done and what the new benefits are? It is very difficult for me to piece that all together just from the comments and cross references to other tickets. Thanks! :)

In general, I would wait until the dependent tickets are merged. However, I am a bad role model due to my impatience.

Last edited 3 months ago by gh-mjungmath (previous) (diff)

comment:51 Changed 3 months ago by egourgoulhon

  • Status changed from new to needs_review

comment:52 in reply to: ↑ 49 Changed 3 months ago by egourgoulhon

Replying to mkoeppe:

Should we set this ticket to "needs_review"?

Yes, indeed.

comment:53 in reply to: ↑ 50 ; follow-up: Changed 3 months ago by egourgoulhon

Replying to gh-mjungmath:

Can you please briefly summarize in the description what has actually been done and what the new benefits are? It is very difficult for me to piece that all together just from the comments and cross references to other tickets. Thanks! :)

The code added in this ticket is described in comment:22 and comment:23. The other code in the ticket branch arises from the dependencies tickets #30065 and #30074.

comment:54 in reply to: ↑ 53 Changed 3 months ago by gh-mjungmath

Replying to egourgoulhon:

The code added in this ticket is described in comment:22 and comment:23. The other code in the ticket branch arises from the dependencies tickets #30065 and #30074.

Thanks!

comment:55 Changed 3 months ago by mkoeppe

  • Summary changed from Constructing Euclidean spaces, many variants to Speed up constructing high-dimensional Euclidean spaces

comment:56 Changed 3 months ago by mkoeppe

  • Description modified (diff)

comment:57 Changed 3 months ago by mkoeppe

  • Description modified (diff)

comment:58 Changed 3 months ago by mkoeppe

  • Description modified (diff)

comment:59 Changed 3 months ago by mkoeppe

  • Description modified (diff)

comment:60 Changed 3 months ago by mkoeppe

  • Reviewers set to Matthias Koeppe
  • Status changed from needs_review to positive_review

comment:61 Changed 3 months ago by gh-mjungmath

  • Status changed from positive_review to needs_work

Patchbot not green. The latest branch of #30074 has not been merged yet.

comment:62 follow-up: Changed 3 months ago by tscrim

  • Status changed from needs_work to positive_review

That is not a real error (the other patchbot also came back green). You don't need to have all of the other branches merged in.

comment:63 in reply to: ↑ 62 ; follow-up: Changed 3 months ago by gh-mjungmath

Replying to tscrim:

That is not a real error (the other patchbot also came back green).

It is a real error. Pyflakes is complaining about an unused import.

You don't need to have all of the other branches merged in.

Why is that? Wouldn't that cause a merge conflict later on?

Last edited 3 months ago by gh-mjungmath (previous) (diff)

comment:64 in reply to: ↑ 63 ; follow-up: Changed 3 months ago by tscrim

Replying to gh-mjungmath:

Replying to tscrim:

That is not a real error (the other patchbot also came back green).

It is a real error. Pyflakes is complaining about an unused import.

No, it is not. Pyflakes are not real errors (but they are good to clear).

You don't need to have all of the other branches merged in.

Why is that? Wouldn't that cause in a merge conflict later on?

There is no need to add unnecessary extra merge commits when there is no conflict or a need to have the latest branch. If later changes on the first branch are done to unrelated parts to the changes on this branch, there is no conflict.

Addendum - The pyflakes is handled on the later commits from #30074. So changing that here would cause a merge conflict and resolve a problem already fixed.

Last edited 3 months ago by tscrim (previous) (diff)

comment:65 in reply to: ↑ 64 ; follow-up: Changed 3 months ago by gh-mjungmath

Replying to tscrim:

Replying to gh-mjungmath:

Replying to tscrim:

That is not a real error (the other patchbot also came back green).

It is a real error. Pyflakes is complaining about an unused import.

No, it is not. Pyflakes are not real errors (but they are good to clear).

Okay, I see what you mean.

Replying to tscrim:

Addendum - The pyflakes is handled on the later commits from #30074. So changing that here would cause a merge conflict and resolve a problem already fixed.

Mh. Sorry, I don't get it. The #30074 will be merged into develop first because of the dependence, right? Afterwards, this ticket will be merged. If this ticket consists of an older version than #30074, this would cause a merge conflict, wouldn't it? Or, at least, the pyflakes error would be added again.

Last edited 3 months ago by gh-mjungmath (previous) (diff)

comment:66 in reply to: ↑ 65 ; follow-up: Changed 3 months ago by tscrim

Replying to gh-mjungmath:

Replying to tscrim:

Addendum - The pyflakes is handled on the later commits from #30074. So changing that here would cause a merge conflict and resolve a problem already fixed.

Mh. Sorry, I don't get it. The #30074 will be merged into develop first because of the dependence, right? Afterwards, this ticket will be merged. If this ticket consists of an older version than #30074, this would cause a merge conflict, wouldn't it? Or, at least, the pyflakes error would be added again.

No. The result will look like this:

B--0--X--T----M
    \        /
     0--0--H

where T is this branch and H is the head of #30074. So the #30074 branch will be merged first (all commits marked 0 starting from the latest beta B and H), then the commits from this branch (labelled X and T), which will result in the merge commit M. This would be exactly the same as if you merged in the latest commits from #30074 because that would be doing the same thing. The only difference is you have an extra merge commit in where this branch is pointing to.

comment:67 in reply to: ↑ 66 ; follow-up: Changed 3 months ago by gh-mjungmath

Replying to tscrim:

No. The result will look like this:

B--0--X--T----M
    \        /
     0--0--H

where T is this branch and H is the head of #30074. So the #30074 branch will be merged first (all commits marked 0 starting from the latest beta B and H), then the commits from this branch (labelled X and T), which will result in the merge commit M. This would be exactly the same as if you merged in the latest commits from #30074 because that would be doing the same thing. The only difference is you have an extra merge commit in where this branch is pointing to.

Ah, I see. This is because the commit in #30074 is more recent, right?

comment:68 in reply to: ↑ 67 ; follow-up: Changed 3 months ago by tscrim

Replying to gh-mjungmath:

Ah, I see. This is because the commit in #30074 is more recent, right?

Yes, but it is still based off a common base commit.

comment:69 in reply to: ↑ 68 Changed 3 months ago by gh-mjungmath

Replying to tscrim:

Replying to gh-mjungmath:

Ah, I see. This is because the commit in #30074 is more recent, right?

Yes, but it is still based off a common base commit.

Ah well. Thank you for explaining this, and sorry for delaying the process. This is good to know for future commits of mine.

comment:70 Changed 3 months ago by vbraun

  • Branch changed from public/manifolds/init_module_basis-30061 to 3ce0c15fc9b3d71842946ad6a0b0449771b1ada7
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.