text:
using Parameters

@with_kw mutable struct NSGAII <: Algorithm
    problem::Problem
    eval_fn::Function
    population_size::Int64
    n_iters::Int64
    archive::Archive
    archive_frequency::Int64
    maxtime::Float64            # in seconds
end

"""
    multistart(algo, nstarts; seedpop)

Runs the algorithm nstarts times and returns the pareto pool from each of the
runs.
"""
function multistart(algo, nstarts; seedpop::Vector{Solution}=Vector{Solution}())
    results = Vector{Vector{Solution}}()
    for i in 1:nstarts
        push!(results, garun(algo, seedpop=seedpop))
    end
    return results
end

"Convenience constructor with no archiving, kwarg for setting maxtime"
function NSGAII(problem::Problem, eval_fn::Function, pop_size::Int64, n_iters::Int64, maxtime=0.0)
    archive = Archive(compare_pareto_dominance, Vector{Solution}())
    archive_frequency = 0
    return NSGAII(problem, eval_fn, pop_size, n_iters, archive, archive_frequency, maxtime)
end

"Run the NSGAII algorithm"
function garun(algo::NSGAII;
               starting_gen=1,  # set a different value for restarts
               seedpop::Vector{Solution}=Vector{Solution}(),
               logging_frequency=0,
               logging_destination=Nothing)
    population = init_pop(algo; seedpop=seedpop)

    start_time = time()

    for gen in starting_gen:algo.n_iters
        population = iter_generation(algo, population, gen)
        if (algo.archive_frequency > 0) && (gen % algo.archive_frequency == 0)
            insert_solutions!(algo.archive, population)
        end
        if (logging_frequency > 0) && (gen % logging_frequency == 0)
            log_population(population, gen, logging_destination)
        end
        if (algo.maxtime > 0) && (time() - start_time > algo.maxtime)
            break
        end
    end

    return population
end

function init_pop(algo::NSGAII; seedpop::Vector{Solution}=Vector{Solution}())
    N = algo.population_size
    population = vcat(seedpop,
                      [random_candidate(algo.problem) for i in 1:(N-length(seedpop))])

    for p in population
        evaluate!(p)  # distribute
    end

    nondominated_sort(population)
    return population
end

"""
    iter_generation(algo, population, gen)

Perform one generation of the GA. Returns the new child population.
"""
function iter_generation(algo::NSGAII, population::Vector{Solution}, gen::Int64)
    # create N new indivs: binary tournament -> SBX/PM
    N = algo.population_size
    for i = 1:(N/2)
        p1 = tournament_selector(population, 2, dominance=nondominated_cmp)
        p2 = tournament_selector(population, 2, dominance=nondominated_cmp)
        c1, c2 = [PM(c) for c in SBX(p1, p2)]
        for c in [c1, c2]
            evaluate!(c)
            c.generation = gen
        end
        push!(population, c1, c2)
    end

    # assign ranks and crowding_distance
    nondominated_sort(population)

    # select best N as next generation
    population = nondominated_truncate(population, N)

    return population
end

# TODO: add Archive to NSGAIII
@with_kw mutable struct NSGAIII <: Algorithm
    problem::Problem
    eval_fn::Function
    nobjs::Int64
    ndivs::Int64
    population_size::Int64
    n_iters::Int64
    reference_points::Array{Float64, 2}
end

function NSGAIII(problem::Problem, evalfn::Function, ndivs::Int, niters::Int)
    refpts = generate_regular_reference_points(problem.nobjs, ndivs)
    return NSGAIII(
        problem,
        evalfn,
        problem.nobjs,
        ndivs,
        Int(ceil(size(refpts, 1)/4)*4),
        niters,
        refpts
    )
end

function garun(algo::NSGAIII; seedpop::Vector{Solution}=Vector{Solution}())
    println("Echidna.garun")

    population = init_pop(algo; seedpop=seedpop)

    for gen in 1:algo.n_iters
        println("")
        println("generation: ", gen)
        population = iter_generation(algo, population, gen)
    end

    return population
end

function init_pop(algo::NSGAIII; seedpop::Vector{Solution}=Vector{Solution}())
    N = algo.population_size
    population = vcat(seedpop,
                      [random_candidate(algo.problem) for i in 1:(N-length(seedpop))])

    for p in population
        evaluate!(p)
    end

    nondominated_sort(population)

    return population
end

function iter_generation(algo::NSGAIII, population::Vector{Solution}, gen::Int64)
    # TODO: this is copied from NSGAII - make sure there shouldn't be other changes
    # create N new indivs: binary tournament -> SBX/PM
    N = algo.population_size
    for i = 1:(N/2)
        p1 = tournament_selector(population, 2, dominance=nondominated_cmp)
        p2 = tournament_selector(population, 2, dominance=nondominated_cmp)
        c1, c2 = [PM(c) for c in SBX(p1, p2)]
        for c in [c1, c2]
            evaluate!(c)
            c.generation = gen
        end
        push!(population, c1, c2)
    end

    # assign ranks and crowding_distance
    nondominated_sort(population)

    # select best N as next generation
    population = reference_point_truncate(population, N, algo.reference_points)

    return population
end

max_stars_repo_name: phawthorne/echidna

max_issues_repo_name: phawthorne/echidna

max_forks_repo_name: phawthorne/echidna


text:
# INSIDE MAIN

#=
Main talks to user
=#

#=
function 1:
get data from user
=#

#=
function 2:
send data to user
=#

#=
function 3:
update data
=#

# INSIDE MASTER

#=
function 1:
get data from user
=#

#=
function 2:
send data from user
=#

# INSIDE DATA

#=
function 1:
get data from user
=#

#=
function 2:
send data from user
=#

ext: jl

max_stars_repo_name: cidazc/caDDS

max_issues_repo_name: cidazc/caDDS

max_forks_repo_name: cidazc/caDDS


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule libLLVM_assert_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("libLLVM_assert")
JLLWrappers.@generate_main_file("libLLVM_assert", UUID("393490a8-41a6-5dea-8ab1-81a396b80891"))
end  # module libLLVM_assert_jll


text:
###########################################################################
#
#  Serial Tempering Monte-Carlo
#
#  Takes MCMCTasks and performs Serial Tempering on a set of MCMCTasks
#
#
#  Ref : "Bayes Factors via Serial Tempering" C.Geyer
#
###########################################################################

export SerialTempMC

println("Loading SerialTempMC(steps, burnin, swapPeriod) runner")

immutable SerialTempMC <: MCMCRunner
  steps::Int
  burnin::Int
  swapPeriod::Int

  function SerialTempMC(steps::Int, burnin::Int, swapPeriod::Int)
	  @assert burnin >= 0 "Burnin rounds ($burnin) should be >= 0"
	  @assert steps > burnin "Steps ($steps) should be > to burnin ($burnin)"

    new(steps, burnin, swapPeriod)
  end
end

SerialTempMC(; steps::Int=1, burnin::Int=0, swapPeriod::Int=5) = SerialTempMC(steps, burnin, swapPeriod)

function run_serialtempmc(tasks::Array{MCMCTask})

	local nmods = length(tasks)
	local tsize = tasks[end].model.size
  local steps = tasks[end].runner.steps
  local burnin = tasks[end].runner.burnin
  local swapPeriod = tasks[end].runner.swapPeriod

	@assert all(map(t->t.model.size, tasks) .== tsize) "Models do not have the same parameter vector size"

	tic() # start timer

	# initialize all tasks
	map(t -> consume(t.task), tasks)

	# initialize MCMCChain result (one Chain only)
  res = MCMCChain(1:1:2, DataFrame(tsize, steps-burnin), tasks[end])

	local logW = zeros(nmods)  # log of task weights that will be adapted
	local at = 1  # pick starting task
	local s = consume(tasks[at].task)
	local ppars, logtarget, pars
	ppars, logtarget, pars = s.ppars, s.logtarget, s.pars

	for i in 1:steps  # i = 1

		if i % swapPeriod == 0  # attempt a task switch
			# pick another task at random
			at2 = rand(1:(nmods-1))
			at2 = at2 >= at ? at2+1 : at2

			reset(tasks[at2], pars)
			s2 = consume(tasks[at2].task)
			if rand() < exp(logtarget - s2.logtarget + logW[at2] - logW[at]) # accept swap ?
				at, s = at2, s2
			end
		else
			s = consume(tasks[at].task)
		end

		# TODO : add logW adaptation

		ppars, logtarget, pars = s.ppars, s.logtarget, s.pars

		if i > burnin # store ppars and the model we're on
			pos = i-burnin
			res.samples[:, pos] = ppars
			#res.misc[:mod][pos] = at
		end
	end

	res.runTime = toq()
	res
end

# TODO: check that all elements of array contain MCMCTasks of the same type
function resume_serialtempmc(tasks::Array{MCMCTask}; steps::Int=100)
	run(MCMCTask[tasks[i].model * tasks[i].sampler * SerialTempMC(steps=steps, swapPeriod=tasks[i].runner.swapPeriod)
		for i in 1:length(tasks)])
end


text:
function find_stdlib_deps(ctx::Context, path::String)
    stdlib_deps = Dict{UUID, String}()
    regexps = [Regex("\\b(import|using)\\s+((\\w|\\.)+\\s*,\\s*)*$lib\\b") for lib in values(ctx.stdlibs)]
    for (root, dirs, files) in walkdir(path; onerror = x->nothing)
        for file in files
            endswith(file, ".jl") || continue
            filecontent = try read(joinpath(root, file), String)
                catch e
                    e isa SystemError || rethrow()
                    ""
                end
            for ((uuid, stdlib), r) in zip(ctx.stdlibs, regexps)
                if occursin(r, filecontent)
                    stdlib_deps[uuid] = stdlib
                end
            end
        end
    end
    return stdlib_deps
end

# Backwards compatibility with Pkg2 REQUIRE format
function collect_require!(ctx::Context, pkg::PackageSpec, path::String,
                          fix_deps_map::Dict{UUID,Vector{PackageSpec}})
    fix_deps = PackageSpec[]
    reqfile = joinpath(path, "REQUIRE")
    # Checked out "old-school" packages have by definition a version higher than all registered.
    set_maximum_version_registry!(ctx, pkg)
    !haskey(fix_deps_map, pkg.uuid) && (fix_deps_map[pkg.uuid] = valtype(fix_deps_map)())
    if isfile(reqfile)
        for r in Pkg2.Reqs.read(reqfile)
            r isa Pkg2.Reqs.Requirement || continue
            pkg_name, vspec = r.package, VersionSpec(VersionRange[r.versions.intervals...])
            if pkg_name == "julia"
                if !(VERSION in vspec)
                    @warn("julia version requirement for package $(pkg.name) not satisfied")
                end
            else
                deppkg = PackageSpec(pkg_name, vspec)
                push!(fix_deps_map[pkg.uuid], deppkg)
                push!(fix_deps, deppkg)
            end
        end

        # Packages from REQUIRE files need to get their UUID from the registry
        registry_resolve!(ctx, fix_deps)
        project_deps_resolve!(ctx, fix_deps)
        ensure_resolved(ctx, fix_deps; registry=true)
    end

    # And collect the stdlibs
    stdlibs = find_stdlib_deps(ctx, path)
    for (uuid, name) in stdlibs
        deppkg = PackageSpec(name, uuid)
        push!(fix_deps_map[pkg.uuid], deppkg)
        push!(fix_deps, deppkg)
    end

    return
end

# Pkg2 test/REQUIRE compatibility
function pkg2_test_target_compatibility!(ctx, path, pkgs)
    test_reqfile = joinpath(path, "test", "REQUIRE")
    if isfile(test_reqfile)
        for r in Pkg2.Reqs.read(test_reqfile)
            r isa Pkg2.Reqs.Requirement || continue
            pkg_name, vspec = r.package, VersionSpec(VersionRange[r.versions.intervals...])
            pkg_name == "julia" && continue
            push!(pkgs, PackageSpec(pkg_name, vspec))
        end
        registry_resolve!(ctx, pkgs)
        project_deps_resolve!(ctx, pkgs)
        ensure_resolved(ctx, pkgs; registry=true)
    end
    return nothing
end


text:
# Make sure PERK.krr computes Eq. 11 in Nataraj et al. for exact kernels
function test_krr_1()

    Random.seed!(0)
    T = 100
    xtrain = LinRange(0, 1, T)
    ytrain = xtrain .+ 0.1 .* randn(T)
    kernel = EuclideanKernel()
    ρ = 0.01
    trainData = PERK.krr_train(xtrain, ytrain, kernel, ρ)
    ytest = 0.6
    xhat = PERK.krr(ytest, trainData, kernel)

    K = ytrain * ytrain'
    M = I - ones(T, T) / T
    k = ytest * conj(ytrain) - K * ones(T) / T
    x_correct = xtrain' * (ones(T) / T + M * ((M * K * M + T * ρ * I) \ k))

    # For line coverage, access Q, L, and T in trainData as well
    return xhat ≈ x_correct &&
           trainData.Q == length(ytest) &&
           trainData.L == length(xhat) &&
           trainData.T == 100

end

function test_krr_2()

    Random.seed!(0)
    T = 100
    D = 3
    xtrain = LinRange(0, 1, T)
    ytrain = xtrain' .+ 0.1 .* randn(D, T)
    Λ = ones(D)
    kernel = GaussianKernel(Λ)
    ρ = 0.01
    trainData = PERK.krr_train(xtrain, ytrain, kernel, ρ)
    ytest = transpose([0.6 0.57 0.67])
    xhat = PERK.krr(ytest, trainData, kernel)

    gauss = (p, q) -> exp(-0.5 * norm((1 ./ Λ) .* (p - q))^2)
    K = [gauss(ytrain[:,m], ytrain[:,n]) for m = 1:T, n = 1:T]
    M = I - ones(T, T) / T
    k = [gauss(vec(ytest), ytrain[:,n]) for n = 1:T] - K * ones(T) / T
    x_correct = xtrain' * (ones(T) / T + M * ((M * K * M + T * ρ * I) \ k))

    return xhat[] ≈ x_correct

end

# Make sure PERK.krr with EuclideanKernel produces the same result as ridge
# regression
function test_krr_3()

    Random.seed!(0)
    T = 100
    xtrain = LinRange(0, 1, T)
    ytrain = xtrain .+ 0.1 .* randn(T)
    kernel = EuclideanKernel()
    ρ = 0.01
    trainData = PERK.krr_train(xtrain, ytrain, kernel, ρ)
    N = 10
    ytest = LinRange(0.1, 0.9, N)
    xhat = PERK.krr(ytest, trainData, kernel)

    xtilde = [xtrain .- mean(xtrain); 0]
    Ytilde = [ytrain .- mean(ytrain); sqrt(T * ρ)]
    w = Ytilde \ xtilde
    b = mean(xtrain) - w * mean(ytrain)
    x_correct = w .* ytest .+ b

    return xhat ≈ x_correct

end

function test_krr_4()

    Random.seed!(0)
    T = 100
    D = 3
    xtrain = LinRange(0, 1, T)
    ytrain = xtrain' .+ 0.1 .* randn(D, T)
    kernel = EuclideanKernel()
    ρ = 0.01
    trainData = PERK.krr_train(xtrain, ytrain, kernel, ρ)
    N = 10
    ytest = LinRange(0.1, 0.9, N)' .+ 0.1 .* randn(D, N)
    xhat = PERK.krr(ytest, trainData, kernel)

    xtilde = [xtrain .- mean(xtrain); zeros(D)]
    Ytilde = [(ytrain .- mean(ytrain, dims = 2))'; sqrt(T * ρ) * I]
    w = Ytilde \ xtilde
    b = mean(xtrain) - w' * dropdims(mean(ytrain, dims = 2), dims = 2)
    x_correct = [w' * ytest[:,n] + b for n = 1:N]

    return xhat ≈ x_correct

end

# Other tests
function test_krr_5()

    rng = StableRNG(0)
    f = x -> exp(-30 / x)
    xtrue = 100
    y = f(xtrue)
    T = 200
    xDists = Uniform(10, 500)
    noiseDist = Normal(0, 0.01)
    signalModels = f
    λ = 2.0^-1.5
    H = 100
    kernel = GaussianRFF(H, λ * mean(y))
    ρ = 2.0^-20
    (ytrain, xtrain) = generatenoisydata(rng, T, xDists, noiseDist, signalModels)
    trainData = PERK.krr_train(xtrain, ytrain, kernel, ρ, randn(rng, H),
        rand(rng, H))
    xhat = PERK.krr(y, trainData, kernel)

    # For code coverage, call
    # - krr_train without rng and freq and phase
    # - generatenoisydata without rng
    # - addnoise! without rng
    # - train without rng
    # (results not used)
    PERK.krr_train(xtrain, ytrain, kernel, ρ)
    generatenoisydata(T, xDists, noiseDist, signalModels)
    PERK.addnoise!([y], noiseDist)
    PERK.train(T, xDists, noiseDist, signalModels, kernel, ρ)

    error_rel = abs(xhat - xtrue) / xtrue
    # For line coverage, access Q, L, and H in trainData as well
    return isapprox(error_rel, 0.07054474887124201, atol = 1e-7) &&
           trainData.Q == length(y) &&
           trainData.L == length(xtrue) &&
           trainData.H == H

end

function test_krr_6()

    rng = StableRNG(0)
    f = x -> exp(-30 / x)
    xtrue = 100
    y = f(xtrue)
    T = 200
    xDists = Uniform(10, 500)
    noiseDist = Normal(0, 0.01)
    signalModels = f
    λ = 2.0^-1.5
    H = 100
    kernel = GaussianRFF(H, λ * mean(y))
    ρ = 2.0^-20
    (ytrain, xtrain) = generatenoisydata(rng, T, xDists, noiseDist, signalModels)
    trainData = PERK.krr_train(xtrain, ytrain, kernel, ρ, randn(rng, H, 1),
        rand(rng, H))
    xhat = PERK.krr(y, trainData, kernel)

    error_rel = abs(xhat - xtrue) / xtrue
    return isapprox(error_rel, 0.07054474887124201, atol = 1e-7)

end

@testset "Kernel Ridge Regression" begin

    @test test_krr_1()
    @test test_krr_2()
    @test test_krr_3()
    @test test_krr_4()
    @test test_krr_5()
    @test test_krr_6()

end


text:
"""
    add_sub_plot3!(f)

Internal function to add a subroutine to the GLE script to enable plot3.
"""
function add_sub_plot3!(
            f::Figure
         )::Nothing

    "plot3" ∈ keys(f.subroutines) && return
    xs = "(xo-xmin)/xspan"
    ys = "(yo-ymin)/yspan"
    f.subroutines["plot3"] = """
        sub plot3 data\$ xmin xspan ymin yspan showline showmarker marker\$ mscale
            xo = 0
            yo = 0
            zo = 0
            io = 0
            fopen data\$ file read
            until feof(file)
                fread file x y z
                xo = x
                yo = y
                zo = z
                io = io+1
                if (io<=1) then
                    amove xg3d($xs,$ys,zo) yg3d($xs,$ys,zo)
                    if (showmarker > 0) then
                        marker marker\$ mscale
                    end if
                else
                    if (showline >= 1) then
                        aline xg3d($xs,$ys,zo) yg3d($xs,$ys,zo)
                    else
                        amove xg3d($xs,$ys,zo) yg3d($xs,$ys,zo)
                    end if
                    if (showmarker > 0) then
                        marker marker\$ mscale
                    end if
                end if
            next
            fclose file
        end sub
        """
    return nothing
end


text:
include("models.jl")
include("framework.jl")
include("fit.jl")


text:
module Particles

export Simulation, run!

include("simulation.jl")

end # module

max_stars_repo_name: michakraus/JuliaAtHomeAssignments

max_issues_repo_name: michakraus/JuliaAtHomeAssignments

max_forks_repo_name: michakraus/JuliaAtHomeAssignments


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta

const Normalizer_NormType = (;[
    Symbol("LMax") => Int32(0),
    Symbol("L1") => Int32(1),
    Symbol("L2") => Int32(2),
]...)

mutable struct Normalizer <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function Normalizer(; kwargs...)
        obj = new(meta(Normalizer), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct Normalizer
const __meta_Normalizer = Ref{ProtoMeta}()
function meta(::Type{Normalizer})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Normalizer)
            __meta_Normalizer[] = target = ProtoMeta(Normalizer)
            allflds = Pair{Symbol,Union{Type,String}}[:normType => Int32]
            meta(target, Normalizer, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_Normalizer[]
    end
end
function Base.getproperty(obj::Normalizer, name::Symbol)
    if name === :normType
        return (obj.__protobuf_jl_internal_values[name])::Int32
    else
        getfield(obj, name)
    end
end

export Normalizer_NormType, Normalizer


text:
### reinterpret
#
## Color->T
Base.reinterpret{CV1<:Colorant,CV2<:Colorant}(::Type{CV1}, A::Array{CV2,1}) = _reinterpret_cvarray(CV1, A)
Base.reinterpret{CV1<:Colorant,CV2<:Colorant}(::Type{CV1}, A::Array{CV2})   = _reinterpret_cvarray(CV1, A)
Base.reinterpret{T<:Number,CV<:Colorant}(::Type{T}, A::Array{CV,1}) = _reinterpret_cvarray(T, A)
Base.reinterpret{T<:Number,CV<:Colorant}(::Type{T}, A::Array{CV})   = _reinterpret_cvarray(T, A)

function _reinterpret_cvarray{CV1<:Colorant,CV2}(::Type{CV1}, A::Array{CV2})
    reinterpret(CV1, A, size(A))
end
@inline function _reinterpret_cvarray{T<:Number,CV<:Colorant}(::Type{T}, A::Array{CV})
    reinterpret(T, A, (_len(CV), size(A)...))
end
if squeeze1
    function _reinterpret_cvarray{T<:Number,CV<:Color1}(::Type{T}, A::Array{CV})
        reinterpret(T, A, size(A))
    end
end
_len{C}(::Type{C}) = _len(C, eltype(C))
_len{C}(::Type{C}, ::Type{Any}) = error("indeterminate type")
_len{C,T}(::Type{C}, ::Type{T}) = sizeof(C) ÷ sizeof(T)

## T->Color
# We have to distinguish two forms of call:
#   form 1: reinterpret(RGB, img)
#   form 2: reinterpret(RGB{UFixed8}, img)
Base.reinterpret{T<:Number,CV<:Colorant}(::Type{CV}, A::Array{T,1}) = _reinterpret(CV, eltype(CV), A)
Base.reinterpret{T<:Number,CV<:Colorant}(::Type{CV}, A::Array{T})   = _reinterpret(CV, eltype(CV), A)

function _reinterpret{T,CV<:Colorant}(::Type{CV}, ::Type{Any}, A::Array{T})
    # form 1 (turn into a form 2 call by filling in the element type of the array)
    _reinterpret_array_cv(CV{T}, A)
end
function _reinterpret{T,CV<:Colorant}(::Type{CV}, TT::Type, A::Array{T})
    # form 2
    _reinterpret_array_cv(CV, A)
end

function _reinterpret_array_cv{T,CV<:Colorant}(::Type{CV}, A::Array{T})
    reinterpret(CV, A, tail(size(A)))
end
if squeeze1
    function _reinterpret_array_cv{T,CV<:Color1}(::Type{CV}, A::Array{T})
        reinterpret(CV, A, size(A))
    end
end

# This version is used by the deserializer to convert UInt8 buffers
# back to their original type. Fixes Images#287.
# _reinterpret_array_cv{CV<:Colorant}(::Type{CV}, A::Vector{UInt8}) =
#     reinterpret(CV, A, (div(length(A), sizeof(CV)),))

### convert
#
# The main contribution here is "concretizing" the colorant type: allow
#    convert(RGB, a)
# rather than requiring
#    convert(RGB{U8}, a)
# Where possible the raw element type of the source is retained.
Base.convert{C<:Colorant,n}(::Type{Array{C}},   img::Array{C,n}) = img
Base.convert{C<:Colorant,n}(::Type{Array{C,n}}, img::Array{C,n}) = img

function Base.convert{Cdest<:Colorant,n,Csrc<:Colorant}(::Type{Array{Cdest}},
                                                        img::AbstractArray{Csrc,n})
    convert(Array{Cdest,n}, img)
end

function Base.convert{Cdest<:Colorant,n,Csrc<:Colorant}(::Type{Array{Cdest,n}},
                                                        img::AbstractArray{Csrc,n})
    copy!(Array{ccolor(Cdest, Csrc)}(size(img)), img)
end

# float32, float64, etc. Used for conversions like
#     Array{RGB{U8}} -> Array{RGB{Float32}},
# since
#    convert(Array{RGB{Float32}}, A)
# is annoyingly verbose for such a common operation.
for (fn,T) in (#(:float16, Float16),   # Float16 currently has promotion problems
               (:float32, Float32), (:float64, Float64),
               (:ufixed8, UFixed8), (:ufixed10, UFixed10),
               (:ufixed12, UFixed12), (:ufixed14, UFixed14), (:ufixed16, UFixed16))
    @eval begin
        ($fn){C<:Colorant}(::Type{C}) = base_colorant_type(C){$T}
        ($fn){S<:Number  }(::Type{S}) = $T
        ($fn)(c::Colorant) = convert(($fn)(typeof(c)), c)
        ($fn)(n::Number)   = convert(($fn)(typeof(n)), n)
        @deprecate ($fn){C<:Colorant}(A::AbstractArray{C}) ($fn).(A)
        fname = $(Expr(:quote, fn))
        Tname = $(Expr(:quote, T))
@doc """
    $fname.(img)

converts the raw storage type of `img` to `$Tname`, without changing the color space.
""" $fn

    end
end
const u8 = ufixed8
const u16 = ufixed16


text:
# Similar to ExtrapDimSpec but for only a single dimension
const ExtrapSpec = Union{BoundaryCondition,Tuple{BoundaryCondition,BoundaryCondition}}

# Union over all BoundaryCondition that yield an infinite sequence of knots
# Note: Extrapolation must result in extrapolated knots, simply evaluating
# outside of the interpolation's bounds is not sufficient
const RepeatKnots = Union{Periodic,Reflect}

"""
    KnotIterator{T,ET}(k::AbstractArray{T}, bc::ET)

Defines an iterator over the knots in `k` based on the boundary conditions `bc`.

# Fields
- `knots::Vector{T}` The interpolated knots of the axis to iterate over
- `bc::ET` The Boundary Condition for the axis
- `nknots::Int` The number of interpolated knots (ie. `length(knots)`)

`ET` is `Union{BoundaryCondition,Tuple{BoundaryCondition,BoundaryCondition}}`

# Iterator Interface
The following methods defining Julia's iterator interface have been defined

`IteratorSize(::Type{KnotIterator})` -> Will return one of the following
- `Base.IsInfinite` if the iteration will produces an infinite sequence of knots
- `Base.HasLength` if iteration will produce a finite sequence of knots
- `Base.SizeUnknown` if we can't decided from only the type information

`length` and `size` -> Are defined if IteratorSize is HasLength, otherwise will
raise a MethodError.

`IteratorEltype` will always return `HasEltype`, as we always track the data
types of the knots

`eltype` will return the data type of the knots

`iterate` Defines iteration over the knots starting from the first one and
moving in the forward direction along the axis.

# Knots for Multi-dimensional Interpolants
Iteration over the knots of a multi-dimensional interpolant is done by wrapping
multiple KnotIterator within `Iterators.product`.

# Indexing
`KnotIterator` provides limited support for accessing knots via indexing
- `getindex` is provided for `KnotIterator` but does not support Multidimensional
  interpolations (As wrapped by `ProductIterator`) or non-Int indexes.
- A `BoundsError` will be raised if out of bounds and `checkbounds` has been
  implemented for `KnotIterator`

```jldoctest
julia> using Interpolations;

julia> etp = LinearInterpolation([1.0, 1.2, 2.3, 3.0], rand(4); extrapolation_bc=Periodic());

julia> kiter = knots(etp);

julia> kiter[4]
3.0

julia> kiter[36]
24.3
```

"""
struct KnotIterator{T,ET <: ExtrapSpec}
    knots::Vector{T}
    bc::ET
    nknots::Int
    KnotIterator{T,ET}(k::AbstractArray{T}, bc::ET) where {T,ET} = new(k, bc, length(k))
end

# Dispatch on outputs of k = getknots and bc = etpflag to create one
# KnotIterator per interpolation axis
function KnotIterator(k::NTuple{N,AbstractArray}, bc::NTuple{N, ExtrapSpec}) where {N}
    # One ExtrapSpec per Axis
    map(KnotIterator, k, bc)
end
function KnotIterator(k::NTuple{N,AbstractArray}, bc::BoundaryCondition) where {N}
    # All Axes use the same BoundaryCondition
    map(x -> KnotIterator(x, bc), k)
end
function KnotIterator(k::AbstractArray{T}, bc::ET) where {T,ET <: ExtrapDimSpec}
    # Prior dispatches will end up here: Axis knots: k and boundary condition bc
    KnotIterator{T,ET}(k, bc)
end

IteratorSize(::Type{KnotIterator}) = SizeUnknown()
IteratorSize(::Type{KnotIterator{T}}) where {T} = SizeUnknown()

IteratorSize(::Type{KnotIterator{T,ET}}) where {T,ET} = _knot_iter_size(ET)
_knot_iter_size(::Type{<:BoundaryCondition}) = HasLength()
_knot_iter_size(::Type{<:RepeatKnots}) = IsInfinite()
_knot_iter_size(::Type{Tuple{RevBC, FwdBC}}) where {RevBC,FwdBC} = _knot_iter_size(FwdBC)

length(iter::KnotIterator) = _knot_length(iter, IteratorSize(iter))
_knot_length(iter::KnotIterator, ::HasLength) = iter.nknots
size(iter::KnotIterator) = (length(iter),)

IteratorEltype(::Type{KnotIterator}) = EltypeUnknown()
IteratorEltype(::Type{<:KnotIterator{T}}) where {T} = HasEltype()
eltype(::Type{<:KnotIterator{T}}) where {T} = T

"""
    knots(itp::AbstractInterpolation)
    knots(etp::AbstractExtrapolation)

Returns an iterator over knot locations for an AbstractInterpolation or
AbstractExtrapolation.

Iterator will yield scalar values for interpolations over a single dimension,
and tuples of coordinates for higher dimension interpolations. Iteration over
higher dimensions is taken as the product of knots along each dimension.

ie. Iterator.product(knots on first dim, knots on 2nd dim,...)

Extrapolations with Periodic or Reflect boundary conditions, will produce an
infinite sequence of knots.

# Example
```jldoctest
julia> using Interpolations;

julia> etp = LinearInterpolation([1.0, 1.2, 2.3, 3.0], rand(4); extrapolation_bc=Periodic());

julia> Iterators.take(knots(etp), 5) |> collect
5-element $(Array{typeof(1.0),1}):
 1.0
 1.2
 2.3
 3.0
 3.2

```
"""
function knots(itp::AbstractInterpolation)
    # Construct separate KnotIterator for each dimension, and combine them using
    # Iterators.product
    k = getknots(itp)
    bc = Throw()
    iter = KnotIterator(k, bc)
    length(iter) == 1 ? iter[1] : Iterators.product(iter...)
end

function knots(etp::AbstractExtrapolation)
    # Construct separate KnotIterator for each dimension, and combine them using
    # Iterators.product
    k = getknots(etp)
    bc = etpflag(etp)
    iter = KnotIterator(k, bc)
    length(iter) == 1 ? iter[1] : Iterators.product(iter...)
end

# By default start iteration from the first knot, to start at a different index
# use knotsbetween
firstindex(::KnotIterator) = 1
iterate(iter::KnotIterator) = iterate(iter, firstindex(iter))

# Iterate over knots while checkbounds indicates the idx is still valid
function iterate(iter::KnotIterator{T,ET}, idx) where {T, ET}
    if checkbounds(Bool, iter, idx)
        iter[idx], idx+1
    else
        nothing
    end
end

# Split the ExtrapDimSpec into a tuple of left and right boundary conditions
splitExtrapDimSpec(::Type{ET}) where {ET <: BoundaryCondition} = (ET, ET)
splitExtrapDimSpec(::Type{<:Tuple{L,R}}) where {L,R} = (L,R)

# Checks if idx is a valid for the KnotIterator given it's boundary conditions
function checkbounds(::Type{Bool}, iter::KnotIterator{T,ET}, idx::Int) where {T,ET<:ExtrapDimSpec}
    # Get Left/Right Extrapolation Boundary Conditions for the KnotIterator
    left, right = splitExtrapDimSpec(ET)

    # Check Left/Right Boundary Limits -> If Left/Right is a repeating knot, the
    # checks is true as repeated knots generates an infinite sequence of knots
    leftcheck = left <: RepeatKnots ? true : 1 <= idx
    rightcheck = right <: RepeatKnots ? true : idx <= iter.nknots
    leftcheck && rightcheck
end

# Returns the knot given by idx, or raise a BoundsError
# Despatches to getknotindex(ET, iter, idx) if idx is an extrapolated knot,
# where ET is the boundary conditions for the relevant side
function getindex(iter::KnotIterator{T, ET}, idx::Int) where {T, ET <: ExtrapDimSpec}
    # Get Left/Right Extrapolation Boundary Conditions for the KnotIterator
    left, right = splitExtrapDimSpec(ET)

    # Dispatch to getknotindex if idx is an extrapolated knot, otherwise pull
    # from iter.knots
    if idx < 1
        getknotindex(left, iter, idx)::T
    elseif idx <= iter.nknots
        iter.knots[idx]
    else
        getknotindex(right, iter, idx)::T
    end
end

# For Non-repeating Knots -> Raise BoundsError
function getknotindex(::Type{<:BoundaryCondition}, iter, idx)
    throw(BoundsError(iter, idx))
end

# Get Knots using a Periodic Boundary Condition
@inline function getknotindex(::Type{<:Periodic}, iter, idx)
    knotrange = iter.knots[end] - iter.knots[1]
    offset = knotrange * (cld(idx, iter.nknots-1)-1)
    iter.knots[periodic(idx, 1, iter.nknots)] + offset
end

# Get Knots using a Reflect Boundary Condition
@inline function getknotindex(::Type{<:Reflect}, iter, idx)
    # Get the current knot
    knotrange = iter.knots[end] - iter.knots[1]
    offset = 2*knotrange * (cld(idx, 2*iter.nknots-2)-1)
    cycleidx = mod(idx-1, 2*iter.nknots-2)+1

    if iter.nknots < cycleidx
        offset + 2*iter.knots[end] - iter.knots[2*iter.nknots-cycleidx]
    else
        iter.knots[cycleidx] + offset
    end
end

"""
    nextknotidx(iter::KnotIterator, x)

Returns the index of the first knot such that `x < k` or `nothing` if no such
knot exists.

New boundary conditions should define:

    nextknotidx(::Type{<:NewBoundaryCondition}, knots::Vector, x)

Where `knots` is `iter.knots` and `NewBoundaryCondition` is the new boundary
conditions. This method is expected to handle values of `x` that are both inbounds
or extrapolated.
"""
function nextknotidx(iter::KnotIterator{T,ET}, x) where {T, ET}
    left, right = splitExtrapDimSpec(ET)
    if x < iter.knots[end]
        # Next knot is inbounds or to the left
        nextknotidx(left, iter.knots, x)
    else
        # Next knot is to the right
        nextknotidx(right, iter.knots, x)
    end
end

"""
    priorknotidx(iter::KnotIterator, x)

Returns the index of the last knot such that `k < x` or `nothing` ig not such
knot exists.

New boundary conditions should define

    priorknotidx(::Type{<:NewBoundaryCondition}, knots::Vector, x)

Where knots is `iter.knots` and `NewBoundaryCondition` is the new boundary
condition. This method is expected to handle values of `x` that are both inbounds
or extrapolated.
"""
function priorknotidx(iter::KnotIterator{T,ET}, x) where {T, ET}
    left, right = splitExtrapDimSpec(ET)
    if x <= iter.knots[1]
        # Prior knot is to the left
        priorknotidx(left, iter.knots, x)
    else
        # prior knot is inbounds or to the right
        priorknotidx(right, iter.knots, x)
    end
end

# priorknotidx and nextknotidx for non-repeating knots
function nextknotidx(::Type{<:BoundaryCondition}, knots::Vector, x)
    knots[end] <= x ? nothing : findfirst(x .< knots)
end
function priorknotidx(::Type{<:BoundaryCondition}, knots::Vector, x)
    x <= knots[1] ? nothing : findlast(knots .< x)
end

# Starting iterator state for a periodic knots
function nextknotidx(::Type{<:Periodic}, knots::Vector, start)
    # Find starting offset
    knotrange = knots[end] - knots[1]
    cycle = floor(Int, (start-knots[1])/ knotrange)

    # Find starting index
    cyclepos = periodic(start, knots[1], knots[end])
    cycleidx = findfirst(cyclepos .< knots)::Int
    if cycleidx === length(knots)
        cycleidx = 1
        cycle += 1
    end

    idx = (length(knots)-1) * cycle + cycleidx

    return idx
end
function priorknotidx(::Type{<:Periodic}, knots::Vector, stop)
    # Find stopping offset
    knotrange = knots[end] - knots[1]
    cycle = floor(Int, (stop-knots[1])/ knotrange)

    # Find stopping index
    cyclepos = periodic(stop, knots[1], knots[end])
    if cyclepos == knots[1]
        cycleidx = length(knots) -1
        cycle -= 1
    else
        cycleidx = findlast(knots .< cyclepos)::Int
    end
    idx = (length(knots)-1) * cycle + cycleidx

    return idx
end

# Starting iterator state for a reflecting knots
function nextknotidx(::Type{<:Reflect}, knots::Vector, start)
    refknots = knots[end] .+ (cumsum∘reverse∘diff)(knots)
    knots = vcat(knots, refknots)
    nextknotidx(Periodic, knots, start)
end
function priorknotidx(::Type{<:Reflect}, knots::Vector, stop)
    refknots = knots[end] .+ (cumsum∘reverse∘diff)(knots)
    knots = vcat(knots, refknots)
    priorknotidx(Periodic, knots, stop)
end

"""
    KnotRange(iter::KnotIterator{T}, start, stop)

Defines an iterator over a range of knots such that `start < k < stop`.

# Fields
- `iter::KnotIterator{T}` Underlying `KnotIterator` providing the knots iterated
- `range::R` Iterator defining the range of knot indices iterated. Where
  `R <: Union{Iterators.Count, UnitRange}`

# Iterator Interface
The following methods defining the Julia's iterator interface have been defined

`Base.IteratorSize` -> Will return one of the following:
- `Base.HasLength` if `range` is of finite length
- `Base.IsInfinite` if `range` is of infinite length
- `Base.SizeUnknown` if the type of `range` is unspecified

`Base.IteratorEltype` -> Returns `Base.EltypeUnknown` if type parameter not
provided, otherwise `Base.HasEltype`

`length` and `size` -> Returns the number of knots to be iterated if
`IteratorSize !== IsInfinite`, otherwise will raise `MethodError`

# Multidimensional Interpolants
Iteration over the knots of a multi-dimensional interpolant is done by wrapping
multiple `KnotRange` iterators within `Iterators.product`.
"""
struct KnotRange{T, R <: Union{Iterators.Count, UnitRange}}
    iter::KnotIterator{T}
    range::R
    KnotRange(iter::KnotIterator{T}, range::R) where {T,R} = new{T,R}(iter, range)
end

function KnotRange(iter::KnotIterator, start, stop)
    # Get the index of the first knot larger than start or firstindex(iter) iff
    # start is missing
    sdx = start === nothing ? firstindex(iter) : nextknotidx(iter, start)

    # Generate knot index range from provided stop / start
    if sdx === nothing
        # Empty Iterator
        range = 1:0
    elseif stop === nothing
        if IteratorSize(iter) === IsInfinite()
            # Stop is missing and iter generates infinite knots
            range = Iterators.countfrom(sdx)
        else
            # Stop is missing and iter generates a finite number of knots
            range = sdx:iter.nknots
        end
    else
        # Stop is provided -> Iterate from start to stop
        edx = priorknotidx(iter, stop)
        range = edx === nothing ? (1:0) : sdx:edx
    end
    KnotRange(iter, range)
end

# Iterator Size is Unknown until we knot range's type
IteratorSize(::Type{KnotRange}) = SizeUnknown()
IteratorSize(::Type{KnotRange{T}}) where {T} = SizeUnknown()
IteratorSize(::Type{KnotRange{T,R}}) where {T, R <: UnitRange} = HasLength()
IteratorSize(::Type{KnotRange{T,R}}) where {T, R <: Iterators.Count} = IsInfinite()

# If type parameter is not provided -> Eltype is Unknown otherwise known
IteratorEltype(::Type{<:KnotRange}) = EltypeUnknown()
IteratorEltype(::Type{<:KnotRange{T}}) where {T} = HasEltype()
eltype(::Type{<:KnotRange{T}}) where {T} = T

# Dispatch length and size to range
length(iter::KnotRange) = length(iter.range)
size(iter::KnotRange) = size(iter.range)

"""
    knotsbetween(iter; start, stop)
    knotsbetween(iter, start, stop)

Iterate over all knots of `iter` such that `start < k < stop`.

`iter` can be an `AbstractInterpolation`, or the output of `knots`
(ie. a `KnotIterator` or `ProductIterator` wrapping `KnotIterator`)

If `start` is not provided, iteration will start from the first knot. An
`ArgumentError` will be raised if both `start` and `stop` are not provided.

If no such knots exists will return a KnotIterator with length 0

# Example
```jldoctest
julia> using Interpolations;

julia> etp = LinearInterpolation([1.0, 1.2, 2.3, 3.0], rand(4); extrapolation_bc=Periodic());

julia> knotsbetween(etp; start=38, stop=42) |> collect
6-element $(Array{typeof(1.0),1}):
 38.3
 39.0
 39.2
 40.3
 41.0
 41.2

```
"""
knotsbetween(itp; start=nothing, stop=nothing) = knotsbetween(itp, start, stop)

# If AbstractInterpolation -> Get Knots -> Then Convert to Knot Range
knotsbetween(itp::AbstractInterpolation, start, stop) = knotsbetween(knots(itp), start, stop)

knotsbetween(::KnotIterator, ::Nothing, ::Nothing) =
    throw(ArgumentError("At least one of `start` or `stop` must be specified"))

knotsbetween(iter::KnotIterator, start, stop) = KnotRange(iter, start, stop)

# Expand Nothing on Start/Stop if tuples
function knotsbetween(iter::Iterators.ProductIterator, start::Union{Nothing, Tuple}, stop::Union{Nothing, Tuple})
    if start === stop === nothing
        throw(ArgumentError("At least one of `start` or `stop` must be specified"))
    else
        kiter = map_ranges(iter.iterators, start, stop)
        Iterators.product(kiter...)
    end
end

# Map iter, start and stop to N KnotRanges, handling start/stop being nothing
function map_ranges(iter, start::Tuple, ::Nothing)
    map((i, s) -> knotsbetween(i, s, nothing), iter, start)
end
function map_ranges(iter, ::Nothing, stop::Tuple)
    map((i, s) -> knotsbetween(i, nothing, s), iter, stop)
end
function map_ranges(iter, start::Tuple, stop::Tuple)
    map(knotsbetween, iter, start, stop)
end

iterate(iter::KnotRange) = iterate(iter, first(iter.range))

function iterate(iter::KnotRange{T,R}, state) where {T, R <: UnitRange}
    last(iter.range) < state && return nothing
    iter.iter[state], state+1
end

function iterate(iter::KnotRange{T,R}, state) where {T, R <: Iterators.Count}
    iter.iter[state], state+1
end

# Dispatch to KnotIterator for nextknotidx and priorknotidx
nextknotidx(iter::KnotRange, x) = nextknotidx(iter.iter, x)
priorknotidx(iter::KnotRange, x) = priorknotidx(iter.iter, x)


text:
# # Vlasov-HMF
#
#md # [notebook](@__NBVIEWER_ROOT_URL__notebooks/vlasov-hmf.ipynb)
#

using LinearAlgebra, QuadGK, Roots, FFTW
using Splittings
using Plots

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

" Compute M₀ by solving F(m) = 0 "
function mag(β, mass)

    F(m) = begin
        g(x, n, m) = (1 / π) * (exp(β * m * cos(x)) * cos(n * x))
        bessel0(x) = g(x, 0, m)
        bessel1(x) = g(x, 1, m)
        mass * quadgk(bessel1, 0, π)[1] / quadgk(bessel0, 0, π)[1] - m
    end

    find_zero(F, (0, mass))
end

#-

function Norm(f::Array{Float64,2}, delta1, delta2)
   return delta1 * sum(delta2 * sum(real(f), dims=1))
end

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

"""
    Compute the electric hamiltonian mean field from a
    2D distribution function
"""
function hmf_poisson!(fᵗ::Array{Complex{Float64},2},
        mesh1::UniformMesh,
        mesh2::UniformMesh,
        ex::Array{Float64})

    n1 = mesh1.length
    rho = mesh2.step .* vec(sum(fᵗ, dims=1))
    kernel = zeros(Float64, n1)
    k = π / (mesh1.stop - mesh1.start)
    kernel[2] = k
    ex .= real(ifft(1im * fft(rho) .* kernel * 4π ))

end

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

function bsl_advection!(f::Array{Complex{Float64},2},
                        mesh1::UniformMesh,
                        mesh2::UniformMesh,
                        v::Array{Float64,1},
                        dt;
                        spline_degree=3)

    fft!(f,1)
    @simd for j in 1:mesh2.length
        alpha = v[j] * dt
        @inbounds f[:,j] .= Splittings.interpolate(spline_degree, f[:,j],
            mesh1.step, alpha)
    end
    ifft!(f,1)
end

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

function push_v!(f, fᵗ, mesh1, mesh2, ex, dt)
    transpose!(fᵗ, f)
    hmf_poisson!(fᵗ, mesh1, mesh2, ex)
    bsl_advection!(fᵗ, mesh2, mesh1, ex, dt)
    transpose!(f, fᵗ)
end

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

function vlasov_hmf_gauss(nbiter = 10000, dt = 0.1)

    mass = 1.0
    T = 0.1
    mesh1 = UniformMesh(-π, π, 64)
    mesh2 = UniformMesh(-8, 8, 64)

    n1, delta1 = mesh1.length, mesh1.step
    n2, delta2 = mesh2.length, mesh2.step
    x, v = mesh1.points, mesh2.points
    X = repeat(x,1,n2)
    V = repeat(v,1,n1)'
    ϵ = 0.1

    b = 1 / T
    m = mag(b, mass)

    w   = sqrt(m)
    f   = zeros(Complex{Float64}, (n1,n2))
    fᵗ  = zeros(Complex{Float64}, (n2,n1))

    f  .= exp.(-b .* ((V.^2 / 2) - m * cos.(X)))
    a   = mass / Norm(real(f), delta1, delta2)
    @.  f =  a * exp(-b * (((V^2) / 2) - m * cos(X))) * (1 + ϵ * cos(X))

    ex = zeros(Float64,n1)
    hmf_poisson!(f, mesh1, mesh2, ex )
    test = copy(f)
    T = Float64[]
    for n in 1:nbiter

        gamma1 = Norm(real(f) .* cos.(X), delta1, delta2)
        push!(T,gamma1)

        @Strang(
            bsl_advection!(f, mesh1, mesh2, v, dt),
            push_v!(f, fᵗ, mesh1, mesh2, ex, dt)
        )

    end

    #Substracting from gamma its long time average

    Gamma1 = Norm(real(f) .*cos.(X), delta1, delta2)
    T .= T .- Gamma1

    range(0., stop=nbiter*deltat, length=nbiter), abs.(T)

end

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

nbiter = 1000
deltat = 0.1
@time t, T = vlasov_hmf_gauss(nbiter, deltat);
plot(t, log.(T), xlabel = "t", ylabel = "|C[f](t)-C[f][T]|")
savefig("vlasov-hmf-plot.png"); nothing # hide

#
#md # ![png](vlasov-hmf-plot.png)
#
using Test                   #src
@testset "Vlasov-HMF" begin  #src
@test length(T) > 0          #src
end                          #src


text:
struct UncertainVariableRef <: JuMP.AbstractVariableRef
    model::RobustModel # `model` owning the variable
    idx::Int       # Index in `model.variables`
end
Base.copy(v::UncertainVariableRef) = v

function Base.:(==)(v::UncertainVariableRef, w::UncertainVariableRef)
    return v.model === w.model && v.idx == w.idx
end

Base.broadcastable(v::UncertainVariableRef) = Ref(v)
JuMP.isequal_canonical(v::UncertainVariableRef, w::UncertainVariableRef) = v == w
function JuMP.add_variable(
    m::RobustModel,
    v::JuMP.AbstractVariable,
    name::String = "",
)
    m.nextvaridx += 1
    vref = UncertainVariableRef(m, m.nextvaridx)
    m.uncertainVariables[vref.idx] = v
    JuMP.set_name(vref, name)
    return vref
end
function JuMP.add_variable(
    model::RobustModel,
    variable::JuMP.VariableConstrainedOnCreation,
    name::String,
)
    var_ref = JuMP.add_variable(model, variable.scalar_variable, name)
    JuMP.add_constraint(model, JuMP.ScalarConstraint(var_ref, variable.set))
    return var_ref
end
function JuMP.add_variable(
    model::RobustModel,
    variable::JuMP.VariablesConstrainedOnCreation,
    names,
)
    var_refs =
        JuMP.add_variable.(
            model,
            variable.scalar_variables,
            JuMP.vectorize(names, variable.shape),
        )
    JuMP.add_constraint(model, JuMP.VectorConstraint(var_refs, variable.set))
    return JuMP.reshape_vector(var_refs, variable.shape)
end

function JuMP.delete(model::RobustModel, vref::UncertainVariableRef)
    @assert JuMP.is_valid(model, vref)
    delete!(model.uncertainVariables, vref.idx)
    return delete!(model.var_to_name, vref.idx)
end
function JuMP.delete(model::RobustModel, vrefs::Vector{UncertainVariableRef})
    return JuMP.delete.(model, vrefs)
end
function JuMP.is_valid(model::RobustModel, vref::UncertainVariableRef)
    return (model === vref.model && vref.idx in keys(model.uncertainVariables))
end
JuMP.num_variables(m::RobustModel) = length(m.uncertainVariables)

variable_info(vref::UncertainVariableRef) = vref.model.uncertainVariables[vref.idx].info
function update_variable_info(vref::UncertainVariableRef, info::JuMP.VariableInfo)
    return vref.model.uncertainVariables[vref.idx] = JuMP.ScalarVariable(info)
end

JuMP.has_lower_bound(vref::UncertainVariableRef) = variable_info(vref).has_lb
function JuMP.lower_bound(vref::UncertainVariableRef)::Float64
    @assert !JuMP.is_fixed(vref)
    return variable_info(vref).lower_bound
end
function JuMP.set_lower_bound(vref::UncertainVariableRef, lower)
    info = variable_info(vref)
    return update_variable_info(
        vref,
        JuMP.VariableInfo(
            true,
            lower,
            info.has_ub,
            info.upper_bound,
            info.has_fix,
            info.fixed_value,
            info.has_start,
            info.start,
            info.binary,
            info.integer,
        ),
    )
end

function JuMP.delete_lower_bound(vref::UncertainVariableRef)
    info = variable_info(vref)
    return update_variable_info(
        vref,
        JuMP.VariableInfo(
            false,
            info.lower_bound,
            info.has_ub,
            info.upper_bound,
            info.has_fix,
            info.fixed_value,
            info.has_start,
            info.start,
            info.binary,
            info.integer,
        ),
    )
end
JuMP.has_upper_bound(vref::UncertainVariableRef) = variable_info(vref).has_ub
function JuMP.upper_bound(vref::UncertainVariableRef)::Float64
    @assert !JuMP.is_fixed(vref)
    return variable_info(vref).upper_bound
end
function JuMP.set_upper_bound(vref::UncertainVariableRef, upper)
    info = variable_info(vref)
    return update_variable_info(
        vref,
        JuMP.VariableInfo(
            info.has_lb,
            info.lower_bound,
            true,
            upper,
            info.has_fix,
            info.fixed_value,
            info.has_start,
            info.start,
            info.binary,
            info.integer,
        ),
    )
end
function JuMP.delete_upper_bound(vref::UncertainVariableRef)
    info = variable_info(vref)
    return update_variable_info(
        vref,
        JuMP.VariableInfo(
            info.has_lb,
            info.lower_bound,
            false,
            info.upper_bound,
            info.has_fix,
            info.fixed_value,
            info.has_start,
            info.start,
            info.binary,
            info.integer,
        ),
    )
end
JuMP.is_fixed(vref::UncertainVariableRef) = variable_info(vref).has_fix
function JuMP.fix_value(vref::UncertainVariableRef)::Float64
    return variable_info(vref).fixed_value
end
function JuMP.fix(vref::UncertainVariableRef, value; force::Bool = false)
    info = variable_info(vref)
    if !force && (info.has_lb || info.has_ub)
        error(
            "Unable to fix $(vref) to $(value) because it has existing bounds.",
        )
    end
    update_variable_info(
        vref,
        JuMP.VariableInfo(
            false,
            info.lower_bound,
            false,
            info.upper_bound,
            true,
            value,
            info.has_start,
            info.start,
            info.binary,
            info.integer,
        ),
    )
    return
end
function JuMP.unfix(vref::UncertainVariableRef)
    info = variable_info(vref)
    return update_variable_info(
        vref,
        JuMP.VariableInfo(
            info.has_lb,
            info.lower_bound,
            info.has_ub,
            info.upper_bound,
            false,
            info.fixed_value,
            info.has_start,
            info.start,
            info.binary,
            info.integer,
        ),
    )
end
function JuMP.start_value(vref::UncertainVariableRef)::Union{Nothing,Float64}
    return variable_info(vref).start
end
function JuMP.set_start_value(vref::UncertainVariableRef, start)
    info = variable_info(vref)
    return update_variable_info(
        vref,
        JuMP.VariableInfo(
            info.has_lb,
            info.lower_bound,
            info.has_ub,
            info.upper_bound,
            info.has_fix,
            info.fixed_value,
            true,
            start,
            info.binary,
            info.integer,
        ),
    )
end
JuMP.is_binary(vref::UncertainVariableRef) = variable_info(vref).binary
function JuMP.set_binary(vref::UncertainVariableRef)
    @assert !JuMP.is_integer(vref)
    info = variable_info(vref)
    return update_variable_info(
        vref,
        JuMP.VariableInfo(
            info.has_lb,
            info.lower_bound,
            info.has_ub,
            info.upper_bound,
            info.has_fix,
            info.fixed_value,
            info.has_start,
            info.start,
            true,
            info.integer,
        ),
    )
end
function JuMP.unset_binary(vref::UncertainVariableRef)
    info = variable_info(vref)
    return update_variable_info(
        vref,
        JuMP.VariableInfo(
            info.has_lb,
            info.lower_bound,
            info.has_ub,
            info.upper_bound,
            info.has_fix,
            info.fixed_value,
            info.has_start,
            info.start,
            false,
            info.integer,
        ),
    )
end
JuMP.is_integer(vref::UncertainVariableRef) = variable_info(vref).integer
function JuMP.set_integer(vref::UncertainVariableRef)
    @assert !JuMP.is_binary(vref)
    info = variable_info(vref)
    return update_variable_info(
        vref,
        JuMP.VariableInfo(
            info.has_lb,
            info.lower_bound,
            info.has_ub,
            info.upper_bound,
            info.has_fix,
            info.fixed_value,
            info.has_start,
            info.start,
            info.binary,
            true,
        ),
    )
end
function JuMP.unset_integer(vref::UncertainVariableRef)
    info = variable_info(vref)
    return update_variable_info(
        vref,
        JuMP.VariableInfo(
            info.has_lb,
            info.lower_bound,
            info.has_ub,
            info.upper_bound,
            info.has_fix,
            info.fixed_value,
            info.has_start,
            info.start,
            info.binary,
            false,
        ),
    )
end


text:
include("nonlinear/forward_pass.jl")
include("nonlinear/reverse_pass.jl")
include("nonlinear/univariate.jl")
include("bound.jl")
include("relax.jl")


text:
using Plots, ProgressMeter
pyplot(leg=false, ticks=nothing) #change to the pyplot backend and define some defaults
x = y = range(-5, stop = 5, length = 40)
zs = zeros(0,40)
n = 100

# create a progress bar for tracking the animation generation
prog = Progress(n,1)

@gif for i in range(0, stop = 2π, length = n)
    f(x,y) = sin(x + 10sin(i)) + cos(y)

    # create a plot with 3 subplots and a custom layout
    l = @layout [a{0.7w} b; c{0.2h}]
    p = plot(x, y, f, st = [:surface, :contourf], layout=l)

    # induce a slight oscillating camera angle sweep, in degrees (azimuth, altitude)
    plot!(p[1], camera=(15*cos(i), 40))

    # add a tracking line
    fixed_x = zeros(40)
    z = map(f, fixed_x, y)
    plot!(p[1], fixed_x, y, z, line = (:black, 5, 0.2))
    vline!(p[2], [0], line = (:black, 5))

    # add to and show the tracked values over time
    global zs = vcat(zs, z')
    plot!(p[3], zs, alpha = 0.2, palette = cgrad(:blues).colors)

    # increment the progress bar
    next!(prog)
end


text:
# Modified from Shreyas Kowshik's implementation.

# Weight initialization
_random_normal(shape...) = map(Float32, rand(Normal(0, 0.1), shape...))
constant_init(shape...) = map(Float32, ones(shape...) * 0.1)
normalize(arr) = (arr .- mean(arr))./(sqrt(var(arr) + 1e-10))

function normalize_across_procs(arr, episode_length)
    arr = reshape(arr, episode_length, NUM_PROCESSES)
    reshape((arr .- mean(arr,dims=2))./(sqrt.(var(arr, dims=2) .+ 1e-10)), 1, episode_length*NUM_PROCESSES)
end

### Utility functions specific to PPO

"""
Returns a Generalized Advantage Estimate for an episode
"""
function gae(policy, states::Array, actions::Array, rewards::Array, next_states::Array, num_steps::Int; γ=0.99, λ=0.95)
    Â = []
    A = 0.0
    for i in reverse(1:length(states))
        if length(states) < num_steps && i == length(states)
            δ = rewards[i] - policy.value_net(states[i]).data[1]
        else
            δ = rewards[i] + γ*policy.value_net(next_states[i]).data[1] - policy.value_net(states[i]).data[1]
        end

        A = δ + (γ*λ*A)
        push!(Â, A)
    end

    return reverse(Â)
end

"""
Returns the cumulative discounted returns for each timestep
"""
function disconunted_returns(rewards::Array, last_val=0; γ=0.99)
    r = 0.0
    returns = []

    for i in reverse(1:length(rewards))
        r = rewards[i] + γ*r  # TODO. γ^i (?)
        if i == length(rewards)
            r = r + last_val
        end
        push!(returns, r)
    end

    return reverse(returns)
end

### Utility functions specific to TRPO

"""
Obtain the gradient vector product. Intermediate utility function, calculates `Σ∇D_kl*x`

`x` : Variable to be estimated using conjugate gradient `(Hx = g); (NUM_PARAMS,1)`
"""
function gvp(policy, states, kl_vars, x)
    model_params = get_policy_params(policy)
    gs = Tracker.gradient(() -> kl_loss(policy, states, kl_vars), model_params; nest=true)

    flat_grads = get_flat_grads(gs, get_policy_net(policy))
    return sum(x' * flat_grads)
end

"""
Computes the Hessian vector product.
Hessian is that of the kl divergence between the old and the new policies w.r.t. the policy parameters.

Returns : `Hx; H = ∇²D_kl`
"""
function Hvp(policy, states, kl_vars, x; damping_coeff=0.1)
    model_params = get_policy_params(policy)
    hessian = Tracker.gradient(() -> gvp(policy, states, kl_vars, x), model_params)
    return get_flat_grads(hessian, get_policy_net(policy)) # .+ (damping_coeff .* x)
end

"""
Compute the conjugate gradient.

`b` : Array of shape `(NUM_PARAMS,1)`

Solves `x` for `Hx = b`
"""
function conjugate_gradients(policy, states, kl_vars, b, nsteps=10, err=1e-10)
    x = zeros(size(b))
    r = copy(b)
    p = copy(b)
    rdotr = r' * r

    for i in 1:nsteps
        hvp = Hvp(policy, states, kl_vars, p).data # Returns array of shape (NUM_PARAMS,1)

        α = rdotr ./ (p' * hvp)

        x = x .+ (α .* p)
        r = r .- (α .* hvp)

        new_rdotr = r' * r
        β = new_rdotr ./ rdotr
        p = r .+ (β .* p)

        rdotr = new_rdotr

        if rdotr[1] < err
            break
        end
    end

    return x
end

"""
Flattens out the gradients and concatenates them.

`models` : An array of models whose parameter `gradients` are to be falttened.

Returns : Tracked Array of shape `(NUM_PARAMS,1)`
"""
function get_flat_grads(gradients, models)

    flat_grads = []

    function flatten!(p)
        if typeof(p) <: TrackedArray
            prod_size = prod(size(p))
            push!(flat_grads, reshape(gradients[p], prod_size))
        end
    end

    for model in models
        mapleaves(flatten!, model)
    end

    flat_grads = cat(flat_grads..., dims=1)
    flat_grads = reshape(flat_grads, length(flat_grads), 1)

    return flat_grads
end

"""
Flattens out the parameters and concatenates them.

`models` : An array of models whose parameters are to be flattened

Returns : Tracked Array of shape `(NUM_PARAMS,1)`
"""
function get_flat_params(models)

    flat_params = []

    function flatten!(p)
        if typeof(p) <: TrackedArray
            prod_size = prod(size(p))
            push!(flat_params, reshape(p, prod_size))
        end
    end

    for model in models
        mapleaves(flatten!, model)
    end

    flat_params = cat(flat_params..., dims=1)
    flat_params = reshape(flat_params, length(flat_params), 1)

    return flat_params
end

"""
Sets values of `parameters` to the `model`.

`parameters` : flattened out array of model parameters.

`models` : an array of models whose parameters are to be set.
"""
function set_flat_params!(parameters, models)
    ptr = 1

    function assign!(p)
        if typeof(p) <: TrackedArray
            prod_size = prod(size(p))

            p.data .= Float32.(reshape(parameters[ptr:ptr + prod_size - 1, :], size(p)...)).data
            ptr += prod_size
        end
    end

    for model in models
        mapleaves(assign!, model)
    end
end


text:
#   This file is part of UnitSystems.jl. It is licensed under the MIT license
#   UnitSystems Copyright (C) 2020 Michael Reed

"""
$(convertext(:charge,"sqrt(planck(U,S)*current(U,S)/electricpotential(U,S))"))

Electric `charge` quantization (C, A⋅s), unit conversion factor.

```Julia
julia> charge(EMU,Metric) # C⋅abC⁻¹
$(charge(EMU,Metric))

julia> charge(EMU,ESU) # stC⋅abC⁻¹
$(charge(EMU,ESU))

julia> charge(ESU,Metric) # C⋅stC⁻¹
$(charge(ESU,Metric))

julia> charge(Metric,SI2019) # C⋅C⁻¹
$(charge(Metric,SI2019))

julia> charge(Hartree,SI2019) # C⋅𝘦⁻¹
$(charge(Hartree,SI2019))
```
"""
@pure charge(U::UnitSystem,S::UnitSystem) = unit(sqrt((planckreduced(S)*permeability(U)*lightspeed(U)*rationalization(U)*lorentz(U)^2)/(planckreduced(U)*permeability(S)*lightspeed(S)*rationalization(S)*lorentz(S)^2)))

"""
$(convertext(:current,"charge(U,S)/time(U,S)"))

Flow of electric `charge` per `time` or `current` (A, C⋅s⁻¹), unit conversion factor.

```Julia
julia> current(EMU,Metric) # A⋅Bi⁻¹
$(current(EMU,Metric))

julia> current(EMU,ESU) # statA⋅Bi⁻¹
$(current(EMU,ESU))

julia> current(ESU,Metric) # A⋅statA⁻¹
$(current(ESU,Metric))

julia> current(Metric,SI2019) # A⋅A⁻¹
$(current(Metric,SI2019))
```
"""
@pure current(U::UnitSystem,S::UnitSystem) = unit(charge(U,S)/time(U,S))

"""
$(convertext(:electricpotential,"energy(U,S)/charge(U,S)"))

Voltage or `electricpotential` or `energy` per `charge` (V, J⋅C⁻¹), unit conversion factor.

```Julia
julia> electricpotential(EMU,Metric) # V⋅abV⁻¹
$(electricpotential(EMU,Metric))

julia> electricpotential(EMU,ESU) # statV⋅abV⁻¹
$(electricpotential(EMU,ESU))

julia> electricpotential(ESU,Metric) # V⋅statV⁻¹
$(electricpotential(ESU,Metric))

julia> electricpotential(Metric,SI2019) # V⋅V⁻¹
$(electricpotential(Metric,SI2019))
```
"""
@pure electricpotential(U::UnitSystem,S::UnitSystem) = unit(energy(U,S)/charge(U,S))
const voltage = electricpotential

"""
$(convertext(:capacitance,"charge(U,S)/electricpotential(U,S)"))

Electrical `capactiance` or `charge` per `electricpotential` (F, C⋅V⁻¹), unit conversion factor.

```Julia
julia> capacitance(EMU,Metric) # F⋅abF⁻¹
$(capacitance(EMU,Metric))

julia> capacitance(ESU,Metric) # F⋅cm⁻¹
$(capacitance(ESU,Metric))

julia> capactiance(Metric,SI2019) # F⋅F⁻¹
$(capacitance(Metric,SI2019))
```
"""
@pure capacitance(U::UnitSystem,S::UnitSystem) = unit(charge(U,S)/electricpotential(U,S))

"""
$(convertext(:resistance,"electricpotential(U,S)/current(U,S)"))

Electrical `resistance` or `electricpotential` per `current` (Ω, S⁻¹, V⋅A⁻¹), unit conversion factor.

```Julia
julia> resistance(EMU,Metric) # Ω⋅abΩ⁻¹
$(resistance(EMU,Metric))

julia> resistance(ESU,Metric) # Ω⋅statΩ⁻¹
$(resistance(ESU,Metric))

julia> resistance(Metric,SI2019) # Ω⋅Ω⁻¹
$(resistance(Metric,SI2019))
```
"""
@pure resistance(U::UnitSystem,S::UnitSystem) = unit(electricpotential(U,S)/current(U,S))

"""
$(convertext(:conductance,"current(U,S)/electricpotential(U,S)"))

Electrical `conductance` or `current` per `electricpotential` (S, Ω⁻¹, A⋅V⁻¹), unit conversion factor.

```Julia
julia> conductance(EMU,Metric) # S⋅abS⁻¹
$(conductance(EMU,Metric))

julia> conductance(ESU,Metric) # S⋅statS⁻¹
$(conductance(ESU,Metric))

julia> conductance(Metric,SI2019) # S⋅S⁻¹
$(conductance(Metric,SI2019))
```
"""
@pure conductance(U::UnitSystem,S::UnitSystem) = unit(current(U,S)/electricpotential(U,S))

"""
$(convertext(:magneticflux,"energy(U,S)/lorentz(U,S)/current(U,S)"))

Surface `magneticflux` or `energy` per `current` (Wb, J⋅A⁻¹, V⋅s), unit conversion factor.

```Julia
julia> magneticflux(EMU,Metric) # Wb⋅Mx⁻¹
$(magneticflux(EMU,Metric))

julia> magneticflux(ESU,Metric) # Wb⋅statWb⁻¹
$(magneticflux(ESU,Metric))

julia> magneticflux(Metric,SI2019) # Wb⋅Wb⁻¹
$(magneticflux(Metric,SI2019))
```
"""
@pure magneticflux(U::UnitSystem,S::UnitSystem) = unit(energy(U,S)/lorentz(U,S)/current(U,S))

"""
$(convertext(:magneticfluxdensity,"mass(U,S)/lorentz(U,S)/current(U,S)/time(U,S)^2"))

Magnetic induction or surface `magneticfluxdensity` (T, Wb⋅m⁻²), unit conversion factor.

```Julia
julia> magneticfluxdensity(EMU,Metric) # T⋅G⁻¹
$(magneticfluxdensity(EMU,Metric))

julia> magneticfluxdensity(EMU,ESU) # statT⋅G⁻¹
$(magneticfluxdensity(EMU,ESU))

julia> magneticfluxdensity(Metric,SI2019) # T⋅T⁻¹
$(magneticfluxdensity(Metric,SI2019))
```
"""
@pure magneticfluxdensity(U::UnitSystem,S::UnitSystem) = unit(mass(U,S)/lorentz(U,S)/current(U,S)/time(U,S)^2)

"""
$(convertext(:inductance,"mass(U,S)*area(U,S)/charge(U,S)"))

Electro-`magneticflux` per `current` or `inductance` (H, Ω⋅s, Wb⋅A⁻¹), unit conversion factor.

```Julia
julia> inductance(EMU,Metric) # H⋅abH⁻¹
$(inductance(EMU,Metric))

julia> inductance(ESU,Metric) # H⋅statH⁻¹
$(inductance(ESU,Metric))

julia> inductance(Metric,SI2019) # H⋅H⁻¹
$(inductance(Metric,SI2019))
```
"""
@pure inductance(U::UnitSystem,S::UnitSystem) = unit(mass(U,S)*area(U,S)/charge(U,S)^2)

# electromagnetics

"""
$(convertext(:electricfluxdensity,"charge(U,S)*rationalization(U,S)/area(U,S)"))

Electric field displacement or surface `electricfluxdensity` (C⋅m⁻²), unit conversion factor.

```Julia
julia> electricfluxdensity(EMU,Metric) # C⋅cm²⋅abC⁻¹⋅m⁻²
$(electricfluxdensity(EMU,Metric))

julia> electricfluxdensity(ESU,Metric) # C⋅cm²⋅statC⁻¹⋅m⁼²
$(electricfluxdensity(ESU,Metric))

julia> electricfluxdensity(Metric,SI2019) # C⋅C⁻¹
$(electricfluxdensity(Metric,SI2019))
```
"""
@pure electricfluxdensity(U::UnitSystem,S::UnitSystem) = unit(charge(U,S)*rationalization(U,S)/area(U,S))

"""
$(convertext(:chargedensity,"charge(U,S)/volume(U,S)"))

Volume `chargedensity` or `charge` per `volume` (C⋅m⁻³), unit conversion factor.

```Julia
julia> chargedensity(EMU,Metric) # C⋅cm³⋅abC⁻¹⋅m⁻³
$(chargedensity(EMU,Metric))

julia> chargedensity(ESU,Metric) # C⋅cm³⋅statC⁻¹⋅m⁼³
$(chargedensity(ESU,Metric))

julia> chargedensity(Metric,SI2019) # C⋅C⁻¹
$(chargedensity(Metric,SI2019))
```
"""
@pure chargedensity(U::UnitSystem,S::UnitSystem) = unit(charge(U,S)/volume(U,S))

"""
$(convertext(:currentdensity,"current(U,S)/area(U,S)"))

Cross-section `currentdensity` or `current` per `area` (A⋅m⁻²), unit conversion factor.

```Julia
julia> currentdensity(EMU,Metric) # A⋅cm²⋅Bi⁻¹⋅m⁻²
$(currentdensity(EMU,Metric))

julia> currentdensity(ESU,Metric) # A⋅cm²⋅statA⁻¹⋅m⁼²
$(currentdensity(ESU,Metric))

julia> currentdensity(Metric,SI2019) # A⋅A⁻¹
$(currentdensity(Metric,SI2019))
```
"""
@pure currentdensity(U::UnitSystem,S::UnitSystem) = unit(current(U,S)/area(U,S))

"""
$(convertext(:conductivity,"conductance(U,S)/length(U,S)"))

Reciprocal `resistivity` or electrical `conductivity` (S⋅m⁻¹), unit conversion factor.

```Julia
julia> conductivity(EMU,Metric) # S⋅cm⋅abS⁻¹⋅m⁻¹
$(conductivity(EMU,Metric))

julia> conductivity(ESU,Metric) # S⋅cm⋅statS⁻¹⋅m⁼¹
$(conductivity(ESU,Metric))

julia> conductivity(Metric,SI2019) # S⋅S⁻¹
$(conductivity(Metric,SI2019))
```
"""
@pure conductivity(U::UnitSystem,S::UnitSystem) = unit(conductance(U,S)/length(U,S))

"""
$(convertext(:permittivity,"capacitance(U,S)*rationalization(U,S)/length(U,S)"))

Absolute `permittivity` or `capacitance` per `length` (F⋅m⁻¹), unit conversion factor.

```Julia
julia> permittivity(EMU,Metric) # F⋅cm⋅abF⁻¹⋅m⁻¹
$(permittivity(EMU,Metric))

julia> permittivity(ESU,Metric) # F⋅m⁼¹
$(permittivity(ESU,Metric))

julia> permittivity(Metric,SI2019) # F⋅F⁻¹
$(permittivity(Metric,SI2019))
```
"""
@pure permittivity(U::UnitSystem,S::UnitSystem) = unit(capacitance(U,S)*rationalization(U,S)/length(U,S))

@doc """
$(convertext(:permeability,"permeability(S)/permeability(U)"))

Magnetic `permeability` or `inductance` per `length` (H⋅m⁻¹), unit conversion factor.

```Julia
julia> permeability(EMU,Metric) # H⋅cm⋅abH⁻¹⋅m⁻¹
$(permeability(EMU,Metric))

julia> permeability(ESU,Metric) # H⋅cm⋅statH⁻¹⋅m⁼¹
$(permeability(ESU,Metric))

julia> permeability(Metric,SI2019) # H⋅H⁻¹
$(permeability(Metric,SI2019))
```
""" permeability(::UnitSystem,::UnitSystem)

"""
$(convertext(:electricfield,"electricpotential(U,S)/length(U,S)"))

The `electricpotential` per `length` or `electricfield` (V⋅m⁻¹), unit conversion factor.

```Julia
julia> electricfield(EMU,Metric) # V⋅cm⋅abV⁻¹⋅m⁻¹
$(electricfield(EMU,Metric))

julia> electricfield(EMU,ESU) # statV⋅abV⁻¹
$(electricfield(EMU,ESU))

julia> electricfield(ESU,Metric) # V⋅cm⋅statV⁻¹⋅m⁻¹
$(electricfield(ESU,Metric))

julia> electricfield(Metric,SI2019) # V⋅V⁻¹
$(electricfield(Metric,SI2019))
```
"""
@pure electricfield(U::UnitSystem,S::UnitSystem) = unit(electricpotential(U,S)/length(U,S))

"""
$(convertext(:magneticfield,"current(U,S)*rationalization(U,S)*lorentz(U,S)/length(U,S)"))

Magnetization or `magneticfield` or `current` per `length` (A⋅m⁻¹), unit conversion factor.

```Julia
julia> magneticfield(EMU,Metric) # A⋅m⁻¹⋅Oe⁻¹
$(magneticfield(EMU,Metric))

julia> magneticfield(ESU,Metric) # A⋅cm⋅m⁻¹⋅statA⁻¹
$(magneticfield(ESU,Metric))

julia> magneticfield(Metric,SI2019) # A⋅A⁻¹
$(magneticfield(Metric,SI2019))
```
"""
@pure magneticfield(U::UnitSystem,S::UnitSystem) = unit(current(U,S)*rationalization(U,S)*lorentz(U,S)/length(U,S))

"""
$(convertext(:exposure,"charge(U,S)/mass(U,S)"))

Ionizing radiation `exposure` or `charge` per `mass` (C⋅kg⁻¹), unit conversion factor.

```Julia
julia> exposure(EMU,Metric) # C⋅g⋅abC⁻¹⋅kg
$(exposure(EMU,Metric))

julia> exposure(EMU,ESU) # statC⋅abC⁻¹
$(exposure(EMU,ESU))

julia> expsure(ESU,Metric) # C⋅g⋅statC⁻¹⋅kg
$(exposure(ESU,Metric))

julia> exposure(Metric,SI2019) # C⋅C⁻¹
$(exposure(Metric,SI2019))
```
"""
@pure exposure(U::UnitSystem,S::UnitSystem) = unit(charge(U,S)/mass(U,S))

"""
$(convertext(:resistivity,"resistance(U,S)*length(U,S)"))

Electrical `resistivity` or `resistance` by `length` (Ω⋅m), unit conversion factor.

```Julia
julia> resistance(EMU,Metric) # Ω⋅m⋅abΩ⁻¹⋅cm⁻¹
$(resistance(EMU,Metric))

julia> resistance(ESU,Metric) # Ω⋅m⋅statΩ⁻¹⋅cm⁻¹
$(resistance(ESU,Metric))

julia> resistance(Metric,SI2019) # Ω⋅Ω⁻¹
$(resistance(Metric,SI2019))
```
"""
@pure resistivity(U::UnitSystem,S::UnitSystem) = unit(resistance(U,S)*length(U,S))

"""
$(convertext(:linearchargedensity,"charge(U,S)/length(U,S)"))

Amount of `linearchargedensity` or `charge` per `length` (C⋅m⁻¹), unit conversion factor.

```Julia
julia> linearchargedensity(EMU,Metric) # C⋅cm⋅abC⁻¹⋅m⁻¹
$(linearchargedensity(EMU,Metric))

julia> linearchargedensity(ESU,Metric) # C⋅cm⋅statC⁻¹⋅m⁼¹
$(linearchargedensity(ESU,Metric))

julia> linearchargedensity(Metric,SI2019) # C⋅C⁻¹
$(linearchargedensity(Metric,SI2019))
```
"""
@pure linearchargedensity(U::UnitSystem,S::UnitSystem) = unit(charge(U,S)/length(U,S))

"""
$(convertext(:magneticdipolemoment,"current(U,S)*lorentz(U,S)/area(U,S)"))

Magnetic dipole moment or `magneticdipolemoment` (J⋅T⁻¹, A⋅m²), unit conversion factor.

```Julia
julia> magneticdipolemoment(EMU,Metric) # J⋅G⋅T⁻¹⋅erg⁻¹
$(magneticdipolemoment(EMU,Metric))

julia> magneticdipolemoment(ESU,Metric) # J⋅statT⋅T⁻¹⋅erg⁼¹
$(magneticdipolemoment(ESU,Metric))

julia> magneticdipolemoment(Metric,SI2019) # A⋅A⁻¹⋅
$(magneticdipolemoment(Metric,SI2019))
```
"""
@pure magneticdipolemoment(U::UnitSystem,S::UnitSystem) = unit(current(U,S)*lorentz(U,S)*area(U,S))

"""
$(convertext(:mobility,"charge(U,S)*time(U,S)/mass(U,S)"))

Electron `mobility` in solid state physics (m²⋅V⁻¹⋅s⁻¹, A⋅s²⋅kg⁻¹), unit conversion factor.

```Julia
julia> mobility(EMU,Metric) # C⋅g⋅abC⁻¹⋅kg
$(mobility(EMU,Metric))

julia> mobility(EMU,ESU) # statC⋅abC⁻¹
$(mobility(EMU,ESU))

julia> mobility(ESU,Metric) # C⋅g⋅statC⁻¹⋅kg
$(mobility(ESU,Metric))

julia> mobility(Metric,SI2019) # C⋅C⁻¹
$(mobility(Metric,SI2019))
```
"""
@pure mobility(U::UnitSystem,S::UnitSystem) = unit(charge(U,S)*time(U,S)/mass(U,S))

"""
$(convertext(:reluctance,"rationalization(U,S)*lorentz(U,S)^2/inductance(U,S)"))

Magnetic `reluctance` or magnetic resistance (H⁻¹, Gb⋅Mx⁻¹), unit conversion factor.

```Julia
julia> reluctance(EMU,Metric) # abH⋅H⁻¹
$(reluctance(EMU,Metric))

julia> reluctance(ESU,Metric) # statH⋅H⁻¹
$(reluctance(ESU,Metric))

julia> reluctance(Metric,SI2019) # H⋅H⁻¹
$(reluctance(Metric,SI2019))
```
""" # reciprocal: permeance -- different concept from inductace but same units
@pure reluctance(U::UnitSystem,S::UnitSystem) = unit(rationalization(U,S)*lorentz(U,S)^2/inductance(U,S))

"""
$(convertext(:vectorpotential,"magneticflux(U,S)/length(U,S)"))

Magnetic `vectorpotential` (Wb⋅m⁻¹), unit conversion factor.

```julia
julia> vectorpotential(EMU,Metric) # Wb⋅cm⋅Mx⁻¹⋅m⁻¹
$(vectorpotential(EMU,Metric))

julia> vectorpotential(ESU,Metric) # Wb⋅cm⋅statWb⁻¹⋅m⁻¹
$(vectorpotential(ESU,Metric))

julia> vectorpotential(Metric,SI2019) # Wb⋅Wb⁻¹
$(vectorpotential(Metric,SI2019))
```
"""
@pure vectorpotential(U::UnitSystem,S::UnitSystem) = unit(magneticflux(U,S)/length(U,S))

"""
$(convertext(:magneticmoment,"magneticflux(U,S)*length(U,S)"))

Amount of `magneticmoment` or `magneticflux` by `length` (Wb⋅m), unit conversion factor.

```Julia
julia> magneticmoment(EMU,Metric) # Wb⋅m⋅Mx⁻¹⋅cm⁻¹
$(magneticmoment(EMU,Metric))

julia> magneticmoment(ESU,Metric) # Wb⋅m⋅statWb⁻¹⋅cm⁻¹
$(magneticmoment(ESU,Metric))

julia> magneticmoment(Metric,SI2019) # Wb⋅Wb⁻¹
$(magneticmoment(Metric,SI2019))
```
"""
@pure magneticmoment(U::UnitSystem,S::UnitSystem) = unit(magneticflux(U,S)*length(U,S))

"""
$(convertext(:rigidity,"magneticfluxdensity(U,S)*length(U,S)"))

Electromagnetic `rigidity` (T⋅m), unit conversion factor.

```Julia
julia> rigidity(EMU,Metric) # T⋅m⋅G⁻¹⋅cm⁻¹
$(rigidity(EMU,Metric))

julia> rigidity(EMU,ESU) # statT⋅G⁻¹
$(rigidity(EMU,ESU))

julia> rigidity(Metric,SI2019) # T⋅T⁻¹
$(rigidity(Metric,SI2019))
```
"""
@pure rigidity(U::UnitSystem,S::UnitSystem) = unit(magneticfluxdensity(U,S)*length(U,S))

"""
$(convertext(:susceptibility,"1/rationalization(U,S)"))

Magnetic/electric volume `susceptibility` (dimensionless), unit conversion factor.

```Julia
julia> susceptibility(EMU,Metric)
$(susceptibility(EMU,Metric))

julia> susceptibility(ESU,Metric)
$(susceptibility(ESU,Metric))

julia> susceptibility(Metric,SI2019)
$(susceptibility(Metric,SI2019))
```
""" # magneticdipolemoment(U,S)/magneticfield(U,S)/volume(U,S)
@pure susceptibility(U::UnitSystem,S::UnitSystem) = unit(rationalization(S,U))

# WARNING unchecked: rigidity, magneticmoment, vectorpotential, mobility, linearchargedensity, exposure

# CGS extra: polarizability, permeance, magnetic-current? + density, magneticresistance

"""
$(convertext(:electricflux,"electricpotential(U,S)*length(U,S)"))

Amount of `electricflux` or `electricpotential` by `length` (V⋅m), unit conversion factor.

```Julia
julia> electricflux(EMU,Metric) # V⋅m⋅abV⁻¹⋅cm⁻¹
$(electricflux(EMU,Metric))

julia> electricflux(EMU,ESU) # statV⋅abV⁻¹
$(electricflux(EMU,ESU))

julia> electricflux(ESU,Metric) # V⋅m⋅statV⁻¹⋅cm⁻¹
$(electricflux(ESU,Metric))

julia> electricflux(Metric,SI2019) # V⋅V⁻¹
$(electricflux(Metric,SI2019))
```
"""
@pure electricflux(U::UnitSystem,S::UnitSystem) = unit(electricpotential(U,S)*length(U,S))

"""
$(convertext(:electricdipolemoment,"charge(U,S)*length(U,S)"))

Electric dipole moment or `electricdipolemoment` (C⋅m), unit conversion factor.

```Julia
julia> electricdipolemoment(EMU,Metric) # C⋅m⋅abC⁻¹⋅cm⁻¹
$(electricdipolemoment(EMU,Metric))

julia> electricdipolemoment(ESU,Metric) # C⋅m⋅statC⁻¹⋅cm⁼¹
$(electricdipolemoment(ESU,Metric))

julia> electricdipolemoment(Metric,Gauss)/1e-18 # D⋅C⁻¹⋅m⁻¹
$(electricdipolemoment(Metric,Gauss)/1e-18)

julia> electricdipolemoment(Metric,SI2019) # C⋅C⁻¹
$(electricdipolemoment(Metric,SI2019))
```
"""
@pure electricdipolemoment(U::UnitSystem,S::UnitSystem) = unit(charge(U,S)*length(U,S))

"""
$(convertext(:magneticpotential,"magneticflux(U,S)*reluctance(U,S)"))

Magnetomotive force or `magneticpotential` (A, Gb), unit conversion factor.

```Julia
julia> magneticpotential(EMU,Metric) # A⋅Gb⁻¹
$(magneticpotential(EMU,Metric))

julia> magneticpotential(Metric,SI2019) # A⋅A⁻¹
$(magneticpotential(Metric,SI2019))
```
"""
@pure magneticpotential(U::UnitSystem,S::UnitSystem) = unit(magneticflux(U,S)*reluctance(U,S))

"""
$(convertext(:polestrength,"magneticdipolemoment(U,S)/length(U,S)"))

Magnetic `polestrength` is analogous to `charge` (A⋅m⁻¹), unit conversion factor.

```Julia
julia> polestrength(EMU,Metric) # A⋅m⁻¹⋅pole⁻¹
$(polestrength(EMU,Metric))

julia> polestrength(Metric,SI2019) # A⋅A⁻¹⋅
$(polestrength(Metric,SI2019))
```
"""
@pure polestrength(U::UnitSystem,S::UnitSystem) = unit(magneticdipolemoment(U,S)/length(U,S))

"""
$(convertext(:permeance,"1/reluctance(U,S)"))

Magnetic `permeance` or magnetic conductance (H, Mx⋅Gb⁻¹), unit conversion factor.

```Julia
julia> permeance(EMU,Metric) # abH⋅H⁻¹
$(permeance(EMU,Metric))

julia> permeance(ESU,Metric) # statH⋅H⁻¹
$(permeance(ESU,Metric))

julia> permeance(Metric,SI2019) # H⋅H⁻¹
$(permeance(Metric,SI2019))
```
"""
@pure permeance(U::UnitSystem,S::UnitSystem) = reluctance(S,U)

"""
$(convertext(:specificsusceptibility,"susceptibility(U,S)/density(U,S)"))

Magnetic/electric mass specific `susceptibility` (m³⋅kg⁻¹), unit conversion factor.

```Julia
julia> specificsusceptibility(EMU,Metric) # m³⋅g⋅kg⁻¹⋅cm⁻³
$(specificsusceptibility(EMU,Metric))

julia> specificsusceptibility(ESU,Metric) # m³⋅g⋅kg⁻¹⋅cm⁻³
$(specificsusceptibility(ESU,Metric))

julia> specificsusceptibility(Metric,SI2019) # m³⋅kg⋅kg⁻¹⋅m⁻³
$(specificsusceptibility(Metric,SI2019))
```
"""
@pure specificsusceptibility(U::UnitSystem,S::UnitSystem) = unit(magneticdipolemoment(U,S)/magneticfield(U,S)/mass(U,S))

"""
$(convertext(:magnetizability,"magneticmoment(U,S)/magneticfluxdensity(U,S)"))

Quantity of `magneticmoment` per `magneticfluxdensity` (m⁻¹), unit conversion factor.

```Julia
julia> magnetizability(EMU,Metric) # cm⋅m⁻¹
$(magnetizability(EMU,Metric))

julia> magnetizability(ESU,Metric) # cm⋅m⁻¹
$(magnetizability(ESU,Metric))

julia> magnetizability(Metric,SI2019) # m⋅m⁻¹
$(magnetizability(Metric,SI2019))
```
"""
@pure magnetizability(U::UnitSystem,S::UnitSystem) = unit(magneticmoment(U,S)/magneticfluxdensity(U,S))

"""
$(convertext(:electricpolarizability,"electricdipolemoment(U,S)/electricfield(U,S)"))

Polarizability or `electricdipolemoment` per `electricfield` (C⋅m²⋅V⁻¹), unit conversion factor.

```Julia
julia> electricpolarizability(EMU,Metric) # C⋅m²⋅abV⋅abC⁻¹⋅cm⁻²⋅V⁻¹
$(electricpolarizability(EMU,Metric))

julia> electricpolarizability(ESU,Metric) # C⋅m²⋅statV⋅statC⁻¹⋅cm⁼²⋅V⁻¹
$(electricpolarizability(ESU,Metric))

julia> electricpolarizability(Metric,Gauss) # D⋅cm²⋅V⁻¹⋅C⁻¹⋅m⁻²⋅abV⁻¹
$(electricpolarizability(Metric,Gauss))

julia> electricpolarizability(Metric,SI2019) # C⋅V⋅C⁻¹⋅V⁻¹
$(electricpolarizability(Metric,SI2019))
```
"""
@pure electricpolarizability(U::UnitSystem,S::UnitSystem) = unit(electricdipolemoment(U,S)/electricfield(U,S))

"""
$(convertext(:magneticpolarizability,"magneticdipolemoment(U,S)/magneticfield(U,S)"))

Polarizability or `magneticdipolemoment` per `magneticfield` (m³), unit conversion factor.

```Julia
julia> electricpolarizability(EMU,Metric) # m³⋅cm⁻³
$(electricpolarizability(EMU,Metric))

julia> electricpolarizability(ESU,Metric) # m³⋅cm⁼³
$(electricpolarizability(ESU,Metric))

julia> electricpolarizability(Metric,Gauss) # cm³⋅m⁻³
$(electricpolarizability(Metric,Gauss))

julia> electricpolarizability(Metric,SI2019)
$(electricpolarizability(Metric,SI2019))
```
"""
@pure magneticpolarizability(U::UnitSystem,S::UnitSystem) = unit(magneticdipolemoment(U,S)/magneticfield(U,S))

"""
$(convertext(:magnetization,"magneticmoment(U,S)/volume(U,S)"))

Amount of `magneticmoment` per `volume` (Wb⋅m⁻²), unit conversion factor.

```Julia
julia> magnetization(EMU,Metric) # Wb⋅cm²⋅Mx⁻¹⋅m⁻²
$(magnetization(EMU,Metric))

julia> magnetization(ESU,Metric) # Wb⋅cm²⋅statWb⁻¹⋅m⁻²
$(magnetization(ESU,Metric))

julia> magnetization(Metric,SI2019) # Wb⋅Wb⁻¹
$(magnetization(Metric,SI2019))
```
"""
@pure magnetization(U::UnitSystem,S::UnitSystem) = unit(magneticmoment(U,S)/volume(U,S))

# specificmagnetization, mass magnetization = 1?
# magneticfluxdensity(Metric,EMU)/density(Metric,EMU)
"""
$(convertext(:specificmagnetization,"magneticmoment(U,S)/mass(U,S)"))

Amount of `magneticmoment` per `mass` (Wb⋅m⋅kg⁻¹), unit conversion factor.

```Julia
julia> specificmagnetization(EMU,Metric) # Wb⋅m⋅g⋅Mx⁻¹⋅cm⁻¹⋅kg⁻¹
$(specificmagnetization(EMU,Metric))

julia> specificmagnetization(ESU,Metric) # Wb⋅m⋅g⋅statWb⁻¹⋅cm⁻¹⋅kg⁻¹
$(specificmagnetization(ESU,Metric))

julia> specificmagnetization(Metric,SI2019) # Wb⋅Wb⁻¹
$(specificmagnetization(Metric,SI2019))
```
"""
@pure specificmagnetization(U::UnitSystem,S::UnitSystem) = unit(magneticmoment(S,U)/mass(S,U))

"""
$(convertext(:demagnetizingfactor,"1/susceptibility(U,S)"))

Quantitiy of `demagnetizingfactor` (dimensionless), unit conversion factor.

```Julia
julia> demagnetizingfactor(EMU,Metric)
$(demagnetizingfactor(EMU,Metric))

julia> demagnetizingfactor(ESU,Metric)
$(demagnetizingfactor(ESU,Metric))

julia> demagnetizingfactor(Metric,SI2019)
$(demagnetizingfactor(Metric,SI2019))
```
"""
@pure demagnetizingfactor(U::UnitSystem,S::UnitSystem) = unit(rationalization(U,S))

# Gyrator-capacitor model alternative:
#@pure magneticfluxrate(U::UnitSystem,S::UnitSystem) = unit(magneticflux(U,S)/time(U,S))
#@pure magneticfluxratedensity(U::UnitSystem,S::UnitSystem) = unit(magneticfluxrate(U,S)/area(U,S))
#@pure magneticresistance(U::UnitSystem,S::UnitSystem) = unit(magneticpotential(U,S)/magneticfluxrate(U,S)) # not reluctance


text:
const testdata = ["127.0.0.1", "127.0.0.1:80", "::1", "[::1]:80",
                  "2605:2700:0:3::4713:93e3", "[2605:2700:0:3::4713:93e3]:80",
                  "::ffff:192.168.173.22", "[::ffff:192.168.173.22]:80",
                  "1::", "[1::]:80", "::", "[::]:80"]

maybev4(ip) = search(ip, '.') > 0 && length(matchall(r":", ip)) < 2
maybev6(ip) = length(matchall(r":", ip)) > 1

function parseip(ip)
    if (mat = match(r"^\[([:.\da-fA-F]+)\]:(\d+)$", ip))!= nothing ||
       (mat = match(r"^([\d.]+)[:/](\d+)$", ip)) != nothing
        port = mat.captures[2]
        ip = mat.captures[1]
    else
        port = "none"
    end
    if maybev4(ip)
        println("Processing ip v4 $ip")
        iphex = hex(Int(Base.IPv4(ip)))
        addresspace = "IPv4"
    elseif maybev6(ip)
        println("Processing ip v6 $ip")
        iphex = hex(UInt128(Base.IPv6(ip)))
        addresspace = "IPv6"
    else
        throw("Bad IP address argument $ip")
    end
    iphex, addresspace, port
end

for ip in testdata
    hx, add, por = parseip(ip)
    println("For input $ip, IP in hex is $hx, address space $add, port $por.")
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module Authorization

using HTTP
using SHA

function login()
    win = Sys.iswindows()

    un = prompt("username/email: ")
    if sha2_512(prompt("password: ")) == USERS[un][:sha_pass]
        println("continue login")
    else
        println("invalid credentials, please try again.")
    end
end

function register_credentials()
    win = Sys.iswindows()

    username = Base.prompt("Enter a username")
    email = Base.prompt("Enter an email address")
    pass = Base.getpass("Enter a password")
    sha_pass = sha2_512(pass)

    Base.shred!(pass)

    while true
        break
        #=
        if length(read(user_pass, String)) > 7
            break
        else
            println("Invalid password. Passwords must contain at least 8 characters. Please try again.")
        end
        =#
    end

    Dict(
        :username=>username,
        :email=>email,
        :sha_pass=>sha_pass,
        )
end

end # module

max_stars_repo_name: dustyirwin/PlaymatSimulator

max_issues_repo_name: dustyirwin/PlaymatSimulator

max_forks_repo_name: dustyirwin/PlaymatSimulator


text:
module SunlightAPIs

import Base.print

using JSON
using HttpCommon
using Requests

# Entity types (organization, individual, politician)
export Entity, Org, Ind, Pol

# Influence Explorer API
export entity_search,
       entity_info,
       top_politicians,
       top_contributors,
       top_industries,
       unknown_industries,
       top_sectors,
       local_breakdown,
       contributor_breakdown,
       fec_summary,
       fec_independent_expenditures,
       top_individuals,
       top_recipient_organizations,
       top_recipient_politicians,
       party_breakdown,
       lobbying_registrants,
       lobbying_clients,
       lobbying_issues,
       top_organizations,
       top_recipients,
       pac_recipients,
       state_federal_breakdown,
       lobbing_registrants,
       bills,
       lobbyists,
       registrant_clients,
       registrant_issues,
       registrant_bills,
       registrant_lobbyists,
       mentions_in_regulations,
       regulatory_comment_submissions,
       faca_memberships,
       top_industries,
       top_industry_organizations,
       get_metadata_overview,
       get_state_metadata,
       bill_search,
       bill_detail,
       legislator_search,
       legislator_detail,
       geo_lookup,
       committee_search,
       committee_detail,
       event_search,
       event_detail,
       district_search,
       district_boundary_lookup

include("utils.jl")
include("influence-explorer.jl")
include("open-states.jl")

end


text:
module TestLSM

using Base.Test

using Bio.Util

@testset "LabelledSquareMatrices" begin

    integral_matrix = """\
    \tA\tB\tC
    A\t1\t2\t4
    B\t2\t1\t2
    C\t4\t2\t1
    """
    integral_matrix_expt = Float64[1 2 4; 2 1 2; 4 2 1;]
    commented_matrix = """\
    # This is a comment that gets skipped. The next blank line should also get
    # skipped.

    \tA\tB\tC
    A\t1\t2\t4
    B\t2\t1\t2
    C\t4\t2\t1
    """
    spacy_matrix = """\
    # This matrix has heaps of extra spaces in the road, to try and trip up the
    # parser.

       A  B  C
    A 1 2          4
    B  2  1  2
    C  4  2  1
    """

    float_matrix = """\
    \tA\tB\tC
    A\t1.0\t2.5\t4.5
    B\t2.5\t1.0\t2.5
    C\t4.5\t2.5\t1.0
    """
    float_matrix_expt = Float64[1. 2.5 4.5; 2.5 1. 2.5; 4.5 2.5 1.;]

    good_matricies = [(float_matrix, float_matrix_expt),
                      (integral_matrix, integral_matrix_expt),
                      (commented_matrix, integral_matrix_expt),
                      (spacy_matrix, integral_matrix_expt),
                     ]

    bad_labels = """\
    \tA\tB\tC
    X\t1\t2\t4
    Y\t2\t1\t2
    Z\t4\t2\t1
    """

    not_square = """\
    \tA\tB\tC
    A\t1\t2\t4
    """

    missing_vals = """\
    \tA\tB\tC
    A\t1\t2\t4
    B\t2\t1\t2
    C\t4
    """
    bad_matricies = [bad_labels,
                     not_square,
                     missing_vals,
                     ]

    file_matricies = [("ints.tab", integral_matrix_expt),
                      ("commented.tab", integral_matrix_expt),
                     ]

    # Maps filename to labels
    NUC44 = ["A", "T", "G", "C", "S", "W", "R", "Y", "K", "M", "B", "V", "H", "D", "N"]
    BLOSUM = ["A", "R", "N", "D", "C", "Q", "E", "G", "H", "I", "L", "K",
              "M", "F", "P", "S", "T", "W", "Y", "V", "B", "Z", "X", "*"]
    subst_matricies = [("BLOSUM62", BLOSUM),
                       ("NUC.4.4", NUC44),
                     ]
    all_subst_matrices = ["BLOSUM45", "BLOSUM50", "BLOSUM62", "BLOSUM80",
                          "BLOSUM90", "NUC.4.4", "PAM250", "PAM30", "PAM70"]

    @testset "types" for T in [Int64, Float64]
        m, l = readlsm(T, IOBuffer(integral_matrix))
        expt_m = Array{T}(integral_matrix_expt)
        @test m == expt_m
        @test eltype(m) == T
        @test l == ["A", "B", "C"]
    end

    @testset "good_matricies" for (mat, expected) in good_matricies
        m, l = readlsm(IOBuffer(mat))
        @test m == expected
        @test l == ["A", "B", "C"]
    end

    @testset "bad_matricies" for mat in bad_matricies
        @test_throws Exception readlsm(IOBuffer(mat))
    end

    @testset "files" for (fname, expected) in file_matricies
        fname = joinpath(dirname(@__FILE__), "..", "BioFmtSpecimens", "LSM", fname)
        m, l = readlsm(fname)
        @test m == expected
        @test l == ["A", "B", "C"]
    end

    @testset "subst_matrices" for (fname, labels) in subst_matricies
        fname = joinpath(dirname(@__FILE__), "..", "BioFmtSpecimens", "LSM", fname)
        N = length(labels)
        m, l = readlsm(Int64, fname)
        @test size(m) == (N, N)
        @test l == labels
    end

    @testset "all_subst_matrices" for fname in all_subst_matrices
        # Check that we can parse all matrices under the alignm module
        fname = joinpath(dirname(@__FILE__), "..", "..", "src", "align", "data", "submat", fname)
        m, l = readlsm(Int64, fname)
        @test eltype(m) === Int64
    end
end

end


text:
include("./exportData.jl")
using .exportData, DifferentialEquations, DiffEqFlux, PyPlot, DelimitedFiles

const N = Int(round(3/0.015 * sqrt(2)))
const D1= 0.0010/(0.015*0.015)
const D2= 0.0002/(0.015*0.015)

# define fenkar system
function H(x;k=100.0)
	return 0.5*(1.0 + tanh(k*x))
end
function fenkar!(dx, x, p, t)

	# parameters
	@views tsi,tv1m,tv2m,tvp,twm,twp,td,to,tr,xk,uc,uv,ucsi,θ,t0,TI,IA = p[1:17]

	@inbounds for n=1:N

		# fenkar dynamics
		dx[n,1] =-(x[n,1]*H(uc-x[n,1])/to + H(x[n,1]-uc)/tr - x[n,2]*H(x[n,1]-uc)*(1.0-x[n,1])*(x[n,1]-uc)/td - x[n,3]*H(x[n,1]-ucsi;k=xk)/tsi)
		dx[n,2] = H(uc-x[n,1])*(1.0-x[n,2])/(tv1m*H(uv-x[n,1]) + tv2m*H(x[n,1]-uv)) - H(x[n,1]-uc)*x[n,2]/tvp
		dx[n,3] = H(uc-x[n,1])*(1.0-x[n,3])/twm - H(x[n,1]-uc)*x[n,3]/twp

		# no-flux with angular dependence on two diffusions
		if n==1
			dx[n,1] = dx[n,1] + (D1*cos(θ)^2 + D2*sin(θ)^2)*(2*x[2,1] - 2*x[n,1])
		elseif n==N
			dx[n,1] = dx[n,1] + (D1*cos(θ)^2 + D2*sin(θ)^2)*(2*x[N-1,1] - 2*x[n,1])
		else
			dx[n,1] = dx[n,1] + (D1*cos(θ)^2 + D2*sin(θ)^2)*(x[n-1,1] + x[n+1,1] - 2*x[n,1])
		end
	end

	@inbounds for n=1:10
		dx[n,1] = dx[n,1] + H(t-t0;k=1.0).*IA*sin(pi*(t-t0)/TI)^500
	end

	return nothing
end

# fenkar system parameters in p[1:13] (BR; from https://doi.org/10.1063/1.166311)
# and then p[14:14] = [fiber_angle]
# and then p[15:17] = [offset_time, stimulus_period, stimulus_amplitude]
   p = [29.0, 19.6, 1250.0, 3.33, 41.0, 870.0, 0.25, 12.5, 33.3, 10.0, 0.13, 0.04, 0.85, 0.0,  1.0,   125.0, 0.1]

tspan = (0.0, 3000.0)
u0 = zeros(Float64,N,3)

#solve and take data
prob = ODEProblem(fenkar!, u0, tspan, p)
sol = solve(prob, Tsit5(); saveat=tspan[begin]:5.0:tspan[end])

x = collect(0:(N-1)).*0.015;
plt.pcolormesh(sol.t, x, sol[:,1,:], rasterized=true, snap=true, vmin=0.0, vmax=1.0)
plt.colorbar()
plt.savefig("./pde/flux/baseline.pdf")
plt.close()

# getExpData
ind = 78
t, data = getExpData(ind, tInds=1:1000)

# remake problem
prob = remake(prob, tspan=(0.0, t[end]), p=p)

# use the save_idxs to select the mid-point of the u-variable
sol = solve(prob, Tsit5(); saveat=t, save_idxs=Int(round(N//2)):Int(round(N//2)))

function loss(p)
  sol = solve(prob, Tsit5(), p=p, saveat = t, save_idxs=Int(round(N//2)):Int(round(N//2)), sensealg=InterpolatingAdjoint(autojacvec=ReverseDiffVJP(true))) #sensealg=ForwardDiffSensitivity()
  if any((s.retcode != :Success for s in sol))
     loss = Inf
  else
     loss = sum(abs2, sol.-data)
  end
  return loss, sol
end
l0, pred = loss(p)

callback = function (p, l, pred; save=false)
  @show l
  plt.cla()
  plt.plot(sol.t, sol[1,:], "-", linewidth=1)
  plt.plot(t, data[1,:], "--")
  plt.plot(pred.t, pred[1,:], ":")
  # Tell sciml_train to not halt the optimization. If return true, then
  # optimization stops.
  if save==true
  	plt.savefig("./pde/flux/$(ind)_fittings.pdf")
  end
  return false
end

fig = plt.figure(figsize=(4,2))
callback(p, l0, pred; save=true)

result = DiffEqFlux.sciml_train(loss, p, maxiters=300, cb = callback)

l0, pred = loss(result)
callback(result, l, pred; save=true)
plt.close("all")

open("./pde/flux/$(ind)_params.txt", "w") do io
	writedlm(io, ["tsi" "tv1m" "tv2m" "tvp" "twm" "twp" "td" "to" "tr" "xk" "uc" "uv" "ucsi" "θ" "t0" "TI" "IA"])
	writedlm(io, "\n")
	writedlm(io, transpose(p))
	writedlm(io, "\n")
	writedlm(io, transpose(result))
end


text:
using CineFiles
using ImageCore
using Images: load
using TiffImages
using Test
using Glob
using XMLDict: parse_xml
using FixedPointNumbers: nbitsfrac

# In order to test against arbitary cine files, export tif image of first frame with gamma = 1
# and xml file using PCC. Ensure the cine file, tif image and xml file have the same name and
# place in either data or proprietary_data folders.

struct cine_test_files
    cine_path::String
    xml_path::String
    first_frame_path::String
end

function cine_test_files(cine_path::String)
    base_name = splitext(cine_path)[1]
    xml_path = base_name * ".xml"
    first_frame_path = base_name * ".tif"
    paths = (cine_path, xml_path, first_frame_path)
    exists = isfile.(paths)
    if all(exists)
        cine_test_files(paths...)
    else
        @error "Missing file(s) $(paths[[.~exists...]])"
    end
end

function compare_images(cine_file, tiff_image, frame_no = 1)
    cine_frame = cine_file[frame_no]
    cine_frame_gamma = (cine_frame) .+ cine_file.header.setup.fOffset
    # Have not included gamma adjustment yet
    cine_frame_gamma[cine_frame_gamma.<0] .= 0
    return all(isapprox.(cine_frame_gamma, tiff_image; atol = 0.01))
end

cine_file_paths = cine_test_files.(glob("*.cine", "data"))
append!(cine_file_paths, cine_test_files.(glob("*.cine", "proprietary_data")))

@testset "CineFiles.jl" begin
    @testset "Dependency internals" begin
        # nbitsfrac is unexported, so we'll test to make sure it doesn't break
        @test nbitsfrac(N0f8)  == 8
        @test nbitsfrac(N4f12) == 12
        @test nbitsfrac(N0f16) == 16
    end

    @time cf8 = CineFile(joinpath("data", "8bpp.cine"))
    @testset "8 bit grayscale" begin
        @test length(cf8) == 202
        @test size(cf8[1]) == (16, 128)
        @test eltype(cf8) == Matrix{Gray{N0f8}}
    end

    cf12 = CineFile(joinpath("data", "12bpp.cine"))
    @testset "12 bit grayscale" begin
        @test length(cf12) == 202
        @test size(cf12[1]) == (16, 128)
        @test eltype(cf12) == Matrix{Gray{N4f12}}
    end

    @test cf12[1] ≈ cf8[1]

    first_load_time = (@elapsed cf12[2])
    cached_load_time = (@elapsed cf12[2])
    @show first_load_time, cached_load_time
    @test first_load_time > cached_load_time

    @testset "Tiff compare first frame" begin
        for file in cine_file_paths
            @info file.cine_path
            cf = CineFile(file.cine_path)
            tiff_img = load(file.first_frame_path)
            @test compare_images(cf, tiff_img)
        end
    end

    @testset "XML compare camera version" begin
        for file in cine_file_paths
            @info file.cine_path
            cf = CineFile(file.cine_path)
            xml_data = parse_xml(read(file.xml_path, String))
            @test Int(cf.header.setup.CameraVersion) ==
                  parse(Int, xml_data["CameraSetup"]["CameraVersion"])
        end
    end
end