Merge pull request #2 from numericalEFT/treedivBZ

Treediv bz

.gitignore

@@ -1,4 +1,32 @@
-*.jl.*.cov
+
+
+# Files generated by invoking Julia with --code-coverage
 *.jl.cov
+*.jl.*.cov
+
+# Files generated by invoking Julia with --track-allocation
 *.jl.mem
-/docs/build/
+*.jl.*.mem
+
+# System-specific files and directories generated by the BinaryProvider and BinDeps packages
+# They contain absolute paths specific to the host computer, and so should not be committed
+deps/deps.jl
+deps/build.log
+deps/downloads/
+deps/usr/
+deps/src/
+
+# Build artifacts for creating documentation generated by the Documenter package
+docs/build/
+docs/site/
+/run/
+
+# File generated by Pkg, the package manager, based on a corresponding Project.toml
+# It records a fixed state of all packages used by the project. As such, it should not be
+# committed for packages, but should be committed for applications that require a static
+# environment.
+Manifest.toml
+
+*.DS_Store
+
+.ipynb_checkpoints

Project.toml

@@ -3,6 +3,11 @@ uuid = "3d4565c3-192b-4227-89f7-5476f913ab49"
 authors = ["Tao Wang", "Xiansheng Cai"]
 version = "0.1.0"
 
+[deps]
+AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+
 [compat]
 julia = "1.6"
 

example/testabstracttree.jl

@@ -0,0 +1,23 @@
+using AbstractTrees
+
+struct MyNode
+    val::Int
+    children::Vector{MyNode}
+end
+
+AbstractTrees.children(node::MyNode) = node.children
+
+function MyNode(val::Int)
+    val = val
+    children = Vector{MyNode}()
+    for i in 1:val
+        push!(children, MyNode(val-1))
+    end
+
+    return MyNode(val, children)
+end
+
+tree = MyNode(3)
+
+print_tree(tree)
+

example/tree.jl

@@ -0,0 +1,153 @@
+
+#--- test tree divided BZ ---
+
+using Random
+using LinearAlgebra
+using Plots
+
+struct SquareNod
+    gridpoints::Matrix{Float64}
+    sons::Vector{SquareNod}
+end
+
+function minarea(sn::SquareNod)
+    return (isempty(sn.sons)) ? area(sn) : min([minarea(ssn) for ssn in sn.sons]...)
+end
+
+function SquareNod!(gridpoints, weightcap::Float64, gpcontainer; treelevelcap = 5)
+
+    K00, K01, K10, K11 = Matrix(gridpoints[1,:]'), Matrix(gridpoints[2,:]'), Matrix(gridpoints[3,:]'), Matrix(gridpoints[4,:]')
+    K0h = 0.5 .* (K00 .+ K01)
+    Kh0 = 0.5 .* (K00 .+ K10)
+    K1h = 0.5 .* (K10 .+ K11)
+    Kh1 = 0.5 .* (K01 .+ K11)
+    Khh = 0.5 .* (K00 .+ K11)
+
+    push!(gpcontainer, K0h)
+    push!(gpcontainer, Kh0)
+    # push!(gpcontainer, K1h)
+    # push!(gpcontainer, Kh1)
+    push!(gpcontainer, Khh)
+
+    sons = []
+    if treelevelcap > 1 && weight(gridpoints) > weightcap
+        sn00 = SquareNod!(Matrix([K00; K0h; Kh0; Khh]), weightcap, gpcontainer;treelevelcap= treelevelcap-1)
+        sn01 = SquareNod!(Matrix([K0h; K01; Khh; Kh1]), weightcap, gpcontainer;treelevelcap= treelevelcap-1)
+        sn10 = SquareNod!(Matrix([Kh0; Khh; K10; K1h]), weightcap, gpcontainer;treelevelcap= treelevelcap-1)
+        sn11 = SquareNod!(Matrix([Khh; Kh1; K1h; K11]), weightcap, gpcontainer;treelevelcap= treelevelcap-1)
+
+        push!(sons, sn00)
+        push!(sons, sn01)
+        push!(sons, sn10)
+        push!(sons, sn11)
+    end
+
+    return SquareNod(gridpoints, sons)
+end
+
+function SquareNod(gridpoints, weightcap::Float64; treelevelcap = 5)
+
+    K00, K01, K10, K11 = Matrix(gridpoints[1,:]'), Matrix(gridpoints[2,:]'), Matrix(gridpoints[3,:]'), Matrix(gridpoints[4,:]')
+    K0h = 0.5 .* (K00 .+ K01)
+    Kh0 = 0.5 .* (K00 .+ K10)
+    K1h = 0.5 .* (K10 .+ K11)
+    Kh1 = 0.5 .* (K01 .+ K11)
+    Khh = 0.5 .* (K00 .+ K11)
+
+    sons = []
+    if treelevelcap > 1 && weight(gridpoints) > weightcap
+        sn00 = SquareNod(Matrix([K00; K0h; Kh0; Khh]), weightcap;treelevelcap= treelevelcap-1)
+        sn01 = SquareNod(Matrix([K0h; K01; Khh; Kh1]), weightcap;treelevelcap= treelevelcap-1)
+        sn10 = SquareNod(Matrix([Kh0; Khh; K10; K1h]), weightcap;treelevelcap= treelevelcap-1)
+        sn11 = SquareNod(Matrix([Khh; Kh1; K1h; K11]), weightcap;treelevelcap= treelevelcap-1)
+
+        push!(sons, sn00)
+        push!(sons, sn01)
+        push!(sons, sn10)
+        push!(sons, sn11)
+    end
+
+    return SquareNod(gridpoints, sons)
+end
+
+@inline function density(K::AbstractVector)
+    me = 0.5
+    T = 0.01
+
+    # μ = 1.0
+    # k = norm(K)
+    # ϵ = k^2 / (2me)
+    # dϵdk = k / me
+
+    μ = 0.5
+    ϵ = -(cos(π*K[1]) + cos(π*K[2])) / (2me)
+    dϵdk = sqrt(sin(π*K[1])^2 + sin(π*K[2])^2) / (2me)
+
+    return (exp((ϵ - μ) / T) / T)/(exp((ϵ - μ) / T) + 1.0)^2 * dϵdk
+end
+
+function area(gridpoints)
+    a = norm(gridpoints[2,:]-gridpoints[1,:])
+    b = norm(gridpoints[3,:]-gridpoints[1,:])
+    return a*b
+end
+
+area(sn::SquareNod) = area(sn.gridpoints)
+
+function weight(gridpoints; N=400)
+    result = 0.0
+    a = area(gridpoints)
+
+    K1, K2 = gridpoints[2,:]-gridpoints[1,:], gridpoints[3,:]-gridpoints[1,:]
+    for i in 1:N
+        K = gridpoints[1, :] .+ rand() .* K1 .+ rand() .* K2
+        result += density(K)
+    end
+
+    return result * a
+end
+
+weight(sn::SquareNod; N=100) = weight(sn.gridpoints; N=N)
+
+gridpoints = Matrix([-1.0 -1.0 ; -1 1 ; 1 -1 ; 1 1])
+K00, K01, K10, K11 = Matrix(gridpoints[1,:]'), Matrix(gridpoints[2,:]'), Matrix(gridpoints[3,:]'), Matrix(gridpoints[4,:]')
+gpcontainer = [K00,]
+
+sn = SquareNod!(gridpoints, weight(gridpoints)/128^2, gpcontainer; treelevelcap = 12)
+
+println(area(sn))
+println(weight(sn))
+
+function Base.print(sn::SquareNod; treelevel = 1)
+    for i in 1:treelevel
+        print(" ")
+    end
+    print("area=$(area(sn)),")
+    print("weight=$(weight(sn)),")
+    print("gridpoints:",sn.gridpoints)
+    print("\n")
+
+    for ssn in sn.sons
+        print(ssn, ;treelevel = treelevel + 1)
+    end
+end
+
+# print(sn)
+
+println("minarea = $(minarea(sn))")
+# print(gpcontainer)
+println("length:",length(gpcontainer))
+
+println("compress:", minarea(sn)*length(gpcontainer)/4)
+
+X, Y = zeros(Float64, length(gpcontainer)), zeros(Float64, length(gpcontainer))
+for (Ki, K) in enumerate(gpcontainer)
+    X[Ki], Y[Ki] = K[1], K[2]
+end
+
+# p = plot(legend = false, size = (1024, 1024))
+# scatter!(p, X, Y, marker = :cross, markersize = 2)
+# savefig(p, "run/grid.pdf")
+# display(p)
+# readline()
+

example/treegrid.jl

@@ -0,0 +1,98 @@
+#--- test tree divided BZ ---
+
+using Random
+using LinearAlgebra
+using Plots
+
+struct MeshNod
+    depth::Int
+    pos::Vector{Int}
+    sons::Vector{MeshNod}
+end
+
+struct TreeMesh
+    dim::Int
+    root::MeshNod
+    latticevectors::Matrix{Float64}
+    gridpoints::Matrix{Float64}
+end
+
+function TreeMesh(latticevectors::Matrix, weightcap::Float64; depthcap = 5)
+
+end
+
+@inline function density(K::AbstractVector)
+    me = 0.5
+    T = 0.01
+
+    μ = 1.0
+    k = norm(K)
+    ϵ = k^2 / (2me)
+    dϵdk = k / me
+
+    # μ = 0.5
+    # ϵ = -(cos(π*K[1]) + cos(π*K[2])) / (2me)
+    # dϵdk = sqrt(sin(π*K[1])^2 + sin(π*K[2])^2) / (2me)
+
+    return (exp((ϵ - μ) / T) / T)/(exp((ϵ - μ) / T) + 1.0)^2 * dϵdk
+end
+
+function area(gridpoints)
+    a = norm(gridpoints[2,:]-gridpoints[1,:])
+    b = norm(gridpoints[3,:]-gridpoints[1,:])
+    return a*b
+end
+
+area(sn::SquareNod) = area(sn.gridpoints)
+
+function weight(gridpoints; N=100)
+    result = 0.0
+    a = area(gridpoints)
+
+    K1, K2 = gridpoints[2,:]-gridpoints[1,:], gridpoints[3,:]-gridpoints[1,:]
+    for i in 1:N
+        K = gridpoints[1, :] .+ rand() .* K1 .+ rand() .* K2
+        result += density(K)
+    end
+
+    return result * a
+end
+
+weight(sn::SquareNod; N=100) = weight(sn.gridpoints; N=N)
+
+gridpoints = Matrix([-1.0 -1.0 ; -1 1 ; 1 -1 ; 1 1])
+K00, K01, K10, K11 = Matrix(gridpoints[1,:]'), Matrix(gridpoints[2,:]'), Matrix(gridpoints[3,:]'), Matrix(gridpoints[4,:]')
+gpcontainer = [K00, K01, K10, K11]
+
+sn = SquareNod!(gridpoints, weight(gridpoints)/500, gpcontainer; treelevelcap = 10)
+
+println(area(sn))
+println(weight(sn))
+
+function Base.print(sn::SquareNod; treelevel = 1)
+    for i in 1:treelevel
+        print(" ")
+    end
+    print("area=$(area(sn)),")
+    print("weight=$(weight(sn)),")
+    print("gridpoints:",sn.gridpoints)
+    print("\n")
+
+    for ssn in sn.sons
+        print(ssn, ;treelevel = treelevel + 1)
+    end
+end
+
+print(sn)
+# print(gpcontainer)
+println("length:",length(gpcontainer))
+
+X, Y = zeros(Float64, length(gpcontainer)), zeros(Float64, length(gpcontainer))
+for (Ki, K) in enumerate(gpcontainer)
+    X[Ki], Y[Ki] = K[1], K[2]
+end
+
+p = scatter(X, Y)
+display(p)
+readline()
+

example/treenode.jl

@@ -0,0 +1,66 @@
+using SpaceGrid
+using SpaceGrid.AbstractTrees, SpaceGrid.GridTree, SpaceGrid.BaseMesh
+using Plots
+using LinearAlgebra
+
+mesh = UniformMesh{2, 4}([0.0, 0.0], [0 1 ; 1 0])
+show(mesh)
+
+naiveisfine(depth, pos) = depth >= 2
+
+# tree = GridNode{2}(naiveisfine)
+# print_tree(tree)
+
+# for node in PostOrderDFS(tree)
+#     if isempty(node.children)
+#         println(node.pos)
+#     end
+# end
+
+function dispersion(k)
+    me = 0.5
+    μ = 1.0
+    return norm(k)^2/(2me)-μ
+
+    # t = 0.5
+    # μ = 0.0
+    # return sum([t * cos(π*p) for p in k]) + μ
+ end
+
+function density(K, latvec)
+    DIM = length(K)
+
+    T = 0.001
+
+    kpoints = [K, ]
+    for i in 1:DIM
+        push!(kpoints, K .+ latvec[i, :])
+        push!(kpoints, K .- latvec[i, :])
+    end
+
+    ϵ = minimum([dispersion(k) for k in kpoints])
+
+    return 1 / (exp((ϵ) / T) + 1.0)
+    # return 1 / ((π * T)^2 + ϵ^2)
+end
+
+latvec = [2 0; 1 sqrt(3)]
+# latvec = [2 0; 0 2]
+# tg = uniformtreegrid(naiveisfine, latvec; N = 3)
+# tg = treegridfromdensity(density, latvec; rtol = 1e-4, maxdepth = 5)
+tg = treegridfromdensity(k->density(k, latvec), latvec; rtol = 1e-1, maxdepth = 6)
+
+X, Y = zeros(Float64, size(tg)), zeros(Float64, size(tg))
+for (pi, p) in enumerate(tg)
+    X[pi], Y[pi] = p[1], p[2]
+    # println(p)
+end
+
+println("size:$(size(tg))")
+println("length:$(length(tg))")
+println("efficiency:$(efficiency(tg))")
+
+# p = plot(legend = false, size = (1024, 1024), xlim = (-1, 1), ylim = (-1, 1))
+p = plot(legend = false, size = (1024, 1024), xlim = (-1.5, 1.5), ylim = (-1.5, 1.5))
+scatter!(p, X, Y, marker = :cross, markersize = 2)
+savefig(p, "run/tg.pdf")