chakravala / grassmann.jl Goto Github PK

I notice that certain constructors for Chain work when called generically, but not when called with the specific type:

julia> Chain{V,1}(1,2)
1v₁ + 2v₂

julia> Chain{V,1,Int}(1,2)
ERROR: MethodError: no method matching Chain{⟨++⟩,1,Int64,253} where 253(::Int64, ::Int64)

Similarly:

julia> Chain(B.v1)
1v₁ + 0v₂

julia> Chain{Int}(B.v1)
1v₁ + 0v₂

but:

julia> Chain{V}(B.v1)
ERROR: MethodError: no method matching SubManifold{⟨++⟩,0,0x0000000000000000}(::Simplex{⟨++⟩,0,v,Int64})

julia> Chain{V,1}(B.v1)
ERROR: DimensionMismatch("No precise constructor for StaticArrays.SArray{Tuple{2},SubManifold{⟨++⟩,1,0x0000000000000001},1,2} found. Length of input was 1.")

julia> Chain{V,1,Int}(B.v1)
ERROR: MethodError: Cannot `convert` an object of type SubManifold{⟨++⟩,1,0x0000000000000001} to an object of type StaticArrays.SArray{Tuple{2},Int64,1,2}

It would be convenient if one could always call a specific constructor. That is, if a constructor creates a certain type Chain{V,1,T}, then I'd prefer if I was always able to call Chain{V,1,T}, and not have to use Chain{V,1} or Chain{T} instead.

This is not urgent.

Which definition of the inner product between multivectors do you use? I'm asking because I was confused by the fact that v1 ⋅ v12 is giving me 0 as a result, and I would have expected some multiple of v1 (which v12 ⋅ v1 gives me).

For property testing I like to use random vectors (simplices, chains, multivectors). It would be convenient to have a method for rand() that does that, called e.g. as rand(Chain{V,1,Float64}).

This is not urgent.

I encountered this problem:

julia> using Grassmann

julia> D = 2
2

julia> S = Signature(D)
⟨++⟩

julia> V = SubManifold(S)
⟨++⟩

julia> x = Chain{V,1}(1,2)
1v₁ + 2v₂

julia> xs = [x]
1-element Array{Chain{⟨++⟩,1,Int64,2},1}:
 1v₁ + 2v₂

julia> @show xs
ERROR: _show_nonempty(::IO, ::AbstractVector, ::String) is not implemented
Stacktrace:
 [1] error(::String) at ./error.jl:33
 [2] _show_nonempty(::IOContext{Base.GenericIOBuffer{Array{UInt8,1}}}, ::Array{Chain{⟨++⟩,1,Int64,2},1}, ::String) at ./arrayshow.jl:420
 [3] show(::Base.GenericIOBuffer{Array{UInt8,1}}, ::Array{Chain{⟨++⟩,1,Int64,2},1}) at ./arrayshow.jl:437
 [4] sprint(::Function, ::Array{Chain{⟨++⟩,1,Int64,2},1}; context::Nothing, sizehint::Int64) at ./strings/io.jl:105
 [5] #repr#353 at ./strings/io.jl:227 [inlined]
 [6] repr(::Array{Chain{⟨++⟩,1,Int64,2},1}) at ./strings/io.jl:227
 [7] top-level scope at show.jl:634

The problem seems to be that Grassmann defines (see multivectors.jl, line 139) the method

@pure Base.ndims(::Vector{Chain{V,G,T,X}} where {G,T,X}) where V = ndims(V)

which confuses the base library when showing vectors.

I think it's wrong to overload ndims for any Vector{T}. All Vector{T} are one-dimensional arrays. The corresponding definition for Base.parent probably also shouldn't be there.

This breaks my code since using @show is important for debugging for me.

This fails:

julia> @basis 0
(⟨⟩, ⟨⟩)

julia> Chain{V,1}(v)
ERROR: DimensionMismatch("No precise constructor for SArray{Tuple{0},SubManifold{⟨⟩,0,0x0000000000000000},1,0} found. Length of input was 1.")
Stacktrace:
 [1] SArray{Tuple{0},SubManifold{⟨⟩,0,0x0000000000000000},1,0}(::Tuple{Tuple{Tuple{Tuple{SubManifold{⟨⟩,0,0x0000000000000000}}}}}) at /Users/eschnett/.julia/packages/StaticArrays/MB5hl/src/convert.jl:1
 [2] StaticArray at /Users/eschnett/.julia/packages/StaticArrays/MB5hl/src/convert.jl:4 [inlined] (repeats 4 times)
 [3] Chain{⟨⟩,1,𝕂,253} where 253 where 𝕂(::SubManifold{⟨⟩,0,0x0000000000000000}) at /Users/eschnett/.julia/packages/Grassmann/Of0BW/src/multivectors.jl:33
 [4] top-level scope at REPL[209]:1
 [5] eval(::Module, ::Any) at ./boot.jl:331
 [6] eval_user_input(::Any, ::REPL.REPLBackend) at /Users/eschnett/src/julia-1.4/usr/share/julia/stdlib/v1.4/REPL/src/REPL.jl:86
 [7] run_backend(::REPL.REPLBackend) at /Users/eschnett/.julia/packages/Revise/0meWR/src/Revise.jl:1075
 [8] top-level scope at none:0

This is of low priority to me.

The following is a quote of @enkimute from the discussion in #3

check the way he defined the complement and figure out if this may provide a method that works for degenerate metrics in the same unified fashion..

Opening this issue is meant to track ongoing discussions regarding degenerate metrics and complements

Please see comments in the pdf of your JuliaCon submission.
grassmann-juliacon-2019.pdf
Summary of Comments on grassmann-juliacon-2019.pdf

This paper is a description of Grassmann.jl and its supporting packages, which are appear quite impressive and extensive. Unfortunately, it
assumes extensive knowledge of terminology and notation from the reader, without ever defining what that reader should know, whose
notations they should be familiar with, and what that reader can expect to gain from reading.

As someone who has written a lot of geometric and topological code in low dimensions, and been interested in Clifford algebra enough
to read the book of Dorst, Fontijne, & Mann and some of Hestenes, I was eager to review this. I was disappointed because it does not
clearly say what the package does, nor what decisions were made between alternative implementations. (There are many intriguing hints,
e.g., Due to the design of the VectorBundle type
system, these operations make code optimization easier at compile-time by evaluating the bit parameters.)

As Fred Brooks has said, "Since anyone can create something new [in a synthetic field], that alone does not establish a contribution.
Rather, one must show that the creation is better." (Quoted from https://cra.org/resources/best-practice-memos/evaluating-computerscientists-and-engineers-for-promotion-and-tenure/) A paper should do more than quote from the documentation of a package; it
should provide context so the reader (from the author's desired audience) can understand the decisions made by the package author and
why they are better than alternatives.

I wanted to do a test your package by implementing Quaternions and tried to create a quaternion q for a rotation around a vector v = (a, b, c) by an angle alpha using exponentials:
q = exp(alpha/2*(ai + bj + ck))

However the exponential function doesn't seem to be fully implemented yet for multivectors? At least i get an error in the following MWE:

julia> using Grassmann

julia> i, j, k = hyperplanes(ℝ^3)
3-element Array{SValue{⟨+++⟩,2,B,Int64} where B,1}:
 -1v₂₃
 1v₁₃ 
 -1v₁₂

julia> alpha = 0.5π
1.5707963267948966

julia> exp(alpha/2*(i))
0.7071067811865476 - 0.7071067811865475v₂₃

julia> a, b, c = 1/sqrt(2) * [1, 1, 0]
3-element Array{Float64,1}:
 0.7071067811865475
 0.7071067811865475
 0.0               

julia> exp(alpha/2*(a*i + b*j + c*k))
ERROR: MethodError: no method matching abs(::StaticArrays.MArray{Tuple{8},Float64,1,8})
Closest candidates are:
  abs(::Bool) at bool.jl:91
  abs(::Float16) at float.jl:520
  abs(::Float32) at float.jl:521
  ...
Stacktrace:
 [1] abs(::MultiVector{Float64,⟨+++⟩,8}) at .julia/packages/AbstractTensors/fUnNY/src/AbstractTensors.jl:55
 [2] _broadcast_getindex_evalf at ./broadcast.jl:578 [inlined]
 [3] _broadcast_getindex at ./broadcast.jl:551 [inlined]
 [4] getindex at ./broadcast.jl:511 [inlined]
 [5] copyto_nonleaf!(::Array{Float64,1}, ::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1},Tuple{Base.OneTo{Int64}},typeof(abs),Tuple{Base.Broadcast.Extruded{Array{TensorAlgebra{⟨+++⟩},1},Tuple{Bool},Tuple{Int64}}}}, ::Base.OneTo{Int64}, ::Int64, ::Int64) at ./broadcast.jl:928
 [6] copy at ./broadcast.jl:791 [inlined]
 [7] materialize at ./broadcast.jl:753 [inlined]
 [8] exp(::SBlade{Float64,⟨+++⟩,2,3}) at .julia/packages/Grassmann/wWv7E/src/algebra.jl:961
 [9] top-level scope at none:0

Applying abs2 to a Chain loses type information

julia> S = S"∞∅++"
⟨∞∅++⟩
julia> V = SubManifold(S)
⟨∞∅++⟩
julia> B = Λ(S)
DirectSum.Basis{⟨∞∅++⟩,16}(v, v∞, v∅, v₁, v₂, v∞∅, v∞₁, v∞₂, v∅₁, v∅₂, v₁₂, v∞∅₁, v∞∅₂, v∞₁₂, v∅₁₂, v∞∅₁₂)
julia> x = Chain(B.v1)
0v∞ + 0v∅ + 1v₁ + 0v₂
julia> x |> typeof
Chain{⟨∞∅++⟩,1,Int64,4}
julia> abs2(x) |> typeof
MultiVector{⟨∞∅++⟩,Any,16}

While x has element type Int, abs2(x) has only Any.

When I add two multivectors A and B everything works like expected:

julia> using Grassmann

julia> @basis "2"
(⟨++⟩, v, v₁, v₂, v₁₂)

julia> A = 2v1 + v2
2v₁ + 1v₂

julia> B = v1 + v2
1v₁ + 1v₂

julia> A + B
3v₁ + 2v₂

julia> A
2v₁ + 1v₂

julia> B
1v₁ + 1v₂

However when I add one of the basis vectors (v1 or v2) that are not orthogonal to the multivector, it changes the multivector:

julia> v1 + A
3v₁ + 1v₂

julia> A
3v₁ + 1v₂

julia> v1 + A
4v₁ + 1v₂

julia> v1 + A
5v₁ + 1v₂

In geometric algebra at least the 0-vector wedged with another 0-vector should be non-zero if both scalars are non-zero.

julia> 2v∧3v
0v

The docstring for >>> seems to document \oslash. I therefore originally assumed these are just different names for the same function, but that doesn't seem to be so.

I'm trying to implement sqrt for pure scalar multivectors (so for Basis{V, 0, 0}, MValue{V, 0, v, T} and SValue{V, 0, v, T}), but I have a hard time to find the right syntax for the methods type declaration. It is easy for Basis:

sqrt(t::Basis{B, 0, 0x0000000000000000}) where {B} = 1

But I'm not sure how to define the type of MValue and SValue for scalar basis v. I can get v into the scope and write

julia> basis"2"
(⟨++⟩, v, v₁, v₂, v₁₂)

julia> sqrt(x::SValue{B, 0, v, T}) where {B, T} = sqrt(x.v)
sqrt (generic function with 23 methods)

julia> sqrt(13*v)
3.605551275463989

but that seems ugly and does not make sure that v matches the scalar basis for any B.
I guess I have to use the function indexbasis but I'm not sure about its syntax.

So this is more like a question than an actual issue (except maybe a documentation issue).

Btw: What is the difference between SValue and MValue?

PS: I hope that its clear, but I'm saying it anyway: I really like this package, thanks for all the work you are putting into it.

julia> using Grassmann

julia> @basis tangent(ℝ^0, 2, 2)
(T²⟨₁₂⟩, v, ∂₁, ∂₂, ∂₁₂)

julia> [1∂1, 2∂2]
Bus error: 10

This is with

julia> versioninfo()
Julia Version 1.4.2-pre.0
Commit ef4fe83698* (2020-04-15 16:24 UTC)
Platform Info:
  OS: macOS (x86_64-apple-darwin19.3.0)
  CPU: Intel(R) Core(TM) i7-8850H CPU @ 2.60GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-8.0.1 (ORCJIT, skylake)

and

(GrassmannPlayground) pkg> status
Project GrassmannPlayground v0.1.0
Status `~/src/jl/GrassmannPlayground/Project.toml`
  [c3fe647b] AbstractAlgebra v0.9.1
  [537997a7] AbstractPlotting v0.10.1
  [159f3aea] Cairo v1.0.3
  [a81c6b42] Compose v0.8.2
  [22fd7b30] DirectSum v0.5.6
  [8d0d7f98] GaloisFields v1.0.1
  [5c1252a2] GeometryBasics v0.2.4
  [a2cc645c] GraphPlot v0.4.2
  [4df31cd9] Grassmann v0.5.7
  [edad4870] Leibniz v0.0.5
  [093fc24a] LightGraphs v1.3.1
  [ee78f7c6] Makie v0.10.0
  [93e0c654] Reduce v1.2.5

Hi @chakravala,

In your JuliaCon paper, I noticed the statement "Mixed-symmetry algebra with Leibniz.jl and Grassmann.jl, having the geometric algebraic product chain rule, yields automatic differ- entiation and...", and I found a related simple example in the README.md file, section "Differential forms and Leibniz tangent algebra".

Would you mind providing one or two concrete examples on symbolic differentiation using Grassmann.jl? For example, with respect to a rotor/versor, or a bivector.
I'm thinking that the key part is about a proper setup of a @mixedbasis tangent .... Am I right?

I'm reading Chapter 8 of the book "Geometric Algebra for Computer Science", but there are just a few examples of basic functions. I would like to apply vector/multivector differentiation on some more functions, and verify my results.
I'm trying to differentiate a function like \sum_i weight_i || V Bivector_0 /V \wedge v_i ||^2 with respect to the versor V.

Thank you!

For some reason, I'm getting the following behavior for exponentials:

exp(0 * (v1 + v2)) # 1v⃖ (as expected)
exp(0 * v1) # 1.0 - NaNv₁ (unexpected)

The problem appears to be specific to zero multiples of basis elements.

I am testing Grassmann.jl with conformal projective geometry, but the geometric product seems to be broken for those spaces.

Easiest way to reproduce is with your example for the null-basis conformal split in README.md. That gives me:

julia> using Grassmann; @basis S"∞∅++"
[ Info: Recompiling stale cache file /home/karl/.julia/compiled/v1.1/Grassmann/i5fGo.ji for Grassmann [4df31cd9-4c27-5bea-88d0-e6a7146666d8]
(⟨∞∅++⟩, v, v∞, v∅, v₁, v₂, v∞∅, v∞₁, v∞₂, v∅₁, v∅₂, v₁₂, v∞∅₁, v∞∅₂, v∞₁₂, v∅₁₂, v∞∅₁₂)

julia> v∞^2, v∅^2, v1^2, v2^2
ERROR: MethodError: no method matching conformal(::Signature{4,3,0x0000000000000002,0}, ::UInt64, ::UInt64)
Closest candidates are:
  conformal(::Any, ::Any, ::W<:(Signature{n,m,s,Diff} where Diff)) where {n, m, s, W<:(Signature{n,m,s,Diff} where Diff)} at .julia/packages/Grassmann/wWv7E/src/parity.jl:164
  conformal(::Any, ::Any, ::W<:(DiagonalForm{n,m,s,Diff} where Diff)) where {n, m, s, W<:(DiagonalForm{n,m,s,Diff} where Diff)} at .julia/packages/Grassmann/wWv7E/src/parity.jl:164
Stacktrace:
 [1] *(::Basis{⟨∞∅++⟩,1,0x0000000000000004}, ::Basis{⟨∞∅++⟩,1,0x0000000000000004}) at .julia/packages/Grassmann/wWv7E/src/algebra.jl:148
 [2] ^(::Basis{⟨∞∅++⟩,1,0x0000000000000004}, ::Int64) at .julia/packages/Grassmann/wWv7E/src/algebra.jl:919
 [3] macro expansion at ./none:0 [inlined]
 [4] literal_pow(::typeof(^), ::Basis{⟨∞∅++⟩,1,0x0000000000000004}, ::Val{2}) at ./none:0
 [5] top-level scope at none:0

It worked in 0.1.4, but doesn't in 0.1.5

http://lomont.org/math/geometric-algebra/Universal%20Geometric%20Algebra%20-%20Hestenes%20-%201988.pdf is the current url for Hestenes' Universal Geometric Algebra paper

The readme contains this example

using Grassmann, Makie; @basis S"∞+++"
streamplot(vectorfield(exp((π/4)*(v12+v∞3)),V(2,3,4),V(1,2,3)),-1.5..1.5,-1.5..1.5,-1.5..1.5,gridsize=(10,10))

With just the statements above, vectorfield is undefined. If I also ] add GeometryTypes and using GeometryTypes, then vectorfield is found, but I'm then stuck at the next error:

julia> streamplot(vectorfield(exp((π/4)*(v12+v∞3)),V(2,3,4),V(1,2,3)),-1.5..1.5,-1.5..1.5,-1.5..1.5,gridsize=(10,10))
ERROR: MethodError: no method matching ptype(::GeometryBasics.Point{3,Float64})
Closest candidates are:
  ptype(::GeometryTypes.Point{N,T} where N) where T at /Users/eschnett/.julia/packages/Grassmann/9buUZ/src/Grassmann.jl:307
Stacktrace:
 [1] (::Grassmann.var"#707#709"{MultiVector{⟨∞+++⟩,Float64,16},SubManifold{⟨∞+++⟩,3,0x000000000000000e},SubManifold{⟨∞+++⟩,3,0x0000000000000007}})(::GeometryBasics.Point{3,Float64}) at /Users/eschnett/.julia/packages/Grassmann/9buUZ/src/Grassmann.jl:310
 [2] (::AbstractPlotting.var"#apply_f#505"{Grassmann.var"#707#709"{MultiVector{⟨∞+++⟩,Float64,16},SubManifold{⟨∞+++⟩,3,0x000000000000000e},SubManifold{⟨∞+++⟩,3,0x0000000000000007}}})(::GeometryBasics.Point{3,Float64}, ::Type{T} where T) at /Users/eschnett/.julia/packages/AbstractPlotting/7mERO/src/basic_recipes/basic_recipes.jl:1061
 [3] streamplot_impl(::Type{T} where T, ::Grassmann.var"#707#709"{MultiVector{⟨∞+++⟩,Float64,16},SubManifold{⟨∞+++⟩,3,0x000000000000000e},SubManifold{⟨∞+++⟩,3,0x0000000000000007}}, ::GeometryBasics.HyperRectangle{3,Float32}, ::Tuple{Int64,Int64}, ::Float64, ::Int64, ::Float64) at /Users/eschnett/.julia/packages/AbstractPlotting/7mERO/src/basic_recipes/basic_recipes.jl:1081
 [4] (::AbstractPlotting.var"#508#512")(::Function, ::GeometryBasics.HyperRectangle{3,Float32}, ::Tuple{Int64,Int64}, ::Float64, ::Int64, ::Float64) at /Users/eschnett/.julia/packages/AbstractPlotting/7mERO/src/basic_recipes/basic_recipes.jl:1138
 [5] lift(::Function, ::Observables.Observable{Grassmann.var"#707#709"{MultiVector{⟨∞+++⟩,Float64,16},SubManifold{⟨∞+++⟩,3,0x000000000000000e},SubManifold{⟨∞+++⟩,3,0x0000000000000007}}}, ::Observables.Observable{GeometryBasics.HyperRectangle{3,Float32}}, ::Vararg{Any,N} where N) at /Users/eschnett/.julia/packages/AbstractPlotting/7mERO/src/interaction/nodes.jl:8
 [6] plot!(::StreamPlot{...}) at /Users/eschnett/.julia/packages/AbstractPlotting/7mERO/src/basic_recipes/basic_recipes.jl:1132
 [7] plot!(::Scene, ::Type{StreamPlot{...}}, ::Attributes, ::Tuple{Observables.Observable{Grassmann.var"#707#709"{MultiVector{⟨∞+++⟩,Float64,16},SubManifold{⟨∞+++⟩,3,0x000000000000000e},SubManifold{⟨∞+++⟩,3,0x0000000000000007}}},Observables.Observable{IntervalSets.Interval{:closed,:closed,Float64}},Observables.Observable{IntervalSets.Interval{:closed,:closed,Float64}},Observables.Observable{IntervalSets.Interval{:closed,:closed,Float64}}}, ::Observables.Observable{Tuple{Grassmann.var"#707#709"{MultiVector{⟨∞+++⟩,Float64,16},SubManifold{⟨∞+++⟩,3,0x000000000000000e},SubManifold{⟨∞+++⟩,3,0x0000000000000007}},GeometryBasics.HyperRectangle{3,Float32}}}) at /Users/eschnett/.julia/packages/AbstractPlotting/7mERO/src/interfaces.jl:633
 [8] plot!(::Scene, ::Type{StreamPlot{...}}, ::Attributes, ::Function, ::Vararg{Any,N} where N; kw_attributes::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at /Users/eschnett/.julia/packages/AbstractPlotting/7mERO/src/interfaces.jl:563
 [9] plot!(::Scene, ::Type{StreamPlot{...}}, ::Attributes, ::Function, ::IntervalSets.Interval{:closed,:closed,Float64}, ::IntervalSets.Interval{:closed,:closed,Float64}, ::Vararg{IntervalSets.Interval{:closed,:closed,Float64},N} where N) at /Users/eschnett/.julia/packages/AbstractPlotting/7mERO/src/interfaces.jl:532
 [10] streamplot(::Function, ::Vararg{Any,N} where N; attributes::Base.Iterators.Pairs{Symbol,Tuple{Int64,Int64},Tuple{Symbol},NamedTuple{(:gridsize,),Tuple{Tuple{Int64,Int64}}}}) at /Users/eschnett/.julia/packages/AbstractPlotting/7mERO/src/recipes.jl:15
 [11] top-level scope at REPL[8]:1
 [12] eval(::Module, ::Any) at ./boot.jl:331
 [13] eval_user_input(::Any, ::REPL.REPLBackend) at /Users/eschnett/src/julia-1.4/usr/share/julia/stdlib/v1.4/REPL/src/REPL.jl:86
 [14] run_backend(::REPL.REPLBackend) at /Users/eschnett/.julia/packages/Revise/C272c/src/Revise.jl:1075
 [15] top-level scope at none:0

The inverse of a vector a works like expected for vectors with positive square a² > 0 in the way that baa⁻¹ = b:

julia> using Grassmann

julia> basis"-+++"
(⟨-+++⟩, v, v₁, v₂, v₃, v₄, v₁₂, v₁₃, v₁₄, v₂₃, v₂₄, v₃₄, v₁₂₃, v₁₂₄, v₁₃₄, v₂₃₄, v₁₂₃₄)

julia> v2^2
v

julia> v1*v2*inv(v2), v1
(v₁, v₁)

It however does not give the correct result, I think, for vectors with negative square a² < 0:

julia> v1^2
-1v

julia> v2*v1*inv(v1), v2
(-1v₂, v₂)

Thanks for making such a great GA lib available.

How do you visualize the results?
I'm trying to analyze geometric structure in GA.
Can you recommend a workflow or any libraries (if you are already aware of by change) ?

Thanks!
Best,
Hongying

P.S. Previously I use "gable" by Leo Dorst, and also noticed GAlgebra (in Python) and its Julia wrapper).

The equality check between Multivectors and Numbers currently does not work correctly for numbers that are not zero:

julia> using Grassmann

julia> basis"2"
(⟨++⟩, v, v₁, v₂, v₁₂)

julia> a = v + v1 - v1
1v⃖

julia> typeof(a)
MultiVector{Int64,⟨++⟩,4}

julia> a == 1
false

julia> b = a - 1
0v⃖

julia> b == 0
true

julia> a - 1 == 0
true

A simple fix would probably be changing lines 309 and 310 in multivectors.jl to:

==(a::Number,b::MultiVector{S,V,G} where {S,V}) where G = prod(0 .== value(b - a))
==(a::MultiVector{S,V,G} where {S,V},b::Number) where G = prod(0 .== value(a - b))

This issue is to track and discuss the new updates related to differential SubManifold algebra.

chakravala/DirectSum.jl@e76bea7
chakravala/DirectSum.jl@1009e80
chakravala/DirectSum.jl@60c0e27
chakravala/AbstractTensors.jl@b7d2ff6

This is related to the finalization of the JuliaCon submission from JuliaCon/proceedings-review#18

It is sometimes convenient to call a multiplication operator with a single argument, as in \wedge(xs...) where xs is a collection of chains. For example, * returns its argument when called with a single argument. It would be convenient if \wedge did the same; this would avoid code such as

            if length(xs) == 1
                vol = abs(xs[1])
            else
                vol = abs(∧(xs...))
            end

Wanted to add an example for barycentric coordinates to your tutorial section ( as suggested on bivector.net). Ran into a problem.

julia> basis"3"

julia> a = 0v1
0v₁

julia> b = 1v1
1v₁

julia> c = 1v2
1v₂

julia> ab = b-a
1v₁

julia> ca = a-c
0v₁ - 1v₂ + 0v₃

julia> bc = c-b
-1v₁ + 1v₂ + 0v₃

julia> A = -ab∧ca
1v₁₂ + 0v₁₃ + 0v₂₃

julia> bc∧ca
1v₁₂ + 0v₁₃ + 0v₂₃

julia> bc∧ca/A
ERROR: DimensionMismatch("No precise constructor for SArray{Tuple{8},Int64,1,8} found. Length of input was 3.")
Stacktrace:
[1] SArray{Tuple{8},Int64,1,8}(::Tuple{Tuple{Tuple{Tuple{Int64,Int64,Int64}}}}) at /Users/mewert/.julia/packages/StaticArrays/1g9bq/src/convert.jl:1
[2] StaticArray at /Users/mewert/.julia/packages/StaticArrays/1g9bq/src/convert.jl:4 [inlined] (repeats 3 times)
[3] macro expansion at /Users/mewert/.julia/packages/StaticArrays/1g9bq/src/convert.jl:4 [inlined]
[4] reverse(::Chain{⟨+++⟩,2,Int64,3}) at /Users/mewert/.julia/packages/Grassmann/OeTVc/src/products.jl:531
[5] ~ at /Users/mewert/.julia/packages/DirectSum/Ofduh/src/generic.jl:242 [inlined]
[6] inv(::Chain{⟨+++⟩,2,Int64,3}) at /Users/mewert/.julia/packages/Grassmann/OeTVc/src/algebra.jl:671
[7] /(::Chain{⟨+++⟩,2,Int64,3}, ::Chain{⟨+++⟩,2,Int64,3}) at /Users/mewert/.julia/packages/AbstractTensors/cHS7X/src/AbstractTensors.jl:156
[8] top-level scope at REPL[31]:1

Seems to be the zero elements

julia> A
1v₁₂ + 0v₁₃ + 0v₂₃

julia> (bc∧ca)/1v12
1.0v⃖

Am I wrong or is it a bug ?

using Reduce
using Grassmann
@basis S"∞∅+++"
p = :x*v1 + :y*v2 + :z*v3
P = v∅ + p + 0.5 * p^2 * v∞

ERROR: LoadError: MethodError: no method matching +(::Expr, ::Float64)
Closest candidates are:
+(::Any, ::Any, !Matched::Any, !Matched::Any...) at operators.jl:529
+(!Matched::Bool, ::T<:AbstractFloat) where T<:AbstractFloat at bool.jl:104
+(!Matched::Float64, ::Float64) at float.jl:395
...
Stacktrace:
[1] macro expansion at .julia/packages/StaticArrays/3KEjZ/src/mapreduce.jl:41 [inlined]
[2] _map at .julia/packages/StaticArrays/3KEjZ/src/mapreduce.jl:21 [inlined]
[3] map at .julia/packages/StaticArrays/3KEjZ/src/mapreduce.jl:14 [inlined]
[4] +(::Array{Any,1}, ::StaticArrays.SizedArray{Tuple{5},Any,1,1}) at /home/gacquewi/.julia/packages/StaticArrays/3KEjZ/src/linalg.jl:11
[5] +(::MChain{Any,⟨∞∅+++⟩,1,5}, ::MultiVector{Any,⟨∞∅+++⟩,32}) at .julia/packages/Grassmann/cEoGT/src/algebra.jl:1310
[6] +(::Basis{⟨∞∅+++⟩,1,0x0000000000000002}, ::MChain{Any,⟨∞∅+++⟩,1,5}, ::MultiVector{Any,⟨∞∅+++⟩,32}) at ./operators.jl:529
[7] top-level scope at wgCGAbug.jl:5
[8] include at ./boot.jl:328 [inlined]
[9] include_relative(::Module, ::String) at ./loading.jl:1094
[10] include(::Module, ::String) at ./Base.jl:31
[11] exec_options(::Base.JLOptions) at ./client.jl:295
[12] _start() at ./client.jl:464
in expression starting at wgCGAbug.jl:5

May be helpful to add AbstractPlotting.inline!(true) to Makie plot examples for newcomers so they're not confused (like I was!). MakieOrg/Makie.jl#392

I'm trying to calculate the inverse of a blade with Grassmann v0.1.3:

julia> basis"3"
(⟨+++⟩, v, v₁, v₂, v₃, v₁₂, v₁₃, v₂₃, v₁₂₃)
julia> inv(v1+v2)
ERROR: MethodError: no method matching inv(::MBlade{Int64,⟨+++⟩,1,3})

Is this function missing, or am I looking in the wrong place?

This fails:

julia> using Grassmann
julia> @basis tangent(ℝ^0, 1, 1)
(T¹⟨₁⟩, v, ∂₁)
julia> ∂1, ∂1 ⊗ ∂1
(∂₁, 0v)

julia> ∂1 ⊗ ∂1 ⊗ ∂1
ERROR: StackOverflowError:
Stacktrace:
 [1] interop(::typeof(⊗), ::Simplex{T¹⟨₁⟩,0,v,Int64}, ::SubManifold{T¹⟨₁⟩,1,0x0000000000000001}) at /Users/eschnett/.julia/packages/AbstractTensors/zchWS/src/AbstractTensors.jl:115
 [2] ⊗(::Simplex{T¹⟨₁⟩,0,v,Int64}, ::SubManifold{T¹⟨₁⟩,1,0x0000000000000001}) at /Users/eschnett/.julia/packages/AbstractTensors/zchWS/src/AbstractTensors.jl:148
 ... (the last 2 lines are repeated 39990 more times)
 [79983] interop(::typeof(⊗), ::Simplex{T¹⟨₁⟩,0,v,Int64}, ::SubManifold{T¹⟨₁⟩,1,0x0000000000000001}) at /Users/eschnett/.julia/packages/AbstractTensors/zchWS/src/AbstractTensors.jl:115

It seems that, if ∂1^n is zero, then ∂1^(n+1) fails with a stack overflow.

According to the documentation of \oslash, this should hold:

help?> ⊘
"⊘" can be typed by \oslash<tab>

search: ⊘

  ⊘(ω::TensorAlgebra,η::TensorAlgebra)

  Sandwich product: ω⊘η = (~η)⊖ω⊖η

julia> v⊘v1 == (~v1)⊖v⊖v1
false

The source code for this library has a very large number of functions marked @pure. My understanding is that the use of this macro is emphatically discouraged. In this issue JuliaLang/julia#20704, Keno said

@pure means I am @vtjnash and I certify that this function doesn't break any of the implicit assumptions I use in the compiler ;)

In this thread: https://discourse.julialang.org/t/pure-macro/3871/3, Matt says that

Improper @pure annotations can introduce bugs. The optimizations it enables rely on an extremely strict definition of pure. It really should be named something like @hyperpure . Some of the restrictions include:

• It must always return exactly ( === ) the same result for a given input. Watch out for mutable types. I think constant globals are okay, though.
• The function it’s used on cannot be further extended by other methods after it gets called.
• It cannot recurse.
• It’s undocumented and not exported (for good reason), but this means the complete list of preconditions is really only in a few people’s heads.

For instance,

Furthermore, adding methods to @pure functions is said to not be kosher and this repository is replete with @pure functions with methods added on.

In this thread https://discourse.julialang.org/t/can-pure-functions-throw-an-error/18459, Jameson advises that @pure functions throwing errors is "Probably, a very very bad idea", but I see that there are functions that throw in multivectors.jl and Grassmann.jl.

Given this, I think cutting back on most if not all uses of the uses of @pure in this package is probably a good idea.

Perhaps this deserves a separate issue, but I noticed that in utilities.jl there are struct definitions

struct Grade{G}
    @pure Grade{G}() where G = new{G}()
end

## Dimension{N}

struct Dimension{N}
    @pure Dimension{N}() where N = new{N}()
end

I believe the above inner constructors are no-ops that don't need to exist. For instance:

julia> struct MyStruct1{T} end

julia> struct MyStruct2{T}
           MyStruct2{T}() where {T} = new{T}()
       end

julia> @code_lowered MyStruct1{1}()
CodeInfo(
1 ─ %1 = (Core.apply_type)(Main.MyStruct1, $(Expr(:static_parameter, 1)))
│   %2 = %new(%1)
└──      return %2
)

julia> @code_lowered MyStruct2{1}()
CodeInfo(
1 ─ %1 = (Core.apply_type)(Main.MyStruct2, $(Expr(:static_parameter, 1)))
│   %2 = %new(%1)
└──      return %2
)

I see that you've made it possible to use REDUCE for symbolic manipulations, but I was hoping to be able to use sympy. I tried the obvious thing:

using SymPy, Grassmann
@vars x y  # This defines x and y as SymPy.Sym objects
@basis ℝ^3
x*v + y*v12

But the last line causes an error:

ERROR: setindex!() with non-isbitstype eltype is not supported by StaticArrays. Consider using SizedArray.
Stacktrace:
 [1] error(::String) at ./error.jl:33
 [2] setindex! at /Users/me/.julia/packages/StaticArrays/mcf7t/src/MArray.jl:134 [inlined]
 [3] setmulti! at /Users/me/.julia/dev/Grassmann/src/algebra.jl:31 [inlined]
 [4] +(::SValue{⟨+++⟩,0,v,Sym}, ::SValue{⟨+++⟩,2,v₁₂,Sym}) at /Users/boyle/.julia/dev/Grassmann/src/algebra.jl:796
 [5] top-level scope at none:0

I tried various other combinations, and it seems that whenever the sympy scalars multiply different Grassmann basis elements, this error occurs. But even if they multiply the same basis element, something weird happens:

julia> (x*v1) + (y*v1)
x + yv₁

This sentence from the README

By default, the coefficients are required to be <:Number.

suggests that is should work right out of the box, because

julia> SymPy.Sym <: Number
true

I also tried

Grassmann.generate_product_algebra(SymPy.Sym)

and

Grassmann.generate_product_algebra(SymPy.Sym,:(SymPy.:*),:(SymPy.:+),:(SymPy.:-))

to no avail.

I'm sorry that my julia skills are insufficient to let me figure this out, but this is my first day with this language (being attracted mostly by this exciting package). I would appreciate any help you can give.

This is fundamentally tied to the Hodge star dual operation

https://en.wikipedia.org/wiki/Hodge_star_operator

Help or discussion is welcome on this topic

I am experimenting with mixedbasis to represent non-trivial metric tensors. (Is there a better way to represent them?) I receive these errors:

julia> using Grassmann
julia> S = S"-+++"
julia> @mixedbasis S

julia> v1w1(v1)
ERROR: MethodError: no method matching evaluate2(::SubManifold{v₁₂₃₄w¹²³⁴,2,0x0000000000000011}, ::SubManifold{⟨-++++---⟩*,1,0x0000000000000001})
Closest candidates are:
  evaluate2(::A, ::B) where {V, A<:TensorTerm{V,1}, B<:TensorTerm{V,1}} at /Users/eschnett/.julia/packages/DirectSum/Ng7R9/src/operations.jl:229
Stacktrace:
 [1] (::SubManifold{v₁₂₃₄w¹²³⁴,2,0x0000000000000011})(::SubManifold{⟨-++++---⟩*,1,0x0000000000000001}) at /Users/eschnett/.julia/packages/DirectSum/Ng7R9/src/operations.jl:249
 [2] interform at /Users/eschnett/.julia/packages/AbstractTensors/zchWS/src/AbstractTensors.jl:124 [inlined]
 [3] (::SubManifold{⟨-++++---⟩*,2,0x0000000000000011})(::SubManifold{⟨-+++⟩,1,0x0000000000000001}) at /Users/eschnett/.julia/packages/DirectSum/Ng7R9/src/operations.jl:255
 [4] top-level scope at REPL[10]:1
 [5] eval(::Module, ::Any) at ./boot.jl:331
 [6] eval_user_input(::Any, ::REPL.REPLBackend) at /Users/eschnett/src/julia-1.4/usr/share/julia/stdlib/v1.4/REPL/src/REPL.jl:86
 [7] run_backend(::REPL.REPLBackend) at /Users/eschnett/.julia/packages/Revise/C272c/src/Revise.jl:1075
 [8] top-level scope at none:0

It seems that operations are not defined on pure basis vectors, only on linear combinations of them. However:

julia> (-v1w1+v2w2)(v1)
ERROR: InexactError: Bool(-1)
Stacktrace:
 [1] Bool at ./float.jl:73 [inlined]
 [2] convert at ./number.jl:7 [inlined]
 [3] convert at ./essentials.jl:310 [inlined] (repeats 2 times)
 [4] setindex! at ./array.jl:826 [inlined]
 [5] dualform(::SubManifold{⟨-++++---⟩*,8,0x00000000000000ff}) at /Users/eschnett/.julia/dev/Grassmann/src/forms.jl:142
 [6] (::Chain{⟨-++++---⟩*,2,Int64,28})(::SubManifold{⟨-++++---⟩*,1,0x0000000000000001}) at /Users/eschnett/.julia/dev/Grassmann/src/forms.jl:187
 [7] interform at /Users/eschnett/.julia/packages/AbstractTensors/zchWS/src/AbstractTensors.jl:124 [inlined]
 [8] (::Chain{⟨-++++---⟩*,2,Int64,28})(::SubManifold{⟨-+++⟩,1,0x0000000000000001}) at /Users/eschnett/.julia/dev/Grassmann/src/forms.jl:176
 [9] top-level scope at REPL[11]:1
 [10] eval(::Module, ::Any) at ./boot.jl:331
 [11] eval_user_input(::Any, ::REPL.REPLBackend) at /Users/eschnett/src/julia-1.4/usr/share/julia/stdlib/v1.4/REPL/src/REPL.jl:86
 [12] run_backend(::REPL.REPLBackend) at /Users/eschnett/.julia/packages/Revise/C272c/src/Revise.jl:1075
 [13] top-level scope at none:0

The latter seems to happen in line 142 of forms.jl:

            @inbounds mV[Q] = (intlog(Y)+1,V[intlog(x)+1])

Should the second tuple element have an != 0 added to handle negative signatures?

I am using projective or conformal bases, such as

julia> S = S"∞++"
julia> B = Λ(S)

I want to write dimension-generic code, and want to access the i-th basis vector. Unfortunately, Grassmann counts starting from v∞, so that v₁ is B.v(2) etc. It would be more convenient if B.v(1) returned v1, independent of whether v∞ or v∅ were present. I currently have to write B.v(i+1) or B.v(i+2), depending on the basis.

I realize this would be a breaking change, and obviously this is not an urgent feature request.

julia> basis"+++"
(⟨+++⟩, v, v₁, v₂, v₃, v₁₂, v₁₃, v₂₃, v₁₂₃)

julia> v1⨼(v1*v2)
-1v₂

julia> v1*(v1*v2)
v₂

Hi,

• On notation:
  What is the difference between every pair of the three: @basis S"++++-", @basis V"++++-" and @basis "++++-"?
  I'm a little confused in trying the examples in README.md. It seems they all give me the same results.

• On algebra
  In Grassmann.js, is V"++++-" also the conformal model G^{4,1} of R^3?
  (I noticed in the section "Null-basis of the conformal split" that we can declare a space with signature ∞ and ∅, but did not figure out how I can construct a signature such that it is in correspondence with "++++-".

Thank you,
Hongying

I get an error for the geometric product of v∅ and -v∞:

julia> using Grassmann

julia> @basis S"∞∅+"
(⟨∞∅+⟩, v, v∞, v∅, v₁, v∞∅, v∞₁, v∅₁, v∞∅₁)

julia> v∅*v∞
-1 - 1v∞∅

julia> v∅*(-v∞)
ERROR: MethodError: no method matching SValue{⟨∞∅+⟩,G,B,T} where T where B where G(::Int64, ::MultiVector{Int64,⟨∞∅+⟩,8})
Closest candidates are:
  SValue{⟨∞∅+⟩,G,B,T} where T where B where G(::T) where {V, T} at .julia/packages/Grassmann/H9Zog/src/multivectors.jl:97
  SValue{⟨∞∅+⟩,G,B,T} where T where B where G(::Any, ::SValue{V,G,B,T} where T where B) where {V, G} at .julia/packages/Grassmann/H9Zog/src/multivectors.jl:90
  SValue{⟨∞∅+⟩,G,B,T} where T where B where G(::Any, ::MValue{V,G,B,T} where T where B) where {V, G} at .julia/packages/Grassmann/H9Zog/src/multivectors.jl:91
  ...
Stacktrace:
 [1] *(::Basis{⟨∞∅+⟩,1,0x0000000000000002}, ::SValue{⟨∞∅+⟩,1,v∞,Int64}) at .julia/packages/Grassmann/H9Zog/src/algebra.jl:157
 [2] top-level scope at none:0

(v1.3) pkg> add Grassmann
Updating registry at ~/.julia/registries/General
Updating git-repo https://github.com/JuliaRegistries/General.git
Resolving package versions...
Installed ReplMaker ─────── v0.2.3
Installed Requires ──────── v1.0.1
Installed AbstractTensors ─ v0.4.3
Installed DirectSum ─────── v0.5.3
Installed Reduce ────────── v1.2.4
Installed Grassmann ─────── v0.5.0
Updating ~/.julia/environments/v1.3/Project.toml
[4df31cd9] + Grassmann v0.5.0
Updating ~/.julia/environments/v1.3/Manifest.toml
[398f06c4] + AbstractLattices v0.1.2
[a8e43f4a] + AbstractTensors v0.4.3
[861a8166] + Combinatorics v1.0.0
[459fdd68] + ComputedFieldTypes v0.1.0
[22fd7b30] + DirectSum v0.5.3
[9dda63f9] + ForceImport v0.0.3
[4df31cd9] + Grassmann v0.5.0
[edad4870] + Leibniz v0.0.3
[f27b6e38] + Polynomials v0.6.0
[3cdcf5f2] + RecipesBase v0.7.0
[93e0c654] + Reduce v1.2.4
[b873ce64] + ReplMaker v0.2.3
[ae029012] + Requires v1.0.1
[90137ffa] + StaticArrays v0.12.1
[a4af3ec5] + SyntaxTree v1.0.1
[2a0f44e3] + Base64
[8ba89e20] + Distributed
[b77e0a4c] + InteractiveUtils
[8f399da3] + Libdl
[37e2e46d] + LinearAlgebra
[56ddb016] + Logging
[d6f4376e] + Markdown
[3fa0cd96] + REPL
[9a3f8284] + Random
[ea8e919c] + SHA
[9e88b42a] + Serialization
[6462fe0b] + Sockets
[2f01184e] + SparseArrays
[10745b16] + Statistics
[8dfed614] + Test
[cf7118a7] + UUIDs
[4ec0a83e] + Unicode
Building Reduce → ~/.julia/packages/Reduce/7af5f/deps/build.log

julia> using Grassmann
[ Info: Precompiling Grassmann [4df31cd9-4c27-5bea-88d0-e6a7146666d8]
ERROR: LoadError: UndefVarError: VectorBundle not defined
Stacktrace:
[1] top-level scope at /Users/mewert/.julia/packages/Leibniz/XlDhG/src/Leibniz.jl:53
[2] include at ./boot.jl:328 [inlined]
[3] include_relative(::Module, ::String) at ./loading.jl:1105
[4] include(::Module, ::String) at ./Base.jl:31
[5] top-level scope at none:2
[6] eval at ./boot.jl:330 [inlined]
[7] eval(::Expr) at ./client.jl:425
[8] top-level scope at ./none:3
in expression starting at /Users/mewert/.julia/packages/Leibniz/XlDhG/src/Leibniz.jl:53
ERROR: LoadError: Failed to precompile Leibniz [edad4870-8a01-11e9-2d75-8f02e448fc59] to /Users/mewert/.julia/compiled/v1.3/Leibniz/wNkI4_RK5IG.ji.
Stacktrace:
[1] error(::String) at ./error.jl:33
[2] compilecache(::Base.PkgId, ::String) at ./loading.jl:1283
[3] _require(::Base.PkgId) at ./loading.jl:1024
[4] require(::Base.PkgId) at ./loading.jl:922
[5] require(::Module, ::Symbol) at ./loading.jl:917
[6] include at ./boot.jl:328 [inlined]
[7] include_relative(::Module, ::String) at ./loading.jl:1105
[8] include(::Module, ::String) at ./Base.jl:31
[9] top-level scope at none:2
[10] eval at ./boot.jl:330 [inlined]
[11] eval(::Expr) at ./client.jl:425
[12] top-level scope at ./none:3
in expression starting at /Users/mewert/.julia/packages/Grassmann/5163y/src/Grassmann.jl:46
ERROR: Failed to precompile Grassmann [4df31cd9-4c27-5bea-88d0-e6a7146666d8] to /Users/mewert/.julia/compiled/v1.3/Grassmann/i5fGo_RK5IG.ji.
Stacktrace:
[1] error(::String) at ./error.jl:33
[2] compilecache(::Base.PkgId, ::String) at ./loading.jl:1283
[3] _require(::Base.PkgId) at ./loading.jl:1024
[4] require(::Base.PkgId) at ./loading.jl:922
[5] require(::Module, ::Symbol) at ./loading.jl:917

julia>

The documentation for \oslash reads:

Sandwich product: ω⊘η = (~ω)⊖η⊖ω

However, it seems to implement η⊘ω instead (with arguments commuted):

julia> R = exp(π/4*v12)
0.7071067811865476 + 0.7071067811865475v₁₂

julia> (~R) ⊖ v1 ⊖ R
0.0 + 2.220446049250313e-16v₁ + 1.0v₂

julia> R ⊘ v1
-0.7071067811865476 + 0.7071067811865475v₁₂

julia> v1 ⊘ R
0.0 + 2.220446049250313e-16v₁ + 1.0v₂

Hi.
Thanks for helping to push geometric algebra into the mainstream with your hard work.
I am reading introductory material on GA and using Clifford (Python) to add, multiply, reflect, etc.
For the most part, I bypassed linear algebra and just took up with GA because it looks interesting.
I would like to use Julia (and Grassmann) to learn GA, but I don't understand most of the README, yet.
In particular, I couldn't tell how to simply add, multiply or reflect a vector in R3.
I was able to use Clifford because I could follow their docs, but the Grassmann README seems geared to people who already know what they are doing...
Thanks, again.

The result type of abs2 seems to depend on argument values:

julia> @basis 2
(⟨++⟩, v, v₁, v₂, v₁₂)

julia> abs2(MultiVector(v+v1)) |> typeof
MultiVector{⟨++⟩,Int64,4}

julia> abs2(MultiVector(v)) |> typeof
Simplex{⟨++⟩,0,v,Int64}

The return type of a function should usually only depend on the argument types.

I am currently not running benchmarks, so this does currently not affect me.

works on a Multivector and Basis, but doesn't seem to give correct result for Blades.

julia> a
1.0v₁ + 2.0v₂ + 0.0v₃

julia> b
3.0v₁

julia> ab = a*b
3.0 - 6.0v₁₂

julia> reverse(ab)
3.0 + 6.0v₁₂

julia> reverse(ab(2))
6.0v₁ + 0.0v₂ + 0.0v₃

julia> typeof(ab(2))
SBlade{Float64,⟨+++⟩,2,3}

I find that the uparrow and downarrow operators for conformal geometries convert rational numbers to floating-point numbers:

julia> using Grassmann

julia> S = Signature(4, 1, 1)
⟨∞∅++⟩

julia> V = SubManifold(S)
⟨∞∅++⟩

julia> B = Λ(S)
DirectSum.Basis{⟨∞∅++⟩,16}(v, v∞, v∅, v₁, v₂, v∞∅, v∞₁, v∞₂, v∅₁, v∅₂, v₁₂, v∞∅₁, v∞∅₂, v∞₁₂, v∅₁₂, v∞∅₁₂)

julia> Chain{V,1}(1//1,2//1,3//1,4//1)
(1//1)v∞ + (2//1)v∅ + (3//1)v₁ + (4//1)v₂

julia> x = Chain{V,1}(1//1,2//1,3//1,4//1)
(1//1)v∞ + (2//1)v∅ + (3//1)v₁ + (4//1)v₂

julia> ↑(x)
-0.0 + 11.5v∞ + 3.0v∅ + 3.0v₁ + 4.0v₂

julia> ↓(x)
0.0 + 1.5v₁ + 2.0v₂

The problem is caused by the terms G.v∞/2 and inv(G.v∞∅), which don't have any type information and thus default to Float64. Multiplying the basis vectors by one(valuetype(\omega)) before the division/inversion solves the problem.

This breaks my code because I like to use rational numbers for testing to avoid round-off errors. I am currently using a corrected copy of these functions that I copied into my code.

I find this:

julia> using Grassmann
julia> V=SubManifold(Signature(1))

julia> Chain{V,1}(SVector(BigInt(10))) |> typeof
Chain{⟨+⟩,1,BigInt,1}

julia> +Chain{V,1}(SVector(BigInt(10))) |> typeof
Chain{⟨+⟩,1,Any,1}

Applying a unary + or - to a chain converts the type to Any.

I encounter this error:

julia> basis"+++"
julia> (v1+v2) + (v1+v2)*(v1+v2)
ERROR: DimensionMismatch("No precise constructor for StaticArrays.MArray{Tuple{3},Int64,1,3} found. Length of input was 8.")
Stacktrace:
 [1] StaticArrays.MArray{Tuple{3},Int64,1,3}(::Tuple{Tuple{Tuple{NTuple{8,Int64}}}}) at /Users/eschnett/.julia/packages/StaticArrays/VyRz3/src/convert.jl:1
 [2] Type at /Users/eschnett/.julia/packages/StaticArrays/VyRz3/src/convert.jl:4 [inlined] (repeats 3 times)
 [3] convert at /Users/eschnett/.julia/packages/StaticArrays/VyRz3/src/convert.jl:10 [inlined]
 [4] value at /Users/eschnett/.julia/packages/Grassmann/LeGo1/src/multivectors.jl:327 [inlined]
 [5] +(::MBlade{Int64,⟨+++⟩,1,3}, ::MultiVector{Int64,⟨+++⟩,8}) at /Users/eschnett/.julia/packages/Grassmann/LeGo1/src/algebra.jl:1045
 [6] top-level scope at none:0

There seems to be a bug in the constructor when creating a MultiVector from a blade:

julia> using Grassmann

julia> basis"++"
(⟨++⟩, v, v₁, v₂, v₁₂)

julia> a = v1 + v2
1v₁ + 1v₂

julia> typeof(a)
MBlade{Int64,⟨++⟩,1,2}

julia> MultiVector(a)
ERROR: type MArray has no field v
Stacktrace:
 [1] getproperty(::Any, ::Symbol) at ./sysimg.jl:18
 [2] MultiVector(::MBlade{Int64,⟨++⟩,1,2}) at .julia/packages/Grassmann/H9Zog/src/multivectors.jl:261
 [3] top-level scope at none:0