cleanup + makieplot fixes

pablosanjose · Nov 3, 2020 · 21dce3e · 21dce3e
1 parent 6a6657f
commit 21dce3e
Show file tree

Hide file tree

Showing 5 changed files with 6 additions and 168 deletions.
diff --git a/src/diagonalizer.jl b/src/diagonalizer.jl
@@ -94,154 +94,4 @@ method_matrixtype(::AbstractDiagonalizeMethod, h) = flatten
 # Type of states in bandstructure. Should be Matrix view, because eigvecs are dense and we need to support degeneracies
 method_matviewtype(::AbstractDiagonalizeMethod, ::AbstractMatrix{T}) where {T} = subarray_matrix_type(orbitaltype(T))
 
-subarray_matrix_type(::Type{M}) where {M} = typeof(view(Matrix{eltype(M)}(undef, 2, 2), :, 1:0))
-
-# #######################################################################
-# # shift and invert methods
-# #######################################################################
-
-# function shiftandinvert(matrix::AbstractMatrix{Tv}, origin) where {Tv}
-#     cols, rows = size(matrix)
-#     # Shift away from real axis to avoid pivot point error in factorize
-#     matrix´ = diagshift!(parent(matrix), origin + im)
-#     F = factorize(matrix´)
-#     lmap = LinearMap{Tv}((x, y) -> ldiv!(x, F, y), cols, rows,
-#                          ismutating = true, ishermitian = false)
-#     return lmap
-# end
-
-# function diagshift!(matrix::AbstractMatrix, origin)
-#     matrix´ = parent(matrix)
-#     vals = nonzeros(matrix´)
-#     rowval = rowvals(matrix´)
-#     for col in 1:size(matrix, 2)
-#         found_diagonal = false
-#         for ptr in nzrange(matrix´, col)
-#             if col == rowval[ptr]
-#                 found_diagonal = true
-#                 vals[ptr] -= origin * I  # To respect non-scalar eltypes
-#                 break
-#             end
-#         end
-#         found_diagonal || throw(error("Sparse work matrix must include the diagonal. Possible bug in `similarmatrix`."))
-#     end
-#     return matrix
-# end
-
-# function invertandshift(ϵ::Vector{T}, origin) where {T}
-#     ϵ´ = similar(ϵ, real(T))
-#     ϵ´ .= real(inv.(ϵ) .+ (origin + im))  # Caution: we assume a real spectrum
-#     return ϵ´
-# end
-
-# #######################################################################
-# # Codiagonalizer
-# #######################################################################
-
-# ## Codiagonalizer
-# ## Uses velocity operators along different directions. If not enough, use finite differences
-# ## along mesh directions
-# struct Codiagonalizer{T,F<:Function}
-#     comatrix::F
-#     matrixindices::UnitRange{Int}
-#     degtol::T
-#     rangesA::Vector{UnitRange{Int}} # Prealloc buffer for degeneray ranges
-#     rangesB::Vector{UnitRange{Int}} # Prealloc buffer for degeneray ranges
-#     perm::Vector{Int}               # Prealloc for sortperm!
-# end
-
-# # mapping = missing is assumed when h is a Function that generates matrices, instead of a Hamiltonian or ParametricHamiltonian
-# function codiagonalizer(h, matrix::AbstractMatrix{T}, mesh, mapping) where {T}
-#     sample_phiparams = map_phiparams(mapping, first(vertices(mesh)))
-#     dirs = codiag_directions(sample_phiparams)
-#     degtol = sqrt(eps(real(eltype(T))))
-#     delta = meshdelta(mesh)
-#     delta = iszero(delta) ? degtol : delta
-#     comatrix, matrixindices = codiag_function(h, matrix, mapping, dirs, delta)
-#     return Codiagonalizer(comatrix, matrixindices, degtol, UnitRange{Int}[], UnitRange{Int}[], Int[])
-# end
-
-# function codiag_function(h::Union{Hamiltonian,ParametricHamiltonian}, matrix, mapping, dirs, delta)
-#     hdual = dual_if_parametric(h)
-#     matrixdual = dualarray(matrix)
-#     anyold = anyoldmatrix(matrix)
-#     ndirs = length(dirs)
-#     matrixindices = 1:(ndirs + ndirs + 1)
-#     comatrix(meshϕs, n) =
-#         if n <= ndirs # automatic differentiation using dual numbers
-#             ϕsparams = dual_phisparams(map_phiparams(mapping, meshϕs), dirs[n])
-#             dualpart.(bloch!(matrixdual, hdual, ϕsparams))
-#         elseif n - ndirs <= ndirs # resort to finite differences
-#             ϕsparams = delta_phisparams(map_phiparams(mapping, meshϕs), delta * dirs[n - ndirs])
-#             bloch!(matrix, h, ϕsparams)
-#         else # use a fixed arbitrary matrix
-#             anyold
-#         end
-#     return comatrix, matrixindices
-# end
-
-# dual_if_parametric(ph::ParametricHamiltonian) = Dual(ph)
-# dual_if_parametric(h::Hamiltonian) = h
-
-# # In the Function case we cannot know what directions to scan (arguments of matrixf). Also,
-# # we cannot be sure that dual numbers propagate. We thus restrict to finite differences in the mesh
-# # Note that mapping is already wrapped into matrixf in the calling bandstructure(::Function,...)
-# function codiag_function(matrixf::Function, matrix, mapping::Missing, meshdirs, delta)
-#     anyold = anyoldmatrix(matrix)
-#     ndirs = length(meshdirs)
-#     matrixindices = 1:(ndirs + 1)
-#     comatrix(meshϕs, n) =
-#         if n <= ndirs # finite differences
-#             matrixf(meshϕs + delta * meshdirs[n])
-#         else # use a fixed arbitrary matrix
-#             anyold
-#         end
-#     return comatrix, matrixindices
-# end
-
-# val_length(::SVector{N}, nt::NamedTuple) where {N} = Val(N+length(nt))
-
-# codiag_directions(phiparams) = codiag_directions(val_length(phiparams...), phiparams)
-
-# function codiag_directions(::Val{L}, phiparams, direlements = 0:1) where {L}
-#     directions = vec(SVector{L,Int}.(Iterators.product(ntuple(_ -> direlements, Val(L))...)))
-#     mask_dirs!(directions, phiparams)
-#     filter!(ispositive, directions)
-#     unique!(directions)
-#     sort!(directions, by = norm)
-#     return directions
-# end
-
-# # Zeros out any direction that cannot modify a param because it is not a number
-# function mask_dirs!(dirs::Vector{S}, pp) where {L,S<:SVector{L}}
-#     valparams = values(last(pp))
-#     valids = valparams .!= maybe_sum.(valparams, 1)
-#     n = length(first(pp))
-#     mask = SVector(ntuple(i -> i <= n || valids[i - n] , Val(L)))
-#     map!(dir -> mask .* dir, dirs, dirs)
-#     return dirs
-# end
-
-# dual_phisparams(ϕs::SVector{N}, dir) where {N} = Dual.(ϕs, frontsvec(dir, Val(N)))
-# dual_phisparams(params::NamedTuple, dir) = NamedTuple{keys(params)}(maybe_dual.(values(params), tailtuple(dir, Val(length(params)))))
-# dual_phisparams((ϕs, params)::Tuple, dir) = (dual_phisparams(ϕs, dir), dual_phisparams(params, dir))
-
-# maybe_dual(param::Number, ε) = Dual(param, ε)
-# maybe_dual(param, ε) = param
-
-# delta_phisparams(ϕs::SVector{N}, dir) where {N} = ϕs + frontsvec(dir, Val(N))
-# delta_phisparams(params::NamedTuple, dir) = NamedTuple{keys(params)}(maybe_sum.(values(params), tailtuple(dir, Val(length(params)))))
-# delta_phisparams((ϕs, params)::Tuple, dir) = (delta_phisparams(ϕs, dir), delta_phisparams(params, dir))
-
-# maybe_sum(param::Number, ε) = param + ε
-# maybe_sum(param, ε) = param
-
-# meshdelta(mesh::BandMesh{<:Any,T}) where {T} = T(0.1) * norm(first(minmax_edge(mesh)))
-
-# function anyoldmatrix(matrix::SparseMatrixCSC, rng = MersenneTwister(1))
-#     s = copy(matrix)
-#     rand!(rng, nonzeros(s))
-#     return s
-# end
-
-# anyoldmatrix(m::DenseArray, rng = MersenneTwister(1)) = rand!(rng, copy(m))
+subarray_matrix_type(::Type{M}) where {M} = typeof(view(Matrix{eltype(M)}(undef, 2, 2), :, 1:0))
diff --git a/src/hamiltonian.jl b/src/hamiltonian.jl
@@ -724,12 +724,6 @@ SparseArrays.issparse(h::Hamiltonian{LA,L,M,A}) where {LA,L,M,A} = false
 
 Base.parent(h::Hamiltonian) = h
 
-# Dual numbers #
-
-DualNumbers.Dual(h::Hamiltonian) = Hamiltonian(h.lattice, Dual.(h.harmonics), h.orbitals)
-
-DualNumbers.Dual(h::HamiltonianHarmonic) = HamiltonianHarmonic(h.dn, dualarray(h.h))
-
 # Iterators #
 
 function nonzero_indices(h::Hamiltonian, rowrange = 1:size(h, 1), colrange = 1:size(h, 2))

diff --git a/src/parametric.jl b/src/parametric.jl
@@ -201,9 +201,6 @@ Base.size(ph::ParametricHamiltonian, n...) = size(ph.h, n...)
 
 Base.eltype(ph::ParametricHamiltonian) = eltype(ph.h)
 
-DualNumbers.Dual(p::ParametricHamiltonian) =
-    ParametricHamiltonian(Dual(p.baseh), Dual(p.h), p.modifiers, p.ptrdata, p.allptrs, p.parameters)
-
 #######################################################################
 # bloch! for parametric
 #######################################################################

diff --git a/src/plot_makie.jl b/src/plot_makie.jl
@@ -263,7 +263,7 @@ end
 
 @recipe(BandPlot2D, bandstructure) do scene
     Theme(
-    linewidth = 3.0,
+    linethickness = 3.0,
     wireframe = true,
     colors = map(t -> RGBAf0((0.8 .* t)...),
         ((0.973, 0.565, 0.576), (0.682, 0.838, 0.922), (0.742, 0.91, 0.734),
@@ -281,15 +281,16 @@ function plot!(plot::BandPlot2D)
         vertices = band.verts
         simplices = band.simpinds
         linesegments!(plot, (t -> vertices[first(t)] => vertices[last(t)]).(simplices),
-                      linewidth = plot[:linewidth], color = color)
+                      linewidth = plot[:linethickness][], color = color)
     end
     return plot
  end
 
 @recipe(BandPlot3D, bandstructure) do scene
     Theme(
-    linewidth = 1.0,
+    linethickness = 1.0,
     wireframe = true,
+    linedarken = 0.5,
     ssao = true, ambient = Vec3f0(0.55), diffuse = Vec3f0(0.4),
     colors = map(t -> RGBAf0(t...),
         ((0.973, 0.565, 0.576), (0.682, 0.838, 0.922), (0.742, 0.91, 0.734),
@@ -313,7 +314,7 @@ function plot!(plot::BandPlot3D)
                 ssao = plot[:ssao][], ambient = plot[:ambient][], diffuse = plot[:diffuse][])
             if plot[:wireframe][]
                 edgevertices = collect(Quantica.edgevertices(band))
-                linesegments!(plot, edgevertices, color = darken(color, 0.4), linewidth = plot[:linewidth])
+                linesegments!(plot, edgevertices, color = darken(color, plot[:linedarken][]), linewidth = plot[:linethickness][])
             end
         end
     end

diff --git a/src/tools.jl b/src/tools.jl
@@ -264,10 +264,6 @@ function append_slice!(dest::AbstractArray, src::AbstractArray{T,N}, Rsrc::Carte
     return dest
 end
 
-dualarray(a::DenseMatrix) = map(x->Dual.(x, 0), a)
-# Need to preserve stored zeros, so we have to treat sparse case as special
-dualarray(s::SparseMatrixCSC) = SparseMatrixCSC(s.m, s.n, s.colptr, s.rowval, map(x->Dual.(x, 0), s.nzval))
-
 ######################################################################
 # convert a matrix/number block to a matrix/inlinematrix string
 ######################################################################