rjrosati / symbolictensors.jl Goto Github PK

Manipulate tensors symbolically in Julia! Currently needs a SymPy dependency, but work is ongoing to change the backend to SymbolicUtils.jl

License: MIT License

Julia 100.00%

julia physics symbolic-computation tensor-algebra symbolic-math symbolic-manipulation

symbolictensors.jl's Introduction

SymbolicTensors

Many numerical tensor manipulation packages exist (e.g. Einsum.jl), but treating tensors at a purely numeric level throws away a lot of potential optimizations. Often, it's possible to exploit the symmetries of a problem to dramatically reduce the calculation steps necessary, or perform some tensor contractions symbolically rather than numerically.

SymbolicTensors.jl is designed to exploit these simplifications to perform symbolic calculations and generate more efficient input into numeric tensor packages than you would write by hand. It based on SymPy.jl, sympy.tensor.tensor, xTensor, and ITensors.jl.

See the talk about this package given at JuliaCon 2020: https://www.youtube.com/watch?v=_b4JIv044GY

Example calculations

using SymbolicTensors
using SymPy

spacetime = TensorIndexType("spacetime","f")
@indices spacetime μ ν σ ρ η
x = TensorHead("x",[spacetime])
δ = spacetime.delta
# one way to write the metric on a sphere
g = 4*δ(-μ,-ν)/(1+x(μ)*x(ν)*δ(-μ,-ν))^2

# compute the christoffel symbols
Γ = (diff(g(-μ,-ν),x(σ)) - diff(g(-ν,-σ),x(μ)) + diff(g(-σ,-μ),x(ν)))/2
Γ = factor(contract_metric(canon_bp(Γ),spacetime.metric))

At this point Γ is a symbolic tensor expression, written in Einstein notation. It could represent 2 dimensional or 10,000 dimensional tensors in the same form -- this is one of the great advantages of working at the symbolic tensor level.

We could manipulate Γ further, or convert it into an array.

# convert Γ to Array{Sym}
xarr = symbols("x y",real=true)
garr = replace_with_arrays(g,Dict(x(ρ) => xarr, spacetime.delta(-ρ,-η) => [1 0; 0 1]))
Γarr = replace_with_arrays(Γ,Dict(x(ρ) => xarr, spacetime.delta(-ρ,-η) => [1 0; 0 1], spacetime => garr))

Now that we have an array in Γarr, Quote allows us to convert it into a native Julia function.

# Then Quote it and eval to a Julia function
quot = Quote("Christoffel",Γarr,[ x for x in xarr])
Christ = eval(quot)
Christ(0,1)

Rough TODO

Nicer errors
testing coverage
add conversion to ITensors or equivalent
documentation
symbolic derivatives?

known bugs

check issues

symbolictensors.jl's People

Contributors

Stargazers

Watchers

Forkers

miguelraz coligne potatoasad yinxx

symbolictensors.jl's Issues

Cannot make @heads work

Cannot figure out quickly what's wrong with that

lorenz = TensorIndexType("lorenz","L",4)
@heads [lorenz] pX p3 p1 p2

the first argument is not recognized as AbstractArray{SymbolicTensors.TensorIndexType}. Instead, Julia think that it is Expr...(?):

MethodError: no method matching var"@heads"(::LineNumberNode, ::Module, ::Expr, ::Symbol, ::Symbol, ::Symbol, ::Symbol)

Closest candidates are:

var"@heads"(::LineNumberNode, ::Module, !Matched::AbstractArray{SymbolicTensors.TensorIndexType}, ::Any...) at \.julia\packages\SymbolicTensors\DPYJA\src\tensor.jl:244

#macroexpand#[email protected]:115[inlined]
macroexpand@expr.jl:114[inlined]
try_macroexpand(::Module, ::Base.UUID, ::Expr)@PlutoRunner.jl:248
var"#run_expression#25"(::Bool, ::typeof(Main.PlutoRunner.run_expression), ::Module, ::Expr, ::Base.UUID, ::Nothing, ::Nothing)@PlutoRunner.jl:477
top-level scope@none:1

the @head macro looks ok, but different to @indices.

Indexing a tensor should also create a scalar

Right now, B(i,-i) isn't treated as a scalar.

Kronecker delta behavior

At the moment, the upstream contract_metric will transform KD(-i,-j) metric(j,k) into KD(-i,k).
This isn't what the Kronecker delta means to me (i.e. always the identity matrix). Is this actually standard somewhere?

Okay, even avoiding this issue, any expression with Kronecker deltas, after enough manipulation, will be a mess of several deltas contracted with each other. canon_bp and contract_metric don't seem to simplify this at all. Should we have a new function, contract_delta or equivalent that reduces these expressions?

For example, here is the Ricci scalar I calculated from the metric 6*α*delta(-i,-j) / (1-F(i)*F(j)*delta(-i,-j))^2:

(-24)*α*(-1)*(ts1 - 1)**(-3)*KD(f_0, -f_0)*KD(f_1, -f_1) + 24*α*(-1)*(ts1 - 1)**(-3)*KD(f_0, f_1)*KD(-f_0, -f_1) + (-864)*α**2*(ts1 - 1)**(-6)*F(f_0)*F(f_1)*KD(-f_0, f_2)*KD(-f_1, f_3)*KD(-f_2, f_4)*KD(-f_3, -f_4) + (-144)*α*(ts1 - 1)**2*1/(ts1**6 - 6*ts1**5 + 15*ts1**4 - 20*ts1**3 + 15*ts1**2 - 6*ts1 + 1)*F(f_0)*F(f_1)*KD(-f_0, f_2)*KD(-f_1, -f_2)*KD(f_3, -f_3) + (-144)*α**2*(ts1 - 1)**(-6)*F(f_0)*F(f_1)*KD(-f_0, f_2)*KD(-f_1, -f_2)*KD(f_3, -f_3)*KD(f_4, -f_4) + 144*α*(ts1 - 1)**2*1/(ts1**6 - 6*ts1**5 + 15*ts1**4 - 20*ts1**3 + 15*ts1**2 - 6*ts1 + 1)*F(f_0)*F(f_1)*KD(-f_0, f_2)*KD(-f_1, f_3)*KD(-f_2, -f_3) + 432*α**2*(ts1 - 1)**(-6)*F(f_0)*F(f_1)*KD(-f_0, f_2)*KD(-f_1, -f_2)*KD(f_3, f_4)*KD(-f_3, -f_4) + 576*α**2*(ts1 - 1)**(-6)*F(f_0)*F(f_1)*KD(-f_0, f_2)*KD(-f_1, f_3)*KD(-f_2, -f_3)*KD(f_4, -f_4))

The correct answer should only depend on α and the TensorIndexType dimension.

Calling epsilon throws `invalid name` error

Just started with the packages (thanks for that), and quickly stumbled upon several issues.
Let me drop it here. Maybe I can cover the issues to tests once it is resolved

lorenz = TensorIndexType("lorenz","L",4)
@indices lorenz μ ν σ ρ
ϵ = lorenz.epsilon
ϵ(μ,ν,σ,ρ)

gives the error:

PyError ($(Expr(:escape, :(ccall(#= C:\Users\Mikhasenko\.julia\packages\PyCall\L0fLP\src\pyfncall.jl:43 =# @pysym(:PyObject_Call), PyPtr, (PyPtr, PyPtr, PyPtr), o, pyargsptr, kw))))) <class 'ValueError'>

ValueError('invalid name')
File "\.julia\conda\3\lib\site-packages\sympy\core\basic.py", line 953, in subs
rv = rv._subs(old, new, **kwargs)
File "\.julia\conda\3\lib\site-packages\sympy\core\cache.py", line 72, in wrapper
retval = cfunc(*args, **kwargs)
File "\.julia\conda\3\lib\site-packages\sympy\core\basic.py", line 1067, in _subs
rv = fallback(self, old, new)
File "C:\Users\Mikhasenko\.julia\conda\3\lib\site-packages\sympy\core\basic.py", line 1044, in fallback
rv = self.func(*args)
File "\.julia\conda\3\lib\site-packages\sympy\tensor\tensor.py", line 1735, in __new__
raise ValueError("invalid name")

pyerr_check@exception.jl:62[inlined]
pyerr_check@exception.jl:66[inlined]
_handle_error(::String)@exception.jl:83
macro expansion@exception.jl:97[inlined]
#[email protected]:43[inlined]
disable_sigint@c.jl:458[inlined]
__pycall!@pyfncall.jl:42[inlined]
_pycall!(::PyCall.PyObject, ::PyCall.PyObject, ::Tuple{SymPy.Sym, SymbolicTensors.TensorIndex}, ::Int64, ::Ptr{Nothing})@pyfncall.jl:29
_pycall!(::PyCall.PyObject, ::PyCall.PyObject, ::Tuple{SymPy.Sym, SymbolicTensors.TensorIndex}, ::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})@pyfncall.jl:11
var"#_#114"(::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, ::PyCall.PyObject, ::SymPy.Sym, ::Vararg{Any})@pyfncall.jl:86
(::PyCall.PyObject)(::SymPy.Sym, ::Vararg{Any})@pyfncall.jl:86
(::SymPy.Sym)(::SymbolicTensors.TensorIndex, ::Vararg{SymbolicTensors.TensorIndex})@utils.jl:54
top-level scope@Local: 1[inlined]

self-contracted tensors aren't properly treated as scalars

In practice, these are probably avoidable unless the user inputs them or calls contract_metric directly.
The resulting objects A(f_0,-f_0) are currently of type Sym and aren't handled correctly in diff or canon_bp.

Two equations* do not render in Pluto

What is about align?

Switch backend to SymbolicUtils?

SymbolicUtils.jl is pure Julia, and is extendable: https://github.com/JuliaSymbolics/SymbolicUtils.jl

It could replace the sympy backend.

`replace_with_arrays` is too picky about up/down indices

replace_with_arrays has sufficient information to work in cases like:

replace_with_arrays(A(i)*A(j)*field.metric(-i,-j), Dict(A(-i) => [0, 1], field.metric(i,j) => [1 0; 0 1]) )

This function should internally know that field.metric(-i,-j) is just inv([1 0; 0 1]), and that A(i) is just A(-j)*field.metric(i,j) == [0,1]' * [1 0; 0 1], but it gives an unspecified metric error.

This error is upstream and also exists in sympy.tensor.tensor. Should we rewrite replace_with_arrays?