text:
module MockFBA

#using LoopVectorization
using FITSIO
using YAML
using CSV
using JLD2
using Formatting
using Printf
using DelimitedFiles
using DataInterpolations
using StaticArrays
using LoopVectorization
using DataFrames
using CodecLz4
using Random
using Statistics
using Dates

export preprocess_tracer

include("FBA_utility.jl")
include("Hardware.jl")
include("Tiles.jl")
include("Geometry.jl")
include("Positioner.jl")
include("ZonesTargets.jl")
include("PostProcess.jl")

#Should be load for development but not deployment
include("PlotLibrary.jl")

"""This goes through the full proe-process step for a particular tracer
 Steps:
  1) call the pre_process_FITS2JLD_tracer which splits the tracer in NGC/SGC and 20x20 deg zones
  2) Assign TILEID per_pass for the objects by load the tile file
  3) Generate indexing fro RA,DEC
  4) Generate the listing for TILEID across all passes

   # Arguments
   - `fname::String`: Name of the fits file containing tracer objects, Only needs RA,DEC columns
   - `priority_list::Array[Integer]`: The list of priority to be assigned randomly
   - `priority_frac_list::Array[Float]`: Fraction of objects in each priority
   - `outfile::String`: Filename for JLD2 file to be created for efficient I/O
   - `tile_file::String`: fits file containing list tiles

"""
function preprocess_tracer(fname,priority_list,priority_frac_list,outfile,
        tile_file,tile_radius_indeg,npass,program)

    #generates the JLD2 file
    pre_process_FITS2JLD_tracer(fname,priority_list,priority_frac_list,outfile)

    #Assign TILEID to the object and ad to JLD2 file for each pass
    overwrite=false #To not overwrite the tileid field if it alredy exists
    #For SGC first and then NGC
    for group in ["NGC","SGC"]
        Append_TileID(outfile,tile_file,group,program,npass,tile_radius_indeg,overwrite)
        MockFBA.Append_column(outfile,group,["XYZ_UNIT"]) #Adding XYZ_UNIT
    end

    #Add the indexing to the file for RA and DEC and listing for tileid
    tileid_cols=[]
    for pass in 0:npass-1
        append!(tileid_cols,["TILEID_PASS$(pass)"])
    end

    for group in ["NGC","SGC",""]
        Indexing_column(outfile,group,["RA","DEC"],nothing)
        udic=LISTING_column(outfile,group,"TILEID",tileid_cols)
    end
end

end #end of module

max_stars_repo_name: shadaba/MockFBA

max_issues_repo_name: shadaba/MockFBA

max_forks_repo_name: shadaba/MockFBA


text:
function testme(MAXTRIES=3; tol=0.1)

    # Verify symmetric kernel calculation by comparing with numerical integration results
    for i = 1:MAXTRIES
        @assert(test_GKGT_random(tol=tol))
    end

    # Verify integrals

    for i = 1:MAXTRIES

        lowery = rand() * 60 - 30
        uppery = lowery + rand() * 60

        lowerx = rand() * 60 - 30
        upperx = lowerx + rand() * 60

        ℓ = rand() * (5.0 - 0.01) + 0.01

        a = lowerx
        b = upperx
        c = lowery
        d = uppery

        #####
        # 1 #
        #####
        @printf("verify_1\n")
        local out = verify_1(lowery = lowery, uppery = uppery, lowerx = lowerx, upperx = upperx, ℓ=ℓ)
        display(out)
        @assert(abs(out[1]-out[3]) < tol)

        @printf("\t verify_sub_1\n\n")
        out = verify_sub_1(b=b, c=c, d=d, ℓ=ℓ)
        display(out)
        @assert(abs(out[1]-out[2]) < tol)
        @assert(abs(out[1]-out[3]) < tol)

        #####
        # 2 #
        #####
        @printf("verify_2\n")
        out = verify_2(lowery = lowery, uppery = uppery, lowerx = lowerx, upperx = upperx, ℓ=ℓ)
        display(out)
        @assert(abs(out[1]-out[3]) < tol)

        @printf("\t verify_sub_2a\n")
        out = verify_sub_2a(a=a,c=c,d=d, ℓ=ℓ)
        display(out)
        @assert(abs(out[1]-out[2]) < tol)

        @printf("\t verify_sub_2b\n\n")
        out = verify_sub_2b(a=a,c=c,d=d, ℓ=ℓ)
        display(out)
        @assert(abs(out[1]-out[2]) < tol)

        #####
        # 3 #
        #####
        @printf("verify_3\n")
        out = verify_3(lowery = lowery, uppery = uppery, lowerx = lowerx, upperx = upperx, ℓ=ℓ)
        display(out)
        @assert(abs(out[1]-out[3]) < tol)

        @printf("\t verify_sub_3a\n")
        out = verify_sub_3a(a=a,c=c,d=d, ℓ=ℓ)
        display(out)
        @assert(abs(out[1]-out[2]) < tol)

        @printf("\t verify_sub_3b\n\n")
        out = verify_sub_3b(a=a,c=c,d=d, ℓ=ℓ)
        display(out)
        @assert(abs(out[1]-out[2]) < tol)

    end

    @info("if you see this message, then tests must be succesful")

end

ext: jl


text:
if ENV["HOME"] == "/mnt/juliabox"
    Pkg.dir(path...)=joinpath("/home/jrun/.julia/v0.6",path...)
else
    for p in ("Knet","JLD","JSON","Images") # ,"WordTokenizers")
        Pkg.installed(p) == nothing && Pkg.add(p)
    end
end
using Images,JLD,Knet # ,WordTokenizers
global atype = gpu()<0 ? Array{Float32}:KnetArray{Float32}

server="people.csail.mit.edu/deniz/"
if !isdir("data/demo")
    info("Downloading sample questions and images from CLEVR dataset...")
    download(server*"data/mac-network/demo.tar.gz","demo.tar.gz")
    run(`tar -xzf demo.tar.gz`)
    rm("demo.tar.gz")
end
if !isfile("models/macnet.jld")
    info("Downloading pre-trained model from our servers...")
    download(server*"models/mac-network/demo_model.jld","models/macnet.jld")
end

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
"""
2520 is the smallest number that can be divided by each of the numbers from 1 to 10 without any remainder.

What is the smallest positive number that is evenly divisible by all of the numbers from 1 to 20?
"""

ismultiple(x::Int,r::UnitRange) = all([x%i == 0 for i in r])

function smallestmutiple(nmin,nmax)
    r = collect(nmin:nmax)
    multiples = []

    for i in eachindex(r)   # store only the common divisor
        t = r[i]
        t == 1 && continue
        push!(multiples,r[i])

        for j in i:length(r)
            if r[j] % t == 0
                r[j] = r[j] / t
            end
        end
    end
    return prod(multiples)
end

@assert smallestmutiple(1,10) == 2520
@show smallestmutiple(1,20)


text:
using PastaQ
using Random
using ITensors
using Test
using LinearAlgebra

#
# Helper functions for tests
#

function state_to_int(state::Array)
  index = 0
  for j in 1:length(state)
    index += 2^(j - 1) * state[length(state) + 1 - j]
  end
  return index
end

function empiricalprobability(samples::Matrix)
  prob = zeros((1 << size(samples)[2]))
  for n in 1:size(samples)[1]
    sample = samples[n, :]
    index = state_to_int(sample)
    prob[index + 1] += 1
  end
  prob = prob / size(samples)[1]
  return prob
end

@testset "generation of preparation states" begin
  N = 4
  nstates = 100
  states = PastaQ.randompreparations(N, nstates)
  @test size(states)[1] == nstates
  @test size(states)[2] == N
end

@testset "generation of measurement bases" begin
  N = 4
  nbases = 100
  bases = randombases(N, nbases)
  @test size(bases)[1] == nbases
  @test size(bases)[2] == N

  bases = fullbases(N)
  @test bases isa Matrix{String}
  @test size(bases) == (3^N,N)
end

@testset "measurements" begin
  Random.seed!(1234)
  N = 3
  depth = 4
  ψ0 = productstate(N)
  gates = randomcircuit(N; depth =  depth)
  ψ = runcircuit(ψ0, gates)
  ψ_vec = PastaQ.array(ψ)
  probs = abs2.(ψ_vec)

  nshots = 50000
  samples = PastaQ.getsamples(ψ, nshots)
  @test size(samples)[1] == nshots
  @test size(samples)[2] == N
  data_prob = empiricalprobability(samples)
  @test probs ≈ data_prob atol = 1e-2

  ρ = runcircuit(N, gates; noise=("amplitude_damping", (γ=0.01,)))
  ρ_mat = PastaQ.array(ρ)
  probs = real(diag(ρ_mat))

  samples = PastaQ.getsamples(ρ, nshots)
  @test size(samples)[1] == nshots
  @test size(samples)[2] == N
  data_prob = empiricalprobability(samples)
  @test probs ≈ data_prob atol = 1e-2
end

@testset "basis rotations" begin
  s = siteinds("Qubit",1)
  #ϕ = qubits(1)
  ψ0 = state("X+",s[1])
  gates = [("basisX",1,(dag=true,))]
  ψ = runcircuit(ψ0,gates)
  @test PastaQ.array(ψ) ≈ state("Z+")
  ψ0 = state("X-",s[1])
  gates = [("basisX",1,(dag=true,))]
  ψ = runcircuit(ψ0,gates)
  @test PastaQ.array(ψ) ≈ state("Z-")

  ψ0 = state("Y+",s[1])
  gates = [("basisY",1,(dag=true,))]
  ψ = runcircuit(ψ0,gates)
  @test PastaQ.array(ψ) ≈ state("Z+")
  ψ0 = state("Y-",s[1])
  gates = [("basisY",1,(dag=true,))]
  ψ = runcircuit(ψ0,gates)
  @test PastaQ.array(ψ) ≈ state("Z-")
end

@testset "project unitary MPO" begin
  N = 4
  ntrial = 100
  gates = randomcircuit(N; depth =  4, layered=false)

  U = runcircuit(N, gates; process=true)

  bases = randombases(N, ntrial)
  preps = randompreparations(N, ntrial)

  for n in 1:ntrial
    mgates = PastaQ.measurementgates(bases[n, :])
    ψ_in = productstate(N, preps[n, :])
    ψ_out = runcircuit(ψ_in, gates)

    Ψ_out = PastaQ.projectchannel(U, preps[n, :])
    @test PastaQ.array(ψ_out) ≈ PastaQ.array(Ψ_out)

    ψ_m = runcircuit(ψ_out, mgates)
    Ψ_m = runcircuit(Ψ_out, mgates)
    @test PastaQ.array(ψ_m) ≈ PastaQ.array(Ψ_m)
  end
end

@testset "project unitary ITensor" begin
  N = 4
  ntrial = 100
  gates = randomcircuit(N; depth =  4, layered=false)

  U = runcircuit(N, gates; process=true, full_representation = true)

  bases = randombases(N, ntrial)
  preps = randompreparations(N, ntrial)

  for n in 1:ntrial
    mgates = PastaQ.measurementgates(bases[n, :])
    ψ_in = productstate(N, preps[n, :])
    ψ_out = runcircuit(ψ_in, gates)

    Ψ_out = PastaQ.projectchannel(U, preps[n, :])
    @test PastaQ.array(ψ_out) ≈ PastaQ.array(Ψ_out)

    ψ_m = runcircuit(ψ_out, mgates)
    Ψ_m = runcircuit(Ψ_out, mgates)
    @test PastaQ.array(ψ_m) ≈ PastaQ.array(Ψ_m)
  end
end

@testset "project Choi MPO" begin
  N = 4
  ntrial = 100
  gates = randomcircuit(N; depth =  4, layered=false)

  Λ = runcircuit(N, gates; process=true, noise=("amplitude_damping", (γ=0.1,)))

  bases = randombases(N, ntrial)
  preps = PastaQ.randompreparations(N, ntrial)
  for n in 1:ntrial
    mgates = PastaQ.measurementgates(bases[n, :])
    ψ_in = productstate(N, preps[n, :])
    ρ_out = runcircuit(ψ_in, gates; noise=("amplitude_damping", (γ=0.1,)))

    Λ_out = PastaQ.projectchannel(Λ, preps[n, :])
    @test PastaQ.array(ρ_out) ≈ PastaQ.array(Λ_out) atol = 1e-6

    ρ_m = runcircuit(ρ_out, mgates)
    Λ_m = runcircuit(Λ_out, mgates)
    @test PastaQ.array(ρ_m) ≈ PastaQ.array(Λ_m) atol = 1e-6
  end
end

@testset "project Choi ITensor" begin
  N = 4
  ntrial = 100
  gates = randomcircuit(N; depth =  4, layered=false)

  @disable_warn_order begin
    Λ = runcircuit(N, gates; process=true, noise=("amplitude_damping", (γ=0.1,)), full_representation = true)

    bases = randombases(N, ntrial)
    preps = PastaQ.randompreparations(N, ntrial)
    for n in 1:ntrial
      mgates = PastaQ.measurementgates(bases[n, :])
      ψ_in = productstate(N, preps[n, :])
      ρ_out = runcircuit(ψ_in, gates; noise=("amplitude_damping", (γ=0.1,)))

      Λ_out = PastaQ.projectchannel(Λ, preps[n, :])
      @test PastaQ.array(ρ_out) ≈ PastaQ.array(Λ_out) atol = 1e-6

      ρ_m = runcircuit(ρ_out, mgates)
      Λ_m = runcircuit(Λ_out, mgates)
      @test PastaQ.array(ρ_m) ≈ PastaQ.array(Λ_m) atol = 1e-6
    end
  end
end

@testset "getsamples: states" begin
  N = 3
  circuit = randomcircuit(N; depth = 3)

  # quantum states
  ψ = runcircuit(N, circuit)
  ρ = runcircuit(N, circuit; noise=("amplitude_damping", (γ=0.1,)))

  nbases = 11
  bases = randombases(N, nbases)
  data = getsamples(ψ, bases)
  @test size(data) == (nbases,N)
  nshots = 3
  data = getsamples(ψ, bases, nshots)
  @test size(data) == (nbases*nshots,N)
  for b in 1:nbases
    for k in 1:nshots
      @test first.(data[(b-1)*nshots+k,:]) == bases[b,:]
    end
  end

  data = getsamples(ρ, bases)
  @test size(data) == (nbases,N)
  nshots = 3
  data = getsamples(ρ, bases, nshots)
  @test size(data) == (nbases*nshots,N)
  for b in 1:nbases
    for k in 1:nshots
      @test first.(data[(b-1)*nshots+k,:]) == bases[b,:]
    end
  end

  # ITensors quantum states
  ψ = runcircuit(N, circuit; full_representation = true)
  ρ = runcircuit(N, circuit; noise=("amplitude_damping", (γ=0.1,)), full_representation = true)

  nbases = 11
  bases = randombases(N, nbases)
  data = getsamples(ψ, bases)
  @test size(data) == (nbases,N)
  nshots = 3
  data = getsamples(ψ, bases, nshots)
  @test size(data) == (nbases*nshots,N)
  for b in 1:nbases
    for k in 1:nshots
      @test first.(data[(b-1)*nshots+k,:]) == bases[b,:]
    end
  end

  data = getsamples(ρ, bases)
  @test size(data) == (nbases,N)
  nshots = 3
  data = getsamples(ρ, bases, nshots)
  @test size(data) == (nbases*nshots,N)
  for b in 1:nbases
    for k in 1:nshots
      @test first.(data[(b-1)*nshots+k,:]) == bases[b,:]
    end
  end
end

@testset "getsamples: channels" begin

  N = 3
  circuit = randomcircuit(N; depth = 3)

  # quantum processes
  U = runcircuit(N, circuit; process = true)
  Λ = runcircuit(N, circuit; noise=("amplitude_damping", (γ=0.1,)), process = true)

  npreps = 5
  nbases = 7
  preps = randompreparations(N, npreps)
  bases = randombases(N, nbases)

  data = getsamples(U, preps, bases)
  @test size(data) == (npreps*nbases,N)
  nshots = 3
  data = getsamples(U, preps, bases, nshots)
  @test size(data) == (npreps*nbases*nshots,N)

  for p in 1:npreps
    for b in 1:nbases
      for k in 1:nshots
        pdata = first.(data[(p-1)*nbases*nshots+(b-1)*nshots+k,:])
        bdata = first.(last.(data[(p-1)*nbases*nshots+(b-1)*nshots+k,:]))
        @test pdata == preps[p,:]
        @test bdata == bases[b,:]
      end
    end
  end

  npreps = 5
  nbases = 5
  preps = randompreparations(N, npreps)
  bases = randombases(N, nbases)

  data = getsamples(U, preps .=> bases)
  @test size(data) == (npreps,N)

  data = getsamples(Λ, preps, bases)
  @test size(data) == (npreps*nbases,N)
  nshots = 3
  data = getsamples(Λ, preps, bases, nshots)
  @test size(data) == (npreps*nbases*nshots,N)

  for p in 1:npreps
    for b in 1:nbases
      for k in 1:nshots
        pdata = first.(data[(p-1)*nbases*nshots+(b-1)*nshots+k,:])
        bdata = first.(last.(data[(p-1)*nbases*nshots+(b-1)*nshots+k,:]))
        @test pdata == preps[p,:]
        @test bdata == bases[b,:]
      end
    end
  end

  npreps = 5
  nbases = 5
  preps = randompreparations(N, npreps)
  bases = randombases(N, nbases)

  data = getsamples(Λ, preps .=> bases)
  @test size(data) == (npreps,N)

  # full representation
  U = runcircuit(N, circuit; process = true, full_representation = true)
  Λ = runcircuit(N, circuit; noise=("amplitude_damping", (γ=0.1,)), process = true,full_representation = true)

  npreps = 5
  nbases = 7
  preps = randompreparations(N, npreps)
  bases = randombases(N, nbases)

  data = getsamples(U, preps, bases)
  @test size(data) == (npreps*nbases,N)
  nshots = 3
  data = getsamples(U, preps, bases, nshots)
  @test size(data) == (npreps*nbases*nshots,N)

  for p in 1:npreps
    for b in 1:nbases
      for k in 1:nshots
        pdata = first.(data[(p-1)*nbases*nshots+(b-1)*nshots+k,:])
        bdata = first.(last.(data[(p-1)*nbases*nshots+(b-1)*nshots+k,:]))
        @test pdata == preps[p,:]
        @test bdata == bases[b,:]
      end
    end
  end

  data = getsamples(Λ, preps, bases)
  @test size(data) == (npreps*nbases,N)
  nshots = 3
  data = getsamples(Λ, preps, bases, nshots)
  @test size(data) == (npreps*nbases*nshots,N)

  for p in 1:npreps
    for b in 1:nbases
      for k in 1:nshots
        pdata = first.(data[(p-1)*nbases*nshots+(b-1)*nshots+k,:])
        bdata = first.(last.(data[(p-1)*nbases*nshots+(b-1)*nshots+k,:]))
        @test pdata == preps[p,:]
        @test bdata == bases[b,:]
      end
    end
  end
end


text:
# Testcase from example given in Mocking.jl's README
@testset "readme" begin
    # Note: Function only works in UNIX environments.
    function randdev(n::Integer)
        @mock open("/dev/urandom") do fp
            reverse(read(fp, n))
        end
    end

    n = 10
    if Sys.isunix()
        result = randdev(n)  # Reading /dev/urandom only works on UNIX environments
        @test eltype(result) == UInt8
        @test length(result) == n
    end

    # Produces a string with sequential UInt8 values from 1:n
    data = unsafe_string(pointer(convert(Array{UInt8}, 1:n)))

    # Generate a alternative method of `open` which call we wish to mock
    patch = @patch open(fn::Function, f::AbstractString) = fn(IOBuffer(data))

    # Apply the patch which will modify the behaviour for our test
    apply(patch) do
        @test randdev(n) == convert(Array{UInt8}, n:-1:1)
    end

    if Sys.isunix()
        # Outside of the scope of the patched environment `@mock` is essentially a no-op
        @test randdev(n) != convert(Array{UInt8}, n:-1:1)
    end
end


text:
"""
	violates_bapat_sunder(A,B, tol = ATOL)

Checks if matrices `A` and `B` violate the Bapat-Sunder conjecture, see
[Boson bunching is not maximized by indistinguishable particles](https://arxiv.org/abs/2203.01306)
"""
function violates_bapat_sunder(A,B, tol = ATOL)

	diagonal_el_product = prod([B[i,i] for i in 1:size(B)[1]])

	# @test (imag(diagonal_el_product) ≈ 0.)

	if !(imag(diagonal_el_product) ≈ 0.)
		throw(Exception("product of diagonal elements of B is not real"))
	elseif !(is_positive_semidefinite(A) && is_positive_semidefinite(B))
		throw(Exception("A or B not semidefinitepositive"))
	else
		return abs(ryser(A .* B)) / (abs(ryser(A)) * real(prod([B[i,i] for i in 1:size(B)[1]]))) > 1+tol
	end

	return nothing

end


text:
export RK1Step, RK4Step

struct RK1Step
end
struct RK4Step
end

# changes: added d as input and
function step(stepper::RK1Step,odefun,J,U::AbstractArray{R},Θ,tk::R,tkp1::R) where R <: Real
    U  += (tkp1-tk) .* odefun(J, U, Θ, tk)
    return U
end

function step(stepper::RK4Step,odefun,J,U0::AbstractArray{R},Θ,tk::R,tkp1::R) where R<:Real
    h    = tkp1 - tk

    K  = h .* odefun(J, U0, Θ, tk)
    U  = U0 + R.(1/6) .* K

    K  = h .* odefun(J, U0 + 0.5 .* K, Θ, tk+h/2)
    U  = U + R.(2/6) .* K

    K  = h .* odefun(J, U0 + 0.5 .* K, Θ, tk+h/2)
    U  = U + R.(2/6) .* K

    K  = h .* odefun(J, U0 + K, Θ, tk+h)
    U  = U + R.(1/6) .* K

    return U
end

"""
numerical integration of ODE
"""
function integrate(stepper,odefun,J,U::AbstractArray{R},Θ,tspan::AbstractArray{R},N::Int) where R <: Real
    h   = (tspan[2]-tspan[1])/N
    tk  = tspan[1]
    for k=1:N
        U = step(stepper,odefun,J, U, Θ,tk,tk+h)
        tk += h
    end
    return U
end

"""
numerical integration of ODE, with storing intermediate states
"""
function integrate2(stepper,odefun,J,U0::AbstractArray{R},Θ,tspan::AbstractArray{R},N::Int) where R <: Real
    h = (tspan[2]-tspan[1])/N
    tk = tspan[1]
    UArray = zeros(R,tuple([size(U0)... N+1]...))
    UArray[:,:,1] = U0
    d   = size(U0,1)-4
    X0  = (U0[1:d,:])
    for k=1:N
        temp = step(stepper,odefun, J, UArray[:,:,k], Θ , tk, tk+h)
        UArray[:,:,k+1] = temp
		tk += h
    end
    return UArray
end

max_stars_repo_name: EmoryMLIP/MFGnet

max_issues_repo_name: EmoryMLIP/MFGnet

max_forks_repo_name: EmoryMLIP/MFGnet


text:
# broadcasting

using Base.Broadcast: BroadcastStyle, Broadcasted

struct CuArrayStyle{N} <: AbstractGPUArrayStyle{N} end
CuArrayStyle(::Val{N}) where N = CuArrayStyle{N}()
CuArrayStyle{M}(::Val{N}) where {N,M} = CuArrayStyle{N}()

BroadcastStyle(::Type{CuArray{T,N}}) where {T,N} = CuArrayStyle{N}()

Base.similar(bc::Broadcasted{CuArrayStyle{N}}, ::Type{T}) where {N,T} =
    similar(CuArray{T}, axes(bc))

Base.similar(bc::Broadcasted{CuArrayStyle{N}}, ::Type{T}, dims) where {N,T} =
    CuArray{T}(undef, dims)

## replace base functions with libdevice alternatives

cufunc(f) = f
cufunc(::Type{T}) where T = (x...) -> T(x...) # broadcasting type ctors isn't GPU compatible

Broadcast.broadcasted(::CuArrayStyle{N}, f, args...) where {N} =
  Broadcasted{CuArrayStyle{N}}(cufunc(f), args, nothing)

const device_intrinsics = :[
  cos, cospi, sin, sinpi, tan, acos, asin, atan,
  cosh, sinh, tanh, acosh, asinh, atanh,
  log, log10, log1p, log2, logb, ilogb,
  exp, exp2, exp10, expm1, ldexp,
  erf, erfinv, erfc, erfcinv, erfcx,
  brev, clz, ffs, byte_perm, popc,
  isfinite, isinf, isnan, nearbyint,
  nextafter, signbit, copysign, abs,
  sqrt, rsqrt, cbrt, rcbrt, pow,
  ceil, floor, saturate,
  lgamma, tgamma,
  j0, j1, jn, y0, y1, yn,
  normcdf, normcdfinv, hypot,
  fma, sad, dim, mul24, mul64hi, hadd, rhadd, scalbn].args

for f in device_intrinsics
  isdefined(Base, f) || continue
  @eval cufunc(::typeof(Base.$f)) = $f
end

# broadcast ^

culiteral_pow(::typeof(^), x::T, ::Val{0}) where {T<:Real} = one(x)
culiteral_pow(::typeof(^), x::T, ::Val{1}) where {T<:Real} = x
culiteral_pow(::typeof(^), x::T, ::Val{2}) where {T<:Real} = x * x
culiteral_pow(::typeof(^), x::T, ::Val{3}) where {T<:Real} = x * x * x
culiteral_pow(::typeof(^), x::T, ::Val{p}) where {T<:Real,p} = pow(x, Int32(p))

cufunc(::typeof(Base.literal_pow)) = culiteral_pow
cufunc(::typeof(Base.:(^))) = pow

using MacroTools

const _cufuncs = [copy(device_intrinsics); :^]
cufuncs() = (global _cufuncs; _cufuncs)

_cuint(x::Int) = Int32(x)
_cuint(x::Expr) = x.head == :call && x.args[1] == :Int32 && x.args[2] isa Int ? Int32(x.args[2]) : x
_cuint(x) = x

function _cupowliteral(x::Expr)
  if x.head == :call && x.args[1] == :(CUDA.cufunc(^)) && x.args[3] isa Int32
    num = x.args[3]
    if 0 <= num <= 3
      sym = gensym(:x)
      new_x = Expr(:block, :($sym = $(x.args[2])))

      if iszero(num)
        push!(new_x.args, :(one($sym)))
      else
        unroll = Expr(:call, :*)
        for x = one(num):num
          push!(unroll.args, sym)
        end
        push!(new_x.args, unroll)
      end

      x = new_x
    end
  end
  x
end
_cupowliteral(x) = x

function replace_device(ex)
  global _cufuncs
  MacroTools.postwalk(ex) do x
    x = x in _cufuncs ? :(CUDA.cufunc($x)) : x
    x = _cuint(x)
    x = _cupowliteral(x)
    x
  end
end

macro cufunc(ex)
  global _cufuncs
  def = MacroTools.splitdef(ex)
  f = def[:name]
  def[:name] = Symbol(:cu, f)
  def[:body] = replace_device(def[:body])
  push!(_cufuncs, f)
  quote
    $(esc(MacroTools.combinedef(def)))
    CUDA.cufunc(::typeof($(esc(f)))) = $(esc(def[:name]))
  end
end


text:
using DrWatson
quickactivate(@__DIR__)
@show projectdir()
##
using CovidSurvey
using CSV
using DataFrames
using Dates
# using RCall
using Turing
using StatsPlots
using Plots.PlotMeasures
using Random
using BSON
using PrettyTables
plotlyjs()
scalefontsizes(1.5)
## ===========================================================================
# load data and model
# ============================================================================
fdir = projectdir("reports/regional/")
fname = "cases_start=2020-06-01_chains=10_epidemic_start=30_folder=regional_model=hospit_num_impute=6_observations_end=2021-01-13_prefix=_rw_step=1_steps=4000_warmup=2000"
PREFIX = ""
#-----------------------------------------------------------------------------
# read params
_, parsed_args, _ = parse_savename(fname, parsetypes = (Int, Float64, Bool))
!("preds" in keys(parsed_args)) && (parsed_args["preds"] = nothing)
# ps = NamedTuple{Tuple(Symbol.(keys(ps)))}(values(ps))

#------------------------------------------------------------------------------
# read data
data_params = (
      observations_end  = parsed_args["observations_end"]
    , predictors        = parsed_args["preds"]|> CovidSurvey.parse_predictors
    , cases_start       = parsed_args["start"]
    , rw_step           = parsed_args["prefix=_rw_step"]
    , epidemic_start    = parsed_args["epidemic_start"]
    , num_impute        = parsed_args["num_impute"]
    , link              = KLogistic(3.)
    , invlink           = KLogit(3.)
    , lockdown          = "2020-03-18"
    , covariates_kwargs = Dict(
      :fname => normpath( homedir(), "data/covidsurvey/smoothed_contacts.csv" ),
      :shift => -1,
      :startdate => "2020-11-10", # >= 1000 cases / day
    )
)

data        = Regional.load_data(; data_params... )
turing_data = data.turing_data;
model       = Regional.name2model[parsed_args["model"]]
#-----------------------------------------------------------------------------
# read chain
function not(c::Chains, i)
    n = size(c,3)
    s = filter(x->x ∉ i, 1:n)
    return c[:,:,s]
end
chainraw = read(fdir*PREFIX*"CHAIN_"*fname*".jls", Chains)
chainraw = chainraw[parsed_args["warmup"]+1:end,:,:]
chainraw = chainraw[1:10:end,:,:]
chain    = not(chainraw, [7,8])
# chain = chain[1:10:end,:,[1,2,3,5,6,7,9,10]]
# fname = normpath( fdir, savename(PREFIX*"CHAIN", ps, "jls") )
# safesave( fname, chain)

#-----------------------------------------------------------------------------
# plot chain
n = filter( x->!occursin(r"latent_Rts", x), String.(names(chain)))
let
    p = plot(chain[n]);
    saveto = normpath( fdir, PREFIX*"CHAINSPLOT_"*fname*".html")
    savefig(p, saveto )
    run(`firefox $(saveto)`, wait=false)
end

## ===========================================================================
# plot results
# ============================================================================
@info "make predictions"
generated_posterior = Regional.posterior(model, turing_data, chain)
let
    dic = Dict( zip( keys(generated_posterior), values(generated_posterior) ) )
    bson( normpath( fdir, "GENERATED-QUANTITIES_"*fname*".bson" ) ,  dic )
end
## ==========================================================================
# pgfplotsx()
for r in Regional.regions
    recipe = Regional.RegionPlottingRecipe(data, generated_posterior, r)
    p = plot(recipe)
    savefig( p, normpath( fdir, "FIG-$(uppercase(r))_"*fname*".html" ))
end
##
let
    r = Regional.RtsPlottingRecipe(data, generated_posterior)
    p = plot(r)
    savefig( p, normpath( fdir, "FIG-RT_"*fname*".html" ) )
end
##
let
    r = Regional.HospitsPlottingRecipe(data, generated_posterior)
    p = plot(r)
    savefig( p, normpath( fdir, "FIG-HOSPIT_"*fname*".html" ) )
end

#-----------------------------------------------------------------------------
# @info "store reproduction number"
# rt = let
#     rt = let
#         ib = findfirst(==(Date("2020-05-15")), data.dates)
#         ie = findfirst(==(Date(parsed_args["observations_end"])), data.dates)
#         Rt_array = hcat(generated_posterior.Rts...)[ib:ie,:]
#
#         qs = [quantile(v, [0.025, 0.25, 0.5, 0.75, 0.975]) for v in eachrow(Rt_array)]
#         llq, lq, mq, uq, uuq = (eachrow(hcat(qs...))..., )
#
#         date = data.dates[ib:ie]
#         DataFrame((;date, llq, lq, mq, uq, uuq))
#     end
#     saveto =  normpath( fdir, "Rt_"*fname*".csv" )
#     save(saveto, rt)
#     rt
# end
#-----------------------------------------------------------------------------
# perform diagnostics
if parsed_args["chains"] > 1
    @info "gelman diagnostics"
    diagnostics = gelmandiag(chain)
    saveto =  normpath( fdir, PREFIX*"GELMANDIAG_"*fname*".csv" )
    save( saveto, diagnostics)
    pretty_table(diagnostics; crop=:none)
end
#-----------------------------------------------------------------------------
@info "meanplot"
let
    p = meanplot(chain[n]);
    saveto =  normpath( fdir, PREFIX*"MEANPLOT_"*fname*".html" )
    savefig(p, saveto )
    run(`firefox $(saveto)`, wait=false)
end

max_stars_repo_name: andreaskoher/CovidSurvey

max_issues_repo_name: andreaskoher/CovidSurvey

max_forks_repo_name: andreaskoher/CovidSurvey


text:
# NIfTI.jl
# Methods for reading NIfTI MRI files in Julia

module NIfTI

using CodecZlib, Mmap, MappedArrays, TranscodingStreams

import Base.getindex, Base.size, Base.ndims, Base.length, Base.write, Base64
export NIVolume, niread, niwrite, voxel_size, time_step, vox, getaffine, setaffine

function define_packed(ty::DataType)
    packed_offsets = cumsum([sizeof(x) for x in ty.types])
    sz = pop!(packed_offsets)
    pushfirst!(packed_offsets, 0)

    @eval begin
        function Base.read(io::IO, ::Type{$ty})
            bytes = read!(io, Array{UInt8}(undef, $sz...))
            hdr = $(Expr(:new, ty, [:(unsafe_load(convert(Ptr{$(ty.types[i])}, pointer(bytes)+$(packed_offsets[i])))) for i = 1:length(packed_offsets)]...,))
            if hdr.sizeof_hdr == ntoh(Int32(348))
                return byteswap(hdr), true
            end
            hdr, false
        end
        function Base.write(io::IO, x::$ty)
            bytes = UInt8[]
            for name in fieldnames($ty)
                append!(bytes, reinterpret(UInt8, [getfield(x,name)]))
            end
            write(io, bytes)
            $sz
        end
    end
    nothing
end

mutable struct NIfTI1Header
    sizeof_hdr::Int32

    data_type::NTuple{10,UInt8}
    db_name::NTuple{18,UInt8}
    extents::Int32
    session_error::Int16
    regular::Int8

    dim_info::Int8
    dim::NTuple{8,Int16}
    intent_p1::Float32
    intent_p2::Float32
    intent_p3::Float32
    intent_code::Int16
    datatype::Int16
    bitpix::Int16
    slice_start::Int16
    pixdim::NTuple{8,Float32}
    vox_offset::Float32
    scl_slope::Float32
    scl_inter::Float32
    slice_end::Int16
    slice_code::Int8
    xyzt_units::Int8
    cal_max::Float32
    cal_min::Float32
    slice_duration::Float32
    toffset::Float32

    glmax::Int32
    glmin::Int32

    descrip::NTuple{80,UInt8}
    aux_file::NTuple{24,UInt8}

    qform_code::Int16
    sform_code::Int16
    quatern_b::Float32
    quatern_c::Float32
    quatern_d::Float32
    qoffset_x::Float32
    qoffset_y::Float32
    qoffset_z::Float32

    srow_x::NTuple{4,Float32}
    srow_y::NTuple{4,Float32}
    srow_z::NTuple{4,Float32}

    intent_name::NTuple{16,UInt8}

    magic::NTuple{4,UInt8}
end
define_packed(NIfTI1Header)

# byteswapping

function byteswap(hdr::NIfTI1Header)
    for fn in fieldnames(typeof(hdr))
        val = getfield(hdr, fn)
        if isa(val, Number) && sizeof(val) > 1
            setfield!(hdr, fn, ntoh(val))
        elseif isa(val, NTuple) && sizeof(eltype(val)) > 1
            setfield!(hdr, fn, map(ntoh, val))
        end
    end
    hdr
end

const SIZEOF_HDR = Int32(348)

const NIfTI_DT_BITSTYPES = Dict{Int16,Type}([
    (Int16(1), Bool),
    (Int16(2), UInt8),
    (Int16(4), Int16),
    (Int16(8), Int32),
    (Int16(16), Float32),
    (Int16(32), ComplexF32),
    (Int16(64), Float64),
    (Int16(256), Int8),
    (Int16(512), UInt16),
    (Int16(768), UInt32),
    (Int16(1024), Int64),
    (Int16(1280), UInt64),
    (Int16(1792), ComplexF64)
])
const NIfTI_DT_BITSTYPES_REVERSE = Dict{Type,Int16}()
for (k, v) in NIfTI_DT_BITSTYPES
    NIfTI_DT_BITSTYPES_REVERSE[v] = k
end

const NP1_MAGIC = (0x6e,0x2b,0x31,0x00)
const NI1_MAGIC = (0x6e,0x69,0x31,0x00)

function string_tuple(x::String, n::Int)
    a = codeunits(x)
    padding = zeros(UInt8, n-length(a))
    (a..., padding...)
end
string_tuple(x::AbstractString) = string_tuple(bytestring(x))

mutable struct NIfTI1Extension
    ecode::Int32
    edata::Vector{UInt8}
end

mutable struct NIVolume{T<:Number,N,R} <: AbstractArray{T,N}
    header::NIfTI1Header
    extensions::Vector{NIfTI1Extension}
    raw::R

end
NIVolume(header::NIfTI1Header, extensions::Vector{NIfTI1Extension}, raw::R) where {R}=
    niupdate(new(header, extensions, raw))

NIVolume(header::NIfTI1Header, extensions::Vector{NIfTI1Extension}, raw::AbstractArray{T,N}) where {T<:Number,N} =
    NIVolume{typeof(one(T)*1f0+1f0),N,typeof(raw)}(header, extensions, raw)
NIVolume(header::NIfTI1Header, raw::AbstractArray{T,N}) where {T<:Number,N} =
    NIVolume{typeof(one(T)*1f0+1f0),N,typeof(raw)}(header, NIfTI1Extension[], raw)
NIVolume(header::NIfTI1Header, extensions::Vector{NIfTI1Extension}, raw::AbstractArray{Bool,N}) where {N} =
    NIVolume{Bool,N,typeof(raw)}(header, extensions, raw)
NIVolume(header::NIfTI1Header, raw::AbstractArray{Bool,N}) where {N} =
    NIVolume{Bool,N,typeof(raw)}(header, NIfTI1Extension[], raw)

# Conversion factors to mm/ms
# http://nifti.nimh.nih.gov/nifti-1/documentation/nifti1fields/nifti1fields_pages/xyzt_units.html
const SPATIAL_UNIT_MULTIPLIERS = [
    1000,   # 1 => NIfTI_UNITS_METER
    1,      # 2 => NIfTI_UNITS_MM
    0.001   # 3 => NIfTI_UNITS_MICRON
]
const TIME_UNIT_MULTIPLIERS = [
    1000,   # NIfTI_UNITS_SEC
    1,      # NIfTI_UNITS_MSEC
    0.001,  # NIfTI_UNITS_USEC
    1,      # NIfTI_UNITS_HZ
    1,      # NIfTI_UNITS_PPM
    1       # NIfTI_UNITS_RADS
]

# Always in mm
voxel_size(header::NIfTI1Header) =
    [header.pixdim[i] * SPATIAL_UNIT_MULTIPLIERS[header.xyzt_units & Int8(3)] for i = 2:min(header.dim[1], 3)+1]

# Always in ms
time_step(header::NIfTI1Header) =
    header.pixdim[5] * TIME_UNIT_MULTIPLIERS[header.xyzt_units >> 3]

function to_dim_info(dim_info::Tuple{Integer,Integer,Integer})
    if dim_info[1] > 3 || dim_info[1] < 0
        error("Invalid frequency dimension $(dim_info[1])")
    elseif dim_info[2] > 3 || dim_info[2] < 0
        error("Invalid phase dimension $(dim_info[2])")
    elseif dim_info[3] > 3 || dim_info[3] < 0
        error("Invalid slice dimension $(dim_info[3])")
    end

    Int8(dim_info[1] | (dim_info[2] << 2) | (dim_info[3] << 4))
end

# Returns or sets dim_info as a tuple whose values are the frequency, phase, and slice dimensions
dim_info(header::NIfTI1Header) = (header.dim_info & int8(3), (header.dim_info >> 2) & int8(3),
    (header.dim_info >> 4) & int8(3))
dim_info(header::NIfTI1Header, dim_info::Tuple{T, T, T}) where {T<:Integer} =
    header.dim_info = to_dim_info(dim_info)

# Gets dim to be used in header
function nidim(x::AbstractArray)
    dim = ones(Int16, 8)
    dim[1] = ndims(x)
    dim[2:dim[1]+1] = [size(x)...]
    (dim...,)
end

# Gets datatype to be used in header
function nidatatype(t::Type)
    t = get(NIfTI_DT_BITSTYPES_REVERSE, t, nothing)
    if t == nothing
        error("Unsupported data type $t")
    end
    t
end

# Gets the size of a type in bits
nibitpix(t::Type) = Int16(sizeof(t)*8)
nibitpix(::Type{Bool}) = Int16(1)

# Convert a NIfTI header to a 4x4 affine transformation matrix
function getaffine(h::NIfTI1Header)
    pixdim = convert(Tuple{Vararg{Float64}}, h.pixdim)
    # See documentation at
    # http://nifti.nimh.nih.gov/nifti-1/documentation/nifti1fields/nifti1fields_pages/qsform.html
    if h.sform_code > 0
        # Method 3
        return Float64[
            collect(h.srow_x)'
            collect(h.srow_y)'
            collect(h.srow_z)'
            0 0 0 1
        ]
    elseif h.qform_code > 0
        # Method 2
        # http://nifti.nimh.nih.gov/nifti-1/documentation/nifti1fields/nifti1fields_pages/quatern.html
        b::Float64 = h.quatern_b
        c::Float64 = h.quatern_c
        d::Float64 = h.quatern_d
        a = sqrt(1 - b*b - c*c - d*d)
        A = [
            a*a+b*b-c*c-d*d   2*b*c-2*a*d       2*b*d+2*a*c
            2*b*c+2*a*d       a*a+c*c-b*b-d*d   2*c*d-2*a*b
            2*b*d-2*a*c       2*c*d+2*a*b       a*a+d*d-c*c-b*b
        ]
        B = [
            pixdim[2]
            pixdim[3]
            pixdim[1]*pixdim[4]
        ]
        return vcat(hcat(A*B, [
            h.qoffset_x
            h.qoffset_y
            h.qoffset_z
        ]), [0 0 0 1])
    elseif h.qform_code == 0
        # Method 1
        return Float64[
             pixdim[2] 0         0          0 0
                       0 pixdim[3]          0 0
                       0         0  bitpix[4] 0
                       0         0          0 1
        ]
    end
end

# Set affine matrix of NIfTI header
function setaffine(h::NIfTI1Header, affine::Array{T,2}) where {T}
    size(affine, 1) == size(affine, 2) == 4 ||
        error("affine matrix must be 4x4")
    affine[4, 1] == affine[4, 2] == affine[4, 3] == 0 && affine[4, 4] == 1 ||
        error("last row of affine matrix must be [0 0 0 1]")
    h.qform_code = one(Int16)
    h.sform_code = one(Int16)
    h.pixdim = (zero(Float32), h.pixdim[2:end]...,)
    h.quatern_b = zero(Float32)
    h.quatern_c = zero(Float32)
    h.quatern_d = zero(Float32)
    h.qoffset_x = zero(Float32)
    h.qoffset_y = zero(Float32)
    h.qoffset_z = zero(Float32)
    h.srow_x = convert(Tuple{Vararg{Float32}}, (affine[1, :]...,))
    h.srow_y = convert(Tuple{Vararg{Float32}}, (affine[2, :]...,))
    h.srow_z = convert(Tuple{Vararg{Float32}}, (affine[3, :]...,))
    h
end

# Constructor
function NIVolume(
    # Optional MRI volume; if not given, an empty volume is used
    raw::AbstractArray{T}=Int16[],
    extensions::Union{Vector{NIfTI1Extension},Nothing}=nothing;

    # Fields specified as UNUSED in NIfTI1 spec
    data_type::AbstractString="", db_name::AbstractString="", extents::Integer=Int32(0),
    session_error::Integer=Int16(0), regular::Integer=Int8(0), glmax::Integer=Int32(0),
    glmin::Integer=Int16(0),

    # The frequency encoding, phase encoding and slice dimensions.
    dim_info::NTuple{3, Integer}=(0, 0, 0),
    # Describes data contained in the file; for valid values, see
    # http://nifti.nimh.nih.gov/nifti-1/documentation/nifti1fields/nifti1fields_pages/group__NIfTI1__INTENT__CODES.html
    intent_p1::Real=0f0, intent_p2::Real=0f0, intent_p3::Real=0f0,
    intent_code::Integer=Int16(0), intent_name::AbstractString="",
    # Information about which slices were acquired, in case the volume has been padded
    slice_start::Integer=Int16(0), slice_end::Integer=Int16(0), slice_code::Int8=Int8(0),
    # The size of each voxel and the time step. These are formulated in mm unless xyzt_units is
    # explicitly specified
    voxel_size::NTuple{3, Real}=(1f0, 1f0, 1f0), time_step::Real=0f0, xyzt_units::Int8=Int8(18),
    # Slope and intercept by which volume shoudl be scaled
    scl_slope::Real=1f0, scl_inter::Real=0f0,
    # These describe how data should be scaled when displayed on the screen. They are probably
    # rarely used
    cal_max::Real=0f0, cal_min::Real=0f0,
    # The amount of time taken to acquire a slice
    slice_duration::Real=0f0,
    # Indicates a non-zero start point for time axis
    toffset::Real=0f0,

    # "Any text you like"
    descrip::AbstractString="",
    # Name of auxiliary file
    aux_file::AbstractString="",

    # Parameters for Method 2. See the NIfTI spec
    qfac::Float32=0f0, quatern_b::Real=0f0, quatern_c::Real=0f0, quatern_d::Real=0f0,
    qoffset_x::Real=0f0, qoffset_y::Real=0f0, qoffset_z::Real=0f0,
    # Orientation matrix for Method 3
    orientation::Union{Matrix{Float32},Nothing}=nothing) where {T <: Number}
    local t
    if isempty(raw)
        raw = Int16[]
        t = Int16
    else
        t = T
    end

    if extensions == nothing
        extensions = NIfTI1Extension[]
    end

    method2 = qfac != 0 || quatern_b != 0 || quatern_c != 0 || quatern_d != 0 || qoffset_x != 0 ||
        qoffset_y != 0 || qoffset_z != 0
    method3 = orientation != nothing

    if method2 && method3
        error("Orientation parameters for Method 2 and Method 3 are mutually exclusive")
    end

    if method3
        if size(orientation) != (3, 4)
            error("Orientation matrix must be of dimensions (3, 4)")
        end
    else
        orientation = zeros(Float32, 3, 4)
    end

    if slice_start == 0 && slice_end == 0 && dim_info[3] != 0
        slice_start = 0
        slice_end = size(raw, dim_info[3]) - 1
    end

    NIVolume(NIfTI1Header(SIZEOF_HDR, string_tuple(data_type, 10), string_tuple(db_name, 18), extents, session_error,
        regular, to_dim_info(dim_info), nidim(raw), intent_p1, intent_p2,
        intent_p3, intent_code, nidatatype(t), nibitpix(t),
        slice_start, (qfac, voxel_size..., time_step, 0, 0, 0), 352,
        scl_slope, scl_inter, slice_end, slice_code,
        xyzt_units, cal_max, cal_min, slice_duration,
        toffset, glmax, glmin, string_tuple(descrip, 80), string_tuple(aux_file, 24), (method2 || method3),
        method3, quatern_b, quatern_c, quatern_d,
        qoffset_x, qoffset_y, qoffset_z, (orientation[1, :]...,),
        (orientation[2, :]...,), (orientation[3, :]...,), string_tuple(intent_name, 16), NP1_MAGIC), extensions, raw)
end

# Calculates the size of a NIfTI extension
esize(ex::NIfTI1Extension) = 8 + ceil(Int, length(ex.edata)/16)*16

# Validates the header of a volume and updates it to match the volume's contents
function niupdate(vol::NIVolume{T}) where {T}
    vol.header.sizeof_hdr = SIZEOF_HDR
    vol.header.dim = nidim(vol.raw)
    vol.header.datatype = nidatatype(T)
    vol.header.bitpix = nibitpix(T)
    vol.header.vox_offset = isempty(vol.extensions) ? Int32(352) :
        Int32(mapreduce(esize, +, vol.extensions) + SIZEOF_HDR)
    vol
end

# Avoid method ambiguity
write(io::Base64.Base64EncodePipe, vol::NIVolume{UInt8,1}) = invoke(write, (IO, NIVolume{UInt8,1}), io, vol)

# Write a NIfTI file
function write(io::IO, vol::NIVolume)
    write(io, niupdate(vol).header)
    if isempty(vol.extensions)
        write(io, Int32(0))
    else
        for ex in vol.extensions
            sz = esize(ex)
            write(io, Int32(sz))
            write(io, Int32(ex.ecode))
            write(io, ex.edata)
            write(io, zeros(UInt8, sz - length(ex.edata)))
        end
    end
    if eltype(vol.raw) == Bool
        write(io, BitArray(vol.raw))
    else
        write(io, vol.raw)
    end
end

# Convenience function to write a NIfTI file given a path
function niwrite(path::AbstractString, vol::NIVolume)
    if split(path,".")[end] == "gz"
        io = open(path, "w9")
        stream = GzipCompressorStream(io)
        write(stream, vol)
        close(stream)
        close(io)
    else
        io = open(path, "w")
        write(io, vol)
        close(io)
    end
end

# Read header from a NIfTI file
function read_header(io::IO)
    header, swapped = read(io, NIfTI1Header)
    if header.sizeof_hdr != SIZEOF_HDR
        error("This is not a NIfTI-1 file")
    end
    header, swapped
end

# Read extension fields following NIfTI header
function read_extensions(io::IO, header::NIfTI1Header)
    if eof(io)
        return NIfTI1Extension[]
    end

    extension = read!(io, Array{UInt8}(undef, 4))
    if extension[1] != 1
        return NIfTI1Extension[]
    end

    extensions = NIfTI1Extension[]
    should_bswap = header.dim[1] > 7
    while header.magic == NP1_MAGIC ? position(io) < header.vox_offset : !eof(io)
        esize = read(io, Int32)
        ecode = read(io, Int32)

        if should_bswap
            esize = bswap(esize)
            ecode = bswap(ecode)
        end

        push!(extensions, NIfTI1Extension(ecode, read!(io, Array{UInt8}(undef, esize-8))))
    end
    extensions
end

# Look for a gzip header in an IOStream
function isgz(io::IO)
    ret = read(io, UInt8) == 0x1F && read(io, UInt8) == 0x8B
    seek(io, 0)
    ret
end

# Read a NIfTI file. The optional mmap argument determines whether the contents are read in full
# (if false) or mmaped from the disk (if true).
"""

"""
function niread(file::AbstractString; mmap::Bool=false, mode::AbstractString="r")
    file_io = open(file, mode)
    header_gzipped = isgz(file_io)
    header_io = header_gzipped ? GzipDecompressorStream(file_io) : file_io
    header, swapped = read_header(header_io)
    extensions = read_extensions(header_io, header)
    dims = convert(Tuple{Vararg{Int}}, header.dim[2:header.dim[1]+1])

    if !haskey(NIfTI_DT_BITSTYPES, header.datatype)
        error("data type $(header.datatype) not yet supported")
    end
    dtype = NIfTI_DT_BITSTYPES[header.datatype]

    ArrayType = if dtype == Bool
            BitArray{length(dims)}
        else
            Array{dtype, length(dims)}
        end

    local volume
    if header.magic == NP1_MAGIC
        if mmap
            if header_gzipped
                close(header_io)
                close(file_io)
                error("cannot mmap a gzipped NIfTI file")
            else
                volume = Mmap.mmap(header_io, ArrayType, dims, Int(header.vox_offset))
            end
        else
            seekstart(header_io)
            read(header_io, Int(header.vox_offset))
            volume = read!(header_io, ArrayType(undef, dims))
            if !eof(header_io)
                println("Warning! File size does not match length of data; some data may be ignored")
            end
            close(header_io)
            !header_gzipped || close(file_io)
        end
    elseif header.magic == NI1_MAGIC
        close(header_io)
        !header_gzipped || close(file_io)

        volume_name = replace(file, r"\.\w+(\.gz)?$" => "")*".img"
        if !isfile(volume_name)
            if isfile(volume_name*".gz")
                volume_name *= ".gz"
            else
                error("NIfTI file is dual file storage, but $volume_name does not exist")
            end
        end

        volume_io = open(volume_name, mode)
        volume_gzipped = isgz(volume_io)
        if mmap
            if volume_gzipped
                close(volume_io)
                error("cannot mmap a gzipped NIfTI file")
            else
                volume = Mmap.mmap(volume_io, ArrayType, dims)
            end
        else
            if volume_gzipped
                volume_gz_io = GzipDecompressorStream(volume_io)
                volume = read!(volume_gz_io, ArrayType(undef, dims))
                close(volume_gz_io)
            else
                volume = read!(volume_io, ArrayType(undef, dims))
            end
            close(volume_io)
        end
    end
    if swapped && sizeof(eltype(volume)) > 1
        volume = mappedarray(ntoh, hton, volume)
    end

    return NIVolume(header, extensions, volume)
end

# Allow file to be indexed like an array, but with indices yielding scaled data
@inline getindex(f::NIVolume{T}, idx::Vararg{Int}) where {T} =
    getindex(f.raw, idx...,) * f.header.scl_slope + f.header.scl_inter

add1(x::Union{AbstractArray{T},T}) where {T<:Integer} = x + 1
add1(::Colon) = Colon()
@inline vox(f::NIVolume, args...,) = getindex(f, map(add1, args)...,)
size(f::NIVolume) = size(f.raw)
size(f::NIVolume, d) = size(f.raw, d)
ndims(f::NIVolume) = ndims(f.raw)
length(f::NIVolume) = length(f.raw)
lastindex(f::NIVolume) = lastindex(f.raw)

end


text:
module EodData
#=
    References
     ----------

    Julia language documentation
    - http://julia.readthedocs.org/en/latest/manual/
    - http://julia.readthedocs.org/en/latest/stdlib/base/
    - http://learnxinyminutes.com/docs/julia/
    - http://www.scolvin.com/juliabyexample/
    - http://perldoc.perl.org/perlre.html (regular expressions : regex)

    Hypertext Transfer Protocol
    - http://en.wikipedia.org/wiki/Hypertext_Transfer_Protocol#Example_session

     EODData now offers a complete end of day Web Service ideal for trading
    applications, web sites, portfolio management systems, etc.
    - http://ws.eoddata.com/data.asmx
    - http://ws.eoddata.com/data.asmx?wsdl
    - http://eoddata.com/products/webservice.aspx

     Package: HTTPClient
    - https://github.com/JuliaWeb/HTTPClient.jl
    - https://github.com/JuliaWeb/HTTPClient.jl/blob/master/src/HTTPC.jl
    - https://github.com/JuliaWeb/HTTPClient.jl/blob/master/test/tests.jl

    Package: LibExpat
    - https://github.com/amitmurthy/LibExpat.jl
    - http://nbviewer.ipython.org/github/amitmurthy/LibExpat.jl/blob/master/libexpat_test.ipynb
=#
using DataStructures
if VERSION < v"0.4-"
    using Dates
else
    using Base.Dates
end
using HTTPClient.HTTPC
using LibExpat

# ==================
# Internal constants
const DATETIMEFORMAT_SS = "yyyy-mm-ddTHH:MM:SS"
const DATETIMEFORMAT_MS = "yyyy-mm-ddTHH:MM:SS.ss"
const HEADER_DELIMITERS = [',', ';', ' ']

# ============================
# Make EodData types available
export DataFormat, DataFormatColumn, Exchange, Fundamental, LoginResponse,
        Quote, Quote_2, Split, TickerChange, Ticker

# =============
# EodData Types

# ::Country
type Country
    code::String
    name::String

    # Default constructor
    Country(code::String, name::String) = new(code, name)

    # Construct from XML ::ETree
    function Country(c::ETree)
        code::String = strip(get(c.attr,"Code",""))
        name::String = strip(get(c.attr,"Name",""))

        return new(code, name)
    end
end

# ::DataFormatColumn
type DataFormatColumn
    column_header::String
    sort_order::Int
    data_format_code::String
    data_format_name::String
    column_code::String
    column_name::String
    column_type_id::Int
    column_type::String

    # Default constructor
    DataFormatColumn(column_header::String, sort_order::Int, data_format_code::String,
                     data_format_name::String, column_code::String, column_name::String,
                     column_type_id::Int, column_type::String) =
        new(column_header, sort_order, data_format_code, data_format_name, column_code,
            column_name, column_type_id, column_type)

    # Construct from XML ::ETree
    function DataFormatColumn(col::ETree)
        column_header::String = strip(get(col.attr,"Header",""))
        sort_order::Int = int(strip(get(col.attr,"SortOrder","")))
        data_format_code::String = strip(get(col.attr,"Code",""))
        data_format_name::String = strip(get(col.attr,"DataFormat",""))
        column_code::String = strip(get(col.attr,"ColumnCode",""))
        column_name::String = strip(get(col.attr,"ColumnName",""))
        column_type_id::Int = int(strip(get(col.attr,"ColumnTypeId","")))
        column_type::String = strip(get(col.attr,"ColumnType",""))

        return new(column_header, sort_order, data_format_code, data_format_name, column_code,
                   column_name, column_type_id, column_type)
    end
end

# ::DataFormat
type DataFormat
    code::String
    name::String
    header_format::Vector{String}
    date_format::String
    extension::String
    include_suffix::Bool
    tab_column_seperator::Bool
    column_seperator::String
    text_qualifier::String
    filename_prefix::String
    filename_exchange_code::Bool
    filename_date::Bool
    include_header_row::Bool
    hour_format::String
    datetime_seperator::String
    exchange_filename_format_date::String
    exchange_filename_format_date_range::String
    ticker_filename_format_date::String
    ticker_filename_format_date_range::String
    columns::OrderedDict{Int, DataFormatColumn}

    # Default constructor
    DataFormat(code::String, name::String, header_format::Vector{String}, date_format::String,
               extension::String, include_suffix::Bool, tab_column_seperator::Bool,
               column_seperator::String, text_qualifier::String, filename_prefix::String,
               filename_exchange_code::String, filename_date::Bool, include_header_row::Bool,
               hour_format::String, datetime_seperator::String, exchange_filename_format_date::String,
               exchange_filename_format_date_range::String, ticker_filename_format_date::String,
               ticker_filename_format_date_range::String, columns::OrderedDict{Int, DataFormatColumn}) =
        new(code, name, header_format, date_format, extension, include_suffix, tab_column_seperator,
            column_seperator, text_qualifier, filename_prefix, filename_exchange_code, filename_date,
            include_header_row, hour_format, datetime_seperator, exchange_filename_format_date,
            exchange_filename_format_date_range, ticker_filename_format_date,
            ticker_filename_format_date_range, columns)

    # Construct from XML ::ETree
    function DataFormat(df::ETree)
        code::String = strip(get(df.attr,"Code",""))
        name::String = strip(get(df.attr,"Name",""))
        header_format::Vector{String} = convert(Array{String}, split(strip(get(df.attr,"Header","")), HEADER_DELIMITERS))
        date_format::String = strip(get(df.attr,"DateFormat",""))
        extension::String = strip(get(df.attr,"Extension",""))
        include_suffix::Bool = lowercase(strip(get(df.attr,"IncludeSuffix",""))) == "true" ? true : false
        tab_column_seperator::Bool = lowercase(strip(get(df.attr,"TabColumnSeperator","")))  == "true" ? true : false
        column_seperator::String = strip(get(df.attr,"ColumnSeperator",""))
        text_qualifier::String = strip(get(df.attr,"TextQualifier",""))
        filename_prefix::String = strip(get(df.attr,"FilenamePrefix",""))
        filename_exchange_code::Bool = lowercase(strip(get(df.attr,"FilenameExchangeCode",""))) == "true" ? true : false
        filename_date::Bool = lowercase(strip(get(df.attr,"FilenameDate",""))) == "true" ? true : false
        include_header_row::Bool = lowercase(strip(get(df.attr,"IncludeHeaderRow",""))) == "true" ? true : false
        hour_format::String = strip(get(df.attr,"HourFormat",""))
        datetime_seperator::String = strip(get(df.attr,"DateTimeSeperator",""))
        exchange_filename_format_date::String = strip(get(df.attr,"ExchangeFilenameFormatDate",""))
        exchange_filename_format_date_range::String = strip(get(df.attr,"ExchangeFilenameFormatDateRange",""))
        ticker_filename_format_date::String = strip(get(df.attr,"SymbolFilenameFormatDate",""))
        ticker_filename_format_date_range::String = strip(get(df.attr,"SymbolFilenameFormatDateRange",""))

        columns = OrderedDict{Int, DataFormatColumn}()
        for col_xml in find(df, "COLUMNS/DATAFORMAT_COLUMN")
            columns[col.sort_order] = col = DataFormatColumn(col_xml)
        end

        return new(code, name, header_format, date_format, extension, include_suffix, tab_column_seperator,
                   column_seperator, text_qualifier, filename_prefix, filename_exchange_code, filename_date,
                   include_header_row, hour_format, datetime_seperator, exchange_filename_format_date,
                   exchange_filename_format_date_range, ticker_filename_format_date,
                   ticker_filename_format_date_range, columns)
    end
end

# ::Exchange
type Exchange
    code::String
    name::String
    last_trade_date_time::DateTime
    country_code::String
    currency_code::String
    advances::Float64
    declines::Float64
    suffix::String
    time_zone::String
    is_intraday::Bool
    intraday_start_date::DateTime
    has_intraday_product::Bool

    # Default constructor
    Exchange(code::String, name::String, last_trade_date_time::DateTime, country_code::String,
             currency_code::String, advances::Float64, declines::Float64,suffix::String,
             time_zone::String, is_intraday::Bool, intraday_start_date::DateTime,
             has_intraday_product::Bool) =
        new(code, name, last_trade_date_time, country_code, currency_code, advances, declines,
            suffix, time_zone, is_intraday, intraday_start_date, has_intraday_product)

    # Construct from XML ::ETree
    function Exchange(ex::ETree)
        code::String = strip(get(ex.attr,"Code",""))
        name::String = strip(get(ex.attr,"Name",""))
        last_trade_date_time::DateTime  = DateTime(strip(get(ex.attr,"LastTradeDateTime",""))[1:19], DATETIMEFORMAT_SS)
        country_code::String = strip(get(ex.attr,"Country",""))
        currency_code::String = strip(get(ex.attr,"Currency",""))
        advances::Float64 = float(strip(get(ex.attr,"Advances","")))
        declines::Float64 = float(strip(get(ex.attr,"Declines","")))
        suffix::String = strip(get(ex.attr,"Suffix",""))
        time_zone::String = strip(get(ex.attr,"TimeZone",""))
        is_intraday::Bool = lowercase(strip(get(ex.attr,"IsIntraday",""))) == "true" ? true : false
        intraday_start_date::DateTime = DateTime(strip(get(ex.attr,"IntradayStartDate",""))[1:19], DATETIMEFORMAT_SS)
        has_intraday_product::Bool = lowercase(strip(get(ex.attr,"HasIntradayProduct",""))) == "true" ? true : false

        return new(code, name, last_trade_date_time, country_code, currency_code, advances,
                   declines,suffix, time_zone, is_intraday, intraday_start_date, has_intraday_product)
    end
end

# ::Fundamental
type Fundamental
    ticker_code::String
    name::String
    description::String
    date_time::DateTime
    industry::String
    sector::String
    share_float::Float64
    market_cap::Float64
    pe_ratio::Float64
    earnings_per_share::Float64
    net_tangible_assets::Float64
    dividend_yield::Float64
    dividend::Float64
    dividend_date::DateTime
    dividend_per_share::Float64
    imputation_credits::Float64
    ebitda::Float64
    peg_ratio::Float64
    ps_ratio::Float64
    pb_ratio::Float64
    yield::Float64

    # Default constructor
    Fundamental(ticker_code::String, name::String, description::String, date_time::DateTime,
                industry::String, sector::String, share_float::Float64, market_cap::Float64,
                pe_ratio::Float64, earnings_per_share::Float64, net_tangible_assets::Float64,
                dividend_yield::Float64, dividend::Float64, dividend_date::DateTime,
                dividend_per_share::Float64, imputation_credits::Float64, ebitda::Float64,
                peg_ratio::Float64, ps_ratio::Float64, pb_ratio::Float64, yield::Float64) =
        new(ticker_code, name, description, date_time, industry, sector, share_float,
            market_cap, pe_ratio, earnings_per_share, net_tangible_assets, dividend_yield,
            dividend, dividend_date, dividend_per_share, imputation_credits, ebitda,
            peg_ratio, ps_ratio, pb_ratio, yield)

    # Construct from XML ::ETree
    function Fundamental(fl::ETree)
        ticker_code::String = strip(get(fl.attr,"Symbol",""))
        name::String = strip(get(fl.attr,"Name",""))
        description::String = strip(get(fl.attr,"Description",""))
        date_time::DateTime = DateTime(strip(get(fl.attr,"DateTime",""))[1:19], DATETIMEFORMAT_SS)
        industry::String = strip(get(fl.attr,"Industry",""))
        sector::String = strip(get(fl.attr,"Sector",""))
        share_float::Float64 = float(strip(get(fl.attr,"Shares","")))
        market_cap::Float64 = float(strip(get(fl.attr,"MarketCap","")))
        pe_ratio::Float64 = float(strip(get(fl.attr,"PE","")))
        earnings_per_share::Float64 = float(strip(get(fl.attr,"EPS","")))
        net_tangible_assets::Float64 = float(strip(get(fl.attr,"NTA","")))
        dividend_yield::Float64 = float(strip(get(fl.attr,"DivYield","")))
        dividend::Float64 = float(strip(get(fl.attr,"Dividend","")))
        dividend_date::DateTime = DateTime(strip(get(fl.attr,"DividendDate",""))[1:19], DATETIMEFORMAT_SS)
        dividend_per_share::Float64 = float(strip(get(fl.attr,"DPS","")))
        imputation_credits::Float64 = float(strip(get(fl.attr,"ImputationCredits","")))
        ebitda::Float64 = float(strip(get(fl.attr,"EBITDA","")))
        peg_ratio::Float64 = float(strip(get(fl.attr,"PEG","")))
        ps_ratio::Float64 = float(strip(get(fl.attr,"PtS","")))
        pb_ratio::Float64 = float(strip(get(fl.attr,"PtB","")))
        yield::Float64 = float(strip(get(fl.attr,"Yield","")))

        return new(ticker_code, name, description, date_time, industry, sector, share_float,
                   market_cap, pe_ratio, earnings_per_share, net_tangible_assets, dividend_yield,
                   dividend, dividend_date, dividend_per_share, imputation_credits, ebitda,
                   peg_ratio, ps_ratio, pb_ratio, yield)
    end
end

# ::LoginResponse
type LoginResponse
    message::String
    token::String

    # Default constructor
    LoginResponse(message::String, token::String) = new(message, token)

    # Construct from XML ::ETree
    function LoginResponse(lr::ETree)
        message = strip(find(lr, "/LOGINRESPONSE[1]{Message}"))
        token = strip(find(lr, "/LOGINRESPONSE[1]{Token}"))

        return new(message, token)
    end
end

# ::Quote
type Quote
    ticker_code::String
    description::String
    name::String
    date_time::DateTime
    open::Float64
    high::Float64
    low::Float64
    close::Float64
    volume::Float64
    open_interest::Float64
    previous::Float64
    change::Float64
    simple_return::Float64
    bid::Float64
    ask::Float64
    previous_close::Float64
    next_open::Float64
    modified::DateTime

    # Default constructor
    Quote(ticker_code::String, description::String, name::String, date_time::DateTime, open::Float64,
          high::Float64, low::Float64, close::Float64, volume::Float64, open_interest::Float64,
          previous::Float64, change::Float64, simple_return::Float64, bid::Float64, ask::Float64,
          previous_close::Float64, next_open::Float64, modified::DateTime) =
        new(ticker_code, description, name, date_time, open, high, low, close, volume, open_interest,
                   previous, change, simple_return, bid, ask, previous_close, next_open, modified)

    # Construct from XML ::ETree
    function Quote(qt::ETree)
        ticker_code::String = strip(get(qt.attr,"Symbol",""))
        description::String = strip(get(qt.attr,"Description",""))
        name::String = strip(get(qt.attr,"Name",""))
        date_time::DateTime = DateTime(strip(get(qt.attr,"DateTime",""))[1:19], DATETIMEFORMAT_SS)
        open::Float64 = float(strip(get(qt.attr,"Open","")))
        high::Float64 = float(strip(get(qt.attr,"High","")))
        low::Float64 = float(strip(get(qt.attr,"Low","")))
        close::Float64 = float(strip(get(qt.attr,"Close","")))
        volume::Float64 = float(strip(get(qt.attr,"Volume","")))
        open_interest::Float64 = float(strip(get(qt.attr,"OpenInterest","")))
        previous::Float64 = float(strip(get(qt.attr,"Previous","")))
        change::Float64 = float(strip(get(qt.attr,"Change","")))
        simple_return::Float64 = change / previous
        bid::Float64 = float(strip(get(qt.attr,"Bid","")))
        ask::Float64 = float(strip(get(qt.attr,"Ask","")))
        previous_close::Float64 = float(strip(get(qt.attr,"PreviousClose","")))
        next_open::Float64 = float(strip(get(qt.attr,"NextOpen","")))
        modified::DateTime = DateTime(strip(get(qt.attr,"Modified",""))[1:19], DATETIMEFORMAT_SS)

        return new(ticker_code, description, name, date_time, open, high, low, close, volume,
                   open_interest, previous, change, simple_return, bid, ask, previous_close,
                   next_open, modified)
    end
end

# ::Quote_2
type Quote_2
    ticker_code::String
    date_time::DateTime
    open::Float64
    high::Float64
    low::Float64
    close::Float64
    volume::Float64
    open_interest::Float64
    bid::Float64
    ask::Float64

    # Default constructor
    Quote_2(ticker_code::String, date_time::DateTime, open::Float64, high::Float64, low::Float64,
            close::Float64, volume::Float64, open_interest::Float64, bid::Float64, ask::Float64) =
        new(ticker_code, date_time, open, high, low, close, volume, open_interest, bid, ask)

    # Construct from XML ::ETree
    function Quote_2(qt_2::ETree)
        ticker_code::String = strip(get(qt_2.attr,"s",""))
        date_time::DateTime = DateTime(strip(get(qt_2.attr,"d",""))[1:19], DATETIMEFORMAT_SS)
        open::Float64 = float(strip(get(qt_2.attr,"o","")))
        high::Float64 = float(strip(get(qt_2.attr,"h","")))
        low::Float64 = float(strip(get(qt_2.attr,"l","")))
        close::Float64 = float(strip(get(qt_2.attr,"c","")))
        volume::Float64 = float(strip(get(qt_2.attr,"v","")))
        open_interest::Float64 = float(strip(get(qt_2.attr,"i","")))
        bid::Float64 = float(strip(get(qt_2.attr,"b","")))
        ask::Float64 = float(strip(get(qt_2.attr,"a","")))

        return new(ticker_code, date_time, open, high, low, close, volume, open_interest, bid, ask)
    end
end

# ::Split
type Split
    exchange_code::String
    ticker_code::String
    date_time::DateTime
    ratio::String
    price_multiplier::Float64
    share_float_multiplier::Float64
    is_reverse_split::Bool

    # Default constructor
    Split(exchange_code::String, ticker_code::String, date_time::DateTime, ratio::String,
          price_multiplier::Float64, share_float_multiplier::Float64, is_reverse_split::Bool) =
        new(exchange_code, ticker_code, date_time, ratio, price_multiplier,
            share_float_multiplier, is_reverse_split)

    # Construct from XML ::ETree
    function Split(sp::ETree)
        exchange_code::String = strip(get(sp.attr,"Exchange",""))
        ticker_code::String = strip(get(sp.attr,"Symbol",""))
        date_time::DateTime = DateTime(strip(get(sp.attr,"DateTime",""))[1:19], DATETIMEFORMAT_SS)
        ratio::String = strip(get(sp.attr,"Ratio",""))
        price_multiplier::Float64 = float(strip(split(ratio, "-")[2])) / float(strip(split(ratio, "-")[1]))
        share_float_multiplier::Float64 = float(strip(split(ratio, "-")[1])) / float(strip(split(ratio, "-")[2]))
        is_reverse_split::Bool = price_multiplier > 1.0 ? true : false

        new(exchange_code, ticker_code, date_time, ratio, price_multiplier,
            share_float_multiplier, is_reverse_split)
    end
end

# ::Technical
type Technical
    ticker_code::String
    name::String
    description::String
    date_time::DateTime
    previous::Float64
    change::Float64
    ma_1::Float64
    ma_2::Float64
    ma_5::Float64
    ma_20::Float64
    ma_50::Float64
    ma_100::Float64
    ma_200::Float64
    ma_percent::Float64
    ma_return::Float64
    volume_change::Float64
    three_month_change::Float64
    six_month_change::Float64
    week_high::Float64
    week_low::Float64
    week_change::Float64
    avg_week_change::Float64
    avg_week_volume::Float64
    week_volume::Float64
    month_high::Float64
    month_low::Float64
    month_change::Float64
    avg_month_change::Float64
    avg_month_volume::Float64
    month_volume::Float64
    year_high::Float64
    year_low::Float64
    year_change::Float64
    avg_year_change::Float64
    avg_year_volume::Float64
    ytd_change::Float64
    rsi_14::Float64
    sto_9::Float64
    wpr_14::Float64
    mtm_14::Float64
    roc_14::Float64
    ptc::Float64
    sar::Float64
    volatility::Float64
    liquidity::Float64

    # Default constructor
    Technical(ticker_code::String, name::String, description::String, date_time::DateTime,
              previous::Float64, change::Float64, ma_1::Float64, ma_2::Float64, ma_5::Float64,
              ma_20::Float64, ma_50::Float64, ma_100::Float64, ma_200::Float64, ma_percent::Float64,
              ma_return::Float64, volume_change::Float64, three_month_change::Float64,
              six_month_change::Float64, week_high::Float64, week_low::Float64, week_change::Float64,
              avg_week_change::Float64, avg_week_volume::Float64, week_volume::Float64,
              month_high::Float64, month_low::Float64, month_change::Float64, avg_month_change::Float64,
              avg_month_volume::Float64, month_volume::Float64, year_high::Float64, year_low::Float64,
              year_change::Float64, avg_year_change::Float64, avg_year_volume::Float64,
              ytd_change::Float64, rsi_14::Float64, sto_9::Float64, wpr_14::Float64, mtm_14::Float64,
              roc_14::Float64, ptc::Float64, sar::Float64, volatility::Float64, liquidity::Float64) =
        new(ticker_code, name, description, date_time, previous, change, ma_1, ma_2, ma_5,
            ma_20, ma_50, ma_100, ma_200, ma_percent, ma_return, volume_change,
            three_month_change, six_month_change, week_high, week_low, week_change,
            avg_week_change, avg_week_volume, week_volume, month_high, month_low,
            month_change, avg_month_change, avg_month_volume, month_volume, year_high,
            year_low, year_change, avg_year_change, avg_year_volume, ytd_change, rsi_14,
            sto_9, wpr_14, mtm_14, roc_14, ptc, sar, volatility, liquidity)

    # Construct from XML ::ETree
    function Technical(tl::ETree)
        ticker_code::String = strip(get(tl.attr,"Symbol",""))
        name::String = strip(get(tl.attr,"Name",""))
        description::String = strip(get(tl.attr,"Description",""))
        date_time::DateTime = DateTime(strip(get(tl.attr,"DateTime",""))[1:19], DATETIMEFORMAT_SS)
        previous::Float64 = float(strip(get(tl.attr,"Previous","")))
        change::Float64 = float(strip(get(tl.attr,"Change","")))
        ma_1::Float64 = float(strip(get(tl.attr,"MA1","")))
        ma_2::Float64 = float(strip(get(tl.attr,"MA2","")))
        ma_5::Float64 = float(strip(get(tl.attr,"MA5","")))
        ma_20::Float64 = float(strip(get(tl.attr,"MA20","")))
        ma_50::Float64 = float(strip(get(tl.attr,"MA50","")))
        ma_100::Float64 = float(strip(get(tl.attr,"MA100","")))
        ma_200::Float64 = float(strip(get(tl.attr,"MA200","")))
        ma_percent::Float64 = float(strip(get(tl.attr,"MAPercent","")))
        ma_return::Float64 = float(strip(get(tl.attr,"MAReturn","")))
        volume_change::Float64 = float(strip(get(tl.attr,"VolumeChange","")))
        three_month_change::Float64 = float(strip(get(tl.attr,"ThreeMonthChange","")))
        six_month_change::Float64 = float(strip(get(tl.attr,"SixMonthChange","")))
        week_high::Float64 = float(strip(get(tl.attr,"WeekHigh","")))
        week_low::Float64 = float(strip(get(tl.attr,"WeekLow","")))
        week_change::Float64 = float(strip(get(tl.attr,"WeekChange","")))
        avg_week_change::Float64 = float(strip(get(tl.attr,"AvgWeekChange","")))
        avg_week_volume::Float64 = float(strip(get(tl.attr,"AvgWeekVolume","")))
        week_volume::Float64 = float(strip(get(tl.attr,"WeekVolume","")))
        month_high::Float64 = float(strip(get(tl.attr,"MonthHigh","")))
        month_low::Float64 = float(strip(get(tl.attr,"MonthLow","")))
        month_change::Float64 = float(strip(get(tl.attr,"MonthChange","")))
        avg_month_change::Float64 = float(strip(get(tl.attr,"AvgMonthChange","")))
        avg_month_volume::Float64 = float(strip(get(tl.attr,"AvgMonthVolume","")))
        month_volume::Float64 = float(strip(get(tl.attr,"MonthVolume","")))
        year_high::Float64 = float(strip(get(tl.attr,"YearHigh","")))
        year_low::Float64 = float(strip(get(tl.attr,"YearLow","")))
        year_change::Float64 = float(strip(get(tl.attr,"YearChange","")))
        avg_year_change::Float64 = float(strip(get(tl.attr,"AvgYearChange","")))
        avg_year_volume::Float64 = float(strip(get(tl.attr,"AvgYearVolume","")))
        ytd_change::Float64 = float(strip(get(tl.attr,"YTDChange","")))
        rsi_14::Float64 = float(strip(get(tl.attr,"RSI14","")))
        sto_9::Float64 = float(strip(get(tl.attr,"STO9","")))
        wpr_14::Float64 = float(strip(get(tl.attr,"WPR14","")))
        mtm_14::Float64 = float(strip(get(tl.attr,"MTM14","")))
        roc_14::Float64 = float(strip(get(tl.attr,"ROC14","")))
        ptc::Float64 = float(strip(get(tl.attr,"PTC","")))
        sar::Float64 = float(strip(get(tl.attr,"SAR","")))
        volatility::Float64 = float(strip(get(tl.attr,"Volatility","")))
        liquidity::Float64 = float(strip(get(tl.attr,"Liquidity","")))

        return new(ticker_code, name, description, date_time, previous, change, ma_1, ma_2,
                   ma_5, ma_20, ma_50, ma_100, ma_200, ma_percent, ma_return, volume_change,
                   three_month_change, six_month_change, week_high, week_low, week_change,
                   avg_week_change, avg_week_volume, week_volume, month_high, month_low,
                   month_change, avg_month_change, avg_month_volume, month_volume, year_high,
                   year_low, year_change, avg_year_change, avg_year_volume, ytd_change,
                   rsi_14, sto_9, wpr_14, mtm_14, roc_14, ptc, sar, volatility, liquidity)
    end
end

# ::Ticker
type Ticker
    code::String
    name::String
    long_name::String
    date_time::DateTime

    # Default constructor
    Ticker(code::String, name::String, long_name::String, date_time::DateTime) =
        new(code, name, long_name, date_time)

    # Construct from XML ::ETree
    function Ticker(tk::ETree)
        code::String = strip(get(tk.attr,"Code",""))
        name::String = strip(get(tk.attr,"Name",""))
        long_name::String = strip(get(tk.attr,"LongName",""))
        date_time::DateTime = DateTime(strip(get(tk.attr,"DateTime",""))[1:19], DATETIMEFORMAT_SS)

        return new(code, name, long_name, date_time)
    end
end

# ::Ticker_2
type Ticker_2
    code::String
    name::String

    # Default constructor
    Ticker_2(code::String, name::String) = new(code, name)

    # Construct from XML ::ETree
    function Ticker_2(tk_2::ETree)
        code::String = strip(get(tk_2.attr,"Code",""))
        name::String = strip(get(tk_2.attr,"Name",""))

        return new(code, name)
    end
end

# ::TickerChange
type TickerChange
    old_exchange_code::String
    new_exchange_code::String
    old_ticker_code::String
    new_ticker_code::String
    date_time::DateTime
    is_change_of_exchange_code::Bool
    is_change_of_ticker_code::Bool

    # Default constructor
    TickerChange(old_exchange_code::String, new_exchange_code::String, old_ticker_code::String,
                 new_ticker_code::String, date_time::DateTime, is_change_of_exchange_code::Bool,
                 is_change_of_ticker_code::Bool) =
        new(old_exchange_code, new_exchange_code, old_ticker_code, new_ticker_code, date_time,
            is_change_of_exchange_code, is_change_of_ticker_code)

    # Construct from XML ::ETree
    function TickerChange(tc::ETree)
        old_exchange_code::String = strip(get(tc.attr,"ExchangeCode",""))
        new_exchange_code::String = strip(get(tc.attr,"NewExchangeCode",""))
        old_ticker_code::String = strip(get(tc.attr,"OldSymbol",""))
        new_ticker_code::String = strip(get(tc.attr,"NewSymbol",""))
        date_time::DateTime = DateTime(strip(get(tc.attr,"DateTime",""))[1:19], DATETIMEFORMAT_SS)
        is_change_of_exchange_code::Bool = old_exchange_code != new_exchange_code ? true : false
        is_change_of_ticker_code::Bool = old_ticker_code != new_ticker_code ? true : false

        return    new(old_exchange_code, new_exchange_code, old_ticker_code, new_ticker_code,
                    date_time, is_change_of_exchange_code, is_change_of_ticker_code)
    end
end

# =====
# Files
include("eod_utils_external.jl")
include("eod_utils_internal.jl")
include("eod_ws.jl")

end # module


text:
# Translating
# https://proplot.readthedocs.io/en/latest/subplots.html#Automatic-figure-sizing
# to Julia

# Import proplot using PyCall
using PyCall
proplot = pyimport("proplot")

# Cartopy projections
f, axs = proplot.subplots(ncols=2, nrows=3, proj="robin")
map(x -> x.format(
    land=true, landcolor="k",
    suptitle="Auto figure sizing with grid of cartopy projections"
   ), axs)
f.savefig("figs/Automatic-figure-sizing-1a.svg", transparent=false)

# Images
f, axs = proplot.subplots(ncols=2, nrows=3)
colors = cumsum(rand(10, 20, 3), dims=3)
colors /= maximum(colors)
[ax.imshow(colors) for ax in axs]
[ax.format(
    suptitle="Auto figure sizing with grid of images"
) for ax in axs]
f.savefig("figs/Automatic-figure-sizing-1b.svg", transparent=false)

max_stars_repo_name: natgeo-wong/JuliaProPlotExamples

max_issues_repo_name: natgeo-wong/JuliaProPlotExamples

max_forks_repo_name: natgeo-wong/JuliaProPlotExamples


text:
export DimRedux, GaussDimRedux, SparseDimRedux, SSRFTDimRedux, *

### Dimension reduction operators

"""DimRedux is an abstract dimension reduction map
over the field F, which is either the real or complex numbers"""
abstract type DimRedux{F<:Number} end

"""By default, a DimRedux map reduces a matrix of data via left or right multiplication by its internal variable Ξ.
We extend the multiplication method * to allow multiplication by dimension reduction maps"""
*{F<:Number}(dr::DimRedux{F}, M::Matrix{F}) = dr.Ξ * M
*{F<:Number}(M::Matrix{F}, dr::DimRedux{F}) = M * dr.Ξ' #' nb: different from paper

"""GaussDimRedux is a kind of dimension reduction map
It reduces data over the field F using the Gaussian random matrix Ξ"""
type GaussDimRedux{F<:Number}<:DimRedux{F}
    Ξ::Matrix{F}
end

"""Constructs a GaussDimRedux map which reduces data in F^n to F^k
By default, the field F is a real (floating point) number"""
function GaussDimRedux{F<:Number}(k::Int, n::Int, ::Type{F}=Float64)
    if F<:Real
        Ξ = randn(k,n)
    elseif F<:Complex
        Ξ = randn(k,n) + im*randn(k,n)
    else error("GaussDimRedux not defined over the field $F")
    end
    return GaussDimRedux(Ξ)
end

"""SparseDimRedux is a kind of dimension reduction map
It reduces data over the field F using the sparse random matrix Ξ"""
type SparseDimRedux{F<:Number}<:DimRedux{F}
    Ξ::SparseMatrixCSC{F}
end

"""Constructs a SparseDimRedux map which reduces data in F^n to F^k
By default, the field F is the real numbers
p is the sparsity level, with default level .1 (???)"""
function SparseDimRedux{F<:Number}(k::Int, n::Int, p::Float64=.1, ::Type{F}=Real)
    if F<:Real
        Ξ = sprandn(k,n,p) # nb --- not defined in paper
    elseif F<:Complex
        Ξ = sprandn(k,n,p) + im*sprandn(k,n,p) # nb --- not defined in paper
    else error("SparseDimRedux not defined over the field $F")
    end
    return SparseDimRedux(Ξ)
end

"""SSRFTDimRedux is a kind of dimension reduction map
It reduces data over the field F using a scrambled SRFT
Its internal variables are
the permutations Π1 and Π2, the random signs ϵ1 and ϵ2, and the coordinate restriction R"""
type SSRFTDimRedux{F<:Number}<:DimRedux{F}
    Π1::Array{Int,1} # first permutation
    ϵ1::Array{F,1} # first random signs
    Π2::Array{Int,1} # second permutation
    ϵ2::Array{F,1} # second random signs
    R::Array{Int,1} # restriction onto k coordinates
end

"""Constructs a SSRFTDimRedux map which reduces data in F^n to F^k
By default, the field F is the real numbers"""
function SSRFTDimRedux{F<:Number}(k::Int, n::Int, ::Type{F}=Float64)
    Π1 = randperm(n) # a random permutation of {1,...,n}
    Π2 = randperm(n) # a random permutation of {1,...,n}
    R = randperm(n)[1:k] # a random selection of k elements from {1,...,n}
    if F<:Real
        ϵ1 = one(F) * rand([-1, 1], n) # a random ±1 vector of length n in the field F
        ϵ2 = one(F) * rand([-1, 1], n) # a random ±1 vector of length n in the field F
    elseif F<:Complex
        ϵ1 = one(F) * rand([-1, 1], n) + im*one(F) * rand([-1, 1], n) # a random ±1±i vector of length n in the field F
        ϵ2 = one(F) * rand([-1, 1], n) + im*one(F) * rand([-1, 1], n) # a random ±1±i vector of length n in the field F
    else error("SSRFTDimRedux not defined over the field $F")
    end
    return SSRFTDimRedux(Π1,ϵ1,Π2,ϵ2,R)
end

"""The SSRFTDimRedux map uses an efficient specialized method to reduce a matrix of data,
multiplying it on the left by Ξ = R FFT ϵ2 Π2 FFT ϵ1 Π1"""
function *{F<:Number}(dr::SSRFTDimRedux{F}, M::Matrix{F})
    B = M[dr.Π1,:] # scramble rows
    B .= dr.ϵ1 .* B # scramble signs - equivalent to diagm(dr.ϵ1) * B
    if F<:Real
        B = dct(B)
    elseif F<:Complex
        B = fft(B)
    end
    B = B[dr.Π2,:] # scramble rows
    B .= dr.ϵ2 .* B # scramble signs - equivalent to diagm(dr.ϵ2) * B
    if F<:Real
        B = dct(B)
    elseif F<:Complex
        B = fft(B)
    end
    B = B[dr.R,:] # select rows
    return B
end
"""The SSRFTDimRedux map uses an efficient specialized method to reduce a matrix of data,
multiplying it on the right by Ξ' = (R FFT ϵ2 Π2 FFT ϵ1 Π1)'"""
function *{F<:Number}(M::Matrix{F}, dr::SSRFTDimRedux{F})
  B = M[:,dr.Π1] # scramble columns
  B .= B .* dr.ϵ1' # scramble signs - equivalent to B * diagm(dr.ϵ1)
  if F<:Real
      B = dct(B)
  elseif F<:Complex
      B = ifft(B)
  end
  B = M[:,dr.Π2] # scramble columns
  B .= B .* dr.ϵ2' # scramble signs - equivalent to B * diagm(dr.ϵ1)
  if F<:Real
      B = dct(B)
  elseif F<:Complex
      B = ifft(B)
  end
  B = B[:,dr.R] # select columns
  return B
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Functions related with the command line.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

"""
    _print_cmd_message!(pagerd::Pager, msg::String)

Print the message `msg` in the command line of pager `pagerd`. The string
formatting can be selected using the keyword `crayon`.

"""
function _print_cmd_message!(pagerd::Pager, msg::String; crayon::Crayon = crayon())
    term = pagerd.term
    display_size = pagerd.display_size

    if get(term.out_stream, :color, true)
        _d = string(Crayon(reset = true))
        _h = string(crayon)
    else
        _d = ""
        _h = ""
    end

    # Move the cursor to the last line and print the message.
    _move_cursor(term.out_stream, display_size[1], 0)
    write(term.out_stream, _h)
    write(term.out_stream, msg)
    write(term.out_stream, _d)
    _clear_to_eol(term.out_stream)

    return nothing
end

"""
    _redraw_cmd_line!(pagerd::Pager)

Print the command line of pager `pagerd` to the display.

"""
function _redraw_cmd_line!(pagerd::Pager)
    # Unpack variables.
    term          = pagerd.term
    display_size  = pagerd.display_size
    num_lines     = pagerd.num_lines
    lines_cropped = pagerd.lines_cropped
    mode          = pagerd.mode
    features      = pagerd.features

    if get(term.out_stream, :color, true)::Bool
        _d = string(Crayon(reset = true))
        _g = string(crayon"dark_gray")
    else
        _d = ""
        _g = ""
    end

    # Compute the information considering the current mode.
    if mode == :view
        cmd_help = "(↑ ↓ ← →:move, "
        :help ∈ features && (cmd_help *= "?:help, ")
        cmd_help *= "q:quit)"
    elseif mode == :searching
        active_search_match_id = pagerd.active_search_match_id
        search_matches         = pagerd.search_matches

        num_matches = length(search_matches)

        # Check if there are matches.
        if num_matches > 0
            cmd_help = "(match " * string(active_search_match_id) *
                " of " * string(num_matches) * ")"
        else
            cmd_help = "(no match found)"
        end
    else
        cmd_help = "ERROR"
    end

    # Compute the scroll position
    pos = lpad(string(round(Int, 100 * (1 - lines_cropped / num_lines))), 3)
    cmd_help *= " " * pos * "%"

    lcmd_help = length(cmd_help)

    if display_size[2] > (lcmd_help + 4)
        cmd_aligned = " "^(display_size[2] - lcmd_help - 1) * _g * cmd_help * _d
    else
        cmd_aligned = ""
    end

    # Move the cursor to the last line and print the command line.
    _move_cursor(term.out_stream, display_size[1], 0)
    write(term.out_stream, ":" * cmd_aligned)
    _move_cursor(term.out_stream, display_size[1], 2)

    return nothing
end

"""
    _read_cmd!(pagerd::Pager; prefix::String = "/")

Read a command in the pager `pagerd`.

This function returns a string with the command.

"""
function _read_cmd!(pagerd::Pager; prefix::String = "/")
    # Unpack values.
    term         = pagerd.term
    display_size = pagerd.display_size

    # Initialize variables.
    cmd = ""
    cmd_width = 0
    cursor_pos = 1
    prefix_size = textwidth(prefix)
    redraw = true

    while true
        if redraw
            # Clear command line.
            _move_cursor(term.out_stream, display_size[1], 1)
            _clear_to_eol(term.out_stream)
            write(term.out_stream, prefix * cmd)

            # Restore the cursor position
            _move_cursor(
                term.out_stream,
                display_size[1],
                cursor_pos + prefix_size
            )

            redraw = false
        end

        k = _jlgetch(term.in_stream)

        if k.value == :enter
            break

        elseif k.value isa String
            cmd = first(cmd, (cursor_pos - 1)) *
                k.value *
                last(cmd, cmd_width - (cursor_pos - 1))
            cmd_width += 1
            cursor_pos += 1
            redraw = true
        elseif k.value == :backspace
            if cmd_width > 0
                cmd = first(cmd, cmd_width - 1)
                cmd_width -= 1
                cursor_pos -= 1
                redraw = true
            else
                break
            end
        elseif k.value == :left
            if cursor_pos > 1
                cursor_pos -= 1
                _cursor_back(term.out_stream)
            end
        elseif k.value == :right
            if cursor_pos < cmd_width + 1
                cursor_pos += 1
                _cursor_forward(term.out_stream)
            end
        elseif k.value == :home
            cursor_pos = 1
            _move_cursor(
                term.out_stream,
                display_size[1],
                cursor_pos + prefix_size
            )
        elseif k.value == :end
            cursor_pos = cmd_width + 1
            _move_cursor(
                term.out_stream,
                display_size[1],
                cursor_pos + prefix_size
            )
        end
    end

    return cmd
end


text:
using Pkg.Artifacts, LibGit2, ghr_jll

function generate_artifact(
    datafiles;
    data_dir="dev/artifacts/data",
    root=".",
    artifact_toml=joinpath(root, "Artifacts.toml"),
    deploy=true,
    tag="data"
)

    if deploy && !haskey(ENV, "GITHUB_TOKEN")
        @warn "For automatic github deployment, need GITHUB_TOKEN. Not found in ENV, attemptimg global git config."
    end

    if deploy
        # Where we will put our tarballs
        tempdir = mktempdir()

        function get_git_remote_url(repo_path::String)
            repo = LibGit2.GitRepo(repo_path)
            origin = LibGit2.get(LibGit2.GitRemote, repo, "origin")
            return LibGit2.url(origin)
        end

        # Try to detect where we should upload these weights to (or just override
        # as shown in the commented-out line)
        origin_url = get_git_remote_url(root)
        deploy_repo = "$(basename(dirname(origin_url)))/$(splitext(basename(origin_url))[1])"

    end

    # For each BSON file, generate its own artifact:
    for datafile in datafiles

        # Name for hash/artifact:
        name = splitext(datafile)[1]

        # create_artifact() returns the content-hash of the artifact directory once we're finished creating it
        hash = create_artifact() do artifact_dir
            cp(joinpath(data_dir, datafile), joinpath(artifact_dir, datafile))
        end

        # Spit tarballs to be hosted out to local temporary directory:
        if deploy

            tarball_hash = archive_artifact(hash, joinpath(tempdir, "$(name).tar.gz"))

            # Calculate tarball url
            tarball_url = "https://github.com/$(deploy_repo)/releases/download/$(tag)/$(name).tar.gz"

            # Bind this to an Artifacts.toml file
            @info("Binding $(name) in Artifacts.toml...")
            bind_artifact!(
                artifact_toml, name, hash;
                download_info=[(tarball_url, tarball_hash)], lazy=true, force=true
            )
        end

    end

    if deploy
        # Upload tarballs to a special github release
        @info("Uploading tarballs to $(deploy_repo) tag `$(tag)`")

        ghr() do ghr_exe
            println(readchomp(`$ghr_exe -replace -u $(dirname(deploy_repo)) -r $(basename(deploy_repo)) $(tag) $(tempdir)`))
        end

        @info("Artifacts.toml file now contains all bound artifact names")
    end

end

# generate_artifact()


text:
import Franklin, LiveServer
# Swap out the highlight.js library for prerendering
Franklin.HIGHLIGHTJS[] = get(ENV, "HIGHLIGHTJS", "highlight.js")

function on_write(pg, vars; output_path=nothing)
    output_path === nothing && return
    if endswith(output_path, "highlight-julia/index.html")
        l = length(pg)
        r = r"""^

<code"""m
        s = """

 s)
        l′ = length(pg)
        r = r"""

"""m
        s = """

"""
        pg = replace(pg, r => s)
        l′′ = length(pg)
        @assert div(l - l′, 5) == div(l′ - l′′, 4)
        @assert l - l′ > 0
        write(output_path, pg)
    end
    return nothing
end