text:
using Plots
using simulator.materials.spesifics

# plotly()
gr()

λs = 500e-9:1e-9:1600e-9

fig_n_Au = plot(title="Real part of gold refractive indicies", legend=:outerright, xlabel="Wavelength [nm]", ylabel="n", ticks=:native)
plot!(fig_n_Au, λs .* 1e9, real(Au_Johnson.(λs)), label="Johnson")
plot!(fig_n_Au, λs .* 1e9, real(Au_Werner.(λs)), label="Werner")
plot!(fig_n_Au, λs .* 1e9, real(Au_unloaded.(λs)), label="Palm 2019")
plot!(fig_n_Au, λs .* 1e9, real(Au_McPeak.(λs)), label="McPeak")
plot!(fig_n_Au, λs .* 1e9, real(Au_Babar.(λs)), label="Babar")
plot!(fig_n_Au, λs .* 1e9, real(Au_OlmonEvaporated.(λs)), label="Olmon evaporated")
plot!(fig_n_Au, λs .* 1e9, real(Au_OlmonSingleChrystaline.(λs)), label="Olmon single chrystaline")
plot!(fig_n_Au, λs .* 1e9, real(Au_OlmonTemplateStripped.(λs)), label="Olmon template stripped")
plot!(fig_n_Au, λs .* 1e9, real(Au_11nm_Rosenblatt.(λs)), label="Rosenblatt 11 nm film")
plot!(fig_n_Au, λs .* 1e9, real(Au_21nm_Rosenblatt.(λs)), label="Rosenblatt 21 nm film")
plot!(fig_n_Au, λs .* 1e9, real(Au_44nm_Rosenblatt.(λs)), label="Rosenblatt 44 nm film")
plot!(fig_n_Au, λs .* 1e9, real(SiO2_core_Sellmeier.(λs)), label="SiO2 core")
display(fig_n_Au)
# savefig(fig_n_Au, "../rapport/figures/IndexSourceComparisons/Au_n.pdf")

fig_k_Au = plot(title="Imaginary part of gold refractive indicies", legend=:bottomleft, xlabel="Wavelength [nm]", ylabel="k", ticks=:native)
plot!(fig_k_Au, λs .* 1e9, imag(Au_Johnson.(λs)), label="Johnson")
plot!(fig_k_Au, λs .* 1e9, imag(Au_Werner.(λs)), label="Werner")
plot!(fig_k_Au, λs .* 1e9, imag(Au_unloaded.(λs)), label="Palm 2019")
plot!(fig_k_Au, λs .* 1e9, imag(Au_McPeak.(λs)), label="McPeak")
plot!(fig_k_Au, λs .* 1e9, imag(Au_Babar.(λs)), label="Babar")
plot!(fig_k_Au, λs .* 1e9, imag(Au_OlmonEvaporated.(λs)), label="Olmon evaporated")
plot!(fig_k_Au, λs .* 1e9, imag(Au_OlmonSingleChrystaline.(λs)), label="Olmon single chrystaline")
plot!(fig_k_Au, λs .* 1e9, imag(Au_OlmonTemplateStripped.(λs)), label="Olmon template stripped")
plot!(fig_k_Au, λs .* 1e9, imag(Au_11nm_Rosenblatt.(λs)), label="Rosenblatt 11 nm film")
plot!(fig_k_Au, λs .* 1e9, imag(Au_21nm_Rosenblatt.(λs)), label="Rosenblatt 21 nm film")
plot!(fig_k_Au, λs .* 1e9, imag(Au_44nm_Rosenblatt.(λs)), label="Rosenblatt 44 nm film")
plot!(fig_k_Au, λs .* 1e9, imag(SiO2_core_Sellmeier.(λs)), label="SiO2 core")
display(fig_k_Au)
# savefig(fig_k_Au, "../rapport/figures/IndexSourceComparisons/Au_k.pdf")

fig_n_Pd = plot(title="Real part of palladium refractive indicies", legend=:outerright, xlabel="Wavelength [nm]", ylabel="n", ticks=:native)
plot!(fig_n_Pd, λs .* 1e9, real(Pd_Johnson.(λs)), label="Johnson")
plot!(fig_n_Pd, λs .* 1e9, real(Pd_Werner.(λs)), label="Werner")
plot!(fig_n_Pd, λs .* 1e9, real(Pd_Palm_2018.(λs)), label="Palm 2018")
plot!(fig_n_Pd, λs .* 1e9, real(Pd_unloaded.(λs)), label="Palm 2019 unloaded")
plot!(fig_n_Pd, λs .* 1e9, real(Pd_loaded.(λs)), label="Palm 2019 loaded")
plot!(fig_n_Pd, λs .* 1e9, real(SiO2_core_Sellmeier.(λs)), label="SiO2 core")
display(fig_n_Pd)
# savefig(fig_n_Pd, "../rapport/figures/IndexSourceComparisons/Pd_n.pdf")

fig_k_Pd = plot(title="Imaginary part of palladium refractive indicies", legend=:bottomleft, xlabel="Wavelength [nm]", ylabel="k", ticks=:native)
plot!(fig_k_Pd, λs .* 1e9, imag(Pd_Johnson.(λs)), label="Johnson")
plot!(fig_k_Pd, λs .* 1e9, imag(Pd_Werner.(λs)), label="Werner")
plot!(fig_k_Pd, λs .* 1e9, imag(Pd_Palm_2018.(λs)), label="Palm 2018")
plot!(fig_k_Pd, λs .* 1e9, imag(Pd_unloaded.(λs)), label="Palm 2019 unloaded")
plot!(fig_k_Pd, λs .* 1e9, imag(Pd_loaded.(λs)), label="Palm 2019 loaded")
plot!(fig_k_Pd, λs .* 1e9, imag(SiO2_core_Sellmeier.(λs)), label="SiO2 core")
display(fig_k_Pd)
# savefig(fig_k_Pd, "../rapport/figures/IndexSourceComparisons/Pd_k.pdf")

fig_n_lab = plot(title="Real part of refractive index for samples", legend=:outerright, xlabel="Wavelength [nm]", ylabel="n", ticks=:native)
plot!(fig_n_lab, λs .* 1e9, real(Au_Johnson.(λs)), label="Au Johnson")
plot!(fig_n_lab, λs .* 1e9, real(Pd_Johnson.(λs)), label="Pd Johnson")

plot!(fig_n_lab, λs .* 1e9, real(Pd000S5.(λs)), label="S5 0% Pd")
plot!(fig_n_lab, λs .* 1e9, real(Pd100S20.(λs)), label="S20 100% Pd")
plot!(fig_n_lab, λs .* 1e9, real(Pd070S9to17.(λs)), label="S9-17 70% Pd")

plot!(fig_n_lab, λs .* 1e9, real(Pd041S21.(λs)), label="S21 41% Pd")
plot!(fig_n_lab, λs .* 1e9, real(Pd051S22.(λs)), label="S22 51% Pd")
plot!(fig_n_lab, λs .* 1e9, real(Pd061S23.(λs)), label="S23 61% Pd")
plot!(fig_n_lab, λs .* 1e9, real(Pd069S24.(λs)), label="S24 69% Pd")

plot!(fig_n_lab, λs .* 1e9, real(Pd064S7and18.(λs)), label="S7, 18 64% Pd")
plot!(fig_n_lab, λs .* 1e9, real(Pd053S6.(λs)), label="S6 53% Pd")
plot!(fig_n_lab, λs .* 1e9, real(Pd074S8drude.(λs)), label="S8 drude 74% Pd")
plot!(fig_n_lab, λs .* 1e9, real(Pd061S19.(λs)), label="S19 61% Pd")
plot!(fig_n_lab, λs .* 1e9, real(Pd061S19drude.(λs)), label="S19 drude 61% Pd")
display(fig_n_lab)
# savefig(fig_n_lab, "../rapport/figures/thesisfigures/Samples_n.pdf")

fig_k_lab = plot(title="Imaginary part of refractive index for samples", legend=:bottomleft, xlabel="Wavelength [nm]", ylabel="k", ticks=:native)
plot!(fig_k_lab, λs .* 1e9, imag(Au_Johnson.(λs)), label="Au Johnson")
plot!(fig_k_lab, λs .* 1e9, imag(Pd_Johnson.(λs)), label="Pd Johnson")

plot!(fig_k_lab, λs .* 1e9, imag(Pd000S5.(λs)), label="S5 0% Pd")
plot!(fig_k_lab, λs .* 1e9, imag(Pd100S20.(λs)), label="S20 100% Pd")
plot!(fig_k_lab, λs .* 1e9, imag(Pd070S9to17.(λs)), label="S9-17 70% Pd")

plot!(fig_k_lab, λs .* 1e9, imag(Pd041S21.(λs)), label="S21 41% Pd")
plot!(fig_k_lab, λs .* 1e9, imag(Pd051S22.(λs)), label="S22 51% Pd")
plot!(fig_k_lab, λs .* 1e9, imag(Pd061S23.(λs)), label="S23 61% Pd")
plot!(fig_k_lab, λs .* 1e9, imag(Pd069S24.(λs)), label="S24 69% Pd")

plot!(fig_k_lab, λs .* 1e9, imag(Pd064S7and18.(λs)), label="S7, 18 64% Pd")
plot!(fig_k_lab, λs .* 1e9, imag(Pd053S6.(λs)), label="S6 53% Pd")
plot!(fig_k_lab, λs .* 1e9, imag(Pd074S8drude.(λs)), label="S8 drude 74% Pd")
plot!(fig_k_lab, λs .* 1e9, imag(Pd061S19.(λs)), label="S19 61% Pd")
plot!(fig_k_lab, λs .* 1e9, imag(Pd061S19drude.(λs)), label="S19 drude 61% Pd")
display(fig_k_lab)
# savefig(fig_k_lab, "../rapport/figures/thesisfigures/Samples_k.pdf")

fig_n_comp = plot(title="Real part of refractive index for samples", legend=:outerright, xlabel="Wavelength [nm]", ylabel="n", ticks=:native)
plot!(fig_n_comp, λs .* 1e9, real(SiO2_core_Sellmeier.(λs)), label="SiO2 core")

plot!(fig_n_comp, λs .* 1e9, real(Pd041S21.(λs)), label="S21 41% Pd")
plot!(fig_n_comp, λs .* 1e9, real(Pd051S22.(λs)), label="S22 51% Pd")
plot!(fig_n_comp, λs .* 1e9, real(Pd053S6.(λs)), label="S6 53% Pd")
plot!(fig_n_comp, λs .* 1e9, real(Pd070S9to17.(λs)), label="S9-17 70% Pd")
plot!(fig_n_comp, λs .* 1e9, real(Pd074S8drude.(λs)), label="S8 drude 74% Pd")

plot!(fig_n_comp, λs .* 1e9, real(Pd042_unloaded.(λs)), label="Palm 42% Pd")
plot!(fig_n_comp, λs .* 1e9, real(Pd052_unloaded.(λs)), label="Palm 52% Pd")
plot!(fig_n_comp, λs .* 1e9, real(Pd073_unloaded.(λs)), label="Palm 73% Pd")
display(fig_n_comp)
# savefig(fig_n_comp, "../rapport/figures/thesisfigures/Samples_n_comp.pdf")

fig_k_comp = plot(title="Imaginary part of refractive index for samples", legend=:bottomleft, xlabel="Wavelength [nm]", ylabel="k", ticks=:native)
plot!(fig_k_comp, λs .* 1e9, imag(SiO2_core_Sellmeier.(λs)), label="SiO2 core")

plot!(fig_k_comp, λs .* 1e9, imag(Pd041S21.(λs)), label="S21 41% Pd")
plot!(fig_k_comp, λs .* 1e9, imag(Pd051S22.(λs)), label="S22 51% Pd")
plot!(fig_k_comp, λs .* 1e9, imag(Pd053S6.(λs)), label="S6 53% Pd")
plot!(fig_k_comp, λs .* 1e9, imag(Pd070S9to17.(λs)), label="S9-17 70% Pd")
plot!(fig_k_comp, λs .* 1e9, imag(Pd074S8drude.(λs)), label="S8 drude 74% Pd")

plot!(fig_k_comp, λs .* 1e9, imag(Pd042_unloaded.(λs)), label="Palm 42% Pd")
plot!(fig_k_comp, λs .* 1e9, imag(Pd052_unloaded.(λs)), label="Palm 52% Pd")
plot!(fig_k_comp, λs .* 1e9, imag(Pd073_unloaded.(λs)), label="Palm 73% Pd")
display(fig_k_comp)
# savefig(fig_k_comp, "../rapport/figures/thesisfigures/Samples_k_comp.pdf")


text:
using H3
using H3.Lib
using H3.Lib: H3Index, GeoCoord, Geofence, GeoPolygon
import H3.API.geoToH3
using Geodesy
using DataStructures

using Plots
using HexGeometry

# The goal here is to search a plausible sized grid for the nearest n point.
# The MRMS data has a grid 3500 x 7000 points which are roughly 1km apart
# Here, we will use a test grid that is evenly spaced in latitude from 20-50 degrees
# and longitude from -55 to -120 degrees

struct Location
    lat::Float64
    lon::Float64
end

"""
grid(res=(7000,3500), r=1000, maxResolution=11, minResolution=6)

Creates a grid of points with a specified number of points in east-west and
north-south directions (defaul 7000ew x 3500ns). Points will have a specified
spacing (default 1000m = 1km).

In addition to locating points on a grid, an index will be created
so that points can be found using different scales of H3Index.
"""
function grid(;res=(7000,3500), r=1000.0, maxResolution=11, minResolution=6)
    center = (40.0,-95)
    # planar projection tangent to roughly the center of continental US
    tx = ENUfromLLA(LLA(center...,0), wgs84)
    ix = inv(tx)
    grid = Vector{Location}()
    index = DefaultDict{H3Index, Vector{Int32}}(Vector{Int32})
    x0, y0 = res ./ 2
    id = 0
    grid = Vector{Location}(undef, reduce(*, res .+ 1))
    for x = -x0:x0, y = -y0:y0
        id += 1
        px = ix(ENU(x*r, y*r, 0))
        grid[id] = Location(px.lat, px.lon)
        for resolution in minResolution:maxResolution
            hex = geoToH3(px.lat, px.lon, resolution)
            push!(index[hex], id)
        end
    end
    return grid, index
end

function nearest(query, grid, index, range)
    tx = ENUfromLLA(LLA(query...), wgs84)
    for res in max(range...):-1:min(range...)
        q = geoToH3(query..., res)
        r = [grid[z] for u in map(x -> i[x], H3.API.kRing(q, 1)) for z in u]
        if length(r) > 1
            p = map(x -> tx(LLA(x.lat, x.lon, 0.0)), r)
            distance = map(x->sum(x.^2), p)
            best = findfirst(distance .== reduce(min, distance))
            return [r[best], sqrt(distance[best])]
        end
    end
    return nothing
end

ext: jl

max_stars_repo_name: gunnarnewell/weather-server

max_issues_repo_name: gunnarnewell/weather-server

max_forks_repo_name: gunnarnewell/weather-server


text:
@testset "RCF" begin

  @testset "abelian extensions && related examples" begin

    Qx,x=PolynomialRing(FlintQQ, "x")
    K,a=NumberField(x-1,"a")
    O=maximal_order(K)
    l=Hecke.abelian_normal_extensions(O,[2,2], fmpz(10)^4)
    @test length(l)==47
    l1=collect(Hecke.C22_extensions(10^4))
    @test length(l1)==47

    K,a=NumberField(x^2+1, "a")
    O=maximal_order(K)
    l=Hecke.abelian_normal_extensions(O, Int[2], fmpz(10)^5, with_autos=Val{true})
    @test length(l)==41
    for x in l[1:5]
      K, autos=Hecke._from_relative_to_abs(x[1], x[2])
      @test length(autos)==2
      y=small_generating_set(closure(autos, *))
      @test length(y)==1 || length(y)==2
      if length(y)==2
        @test y[1]*y[2]==y[2]*y[1]
      else
        @test y[1]^2!=Hecke.NfToNfMor(K,K,gen(K))
      end
    end
  end

end


text:
# Represent as a permtuation matrix.
"""
    Matrix(p::Permutation)

Return the permutation matrix for `p`.
"""
function Matrix{T}(p::Permutation) where {T}
    n = length(p)
    A = Matrix{T}(I, n, n)     # A = eye(T,n)   #  int(eye(n))
    return A[:, p.data]
end

Matrix(p::Permutation) = Matrix{Int}(p)
Array(p::Permutation) = Matrix(p)
AbstractMatrix(p::Permutation) = Matrix(p)
AbstractArray(p::Permutation) = Matrix(p)

Array{T}(p::Permutation) where {T} = Matrix{T}(p)
AbstractMatrix{T}(p::Permutation) where {T} = Matrix{T}(p)
AbstractArray{T}(p::Permutation) where {T} = Matrix{T}(p)

function Permutation(M::Matrix)
    n, c = size(M)
    @assert n == c "Matrix must be square"
    v = [findfirst(Bool.(M[:, k])) for k = 1:n]
    return Permutation(v)
end


text:
# mechanism to disallow indexing

const _allowscalar = Ref(true)

allowscalar(flag = true) = (_allowscalar[] = flag)

function assertscalar(op = "operation")
  _allowscalar[] || error("$op is disallowed")
  return
end

macro allowscalar(ex)
    quote
        local prev = _allowscalar[]
        _allowscalar[] = true
        local ret = $(esc(ex))
        _allowscalar[] = prev
        ret
    end
end

macro disallowscalar(ex)
    quote
        local prev = _allowscalar[]
        _allowscalar[] = false
        local ret = $(esc(ex))
        _allowscalar[] = prev
        ret
    end
end

# basic indexing

Base.IndexStyle(::Type{<:GPUArray}) = Base.IndexLinear()

function _getindex(xs::GPUArray{T}, i::Integer) where T
    x = Array{T}(undef, 1)
    copyto!(x, 1, xs, i, 1)
    return x[1]
end

function Base.getindex(xs::GPUArray{T}, i::Integer) where T
    ndims(xs) > 0 && assertscalar("scalar getindex")
    _getindex(xs, i)
end

function _setindex!(xs::GPUArray{T}, v::T, i::Integer) where T
    x = T[v]
    copyto!(xs, i, x, 1, 1)
    return v
end

function Base.setindex!(xs::GPUArray{T}, v::T, i::Integer) where T
    assertscalar("scalar setindex!")
    _setindex!(xs, v, i)
end

Base.setindex!(xs::GPUArray, v, i::Integer) = xs[i] = convert(eltype(xs), v)

# Vector indexing

to_index(a, x) = x
to_index(::A, x::Array{ET}) where {A, ET} = copyto!(similar(A, ET, size(x)), x)
to_index(a, x::UnitRange{<: Integer}) = convert(UnitRange{Int}, x)
to_index(a, x::Base.LogicalIndex) = error("Logical indexing not implemented")

@generated function index_kernel(state, dest::AbstractArray, src::AbstractArray, idims, Is)
    N = length(Is.parameters)
    quote
        i = linear_index(state)
        i > length(dest) && return
        is = gpu_ind2sub(idims, i)
        @nexprs $N i -> @inbounds I_i = Is[i][is[i]]
        @inbounds dest[i] = @ncall $N getindex src i -> I_i
        return
    end
end

function Base._unsafe_getindex!(dest::GPUArray, src::GPUArray, Is::Union{Real, AbstractArray}...)
    if length(Is) == 1 && isa(first(Is), Array) && isempty(first(Is)) # indexing with empty array
        return dest
    end
    idims = map(length, Is)
    gpu_call(index_kernel, dest, (dest, src, idims, map(x-> to_index(dest, x), Is)))
    return dest
end

# FIXME: simple broadcast getindex like function... reuse from Base
@inline bgetindex(x::AbstractArray, i) = x[i]
@inline bgetindex(x, i) = x

@generated function setindex_kernel!(state, dest::AbstractArray, src, idims, Is, len)
    N = length(Is.parameters)
    idx = ntuple(i-> :(Is[$i][is[$i]]), N)
    quote
        i = linear_index(state)
        i > len && return
        is = gpu_ind2sub(idims, i)
        @inbounds setindex!(dest, bgetindex(src, i), $(idx...))
        return
    end
end

# FIXME: this should use adapt
gpu_convert(GPUType, x::GPUArray) = x
function gpu_convert(GPUType, x::AbstractArray)
    isbits(x) ? x : convert(GPUType, x)
end
function gpu_convert(GPUType, x)
    isbits(x) ? x : error("Only isbits types are allowed for indexing. Found: $(typeof(x))")
end

function Base._unsafe_setindex!(::IndexStyle, dest::T, src, Is::Union{Real, AbstractArray}...) where T <: GPUArray
    if length(Is) == 1 && isa(first(Is), Array) && isempty(first(Is)) # indexing with empty array
        return dest
    end
    idims = length.(Is)
    len = prod(idims)
    src_gpu = gpu_convert(T, src)
    gpu_call(setindex_kernel!, dest, (dest, src_gpu, idims, map(x-> to_index(dest, x), Is), len), len)
    return dest
end


text:
using Gadfly, RDatasets, Distributions

set_default_plot_size(6inch, 3inch)

plot(dataset("ggplot2", "diamonds"), x=:Price, Geom.density)

# manual densities
dist = MixtureModel(Normal, [(0.5, 0.2), (1, 0.1)])
xs = rand(dist, 10^5)
plot(layer(x=xs, Geom.density, color=["auto"]),
    layer(x=xs, Geom.density(bandwidth=0.0003), color=["bw=0.0003"]),
    layer(x=xs, Geom.density(bandwidth=0.25), color=["bw=0.25"]),
    Scale.color_discrete_manual("orange", "green", "purple"),
    Guide.colorkey(title="bandwidth"))


text:
struct LowDimArray{D,T,N,A<:DenseArray{T,N}} <: DenseArray{T,N}
    data::A
end
function LowDimArray{D}(data::A) where {D,T,N,A <: AbstractArray{T,N}}
    LowDimArray{D,T,N,A}(data)
end
Base.@propagate_inbounds Base.getindex(A::LowDimArray, i::Vararg{Union{Integer,CartesianIndex},K}) where {K} = getindex(A.data, i...)
@inline Base.size(A::LowDimArray) = Base.size(A.data)
@inline Base.size(A::LowDimArray, i) = Base.size(A.data, i)
@inline Base.strides(A::LowDimArray) = strides(A.data)
@inline ArrayInterface.parent_type(::Type{LowDimArray{D,T,N,A}}) where {T,D,N,A} = A
@inline ArrayInterface.strides(A::LowDimArray) = ArrayInterface.strides(A.data)
@inline ArrayInterface.device(A::LowDimArray) = ArrayInterface.CPUPointer()
@generated function ArrayInterface.size(A::LowDimArray{D,T,N}) where {D,T,N}
    t = Expr(:tuple)
    gf = GlobalRef(Core,:getfield)
    for n ∈ 1:N
        if n > length(D) || D[n]
            push!(t.args, Expr(:call, gf, :s, n, false))
        else
            push!(t.args, Expr(:call, Expr(:curly, lv(:StaticInt), 1)))
        end
    end
    Expr(:block, Expr(:meta,:inline), :(s = ArrayInterface.size(parent(A))), t)
end
Base.parent(A::LowDimArray) = getfield(A, :data)
Base.unsafe_convert(::Type{Ptr{T}}, A::LowDimArray{D,T}) where {D,T} = pointer(parent(A))
ArrayInterface.contiguous_axis(A::LowDimArray) = ArrayInterface.contiguous_axis(parent(A))
ArrayInterface.contiguous_batch_size(A::LowDimArray) = ArrayInterface.contiguous_batch_size(parent(A))
ArrayInterface.stride_rank(A::LowDimArray) = ArrayInterface.stride_rank(parent(A))
ArrayInterface.offsets(A::LowDimArray) = ArrayInterface.offsets(parent(A))

@generated function _lowdimfilter(::Val{D}, tup::Tuple{Vararg{Any,N}}) where {D,N}
    t = Expr(:tuple)
    gf = GlobalRef(Core,:getfield)
    for n ∈ 1:N
        if n > length(D) || D[n]
            push!(t.args, Expr(:call, gf, :tup, n, false))
        end
    end
    Expr(:block, Expr(:meta,:inline), t)
end

struct ForBroadcast{T,N,A<:AbstractArray{T,N}} <: AbstractArray{T,N}
    data::A
end
@inline Base.parent(fb::ForBroadcast) = getfield(fb, :data)
Base.@propagate_inbounds Base.getindex(A::ForBroadcast, i::Vararg{Any,K}) where {K} = parent(A)[i...]
const LowDimArrayForBroadcast{D,T,N,A} = ForBroadcast{T,N,LowDimArray{D,T,N,A}}
@inline function VectorizationBase.contiguous_axis(fb::LowDimArrayForBroadcast{D,T,N,A}) where {D,T,N,A}
    _contiguous_axis(Val{D}(), VectorizationBase.contiguous_axis(parent(parent(fb))))
end
@inline forbroadcast(A::AbstractArray) = ForBroadcast(A)
@inline forbroadcast(A::AbstractRange) = A
@inline forbroadcast(A) = A
@inline forbroadcast(A::Adjoint) = forbroadcast(parent(A))
@inline forbroadcast(A::Transpose) = forbroadcast(parent(A))
Base.size(fb::ForBroadcast) = size(parent(fb))

# @inline function VectorizationBase.contiguous_batch_size(fb::LowDimArrayForBroadcast{D,T,N,A}) where {D,T,N,A}
#     _contiguous_axis(Val{D}(), VectorizationBase.contiguous_batch_size(parent(parent(fb))))
# end
@generated function _contiguous_axis(::Val{D}, ::StaticInt{C}) where {D,C}
    Dlen = length(D)
    (C > 0) || return Expr(:block,Expr(:meta,:inline), staticexpr(C))
    if C ≤ Dlen
        D[C] || return Expr(:block,Expr(:meta,:inline), staticexpr(-1))
    end
    Cnew = 0
    for n ∈ 1:C
        Cnew += ((n > Dlen)) || D[n]
    end
    Expr(:block,Expr(:meta,:inline), staticexpr(Cnew))
end
@inline function VectorizationBase.val_stride_rank(fb::LowDimArrayForBroadcast{D}) where {D}
    VectorizationBase.asvalint(_lowdimfilter(Val(D), ArrayInterface.stride_rank(parent(parent(fb)))))
end
@inline function VectorizationBase.val_dense_dims(fb::LowDimArrayForBroadcast{D}) where {D}
    VectorizationBase.asvalbool(_lowdimfilter(Val(D), ArrayInterface.dense_dims(parent(parent(fb)))))
end
@inline function VectorizationBase.bytestrides(fb::LowDimArrayForBroadcast{D}) where {D}
    p = parent(parent(fb))
    s = _bytestrides(ArrayInterface.size(p), ArrayInterface.strides(p), p)
    _lowdimfilter(Val(D), s)
end
@inline function ArrayInterface.offsets(fb::LowDimArrayForBroadcast{D}) where {D}
    _lowdimfilter(Val(D), ArrayInterface.offsets(parent(parent(fb))))
end

for f ∈ [ # groupedstridedpointer support
    :(VectorizationBase.memory_reference),
    :(ArrayInterface.contiguous_axis),
    :(ArrayInterface.contiguous_batch_size),
    :(ArrayInterface.device),
    :(VectorizationBase.val_stride_rank),
    :(VectorizationBase.val_dense_dims),
    :(ArrayInterface.offsets)
]
    @eval @inline $f(fb::ForBroadcast) = $f(getfield(fb, :data))
end
@inline _bytestrides(s,paren) = VectorizationBase.bytestrides(paren)
@generated function _bytestrides(s::Tuple{Vararg{Integer,N}}, x::Tuple{Vararg{Integer,N}}, paren::AbstractArray{T,N}) where {T,N}
    q = Expr(:block, Expr(:meta,:inline))
    strd_tup = Expr(:tuple)
    gf = GlobalRef(Core, :getfield)
    ifel = GlobalRef(Core, :ifelse)
    st = staticexpr(sizeof(T))
    for n ∈ 1:N
        s_type = s.parameters[n]
        if s_type <: Static
            if s_type === One
                push!(strd_tup.args, Expr(:call, lv(:Zero)))
            else
                push!(strd_tup.args, :($gf(x, $n, false) * $st))
            end
        else
            Xₙ_type = x.parameters[n]
            if Xₙ_type <: Static # FIXME; what to do here? Dynamic dispatch?
                push!(strd_tup.args, :($gf(x, $n, false)*$st))
            else
                push!(strd_tup.args, :($ifel(isone($gf(s, $n, false)), zero($Xₙ_type), $gf(x, $n, false)*$st)))
            end
        end
    end
    push!(q.args, strd_tup)
    q
end
@inline function VectorizationBase.bytestrides(fb::ForBroadcast)
    p = getfield(fb,:data)
    _bytestrides(ArrayInterface.size(p), ArrayInterface.strides(p), p)
end

struct Product{A,B}
    a::A
    b::B
end
@inline function Base.size(p::Product)
    M = @inbounds size(p.a)[1]
    (M, Base.tail(size(p.b))...)
end
@inline function Base.size(p::Product, i::Integer)
    i == 1 && return @inbounds size(p.a)[1]
    @inbounds size(p.b)[i]
end
@inline Base.length(p::Product) = prod(size(p))
@inline Base.broadcastable(p::Product) = p
@inline numdims(A) = ndims(A) # fallback
@inline numdims(::Type{Product{A,B}}) where {A,B} = numdims(B)
@inline Base.ndims(::Type{Product{A,B}}) where {A,B} = numdims(B)
# This numdims nonsense is a hack to avoid type piracy in defining:
@inline numdims(::Type{B}) where {N, S <: Base.Broadcast.AbstractArrayStyle{N}, B <: Base.Broadcast.Broadcasted{S}} = N

Base.Broadcast._broadcast_getindex_eltype(::Product{A,B}) where {T, A <: AbstractVecOrMat{T}, B <: AbstractVecOrMat{T}} = T
function Base.Broadcast._broadcast_getindex_eltype(p::Product)
    promote_type(
        Base.Broadcast._broadcast_getindex_eltype(p.a),
        Base.Broadcast._broadcast_getindex_eltype(p.b)
    )
end

# recursive_eltype(::Type{A}) where {T, A <: AbstractArray{T}} = T
# recursive_eltype(::Type{NTuple{N,T}}) where {N,T<:Union{Float32,Float64}} = T
# recursive_eltype(::Type{Float32}) = Float32
# recursive_eltype(::Type{Float64}) = Float64
# recursive_eltype(::Type{Tuple{T}}) where {T} = T
# recursive_eltype(::Type{Tuple{T1,T2}}) where {T1,T2} = promote_type(recursive_eltype(T1), recursive_eltype(T2))
# recursive_eltype(::Type{Tuple{T1,T2,T3}}) where {T1,T2,T3} = promote_type(recursive_eltype(T1), recursive_eltype(T2), recursive_eltype(T3))
# recursive_eltype(::Type{Tuple{T1,T2,T3,T4}}) where {T1,T2,T3,T4} = promote_type(recursive_eltype(T1), recursive_eltype(T2), recursive_eltype(T3), recursive_eltype(T4))
# recursive_eltype(::Type{Tuple{T1,T2,T3,T4,T5}}) where {T1,T2,T3,T4,T5} = promote_type(recursive_eltype(T1), recursive_eltype(T2), recursive_eltype(T3), recursive_eltype(T4), recursive_eltype(T5))

# function recursive_eltype(::Type{Broadcasted{S,A,F,ARGS}}) where {S,A,F,ARGS}
#     recursive_eltype(ARGS)
# end

"""
    A *ˡ B

A lazy product of `A` and `B`. While functionally identical to `A * B`, this may avoid the
need for intermediate storage for any computations in `A` or `B`.  Example:

    @avx @. a + B *ˡ (c + d')

which is equivalent to

     a .+ B * (c .+ d')

It should only be used inside an `@avx` block, and to materialize the result it cannot be
the final operation.
"""
@inline *ˡ(a::A, b::B) where {A,B} = Product{A,B}(a, b)
@inline Base.Broadcast.broadcasted(::typeof(*ˡ), a::A, b::B) where {A, B} = Product{A,B}(a, b)
# TODO: Need to make this handle A or B being (1 or 2)-D broadcast objects.
function add_broadcast!(
    ls::LoopSet, mC::Symbol, bcname::Symbol, loopsyms::Vector{Symbol},
    @nospecialize(prod::Type{<:Product}), elementbytes::Int
)
    A, B = prod.parameters
    Krange = gensym!(ls, "K")
    Klen = gensym!(ls, "K")
    mA = gensym!(ls, "Aₘₖ")
    mB = gensym!(ls, "Bₖₙ")
    gf = GlobalRef(Core,:getfield)
    pushprepreamble!(ls, Expr(:(=), mA, Expr(:call, :forbroadcast, Expr(:(.), bcname, QuoteNode(:a)))))
    pushprepreamble!(ls, Expr(:(=), mB, Expr(:call, :forbroadcast, Expr(:(.), bcname, QuoteNode(:b)))))
    pushprepreamble!(ls, Expr(:(=), Klen, Expr(:call, gf, Expr(:call, :size, mB), 1, false)))
    pushpreamble!(ls, Expr(:(=), Krange, Expr(:call, :(:), staticexpr(1), Klen)))
    k = gensym!(ls, "k")
    add_loop!(ls, Loop(k, 1, Klen, 1, Krange, Klen), k)
    m = loopsyms[1];
    if numdims(B) == 1
        bloopsyms = Symbol[k]
        cloopsyms = Symbol[m]
        reductdeps = Symbol[m, k]
        kvec = bloopsyms
    elseif numdims(B) == 2
        n = loopsyms[2];
        bloopsyms = Symbol[k,n]
        cloopsyms = Symbol[m,n]
        reductdeps = Symbol[m, k, n]
        kvec = Symbol[k]
    else
        throw("B must be a vector or matrix.")
    end
    # load A
    # loadA = add_load!(ls, gensym!(ls, :A), productref(A, mA, m, k), elementbytes)
    loadA = add_broadcast!(ls, gensym!(ls, "A"), mA, Symbol[m,k], A, elementbytes)
    # load B
    loadB = add_broadcast!(ls, gensym!(ls, "B"), mB, bloopsyms, B, elementbytes)
    # set Cₘₙ = 0
    # setC = add_constant!(ls, zero(promote_type(recursive_eltype(A), recursive_eltype(B))), cloopsyms, mC, elementbytes)
    # targetC will be used for reduce_to_add
    mCt = gensym!(ls, mC)
    targetC = add_constant!(ls, gensym!(ls, "zero"), cloopsyms, mCt, elementbytes, :numericconstant)
    push!(ls.preamble_zeros, (identifier(targetC), IntOrFloat))
    setC = add_constant!(ls, gensym!(ls, "zero"), cloopsyms, mC, elementbytes, :numericconstant)
    push!(ls.preamble_zeros, (identifier(setC), IntOrFloat))
    setC.reduced_children = kvec
    # compute Cₘₙ += Aₘₖ * Bₖₙ
    instrsym = Base.libllvm_version < v"11.0.0" ? :vfmadd231 : :vfmadd
    reductop = Operation(
        ls, mC, elementbytes, instrsym, compute, reductdeps, kvec, Operation[loadA, loadB, setC]
    )
    reductop = pushop!(ls, reductop, mC)
    reductfinal = Operation(
        ls, mCt, elementbytes, :reduce_to_add, compute, cloopsyms, kvec, Operation[reductop, targetC]
    )
    pushop!(ls, reductfinal, mCt)
end

function extract_all_1_array!(ls::LoopSet, bcname::Symbol, N::Int, elementbytes::Int)
    refextract = gensym!(ls, bcname)
    ref = Expr(:ref, bcname);
    for _ ∈ 1:N
        push!(ref.args, :begin)
    end
    pushprepreamble!(ls, Expr(:(=), refextract, :(@inbounds $ref)))
    return add_constant!(ls, refextract, elementbytes) # or replace elementbytes with sizeof(T) ? u
end
function doaddref!(ls::LoopSet, op::Operation)
    push!(ls.syms_aliasing_refs, name(op))
    push!(ls.refs_aliasing_syms, op.ref)
    op
end

function add_broadcast!(
    ls::LoopSet, destname::Symbol, bcname::Symbol, loopsyms::Vector{Symbol},
    @nospecialize(LDA::Type{LowDimArray{D,T,N,A}}), elementbytes::Int
) where {D,T,N,A}
    # D,T,N::Int,_ = LDA.parameters
    Dlen = length(D)
    if Dlen == N && !any(D) # array is a scalar, as it is broadcasted on all dimensions
        return extract_all_1_array!(ls, bcname, N, elementbytes)
    end
    fulldims = Symbol[loopsyms[n] for n ∈ 1:N if ((Dlen < n) || D[n]::Bool)]
    ref = ArrayReference(bcname, fulldims)
    loadop = add_simple_load!(ls, destname, ref, elementbytes, true, true )::Operation
    doaddref!(ls, loadop)
end
function add_broadcast_adjoint_array!(
    ls::LoopSet, destname::Symbol, bcname::Symbol, loopsyms::Vector{Symbol}, ::Type{A}, elementbytes::Int
) where {T,N,A<:AbstractArray{T,N}}
    # parent = gensym!(ls, "parent")
    # pushprepreamble!(ls, Expr(:(=), parent, Expr(:call, :parent, bcname)))
    # isone(length(loopsyms)) && return extract_all_1_array!(ls, bcname, N, elementbytes)
    ref = ArrayReference(bcname, Symbol[loopsyms[N + 1 - n] for n ∈ 1:N])
    loadop = add_simple_load!( ls, destname, ref, elementbytes, true, true )::Operation
    doaddref!(ls, loadop)
end
function add_broadcast_adjoint_array!(
    ls::LoopSet, destname::Symbol, bcname::Symbol, loopsyms::Vector{Symbol}, ::Type{<:AbstractVector}, elementbytes::Int
)
    # isone(length(loopsyms)) && return extract_all_1_array!(ls, bcname, N, elementbytes)
    # parent = gensym!(ls, "parent")
    # pushprepreamble!(ls, Expr(:(=), parent, Expr(:call, :parent, bcname)))

    ref = ArrayReference(bcname, Symbol[loopsyms[2]])
    loadop = add_simple_load!( ls, destname, ref, elementbytes, true, true )::Operation
    doaddref!(ls, loadop)
end
function add_broadcast!(
    ls::LoopSet, destname::Symbol, bcname::Symbol, loopsyms::Vector{Symbol},
    ::Type{Adjoint{T,A}}, elementbytes::Int
) where {T, A <: AbstractArray{T}}
    add_broadcast_adjoint_array!( ls, destname, bcname, loopsyms, A, elementbytes )
end
function add_broadcast!(
    ls::LoopSet, destname::Symbol, bcname::Symbol, loopsyms::Vector{Symbol},
    ::Type{Transpose{T,A}}, elementbytes::Int
) where {T, A <: AbstractArray{T}}
    add_broadcast_adjoint_array!( ls, destname, bcname, loopsyms, A, elementbytes )
end
function add_broadcast!(
    ls::LoopSet, destname::Symbol, bcname::Symbol, loopsyms::Vector{Symbol},
    ::Type{<:AbstractArray{T,N}}, elementbytes::Int
) where {T,N}
    loadop = add_simple_load!(ls, destname, ArrayReference(bcname, @view(loopsyms[1:N])), elementbytes, true, true)
    doaddref!(ls, loadop)
end
function add_broadcast!(
    ls::LoopSet, ::Symbol, bcname::Symbol, loopsyms::Vector{Symbol}, ::Type{T}, elementbytes::Int
) where {T<:Number}
    add_constant!(ls, bcname, elementbytes) # or replace elementbytes with sizeof(T) ? u
end
function add_broadcast!(
    ls::LoopSet, ::Symbol, bcname::Symbol, loopsyms::Vector{Symbol}, ::Type{Base.RefValue{T}}, elementbytes::Int
) where {T}
    refextract = gensym!(ls, bcname)
    pushprepreamble!(ls, Expr(:(=), refextract, Expr(:ref, bcname)))
    add_constant!(ls, refextract, elementbytes) # or replace elementbytes with sizeof(T) ? u
end
function add_broadcast!(
    ls::LoopSet, destname::Symbol, bcname::Symbol, loopsyms::Vector{Symbol},
    @nospecialize(_::Type{SubArray{T,N,A,S,B}}), elementbytes::Int
) where {T,N,N2,A<:AbstractArray{T,N2},B,N3,S <: Tuple{Int,Vararg{Any,N3}}}
    inds = Vector{Symbol}(undef, N+1)
    inds[1] = DISCONTIGUOUS
    inds[2:end] .= @view(loopsyms[1:N])
    loadop = add_simple_load!(ls, destname, ArrayReference(bcname, inds), elementbytes, true, true)
    doaddref!(ls, loadop)
end
const BroadcastedArray{S<:Broadcast.AbstractArrayStyle,F,A} = Broadcasted{S,Nothing,F,A}
function add_broadcast!(
    ls::LoopSet, destname::Symbol, bcname::Symbol, loopsyms::Vector{Symbol},
    @nospecialize(B::Type{<:BroadcastedArray}),
    elementbytes::Int
)
    S,_,F,A = B.parameters
    instr = get(FUNCTIONSYMBOLS, F) do
        f = gensym!(ls, "func")
        pushprepreamble!(ls, Expr(:(=), f, Expr(:(.), bcname, QuoteNode(:f))))
        Instruction(bcname, f)
    end
    args = A.parameters
    Nargs = length(args)
    bcargs = gensym!(ls, "bcargs")
    pushprepreamble!(ls, Expr(:(=), bcargs, Expr(:(.), bcname, QuoteNode(:args))))
    # this is the var name in the loop
    parents = Operation[]
    deps = Symbol[]
    # reduceddeps = Symbol[]
    gf = GlobalRef(Core, :getfield)
    for (i,arg) ∈ enumerate(args)
        argname = gensym!(ls, "arg")
        pushprepreamble!(ls, Expr(:(=), argname, Expr(:call, :forbroadcast, Expr(:call, gf, bcargs, i, false))))
        # dynamic dispatch
        parent = add_broadcast!(ls, gensym!(ls, "temp"), argname, loopsyms, arg, elementbytes)::Operation
        push!(parents, parent)
        mergesetdiffv!(deps, loopdependencies(parent), reduceddependencies(parent))
    end
    op = Operation(
        length(operations(ls)), destname, elementbytes, instr, compute, deps, NOPARENTS, parents
    )
    pushop!(ls, op, destname)
end

function add_broadcast_loops!(ls::LoopSet, loopsyms::Vector{Symbol}, destsym::Symbol)
    axes_tuple = Expr(:tuple)
    pushpreamble!(ls, Expr(:(=), axes_tuple, Expr(:call, :axes, destsym)))
    for (n,itersym) ∈ enumerate(loopsyms)
        Nrange = gensym!(ls, "N")
        Nlower = gensym!(ls, "N")
        Nupper = gensym!(ls, "N")
        Nlen = gensym!(ls, "N")
        add_loop!(ls, Loop(itersym, Nlower, Nupper, 1, Nrange, Nlen), itersym)
        push!(axes_tuple.args, Nrange)
        pushpreamble!(ls, Expr(:(=), Nlower, Expr(:call, lv(:maybestaticfirst), Nrange)))
        pushpreamble!(ls, Expr(:(=), Nupper, Expr(:call, lv(:maybestaticlast), Nrange)))
        pushpreamble!(ls, Expr(:(=), Nlen, Expr(:call, lv(:maybestaticlength), Nrange)))
    end
end
# size of dest determines loops
# function vmaterialize!(
@generated function vmaterialize!(
    dest::AbstractArray{T,N}, bc::BC, ::Val{Mod}, ::Val{UNROLL}
) where {T <: NativeTypes, N, BC <: Union{Broadcasted,Product}, Mod, UNROLL}
    # 2+1
    # we have an N dimensional loop.
    # need to construct the LoopSet
    # @show typeof(dest)
    ls = LoopSet(Mod)
    inline, u₁, u₂, isbroadcast, W, rs, rc, cls, l1, l2, l3, threads = UNROLL
    set_hw!(ls, rs, rc, cls, l1, l2, l3)
    ls.isbroadcast[] = isbroadcast # maybe set `false` in a DiffEq-like `@..` macro
    loopsyms = [gensym!(ls, "n") for n ∈ 1:N]
    add_broadcast_loops!(ls, loopsyms, :dest)
    elementbytes = sizeof(T)
    add_broadcast!(ls, :dest, :bc, loopsyms, BC, elementbytes)
    storeop = add_simple_store!(ls, :dest, ArrayReference(:dest, loopsyms), elementbytes)
    doaddref!(ls, storeop)
    resize!(ls.loop_order, num_loops(ls)) # num_loops may be greater than N, eg Product
    # return ls
    Expr(:block, Expr(:meta,:inline), setup_call(ls, :(Base.Broadcast.materialize!(dest, bc)), LineNumberNode(0), inline, false, u₁, u₂, threads%Int), :dest)
end
@generated function vmaterialize!(
    dest′::Union{Adjoint{T,A},Transpose{T,A}}, bc::BC, ::Val{Mod}, ::Val{UNROLL}
) where {T <: NativeTypes, N, A <: AbstractArray{T,N}, BC <: Union{Broadcasted,Product}, Mod, UNROLL}
    # we have an N dimensional loop.
    # need to construct the LoopSet
    ls = LoopSet(Mod)
    inline, u₁, u₂, isbroadcast, W, rs, rc, cls, l1, l2, l3, threads = UNROLL
    set_hw!(ls, rs, rc, cls, l1, l2, l3)
    ls.isbroadcast[] = isbroadcast # maybe set `false` in a DiffEq-like `@..` macro
    loopsyms = [gensym!(ls, "n") for n ∈ 1:N]
    pushprepreamble!(ls, Expr(:(=), :dest, Expr(:call, :parent, :dest′)))
    add_broadcast_loops!(ls, loopsyms, :dest′)
    elementbytes = sizeof(T)
    add_broadcast!(ls, :dest, :bc, loopsyms, BC, elementbytes)
    storeop = add_simple_store!(ls, :dest, ArrayReference(:dest, reverse(loopsyms)), elementbytes)
    doaddref!(ls, storeop)
    resize!(ls.loop_order, num_loops(ls)) # num_loops may be greater than N, eg Product
    Expr(:block, Expr(:meta,:inline), setup_call(ls, :(Base.Broadcast.materialize!(dest′, bc)), LineNumberNode(0), inline, false, u₁, u₂, threads%Int), :dest′)
end
# these are marked `@inline` so the `@avx` itself can choose whether or not to inline.
@generated function vmaterialize!(
    dest::AbstractArray{T,N}, bc::Broadcasted{Base.Broadcast.DefaultArrayStyle{0},Nothing,typeof(identity),Tuple{T2}}, ::Val{Mod}, ::Val{UNROLL}
) where {T <: NativeTypes, N, T2 <: Number, Mod, UNROLL}
    inline, u₁, u₂, isbroadcast, W, rs, rc, cls, l1, l2, l3, threads = UNROLL
    quote
        $(Expr(:meta,:inline))
        arg = T(first(bc.args))
        @avx inline=$inline unroll=($u₁,$u₂) thread=$threads for i ∈ eachindex(dest)
            dest[i] = arg
        end
        dest
    end
end
@generated function vmaterialize!(
    dest′::Union{Adjoint{T,A},Transpose{T,A}}, bc::Broadcasted{Base.Broadcast.DefaultArrayStyle{0},Nothing,typeof(identity),Tuple{T2}}, ::Val{Mod}, ::Val{UNROLL}
) where {T <: NativeTypes, N, A <: AbstractArray{T,N}, T2 <: Number, Mod, UNROLL}
    inline, u₁, u₂, isbroadcast, W, rs, rc, cls, l1, l2, l3, threads = UNROLL
    quote
        $(Expr(:meta,:inline))
        arg = T(first(bc.args))
        dest = parent(dest′)
        @avx inline=$inline unroll=($u₁,$u₂) thread=$threads for i ∈ eachindex(dest)
            dest[i] = arg
        end
        dest′
    end
end

@inline function vmaterialize(
    bc::Broadcasted, ::Val{Mod}, ::Val{UNROLL}
) where {Mod,UNROLL}
    ElType = Base.Broadcast.combine_eltypes(bc.f, bc.args)
    dest = similar(bc, ElType)
    vmaterialize!(dest, bc, Val{Mod}(), Val{UNROLL}())
end

vmaterialize!(dest, bc, ::Val, ::Val, ::StaticInt, ::StaticInt, ::StaticInt) = Base.Broadcast.materialize!(dest, bc)


text:
# Tools for dealing with PCI
const DRV_IS_PCI_VENDOR_ID_INTEL = 0x8086
const VENDOR_ID_MASK = 0x0000_FFFF
const DEVICE_ID_MASK = 0xFFFF_0000
const DEVICE_ID_BITSHIFT = 16
const PCI_ENABLE = 0x8000_0000

# These are found in uncore performance monitoring guide
# Under: "Uncore Performance Monitoring State in PCICFG space"
# Device and Function IDS for Skylake based IMCs
const SKYLAKE_IMC_REGISTERS = (
    # IMC 0 - Channels 0, 1, 2
    ((device = 10, fn = 2), (device = 10, fn = 6), (device = 11, fn = 2)),
    # IMC 1 - Channels 0, 1, 2
    ((device = 12, fn = 2), (device = 12, fn = 6), (device = 13, fn = 2)),
)

const IMC_DEVICE_IDS = (0x2042, 0x2046, 0x204A)

# This is taken from the PCM code
#
# TODO: Dynamically find number of sockets
#
# The general idea is that we know the device and function numbers of IMC devices
# Furthermore, we know the device IDS
#
# So, we enumerate all of the PCI busses until
#   1. We find a valid path
#   2. We read the vendor ID and device ID from the device
#       - Vendor ID must match Intel (0x8086)
#       - Device ID must match one of the iMC device IDS
#
# When we've found this, we've found the bus for the socket.
#
# NOTE: I'm assuming that lower bus numbers correspond to lower sockets.
# This appears consistent with what's happening in PCM.
#
# NOTE: There's a lot of business with PCI group numbers.
# I'm ignoring that for now because we only have group 0 on our system.
#
# TODO: We can abstract this for different architectures using dispatch, but I'm not
# too worried about that at the moment.
function findbusses(;
    device = first(first(SKYLAKE_IMC_REGISTERS)).device,
    fn = first(first(SKYLAKE_IMC_REGISTERS)).fn,
    device_ids = IMC_DEVICE_IDS,
)
    socket_to_bus = UInt[]
    bus_numbers = 0:255
    for bus in bus_numbers
        # Check to see if the path exists
        path = pcipath(bus, device, fn)
        ispath(path) || continue
        pci = Handle(path)

        # Read the first value from the bus - compare against vendor
        value = read(pci, UInt32, IndexZero(0))
        vendor_id = value & VENDOR_ID_MASK
        device_id = (value & DEVICE_ID_MASK) >> DEVICE_ID_BITSHIFT

        # Check if this is run by Intel
        vendor_id == DRV_IS_PCI_VENDOR_ID_INTEL || continue
        in(device_id, device_ids) || continue
        push!(socket_to_bus, bus)

        close(pci)
    end

    return socket_to_bus
end

#####
##### Ice Lake Stuff
#####

# # MMIO_BASE found at Bus U0, Device 0, Function 1, offset D0h.
# const ICX_IMC_MMIO_BASE_OFFSET = 0xd0
# const ICX_IMC_MMIO_BASE_MASK = 0x1FFFFFFF
#
# # MEM0_BAR found at Bus U0, Device 0, Function 1, offset D8h.
# const ICX_IMC_MMIO_MEM0_OFFSET = 0xd8
# const ICX_IMC_MMIO_MEM_STRIDE = 0xd04
# const ICX_IMC_MMIO_MEM_MASK = 0x7FF
#
# # Each IMC has two channels. But there is addressing for three. Need to
# # determine which two channels are active on the system.
# # The offset starts from 0x22800 with stride 0x4000
# #
# const ICX_IMC_MMIO_CHN_OFFSET = 0x22800
# const ICX_IMC_MMIO_CHN_STRIDE = 0x4000
# # /* IMC MMIO size*/
# const ICX_IMC_MMIO_SIZE = 0x4000
#
# const SERVER_UBOX0_REGISTER_DEV_ADDR = 0
# const SERVER_UBOX0_REGISTER_FUNC_ADDR = 1
#
# function map_imc_pmon(package_id::Integer, imc_index::Integer, channel_index::Integer; kw...)
#     return map_imc_pmon(convert(UInt, package_id), convert(UInt, imc_index), convert(UInt, channel_index), kw...)
# end
#
# function map_imc_pmon(package_id::UInt64, imc_index::UInt64, channel_index::UInt64; device_ids = (0x3451,))
#     bus_numbers = 0:255
#     device = SERVER_UBOX0_REGISTER_DEV_ADDR
#     fn = SERVER_UBOX0_REGISTER_FUNC_ADDR
#
#     socket_to_bus = UInt[]
#
#     for bus in bus_numbers
#         path = pcipath(bus, device, fn)
#         ispath(path) || continue
#         pci = Handle(path)
#         #seek(pci, IndexZero(0))
#         #value = read(pci, UInt32)
#         value = read(pci, UInt32, IndexZero(0))
#
#         vendor_id = value & VENDOR_ID_MASK
#         device_id = (value & DEVICE_ID_MASK) >> DEVICE_ID_BITSHIFT
#
#         vendor_id == DRV_IS_PCI_VENDOR_ID_INTEL || continue
#         in(device_id, device_ids) || continue
#         push!(socket_to_bus, bus)
#
#         # # Do some things
#         # # Read MEMn addr (51:23) from MMIO_BASE register
#         # seek(pci, IndexZero(ICX_IMC_MMIO_BASE_OFFSET))
#         # pci_uint32 = read(pci, UInt32, IndexZero(ICX_IMC_MMIO_MEM0_OFFSET)
#         # address = (pci_uint32 & ICX_IMC_MMIO_BASE_MASK) << 23
#
#         # # read MEMn addr (22:12) from MEMn_BAR register
#         # mem_offset = ICX_IMC_MMIO_MEM0_OFFSET + imc_index * ICX_IMC_MMIO_MEM_STRIDE;
#         # seek(pci, IndexZero(mem_offset))
#         # pci_uint32 = read(pci, UInt32)
#         # address |= (pci_uint32 & ICX_IMC_MMIO_MEM_MASK) << 12;
#
#         # # IMC PMON registers start from PMONUNITCTRL */
#         # address += ICX_IMC_MMIO_CHN_OFFSET + channel_index * ICX_IMC_MMIO_CHN_STRIDE;
#         # @show address
#         close(pci)
#     end
#     return socket_to_bus
# end
#
# # /*
# # * pkg_id: Socket id
# # * imc_idx: The IMC index
# # * channel_idx: The channel index
# # */
# # Void *map_imc_pmon(int pkg_id, int imc_idx, int channel_idx)
# # {
# #   struct pci_dev *pdev = NULL;
# #   resource_size_t addr;
# #   u32 pci_dword;
# #   void *io_addr;
# #   int mem_offset;
# #   /*
# #   * Device ID of Bus U0, Device 0, Function 1 is 0x3451 */
# #   * Get its pdev on the specific socket.
# #   */
# #   while(1){
# #       pdev = pci_get_device(PCI_VENDOR_ID_INTEL, 0x3451, pdev);
# #       if ((!pdev) || (pdev->bus == UNC_UBOX_package_to_bus_map[pkg_id]))
# #           break;
# #   }
# #   if (!pdev)
# #       return NULL;
# #
# #   /* read MEMn addr (51:23) from MMIO_BASE register */
# #   pci_read_config_dword(pdev, ICX_IMC_MMIO_BASE_OFFSET, &pci_dword);
# #   addr = (pci_dword & ICX_IMC_MMIO_BASE_MASK) << 23;
# #
# #   /* read MEMn addr (22:12) from MEMn_BAR register */
# #   mem_offset = ICX_IMC_MMIO_MEM0_OFFSET + mem_idx * ICX_IMC_MMIO_MEM_STRIDE;
# #   pci_read_config_dword(pdev, mem_offset, &pci_dword);
# #   addr |= (pci_dword & ICX_IMC_MMIO_MEM_MASK) << 12;
#
#   /* IMC PMON registers start from PMONUNITCTRL */
#   addr += ICX_IMC_MMIO_CHN_OFFSET + channel_idx * ICX_IMC_MMIO_CHN_STRIDE;
#   /* map the IMC PMON registers */
#   io_addr = ioremap(addr, ICX_IMC_MMIO_SIZE);
#   return io_addr;
# }


text:
using GLM
using StatsBase: RegressionModel
import Base: show

# Similar to typealiases in GLM
typealias BlasReal Union{Float32,Float64}
typealias FP AbstractFloat
typealias FPVector{T<:FP} DenseArray{T,1}

# Initial code imported from pull request https://github.com/JuliaStats/GLM.jl/pull/70/
# See also https://github.com/JuliaStats/GLM.jl/pull/65

effects(mod::RegressionModel) = effects(mod.model)
effects{T<:BlasReal,V<:FPVector}(
  mod::LinearModel{LmResp{V},DensePredQR{T}}) =
    (mod.pp.qr[:Q]'*mod.rr.y)[1:size(mod.pp.X,2)]

function effects{T<:BlasReal,V<:FPVector,C<:Any}(
  mod::LinearModel{LmResp{V},DensePredChol{T,C}})
  # TODO this is hackish, and recomputes Q1 with every call
  R = cholfact!(mod.pp)[:U]
  Q1 = mod.pp.X / R
  return Q1' * mod.rr.y
end

type AnovaTable{T<:BlasReal}
  responsename
  termnames::Vector{Any}
  DF::Vector{T}
  SS::Vector{T}
  MS::Vector{T}
  F::Vector{T}
  PrF::Vector{T}
end

function anova(mod::RegressionModel)
  eff = effects(mod)
  with_intercept = all(mod.model.pp.X[:,1] .== 1.0)
  T = eltype(mod.model.pp.X)
  response = mod.mf.terms.eterms[1]
  termnames = mod.mf.terms.terms
  assign = mod.mm.assign
  df_residual = dof_residual(mod)
  devnc = deviance(mod)
  return _anova(eff,with_intercept,T,response,termnames,assign,df_residual,devnc)
end
function anova{M<:LinearModel}(mod::M)
  eff = effects(mod)
  with_intercept = all(mod.pp.X[:,1] .== 1.0)
  T = eltype(mod.pp.X)
  response = "y"
  if with_intercept
    termnames = ["x$i" for i in 1:size(mod.pp.X,2)-1]
  else
    termnames = ["x$i" for i in 1:size(mod.pp.X,2)]
  end
  assign = 0:1:size(mod.pp.X,2)-1
  df_residual = dof_residual(mod)
  devnc = deviance(mod)
  return _anova(eff,with_intercept,T,response,termnames,assign,df_residual,devnc)
end

function _anova(eff,with_intercept,T,response,termnames,assign,df_residual,devnc)
  k = with_intercept ? 1 : 0
  unq_assign = unique(assign)
  n = length(unq_assign) + 1 - k
  DF = zeros(T,n)
  SS = zeros(T,n)
  MS = zeros(T,n)
  fstat = zeros(T,n-1)
  pval = zeros(T,n-1)
  DF[end] = df_residual
  SS[end] = devnc
  MS[end] = SS[end]/DF[end]

  @inbounds for i in 1:n-1
    v = unq_assign[i+k]
    mask = assign .== v
    DF[i] = sum(mask)
    SS[i] = sumabs2(eff[mask])
    MS[i] = SS[i]/DF[i]
    fstat[i] = MS[i]/MS[end]
    pval[i] = ccdf(FDist(DF[i], DF[end]), fstat[i])
  end
  return AnovaTable(response,termnames,DF,SS,MS,fstat,pval)
end

function show(io::IO, tab::AnovaTable)
  rou(x) = round(x,3)
  pad(x,w) = string(x,RepString(" ",w-length(x)))

  println(io,"Analysis of Variance Table")
  println(io,"")
  println(io,"Response: ",tab.responsename)
  w = maximum(length.(string.(tab.termnames)))
  resd = "Residuals"
  w = max(w,length(resd))
  println(io,RepString(" ",w+3),'\t',"DF",'\t',"SS",'\t',"MS",'\t',
    "F",'\t',"Pr(>F)")
  for i in 1:endof(tab.termnames)
    println(io,pad(string(tab.termnames[i]),w),'\t',rou(tab.DF[i]),'\t',rou(tab.SS[i]),
      '\t',rou(tab.MS[i]),'\t',rou(tab.F[i]),'\t',rou(tab.PrF[i]))
  end
  i = length(tab.DF)
  println(io,pad(resd,w),'\t',rou(tab.DF[i]),'\t',rou(tab.SS[i]),
    '\t',rou(tab.MS[i]))
end


text:
include("Styles.jl")

include("game.jl")
import .Tetris

logo = raw"""

                    _____ ___ _____ ___ ___ ___
                   |_   _| __|_   _| _ \_ _/ __|
                     | | | _|  | | |   /| |\__ \
        ___ ___  __ _| |_|___| |_|_|_|_\___|___/ ___ ___ ___
        | __| _ \/ _ \|  \/  | |_   _| __|_   _| _ \_ _/ __|
        | _||   / (_) | |\/| |   | | | _|  | | |   /| |\__ \
        |_| |_|_\\___/|_|  |_|   |_| |___| |_| |_|_\___|___/

"""

board_string = """
            ______________________
            |%|     SCORE:  &
            |%|     LEVEL:  !
            |%|     RECORD: @
            |%|
            |%|
            |%|     NEXT TETROMINO:
            |%|      ___________
            |%|      |*|
            |%|      |*|
            |%|      |*|
            |%|      |*|
            |%|      ‾‾‾‾‾‾‾‾‾‾‾
            |%|
            |%|     CONTROLS:
            |%|     left: a
            |%|     down (1 block): s
            |%|     right: d
            |%|     rotate: w
            |%|     down (all blocks): l
            |%|     arrow keys also work
            |%|
            ‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾
"""

const BLOCK_CHAR =  "██"
const WHITE_BLOCK = " ⋅"

#                 blank    I     J      L       O      S         T       z
const COLORS = [CWHITE, CWHITE, CBLUE, CRED, CYELLOW, CGREEN, CVIOLET, CBEIGE]

function format_board_line(board_line)
    line = CBLACKBG * if length(board_line)<=4  WHITE_BLOCK ^ (4-length(board_line)) else "" end

    for col in board_line
        fill_char = if col == 0 WHITE_BLOCK else BLOCK_CHAR end

        line *= "$(COLORS[col+1])$(fill_char)"
    end
    return line * CEND
end

function print_board(args)
    body
end

print_board(gs::Tetris.GameStatus) = print_board(gs.board, gs.next_tetromino, gs.score,gs.level,gs.record)

function print_board(board, next_tetromino=[], score="None", level="None", record="None")
    next_tetromino = next_tetromino.type

    run(`clear`)
    board_line_counter = 1
    next_tetromino_line_counter = 1
    output_text = logo

    for line in split(board_string, '\n')
        if '%' in line
            line = replace(line, '%'=>format_board_line(board[board_line_counter,:]))
            board_line_counter += 1
        end

        if '*' in line
            if next_tetromino_line_counter <= length(next_tetromino[:,1])
                line = replace(line, '*'=>format_board_line(next_tetromino[next_tetromino_line_counter,:]))
                next_tetromino_line_counter += 1
            else
                line = replace(line, '*'=>format_board_line([0,0,0,0]))
            end
        end

        line = replace(line, '&'=>score)
        line = replace(line, '!'=>level)
        line = replace(line, '@'=>record)
        output_text *= line * '\n'
    end
    # write(stdout, output_text)
    # println("PRINT BOARD")
    println(output_text)
end

max_stars_repo_name: GabrielTrettel/tetris

max_issues_repo_name: GabrielTrettel/tetris

max_forks_repo_name: GabrielTrettel/tetris


text:
# ----------- This file implements the abstract one-particle basis interface

function evaluate!(A, basis::OneParticleBasis, cfg::AbstractConfiguration)
   fill!(A, 0)
   for X in cfg
      add_into_A!(A, basis, X)
   end
   return A
end

function evaluate!(A, basis::OneParticleBasis, X::AbstractState)
   fill!(A, 0)
   add_into_A!(A, basis, X)
   return A
end

_check_args_is_sym() = true
_check_args_is_sym(::Symbol) = true

# args... may be empty or a symbol  for partial derivatives
function evaluate_d!(dA, basis::OneParticleBasis,
                     X::Union{AbstractState, AbstractConfiguration},
                     args...)
   @assert _check_args_is_sym(args...)
   A = acquire_B!(basis, X)
   evaluate_ed!(A, dA, basis, X, args...)
   release_B!(basis, A)
   return dA
end

# args... may be empty or a symbol  for partial derivatives
function evaluate_ed!(A, dA, basis::OneParticleBasis,
                      cfg::AbstractConfiguration, args...)
   @assert _check_args_is_sym(args...)
   fill!(A, 0)
   for (j, X) in enumerate(cfg)
      add_into_A_dA!(A, (@view dA[:, j]), basis, X, args...)
   end
   return A, dA
end

# args... may be empty or a symbol for partial derivatives
function evaluate_ed!(A, dA, basis::OneParticleBasis, X::AbstractState, args...)
   @assert _check_args_is_sym(args...)
   fill!(A, 0)
   add_into_A_dA!(A, dA, basis, X, args...)
   return A, dA
end

# fix evaluate_d and evaluate_ed for partial D, i.e. args = a symbol

evaluate_d(basis::ACEBasis, X::Union{AbstractState, AbstractConfiguration}, sym::Symbol) =
      evaluate_d!( acquire_dB!(basis, X), basis, X, sym)

evaluate_ed(basis::ACEBasis, X::Union{AbstractState, AbstractConfiguration}, sym::Symbol) =
      evaluate_ed!( acquire_B!(basis, X), acquire_dB!(basis, X), basis, X, sym )

function add_into_A_dA!(A, dAj, basis::OneParticleBasis, Xj, sym::Symbol)
   if sym in argsyms(basis)
      add_into_A_dA!(A, dAj, basis::OneParticleBasis, Xj)
   else
      add_into_A!(A, basis, Xj)
      fill!(dAj, zero(eltype(dAj)))
   end
end

# if a basis hasn't implemented this ...
function add_into_A_dA!(A, dAj, basis::OneParticleBasis, Xj)
   ϕ = acquire_B!(basis, Xj)
   evaluate_ed!(ϕ, dAj, basis, Xj)
   @. A .+= ϕ
   release_B!(basis, ϕ)
   return nothing
end

function add_into_A!(A, basis::OneParticleBasis, Xj)
   ϕ = acquire_B!(basis, Xj)
   evaluate!(ϕ, basis, Xj)
   @. A += ϕ
   release_B!(basis, ϕ)
   return nothing
end

# -------------------- AD codes

import ChainRules
import ChainRules: rrule, NoTangent, ZeroTangent

function _rrule_evaluate(basis::OneParticleBasis,
                         cfg::ACEConfig,
                         W )
   dXs = [ _rrule_evaluate(basis, X, W) for X in cfg ]
   return DACEConfig(dXs)
end

function rrule(::typeof(evaluate), basis::OneParticleBasis,
                cfg::AbstractConfiguration)
   return evaluate(basis, cfg),
          w -> (NoTangent(), ZeroTangent(),
                _rrule_evaluate(basis, cfg, w) )
end

# -------------------- END AD codes

# function set_Aindices!(basis::OneParticleBasis)
#    NZ = numz(basis)
#    for iz0 = 1:NZ
#       idx = 0
#       for iz = 1:NZ
#          len = length(basis, iz, iz0)
#          basis.Aindices[iz, iz0] = (idx+1):(idx+len)
#          idx += len
#       end
#    end
#    return basis
# end
#
# get_basis_spec(basis::OneParticleBasis, iz0::Integer) =
#       get_basis_spec(basis, i2z(iz0))
#
# get_basis_spec(basis::OneParticleBasis, s::Symbol) =
#       get_basis_spec(basis, atomic_number(s))

# --------------------

"""
`struct One1p` : the 1-p basis where everything just maps to 1
"""
struct One1pBasis <: OneParticleBasis{Bool}
end

Base.length(::One1pBasis) = 1

function evaluate!(B, basis::One1pBasis, args...)
   B[1] = 1
   return B
end

function evaluate_d!(dB, basis::One1pBasis, args...)
   dB[1] = zero(eltype(dB))
   return dB
end

fltype(::One1pBasis) = Bool

symbols(::One1pBasis) = Symbol[]

indexrange(::One1pBasis) = NamedTuple()

isadmissible(b, ::One1pBasis) = true

degree(b, ::One1pBasis) = 0

get_spec(::One1pBasis) = [(1,)]


text:
n, p = 500, 5
((X, y, θ), (X1, y1, θ1)) = generate_continuous(n, p; seed=525)

# adding some outliers (both positive and negative)
Random.seed!(543)
y1a = outlify(y1, 0.1)

@testset "QuantileReg" begin
    δ = 0.5 # effectively LAD regression
    λ = 1.0
    rr = QuantileRegression(δ, lambda=λ, penalize_intercept=true)
    J = objective(rr, X, y1a)
    o = RobustLoss(Quantile(δ)) + λ * L2Penalty()
    @test J(θ1) ≈ o(y1a, X1*θ1, θ1)
    ls = LinearRegression()
    θ_ls    = fit(ls, X, y1a)
    θ_lbfgs = fit(rr, X, y1a, solver=LBFGS())
    θ_iwls  = fit(rr, X, y1a, solver=IWLSCG())
    @test isapprox(J(θ1),      491.94661, rtol=1e-5)
    @test isapprox(J(θ_ls),    614.70403, rtol=1e-5)  # note that LS is crap due to outliers
    @test isapprox(J(θ_lbfgs), 491.65694, rtol=1e-5)
    @test isapprox(J(θ_iwls),  491.65694, rtol=1e-5)

    # NOTE: newton and newton-cg not available because ϕ = 0 identically
    # will throw an error if called.
    @test_throws ErrorException fit(rr, X, y1, solver=Newton())
    @test_throws ErrorException fit(rr, X, y1, solver=NewtonCG())

    # don't penalize intercept
    rr = QuantileRegression(δ, lambda=λ)
    J = objective(rr, X, y1a)
    ls = LinearRegression()
    θ_ls    = fit(ls, X, y1a)
    θ_lbfgs = fit(rr, X, y1a, solver=LBFGS())
    θ_iwls  = fit(rr, X, y1a, solver=IWLSCG())
    @test isapprox(J(θ1),      489.48936, rtol=1e-5)
    @test isapprox(J(θ_ls),    612.16417, rtol=1e-5)  # note that LS is crap due to outliers
    @test isapprox(J(θ_lbfgs), 489.21223, rtol=1e-5)
    @test isapprox(J(θ_iwls),  489.22193, rtol=1e-5)

    if DO_COMPARISONS
        # Compare with R's QuantReg package
        # NOTE: QuantReg doesn't allow for penalties so re-fitting with λ=0
        rr = QuantileRegression(δ, lambda=0)
        J  = objective(rr, X, y1a)
        θ_ls     = fit(LinearRegression(), X, y1a)
        θ_lbfgs  = fit(rr, X, y1a, solver=LBFGS())
        θ_iwls   = fit(rr, X, y1a, solver=IWLSCG())
        θ_qr_br  = rcopy(QUANTREG.rq_fit_br(X1, y1a))[:coefficients]
        θ_qr_fnb = rcopy(QUANTREG.rq_fit_fnb(X1, y1a))[:coefficients]
        # NOTE: we take θ_qr_br as reference point
        @test isapprox(J(θ_ls), 610.41023,  rtol=1e-5)
        @test J(θ_qr_br) ≈      486.36730 # <- ref value
        # Their IP algorithm essentially gives the same answer
        @test (J(θ_qr_fnb) - J(θ_qr_br)) ≤ 1e-10
        # Our algorithms are close enough
        @test isapprox(J(θ_lbfgs), 486.36782, rtol=1e-5)
        @test isapprox(J(θ_iwls),  486.37331, rtol=1e-5)
    end
end

###########################
## With Sparsity penalty ##
###########################

n, p = 500, 100
((X, y, θ), (X1, y1, θ1)) = generate_continuous(n, p;  seed=51112, sparse=0.1)
# pepper with outliers
y1a  = outlify(y1, 0.1)

@testset "LAD+L1" begin
    λ = 5.0
    γ = 10.0
    rr = LADRegression(λ, γ; penalize_intercept=true)
    J  = objective(rr, X, y1a)
    θ_ls    = X1 \ y1a
    θ_fista = fit(rr, X, y1a, solver=FISTA())
    θ_ista  = fit(rr, X, y1a, solver=ISTA())
    @test isapprox(J(θ_ls),    943.33942, rtol=1e-5)
    @test isapprox(J(θ_fista), 390.53177, rtol=1e-5)
    @test isapprox(J(θ_ista),  390.53177, rtol=1e-3)
    @test nnz(θ_fista) == 55
    @test nnz(θ_ista)  == 51 # interesting, things look a bit unstable?

    if DO_COMPARISONS
        # Compare with R's QuantReg package
        # NOTE: QuantReg doesn't apply the penalty on the intercept
        rr = LADRegression(5.0; penalty=:l1)
        J  = objective(rr, X, y1a)
        θ_ls       = X1 \ y1a
        θ_fista    = fit(rr, X, y1a, solver=FISTA())
        θ_ista     = fit(rr, X, y1a, solver=ISTA())
        θ_qr_lasso = rcopy(QUANTREG.rq_fit_lasso(X1, y1a))[:coefficients]
        # NOTE: we take θ_qr_br as reference point
        @test isapprox(J(θ_ls),       799.18992, rtol=1e-5)
        @test isapprox(J(θ_qr_lasso), 373.15568, rtol=1e-5) # <- ref value
        # Our algorithms are close enough
        @test isapprox(J(θ_fista),    372.58504, rtol=1e-5)
        @test isapprox(J(θ_ista),     372.91354, rtol=1e-5)
        # in this case we do a fair bit better
        @test nnz(θ_qr_lasso) == 101
        @test nnz(θ_fista)    == 89
        @test nnz(θ_ista)     == 82
        # in this case fista is best
        @test J(θ_fista) < J(θ_ista) < J(θ_qr_lasso)
    end
end


text:
"""
    HubbardModelAttractive(lattice; params...)
    HubbardModelAttractive(L, dims; params...)
    HubbardModelAttractive(params::Dict)
    HubbardModelAttractive(params::NamedTuple)
    HubbardModelAttractive(; params...)

Defines an attractive (negative `U`) Hubbard model on a given (or derived)
`lattice`. If a linear system size `L` and dimensionality `dims` is given, the
`lattice` will be a Cubic lattice of fitting size.

Additional parameters (keyword arguments) include:
* `l::AbstractLattice = lattice`: The lattice the model uses. The keyword
argument takes precedence over the argument `lattice`.
* `U::Float64 = 1.0 > 0.0` is the absolute value of the Hubbard Interaction.
* `t::Float64 = 1.0` is the hopping strength.
* `mu::Float64` is the chemical potential.

Internally, a discrete Hubbard Stratonovich transformation (Hirsch
transformation) is used in the spin/magnetic channel to enable DQMC. The
resulting Hubbard Stratonovich fiels is real.
Furthermore, we use spin up/down symmetry to speed up the simulation. As a
result the greens matrix is of size (N, N) with N the number of sites, and the
element G[i, j] corresponds to the up-up and down-down element.
"""
@with_kw_noshow struct HubbardModelAttractive{LT<:AbstractLattice} <: HubbardModel
    # user optional
    mu::Float64 = 0.0
    U::Float64 = 1.0
    @assert U >= 0. "U must be positive."
    t::Float64 = 1.0
    l::LT

    # non-user fields
    flv::Int = 1
    # to avoid allocations (TODO always real?)
    IG::Vector{Float64} = Vector{Float64}(undef, length(l))
    G::Vector{Float64} = Vector{Float64}(undef, length(l))
end

HubbardModelAttractive(params::Dict{Symbol}) = HubbardModelAttractive(; params...)
HubbardModelAttractive(params::NamedTuple) = HubbardModelAttractive(; params...)
function HubbardModelAttractive(lattice::AbstractLattice; kwargs...)
    HubbardModelAttractive(l = lattice; kwargs...)
end
function HubbardModelAttractive(L, dims; kwargs...)
    l = choose_lattice(HubbardModelAttractive, dims, L)
    HubbardModelAttractive(l = l; kwargs...)
end

# cosmetics
import Base.summary
import Base.show
Base.summary(model::HubbardModelAttractive) = "attractive Hubbard model"
function Base.show(io::IO, model::HubbardModelAttractive)
    print(io, "attractive Hubbard model, $(length(model.l)) sites")
end
Base.show(io::IO, m::MIME"text/plain", model::HubbardModelAttractive) = print(io, model)

# Convenience
@inline parameters(m::HubbardModelAttractive) = (L = m.L, t = m.t, U = m.U, mu = m.mu)

"""
Calculates the hopping matrix \$T_{i, j}\$ where \$i, j\$ are
site indices.

Note that since we have a time reversal symmetry relating spin-up
to spin-down we only consider one spin sector (one flavor) for the attractive
Hubbard model in the DQMC simulation.

This isn't a performance critical method as it is only used once before the
actual simulation.
"""
function hopping_matrix(mc::DQMC, m::HubbardModelAttractive{L}) where {L<:AbstractLattice}
    N = length(m.l)
    T = diagm(0 => fill(-m.mu, N))

    # Nearest neighbor hoppings
    @inbounds @views begin
        for (src, trg) in neighbors(m.l, Val(true))
            trg == -1 && continue
            T[trg, src] += -m.t
        end
    end

    return T
end

"""
Calculate the interaction matrix exponential `expV = exp(- power * delta_tau * V(slice))`
and store it in `result::Matrix`.

This is a performance critical method.
"""
@inline @bm function interaction_matrix_exp!(mc::DQMC, m::HubbardModelAttractive,
            result::Diagonal, conf::HubbardConf, slice::Int, power::Float64=1.)
    dtau = mc.p.delta_tau
    lambda = acosh(exp(0.5 * m.U * dtau))

    N = size(result, 1)
    @inbounds for i in 1:N
        result.diag[i] = exp(sign(power) * lambda * conf[i, slice])
    end
    nothing
end

@inline @bm function propose_local(
        mc::DQMC, m::HubbardModelAttractive, i::Int, slice::Int, conf::HubbardConf
    )
    # see for example dos Santos Introduction to quantum Monte-Carlo
    greens = mc.s.greens
    dtau = mc.p.delta_tau
    lambda = acosh(exp(m.U * dtau/2))

    @inbounds ΔE_boson = -2. * lambda * conf[i, slice]
    γ = exp(ΔE_boson) - 1
    @inbounds detratio = (1 + γ * (1 - greens[i,i]))^2
    # squared because of two spin sectors

    return detratio, ΔE_boson, γ
end

@inline @bm function accept_local!(
        mc::DQMC, m::HubbardModelAttractive, i::Int, slice::Int, conf::HubbardConf,
        detratio, ΔE_boson, γ)
    greens = mc.s.greens

    # Unoptimized Version
    # u = -greens[:, i]
    # u[i] += 1.
    # greens .-= kron(u * 1./(1 + gamma * u[i]), transpose(gamma * greens[i, :]))
    # conf[i, slice] .*= -1.

    # Optimized
    # copy! and `.=` allocate, this doesn't. Synced loop is marginally faster
    @avx for j in eachindex(m.IG)
        m.IG[j] = -greens[j, i]
        m.G[j] = greens[i, j]
    end
    @inbounds m.IG[i] += 1.0
    # This is way faster for small systems and still ~33% faster at L = 15
    # Also no allocations here
    @inbounds x = γ / (1.0 + γ * m.IG[i])
    @avx for k in eachindex(m.IG), l in eachindex(m.G)
        greens[k, l] -= m.IG[k] * x * m.G[l]
    end
    @inbounds conf[i, slice] *= -1
    nothing
end

"""
Calculate energy contribution of the boson, i.e. Hubbard-Stratonovich/Hirsch field.
"""
@inline function energy_boson(mc::DQMC, m::HubbardModelAttractive)
    hsfiled = conf(mc)
    dtau = mc.p.delta_tau
    lambda = acosh(exp(m.U * dtau/2))
    return lambda * sum(hsfield)
end

function greens(mc::DQMC, model::HubbardModelAttractive)
    G = greens!(mc)
    vcat(hcat(G, zeros(size(G))), hcat(zeros(size(G)), G))
end

function save_model(
        file::JLDFile,
        m::HubbardModelAttractive,
        entryname::String="Model"
    )
    write(file, entryname * "/VERSION", 1)
    write(file, entryname * "/type", typeof(m))

    write(file, entryname * "/mu", m.mu)
    write(file, entryname * "/U", m.U)
    write(file, entryname * "/t", m.t)
    save_lattice(file, m.l, entryname * "/l")
    write(file, entryname * "/flv", m.flv)

    nothing
end

#     load_parameters(data, ::Type{<: DQMCParameters})
#
# Loads a DQMCParameters object from a given `data` dictionary produced by
# `JLD.load(filename)`.
function _load(data, ::Type{T}) where T <: HubbardModelAttractive
    if !(data["VERSION"] == 1)
        throw(ErrorException("Failed to load HubbardModelAttractive version $(data["VERSION"])"))
    end

    l = _load(data["l"], data["l"]["type"])
    data["type"](
        mu = data["mu"],
        U = data["U"],
        t = data["t"],
        l = l,
        flv = data["flv"]
    )
end

################################################################################
### Measurement overloads
################################################################################

# Need some measurement overwrites because nflavors = 1
checkflavors(::HubbardModelAttractive) = nothing

function cdc_kernel(mc, ::HubbardModelAttractive, i, j, G::AbstractArray)
    # spin up and down symmetric, so (i+N, i+N) = (i, i); (i+N, i) drops
    4 * (1 - G[i, i]) * (1 - G[j, j]) + 2 * (I[j, i] - G[j, i]) * G[i, j]
end
function cdc_kernel(mc, ::HubbardModelAttractive, i, j, pg::NTuple{4})
    G00, G0l, Gl0, Gll = pg
    # spin up and down symmetric, so (i+N, i+N) = (i, i); (i+N, i) drops
    4 * (1 - Gll[i,i]) * (1 - G00[j,j]) - 2 * G0l[j,i] * Gl0[i,j]
end

mx_kernel(mc, ::HubbardModelAttractive, i, G::AbstractArray) = 0.0
my_kernel(mc, ::HubbardModelAttractive, i, G::AbstractArray) = 0.0
mz_kernel(mc, ::HubbardModelAttractive, i, G::AbstractArray) = 0.0

sdc_x_kernel(mc, ::HubbardModelAttractive, i, j, G::AbstractArray) = 2(I[j,i] - G[j,i]) * G[i,j]
sdc_y_kernel(mc, ::HubbardModelAttractive, i, j, G::AbstractArray) = 2(I[j,i] - G[j,i]) * G[i,j]
sdc_z_kernel(mc, ::HubbardModelAttractive, i, j, G::AbstractArray) = 2(I[j,i] - G[j,i]) * G[i,j]

sdc_x_kernel(mc, ::HubbardModelAttractive, i, j, pg::NTuple{4}) = -2 * pg[2][j,i] * pg[3][i,j]
sdc_y_kernel(mc, ::HubbardModelAttractive, i, j, pg::NTuple{4}) = -2 * pg[2][j,i] * pg[3][i,j]
sdc_z_kernel(mc, ::HubbardModelAttractive, i, j, pg::NTuple{4}) = -2 * pg[2][j,i] * pg[3][i,j]

function pc_kernel(mc, ::HubbardModelAttractive, src1, trg1, src2, trg2, G::AbstractArray)
    G[src1, src2] * G[trg1, trg2]
end
function pc_kernel(mc, ::HubbardModelAttractive, src1, trg1, src2, trg2, pg::NTuple{4})
    pg[3][src1, src2] * pg[3][trg1, trg2]
end

function cc_kernel(mc, ::HubbardModelAttractive, src1, trg1, src2, trg2, pg::NTuple{4})
    G00, G0l, Gl0, Gll = pg
    N = length(lattice(mc))
    T = mc.s.hopping_matrix

    # up-up counts, down-down counts, mixed only on 11s or 22s
    s1 = src1; t1 = trg1
    s2 = src2; t2 = trg2
    output = 4.0 *
        (T[s1, t1] * Gll[t1, s1] - T[t1, s1] * Gll[s1, t1]) *
        (T[s2, t2] * G00[t2, s2] - T[t2, s2] * G00[s2, t2]) +
        2.0 * T[t1, s1] * T[t2, s2] * (- G0l[s2, t1]) * Gl0[s1, t2] -
        2.0 * T[s1, t1] * T[t2, s2] * (- G0l[s2, s1]) * Gl0[t1, t2] -
        2.0 * T[t1, s1] * T[s2, t2] * (- G0l[t2, t1]) * Gl0[s1, s2] +
        2.0 * T[s1, t1] * T[s2, t2] * (- G0l[t2, s1]) * Gl0[t1, s2]

    output
end


text:
using Faulhaber
using Test

@testset "Faulhaber.jl" begin
    @test_throws ArgumentError faulhaber(10, Val(-7))

    for m in 1:200, p in 0:7, T in [Int, Int128, Float64]
        if T == Float64
            @test faulhaber(float(m), Val(p)) ≈ sum(i->i^p, 1:big(m))
        else
            @test faulhaber(T(m), Val(p)) === T(sum(i->i^p, 1:big(m)))
        end
    end
    for m in 1:13, p in 8:15
        @test faulhaber(m, Val(p)) == sum(i->i^p, 1:big(m))
    end
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiAuditregistrationV1alpha1WebhookClientConfig <: SwaggerModel
    caBundle::Any # spec type: Union{ Nothing, Vector{UInt8} } # spec name: caBundle
    service::Any # spec type: Union{ Nothing, IoK8sApiAuditregistrationV1alpha1ServiceReference } # spec name: service
    url::Any # spec type: Union{ Nothing, String } # spec name: url

    function IoK8sApiAuditregistrationV1alpha1WebhookClientConfig(;caBundle=nothing, service=nothing, url=nothing)
        o = new()
        validate_property(IoK8sApiAuditregistrationV1alpha1WebhookClientConfig, Symbol("caBundle"), caBundle)
        setfield!(o, Symbol("caBundle"), caBundle)
        validate_property(IoK8sApiAuditregistrationV1alpha1WebhookClientConfig, Symbol("service"), service)
        setfield!(o, Symbol("service"), service)
        validate_property(IoK8sApiAuditregistrationV1alpha1WebhookClientConfig, Symbol("url"), url)
        setfield!(o, Symbol("url"), url)
        o
    end
end # type IoK8sApiAuditregistrationV1alpha1WebhookClientConfig

const _property_map_IoK8sApiAuditregistrationV1alpha1WebhookClientConfig = Dict{Symbol,Symbol}(Symbol("caBundle")=>Symbol("caBundle"), Symbol("service")=>Symbol("service"), Symbol("url")=>Symbol("url"))
const _property_types_IoK8sApiAuditregistrationV1alpha1WebhookClientConfig = Dict{Symbol,String}(Symbol("caBundle")=>"Vector{UInt8}", Symbol("service")=>"IoK8sApiAuditregistrationV1alpha1ServiceReference", Symbol("url")=>"String")
Base.propertynames(::Type{ IoK8sApiAuditregistrationV1alpha1WebhookClientConfig }) = collect(keys(_property_map_IoK8sApiAuditregistrationV1alpha1WebhookClientConfig))
Swagger.property_type(::Type{ IoK8sApiAuditregistrationV1alpha1WebhookClientConfig }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiAuditregistrationV1alpha1WebhookClientConfig[name]))}
Swagger.field_name(::Type{ IoK8sApiAuditregistrationV1alpha1WebhookClientConfig }, property_name::Symbol) =  _property_map_IoK8sApiAuditregistrationV1alpha1WebhookClientConfig[property_name]

function check_required(o::IoK8sApiAuditregistrationV1alpha1WebhookClientConfig)
    true
end

function validate_property(::Type{ IoK8sApiAuditregistrationV1alpha1WebhookClientConfig }, name::Symbol, val)
    if name === Symbol("caBundle")
    end
end


text:
export FixedVariablesSet

struct FixedVariablesIdeal{V<:AbstractVariable, T<:Number, MT<:AbstractMonomialLike} <: AbstractPolynomialIdeal
    substitutions::Union{Nothing, Dict{V, T}}
end
function Base.convert(::Type{FixedVariablesIdeal{V, T, MT}}, ideal::FixedVariablesIdeal{V, T, MT}) where {V<:AbstractVariable, T<:Number, MT<:AbstractMonomialLike}
    return ideal
end
function Base.convert(::Type{FixedVariablesIdeal{V, T, MT}}, ideal::FixedVariablesIdeal{V, S, MT}) where {V, S, T, MT}
    subs = ideal.substitutions
    if subs !== nothing
        subs = convert(Dict{V, T}, subs)
    end
    return FixedVariablesIdeal{V, T, MT}(subs)
end

function MP.variables(ideal::FixedVariablesIdeal{V}) where V
    if generate_nonzero_constant(ideal)
        return V[]
    else
        return sort(collect(keys(ideal.substitutions)), rev=true)
    end
end
# In that case, the ideal can generate any polynomial.
function generate_nonzero_constant(I::FixedVariablesIdeal)
    return I.substitutions === nothing
end
function Base.rem(p::AbstractPolynomialLike, I::FixedVariablesIdeal)
    if generate_nonzero_constant(I)
        return zero(p)
    else
        substitutions = [key => value for (key, value) in I.substitutions]
        return subs(p, substitutions...)
    end
end

struct FixedVariablesSet{V, T, MT} <: AbstractAlgebraicSet
    ideal::FixedVariablesIdeal{V, T, MT}
end
MP.changecoefficienttype(::Type{FixedVariablesSet{V, S, MT}}, T::Type) where {V, S, MT} = FixedVariablesSet{V, T, MT}
function Base.convert(::Type{FixedVariablesSet{V, T, MT}}, set::FixedVariablesSet{V, T, MT}) where {V, T, MT}
    return set
end
function Base.convert(::Type{FixedVariablesSet{V, T, MT}}, set::FixedVariablesSet{V, S, MT}) where {V, S, T, MT}
    return FixedVariablesSet(convert(FixedVariablesIdeal{V, T, MT}, set.ideal))
end

ideal(set::FixedVariablesSet, args...) = set.ideal
MP.variables(set::FixedVariablesSet) = MP.variables(set.ideal)
function nequalities(set::FixedVariablesSet)
    if set.ideal.substitutions === nothing
        return 1
    else
        return length(set.ideal.substitutions)
    end
end
function equalities(set::FixedVariablesSet{V, T, MT}) where {V, T, MT}
    if set.ideal.substitutions === nothing
        return [constantterm(one(T), MT)]
    else
        return [key - value for (key, value) in set.ideal.substitutions]
    end
end
function Base.intersect(V::FixedVariablesSet{V1, T1, MT1}, W::FixedVariablesSet{V2, T2, MT2}) where {V1, V2, T1, T2, MT1, MT2}
    # For `DynamicPolynomials`, they have the same type and for
    # `TypedPolynomials`, promoting would give `Monomial`.
    VT = V1 == V2 ? V1 : AbstractVariable
    T = promote_type(T1, T2)
    if V.ideal.substitutions === nothing || W.ideal.substitutions === nothing
        sub = nothing
    else
        has_dup = false
        function combine(a, b)
            if a != b
                has_dup = true
            end
            return a
        end
        sub = Dict{VT, T}()
        merge!(combine, sub, V.ideal.substitutions)
        merge!(combine, sub, W.ideal.substitutions)
        if has_dup
            sub = nothing
        end
    end
    ideal = FixedVariablesIdeal{VT, T, promote_type(MT1, MT2)}(sub)
    return FixedVariablesSet(ideal)
end
function Base.intersect(V::AlgebraicSet, W::FixedVariablesSet)
    return V ∩ algebraicset(W)
end
function Base.intersect(V::FixedVariablesSet, W::AlgebraicSet)
    # Need `W` to be first so that the intersection uses its solver.
    # and not the default one taken in `algebraicset(W)`.
    return W ∩ V
end

iszerodimensional(::FixedVariablesSet) = true
function Base.isempty(V::FixedVariablesSet)
    return generate_nonzero_constant(V.ideal)
end

# Assumes `isempty(V)`
function only_point(V::FixedVariablesSet)
    subs = collect(V.ideal.substitutions)
    sort!(subs, rev = true, by = sub -> sub[1])
    return [sub[2] for sub in subs]
end

function Base.iterate(V::FixedVariablesSet, state=nothing)
    if state === nothing && !isempty(V)
        return only_point(V), true
    else
        return nothing
    end
end
Base.length(V::FixedVariablesSet) = isempty(V) ? 0 : 1

struct FixedVariable{V<:AbstractVariable, T}
    variable::V
    value::T
end
function equality(variable::AbstractVariable, value::Number)
    return FixedVariable(variable, value)
end
function equality(value::Number, variable::AbstractVariable)
    return FixedVariable(variable, value)
end

function Base.intersect(el::FixedVariable; kws...)
    subs = Dict(el.variable => el.value)
    return FixedVariablesSet(FixedVariablesIdeal{
        typeof(el.variable), typeof(el.value), typeof(el.variable)}(subs))
end


text:
using SAT
using Test

cnf = CNF([[1, -5, 4], [-1, 5, 3, 4], [-3, -4]])
@test cnf.M == 3
@test cnf.N == 5
σ = solve(cnf);
@test energy(cnf, σ) == 0

cnf = randomcnf(N=2000, α=9.4, k=4, seed=19)
σ = solve(cnf, rstep=0.001, maxiters=1000);
@test energy(cnf, σ) == 0

N = 1000; n=2
planted = [rand([-1,1], N) for _=1:n]
cnf = randomcnf(N=N, α=8., k=4, seed=18, planted=planted)
for p in planted
    @test energy(cnf, p) == 0
end

# E, σ = solve(N=10000, α=9.6, k=4, seed_cnf=19, rstep=0.0002, maxiters=1000);
# @test E == 0

# DECIMATION NOT WORKING
# E, σ = KSAT.solve(method=:decimation, N=10000,α=9.6, k=4, seed_cnf=19,
#         r=0.02, maxiters=1000);
# @test E == 0


text:
using Benchmark

################################################################################
# Test how matrix operations performance compared with elementwise operations
# when @inbounds annotation is turned on.
#
# | Row | Function         | Average     | Relative | Replications |
# |-----|------------------|-------------|----------|--------------|
# | 1   | "softmax_matrix" | 0.00168359  | 2.63545  | 50           |
# | 2   | "softmax_elem"   | 0.000638823 | 1.0      | 50           |
################################################################################

function softmax_matrix(input, the_output)
  output = copy(input)
  output .-= maximum(output, 3) # subtract max along channel dimension
  output = exp(output)
  output ./= sum(output, 3) # normalize along channel dimension
  copy!(the_output, output)
end

function softmax_elem(input, output)
  width, height, channels, num = size(input)

  for w = 1:width
    for h = 1:height
      for n = 1:num
        maxval = -Inf
        @simd for c = 1:channels
          @inbounds maxval = max(maxval, input[w,h,c,n])
        end
        @simd for c = 1:channels
          @inbounds output[w,h,c,n] = exp(input[w,h,c,n]-maxval)
        end
        the_sum = 0.0
        @simd for c = 1:channels
          @inbounds the_sum += output[w,h,c,n]
        end
        @simd for c = 1:channels
          @inbounds output[w,h,c,n] /= the_sum
        end
      end
    end
  end
end

input = rand(5, 5, 10, 128)
o1 = rand(size(input))
o2 = rand(size(input))
softmax_matrix() = softmax_matrix(input, o1)
softmax_elem() = softmax_elem(input, o2)

softmax_matrix()
softmax_elem()
@assert all(abs(o1-o2) .< 1e-10)

df = compare([softmax_matrix, softmax_elem], 50)
println("$df")


text:
using Graph500
using Test

@testset begin
  include("kronecker_generator.jl")
  include("kernel_1.jl")
  include("kernel_2.jl")
  include("validation.jl")
  include("output.jl")
  include("driver.jl")
end


text:
module CowinAPI
using HTTP, JSON
using DataFrames, CSV
using Dates
include("urls.jl")

export get_state, get_state_ids, get_districts, get_district_id
export find_by_pin, find_by_district, calendar_by_pin, calendar_by_district
export calendar_by_district_brief

export stateID

stateID = Dict("Punjab" => 28,"Telangana" => 32,"Bihar" => 5,"Chandigarh" => 6,"Chhattisgarh" => 7,"Jammu and Kashmir" => 14,"Lakshadweep" => 19,"Nagaland" => 25,"Tripura" => 33,"Gujarat" => 11,"Madhya Pradesh" => 20,"Goa" => 10,"Ladakh" => 18,"Odisha" => 26,"Puducherry" => 27,"Tamil Nadu" => 31,"Rajasthan" => 29,"Mizoram" => 24,"Karnataka" => 16,"Jharkhand" => 15,"Himachal Pradesh" => 13,"Assam" => 4,"Haryana" => 12,"Meghalaya" => 23,"Andhra Pradesh" => 2,"Sikkim" => 30,"Delhi" => 9,"West Bengal" => 36,"Daman and Diu" => 37,"Uttarakhand" => 35,"Dadra and Nagar Haveli" => 8,"Manipur" => 22,"Uttar Pradesh" => 34,"Arunachal Pradesh" => 3,"Maharashtra" => 21,"Kerala" => 17,"Andaman and Nicobar Islands" => 1)

function make_API_call(url::String)
    try
        response = HTTP.get(url)
        return String(response.body)
    catch e
        #return "Error occurred: $(e)"
        return e
    end
end

function get_states()
    url = get_states_url()
    try
        str = make_API_call(url)
        jobj = JSON.Parser.parse(str)
        datadict = jobj["states"]
        df = DataFrame(state_id=Int[], state_name=String[])
        nrows = size(datadict)[1]
        for i in 1:nrows
            push!(df.state_id, datadict[i]["state_id"])
            push!(df.state_name, datadict[i]["state_name"])
        end
        return df
    catch e
        return e
    end
end

function get_state_ids()
    states = []
    try
        states = get_states()
        stateIDs = Dict{String,Int}()
        for istate in 1:size(states)[1]
            stateIDs[states[istate, :state_name]]=states[istate, :state_id]
        end
        return stateIDs
    catch e
        return e
    end
end

#stateID = get_state_ids()

function get_districts(state::String)
    state_id = stateID[state]
    url = get_districts_url(state_id)
    try
        str = make_API_call(url)
        jobj = JSON.Parser.parse(str)
        datadict = jobj["districts"]
        df = DataFrame(district_id=Int[], district_name=String[])
        nrows = size(datadict)[1]
        for i in 1:nrows
            push!(df.district_id, datadict[i]["district_id"])
            push!(df.district_name, datadict[i]["district_name"])
        end
        return df
    catch e
        return e
    end
end

function get_district_id(state::String, district::String)
    try
        df = get_districts(state)

        return df[df[!,:district_name].==district,:district_id][1]
    catch e
        return e
    end
end

function find_by_pin(pin::Int, date::String)
    url = find_by_pin_url(pin, date)
    try
        str = make_API_call(url)
        jobj = JSON.Parser.parse(str)
        datadict = jobj["sessions"]
        nrows = size(datadict)[1]
        df = DataFrame(block_name=String[], lat=Float64[], long=Float64[], center_id=Int[],
                       state_name=String[], min_age_limit=Int[], pin=Int[], address=String[],name=String[],
                       session_id=String[], fee_type=String[], date=String[], available_capacity=Int[], vaccine=String[],
                       fee=String[], district_name=String[], from=String[], to=String[], slots=Array{String,1}[])
        for i in 1:nrows
            push!(df.block_name, datadict[i]["block_name"])
            push!(df.lat, datadict[i]["lat"])
            push!(df.long, datadict[i]["long"])
            push!(df.center_id, datadict[i]["center_id"])
            push!(df.state_name, datadict[i]["state_name"])
            push!(df.min_age_limit, datadict[i]["min_age_limit"])
            push!(df.pin, datadict[i]["pincode"])
            push!(df.address, datadict[i]["address"])
            push!(df.name, datadict[i]["name"])
            push!(df.session_id, datadict[i]["session_id"])
            push!(df.fee_type, datadict[i]["fee_type"])
            push!(df.date, datadict[i]["date"])
            push!(df.available_capacity, datadict[i]["available_capacity"])
            push!(df.vaccine, datadict[i]["vaccine"])
            push!(df.fee, datadict[i]["fee"])
            push!(df.district_name, datadict[i]["district_name"])
            push!(df.from, datadict[i]["from"])
            push!(df.to, datadict[i]["to"])
            push!(df.slots, Array{String,1}(datadict[i]["slots"]))
        end

        return df
    catch e
        return e
    end
end

function find_by_district(district_id::Int, date::String)
    url = find_by_district_url(district_id, date)
    try
        str = make_API_call(url)
        jobj = JSON.Parser.parse(str)
        datadict = jobj["sessions"]
        nrows = size(datadict)[1]
        df = DataFrame(block_name=String[], lat=Float64[], long=Float64[], center_id=Int[],
                       state_name=String[], min_age_limit=Int[], pin=Int[], address=String[],name=String[],
                       session_id=String[], fee_type=String[], date=String[], available_capacity=Int[], vaccine=String[],
                       fee=String[], district_name=String[], from=String[], to=String[], slots=Array{String,1}[])
        for i in 1:nrows
            push!(df.block_name, datadict[i]["block_name"])
            push!(df.lat, datadict[i]["lat"])
            push!(df.long, datadict[i]["long"])
            push!(df.center_id, datadict[i]["center_id"])
            push!(df.state_name, datadict[i]["state_name"])
            push!(df.min_age_limit, datadict[i]["min_age_limit"])
            push!(df.pin, datadict[i]["pincode"])
            push!(df.address, datadict[i]["address"])
            push!(df.name, datadict[i]["name"])
            push!(df.session_id, datadict[i]["session_id"])
            push!(df.fee_type, datadict[i]["fee_type"])
            push!(df.date, datadict[i]["date"])
            push!(df.available_capacity, datadict[i]["available_capacity"])
            push!(df.vaccine, datadict[i]["vaccine"])
            push!(df.fee, datadict[i]["fee"])
            push!(df.district_name, datadict[i]["district_name"])
            push!(df.from, datadict[i]["from"])
            push!(df.to, datadict[i]["to"])
            push!(df.slots, Array{String,1}(datadict[i]["slots"]))
        end
        return df
    catch e
        return e
    end
end

function calendar_by_pin(pin::Int, date::String)
    url = calendar_by_pin_url(pin, date)
    try
        str = make_API_call(url)
        jobj = JSON.Parser.parse(str)
        datadict = jobj["centers"]
        df = DataFrame(block_name=String[], date=String[], available_capacity=Float64[], vaccine=String[],
                       lat=Float64[], long=Float64[], center_id=Int[],
                       state_name=String[], pin=Int[], address=String[], name=String[], session_id=String[],
                       slots=Array{String,1}[], min_age_limit=Int[],
                       fee_type=String[], district_name=String[], to=String[], from=String[])

        nrows = size(datadict)[1]
        for i in 1:nrows
            nsessions = size(datadict[i]["sessions"])[1]
            for isession in 1:nsessions
                sessiondata = datadict[i]["sessions"][isession]
                push!(df.block_name, datadict[i]["block_name"])
                push!(df.lat, datadict[i]["lat"])
                push!(df.long, datadict[i]["long"])
                push!(df.center_id, datadict[i]["center_id"])
                push!(df.state_name, datadict[i]["state_name"])
                push!(df.pin, datadict[i]["pincode"])
                push!(df.address, datadict[i]["address"])
                push!(df.name, datadict[i]["name"])

                # sessiondata
                push!(df.session_id, sessiondata["session_id"])
                push!(df.slots, Array{String,1}(sessiondata["slots"]))
                push!(df.date, sessiondata["date"])
                push!(df.available_capacity, sessiondata["available_capacity"])
                push!(df.min_age_limit, sessiondata["min_age_limit"])
                push!(df.vaccine, sessiondata["vaccine"])
                # sessiondata END

                push!(df.fee_type, datadict[i]["fee_type"])
                push!(df.district_name, datadict[i]["district_name"])
                push!(df.from, datadict[i]["from"])
                push!(df.to, datadict[i]["to"])
            end
        end
        return df
    catch e
        return e
    end
end

function calendar_by_district(district_id::Int, date::String)
    url = calendar_by_district_url(district_id, date)
    try
        str = make_API_call(url)
        jobj = JSON.Parser.parse(str)
        datadict = jobj["centers"]

        df = DataFrame(block_name=String[],date=String[], available_capacity=Float64[], available_capacity_dose1=Float64[],
                       available_capacity_dose2=Float64[], vaccine=String[],
                       lat=Float64[], long=Float64[], center_id=Int[],
                       state_name=String[], pin=Int[], address=String[], name=String[], session_id=String[],
                       slots=Array{String,1}[], min_age_limit=Int[],
                       fee_type=String[], district_name=String[], to=String[], from=String[])
                       nrows = size(datadict)[1]
        for i in 1:nrows
            nsessions = size(datadict[i]["sessions"])[1]
            for isession in 1:nsessions
                sessiondata = datadict[i]["sessions"][isession]
                push!(df.block_name, datadict[i]["block_name"])
                push!(df.lat, datadict[i]["lat"])
                push!(df.long, datadict[i]["long"])
                push!(df.center_id, datadict[i]["center_id"])
                push!(df.state_name, datadict[i]["state_name"])
                push!(df.pin, datadict[i]["pincode"])
                push!(df.address, datadict[i]["address"])
                push!(df.name, datadict[i]["name"])

                # sessiondata
                push!(df.session_id, sessiondata["session_id"])
                push!(df.slots, Array{String,1}(sessiondata["slots"]))
                push!(df.date, sessiondata["date"])
                push!(df.available_capacity, sessiondata["available_capacity"])
                push!(df.available_capacity_dose1, sessiondata["available_capacity_dose1"])
                push!(df.available_capacity_dose2, sessiondata["available_capacity_dose2"])
                push!(df.min_age_limit, sessiondata["min_age_limit"])
                push!(df.vaccine, sessiondata["vaccine"])
                # sessiondata END

                push!(df.fee_type, datadict[i]["fee_type"])
                push!(df.district_name, datadict[i]["district_name"])
                push!(df.from, datadict[i]["from"])
                push!(df.to, datadict[i]["to"])
            end
        end
        return df
    catch e
        return e
    end
end

function calendar_by_district_brief(district_id::Int, date::String)
    url = calendar_by_district_url(district_id, date)
    try
        str = make_API_call(url)
        jobj = JSON.Parser.parse(str)
        datadict = jobj["centers"]

        df = DataFrame(center_id=Int[],
                       state_name=String[], pin=Int[], name=String[],
                       date=Dates.Date[], available_capacity=Float64[], available_capacity_dose1=Float64[],
                       available_capacity_dose2=Float64[], min_age_limit=Int[],
                       vaccine=String[], district_name=String[])
                       nrows = size(datadict)[1]
        for i in 1:nrows
            nsessions = size(datadict[i]["sessions"])[1]
            for isession in 1:nsessions
                sessiondata = datadict[i]["sessions"][isession]
                push!(df.center_id, datadict[i]["center_id"])
                push!(df.state_name, datadict[i]["state_name"])
                push!(df.pin, datadict[i]["pincode"])
                push!(df.name, datadict[i]["name"])

                # sessiondata
                push!(df.date, Date(sessiondata["date"],"dd-mm-yyyy"))
                push!(df.available_capacity, sessiondata["available_capacity"])
                push!(df.available_capacity_dose1, sessiondata["available_capacity_dose1"])
                push!(df.available_capacity_dose2, sessiondata["available_capacity_dose2"])
                push!(df.min_age_limit, sessiondata["min_age_limit"])
                push!(df.vaccine, sessiondata["vaccine"])
                # sessiondata END

                push!(df.district_name, datadict[i]["district_name"])
            end
        end
        return df
    catch e
        return e
    end
end
# function download_certificate(beneficiary_reference_id::String)
#     url = download_certificate_url(beneficiary_reference_id)
#     return url
# end

end # module