Skip to content

Commit

Permalink
generalized Moore-Penrose inverses
Browse files Browse the repository at this point in the history
  • Loading branch information
@chakravala
chakravala committed Aug 14, 2020
1 parent 376795c commit 75ee157
Show file tree
Hide file tree
Showing 3 changed files with 95 additions and 22 deletions.
6 changes: 4 additions & 2 deletions src/algebra.jl
View file
Expand Up @@ -346,13 +346,15 @@ Exterior product as defined by the anti-symmetric quotient Λ≡⊗/~
@inline ∧(a::X,b::Y) where {X<:TensorAlgebra,Y<:TensorAlgebra} = interop(∧,a,b)
@inline ∧(a::TensorAlgebra{V},b::UniformScaling{T}) where {V,T<:Field} = a∧V(b)
@inline ∧(a::UniformScaling{T},b::TensorAlgebra{V}) where {V,T<:Field} = V(a)∧b
@generated ∧(t::T) where T<:SVector = Expr(:call,:∧,[:(t[$k]) for k ∈ 1:length(t)]...)
@generated ∧(t::T) where T<:SizedVector = Expr(:call,:∧,[:(t[$k]) for k ∈ 1:length(t)]...)
@generated ∧(t::T) where T<:SVector{N} where N = wedges([:(t[$i]) for i ∈ 1:N])
@generated ∧(t::T) where T<:SizedVector{N} where N = wedges([:(t[$i]) for i ∈ 1:N])
∧(::SVector{0,<:Chain{V}}) where V = one(V) # ∧() = 1
∧(::SizedVector{0,<:Chain{V}}) where V = one(V)
∧(t::Chain{V,1,<:Chain} where V) = ∧(value(t))
∧(a::X,b::Y,c::Z...) where {X<:TensorAlgebra,Y<:TensorAlgebra,Z<:TensorAlgebra} = ∧(a∧b,c...)

wedges(x,i=length(x)-1) = i ≠ 0 ? Expr(:call,:∧,wedges(x,i-1),x[1+i]) : x[1+i]

export ∧, ∨, ⊗

@pure function ∧(a::SubManifold{V},b::SubManifold{V}) where V
Expand Down
96 changes: 76 additions & 20 deletions src/composite.jl
View file
Expand Up @@ -342,13 +342,42 @@ function Cramer(N::Int,j=0)
return x,y,xy
end

@generated ∧(t::T) where T<:SVector{N} where N = wedges([:(t[$i]) for i ∈ 1:N])
@generated ∧(t::T) where T<:SizedVector{N} where N = wedges([:(t[$i]) for i ∈ 1:N])

@generated function Base.:\(t::SVector{M,<:Chain{V,1}},v::Chain{V,1}) where {M,V}
N = ndims(V)-1 # paste this into the REPL for faster eval
x,y,xy = Grassmann.Cramer(N)
W = M≠ndims(V) ? SubManifold(Manifold(M)) : V; N = M-1
if M == 1 && V == ℝ1
return :(Chain{V,1}(SVector(v[1]/t[1][1])))
elseif M == 2 && V == ℝ2
return quote
(a,A),(b,B),(c,C) = value(t[1]),value(t[2]),value(v)
x1 = (c-C*(b/B))/(a-A*(b/B))
return Chain{V,1}(x1,(C-A*x1)/B)
end
elseif M == 3 && V == ℝ3
return quote
dv = v/∧(t)[1]; c1,c2,c3 = value(t)
return Chain{V,1}(
(c2[2]*c3[3] - c3[2]*c2[3])*dv[1] +
(c3[1]*c2[3] - c2[1]*c3[3])*dv[2] +
(c2[1]*c3[2] - c3[1]*c2[2])*dv[3],
(c3[2]*c1[3] - c1[2]*c3[3])*dv[1] +
(c1[1]*c3[3] - c3[1]*c1[3])*dv[2] +
(c3[1]*c1[2] - c1[1]*c3[2])*dv[3],
(c1[2]*c2[3] - c2[2]*c1[3])*dv[1] +
(c2[1]*c1[3] - c1[1]*c2[3])*dv[2] +
(c1[1]*c2[2] - c2[1]*c1[2])*dv[3])
end
end
N<1 && (return :(inv(t)⋅v))
M > ndims(V) && (return :(tt=_transpose(t,$W); tt⋅(inv(Chain{$W,1}(t)⋅tt)⋅v)))
x,y,xy = Grassmann.Cramer(N) # paste this into the REPL for faster eval
mid = [:($(x[i])∧v∧$(y[end-i])) for i ∈ 1:N-1]
out = Expr(:call,:SVector,:(v∧$(y[end])),mid...,:($(x[end])∧v))
return Expr(:block,:((x1,y1)=(t[1],t[end])),xy...,
:(Chain{V,1}(getindex.($(Expr(:call,:./,out,:(t[1]∧$(y[end])))),1))))
detx = :(detx = (t[1]∧$(y[end])))
return Expr(:block,:((x1,y1)=(t[1],t[end])),xy...,detx,
:(Chain{$W,1}(column($(Expr(:call,:.⋅,out,:detx))./abs2(detx)))))
end

@generated function Base.in(v::Chain{V,1},t::SVector{N,<:Chain{V,1}}) where {V,N}
Expand All @@ -366,35 +395,62 @@ end
end
end

@generated function Base.inv(t::SVector{N,<:Chain{V,1}} where N) where V
N = ndims(V)-1
N<1 && (return :(_transpose(SVector(inv(t[1])))))
@generated function Base.inv(t::SVector{M,<:Chain{V,1}}) where {M,V}
W = M≠ndims(V) ? SubManifold(Manifold(M)) : V; N = M-1
N<1 && (return :(_transpose(SVector(inv(t[1])),$W)))
M > ndims(V) && (return :(tt = _transpose(t,$W); tt⋅inv(Chain{$W,1}(t)⋅tt)))
x,y,xy = Grassmann.Cramer(N)
out = if iseven(N)
val = if iseven(N)
Expr(:call,:SVector,y[end],[:($(y[end-i])∧$(x[i])) for i ∈ 1:N-1]...,x[end])
elseif M≠ndims(V)
Expr(:call,:SVector,y[end],[:($(iseven(i) ? :+ : :-)($(y[end-i])∧$(x[i]))) for i ∈ 1:N-1]...,:(-$(x[end])))
else
Expr(:call,:SVector,:(-$(y[end])),[:($(isodd(i) ? :+ : :-)($(y[end-i])∧$(x[i]))) for i ∈ 1:N-1]...,x[end])
end
return Expr(:block,:((x1,y1)=(t[1],t[end])),xy...,:(_transpose(.⋆($(Expr(:call,:./,out,:((t[1]∧$(y[end]))[1])))))))
out = if M≠ndims(V)
:(vector.($(Expr(:call,:./,val,:((t[1]∧$(y[end])))))))
else
:(.⋆($(Expr(:call,:./,val,:((t[1]∧$(y[end]))[1])))))
end
return Expr(:block,:((x1,y1)=(t[1],t[end])),xy...,:(_transpose($out,$W)))
end

@generated function grad(T::SVector{N,<:Chain{V,1}} where N,W=V) where V
N = ndims(V)-1
@generated function grad(T::SVector{M,<:Chain{V,1}}) where {M,V}
W = M≠ndims(V) ? SubManifold(Manifold(M)) : V; N = ndims(V)-1
M < ndims(V) && (return :(ct = Chain{$W,1}(T); map(↓(V),ct⋅inv(_transpose(T,$W)⋅ct))))
x,y,xy = Grassmann.Cramer(N)
out = if iseven(N)
val = if iseven(N)
Expr(:call,:SVector,[:($(y[end-i])∧$(x[i])) for i ∈ 1:N-1]...,x[end])
elseif M≠ndims(V)
Expr(:call,:SVector,y[end],[:($(iseven(i) ? :+ : :-)($(y[end-i])∧$(x[i]))) for i ∈ 1:N-1]...,:(-$(x[end])))
else
Expr(:call,:SVector,[:($(isodd(i) ? :+ : :-)($(y[end-i])∧$(x[i]))) for i ∈ 1:N-1]...,x[end])
end
return Expr(:block,:(t=_transpose(T,W)),:((x1,y1)=(t[1],t[end])),xy...,:(_transpose(.⋆($(Expr(:call,:./,out,:((t[1]∧$(y[end]))[1])))),↓(V))))
out = if M≠ndims(V)
:(vector.($(Expr(:call,:./,val,:((t[1]∧$(y[end])))))))
else
:(.⋆($(Expr(:call,:./,val,:((t[1]∧$(y[end]))[1])))))
end
return Expr(:block,:(t=_transpose(T,$W)),:((x1,y1)=(t[1],t[end])),xy...,:(_transpose($out,↓(V))))
end

@generated function Base.:\(t::SVector{N,<:Chain{M,1}},v::Chain{V,1}) where {N,M,V}
W = M≠ndims(V) ? SubManifold(Manifold(N)) : V
if ndims(M) > ndims(V)
:(ct=Chain{$W,1}(t); ct⋅(inv(_transpose(t,$W)⋅ct)⋅v))
else # ndims(M) < ndims(V) ? inv(tt⋅t)⋅(tt⋅v) : tt⋅(inv(t⋅tt)⋅v)
:(_transpose(t,$W)\v)
end
end
function inv_approx(t::Chain{M,1,<:Chain{V,1}}) where {M,V}
tt = transpose(t)
ndims(M) < ndims(V) ? (inv(tt⋅t))⋅tt : tt⋅inv(t⋅tt)
end

Base.:\(t::Chain{V,1,<:Chain{V,1}},v::Chain{V,1}) where V = value(t)\v
Base.:\(t::Chain{M,1,<:Chain{V,1}},v::Chain{V,1}) where {M,V} = ndims(M) < ndims(V) ? (tt=transpose(t); ((inv(tt⋅t))⋅tt)⋅v) : value(t)\v
Base.:\(t::Chain{V,1,<:Chain{M,1}},v::Chain{V,1}) where {M,V} = ndims(M) < ndims(V) ? ((inv(t⋅transpose(t)))⋅t)⋅v : value(transpose(t))\v
Base.:\(t::Chain{M,1,<:Chain{W,1}},v::Chain{V,1}) where {M,W,V} = value(t)\v
Base.in(v::Chain{V,1},t::Chain{W,1,<:Chain{V,1}}) where {V,W} = v ∈ value(t)
Base.inv(t::Chain{V,1,<:Chain{V,1}}) where V = inv(value(t))
grad(t::Chain{V,1,<:Chain{V,1}}) where V = grad(value(t))
Base.inv(t::Chain{V,1,<:Chain{W,1}}) where {W,V} = inv(value(t))
grad(t::Chain{V,1,<:Chain{W,1}}) where {V,W} = grad(value(t))
INV(m::Chain{V,1,<:Chain{V,1}}) where V = Chain{V,1,Chain{V,1}}(inv(SMatrix(m)))

export vandermonde
Expand All @@ -420,8 +476,8 @@ _vandermonde(x::Chain,V) = _vandermonde(value(x),V)

function vandermondeinterp(x,y,V,grid) # grid=384
coef = vandermonde(x,y,V) # Vandermonde ((inv(X'*X))*X')*y
xp=[minimum(x):(maximum(x)-minimum(x))/grid:maximum(x)...]
yp=zeros(length(xp)).+coef[1]
minx,maxx = minimum(x),maximum(x)
xp,yp = [minx:(maxx-minx)/grid:maxx...],coef[1]*ones(grid+1)
for d ∈ 1:length(coef)-1
yp += coef[d+1].*xp.^d
end # fill in polynomial terms
Expand Down
15 changes: 15 additions & 0 deletions src/multivectors.jl
View file
Expand Up @@ -36,6 +36,19 @@ Chain{V,G}(val::S) where {V,G,S<:AbstractVector{𝕂}} where 𝕂 = Chain{V,G,
ref = SVector{bg}([:(args[$i]) for i ∈ 1:bg])
:(Chain{V,G}($(Expr(:call,:(SVector{$bg,𝕂}),ref...))))
end

@generated function Chain{V}(args::𝕂...) where {V,𝕂}
bg = ndims(V); ref = SVector{bg}([:(args[$i]) for i ∈ 1:bg])
:(Chain{V,1}($(Expr(:call,:(SVector{$bg,𝕂}),ref...))))
end

@generated function Chain(args::𝕂...) where 𝕂
V = SubManifold(Manifold(length(args))); bg = ndims(V)
ref = SVector{bg}([:(args[$i]) for i ∈ 1:bg])
:(Chain{$V,1}($(Expr(:call,:(SVector{$bg,𝕂}),ref...))))
end


function Chain(val::𝕂,v::SubManifold{V,G}) where {V,G,𝕂}
N = ndims(V)
Chain{V,G}(setblade!(zeros(mvec(N,G,𝕂)),val,bits(v),Val{N}()))
Expand Down Expand Up @@ -76,8 +89,10 @@ Chain{V,1}(m::SMatrix{N,N}) where {V,N} = Chain{V,1}(Chain{V,1}.(getindex.(Ref(m
Chain{V,1,Chain{W,1}}(m::SMatrix{M,N}) where {V,W,M,N} = Chain{V,1}(Chain{W,1}.(getindex.(Ref(m),:,SVector{N}(1:N))))

transpose_row(t::SVector{N,<:Chain{V}},i,W=V) where {N,V} = Chain{W,1}(getindex.(t,i))
transpose_row(t::SizedVector{N,<:Chain{V}},i,W=V) where {N,V} = Chain{W,1}(getindex.(t,i))
transpose_row(t::Chain{V,1,<:Chain},i) where V = transpose_row(value(t),i,V)
@generated _transpose(t::SVector{N,<:Chain{V,1}},W=V) where {N,V} = :(Chain{V,1}(transpose_row.(Ref(t),$(SVector{ndims(V)}(1:ndims(V))),W)))
@generated _transpose(t::SizedVector{N,<:Chain{V,1}},W=V) where {N,V} = :(Chain{V,1}(transpose_row.(Ref(t),$(SVector{ndims(V)}(1:ndims(V))),W)))
Base.transpose(t::Chain{V,1,<:Chain{V,1}}) where V = _transpose(value(t))
Base.transpose(t::Chain{V,1,<:Chain{W,1}}) where {V,W} = _transpose(value(t),V)

Expand Down

0 comments on commit 75ee157

Please sign in to comment.