Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Lowering connection cfs from raising connection cfs #112

Merged
merged 10 commits into from
Jun 24, 2024
Merged

Lowering connection cfs from raising connection cfs #112

Show file tree
Hide file tree
Changes from 3 commits
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
da33c81
lowering works
TSGut Jun 21, 2024
04109cc
fix bugs, add QR lowering too
TSGut Jun 21, 2024
d77c898
change convention, implement lowering
TSGut Jun 21, 2024
dc70b6b
X^1 = X
TSGut Jun 23, 2024
5279d29
Update src/SemiclassicalOrthogonalPolynomials.jl
TSGut Jun 24, 2024
d0b95dd
Update src/SemiclassicalOrthogonalPolynomials.jl
TSGut Jun 24, 2024
24d7dd1
consistency change
TSGut Jun 24, 2024
6bafbd9
bugfix
TSGut Jun 24, 2024
e3be3fb
Update ci.yml
dlfivefifty Jun 24, 2024
f4e9683
increase test coverage
TSGut Jun 24, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
145 changes: 104 additions & 41 deletions src/SemiclassicalOrthogonalPolynomials.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -180,19 +180,27 @@ function semiclassical_jacobimatrix(Q::SemiclassicalJacobi, a, b, c)
cholesky_jacobimatrix(I-Q.X,Q)
elseif iszero(Δa) && iszero(Δb) && isone(Δc)
cholesky_jacobimatrix(Q.t*I-Q.X,Q)
elseif isone(-Δa) && iszero(Δb) && iszero(Δc) # in these cases we currently have to reconstruct
# TODO: This is re-constructing. It should instead use reverse Cholesky (or an alternative)!
semiclassical_jacobimatrix(Q.t,a,b,c)
elseif iszero(Δa) && isone(-Δb) && iszero(Δc)
# TODO: This is re-constructing. It should instead use reverse Cholesky (or an alternative)!
semiclassical_jacobimatrix(Q.t,a,b,c)
elseif iszero(Δa) && iszero(Δb) && isone(-Δc)
# TODO: This is re-constructing. It should instead use reverse Cholesky (or an alternative)!
semiclassical_jacobimatrix(Q.t,a,b,c)
elseif a > Q.a # iterative raising by 1
semiclassical_jacobimatrix(SemiclassicalJacobi(Q.t, Q.a+1, Q.b, Q.c, Q), a, b, c)
elseif b > Q.b
semiclassical_jacobimatrix(SemiclassicalJacobi(Q.t, Q.a, Q.b+1, Q.c, Q), a, b, c)
elseif c > Q.c
semiclassical_jacobimatrix(SemiclassicalJacobi(Q.t, Q.a, Q.b, Q.c+1, Q), a, b, c)
# TODO: Implement lowering via QL/reverse Cholesky or via inverting R?
# elseif b < Q.b # iterative lowering by 1
# semiclassical_jacobimatrix(SemiclassicalJacobi(Q.t, Q.a, Q.b-1, Q.c, Q), a, b, c)
# elseif c < Q.c
# semiclassical_jacobimatrix(SemiclassicalJacobi(Q.t, Q.a, Q.b, Q.c-1, Q), a, b, c)
# elseif a < Q.a
# semiclassical_jacobimatrix(SemiclassicalJacobi(Q.t, Q.a-1, Q.b, Q.c, Q), a, b, c)
elseif a < Q.a # iterative lowering by 1
semiclassical_jacobimatrix(SemiclassicalJacobi(Q.t, Q.a-1, Q.b, Q.c, Q), a, b, c)
elseif b < Q.b
semiclassical_jacobimatrix(SemiclassicalJacobi(Q.t, Q.a, Q.b-1, Q.c, Q), a, b, c)
elseif c < Q.c
semiclassical_jacobimatrix(SemiclassicalJacobi(Q.t, Q.a, Q.b, Q.c-1, Q), a, b, c)
else
error("Not Implemented")
end
Expand All @@ -206,7 +214,6 @@ LanczosPolynomial(P::SemiclassicalJacobi{T}) where T =

gives either a mapped `Jacobi` or `LanczosPolynomial` version of `P`.
"""
# TODO: Use ConvertedOPs for integer special cases?
toclassical(P::SemiclassicalJacobi{T}) where T = iszero(P.c) ? Normalized(jacobi(P.b, P.a, UnitInterval{T}())) : LanczosPolynomial(P)

copy(P::SemiclassicalJacobi) = P
Expand All @@ -231,6 +238,7 @@ jacobimatrix(P::SemiclassicalJacobi) = P.X

"""
op_lowering(Q, y)

Gives the Lowering operator from OPs w.r.t. (x-y)*w(x) to Q
as constructed from Chistoffel–Darboux
"""
Expand Down Expand Up @@ -260,56 +268,110 @@ function semijacobi_ldiv(P::SemiclassicalJacobi, Q)
(P \ R) * _p0(R̃) * (R̃ \ Q)
end

# returns conversion operator from SemiclassicalJacobi P to SemiclassicalJacobi Q in a single step via decomposition.
semijacobi_ldiv_direct(Q::SemiclassicalJacobi, P::SemiclassicalJacobi) = semijacobi_ldiv_direct(Q.t, Q.a, Q.b, Q.c, P)
function semijacobi_ldiv_direct(Qt, Qa, Qb, Qc, P::SemiclassicalJacobi)
(Qt ≈ P.t) && (Qa ≈ P.a) && (Qb ≈ P.b) && (Qc ≈ P.c) && return SymTridiagonal(Ones(∞),Zeros(∞))
Δa = Qa-P.a
Δb = Qb-P.b
Δc = Qc-P.c
M = Diagonal(Ones(∞))
if iseven(Δa) && iseven(Δb) && iseven(Δc)
M = qr(P.X^(Δa÷2)*(I-P.X)^(Δb÷2)*(Qt*I-P.X)^(Δc÷2)).R
"""
semijacobi_ldiv_direct(Q::SemiclassicalJacobi, P::SemiclassicalJacobi)

Returns conversion operator from SemiclassicalJacobi `P` to SemiclassicalJacobi `Q` in a single step via decomposition.
Numerically unstable if the parameter modification is large. Typically one should instead use `P \\ Q` which is equivalent to `semijacobi_ldiv(P,Q)` and proceeds step by step.
"""
function semijacobi_ldiv_direct(Q::SemiclassicalJacobi, P::SemiclassicalJacobi)
(Q.t ≈ P.t) && (Q.a ≈ P.a) && (Q.b ≈ P.b) && (Q.c ≈ P.c) && return SymTridiagonal(Ones(∞),Zeros(∞))
Δa = Q.a-P.a
Δb = Q.b-P.b
Δc = Q.c-P.c
# special case (Δa,Δb,Δc) = (2,0,0)
if isone(Δa/2) && iszero(Δb) && iszero(Δc)
TSGut marked this conversation as resolved.
Show resolved Hide resolved
M = qr(P.X^(Δa÷2)).R
TSGut marked this conversation as resolved.
Show resolved Hide resolved
return ApplyArray(*, Diagonal(sign.(view(M,band(0))).*Fill(abs.(1/M[1]),∞)), M) # match normalization choice P_0(x) = 1
TSGut marked this conversation as resolved.
Show resolved Hide resolved
# special case (Δa,Δb,Δc) = (0,2,0)
elseif iszero(Δa) && isone(Δb/2) && iszero(Δc)
TSGut marked this conversation as resolved.
Show resolved Hide resolved
M = qr((I-P.X)^(Δb÷2)).R
return ApplyArray(*, Diagonal(sign.(view(M,band(0))).*Fill(abs.(1/M[1]),∞)), M)
# special case (Δa,Δb,Δc) = (0,0,2)
elseif iszero(Δa) && iszero(Δb) && isone(Δc/2)
TSGut marked this conversation as resolved.
Show resolved Hide resolved
M = qr((Q.t*I-P.X)^(Δc÷2)).R
return ApplyArray(*, Diagonal(sign.(view(M,band(0))).*Fill(abs.(1/M[1]),∞)), M)
elseif isone(Δa) && iszero(Δb) && iszero(Δc) # special case (Δa,Δb,Δc) = (1,0,0)
# special case (Δa,Δb,Δc) = (-2,0,0)
elseif isone(-Δa/2) && iszero(Δb) && iszero(Δc)
TSGut marked this conversation as resolved.
Show resolved Hide resolved
M = qr(Q.X^(-Δa÷2)).R
return ApplyArray(\, M, Diagonal(sign.(view(M,band(0))).*Fill(abs.(M[1]),∞)))
# special case (Δa,Δb,Δc) = (0,-2,0)
elseif iszero(Δa) && isone(-Δb/2) && iszero(Δc)
M = qr((I-Q.X)^(-Δb÷2)).R
return ApplyArray(\, M, Diagonal(sign.(view(M,band(0))).*Fill(abs.(M[1]),∞)))
# special case (Δa,Δb,Δc) = (0,0,-2)
elseif iszero(Δa) && iszero(Δb) && isone(-Δc/2)
M = qr((Q.t*I-Q.X)^(-Δc÷2)).R
return ApplyArray(\, M, Diagonal(sign.(view(M,band(0))).*Fill(abs.(M[1]),∞)))
# special case (Δa,Δb,Δc) = (1,0,0)
elseif isone(Δa) && iszero(Δb) && iszero(Δc)
M = cholesky(P.X).U
return ApplyArray(*, Diagonal(Fill(1/M[1],∞)), M)
elseif iszero(Δa) && isone(Δb) && iszero(Δc) # special case (Δa,Δb,Δc) = (0,1,0)
# special case (Δa,Δb,Δc) = (0,1,0)
elseif iszero(Δa) && isone(Δb) && iszero(Δc)
M = cholesky(I-P.X).U
return ApplyArray(*, Diagonal(Fill(1/M[1],∞)), M)
elseif iszero(Δa) && iszero(Δb) && isone(Δc) # special case (Δa,Δb,Δc) = (0,0,1)
M = cholesky(Qt*I-P.X).U
# special case (Δa,Δb,Δc) = (0,0,1)
elseif iszero(Δa) && iszero(Δb) && isone(Δc)
M = cholesky(Q.t*I-P.X).U
return ApplyArray(*, Diagonal(Fill(1/M[1],∞)), M)
TSGut marked this conversation as resolved.
Show resolved Hide resolved
# special case (Δa,Δb,Δc) = (-1,0,0)
elseif isone(-Δa) && iszero(Δb) && iszero(Δc)
M = cholesky(Q.X).U
return ApplyArray(\, M, Diagonal(Fill(M[1],∞)))
TSGut marked this conversation as resolved.
Show resolved Hide resolved
# special case (Δa,Δb,Δc) = (0,-1,0)
elseif iszero(Δa) && isone(-Δb) && iszero(Δc)
M = cholesky(I-Q.X).U
return ApplyArray(\, M, Diagonal(Fill(M[1],∞)))
TSGut marked this conversation as resolved.
Show resolved Hide resolved
# special case (Δa,Δb,Δc) = (0,0,-1)
elseif iszero(Δa) && iszero(Δb) && isone(-Δc)
M = cholesky(Q.t*I-Q.X).U
return ApplyArray(\, M, Diagonal(Fill(M[1],∞)))
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
return ApplyArray(\, M, Diagonal(Fill(M[1],∞)))
return inv(M/M[1])

Copy link
Member Author

@TSGut TSGut Jun 24, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

inv doesn't work when applied to this type. ApplyArray(inv,M/M[1]) works similarly to ApplyArray(\, M, Diagonal(Fill(M[1],∞))) though Julia doesn't seem to know it's Triangular so it causes some odd behavior here and there. This feels like something to fix elsewhere and then we can make it just read inv(M/M[1]) here? Altrenatively I can always just return UpperTriangular(ApplyArray(inv,M/M[1]) for now...

For the record, with ApplyArray(inv,M/M[1]) we get

MemoryLayout(R_0)
LazyArrays.InvLayout{TriangularLayout{'U', 'N', LazyArrays.LazyLayout}}()

with inv(M/M[1]) it errors:

julia>         R_0 = SemiclassicalJacobi(t, 1, 2, 2) \ P
ERROR: MethodError: no method matching Matrix{Float64}(::UniformScaling{Bool}, ::Tuple{Infinities.InfiniteCardinal{0}, Infinities.InfiniteCardinal{0}})

Closest candidates are:
  Matrix{T}(::UndefInitializer, ::Tuple{Union{Infinities.InfiniteCardinal{0}, Infinities.Infinity}, Union{Infinities.InfiniteCardinal{0}, Infinities.Infinity}}) where T
   @ InfiniteArrays C:\Users\timon\.julia\packages\InfiniteArrays\heAfT\src\infarrays.jl:7
  Matrix{T}(::UndefInitializer, ::Tuple{Integer, Union{Infinities.InfiniteCardinal{0}, Infinities.Infinity}}) where T
   @ InfiniteArrays C:\Users\timon\.julia\packages\InfiniteArrays\heAfT\src\infarrays.jl:4
  Matrix{T}(::UndefInitializer, ::Tuple{Union{Infinities.InfiniteCardinal{0}, Infinities.Infinity}, Integer}) where T
   @ InfiniteArrays C:\Users\timon\.julia\packages\InfiniteArrays\heAfT\src\infarrays.jl:5
  ...

Stacktrace:
 [1] inv(A::UpperTriangular{Float64, BroadcastMatrix{Float64, typeof(/), Tuple{…}}})
   @ LinearAlgebra C:\Users\timon\AppData\Local\Programs\Julia-1.10.0\share\julia\stdlib\v1.10\LinearAlgebra\src\triangular.jl:797    
 [2] semijacobi_ldiv_direct(Q::SemiclassicalJacobi{Float64}, P::SemiclassicalJacobi{Float64})
   @ SemiclassicalOrthogonalPolynomials c:\Users\timon\OneDrive\Documents\Work projects\Archive\SemiclassicalOrthogonalPolynomials.jl\src\SemiclassicalOrthogonalPolynomials.jl:322
 [3] semijacobi_ldiv(Q::SemiclassicalJacobi{Float64}, P::SemiclassicalJacobi{Float64})
   @ SemiclassicalOrthogonalPolynomials c:\Users\timon\OneDrive\Documents\Work projects\Archive\SemiclassicalOrthogonalPolynomials.jl\src\SemiclassicalOrthogonalPolynomials.jl:352
 [4] \(A::SemiclassicalJacobi{Float64}, B::SemiclassicalJacobi{Float64})
   @ SemiclassicalOrthogonalPolynomials c:\Users\timon\OneDrive\Documents\Work projects\Archive\SemiclassicalOrthogonalPolynomials.jl\src\SemiclassicalOrthogonalPolynomials.jl:411
 [5] top-level scope
   @ c:\Users\timon\OneDrive\Documents\Work projects\Archive\SemiclassicalOrthogonalPolynomials.jl\test\runtests.jl:663
Some type information was truncated. Use `show(err)` to see complete types.

elseif isinteger(Δa) && isinteger(Δb) && isinteger(Δc) && (Δa >= 0) && (Δb >= 0) && (Δc >= 0)
M = cholesky(Symmetric(P.X^(Δa)*(I-P.X)^(Δb)*(Q.t*I-P.X)^(Δc))).U
return ApplyArray(*, Diagonal(Fill(1/M[1],∞)), M)
elseif isinteger(Δa) && isinteger(Δb) && isinteger(Δc)
M = cholesky(Symmetric(P.X^(Δa)*(I-P.X)^(Δb)*(Qt*I-P.X)^(Δc))).U
return ApplyArray(*, Diagonal(Fill(1/M[1],∞)), M) # match normalization choice P_0(x) = 1
else
error("Implement")
error("Implement modification by ($Δa,$Δb,$Δc)")
end
end

# returns conversion operator from SemiclassicalJacobi P to SemiclassicalJacobi Q.
semijacobi_ldiv(Q::SemiclassicalJacobi, P::SemiclassicalJacobi) = semijacobi_ldiv(Q.t, Q.a, Q.b, Q.c, P)
function semijacobi_ldiv(Qt, Qa, Qb, Qc, P::SemiclassicalJacobi)
@assert Qt ≈ P.t
(Qt ≈ P.t) && (Qa ≈ P.a) && (Qb ≈ P.b) && (Qc ≈ P.c) && return SymTridiagonal(Ones(∞),Zeros(∞))
Δa = Qa-P.a
Δb = Qb-P.b
Δc = Qc-P.c
"""
semijacobi_ldiv_direct(Q::SemiclassicalJacobi, P::SemiclassicalJacobi)

Returns conversion operator from SemiclassicalJacobi `P` to SemiclassicalJacobi `Q`. Integer distances are covered by decomposition methods, for non-integer cases a Lanczos fallback is attempted.
"""
function semijacobi_ldiv(Q::SemiclassicalJacobi, P::SemiclassicalJacobi)
@assert Q.t ≈ P.t
(Q.t ≈ P.t) && (Q.a ≈ P.a) && (Q.b ≈ P.b) && (Q.c ≈ P.c) && return SymTridiagonal(Ones(∞),Zeros(∞))
TSGut marked this conversation as resolved.
Show resolved Hide resolved
Δa = Q.a-P.a
Δb = Q.b-P.b
Δc = Q.c-P.c
if isinteger(Δa) && isinteger(Δb) && isinteger(Δc) # (Δa,Δb,Δc) are integers -> use QR/Cholesky iteratively
if ((isone(Δa)||isone(Δa/2)) && iszero(Δb) && iszero(Δc)) || (iszero(Δa) && (isone(Δb)||isone(Δb/2)) && iszero(Δc)) || (iszero(Δa) && iszero(Δb) && (isone(Δc)||isone(Δc/2)))
M = semijacobi_ldiv_direct(Qt, Qa, Qb, Qc, P)
elseif Δa > 0 # iterative modification by 1
M = ApplyArray(*,semijacobi_ldiv_direct(Qt, Qa, Qb, Qc, SemiclassicalJacobi(Qt, Qa-1-iseven(Δa), Qb, Qc, P)),semijacobi_ldiv(Qt, Qa-1-iseven(Δa), Qb, Qc, P))
if ((isone(abs(Δa))||isone(Δa/2)) && iszero(Δb) && iszero(Δc)) || (iszero(Δa) && (isone(abs(Δb))||isone(Δb/2)) && iszero(Δc)) || (iszero(Δa) && iszero(Δb) && (isone(abs(Δc))||isone(Δc/2)))
return semijacobi_ldiv_direct(Q, P)
elseif Δa > 0 # iterative raising by 1
QQ = SemiclassicalJacobi(Q.t, Q.a-1-iseven(Δa), Q.b, Q.c, P)
return ApplyArray(*,semijacobi_ldiv_direct(Q, QQ),semijacobi_ldiv(QQ, P))
elseif Δb > 0
M = ApplyArray(*,semijacobi_ldiv_direct(Qt, Qa, Qb, Qc, SemiclassicalJacobi(Qt, Qa, Qb-1-iseven(Δb), Qc, P)),semijacobi_ldiv(Qt, Qa, Qb-1-iseven(Δb), Qc, P))
QQ = SemiclassicalJacobi(Q.t, Q.a, Q.b-1-iseven(Δb), Q.c, P)
return ApplyArray(*,semijacobi_ldiv_direct(Q, QQ),semijacobi_ldiv(QQ, P))
elseif Δc > 0
M = ApplyArray(*,semijacobi_ldiv_direct(Qt, Qa, Qb, Qc, SemiclassicalJacobi(Qt, Qa, Qb, Qc-1-iseven(Δc), P)),semijacobi_ldiv(Qt, Qa, Qb, Qc-1-iseven(Δc), P))
QQ = SemiclassicalJacobi(Q.t, Q.a, Q.b, Q.c-1-iseven(Δc), P)
return ApplyArray(*,semijacobi_ldiv_direct(Q, QQ),semijacobi_ldiv(QQ, P))
elseif Δa < 0 # iterative lowering by 1
QQ = SemiclassicalJacobi(Q.t, Q.a+1+iseven(Δa), Q.b, Q.c, P)
return ApplyArray(*,semijacobi_ldiv_direct(Q, QQ),semijacobi_ldiv(QQ, P))
elseif Δb < 0
QQ = SemiclassicalJacobi(Q.t, Q.a, Q.b+1+iseven(Δb), Q.c, P)
return ApplyArray(*,semijacobi_ldiv_direct(Q, QQ),semijacobi_ldiv(QQ, P))
elseif Δc < 0
QQ = SemiclassicalJacobi(Q.t, Q.a, Q.b, Q.c+1+iseven(Δc), P)
return ApplyArray(*,semijacobi_ldiv_direct(Q, QQ),semijacobi_ldiv(QQ, P))
end
else # fallback to Lancos
R = SemiclassicalJacobi(P.t, mod(P.a,-1), mod(P.b,-1), mod(P.c,-1))
R̃ = toclassical(R)
return (P \ R) * _p0(R̃) * (R̃ \ SemiclassicalJacobi(Qt, Qa, Qb, Qc))
return (P \ R) * _p0(R̃) * (R̃ \ Q)
end
end

Expand Down Expand Up @@ -596,5 +658,6 @@ function sumquotient(wP::SemiclassicalJacobiWeight{T},wQ::SemiclassicalJacobiWei
end

include("neg1c.jl")
include("deprecated.jl")

end
3 changes: 3 additions & 0 deletions src/deprecated.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,3 @@
# The old ldiv variants are still supported for now but deprecated and implemented using non-efficient constructors
Base.@deprecate semijacobi_ldiv_direct(Qt, Qa, Qb, Qc, P::SemiclassicalJacobi) semijacobi_ldiv_direct(SemiclassicalJacobi(Qt, Qa, Qb, Qc), P)
Base.@deprecate semijacobi_ldiv(Qt, Qa, Qb, Qc, P::SemiclassicalJacobi) semijacobi_ldiv(SemiclassicalJacobi(Qt, Qa, Qb, Qc), P)
35 changes: 35 additions & 0 deletions test/runtests.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -619,5 +619,40 @@ end
@test sprint(show, W) == "Weighted(SemiclassicalJacobi with weight x^2.3 * (1-x)^5.3 * (2.0-x)^0.4 on 0..1)"
end

@testset "Lowering" begin
@testset "Decrements of 1 via Cholesky" begin
# Cholesky lowering
t = 2.2
P = SemiclassicalJacobi(t, 2, 2, 2)
R_0 = SemiclassicalJacobi(t, 1, 2, 2) \ P
R_1 = SemiclassicalJacobi(t, 2, 1, 2) \ P
R_t = SemiclassicalJacobi(t, 2, 2, 1) \ P

R_0inv = P \ SemiclassicalJacobi(t, 1, 2, 2)
R_1inv = P \ SemiclassicalJacobi(t, 2, 1, 2)
R_tinv = P \ SemiclassicalJacobi(t, 2, 2, 1)

@test ApplyArray(inv,R_0inv)[1:20,1:20] ≈ R_0[1:20,1:20]
@test ApplyArray(inv,R_1inv)[1:20,1:20] ≈ R_1[1:20,1:20]
@test ApplyArray(inv,R_tinv)[1:20,1:20] ≈ R_t[1:20,1:20]
end
@testset "Decrements of 2 via QR" begin
# Cholesky lowering
t = 2.2
P = SemiclassicalJacobi(t, 3, 3, 3)
R_0 = SemiclassicalJacobi(t, 1, 3, 3) \ P
R_1 = SemiclassicalJacobi(t, 3, 1, 3) \ P
R_t = SemiclassicalJacobi(t, 3, 3, 1) \ P

R_0inv = P \ SemiclassicalJacobi(t, 1, 3, 3)
R_1inv = P \ SemiclassicalJacobi(t, 3, 1, 3)
R_tinv = P \ SemiclassicalJacobi(t, 3, 3, 1)

@test ApplyArray(inv,R_0inv)[1:20,1:20] ≈ R_0[1:20,1:20]
@test ApplyArray(inv,R_1inv)[1:20,1:20] ≈ R_1[1:20,1:20]
@test ApplyArray(inv,R_tinv)[1:20,1:20] ≈ R_t[1:20,1:20]
end
end

include("test_derivative.jl")
include("test_neg1c.jl")
Loading