text:
module rbc_1
using SparseArrays, LinearAlgebra, DifferentiableStateSpaceModels, LaTeXStrings, ModelingToolkit
const n_y = 3
const n_x = 2
const n = 5
const n_p = 3
const n_ϵ = 1
const n_z = 2
const functions_type() = DenseFunctions()
const η = reshape([0; -1], 2, 1)
const Q = [1.0 0.0 0.0 0.0 0.0; 0.0 0.0 1.0 0.0 0.0]
const Γ! = (var"##MTIIPVar#368", var"##MTKArg#364", var"##MTKArg#365", var"##MTKArg#366")->begin
        @inbounds begin
                let (α, β, ρ, δ, σ, Ω_1) = (var"##MTKArg#364"[1], var"##MTKArg#364"[2], var"##MTKArg#364"[3], var"##MTKArg#365"[1], var"##MTKArg#365"[2], var"##MTKArg#365"[3])
                    var"##MTIIPVar#368"[1] = σ
                end
            end
        nothing
    end
const Γ_p! = (var"##MTIIPVar#375", var"##MTKArg#371", var"##MTKArg#372", var"##MTKArg#373")->begin
        @inbounds begin
                let (α, β, ρ, δ, σ, Ω_1) = (var"##MTKArg#371"[1], var"##MTKArg#371"[2], var"##MTKArg#371"[3], var"##MTKArg#372"[1], var"##MTKArg#372"[2], var"##MTKArg#372"[3])
                end
            end
        nothing
    end
const Ω! = (var"##MTIIPVar#382", var"##MTKArg#378", var"##MTKArg#379", var"##MTKArg#380")->begin
        @inbounds begin
                let (α, β, ρ, δ, σ, Ω_1) = (var"##MTKArg#378"[1], var"##MTKArg#378"[2], var"##MTKArg#378"[3], var"##MTKArg#379"[1], var"##MTKArg#379"[2], var"##MTKArg#379"[3])
                    var"##MTIIPVar#382"[1] = Ω_1
                    var"##MTIIPVar#382"[2] = Ω_1
                end
            end
        nothing
    end
const Ω_p! = (var"##MTIIPVar#389", var"##MTKArg#385", var"##MTKArg#386", var"##MTKArg#387")->begin
        @inbounds begin
                let (α, β, ρ, δ, σ, Ω_1) = (var"##MTKArg#385"[1], var"##MTKArg#385"[2], var"##MTKArg#385"[3], var"##MTKArg#386"[1], var"##MTKArg#386"[2], var"##MTKArg#386"[3])
                end
            end
        nothing
    end
const H! = (var"##MTIIPVar#402", var"##MTKArg#392", var"##MTKArg#393", var"##MTKArg#394", var"##MTKArg#395", var"##MTKArg#396", var"##MTKArg#397", var"##MTKArg#398", var"##MTKArg#399", var"##MTKArg#400")->begin
        @inbounds begin
                let (c_p, q_p, i_p, c, q, i, c_ss, q_ss, i_ss, k_p, z_p, k, z, k_ss, z_ss, α, β, ρ, δ, σ, Ω_1) = (var"##MTKArg#392"[1], var"##MTKArg#392"[2], var"##MTKArg#392"[3], var"##MTKArg#393"[1], var"##MTKArg#393"[2], var"##MTKArg#393"[3], var"##MTKArg#394"[1], var"##MTKArg#394"[2], var"##MTKArg#394"[3], var"##MTKArg#395"[1], var"##MTKArg#395"[2], var"##MTKArg#396"[1], var"##MTKArg#396"[2], var"##MTKArg#397"[1], var"##MTKArg#397"[2], var"##MTKArg#398"[1], var"##MTKArg#398"[2], var"##MTKArg#398"[3], var"##MTKArg#399"[1], var"##MTKArg#399"[2], var"##MTKArg#399"[3])
                    var"##MTIIPVar#402"[1] = (-)((/)(1, c), (*)((/)(β, c_p), (+)((*)((*)(α, (exp)(z_p)), (^)(k_p, (-)(α, 1))), (-)(1, δ))))
                    var"##MTIIPVar#402"[2] = (-)((-)((+)(c, k_p), (*)((-)(1, δ), k)), q)
                    var"##MTIIPVar#402"[3] = (-)(q, (*)((exp)(z), (^)(k, α)))
                    var"##MTIIPVar#402"[4] = (-)(z_p, (*)(ρ, z))
                    var"##MTIIPVar#402"[5] = (-)(i, (-)(k_p, (*)((-)(1, δ), k)))
                end
            end
        nothing
    end
const H_yp! = (var"##MTIIPVar#411", var"##MTKArg#405", var"##MTKArg#406", var"##MTKArg#407", var"##MTKArg#408", var"##MTKArg#409")->begin
        @inbounds begin
                let (c, q, i, k, z, α, β, ρ, δ, σ, Ω_1) = (var"##MTKArg#405"[1], var"##MTKArg#405"[2], var"##MTKArg#405"[3], var"##MTKArg#406"[1], var"##MTKArg#406"[2], var"##MTKArg#407"[1], var"##MTKArg#407"[2], var"##MTKArg#407"[3], var"##MTKArg#408"[1], var"##MTKArg#408"[2], var"##MTKArg#408"[3])
                    var"##MTIIPVar#411"[1] = (*)((+)(1, (*)(-1, δ), (*)((exp)(z), α, (^)(k, (+)(-1, α)))), (^)((inv)(c), 2), β)
                end
            end
        nothing
    end
const H_y! = (var"##MTIIPVar#420", var"##MTKArg#414", var"##MTKArg#415", var"##MTKArg#416", var"##MTKArg#417", var"##MTKArg#418")->begin
        @inbounds begin
                let (c, q, i, k, z, α, β, ρ, δ, σ, Ω_1) = (var"##MTKArg#414"[1], var"##MTKArg#414"[2], var"##MTKArg#414"[3], var"##MTKArg#415"[1], var"##MTKArg#415"[2], var"##MTKArg#416"[1], var"##MTKArg#416"[2], var"##MTKArg#416"[3], var"##MTKArg#417"[1], var"##MTKArg#417"[2], var"##MTKArg#417"[3])
                    var"##MTIIPVar#420"[1] = (*)(-1, (^)((inv)(c), 2))
                    var"##MTIIPVar#420"[2] = 1
                    var"##MTIIPVar#420"[7] = -1
                    var"##MTIIPVar#420"[8] = 1
                    var"##MTIIPVar#420"[15] = 1
                end
            end
        nothing
    end
const H_xp! = (var"##MTIIPVar#429", var"##MTKArg#423", var"##MTKArg#424", var"##MTKArg#425", var"##MTKArg#426", var"##MTKArg#427")->begin
        @inbounds begin
                let (c, q, i, k, z, α, β, ρ, δ, σ, Ω_1) = (var"##MTKArg#423"[1], var"##MTKArg#423"[2], var"##MTKArg#423"[3], var"##MTKArg#424"[1], var"##MTKArg#424"[2], var"##MTKArg#425"[1], var"##MTKArg#425"[2], var"##MTKArg#425"[3], var"##MTKArg#426"[1], var"##MTKArg#426"[2], var"##MTKArg#426"[3])
                    var"##MTIIPVar#429"[1] = (*)((*)((*)((*)(-1, (exp)(z)), (^)((inv)(c), 1)), α), (*)((*)(β, (^)(k, (+)(-2, α))), (+)(-1, α)))
                    var"##MTIIPVar#429"[2] = 1
                    var"##MTIIPVar#429"[5] = -1
                    var"##MTIIPVar#429"[6] = (*)((*)((*)((*)(-1, (exp)(z)), (^)((inv)(c), 1)), α), (*)(β, (^)(k, (+)(-1, α))))
                    var"##MTIIPVar#429"[9] = 1
                end
            end
        nothing
    end
const H_x! = (var"##MTIIPVar#438", var"##MTKArg#432", var"##MTKArg#433", var"##MTKArg#434", var"##MTKArg#435", var"##MTKArg#436")->begin
        @inbounds begin
                let (c, q, i, k, z, α, β, ρ, δ, σ, Ω_1) = (var"##MTKArg#432"[1], var"##MTKArg#432"[2], var"##MTKArg#432"[3], var"##MTKArg#433"[1], var"##MTKArg#433"[2], var"##MTKArg#434"[1], var"##MTKArg#434"[2], var"##MTKArg#434"[3], var"##MTKArg#435"[1], var"##MTKArg#435"[2], var"##MTKArg#435"[3])
                    var"##MTIIPVar#438"[2] = (*)(-1, (+)(1, (*)(-1, δ)))
                    var"##MTIIPVar#438"[3] = (*)(-1, (exp)(z), α, (^)(k, (+)(-1, α)))
                    var"##MTIIPVar#438"[5] = (+)(1, (*)(-1, δ))
                    var"##MTIIPVar#438"[8] = (*)(-1, (exp)(z), (^)(k, α))
                    var"##MTIIPVar#438"[9] = (*)(-1, ρ)
                end
            end
        nothing
    end
const H_yp_p! = (var"##MTIIPVar#447", var"##MTKArg#441", var"##MTKArg#442", var"##MTKArg#443", var"##MTKArg#444", var"##MTKArg#445")->begin
        @inbounds begin
                let (c, q, i, k, z, α, β, ρ, δ, σ, Ω_1) = (var"##MTKArg#441"[1], var"##MTKArg#441"[2], var"##MTKArg#441"[3], var"##MTKArg#442"[1], var"##MTKArg#442"[2], var"##MTKArg#443"[1], var"##MTKArg#443"[2], var"##MTKArg#443"[3], var"##MTKArg#444"[1], var"##MTKArg#444"[2], var"##MTKArg#444"[3])
                    (var"##MTIIPVar#447"[1])[1] = (*)((^)((inv)(c), 2), β, (+)((*)((exp)(z), (^)(k, (+)(-1, α))), (*)((log)(k), (exp)(z), α, (^)(k, (+)(-1, α)))))
                    (var"##MTIIPVar#447"[2])[1] = (*)((+)(1, (*)(-1, δ), (*)((exp)(z), α, (^)(k, (+)(-1, α)))), (^)((inv)(c), 2))
                end
            end
        nothing
    end
const H_y_p! = (var"##MTIIPVar#456", var"##MTKArg#450", var"##MTKArg#451", var"##MTKArg#452", var"##MTKArg#453", var"##MTKArg#454")->begin
        @inbounds begin
                let (c, q, i, k, z, α, β, ρ, δ, σ, Ω_1) = (var"##MTKArg#450"[1], var"##MTKArg#450"[2], var"##MTKArg#450"[3], var"##MTKArg#451"[1], var"##MTKArg#451"[2], var"##MTKArg#452"[1], var"##MTKArg#452"[2], var"##MTKArg#452"[3], var"##MTKArg#453"[1], var"##MTKArg#453"[2], var"##MTKArg#453"[3])
                end
            end
        nothing
    end
const H_xp_p! = (var"##MTIIPVar#465", var"##MTKArg#459", var"##MTKArg#460", var"##MTKArg#461", var"##MTKArg#462", var"##MTKArg#463")->begin
        @inbounds begin
                let (c, q, i, k, z, α, β, ρ, δ, σ, Ω_1) = (var"##MTKArg#459"[1], var"##MTKArg#459"[2], var"##MTKArg#459"[3], var"##MTKArg#460"[1], var"##MTKArg#460"[2], var"##MTKArg#461"[1], var"##MTKArg#461"[2], var"##MTKArg#461"[3], var"##MTKArg#462"[1], var"##MTKArg#462"[2], var"##MTKArg#462"[3])
                    (var"##MTIIPVar#465"[1])[1] = (+)((*)(-1, (exp)(z), (^)((inv)(c), 1), α, β, (^)(k, (+)(-2, α))), (*)(-1, (exp)(z), (^)((inv)(c), 1), β, (^)(k, (+)(-2, α)), (+)(-1, α)), (*)(-1, (log)(k), (exp)(z), (^)((inv)(c), 1), α, (+)(-1, α), β, (^)(k, (+)(-2, α))))
                    (var"##MTIIPVar#465"[1])[6] = (+)((*)(-1, (exp)(z), (^)((inv)(c), 1), β, (^)(k, (+)(-1, α))), (*)(-1, (log)(k), (exp)(z), (^)((inv)(c), 1), α, β, (^)(k, (+)(-1, α))))
                    (var"##MTIIPVar#465"[2])[1] = (*)(-1, (exp)(z), (^)((inv)(c), 1), α, (^)(k, (+)(-2, α)), (+)(-1, α))
                    (var"##MTIIPVar#465"[2])[6] = (*)(-1, (exp)(z), (^)((inv)(c), 1), α, (^)(k, (+)(-1, α)))
                end
            end
        nothing
    end
const H_x_p! = (var"##MTIIPVar#474", var"##MTKArg#468", var"##MTKArg#469", var"##MTKArg#470", var"##MTKArg#471", var"##MTKArg#472")->begin
        @inbounds begin
                let (c, q, i, k, z, α, β, ρ, δ, σ, Ω_1) = (var"##MTKArg#468"[1], var"##MTKArg#468"[2], var"##MTKArg#468"[3], var"##MTKArg#469"[1], var"##MTKArg#469"[2], var"##MTKArg#470"[1], var"##MTKArg#470"[2], var"##MTKArg#470"[3], var"##MTKArg#471"[1], var"##MTKArg#471"[2], var"##MTKArg#471"[3])
                    (var"##MTIIPVar#474"[1])[3] = (+)((*)(-1, (exp)(z), (^)(k, (+)(-1, α))), (*)(-1, (log)(k), (exp)(z), α, (^)(k, (+)(-1, α))))
                    (var"##MTIIPVar#474"[1])[8] = (*)(-1, (log)(k), (exp)(z), (^)(k, α))
                    (var"##MTIIPVar#474"[3])[9] = -1
                end
            end
        nothing
    end
const H_p! = (var"##MTIIPVar#483", var"##MTKArg#477", var"##MTKArg#478", var"##MTKArg#479", var"##MTKArg#480", var"##MTKArg#481")->begin
        @inbounds begin
                let (c, q, i, k, z, α, β, ρ, δ, σ, Ω_1) = (var"##MTKArg#477"[1], var"##MTKArg#477"[2], var"##MTKArg#477"[3], var"##MTKArg#478"[1], var"##MTKArg#478"[2], var"##MTKArg#479"[1], var"##MTKArg#479"[2], var"##MTKArg#479"[3], var"##MTKArg#480"[1], var"##MTKArg#480"[2], var"##MTKArg#480"[3])
                    var"##MTIIPVar#483"[1] = (*)(-1, (^)((inv)(c), 1), β, (+)((*)((exp)(z), (^)(k, (+)(-1, α))), (*)((log)(k), (exp)(z), α, (^)(k, (+)(-1, α)))))
                    var"##MTIIPVar#483"[3] = (*)(-1, (log)(k), (exp)(z), (^)(k, α))
                    var"##MTIIPVar#483"[6] = (*)(-1, (+)(1, (*)(-1, δ), (*)((exp)(z), α, (^)(k, (+)(-1, α)))), (^)((inv)(c), 1))
                    var"##MTIIPVar#483"[14] = (*)(-1, z)
                end
            end
        nothing
    end
const Ψ! = (var"##MTIIPVar#492", var"##MTKArg#486", var"##MTKArg#487", var"##MTKArg#488", var"##MTKArg#489", var"##MTKArg#490")->begin
        @inbounds begin
                let (c, q, i, k, z, α, β, ρ, δ, σ, Ω_1) = (var"##MTKArg#486"[1], var"##MTKArg#486"[2], var"##MTKArg#486"[3], var"##MTKArg#487"[1], var"##MTKArg#487"[2], var"##MTKArg#488"[1], var"##MTKArg#488"[2], var"##MTKArg#488"[3], var"##MTKArg#489"[1], var"##MTKArg#489"[2], var"##MTKArg#489"[3])
                    (var"##MTIIPVar#492"[1])[1] = (*)(-2, (+)(1, (*)(-1, δ), (*)((exp)(z), α, (^)(k, (+)(-1, α)))), (^)((inv)(c), 3), β)
                    (var"##MTIIPVar#492"[1])[7] = (*)((exp)(z), (^)((inv)(c), 2), α, β, (^)(k, (+)(-2, α)), (+)(-1, α))
                    (var"##MTIIPVar#492"[1])[8] = (*)((exp)(z), (^)((inv)(c), 2), α, β, (^)(k, (+)(-1, α)))
                    (var"##MTIIPVar#492"[1])[34] = (*)(2, (^)((inv)(c), 3))
                    (var"##MTIIPVar#492"[1])[61] = (*)((exp)(z), (^)((inv)(c), 2), α, β, (^)(k, (+)(-2, α)), (+)(-1, α))
                    (var"##MTIIPVar#492"[1])[67] = (*)(-1, (exp)(z), (^)((inv)(c), 1), α, (+)(-1, α), β, (^)(k, (+)(-3, α)), (+)(-2, α))
                    (var"##MTIIPVar#492"[1])[68] = (*)(-1, (exp)(z), (^)((inv)(c), 1), α, β, (^)(k, (+)(-2, α)), (+)(-1, α))
                    (var"##MTIIPVar#492"[1])[71] = (*)((exp)(z), (^)((inv)(c), 2), α, β, (^)(k, (+)(-1, α)))
                    (var"##MTIIPVar#492"[1])[77] = (*)(-1, (exp)(z), (^)((inv)(c), 1), α, β, (^)(k, (+)(-2, α)), (+)(-1, α))
                    (var"##MTIIPVar#492"[1])[78] = (*)(-1, (exp)(z), (^)((inv)(c), 1), α, β, (^)(k, (+)(-1, α)))
                    (var"##MTIIPVar#492"[3])[89] = (*)(-1, (exp)(z), α, (^)(k, (+)(-2, α)), (+)(-1, α))
                    (var"##MTIIPVar#492"[3])[90] = (*)(-1, (exp)(z), α, (^)(k, (+)(-1, α)))
                    (var"##MTIIPVar#492"[3])[99] = (*)(-1, (exp)(z), α, (^)(k, (+)(-1, α)))
                    (var"##MTIIPVar#492"[3])[100] = (*)(-1, (exp)(z), (^)(k, α))
                end
            end
        nothing
    end
const H̄! = (var"##MTIIPVar#500", var"##MTKArg#495", var"##MTKArg#496", var"##MTKArg#497", var"##MTKArg#498")->begin
        @inbounds begin
                let (c, q, i, k, z, α, β, ρ, δ, σ, Ω_1) = (var"##MTKArg#495"[1], var"##MTKArg#495"[2], var"##MTKArg#495"[3], var"##MTKArg#495"[4], var"##MTKArg#495"[5], var"##MTKArg#496"[1], var"##MTKArg#496"[2], var"##MTKArg#496"[3], var"##MTKArg#497"[1], var"##MTKArg#497"[2], var"##MTKArg#497"[3])
                    var"##MTIIPVar#500"[1] = (+)((^)((inv)(c), 1), (*)(-1, (+)(1, (*)(-1, δ), (*)((exp)(z), α, (^)(k, (+)(-1, α)))), (^)((inv)(c), 1), β))
                    var"##MTIIPVar#500"[2] = (+)(c, k, (*)(-1, (+)(q, (*)(k, (+)(1, (*)(-1, δ))))))
                    var"##MTIIPVar#500"[3] = (+)((*)(-1, (exp)(z), (^)(k, α)), q)
                    var"##MTIIPVar#500"[4] = (+)((*)(-1, z, ρ), z)
                    var"##MTIIPVar#500"[5] = (+)(i, (*)(-1, (+)((*)(-1, k, (+)(1, (*)(-1, δ))), k)))
                end
            end
        nothing
    end
const H̄_w! = (var"##MTIIPVar#508", var"##MTKArg#503", var"##MTKArg#504", var"##MTKArg#505", var"##MTKArg#506")->begin
        @inbounds begin
                let (c, q, i, k, z, α, β, ρ, δ, σ, Ω_1) = (var"##MTKArg#503"[1], var"##MTKArg#503"[2], var"##MTKArg#503"[3], var"##MTKArg#503"[4], var"##MTKArg#503"[5], var"##MTKArg#504"[1], var"##MTKArg#504"[2], var"##MTKArg#504"[3], var"##MTKArg#505"[1], var"##MTKArg#505"[2], var"##MTKArg#505"[3])
                    var"##MTIIPVar#508"[1] = (+)((*)(-1, (^)((inv)(c), 2)), (*)((+)(1, (*)(-1, δ), (*)((exp)(z), α, (^)(k, (+)(-1, α)))), (^)((inv)(c), 2), β))
                    var"##MTIIPVar#508"[2] = 1
                    var"##MTIIPVar#508"[7] = -1
                    var"##MTIIPVar#508"[8] = 1
                    var"##MTIIPVar#508"[15] = 1
                    var"##MTIIPVar#508"[16] = (*)((*)((*)((*)(-1, (exp)(z)), (^)((inv)(c), 1)), α), (*)((*)(β, (^)(k, (+)(-2, α))), (+)(-1, α)))
                    var"##MTIIPVar#508"[17] = (+)(1, (*)(-1, (+)(1, (*)(-1, δ))))
                    var"##MTIIPVar#508"[18] = (*)(-1, (exp)(z), α, (^)(k, (+)(-1, α)))
                    var"##MTIIPVar#508"[20] = (*)(-1, (+)(1, (*)(-1, (+)(1, (*)(-1, δ)))))
                    var"##MTIIPVar#508"[21] = (*)((*)((*)((*)(-1, (exp)(z)), (^)((inv)(c), 1)), α), (*)(β, (^)(k, (+)(-1, α))))
                    var"##MTIIPVar#508"[23] = (*)(-1, (exp)(z), (^)(k, α))
                    var"##MTIIPVar#508"[24] = (+)(1, (*)(-1, ρ))
                end
            end
        nothing
    end
const ȳ_iv! = (var"##MTIIPVar#515", var"##MTKArg#511", var"##MTKArg#512", var"##MTKArg#513")->begin
        @inbounds begin
                let (α, β, ρ, δ, σ, Ω_1) = (var"##MTKArg#511"[1], var"##MTKArg#511"[2], var"##MTKArg#511"[3], var"##MTKArg#512"[1], var"##MTKArg#512"[2], var"##MTKArg#512"[3])
                    var"##MTIIPVar#515"[1] = (-)((^)((/)((+)((-)((/)(1, β), 1), δ), α), (/)(α, (-)(α, 1))), (*)(δ, (^)((/)((+)((-)((/)(1, β), 1), δ), α), (/)(1, (-)(α, 1)))))
                    var"##MTIIPVar#515"[2] = (^)((/)((+)((-)((/)(1, β), 1), δ), α), (/)(α, (-)(α, 1)))
                    var"##MTIIPVar#515"[3] = (*)(δ, (^)((/)((+)((-)((/)(1, β), 1), δ), α), (/)(1, (-)(α, 1))))
                end
            end
        nothing
    end
const x̄_iv! = (var"##MTIIPVar#522", var"##MTKArg#518", var"##MTKArg#519", var"##MTKArg#520")->begin
        @inbounds begin
                let (α, β, ρ, δ, σ, Ω_1) = (var"##MTKArg#518"[1], var"##MTKArg#518"[2], var"##MTKArg#518"[3], var"##MTKArg#519"[1], var"##MTKArg#519"[2], var"##MTKArg#519"[3])
                    var"##MTIIPVar#522"[1] = (^)((/)((+)((-)((/)(1, β), 1), δ), α), (/)(1, (-)(α, 1)))
                end
            end
        nothing
    end
const ȳ! = (var"##MTIIPVar#529", var"##MTKArg#525", var"##MTKArg#526", var"##MTKArg#527")->begin
        @inbounds begin
                let (α, β, ρ, δ, σ, Ω_1) = (var"##MTKArg#525"[1], var"##MTKArg#525"[2], var"##MTKArg#525"[3], var"##MTKArg#526"[1], var"##MTKArg#526"[2], var"##MTKArg#526"[3])
                    var"##MTIIPVar#529"[1] = (-)((^)((/)((+)((-)((/)(1, β), 1), δ), α), (/)(α, (-)(α, 1))), (*)(δ, (^)((/)((+)((-)((/)(1, β), 1), δ), α), (/)(1, (-)(α, 1)))))
                    var"##MTIIPVar#529"[2] = (^)((/)((+)((-)((/)(1, β), 1), δ), α), (/)(α, (-)(α, 1)))
                    var"##MTIIPVar#529"[3] = (*)(δ, (^)((/)((+)((-)((/)(1, β), 1), δ), α), (/)(1, (-)(α, 1))))
                end
            end
        nothing
    end
const x̄! = (var"##MTIIPVar#536", var"##MTKArg#532", var"##MTKArg#533", var"##MTKArg#534")->begin
        @inbounds begin
                let (α, β, ρ, δ, σ, Ω_1) = (var"##MTKArg#532"[1], var"##MTKArg#532"[2], var"##MTKArg#532"[3], var"##MTKArg#533"[1], var"##MTKArg#533"[2], var"##MTKArg#533"[3])
                    var"##MTIIPVar#536"[1] = (^)((/)((+)((-)((/)(1, β), 1), δ), α), (/)(1, (-)(α, 1)))
                end
            end
        nothing
    end
const ȳ_p! = (var"##MTIIPVar#543", var"##MTKArg#539", var"##MTKArg#540", var"##MTKArg#541")->begin
        @inbounds begin
                let (α, β, ρ, δ, σ, Ω_1) = (var"##MTKArg#539"[1], var"##MTKArg#539"[2], var"##MTKArg#539"[3], var"##MTKArg#540"[1], var"##MTKArg#540"[2], var"##MTKArg#540"[3])
                    var"##MTIIPVar#543"[1] = (+)((*)(-1, δ, (+)((*)(-1, (log)((*)((+)(-1, (^)((inv)(β), 1), δ), (^)((inv)(α), 1))), (^)((*)((+)(-1, (^)((inv)(β), 1), δ), (^)((inv)(α), 1)), (^)((inv)((+)(-1, α)), 1)), (^)((inv)((+)(-1, α)), 2)), (*)((*)((*)((*)(-1, (+)(-1, (^)((inv)(β), 1), δ)), (^)((inv)(α), 2)), (^)((*)((+)(-1, (^)((inv)(β), 1), δ), (^)((inv)(α), 1)), (+)(-1, (^)((inv)((+)(-1, α)), 1)))), (^)((inv)((+)(-1, α)), 1)))), (*)((log)((*)((+)(-1, (^)((inv)(β), 1), δ), (^)((inv)(α), 1))), (+)((^)((inv)((+)(-1, α)), 1), (*)(-1, α, (^)((inv)((+)(-1, α)), 2))), (^)((*)((+)(-1, (^)((inv)(β), 1), δ), (^)((inv)(α), 1)), (*)(α, (^)((inv)((+)(-1, α)), 1)))), (*)((*)((*)((*)(-1, (+)(-1, (^)((inv)(β), 1), δ)), (^)((inv)(α), 1)), (^)((*)((+)(-1, (^)((inv)(β), 1), δ), (^)((inv)(α), 1)), (+)(-1, (*)(α, (^)((inv)((+)(-1, α)), 1))))), (^)((inv)((+)(-1, α)), 1)))
                    var"##MTIIPVar#543"[2] = (+)((*)((log)((*)((+)(-1, (^)((inv)(β), 1), δ), (^)((inv)(α), 1))), (+)((^)((inv)((+)(-1, α)), 1), (*)(-1, α, (^)((inv)((+)(-1, α)), 2))), (^)((*)((+)(-1, (^)((inv)(β), 1), δ), (^)((inv)(α), 1)), (*)(α, (^)((inv)((+)(-1, α)), 1)))), (*)((*)((*)((*)(-1, (+)(-1, (^)((inv)(β), 1), δ)), (^)((inv)(α), 1)), (^)((*)((+)(-1, (^)((inv)(β), 1), δ), (^)((inv)(α), 1)), (+)(-1, (*)(α, (^)((inv)((+)(-1, α)), 1))))), (^)((inv)((+)(-1, α)), 1)))
                    var"##MTIIPVar#543"[3] = (*)(δ, (+)((*)(-1, (log)((*)((+)(-1, (^)((inv)(β), 1), δ), (^)((inv)(α), 1))), (^)((*)((+)(-1, (^)((inv)(β), 1), δ), (^)((inv)(α), 1)), (^)((inv)((+)(-1, α)), 1)), (^)((inv)((+)(-1, α)), 2)), (*)((*)((*)((*)(-1, (+)(-1, (^)((inv)(β), 1), δ)), (^)((inv)(α), 2)), (^)((*)((+)(-1, (^)((inv)(β), 1), δ), (^)((inv)(α), 1)), (+)(-1, (^)((inv)((+)(-1, α)), 1)))), (^)((inv)((+)(-1, α)), 1))))
                    var"##MTIIPVar#543"[4] = (+)((*)(-1, (^)((*)((+)(-1, (^)((inv)(β), 1), δ), (^)((inv)(α), 1)), (+)(-1, (*)(α, (^)((inv)((+)(-1, α)), 1)))), (^)((inv)((+)(-1, α)), 1), (^)((inv)(β), 2)), (*)((*)((*)((*)((^)((inv)(α), 1), δ), (^)((*)((+)(-1, (^)((inv)(β), 1), δ), (^)((inv)(α), 1)), (+)(-1, (^)((inv)((+)(-1, α)), 1)))), (^)((inv)((+)(-1, α)), 1)), (^)((inv)(β), 2)))
                    var"##MTIIPVar#543"[5] = (*)(-1, (^)((*)((+)(-1, (^)((inv)(β), 1), δ), (^)((inv)(α), 1)), (+)(-1, (*)(α, (^)((inv)((+)(-1, α)), 1)))), (^)((inv)((+)(-1, α)), 1), (^)((inv)(β), 2))
                    var"##MTIIPVar#543"[6] = (*)((*)((*)((*)(-1, (^)((inv)(α), 1)), δ), (^)((*)((+)(-1, (^)((inv)(β), 1), δ), (^)((inv)(α), 1)), (+)(-1, (^)((inv)((+)(-1, α)), 1)))), (*)((^)((inv)((+)(-1, α)), 1), (^)((inv)(β), 2)))
                end
            end
        nothing
    end
const x̄_p! = (var"##MTIIPVar#550", var"##MTKArg#546", var"##MTKArg#547", var"##MTKArg#548")->begin
        @inbounds begin
                let (α, β, ρ, δ, σ, Ω_1) = (var"##MTKArg#546"[1], var"##MTKArg#546"[2], var"##MTKArg#546"[3], var"##MTKArg#547"[1], var"##MTKArg#547"[2], var"##MTKArg#547"[3])
                    var"##MTIIPVar#550"[1] = (+)((*)(-1, (log)((*)((+)(-1, (^)((inv)(β), 1), δ), (^)((inv)(α), 1))), (^)((*)((+)(-1, (^)((inv)(β), 1), δ), (^)((inv)(α), 1)), (^)((inv)((+)(-1, α)), 1)), (^)((inv)((+)(-1, α)), 2)), (*)((*)((*)((*)(-1, (+)(-1, (^)((inv)(β), 1), δ)), (^)((inv)(α), 2)), (^)((*)((+)(-1, (^)((inv)(β), 1), δ), (^)((inv)(α), 1)), (+)(-1, (^)((inv)((+)(-1, α)), 1)))), (^)((inv)((+)(-1, α)), 1)))
                    var"##MTIIPVar#550"[3] = (*)((*)((*)((*)(-1, (^)((inv)(α), 1)), (^)((*)((+)(-1, (^)((inv)(β), 1), δ), (^)((inv)(α), 1)), (+)(-1, (^)((inv)((+)(-1, α)), 1)))), (^)((inv)((+)(-1, α)), 1)), (^)((inv)(β), 2))
                end
            end
        nothing
    end
const steady_state! = nothing
const allocate_solver_cache = m->begin
        #= C:\Google Drive\Projects\HMCDSGE\src\mtk_constructor.jl:142 =#
        FirstOrderSolverCache(m)
    end
const select_p_ss_hash = LinearAlgebra.UniformScaling{Bool}(true)
const select_p_f_ss_hash = LinearAlgebra.UniformScaling{Bool}(true)
const select_p_perturbation_hash = LinearAlgebra.UniformScaling{Bool}(true)
const select_p_f_perturbation_hash = LinearAlgebra.UniformScaling{Bool}(true)
end


text:
include("../imports.jl")

@testset ExtendedTestSet "TabularModel Components" begin
    @testset ExtendedTestSet "embeddingbackbone" begin
        embed_szs = [(5, 10), (100, 30), (2, 30)]
        embeds = FastAI.Models.tabular_embedding_backbone(embed_szs, 0.)
        x = [rand(1:n) for (n, _) in embed_szs]

        @test size(embeds(x)) == (70, 1)
    end

    @testset ExtendedTestSet "continuousbackbone" begin
        n = 5
        contback = FastAI.Models.tabular_continuous_backbone(n)
        x = rand(5, 1)
        @test size(contback(x)) == (5, 1)
    end

    @testset ExtendedTestSet "TabularModel" begin
        n = 5
        embed_szs = [(5, 10), (100, 30), (2, 30)]

        embeds = FastAI.Models.tabular_embedding_backbone(embed_szs, 0.)
        contback = FastAI.Models.tabular_continuous_backbone(n)

        x = ([rand(1:n) for (n, _) in embed_szs], rand(5, 1))

        tm = TabularModel(embeds, contback; outsize=4)
        @test size(tm(x)) == (4, 1)

        tm2 = TabularModel(embeds, contback, Chain(Dense(100, 4), x->FastAI.Models.sigmoidrange(x, 2, 5)))
        y2 = tm2(x)
        @test all(y2.> 2) && all(y2.<5)

        cardinalities = [4, 99, 1]
        tm3 = TabularModel(n, 4, [200, 100], cardinalities = cardinalities, size_overrides = (10, 30, 30))
        @test size(tm3(x)) == (4, 1)
    end
end

ext: jl


text:
module MicroBenchmarks

using BenchmarkTools
using Compat
if VERSION >= v"0.7.0-DEV.3449"
    using LinearAlgebra
end
if VERSION >= v"0.7.0-beta.85"
    using Statistics
end

# This module contains the Julia microbenchmarks shown in the language
# comparison table at http://julialang.org/.

include(joinpath(dirname(@__FILE__), "..", "utils", "RandUtils.jl"))
using .RandUtils

include("methods.jl")

const SUITE = BenchmarkGroup(["recursion", "fibonacci", "fib",  "parse", "parseint",
                              "mandel", "mandelbrot", "sort", "quicksort", "pi", "π", "sum",
                              "pisum", "πsum", "rand", "randmatstat", "rand", "randmatmul"])

SUITE["fib"] = @benchmarkable perf_micro_fib(20)
SUITE["parseint"] = @benchmarkable perf_micro_parseint(1000)
SUITE["mandel"] = @benchmarkable perf_micro_mandel()
SUITE["quicksort"] = @benchmarkable perf_micro_quicksort(5000)
SUITE["pisum"] = @benchmarkable perf_micro_pisum()
SUITE["randmatstat"] = @benchmarkable perf_micro_randmatstat(1000)
SUITE["randmatmul"] = @benchmarkable perf_micro_randmatmul(1000)
if (VERSION <= v"0.7.0-DEV.914" ? is_unix() : Sys.isunix())
    SUITE["printfd"] = @benchmarkable perf_printfd(10000)
end

end # module


text:
"""
Data structure and resolution methods of some forms of bilevel flow problems
"""
module BilevelFlowProblems

using SparseArrays: spzeros
using JuMP: @variable, @constraint, @objective
import JuMP

import ..MT
import ..VT

export BilevelFlowProblem

import ..build_blp_model
import ..SOS1Complementarity

"""
    BilevelFlowProblem{M1<:MT,M2<:AbstractMatrix{Bool},M3<:MT,M4<:AbstractArray{<:Real,3}}

Holding the information of a bilevel flow problem:
* Upper-level picks arcs to tax and a discrete tax level
* Lower-level solves a mincost flow problem with `minflow`
units out of the source and arc cost (init_cost + tax_options * chosen_option)
"""
struct BilevelFlowProblem{M1<:MT,M2<:AbstractMatrix{Bool},M3<:MT,M4<:AbstractArray{<:Real,3}}
    init_cost::M1
    taxable_edges::M2
    capacities::M3
    tax_options::M4
    nv::Int
    ne::Int
    minflow::Float64
    function BilevelFlowProblem(init_cost::M1,taxable_edges::M2,capacities::M3,tax_options::M4, minflow) where {M1<:MT,M2<:AbstractMatrix{Bool},M3<:MT,M4<:AbstractArray{<:Real,3}}
        size(init_cost) == size(taxable_edges) == size(capacities) == size(tax_options)[1:end-1] || DimensionMismatch("Matrix dimensions")
        nv = size(init_cost)[1]
        nv == size(init_cost)[2] || DimensionMismatch("Cost matrix should be square")
        ne = sum(1 for i in Base.OneTo(nv) for j in Base.OneTo(nv) if capacities[i,j] > 0. && i!=j)
        new{M1,M2,M3,M4}(init_cost, taxable_edges, capacities, tax_options, nv, ne, minflow)
    end
end

"""
	`build_blp_model(bfp::BilevelFlowProblem, solver; comp_method)`

Build the JuMP model from the data of `bfp` and assign it the passed solver.
Return `(m, r, y, f, λ, s)`
- `m`: `JuMP.Model`
- `y[i,j,k]`: binary variables indicating which option tax option is chosen by the upper-level
- `r[i,j]`: auxiliary variable storing the revenue made on arc `[i,j]`
- `f[i,j]`: flow variables (lower level)
- `λ`: flattened dual vector
- `s`: flattened slack variables for the flow problem
"""
function build_blp_model(bfp::BilevelFlowProblem, solver; comp_method = SOS1Complementarity())
    m = JuMP.Model(solver = solver)
    (nv, _, nopt) = size(bfp.tax_options)
    @variable(m, y[i=1:nv,j=1:nv,k=1:nopt], Bin)
    @variable(m, r[i=1:nv,j=1:nv] >= 0.)
    @variable(m, f[i=1:nv,j=1:nv] >= 0.)
    @constraint(m, revenue_limit[i=1:nv,j=1:nv,k=1:nopt], r[i,j] <= bfp.tax_options[i,j,k] * f[i,j] + maximum(bfp.tax_options[i,j,:]) * bfp.capacities[i,j]*(1-y[i,j,k]))
    @constraint(m, unique_opt[i=1:nv,j=1:nv], sum(y[i,j,k] for k in 1:nopt) == 1.)
    @objective(m, Max,
        sum(r[i,j] for i in 1:nv for j in 1:nv if bfp.taxable_edges[i,j])
    )
    flat_flow = [f[i] for i in eachindex(f)]
    lin_cost  = [bfp.init_cost[i,j] + sum(y[i,j,k] * bfp.tax_options[i,j,k] for k in Base.OneTo(nopt)) for j in Base.OneTo(bfp.nv) for i in Base.OneTo(bfp.nv)]
    (B, b) = flow_constraint_standard(bfp)
    @variable(m,
        s[i=1:length(b)] >= 0.
    )
    @constraint(m, B*flat_flow .+ s .== b)
    (_, λ, _) = build_blp_model(m, B, lin_cost, s, comp_method = comp_method)
    return (m, r, y, f, λ, s)
end

"""
    flow_constraint_standard(bfp::BilevelFlowProblem)

Construct the constraint matrix `B` and right-hand side vector `b`
for the lower-level problem
"""
function flow_constraint_standard(bfp::BilevelFlowProblem)
    B = spzeros( # constraint matrix
        1 +        # min flow from source
        bfp.nv-2 + # flow conservation at all nodes but source & sink
        bfp.nv*bfp.nv,    # capacity at all edges
        bfp.nv*bfp.nv
    )
    b = zeros(size(B)[1])
    b[1] = -bfp.minflow
    if bfp.capacities[1,bfp.nv] > 0.
        B[1,end-bfp.nv+1] = -1
    end
    for k in 2:bfp.nv-1
        if bfp.capacities[1,k] > 0.
            B[1,bfp.nv*(k-1)+1] = -1. # arc s -> k
        end
        for i in Base.OneTo(bfp.nv)
            if i != k
                if bfp.capacities[i,k] > 0.
                    B[k,bfp.nv*(k-1)+i] = 1.
                end
                if bfp.capacities[k,i] > 0.
                    B[k,bfp.nv*(i-1)+k] = -1.
                end
            end
        end
    end
    for i in Base.OneTo(bfp.nv), j in Base.OneTo(bfp.nv)
        B[bfp.nv-1+i+bfp.nv*(j-1),i+bfp.nv*(j-1)] = 1.
        b[bfp.nv-1+i+bfp.nv*(j-1)] = bfp.capacities[i,j]
    end
    return (B, b)
end

end


text:
function set_bond_values_reactive(; session::ServerSession, notebook::Notebook, bound_sym_names::AbstractVector{Symbol}, is_first_value=nothing, kwargs...)::Union{Task,TopologicalOrder}
    if is_first_value !== nothing
        @warn "is_first_value is deprecated, you don't need to set it anymore: https://github.com/fonsp/Pluto.jl/pull/975"
    end
    # filter out the bonds that don't need to be set
    to_set = filter(bound_sym_names) do bound_sym
        new_value = notebook.bonds[bound_sym].value
        is_first_value = notebook.bonds[bound_sym].is_first_value

        variable_exists = is_assigned_anywhere(notebook, notebook.topology, bound_sym)
        if !variable_exists
            # a bond was set while the cell is in limbo state
            # we don't need to do anything
            return false
        end

        # TODO: Not checking for any dependents now
        # any_dependents = is_referenced_anywhere(notebook, notebook.topology, bound_sym)

        # fix for https://github.com/fonsp/Pluto.jl/issues/275
        # if `Base.get` was defined to give an initial value (read more about this in the Interactivity sample notebook), then we want to skip the first value sent back from the bond. (if `Base.get` was not defined, then the variable has value `missing`)
        # Check if the variable does not already have that value.
        # because if the initial value is already set, then we don't want to run dependent cells again.
        eq_tester = :(try !ismissing($bound_sym) && ($bound_sym == $new_value) catch; false end) # not just a === comparison because JS might send back the same value but with a different type (Float64 becomes Int64 in JS when it's an integer.)
        if is_first_value && WorkspaceManager.eval_fetch_in_workspace((session, notebook), eq_tester)
            return false
        end
        return true
    end

    if isempty(to_set)
        return TopologicalOrder(notebook.topology, Cell[], Dict{Cell, ReactivityError}())
    end

    new_values = [notebook.bonds[bound_sym].value for bound_sym in to_set]

    function custom_deletion_hook((session, notebook)::Tuple{ServerSession,Notebook}, old_workspace_name, new_workspace_name, to_delete_vars::Set{Symbol}, methods_to_delete::Set{Tuple{UUID,FunctionName}}, to_reimport::Set{Expr}; to_run::AbstractVector{Cell})
        to_delete_vars = Set([to_delete_vars..., to_set...]) # also delete the bound symbols
        WorkspaceManager.move_vars((session, notebook), old_workspace_name, new_workspace_name, to_delete_vars, methods_to_delete, to_reimport)
        for (bound_sym, new_value) in zip(to_set, new_values)
            WorkspaceManager.eval_in_workspace((session, notebook), :($(bound_sym) = $(new_value)))
        end
    end
    to_reeval = where_referenced(notebook, notebook.topology, Set{Symbol}(to_set))

    run_reactive_async!(session, notebook, to_reeval; deletion_hook=custom_deletion_hook, persist_js_state=true, run_async=false, kwargs...)
end


text:
export PR, auc

struct PR{T <: Real} <: AbstractBinaryMetric
  labels::BitVector
  scores::Vector{T}
  AUC::T
  precision::Vector{T}
  recall::Vector{T}
  OP::Vector{T}
  function PR(labels,scores)
    B = BinaryMetrics(labels, scores)
    OP = op(B)
    p = precision(B, OP)
    r = recall(B,OP)
    AUC = get_auc([0.0;r], [1.0;p])
    new{eltype(p)}(B.labels, B.scores, AUC, p, r, OP)
  end
end

function Base.show(io::IO, B::PR)
  println(io,"Precision-Recall (PR) curve")
  AUC = round(auc(B),digits=2)
  print(io,"Samples:$(length(B)), P:$(positives(B)), N:$(negatives(B)), AUC:$(AUC)")
end

"""
  PR(labels, scores)

Generates data for plotting a precision-recall (PR) curve.

* `labels` must be a boolean array e.g. either containing `true`,`false` or `0`,`1`.
* `scores` must be an array with the same length of `labels` containing the scores of the binary classifier.

You can generate a plot automatically using the following commands:
```julia
pr = PR(labels, scores)

using Plots
plot(pr)
```

The Area Under Curve (AUC) FPR and TPR can be accessed by `auc(pr)`, `fpr(pr)` and `tpr(pr)` respectively.

"""
PR(B::BinaryMetrics) = PR(B.labels, B.scores)

auc(curve::PR) = curve.AUC
ppv(curve::PR) = curve.precision
tpr(curve::PR) = curve.recall
op(curve::PR) = curve.OP

function get_curve(curve::PR)
  return [0.0;curve.recall], [1.0;curve.precision]
end


text:
using SafeTestsets, Test

@safetestset "Parsing Test" begin include("variable_parsing.jl") end
@safetestset "Differentiation Test" begin include("derivatives.jl") end
@safetestset "Simplify Test" begin include("simplify.jl") end
@safetestset "Operation Overloads Test" begin include("operation_overloads.jl") end
@safetestset "Direct Usage Test" begin include("direct.jl") end
@safetestset "Build Function Test" begin include("build_function.jl") end
@safetestset "ODESystem Test" begin include("odesystem.jl") end
@safetestset "Mass Matrix Test" begin include("mass_matrix.jl") end
@safetestset "SDESystem Test" begin include("sdesystem.jl") end
@safetestset "NonlinearSystem Test" begin include("nonlinearsystem.jl") end
@safetestset "OptimizationSystem Test" begin include("optimizationsystem.jl") end
@safetestset "ReactionSystem Test" begin include("reactionsystem.jl") end
@safetestset "Build Targets Test" begin include("build_targets.jl") end
@safetestset "Domain Test" begin include("domains.jl") end
@safetestset "Constraints Test" begin include("constraints.jl") end
@safetestset "PDE Construction Test" begin include("pde.jl") end
@safetestset "Lowering Integration Test" begin include("lowering_solving.jl") end
@safetestset "Test Big System Usage" begin include("bigsystem.jl") end
@testset "Latexify recipes Test" begin include("latexify.jl") end
@testset "Distributed Test" begin include("distributed.jl") end


text:
struct Wild{TAG} end

tagof(x::Sym{Wild}) = DEFAULT_TAG
tagof(x::Sym{Wild{TAG}}) where {TAG} = TAG

iswild(x) = false
iswild(x::Sym{<:Wild}) = true


text:
"""
    MortalityTable(CSV.File)

Read and parse a CSV file in the SOA's usual CSV format. You must import and use CSV.jl before calling this function.

# Examples

```julia-repl
julia> path = "path/to/table.csv"
julia> file = CSV.File(path,header=false) # no real header in the file
julia> MortalityTable(file)
MortalityTable (Insured Lives Mortality):
   Name:
       2015 VBT Female Non-Smoker RR50 ALB
   Fields:
       (:select, :ultimate, :metadata)
   Provider:
       American Academy of Actuaries along with the Society of Actuaries
   mort.SOA.org ID:
       3209
   mort.SOA.org link:
       https://mort.soa.org/ViewTable.aspx?&TableIdentity=3209
   Description:
       2015 VBT Relative Risk Table - Female, 50% Non-Smoker, Age Last Birthday, Select
```

And in one line:
```julia-repl
julia> MortalityTables.MortalityTable(CSV.File( "path/to/table.csv",header=false))
MortalityTable (Insured Lives Mortality):
   Name:
       2015 VBT Female Non-Smoker RR50 ALB
   Fields:
       (:select, :ultimate, :metadata)
   Provider:
       American Academy of Actuaries along with the Society of Actuaries
   mort.SOA.org ID:
       3209
   mort.SOA.org link:
       https://mort.soa.org/ViewTable.aspx?&TableIdentity=3209
   Description:
       2015 VBT Relative Risk Table - Female, 50% Non-Smoker, Age Last Birthday, Select
```

"""
function MortalityTable(lines::CSV.File)
	#what lines the table starts at
	table_starts = findall(line -> ~ismissing(line[1]) && line[1] == "Row\\Column",lines) .+ 1

	# Construct MetaData
	raw_meta = Dict()
	for line in lines
		if ismissing(line[1])
			break
		end
		raw_meta[line[1]] = line[2]
	end

	d = TableMetaData(
		name = get(raw_meta,"Table Name:",nothing),
		id = get(raw_meta,"Table Identity:",nothing),
		provider = get(raw_meta,"Provider Name:",nothing),
		reference = get(raw_meta,"Table Reference:",nothing),
		content_type = get(raw_meta,"Content Type:",nothing),
		description = get(raw_meta,"Table Description:",nothing),
		comments = get(raw_meta,"Comments:",nothing),
	)

# 	scale = get(raw_meta,"Scaling Factor:",nothing)

	# Extract values

	# figure out where table ends
	table_ends= [last_values_line(lines,ts) for ts in table_starts]

	# Parse into table

	if length(table_starts) == 1
		# ultimate only
		ult_start, ult_end = table_starts[1],table_ends[1]
		ult_rates = map(lines[ult_start:ult_end]) do row
			Parsers.parse(Float64,row[2])
		end

		ult = UltimateMortality(ult_rates,start_age=Parsers.parse(Int,lines[ult_start][1]))

		return MortalityTable(ult; metadata=d)
	else
		# select and ultimate
		ult_start, ult_end = table_starts[2],table_ends[2]
		ult_rates = map(lines[ult_start:ult_end]) do row
			Parsers.parse(Float64,row[2])
		end

		ult = UltimateMortality(ult_rates,start_age=Parsers.parse(Int,lines[ult_start][1]))

		ult_end_age = lastindex(ult)

		sel_start, sel_end = table_starts[1],table_ends[1]

		sel_start_age = parsemaybe(Int,lines[sel_start][1])
		sel_rates = OffsetArray(
			map(sel_start:sel_end) do r
				start_age = sel_start_age + r - sel_start
				sel_end_age = length(lines[sel_start]) + start_age - 2
				while ismissing(lines[r][sel_end_age - start_age + 2])
					sel_end_age -= 1
				end

				age_range = start_age:max(ult_end_age,sel_end_age)
				rates = map(age_range) do age
					c = age - start_age + 2
					if age > sel_end_age
						return ult[age]
					else
						return parsemaybe(Float64,lines[r][c])
					end
				end

				return OffsetArray(
					rates,
					start_age - 1
				)
			end,
			sel_start_age - 1
		)

		return MortalityTable(sel_rates,ult,metadata=d)

	end

end

function last_values_line(lines,startline)
	for i in startline:lastindex(lines)
		if ismissing(lines[i][1])
			return i - 1
		end
	end
	return lastindex(lines)
end

# because of the poor standardization of the CSV formatted tables from mort.SOA.org,
# sometimes the value comes through as a string, sometimes as a number when CSV.jl parses it
parsemaybe(t,x) = typeof(x) <: AbstractString ? Parsers.parse(t,x) : x


text:
using LinearAlgebra
using Test

using QuantumOpticsBase

@testset "pauli" begin

@test_throws MethodError PauliBasis(1.4)

# Test conversion of unitary matrices to superoperators.
q2 = PauliBasis(2)
q3 = PauliBasis(3)
CZ = DenseOperator(q2, q2, diagm(0 => [1,1,1,-1]))
CZ_sop = SuperOperator(CZ)

# Test conversion of unitary matrices to superoperators.
@test diag(CZ_sop.data) ==  ComplexF64[1,1,1,-1,1,1,1,-1,1,1,1,-1,-1,-1,-1,1]
@test CZ_sop.basis_l == CZ_sop.basis_r == (q2, q2)

# Test conversion of superoperator to Pauli transfer matrix.
CZ_ptm = PauliTransferMatrix(CZ_sop)

# Test DensePauliTransferMatrix constructor.
@test_throws DimensionMismatch DensePauliTransferMatrix((q2, q2), (q3, q3), CZ_ptm.data)
@test DensePauliTransferMatrix((q2, q2), (q2, q2), CZ_ptm.data) == CZ_ptm

@test all(isapprox.(CZ_ptm.data[[1,30,47,52,72,91,117,140,166,185,205,210,227,256]], 1))
@test all(isapprox.(CZ_ptm.data[[106,151]], -1))

@test CZ_ptm == PauliTransferMatrix(ChiMatrix(CZ))

# Test construction of non-symmetric unitary.
CNOT = DenseOperator(q2, q2, diagm(0 => [1,1,0,0], 1 => [0,0,1], -1 => [0,0,1]))
CNOT_sop = SuperOperator(CNOT)
CNOT_chi = ChiMatrix(CNOT)
CNOT_ptm = PauliTransferMatrix(CNOT)

@test CNOT_sop.basis_l == CNOT_sop.basis_r == (q2, q2)
@test CNOT_chi.basis_l == CNOT_chi.basis_r == (q2, q2)
@test CNOT_ptm.basis_l == CNOT_ptm.basis_r == (q2, q2)

@test all(isapprox.(imag.(CNOT_sop.data), 0))
@test all(isapprox.(imag.(CNOT_chi.data), 0))
@test all(isapprox.(imag.(CNOT_ptm.data), 0))

@test all(isapprox.(CNOT_sop.data[[1,18,36,51,69,86,104,119,141,158,176,191,201,218,236,251]], 1))
@test all(isapprox.(CNOT_chi.data[[1,2,13,17,18,29,193,194,205,222]], 1))
@test all(isapprox.(CNOT_chi.data[[14,30,206,209,210,221]], -1))
@test all(isapprox.(CNOT_ptm.data[[1,18,47,64,70,85,108,138,153,183,205,222,227,244,]], 1))
@test all(isapprox.(CNOT_ptm.data[[123,168]], -1))

# Test DenseChiMatrix constructor.
@test_throws DimensionMismatch DenseChiMatrix((q2, q2), (q3, q3), CNOT_chi.data)
@test DenseChiMatrix((q2, q2), (q2, q2), CNOT_chi.data) == CNOT_chi

# Test equality and conversion among all three bases.
ident = Complex{Float64}[1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1]

IDENT = DenseOperator(q2, ident)

IDENT_sop = SuperOperator(IDENT)
IDENT_chi = ChiMatrix(IDENT)
IDENT_ptm = PauliTransferMatrix(IDENT)

@test ChiMatrix(IDENT_sop) == IDENT_chi
@test ChiMatrix(IDENT_ptm) == IDENT_chi
@test SuperOperator(IDENT_chi) == IDENT_sop
@test SuperOperator(IDENT_ptm) == IDENT_sop
@test PauliTransferMatrix(IDENT_sop) == IDENT_ptm
@test PauliTransferMatrix(IDENT_chi) == IDENT_ptm

# Test approximate equality and conversion among all three bases.
cphase = Complex{Float64}[1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 exp(1im*.6)]

CPHASE = DenseOperator(q2, cphase)

CPHASE_sop = SuperOperator(CPHASE)
CPHASE_chi = ChiMatrix(CPHASE)
CPHASE_ptm = PauliTransferMatrix(CPHASE)

@test isapprox(ChiMatrix(CPHASE_sop), CPHASE_chi)
@test isapprox(ChiMatrix(CPHASE_ptm), CPHASE_chi)
@test isapprox(SuperOperator(CPHASE_chi), CPHASE_sop)
@test isapprox(SuperOperator(CPHASE_ptm), CPHASE_sop)
@test isapprox(PauliTransferMatrix(CPHASE_sop), CPHASE_ptm)
@test isapprox(PauliTransferMatrix(CPHASE_chi), CPHASE_ptm)

# Test composition.
@test isapprox(ChiMatrix(CPHASE) * ChiMatrix(CNOT), ChiMatrix(CPHASE * CNOT))
@test isapprox(PauliTransferMatrix(CPHASE) * PauliTransferMatrix(CNOT), PauliTransferMatrix(CPHASE * CNOT))

end # testset


text:
"""
Structure holding an information
that `what` should be executed by scheduler at time `when`;
`what` should accept one argument of type `Scheduler`.
"""
struct Action{T<:Real}
    what::Function
    when::T
end

"""
    pq_insert!(pq, what, when)

Specialized core Julia code for insertion to priority queue `pq`
Not exported
Put `Action(what, when)` to `pq`.
Return inserted `Action`
"""
function pq_insert!(pq::Vector{Action{T}}, what::Function, when::T) where T<:Real
    p = Action(what, when)
    push!(pq, p)
    i = length(pq)
    @inbounds while (j = div(i, 2)) >= 1
        if when < pq[j].when
            pq[i] = pq[j]
            i = j
        else
            break
        end
    end
    pq[i] = p
end

"""
    pq_remove!(pq)

Specialized core Julia code for getting top element from priority queue `pq`
Not exported
Return `Action`
"""
function pq_remove!(pq::Vector{Action{T}}) where T<:Real
    x = pq[1]
    y = pop!(pq)
    if !isempty(pq)
        i = 1
        len = length(pq)
        @inbounds while (l = 2i) <= len
            r = 2i + 1
            j = r > len || pq[l].when < pq[r].when ? l : r
            if pq[j].when < y.when
                pq[i] = pq[j]
                i = j
            else
                break
            end
        end
        pq[i] = y
    end
    return x
end

"""
Abstract type for holding state of the simulation
"""
abstract type AbstractState end

"""
Simplest concrete type implementing `AbstractState`
that does not hold any data
"""
struct EmptyState <: AbstractState
end

"""
Holds information about current simulation state
Contains the following fields:
* `now`         current simulation time
* `event_queue` priority queue of `Actions` planned to be executed
* `state`       user defined subtype of `AbstractState` of the simulation
* `monitor`     function that is called before event is triggered
                must accept two arguments `Scheduler` and `Δ`, a difference
                between time of event to be executed and time of last executed event

If two `Action`s have identical `when` time in `event_queue` then
the order of their execution is undefined

When `monitor` is executed the event to happen is still on `event_queue`,
but time is updated to time when the event is to be executed (i.e. `monitor`
sees the state of the simulation just before the event is triggered).
Therefore for calculating summary statistics `monitor` may assume that
the simulation spent `Δ` time in this state.
Function `monitor` should not modify `event_queue[1]` as EventSimulation assumes
that the event to be triggered after `monitor` executes will not be modified.
Additionally it it not guaranteed that `event_queue[1]` will be executed after
`monitor` finishes because simulation might terminate earlier.
"""
mutable struct Scheduler{S<:AbstractState, T<:Real}
    now::T
    event_queue::Vector{Action{T}}
    state::S
    monitor::Function
end

# Construct empty `Scheduler`
# with `state` and time defined by type `T`
"""
    Scheduler(state, T)

A convenience constructor of `Scheduler`

Arguments:
* `state=EmptyState()` subtype of `AbstractState` to be used by `Scheduler`
* `T=Float64`          type that will be used to hold time in `Scheduler`
* `monitor=(s,Δ) -> nothing`  monitor function

By default an empty `event_queue` is created, `now` is set to `zero(T)`
and there is idle `monitor`
"""
function Scheduler(state=EmptyState(), T::Type{R}=Float64,
                   monitor::Function=(s,Δ) -> nothing) where R<:Real
    Scheduler(zero(T), Vector{Action{T}}(), state, monitor)
end

"""
    register!(s, what, Δ)

Put `what` at time `s.now+Δ` to `s.event_queue`.
`what` must accept exactly one argument of type `Scheduler`.
The function does not check if `Δ` is a valid (finite) number.
Returns inserted `Action`.
"""
function register!(s::Scheduler{S, T}, what::Function,
                   Δ::T = zero(T)) where {S<:AbstractState, T<:Real}
    when = s.now + Δ
    pq_insert!(s.event_queue, what, when)
end

"""
    repeat_register!(s, what, interval)

Put `what` to `s.event_queue` repeatedly in time intervals specified by
`interval` function, which must accept one argument of type `Scheduler`.
`what` must accept exactly one argument of type `Scheduler`.
`interval` function is called after the previous event was executed.
Returns `nothing`.
Calling `terminate!` in function `interval` will not stop the simulation.
Instead, if `interval` returns `nothing` the action
is not scheduled and `repeat_register` will effectively terminate.
"""
function repeat_register!(s::Scheduler, what::Function, interval::Function)
    function wrap_what(x)
        what(x)
        i = interval(x)
        isa(i, Nothing) || register!(x, wrap_what, i)
    end
    i = interval(s)
    isa(i, Nothing) || register!(s, wrap_what, i)
    nothing
end

"""
    bulk_register!(s, who, what, Δ, randomize)

Put event at time `s.now+Δ` to `s.event_queue`
that will execute `what(scheduler, w)` for all `w` in `who`.
If `randomize` is `false` then `who` is traversed in natural order
otherwise it is traversed in random order.
`what` must accept exactly two arguments of type `Scheduler` and `eltype(who)`.
The function does not check if `Δ` is a valid (finite) number.
Returns inserted bulk `Action`.

Function is designed to efficiently handle case when the same action
has to be executed at the same simulation time by many agents.
"""
function bulk_register!(s::Scheduler{S, T}, who::AbstractVector,
                        what::Function, Δ::T = zero(T),
                        randomize::Bool=false) where {S<:AbstractState, T<:Real}
    function bulk_what(x::Scheduler)
        if randomize
            for i in randperm(length(who))
                what(x, who[i])
            end
        else
            for w in who
                what(x, w)
            end
        end
    end
    register!(s, bulk_what, Δ)
end

"""
    repeat_bulk_register!(s, who, what, interval, randomize)

Repeat `bulk_register!` at time intervals specified by `interval` function,
which must accept `Scheduler` argument.
`interval` function is called after the previous event was executed.
`what` must accept exactly two arguments of type `Scheduler` and `typeof(who)`.
Returns `nothing`.
Calling `terminate!` in function `interval` will not stop the simulation.
Instead, if `interval` returns `nothing` the action
is not scheduled and `repeat_register` will effectively terminate.
"""
function repeat_bulk_register!(s::Scheduler, who::AbstractVector, what::Function,
                               interval::Function, randomize::Bool=false)
    function wrap_bulk_what(x)
        if randomize
            for j in randperm(length(who))
                what(x, who[j])
            end
        else
            for w in who
                what(x, w)
            end
        end
        i = interval(x)
        isa(i, Nothing) || register!(x, wrap_bulk_what, i)
    end
    i = interval(s)
    isa(i, Nothing) || register!(s, wrap_bulk_what, i)
    nothing
end

"""
    interrupt!(s, a)

First occurrence of `Action` `a` is replaced by no-op in event queue.
This way there is no need to fix heap in this operation and it is fast.
Returns `true` if `a` was found in queue and `false` otherwise.
"""
function interrupt!(s::Scheduler, a::Action)
    i = findfirst(isequal(a), s.event_queue)
    isa(i, Nothing) && return false
    s.event_queue[i] = Action(x -> nothing, s.event_queue[i].when)
    return true
end

"""
    terminate!(s)

Empties `s.event_queue` which will lead to termination of simulation
unless it is refilled before execution returns to `go!`.
Useful for event-triggered termination of simulation.
"""
function terminate!(s::Scheduler)
    empty!(s.event_queue)
end

"""
    go!(s, until)

Runs simulation defined by `s` until `s.now` is greater or equal than `until`
or `s.event_queue` is empty (i.e. nothing is left to be done).
By default `until` equals `Inf`.
"""
function go!(s::Scheduler, until::Real=Inf)
    while !isempty(s.event_queue)
        a = s.event_queue[1]
        when = a.when
        if when <= until
            Δ = when - s.now
            s.now = when
            s.monitor(s, Δ)
            pq_remove!(s.event_queue)
            a.what(s)
        else
            Δ = until - s.now
            s.now = until
            Δ > 0 && s.monitor(s, Δ)
            break
        end
    end
end


text:
module test_triangulate
using Triangulate

function test()
    triin=Triangulate.TriangulateIO()
    triin.pointlist=Matrix{Float64}([1.0 0.0 ; 0.0 1.0 ; -1.0 0.0 ; 0.0 -1.0]')
    triin.segmentlist=Matrix{Int32}([1 2 ; 2 3 ; 3 4 ; 4 1 ]')
    triin.segmentmarkerlist=Vector{Int32}([1, 2, 3, 4])
    triin.regionlist=Matrix{Float64}([0.5 0.5 1 0.01;]')

    (triout, vorout)=triangulate("paAqQ", triin)

    size(triout.pointlist,2)==177 && size(triout.trianglelist,2)==319 && size(triout.segmentlist,2)==33
end
end


text:
module TightBinding

    include("./TBfromk2r.jl")
    #using Plots
    #using .Plotfuncs
    using .TBfromk2r
    using LinearAlgebra
    using Requires
    export set_Lattice,add_atoms!,add_hoppings!,add_diagonals!,hamiltonian_k,
    dispersion,get_position,calc_band,get_position_kspace,hamiltonian_k_1d,
    set_Klines,show_neighbors,add_Kpoints!,set_onsite!,set_μ!,
    write_hr,make_supercell,read_wannier
    #,calc_band_plot,plotfuncs,
    #plot_lattice_2d,calc_band_plot_finite,plot_DOS

    export surfaceHamiltonian,σx,σy,σz,σ0
    export get_DOS

    struct Hopping
        amplitude
        ijatoms::Array{Int64,1}
        ijpositions::Array{Float64,1}
    end

    struct Kpoints
        kmin::Array{<:Number,1}
        kmax::Array{<:Number,1}
        nk::Int64
        name_start
        name_end
    end

    mutable struct Klines
        numlines::Int64
        kpoints::Array{Kpoints,1}
    end

    function __init__()
        @require Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" begin
            include("./Plotfuncs.jl")
            export Plotfuncs,plot_lattice_2d,calc_band_plot,plot_DOS,calc_band_plot_finite #From Plotfuncs.jl
            export plot_fermisurface_2D
        end
    end

    function set_Klines()
        kpoint = Array{Kpoints,1}[]
        klines = Klines(0,kpoint)
        return klines
    end

    function add_Kpoints!(klines,kmin,kmax,name_start,name_end;nk=20)
        kpoint = Kpoints(kmin,kmax,nk,name_start,name_end)
        klines.numlines += 1
        push!(klines.kpoints,kpoint)
    end

    mutable struct Lattice
        dim::Int8
        vectors::Array{Array{Float64,1},1}
        numatoms::Int64
        positions::Array{Array{Float64,1},1}
        numhopps::Int64
        hoppings::Array{Hopping,1}
        diagonals::Array{Float64,1}
        rvectors::Array{Array{Float64,1},1} #reciplocal lattice vectors
        μ::Float64
    end

    function Base.display(lattice::Lattice)
        println("Dimension: ", lattice.dim)
        println("Unit vectors ")
        for id=1:lattice.dim
            println(lattice.vectors[id])
        end
        println("Reciplocal lattice vectors ")
        for id=1:lattice.dim
            println(lattice.rvectors[id])
        end
        println("Chemical potential: ", lattice.μ)
        println("num. of atoms ",lattice.numatoms)
        for i=1:lattice.numatoms
            println("$(i)-th atom ",lattice.positions[i], " onsite energy: ", lattice.diagonals[i])
        end

        println("num. of hoppings ",lattice.numhopps)

        if lattice.numhopps != 0
            println("iband jband hopping amplitude")
            for i=1:lattice.numhopps
                hopping = lattice.hoppings[i]
                iband = hopping.ijatoms[1]
                jband = hopping.ijatoms[2]
                println("$iband $jband $(hopping.ijpositions) $(hopping.amplitude)")
            end
        end

        #=

        #println("iband jband hopping amplitude")
        for iatom=1:lattice.numatoms
            for jatom=1:lattice.numatoms
                for i=1:lattice.numhopps
                    hopping = lattice.hoppings[i]
                    iband = hopping.ijatoms[1]
                    jband = hopping.ijatoms[2]
                    if iatom == iband && jatom == jband

                        if jatom < iatom
                            println("$jband $iband $(-hopping.ijpositions) $(hopping.amplitude)")
                        else
                            println("$iband $jband $(hopping.ijpositions) $(hopping.amplitude)")
                        end
                    end
                end
            end
        end
        #println("num. of hoppings ",lattice.numhopps)
        =#

    end

    """
        set_Lattice(dim::Integer,vectors::Array{Array{Float64,1},1})
    Initialize lattice.
    We have to call this before making Hamiltonian.
    dim: Dimension of the system.
    vector:: Primitive vectors.

    Example:

    1D system

    ```julia
    la1 = set_Lattice(1,[[1.0]])
    ```

    2D system

    ```julia
    a1 = [sqrt(3)/2,1/2]
    a2 = [0,1]
    la2 = set_Lattice(2,[a1,a2])
    ```

    3D system

    ```julia
    a1 = [1,0,0]
    a2 = [0,1,0]
    a3 = [0,0,1]
    la2 = set_Lattice(3,[a1,a2,a3])
    ```
    """
    function set_Lattice(dim,vectors)
        positions = Array{Float64,1}[]
        hoppings = Array{Hopping,1}[]
        numatoms = 0
        numhopps=0
        diagonals = Float64[]
        μ = 0.0

        #making primitive vectors
        pvector_1 = zeros(Float64,3)
        pvector_2 = zeros(Float64,3)
        pvector_3 = zeros(Float64,3)
        if dim ==1
            pvector_1[1] = vectors[1][1]
            pvector_2[2] = 1.0 # 0 1 0
            pvector_3[3] = 1.0 # 0 0 1
        elseif dim == 2
            pvector_1[1:2] = vectors[1][1:2]
            pvector_2[1:2] = vectors[2][1:2]
            pvector_3[3] = 1.0 # 0 0 1
        elseif dim == 3
            pvector_1[1:3] = vectors[1][1:3]
            pvector_2[1:3] = vectors[2][1:3]
            pvector_3[1:3] = vectors[3][1:3]
        end
        #making reciplocal lattice vectors
        area = dot(pvector_1,cross(pvector_2,pvector_3))
        rvector_1 = 2π*cross(pvector_2,pvector_3)/area
        rvector_2 = 2π*cross(pvector_3,pvector_1)/area
        rvector_3 = 2π*cross(pvector_1,pvector_2)/area

        if dim ==1
            rvectors = [[rvector_1[1]]]
        elseif dim == 2
            rvectors = [rvector_1[1:2],rvector_2[1:2]]
        elseif dim == 3
            rvectors = [rvector_1[1:3],rvector_2[1:3],rvector_3[1:3]]
        end

        lattice = Lattice(dim,vectors,numatoms,positions,numhopps,hoppings,diagonals,rvectors,μ)
        return lattice
    end

    function add_atoms!(lattice,position::Array{<:Number,1})
        lattice.numatoms += 1
        push!(lattice.positions,position)
        add_diagonals!(lattice,0.0)
    end

    function add_atoms!(lattice,positions::Array{<:Array})
        for i=1:length(positions)
            lattice.numatoms += 1
            push!(lattice.positions,positions[i])
            add_diagonals!(lattice,0.0)
        end
    end

    function set_μ!(lattice,μ)
        lattice.μ = μ
    end

    function set_onsite!(lattice,onsite::Array{<:Number,1})
        lattice.diagonals = onsite[:]
    end

    function add_hoppings!(lattice,amplitude,iatom,jatom,hopping)
        hop = Hopping(amplitude,[iatom,jatom],hopping)
        lattice.numhopps += 1
        push!(lattice.hoppings,hop)
    end

    function add_diagonals!(lattice,diagonal::Number)
        push!(lattice.diagonals,diagonal)
    end

    function add_diagonals!(lattice,diagonals::Array{<:Number,1})
        for i=1:length(diagonals)
            push!(lattice.diagonals,diagonals[i])
        end
    end

    function get_position(lattice,vec)
        position = zeros(Float64,lattice.dim)
        for i=1:lattice.dim
            position += vec[i]*lattice.vectors[i]
        end
        return position
    end

    function get_position_kspace(lattice,vec)
        position = zeros(Float64,lattice.dim)
        for i=1:lattice.dim
            position += vec[i]*lattice.rvectors[i]
        end
        return position
    end

    function get_hopindex(lattice)
        hopindices = []
        for ihop=1:lattice.numhopps
            hopping = lattice.hoppings[ihop]
            i = hopping.ijatoms[1]
            ri = lattice.positions[i]
            rhop = ri+hopping.ijpositions

            positionindex = get_positionindex(lattice,rhop)
            iabhop = positionindex[:]

#            iabhop = round.(Int,rhop).-1
            #j = hopping.ijatoms[2]
            push!(hopindices,iabhop)
        end
        return hopindices
    end

    function hamiltonian_k_1d(lattice,direction;nsites = 20,periodic=true)
        numatoms = lattice.numatoms
        diagonals = zeros(Float64,numatoms)
        vectors = lattice.vectors
        hopindices = get_hopindex(lattice)
        numhopps = lattice.numhopps
        μ = lattice.μ

        for i=1:length(lattice.diagonals)
            diagonals[i] = lattice.diagonals[i]-μ
        end

        function calc_ham(k)
            N = numatoms*nsites
            realpart_ham = zeros(Float64,N,N)
            imagpart_ham = zeros(Float64,N,N)
            for isite=1:nsites
                for i=1:numatoms
                    ii = (isite-1)*numatoms+i
                    realpart_ham[ii,ii] = diagonals[i]
                end
                for ihop=1:numhopps
                    δ = hopindices[ihop][direction]
                    hopping = lattice.hoppings[ihop]
                    i = hopping.ijatoms[1]
                    ii = (isite-1)*numatoms+i
                    jsite = isite + δ
                    if periodic
                        while jsite > nsites || jsite < 1
                            jsite += ifelse(jsite < 1,nsites,0)
                            jsite += ifelse(jsite > nsites,-nsites,0)
                        end
                    end
                    if 1 <= jsite <= nsites
                        j= hopping.ijatoms[2]
                        jj = (jsite-1)*numatoms+j

                        ak = 0.0
                        ik = 0
                        hop = hopping.ijpositions
                        for idim = 1:lattice.dim
                            if idim != direction
                                ik += 1
                                ak += k[ik]*hop[idim]
                            end
                        end

                        ampcos = hopping.amplitude*cos(ak)
                        realpart_ham[ii,jj] +=  ampcos
                        realpart_ham[jj,ii] +=  ampcos
                        ampsin = hopping.amplitude*sin(ak)
                        imagpart_ham[ii,jj] +=  ampsin
                        imagpart_ham[jj,ii] +=  -ampsin

                    end
                end
            end
            if sum(abs.(imagpart_ham)) == 0.0
                ham = realpart_ham[:,:]
            else
                ham = realpart_ham[:,:]+im*imagpart_ham[:,:]
            end
            ham
        end
        k -> calc_ham(k)
    end

    function find_orbital(lattice::Lattice,supercell,position,iorbital)
        position_sp  =position[:]
        for id=1:lattice.dim
            while position_sp[id] >= 1
                position_sp[id] -= 1
            end
            while position_sp[id] < 0
                position_sp[id] += 1
            end
        end

        originalnumatoms = div(lattice.numatoms,prod(supercell))

        position_ori  =position_sp[:]
        positionindex = zeros(Int64,lattice.dim)

        for id=1:lattice.dim
            position_ori[id] *= supercell[id]
            while position_ori[id] >= 1
                position_ori[id] -= 1
                positionindex[id] += 1
            end
            while position_ori[id] < 0
                position_ori[id] += 1
                positionindex[id] += 1
            end
        end
        #println("positionindex $positionindex")

        if lattice.dim == 1
            iatom = ( positionindex[1])*originalnumatoms + iorbital
        elseif lattice.dim == 2
            iatom = ((positionindex[2])*supercell[1] + positionindex[1])*originalnumatoms + iorbital
        elseif lattice.dim == 3
            iatom = ((positionindex[3]*supercell[2]+positionindex[2])*supercell[1] + positionindex[1])*originalnumatoms + iorbital
        end

        #println("iorbital $iorbital iatom $iatom")
        return iatom

        for iatom = 1:lattice.numatoms
            if norm(lattice.positions[iatom] - position_sp) < 1e-15
                push!(iatomcandidate,iatom)
#                return iatom
            end
        end
        if length(iatomcandidate) != 0
            return iatomcandidate[iorbital]
        end

        println(position_sp)
        error("Orbital index was not found")
    end

    function make_supercell(lattice::Lattice,supercell)
        @assert lattice.dim == length(supercell)
        numatoms = lattice.numatoms

        vectors = deepcopy(lattice.vectors)
        for id = 1:lattice.dim
            vectors[id][:] = supercell[id]*vectors[id][:]
        end

        lattice_super = set_Lattice(lattice.dim,vectors)
        diagonals = zeros(Float64,numatoms*prod(supercell))

        set_μ!(lattice_super,lattice.μ)

        if lattice.dim == 1
            icount = 0
            for i1 = 1:supercell[1]
                for iatom =1:numatoms
                    position = lattice.positions[iatom][1] + (i1-1) #*lattice.vectors[1][1]
                    position /= supercell[1]
                    add_atoms!(lattice_super,[position])
#                    icount += 1
                    icount = (i1-1)*numatoms + iatom
                    diagonals[icount] = lattice.diagonals[iatom]
                end
            end

            set_onsite!(lattice_super,diagonals)

            #=
            for iatom=1:lattice_super.numatoms
                jatom =find_orbital(lattice_super,lattice_super.positions[iatom])
                println("$jatom $iatom")
            end
            =#

            for δ=1:lattice.numhopps
                hopping = lattice.hoppings[δ] #Type:
                iband=hopping.ijatoms[1]
                jband=hopping.ijatoms[2]
                hopposition = hopping.ijpositions
                hopposition_sp = hopposition[:]/supercell[1]
                v = hopping.amplitude
                #println("$iband $jband $hopposition $hopposition_sp")
                for i1 = 1:supercell[1]
                    iiband = (i1-1)*numatoms + iband
                    atomposition_i = lattice_super.positions[iiband][:]
                    iatom = find_orbital(lattice_super,supercell,atomposition_i,iband)
#                    iatom = find_orbital(lattice_super,atomposition_i)
                    atomposition_j = atomposition_i[:] + hopposition_sp[:]
#                    jatom =find_orbital(lattice_super,atomposition_j)
                    jatom =find_orbital(lattice_super,supercell,atomposition_j,jband)
                    #println("$iatom $jatom $atomposition_i $atomposition_j")
                    add_hoppings!(lattice_super,v,iatom,jatom,hopposition_sp)
                end
            end

        elseif lattice.dim == 2

            for i2 = 1:supercell[2]
                for i1 = 1:supercell[1]
                    for iatom =1:numatoms
                        position_x = (lattice.positions[iatom][1] + (i1-1))/supercell[1]
                        position_y = (lattice.positions[iatom][2] + (i2-1))/supercell[2]

                        add_atoms!(lattice_super,[position_x,position_y])
    #                    icount += 1
                        icount = ((i2-1)*supercell[1] + (i1-1))*numatoms + iatom
                        diagonals[icount] = lattice.diagonals[iatom]
                    end
                end
            end
            set_onsite!(lattice_super,diagonals)

            for δ=1:lattice.numhopps
                hopping = lattice.hoppings[δ] #Type:
                iband=hopping.ijatoms[1]
                jband=hopping.ijatoms[2]
                hopposition = deepcopy(hopping.ijpositions)
                #println(hopposition)
                hopposition_sp = hopposition[:]
                hopposition_sp[1] /= supercell[1]
                hopposition_sp[2] /= supercell[2]

                v = hopping.amplitude
                #println("$iband $jband $hopposition $hopposition_sp")
                for i2 = 1:supercell[2]
                    for i1 = 1:supercell[1]
                        iiband = ((i2-1)*supercell[1]+(i1-1))*numatoms + iband
                        atomposition_i = lattice_super.positions[iiband][:]
                        iatom = find_orbital(lattice_super,supercell,atomposition_i,iband)
                        #println("iband iatom $iband $iatom")
#                        iatom = find_orbital(lattice_super,atomposition_i)

                        atomposition_j = atomposition_i[:] + hopposition_sp[:]
#                        jatom =find_orbital(lattice_super,atomposition_j)
                        jatom =find_orbital(lattice_super,supercell,atomposition_j,jband)
                        #println("$iatom $jatom $atomposition_i $atomposition_j")
                        add_hoppings!(lattice_super,v,iatom,jatom,hopposition_sp)
                        positionindex = get_positionindex(lattice_super,atomposition_j)
                        #println("positionindex $positionindex")
                    end
                end
            end

        elseif lattice.dim == 3
            for i3 = 1:supercell[3]
                for i2 = 1:supercell[2]
                    for i1 = 1:supercell[1]
                        for iatom =1:numatoms
                            position_x = (lattice.positions[iatom][1] + (i1-1))/supercell[1]
                            position_y = (lattice.positions[iatom][2] + (i2-1))/supercell[2]
                            position_z = (lattice.positions[iatom][3] + (i3-1))/supercell[3]

                            add_atoms!(lattice_super,[position_x,position_y,position_z])
        #                    icount += 1
                            icount = (((i3-1)*supercell[2]+i2-1)*supercell[1] + (i1-1))*numatoms + iatom
                            diagonals[icount] = lattice.diagonals[iatom]
                        end
                    end
                end
            end
            set_onsite!(lattice_super,diagonals)

            for δ=1:lattice.numhopps
                hopping = lattice.hoppings[δ] #Type:
                iband=hopping.ijatoms[1]
                jband=hopping.ijatoms[2]
                hopposition = deepcopy(hopping.ijpositions)
                #println(hopposition)
                hopposition_sp = hopposition[:]
                hopposition_sp[1] /= supercell[1]
                hopposition_sp[2] /= supercell[2]
                hopposition_sp[3] /= supercell[3]

                v = hopping.amplitude
                #println("$iband $jband $hopposition $hopposition_sp")
                for i3 = 1:supercell[3]
                    for i2 = 1:supercell[2]
                        for i1 = 1:supercell[1]
                            iiband = (((i3-1)*supercell[2]+i2-1)*supercell[1]+(i1-1))*numatoms + iband
                            atomposition_i = lattice_super.positions[iiband][:]
                            iatom = find_orbital(lattice_super,supercell,atomposition_i,iband)
#                            iatom = find_orbital(lattice_super,atomposition_i)
                            atomposition_j = atomposition_i[:] + hopposition_sp[:]
                            jatom =find_orbital(lattice_super,supercell,atomposition_j,jband)
#                            jatom =find_orbital(lattice_super,atomposition_j)
                            #println("$iatom $jatom $atomposition_i $atomposition_j")
                            add_hoppings!(lattice_super,v,iatom,jatom,hopposition_sp)
                        end
                    end
                end
            end
        end

        return lattice_super
    end

    function get_positionindex(lattice::Lattice,position)

        positionindex = zeros(Int64,lattice.dim)
        for id=1:lattice.dim
            positionindex[id] = round(Int,position[id],RoundDown)
            #=
            if position[id] < 0
                positionindex[id] = -round(Int,-position[id],RoundDown)-1
            else
                positionindex[id] = round(Int,position[id],RoundDown)
            end
            =#
        end

        return positionindex

        position_sp  =position[:]
        positionindex = zeros(Int64,lattice.dim)
        for id=1:lattice.dim
            while position_sp[id] >= 1
                positionindex[id] += 1
                position_sp[id] -= 1
            end
            while position_sp[id] < 0
                positionindex[id] -= 1
                position_sp[id] += 1
            end
        end
        return positionindex

    end

    function cut_hg(lattice::Lattice)
        numatoms = lattice.numatoms
        la = set_Lattice(lattice.dim,lattice.vectors)
        hopindices = get_hopindex(lattice)

        for iatom = 1:lattice.numatoms
            add_atoms!(la,lattice.positions[iatom][:])
        end
        set_onsite!(la,lattice.diagonals)
        set_μ!(la,lattice.μ)

        if lattice.numhopps != 0
            for δ=1:lattice.numhopps
                hopping = lattice.hoppings[δ] #Type:
                iband=hopping.ijatoms[1]
                jband=hopping.ijatoms[2]
                hoppositions = hopping.ijpositions
                v = hopping.amplitude

                double = false
                if δ > 1
                    for δ2=1:δ-1
                        hopping2 = lattice.hoppings[δ2] #Type:
                        iband2=hopping2.ijatoms[1]
                        jband2=hopping2.ijatoms[2]
                        hoppositions2 = -hopping2.ijpositions
                        if iband == jband2 && jband == iband2 &&  norm(hoppositions2-hoppositions) == 0 && δ2 != δ
                            double = true
                            break
                        end
                    end
                end
                if double == false
                    add_hoppings!(la,v,iband,jband,hoppositions)
                end
            end
        end

        return la

    end

    function read_wannier(lattice::Lattice,filename)
        la_hr = set_Lattice(lattice.dim,lattice.vectors)

        for iatom = 1:lattice.numatoms
            add_atoms!(la_hr,lattice.positions[iatom][:])
        end

        set_μ!(la_hr,lattice.μ)
        diagonals = zeros(Float64,lattice.numatoms)

        data = readlines(pwd()*"/"*filename)
        i = 2
        u = data[i]
        #println(u)
        num_wann = parse(Int64,u)
        @assert num_wann == lattice.numatoms
        i = 3
        nrpts = parse(Int64,data[i])
        factors = zeros(Int64,nrpts)

        num = div(nrpts,15)
        nmod = nrpts % 15
        count = 0
        for j=1:num
            i += 1
            u = split(data[i])
            for k=1:15
                count += 1
                factors[count] = parse(Int64,u[k])
            end
        end
        if nmod != 0
            i += 1
            u = split(data[i])
            for k=1:nmod
                count += 1
                factors[count] = parse(Int64,u[k])
            end
        end

        for k=1:nrpts
            for μ=1:num_wann
                for ν=1:num_wann
                    i += 1
                    u = split(data[i])
                    hops = parse.(Int64,u[1:3])
                    iband = parse(Int64,u[4])
                    jband = parse(Int64,u[5])
                    @assert ν == iband
                    @assert μ == jband
                    if sum(abs.(hops)) == 0 && iband == jband
                        diagonals[iband] = parse(Float64,u[6])
                    else
                        tij = parse(Float64,u[6])-im*parse(Float64,u[7])

                        if abs(tij) != 0
                            #atomposition_i = lattice.positions[iband][:]
                            #atomposition_j = lattice.positions[jband][:]

                            Ri = lattice.positions[iband][:]
                            Rj = lattice.positions[jband][:]
                            for id=1:lattice.dim
                                Rj[id] += hops[id]
                            end
                            dr = Rj-Ri

                            #for id=1:lattice.dim
                            #    atomposition_j[id] += hops[id]#*lattice.vectors[id][:]
                            #end
                            #hopping = atomposition_j- atomposition_i

                            #add_hoppings!(la_hr,tij,jband,iband,-hopping)
                            add_hoppings!(la_hr,tij,iband,jband,dr)
                        end
                    end

                end
            end
        end

        set_onsite!(la_hr,diagonals .+ lattice.μ)

        return cut_hg(la_hr)

    end

    function write_hr(lattice::Lattice;filename="wannier90_hr.dat",maxhop = 10)
        numatoms = lattice.numatoms
        hopindices = get_hopindex(lattice)

        diagonals = zeros(Float64,numatoms)
        vectors = lattice.vectors
        μ = lattice.μ
        for i=1:length(lattice.diagonals)
            diagonals[i] = lattice.diagonals[i]-μ
        end
        hopdata = Array{Int64,1}[]
        valdata = Array{Float64,1}[]
        for i=1:numatoms
            v = diagonals[i]
            iband = i
            jband = i
            hop = [0,0,0,iband,jband]
            val = [real(v),-imag(v)]
            push!(hopdata,hop)
            push!(valdata,val)
        end
        if lattice.numhopps != 0
            for δ=1:lattice.numhopps
                hopping = lattice.hoppings[δ] #Type:
                #positionindex = hopindices[δ][:]
                iband=hopping.ijatoms[1]
                jband=hopping.ijatoms[2]

                ri = lattice.positions[iband][:]
                rhop = ri+hopping.ijpositions
                positionindex = get_positionindex(lattice,rhop)

                #=
                hopping = lattice.hoppings[δ] #Type:
                iband=hopping.ijatoms[1]
                jband=hopping.ijatoms[2]
                hopposition = hopping.ijpositions[:]

                atomposition_i = lattice.positions[iband][:]
                atomposition_j = atomposition_i[:] + hopposition[:]

                positionindex = get_positionindex(lattice,atomposition_j)
                =#
                #=
                Ri = lattice.positions[iband][:]
                Rj = lattice.positions[jband][:]
                for id=1:lattice.dim
                    Rj[id] += positionindex[id]
                end
                dr = Rj-Ri
                if norm(dr - hopping.ijpositions) != 0
                    println("positionindex: $positionindex")
                    println("Ri $Ri Rj $Rj")
                    println("dr $dr ijpositions",hopping.ijpositions)
                end
                =#

                v = hopping.amplitude
                if lattice.dim == 1
                    hop = [positionindex[1],0,0,iband,jband]
                elseif lattice.dim == 2
                    hop = [positionindex[1],positionindex[2],0,iband,jband]
                elseif lattice.dim == 3
                    hop = [positionindex[1],positionindex[2],positionindex[3],iband,jband]
                end
                val = [real(v),imag(v)]
                push!(hopdata,hop)
                push!(valdata,val)

                ri = lattice.positions[jband][:]
                rhop = ri-hopping.ijpositions

                positionindex = get_positionindex(lattice,rhop)

                #=

                Ri = lattice.positions[jband][:]
                Rj = lattice.positions[iband][:]
                for id=1:lattice.dim
                    Rj[id] += positionindex[id]
                end
                dr = Rj-Ri
                if norm(dr + hopping.ijpositions) != 0
                    println("rhop $rhop")
                    println("positionindex: $positionindex")
                    println("Ri $Ri Rj $Rj")
                    println("oriRj ", lattice.positions[iband][:])
                    println("dr $dr ijpositions",-hopping.ijpositions)
                end
                =#

                if lattice.dim == 1
                    hop = [positionindex[1],0,0,jband,iband]
                elseif lattice.dim == 2
                    hop = [positionindex[1],positionindex[2],0,jband,iband]
                elseif lattice.dim == 3
                    hop = [positionindex[1],positionindex[2],positionindex[3],jband,iband]
                end
                val = [real(v),-imag(v)]
                push!(hopdata,hop)
                push!(valdata,val)

            end
        end
        count = length(hopdata)
        if lattice.dim == 1
            hopmatrix = zeros(ComplexF64,2*maxhop+1,1,1,numatoms,numatoms)
        elseif lattice.dim == 2
            hopmatrix = zeros(ComplexF64,2*maxhop+1,2*maxhop+1,1,numatoms,numatoms)
        elseif lattice.dim == 3
            hopmatrix = zeros(ComplexF64,2*maxhop+1,2*maxhop+1,2*maxhop+1,numatoms,numatoms)
        end

        for i=1:count
            hop = hopdata[i]
            val = valdata[i]
            ix = hop[1] + maxhop+1
            iy = hop[2]+ maxhop+1
            iz = hop[3]+ maxhop+1
            μ = hop[4]
            ν = hop[5]
            if lattice.dim == 1
                hopmatrix[ix,1,1,μ,ν] = val[1] + im*val[2]
            elseif lattice.dim == 2
                hopmatrix[ix,iy,1,μ,ν] = val[1] + im*val[2]
            elseif lattice.dim == 3
                hopmatrix[ix,iy,iz,μ,ν] = val[1] + im*val[2]
            end
        end

        hopcount = 0
        if lattice.dim == 1
            for ix=1:2*maxhop+1
                if sum(abs.(hopmatrix[ix,1,1,:,:])) != 0
                    hopcount += 1
                end
            end
        elseif lattice.dim == 2
            for ix=1:2*maxhop+1
                for iy=1:2*maxhop+1
                    if sum(abs.(hopmatrix[ix,iy,1,:,:])) != 0
                        hopcount += 1
                    end
                end
            end
        elseif lattice.dim == 3
            for ix=1:2*maxhop+1
                for iy=1:2*maxhop+1
                    for iz = 1:2*maxhop+1
                        if sum(abs.(hopmatrix[ix,iy,iz,:,:])) != 0
                            hopcount += 1
                        end
                    end
                end
            end
        end

        count = hopcount #*numatoms*numatoms

        fp = open(filename,"w")
        println(fp,"generated by TightBinding.jl")
        println(fp,numatoms)
        println(fp,count)
        num = div(count,15)
        nmod = count % 15
        for j=1:num
            for i=1:15
                print(fp,1,"\t")
            end
            println(fp,"\t")
        end
        if nmod != 0
            for i=1:nmod
                print(fp,1,"\t")
            end
            println(fp,"\t")
        end

        if lattice.dim == 1
            for i1=1:2*maxhop+1
                if sum(abs.(hopmatrix[i1,1,1,:,:])) != 0
                    ix = i1 - maxhop-1
                    iy = 0
                    iz = 0
                    for ν=1:numatoms
                        for μ=1:numatoms
                            val = [real(hopmatrix[i1,1,1,μ,ν]),-imag(hopmatrix[i1,1,1,μ,ν])]
                            println(fp,"$ix\t$iy\t$iz\t$μ\t$ν\t$(val[1])\t$(val[2])")
                        end
                    end
                end
            end
        elseif lattice.dim == 2
            for i1=1:2*maxhop+1
                for i2=1:2*maxhop+1
                    if sum(abs.(hopmatrix[i1,i2,1,:,:])) != 0
                        ix = i1 - maxhop-1
                        iy = i2 - maxhop-1
                        iz = 0

                        for ν=1:numatoms
                            for μ=1:numatoms
                                val = [real(hopmatrix[i1,i2,1,μ,ν]),imag(hopmatrix[i1,i2,1,μ,ν])]
                                println(fp,"$ix\t$iy\t$iz\t$μ\t$ν\t$(val[1])\t$(val[2])")
                            end
                        end
                    end
                end
            end
        elseif lattice.dim == 3
            for i1=1:2*maxhop+1
                for i2=1:2*maxhop+1
                    for i3 = 1:2*maxhop+1
                        if sum(abs.(hopmatrix[i1,i2,i3,:,:])) != 0
                            ix = i1 - maxhop-1
                            iy = i2 - maxhop-1
                            iz = i3 - maxhop-1

                            for ν=1:numatoms
                                for μ=1:numatoms
                                    val = [real(hopmatrix[i1,i2,i3,μ,ν]),imag(hopmatrix[i1,i2,i3,μ,ν])]
                                    println(fp,"$ix\t$iy\t$iz\t$μ\t$ν\t$(val[1])\t$(val[2])")
                                end
                            end
                        end
                    end
                end
            end
        end
        close(fp)
        return

        for i=1:count
            hop = hopdata[i]
            val = valdata[i]
            ix = hop[1]
            iy = hop[2]
            iz = hop[3]
            μ = hop[4]
            ν = hop[5]

            println(fp,"$ix\t$iy\t$iz\t$μ\t$ν\t$(val[1])\t$(val[2])")
        end
        close(fp)

        #println(hopdata)
        #println(valdata)
    end

    function hamiltonian_k(lattice)
        numatoms = lattice.numatoms
        diagonals = zeros(Float64,numatoms)
        vectors = lattice.vectors
        μ = lattice.μ
        for i=1:length(lattice.diagonals)
            diagonals[i] = lattice.diagonals[i]-μ
        end

        function calc_ham(k)
            realpart_ham = zeros(Float64,numatoms,numatoms)
            imagpart_ham = zeros(Float64,numatoms,numatoms)
            for i=1:numatoms
                realpart_ham[i,i] = diagonals[i]
            end
            if lattice.numhopps != 0
                for δ=1:lattice.numhopps
                    hopping = lattice.hoppings[δ] #Type:
                    i=hopping.ijatoms[1]
                    j=hopping.ijatoms[2]
                    hop = hopping.ijpositions
                    vec_a = zeros(Float64,lattice.dim)
                    for i=1:length(lattice.vectors)
                        vec_a += hop[i].*lattice.vectors[i]
                    end
                    ak = 0.0
                    for idim = 1:lattice.dim
                        ak += k[idim]*vec_a[idim]
                    end

                    ampcos = hopping.amplitude*cos(ak)
                    realpart_ham[i,j] +=  ampcos
                    realpart_ham[j,i] +=  ampcos
                    ampsin = hopping.amplitude*sin(ak)
                    imagpart_ham[i,j] +=  ampsin
                    imagpart_ham[j,i] +=  -ampsin

                end
            end
            if sum(abs.(imagpart_ham)) == 0.0
                ham = realpart_ham[:,:]
            else
                ham = realpart_ham[:,:]+im*imagpart_ham[:,:]
            end

            ham

        end
        k -> calc_ham(k)
    end

    function dispersion(dim,n,ham,k)
        if n==1
            return [ham(k)[1]]
        else
            energy = eigen(ham(k)).values
            return energy
        end
    end

    function calc_band(kmin_real,kmax_real,nk,lattice,ham)
        n = lattice.numatoms
        dim = lattice.dim
        energies = zeros(Float64,n,nk)
        vec_k = zeros(Float64,dim,nk)
        for idim=1:dim
            k = range(kmin_real[idim],length=nk,stop=kmax_real[idim])
            vec_k[idim,:] = k[:]
        end
        for ik = 1:nk
            energies[:,ik] = dispersion(dim,n,ham,vec_k[:,ik])
        end

        return vec_k,energies
    end

    function show_neighbors(lattice;maxn = 1,cutoff=2.0)
        dim = lattice.dim
        println("Possible hoppings")
        if dim==1
            for i1 = -maxn:maxn
                for i=1:lattice.numatoms
                    r0 = [i1]
                    for i=1:lattice.numatoms
                        ri = lattice.positions[i]
                        for j=i:lattice.numatoms
                            rj = lattice.positions[j]+r0
                            r = rationalize.(rj-ri)
                            r_pos = get_position(lattice,r)
                            if norm(r_pos) < cutoff
                                println("($i,$j), x:$(r[1])")
                            end
                        end
                    end
                end
            end

        elseif dim == 2
            for i1 = -maxn:maxn
                for i2 = -maxn:maxn
                    r0 = [i1,i2]
                    for i=1:lattice.numatoms
                        ri = lattice.positions[i]
                        for j=i:lattice.numatoms
                            rj = lattice.positions[j]+r0
                            r = rationalize.(rj-ri)
                            r_pos = get_position(lattice,r)
                            if norm(r_pos) < cutoff
                                println("($i,$j), x:$(r[1]), y:$(r[2])")
                            end
                        end
                    end
                end
            end
        elseif dim == 3
            for i1 = -maxn:maxn
                for i2 = -maxn:maxn
                    for i3 = -maxn:maxn
                        r0 = [i1,i2,i3]
                        for i=1:lattice.numatoms
                            ri = lattice.positions[i]
                            for j=i:lattice.numatoms
                                rj = lattice.positions[j]+r0
                                r = rationalize.(rj-ri)
                                r_pos = get_position(lattice,r)
                                if norm(r_pos) < cutoff
                                    println("($i,$j), x:$(r[1]), y:$(r[2]),z:$(r[3])")
                                end
                            end
                        end

                    end
                end
            end
        end
        return

    end

    function get_DOS(lattice,nk;nbins=100)
        dim = lattice.dim
        ham = hamiltonian_k(lattice)
        n = lattice.numatoms

        if dim == 1
            energies = zeros(Float64,n,nk)
            dk = 1/(nk-1)
            for ik = 1:nk
                k = (ik-1)*dk
            end
        elseif dim == 2
            energies = zeros(Float64,n,nk,nk)
            dk1 = 1/(nk-1)
            dk2 = 1/(nk-1)
            for ik1 = 1:nk
                k1 = (ik1-1)*dk1
                for ik2 = 1:nk
                    k2 = (ik2-1)*dk2
                    kx,ky = get_position_kspace(lattice,[k1,k2])
                    energy = dispersion(dim,n,ham,[kx,ky])
                    energies[:,ik1,ik2] = energy[:]

                end
            end
        elseif dim == 3
            energies = zeros(Float64,nk,nk,nk)
            dk1 = 1/(nk-1)
            dk2 = 1/(nk-1)
            dk3 = 1/(nk-1)

            for ik1 = 1:nk
                k1 = (ik1-1)*dk1
                for ik2 = 1:nk
                    k2 = (ik2-1)*dk2
                    for ik3=1:nk
                        k3= (ik3-1)*dk3
                        kx,ky,kz = get_position_kspace(lattice,[k1,k2,k3])
                        energy = dispersion(dim,n,ham,[kx,ky,kz])
                        energies[:,ik1,ik2,ik3] = energy[:]
                    end

                end
            end
        end
        energies = vec(energies)
#            println(energies)
        #pls = histogram(energies,bins=nbins,normalize=true)
        hist = fit(Histogram,energies,nbins=nbins)

        return hist
    end

    include("./Testfuncs.jl")

end # module


text:
using Test

@testset "all" begin

end

max_stars_repo_name: Quantum-Many-Body/TensorNetworks

max_issues_repo_name: Quantum-Many-Body/TensorNetworks

max_forks_repo_name: Quantum-Many-Body/TensorNetworks


text:
export InitialData

"""
Abstract type defining initial data from which initial-value problems can be built,
through `Problem( model :: AbstractModel, initial :: InitialData, param :: NamedTuple )`
"""
abstract type InitialData end

show(io::IO, i::InitialData) =
    try print(io,i.info)
    catch
        print(io,"Initial data: $(i.label)")
    end

function Base.:(==)(i1::InitialData, i2::InitialData)
    x=rand(2)
    i1.η(x) == i2.η(x) && i1.v(x) == i2.v(x) && i1.label == i2.label
end

export Init
"""
    Init(data ; fast, label)

Generate an initial data to be used in the function `Problem`.

`data` should contain either

- a function `η` and a function `v` (in this order);
- a `mesh` (of type `Mesh`) and two vectors representing `η(mesh.x)` and `v(mesh.x)` (in this order);
- an array of collocation points and two vectors representing `η(x)` and `v(x)` (in this order).

In the last two cases, an optional keyword argument `fast` can be set to `true`,
(default is `false`), in which case the algorithm is faster and uses less allocations, but is less precise.

In the last case, the collocation points must be regularly spaced, otherwise an `ErrorException` is raised.

If the keyword `label::String` (used to display information to the output stream)
is not provided, then it is set to the `"user-defined"`.
"""
struct Init <: InitialData

    η :: Any
    v :: Any
    label :: String

    function Init(η , v;  label = "user-defined")
        new( x->η(x) , x->v(x), label )
    end

    function Init( mesh :: Mesh, η , v ; fast = false,  label = "user-defined")
        new( x->interpolate(mesh,η,x;fast=fast) , x->interpolate(mesh,v,x;fast=fast), label )
    end

    function Init( x :: Array, η , v ; fast = false, label = "user-defined")
        if !(x[2:end].-x[2]≈x[1:end-1].-x[1])
            @error("Collocation points must be equally spaced.")
        end
        mesh=Mesh(x)
        Init(mesh, η , v ; fast = fast, label = label)
    end

end


text:
using ParserCombinator

struct SExpr
    vec
end

expry         = Delayed()
floaty       = p"[-+]?[0-9]*\.[0-9]+([eE][-+]?[0-9]+)?" > (x -> parse(Float32, x[1:end]))
opt_ws       = p"\s" | e""
inty         = p"[-+]?\d+" > (x -> parse(Int, x))
raty         = inty + E"/" + inty > Rational{Int}
rooty        = E"(root-obj" + ~opt_ws + expry + ~opt_ws + inty + E")" > ((x,y)->AlgebraicNumber(to_expr(x), y))
booly        = p"(true|false)" > (x -> x == "true" ? true : false)
symboly      = p"[^\d(){}#'`,@~;~\[\]\s][^\s()#'`,@~;{}~\[\]]*" > Symbol
sexpr        = E"(" + ~opt_ws + Repeat(expry + ~opt_ws) + E")" |> (x -> SExpr(x))

expry.matcher =  floaty | raty | inty | booly | rooty | symboly | sexpr

top_level    = Repeat(~opt_ws + expry) + ~opt_ws + Eos()

to_expr(x) = x isa SExpr ? Expr(:call, map(to_expr, x.vec)...) : x

function parse_sexpr(s::String)
  x = parse_one(s, expry)
  x[1]
end


text:
# SingleDecision #
# ========================== #
const FixingConstraint{T} = CI{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}
mutable struct SingleDecisionConstraintBridge{T, S <: MOI.AbstractScalarSet} <: MOIB.Constraint.AbstractBridge
    constraint::CI{MOI.SingleVariable, S}
    fixing_constraint::FixingConstraint{T}
    decision::SingleDecision
end

function MOIB.Constraint.bridge_constraint(::Type{SingleDecisionConstraintBridge{T,S}},
                                           model,
                                           f::SingleDecision,
                                           set::S) where {T, S <: MOI.AbstractScalarSet}
    # Perform the bridge mapping manually
    g = MOIB.bridged_variable_function(model, f.decision)
    mapped_variable = MOI.SingleVariable(only(g.terms).variable_index)
    # Add the bridged constraint
    constraint = MOI.add_constraint(model,
                                    mapped_variable,
                                    set)
    # Check state of decision
    fixing_constraint = if !iszero(g.constant)
        # Decision initially fixed
        fixing_constraint =
            MOI.add_constraint(model,
                               MOI.ScalarAffineFunction{T}([only(g.terms)], zero(T)),
                               MOI.EqualTo(g.constant))
    else
        fixing_constraint = FixingConstraint{T}(0)
    end
    # Save the constraint index and the decision to allow modifications
    return SingleDecisionConstraintBridge{T,S}(constraint, fixing_constraint, f)
end

function MOIB.Constraint.bridge_constraint(::Type{SingleDecisionConstraintBridge{T,FreeDecision}},
                                           model,
                                           f::SingleDecision,
                                           set::FreeDecision) where T
    # Do not need to add constraint, just save bridge handle fixing
    constraint = CI{MOI.SingleVariable, FreeDecision}(0)
    return SingleDecisionConstraintBridge{T,FreeDecision}(constraint, FixingConstraint{T}(0), f)
end

function MOI.supports_constraint(::Type{<:SingleDecisionConstraintBridge{T}},
                                 ::Type{SingleDecision},
                                 ::Type{<:MOI.AbstractScalarSet}) where T
    return true
end
function MOIB.added_constrained_variable_types(::Type{<:SingleDecisionConstraintBridge})
    return Tuple{DataType}[]
end
function MOIB.added_constraint_types(::Type{<:SingleDecisionConstraintBridge{T,S}}) where {T,S}
    return [(MOI.SingleVariable, S), (MOI.ScalarAffineFunction{T}, MOI.EqualTo{T})]
end
function MOIB.Constraint.concrete_bridge_type(::Type{<:SingleDecisionConstraintBridge{T}},
                              ::Type{SingleDecision},
                              S::Type{<:MOI.AbstractScalarSet}) where T
    return SingleDecisionConstraintBridge{T,S}
end

MOI.get(b::SingleDecisionConstraintBridge{T,S}, ::MOI.NumberOfConstraints{MOI.SingleVariable, S}) where {T,S} = 1
MOI.get(b::SingleDecisionConstraintBridge{T}, ::MOI.NumberOfConstraints{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}) where T =
    b.fixing_constraint.value == 0 ? 0 : 1
MOI.get(b::SingleDecisionConstraintBridge{T}, ::MOI.ListOfConstraintIndices{MOI.SingleVariable, S}) where {T,S} = [b.constraint]
MOI.get(b::SingleDecisionConstraintBridge{T}, ::MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}) where {T} =
    b.fixing_constraint.value == 0 ? FixingConstraint{T}[] : [b.fixing_constraint]

function MOI.get(model::MOI.ModelLike, attr::MOI.ConstraintFunction,
                 bridge::SingleDecisionConstraintBridge)
    return bridge.decision
end

function MOI.get(model::MOI.ModelLike, attr::MOI.ConstraintSet,
                 bridge::SingleDecisionConstraintBridge)
    return MOI.get(model, MOI.ConstraintSet(), bridge.constraint)
end

function MOI.get(model::MOI.ModelLike, attr::MOI.ConstraintPrimal,
                 bridge::SingleDecisionConstraintBridge)
    return MOI.get(model, attr, bridge.constraint)
end

function MOI.get(model::MOI.ModelLike, attr::MOI.ConstraintDual,
                 bridge::SingleDecisionConstraintBridge)
    return MOI.get(model, attr, bridge.constraint)
end

function MOI.delete(model::MOI.ModelLike, bridge::SingleDecisionConstraintBridge)
    if bridge.constraint.value != 0
        MOI.delete(model, bridge.constraint)
    end
    return nothing
end

function MOI.set(model::MOI.ModelLike, ::MOI.ConstraintSet,
                 bridge::SingleDecisionConstraintBridge{T,S}, change::S) where {T, S <: MOI.AbstractScalarSet}
    MOI.set(model, MOI.ConstraintSet(), bridge.constraint, change)
    return nothing
end

function MOI.modify(model::MOI.ModelLike, bridge::SingleDecisionConstraintBridge{T,S}, change::DecisionStateChange) where {T,S}
    # Perform the bridge mapping manually
    if bridge.decision.decision != change.decision
        # Decision not in constraint, nothing to do
        return nothing
    end
    # Switch on state transition
    if change.new_state == NotTaken
        if bridge.fixing_constraint.value != 0
            # Remove the fixing constraint
            MOI.delete(model, bridge.fixing_constraint)
            bridge.fixing_constraint = FixingConstraint{T}(0)
        end
    end
    if change.new_state == Taken
        if bridge.fixing_constraint.value != 0
            # Remove any existing fixing constraint
            MOI.delete(model, bridge.fixing_constraint)
        end
        # Perform the bridge mapping manually
        aff = MOIB.bridged_variable_function(model, bridge.decision.decision)
        f = MOI.ScalarAffineFunction{T}([only(aff.terms)], zero(T))
        # Get the decision value
        set = MOI.EqualTo(aff.constant)
        # Add a fixing constraint to ensure that fixed decision is feasible.
        bridge.fixing_constraint = MOI.add_constraint(model, f, set)
    end
    return nothing
end