text:
# Hyper parameters
ωs = [1.0, 0.0, 0.0]
ωp = 1e-7
# PO₄ mean and variance of observations fom WOA18
WOA = WorldOceanAtlasTools
μDIPobs3D, σ²DIPobs3D = WOA.fit_to_grid(grid, 2018, "phosphate", "annual", "1°", "an")
μDIPobs, σ²DIPobs = μDIPobs3D[iwet], σ²DIPobs3D[iwet]
μx = (μDIPobs, missing, missing)
σ²x = (σ²DIPobs, missing, missing)
# generate mismatch functions
f   =   generate_objective(ωs, μx, σ²x, v, ωp, mean_obs(p₀), variance_obs(p₀))
∇ₓf = generate_∇ₓobjective(ωs, μx, σ²x, v, ωp, mean_obs(p₀), variance_obs(p₀))
∇ₚf = generate_∇ₚobjective(ωs, μx, σ²x, v, ωp, mean_obs(p₀), variance_obs(p₀))

#===========================================
Tests
===========================================#

@testset "Mismatch / cost functions" begin
    nt = length(T_all)
    n = nt * nb
    x₀ = p₀.xgeo * ones(n)
    f₀ = f(x₀, p₀)
    ∇ₓf₀ = ∇ₓf(x₀, p₀)
    ∇ₚf₀ = ∇ₚf(x₀, p₀)
    @testset "Mismatch function" begin
        @test f₀ isa Float64
        @test f₀ > 0
    end
    @testset "Derivative w.r.t the state, x" begin
        @test size(∇ₓf₀) == (1, n)
    end
    @testset "Derivative w.r.t the parameters, p" begin
        @test size(∇ₚf₀) == (1, m)
    end
    @testset "Lognormal <-> Normal" begin
        mm, vv = rand(10), rand(10)
        @test AIBECS.LNm(AIBECS.LNμ(mm, vv), AIBECS.LNσ²(mm, vv)) ≈ mm
        @test AIBECS.LNv(AIBECS.LNμ(mm, vv), AIBECS.LNσ²(mm, vv)) ≈ vv
    end
end


text:
using Test, MPI

using CLIMA
using CLIMA.MPIStateArrays

CLIMA.init()
const ArrayType = CLIMA.array_type()
const mpicomm = MPI.COMM_WORLD

@testset "MPIStateArray basics" begin
  Q = MPIStateArray{Float32}(mpicomm, ArrayType, 4, 6, 8)

  @test eltype(Q) == Float32
  @test size(Q) == (4, 6, 8)

  fillval = 0.5f0
  fill!(Q, fillval)

  CLIMA.gpu_allowscalar(true)
  @test Q[1] ==  fillval
  @test Q[2, 3, 4] == fillval
  @test Q[end] == fillval

  @test Array(Q) == fill(fillval, 4, 6, 8)

  Q[2, 3, 4] = 2fillval
  @test Q[2, 3, 4] != fillval
  CLIMA.gpu_allowscalar(false)

  Qp = copy(Q)

  @test typeof(Qp) == typeof(Q)
  @test eltype(Qp) == eltype(Q)
  @test size(Qp) == size(Q)
  @test Array(Qp) == Array(Q)

  Qp = similar(Q)

  @test typeof(Qp) == typeof(Q)
  @test eltype(Qp) == eltype(Q)
  @test size(Qp) == size(Q)

  copyto!(Qp, Q)
  @test Array(Qp) == Array(Q)
end

ext: jl

max_stars_repo_name: favba/CLIMA

max_issues_repo_name: favba/CLIMA

max_forks_repo_name: favba/CLIMA


text:
module GadgetReader
using HDF5
export read_snap

function read_snap(filename)

    T=Float64
    header::Dict = h5readattr(filename, "/Header")
    boxsize::T = header["BoxSize"]
    time::T    = header["Time"]

    N_gas::Int64 = header["NumPart_ThisFile"][1]

    pos_gas::Matrix{T} = h5read(filename, "PartType0/Coordinates");
    vel_gas::Matrix{T} = h5read(filename, "PartType0/Velocities");
    rho::Vector{T}     = h5read(filename, "PartType0/Density");
    u::Vector{T}       = h5read(filename, "PartType0/InternalEnergy");
    m_gas::Vector{T}   = h5read(filename, "PartType0/Masses");
    hsml::Vector{T}    = h5read(filename, "PartType0/SmoothingLength");
    #scal::Vector{T}       = h5read(filename, "PartType0/PassiveScalarField");

    id_gas::Vector{Int64} = h5read(filename, "PartType0/ParticleIDs");

    #abund::Matrix{T} = h5read(filename, "PartType0/ChemicalAbundancesSG");
    #fH2ss::Vector{T}    = h5read(filename, "PartType0/ShieldingFactorH2");
    #fH2dust::Vector{T}    = h5read(filename, "PartType0/ShieldingFactorDust");

    #col_tot::Matrix{T}    = h5read(filename, "PartType0/TreecolColumnDensitiesAll");
    #col_H2::Matrix{T}    = h5read(filename, "PartType0/TreecolColumnDensitiesH2");
    #col_CO::Matrix{T}    = h5read(filename, "PartType0/TreecolColumnDensitiesCO");

    #Tdust::Vector{T}    = h5read(filename, "PartType0/DustTemperature");

    #N_star::Int64 = header["NumPart_ThisFile"][5]

    #pos_star::Matrix{T} = h5read(filename, "PartType4/Coordinates");
    #vel_star::Matrix{T} = h5read(filename, "PartType4/Velocities");
    #m_star::Vector{T}   = h5read(filename, "PartType4/Masses");
    #sftime::Vector{T}      = h5read(filename, "PartType4/StellarFormationTime");
    #id_star::Vector{Int64} = h5read(filename, "PartType4/ParticleIDs");

    return N_gas, pos_gas, vel_gas, rho, u, m_gas, hsml, id_gas,
        #abund, fH2ss, fH2dust, col_tot, col_H2, col_CO, Tdust,
        #N_star, pos_star, vel_star, m_star, sftime, id_star,
        boxsize, time
end

end


text:
using PIMD
using Documenter

makedocs(;
    modules=[PIMD],
    authors="Songchen Tan",
    repo="https://github.com/tansongchen/PIMD.jl/blob/{commit}{path}#L{line}",
    sitename="PIMD.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://tansongchen.com/PIMD.jl",
        assets=String[],
    ),
    pages=[
        "Introduction" => "index.md",
        "API" => "api.md"
    ],
)

deploydocs(;
    repo="github.com/tansongchen/PIMD.jl",
)


text:
module TuringTools
# Tools for Turing

using DataFrames
using Dates
using Mustache

using ..Utils

export @tbenchmark,
    @tbenchmark_expr,
    build_log_data,
    print_log

macro tbenchmark(alg, model, data)
    model = :(($model isa String ? eval(Meta.parse($model)) : $model))
    model_dfn = (data isa Expr && data.head == :tuple) ?
        :(model_f = $model($(data)...)) : model_f = :(model_f = $model($data))
    esc(quote
        $model_dfn
        chain, t_elapsed, mem, gctime, memallocs  = @timed sample(model_f, $alg)
        Dict(
            "engine" => $(string(alg)),
            "time" => t_elapsed,
            "memory" => mem,
        )
        end)
end

macro tbenchmark_expr(engine, expr)
    quote
        chain, t_elapsed, mem, gctime, memallocs  = @timed $(esc(expr))
        Dict(
            "engine" => $(string(engine)),
            "time" => t_elapsed,
            "memory" => mem,
        )
    end
end

const tmpl_log_string = """
/=======================================================================
| Benchmark Result
|-----------------------------------------------------------------------
| Overview
|-----------------------------------------------------------------------
| Inference Engine  : {{{ engine }}}
| Time Used (s)     : {{{ time }}}
{{#time_stan}}
|   -> time by Stan : {{{ time_stan }}}
{{/time_stan}}
| Mem Alloc (bytes) : {{{ memory }}}
{{#turing}}
|-----------------------------------------------------------------------
| Turing Inference Result
|-----------------------------------------------------------------------
{{/turing}}
{{#turing_items}}
| >> {{{ name }}} <<
| mean = {{{ mean }}}
{{#analytic}}
|   -> analytic = {{{ analytic }}}
{{/analytic}}
{{#anal_diff}}
|        |--*-->  diff = {{{ anal_diff }}}
{{/anal_diff}}
{{#stan}}
|   -> Stan     = {{{ stan }}}
{{/stan}}
{{#stan_diff}}
|        |--*--> diff = {{{ anal_diff }}}
{{/stan_diff}}
{{/turing_items}}
{{#turing_strings}}
| {{.}}
{{/turing_strings}}
{{#note}}
|-----------------------------------------------------------------------
| Note:
|   {{{ note }}}
{{/note}}
\\=======================================================================
"""

function stringify_log(logd::Dict, monitor=[])
    data = Dict{Any, Any}(
        # "name" => logd["name"],
        "engine" => logd["engine"],
        "time" => logd["time"],
        "memory" => logd["memory"],
    )
    haskey(logd, "time_stan") && (data["time_stan"] = logd["time_stan"])
    haskey(logd, "note") && (data["note"] = logd["note"])

    if haskey(logd, "turing")
        data["turing"] = true
        data["turing_items"] = []
        data["turing_strings"] = []
        if isa(logd["turing"], String)
            push!(data["turing_strings"], logd["turing"])
        else
            for (v, m) = logd["turing"]
                (!isempty(monitor) && !(v in monitor)) && continue
                item = Dict{Any, Any}("name" => v)
                item["mean"] = round.(m, digits=3)
                if haskey(logd, "analytic") && haskey(logd["analytic"], v)
                    item["analytic"] = round(logd["analytic"][v], digits=3)
                    diff = abs.(m - logd["analytic"][v])
                    if sum(diff) > 0.2
                        item["anal_diff"] = round(diff, digits=3)
                    end
                end
                if haskey(logd, "stan") && haskey(logd["stan"], v)
                    item["stan"] = round.(logd["stan"][v], digits=3)
                    diff = abs.(m - logd["stan"][v])
                    if sum(diff) > 0.2
                        item["stan_diff"] = round.(diff, digits=3)
                    end
                end
                push!(data["turing_items"], item)
            end
        end
    end
    render(tmpl_log_string, data)
end

function stringify_log(logd::DataFrame, monitor=[])
    return string(logd)
end

print_log(log_data::Dict, monitor=[]) = print(stringify_log(log_data, monitor))

end #module


text:
# copied from https://github.com/JuliaLang/julia/pull/30924
# modified to be in compliance with C99: http://port70.net/~nsz/c/c99/n1256.html#6.7.2.2
# Cloned from github on 19/09/2021
# Modified to support StrEnum from Python
module PyEnum

import Core.Intrinsics.bitcast
export @pyenum

function namemap end
function name_value_pairs end

# Support both Integer and Strings
allowedtypes = Union{Integer,String}

struct PyenumValues{T}
    typename::Any
    value::T
end

# basetype(::Type{<:PyenumValues{T}}) where {T<:allowedtypes} = T

# (::Type{T})(x::PyenumValues{T2}) where {T<:allowedtypes,T2<:allowedtypes} = T(bitcast(T2, x.value))::T

Base.cconvert(::Type{T}, x::PyenumValues{T2}) where {T<:allowedtypes,T2<:allowedtypes} = T(x.value)
Base.write(io::IO, x::PyenumValues{T}) where {T<:allowedtypes} = write(io, T(x.value))

# Check this, could lead to potential problems:
Base.read(io::IO, ::Type{T}) where {T<:PyenumValues} = T(read(io, T(x.value)))

Base.isless(x::T, y::T) where {T<:PyenumValues} = isless(x.value, y.value)

Base.Symbol(x::PyenumValues{T}) where {T <: Integer} = get(namemap(x.typename), Integer(x.value), :UnknownMember)
Base.Symbol(x::PyenumValues{T}) where {T <: String} = get(namemap(x.typename), string(x.value), :UnknownMember)

Base.print(io::IO, x::PyenumValues) = print(io, Symbol(x.typename))

function Base.show(io::IO, x::PyenumValues)
    sym = x.value
    if !get(io, :compact, false)
        from = get(io, :module, Main)
        def = typeof(x.typename).name.module
        if from === nothing || !Base.isvisible(sym, def, from)
            show(io, def)
            print(io, ".")
        end
    end
    print(io, sym)
end

function Base.show(io::IO, ::MIME"text/plain", x::PyenumValues)
    print(io, x, "::")
    show(IOContext(io, :compact => true), typeof(x.value))
    print(io, " = ")
    show_pyenum_value(io, x.value)
end

function Base.show(io::IO, ::MIME"text/plain", t::Type{<:PyenumValues})
    print(io, "PyEnum ")
    print(io, typeof(x.value))
    print(io, ":")
    for (s, i) in name_value_pairs(t)
        print(io, "\n", Symbol(s), " = ")
        show_pyenum_value(io, i)
    end
end

function show_pyenum_value(io::IO, x::allowedtypes)
    if x isa Integer
        show(io, Integer(x))
    elseif x isa String
        show(io, String(x))
    else
        show(io, x)
    end
end

# give PyEnum types scalar behavior in broadcasting
Base.broadcastable(x::PyenumValues) = Ref(x)

@noinline enum_argument_error(typename, x) = throw(ArgumentError(string("input value out of range for PyEnum $(typename): $x")))

macro pyenum(T::Union{Symbol,Expr}, syms...)
    if isempty(syms)
        throw(ArgumentError("no arguments given for PyEnum $T"))
    end
    basetype = Int32
    typename = T
    if isa(T, Expr) && T.head == :(::) && length(T.args) == 2 && isa(T.args[1], Symbol)
        typename = T.args[1]
        basetype = Core.eval(__module__, T.args[2])
        if !isa(basetype, DataType)
            throw(ArgumentError("invalid base type for PyEnum $typename, $T=::$basetype"))
        end
    elseif !isa(T, Symbol)
        throw(ArgumentError("invalid type expression for PyEnum $T"))
    end
    seen = Set{Symbol}()
    name_values = Tuple{Symbol,basetype}[]
    namemap = Dict{basetype,Symbol}()
    # lo = hi = 0
    i = Nothing
    if basetype <: Integer
        i = zero(basetype)
    elseif basetype <: String
        i = typemin(basetype)
    else
        throw(ArgumentError("Invalid Datatype for PyEnum $basetype"))
    end

    if length(syms) == 1 && syms[1] isa Expr && syms[1].head == :block
        syms = syms[1].args
    end
    for s in syms
        s isa LineNumberNode && continue
        if isa(s, Symbol)
            # TODO: Check if it works with Strings
            if i == typemin(basetype) && !isempty(name_values)
                throw(ArgumentError("overflow in value \"$s\" of PyEnum $typename"))
            end
        elseif isa(s, Expr) &&
               (s.head == :(=) || s.head == :kw) &&
               length(s.args) == 2 && isa(s.args[1], Symbol)
            i = Core.eval(__module__, s.args[2]) # allow exprs, e.g. uint128, "1"
            i = convert(basetype, i)
            s = s.args[1]
        else
            throw(ArgumentError(string("invalid argument for PyEnum ", typename, ": ", s)))
        end
        if !Base.isidentifier(s)
            throw(ArgumentError("invalid name for PyEnum $typename; \"$s\" is not a valid identifier"))
        end
        haskey(namemap, i) || (namemap[i] = s;)
        if s in seen
            throw(ArgumentError("name \"$s\" in PyEnum $typename is not unique"))
        end
        push!(seen, s)
        push!(name_values, (s,i))
        # if length(name_values) == 1
        #     lo = hi = i
        # else
        #     lo = min(lo, i)
        #     hi = max(hi, i)
        # end
        # i += oneunit(i)

    end
    # print(typename)
    # print(basetype)
    # print(doc)
    # doc = basetype <: Integer ? :(primitive type $(typename) <: PyenumValues{$(basetype)} $(sizeof(basetype) * 8) end) :
    #         :(primitive type $(typename) <: Pyenum{$(basetype)} $(sizeof(Char) * 8) end)
    # println(doc)

    blk = quote
        # enum definition
        # Base.@__doc__($(doc))

        # if $(basetype) <: Integer
        #     function $(esc(typename))(x::Integer)
        #         x ≤ typemax(x) || enum_argument_error($(Expr(:quote, typename)), x)
        #         return bitcast($(esc(typename)), convert($(basetype), x))
        #     end
        # elseif $(basetype) <: String
        #     function $(esc(typename))(x::String)
        #         # Does not work, bitcast cant do anything for String type
        #         # return bitcast($(esc(typename)), convert($(basetype), x))
        # end
        if $(basetype) <: $(esc(allowedtypes))
            function $(esc(typename))(x::String)
                return PyenumValues($(esc(typename)), convert($(basetype), x))
            end
        else
            throw(ArgumentError("Invalid Datatype for PyEnum $basetype"))
        end
        PyEnum.namemap(::Type{$(esc(typename))}) = $(esc(namemap))
        PyEnum.name_value_pairs(::Type{$(esc(typename))}) = $(esc(name_values))
        # Base.typemin(x::Type{$(esc(typename))}) = $(esc(typename))($lo)
        # Base.typemax(x::Type{$(esc(typename))}) = $(esc(typename))($hi)
        let insts = (Any[$(esc(typename))(v[2]) for v in $name_values ]...,)
            Base.instances(::Type{$(esc(typename))}) = insts
        end
    end
    if isa(typename, Symbol)
        for (sym, i) in name_values
            push!(blk.args, :(const $(esc(sym)) = $(esc(typename))($i)))
        end
    end
    push!(blk.args, :nothing)
    blk.head = :toplevel
    return blk
end

include("operators_test.jl")

end # module

max_stars_repo_name: MiguelMarcelino/PyEnum

max_issues_repo_name: MiguelMarcelino/PyEnum

max_forks_repo_name: MiguelMarcelino/PyEnum


text:
"""
    @make_codeinfo begin
        <stmt>::<type>
    end

Create a typed `CodeInfo` object, if field `<type>` is not specified, it will
use `Any`.
"""
macro make_codeinfo(ex)
    esc(make_codeinfo_m(ex))
end

"""
    @make_ircode begin
        <stmt>::<type>
    end

Create a typed `IRCode` object, if field `<type>` is not specified, it will
use `Any`. See also [`make_codeinfo`](@ref).
"""
macro make_ircode(ex)
    quote
        let
            ci = $(make_codeinfo_m(ex))
            Core.Compiler.inflate_ir(ci)
        end
    end |> esc
end

function make_codeinfo_m(ex)
    Meta.isexpr(ex, :block) || error("expect a begin ... end")

    code = Expr(:ref, :Any)
    ssavaluetypes = Expr(:ref, :Any)

    for each in ex.args
        @switch each begin
            @case :($stmt::$type)
                push!(code.args, stmt)
                push!(ssavaluetypes.args, type)
            @case ::LineNumberNode
                continue
            @case _
                push!(code.args, each)
                push!(ssavaluetypes.args, :Any)
        end
    end

    @gensym ci nstmts
    quote
        let ci = (Meta.@lower 1 + 1).args[1]
            ci.code = $code
            nstmts = length(ci.code)
            ci.ssavaluetypes = $ssavaluetypes
            ci.codelocs = fill(Int32(1), nstmts)
            ci.ssaflags = fill(Int32(0), nstmts)
            ci
        end
    end
end

"""
    method_instance(f, tt, world=Base.get_world_counter())

Return the `MethodInstance`, unlike `Base.method_instances`, `tt` must be specified type.
"""
function method_instance(@nospecialize(f), @nospecialize(tt), world=Base.get_world_counter())
    # get the method instance
    meth = which(f, tt)
    sig = Base.signature_type(f, tt)::Type
    (ti, env) = ccall(:jl_type_intersection_with_env, Any, (Any, Any), sig,
        meth.sig)::Core.SimpleVector
    meth = Base.func_for_method_checked(meth, ti, env)
    return ccall(:jl_specializations_get_linfo, Ref{Core.MethodInstance},
        (Any, Any, Any, UInt), meth, ti, env, world)
end

macro test_codeinfo(ci, ex)
    esc(test_codeinfo_m(ci, ex))
end

function rm_code_coverage_effect(ci::CodeInfo)
    # make sure we don't change original one
    new = NewCodeInfo(copy(ci))
    for (v, stmt) in new
        if Meta.isexpr(stmt, :code_coverage_effect)
            delete!(new, v)
        elseif stmt isa GotoIfNot
            # idx always exists since code_coverage_effect was inserted
            idx = findfirst(stmt.dest:length(ci.code)) do k
                !Meta.isexpr(ci.code[k], :code_coverage_effect)
            end
            new[v] = GotoIfNot(stmt.cond, stmt.dest+idx-1)
        elseif stmt isa GotoNode
            # idx always exists since code_coverage_effect was inserted
            idx = findfirst(stmt.label:length(ci.code)) do k
                !Meta.isexpr(ci.code[k], :code_coverage_effect)
            end
            new[v] = GotoNode(stmt.label+idx-1)
        end
    end
    new_ci = finish(new)

    if ci.ssavaluetypes isa Vector
        types = []
        for (v, stmt) in enumerate(ci.code)
            if !Meta.isexpr(stmt, :code_coverage_effect)
                push!(types, ci.ssavaluetypes[v])
            end
        end

        new_ci.ssavaluetypes = types
        new_ci.inferred = ci.inferred
    end
    return new_ci
end

function test_codeinfo_m(ci, ex::Expr)
    Meta.isexpr(ex, :block) || error("expect a begin ... end")
    @gensym nocodedov_ci

    ret = quote
        $nocodedov_ci = $CompilerPluginTools.rm_code_coverage_effect($ci)
    end
    stmt_count = 1
    for each in ex.args
        @switch each begin
            @case :($stmt::$type)
                push!(ret.args, quote
                    @test $MLStyle.@match $nocodedov_ci.code[$stmt_count] begin
                        $stmt => true
                        _ => false
                    end
                end)
                push!(ret.args, quote
                    @test $nocodedov_ci.ssavaluetypes[$stmt_count] == $type
                end)
                stmt_count += 1
            @case ::LineNumberNode
                continue
            @case _
                push!(ret.args, quote
                    @test $MLStyle.@match $nocodedov_ci.code[$stmt_count] begin
                        $each => true
                        _ => false
                    end
                end)
                stmt_count += 1
        end
    end
    return ret
end


text:
using DifferentialEquations
using DataInterpolations

function munge(A::Vector{Vector{Float64}})
    l1 = length(A)
    l2 = length(A[1])
    M = zeros(l2,l1)
    for  i = 1:l1
        M[:,i] = A[i]
    end
    M
end

function lorenz!(du,u,p,t)
    du[1] = p[1]*(u[2]-u[1])
    du[2] = u[1]*(p[2]-u[3]) - u[2]
    du[3] = u[1]*u[2] - p[3]*u[3]
end

u₀ = [1,0,0]
p = [10,12,8/3] # Change d_lorenz
tspan = (0,10)
t = 0:0.01:10
prob = ODEProblem(lorenz!,u₀,tspan,p)
sol = solve(prob,Tsit5(),saveat=t)

U = munge(sol.u)
T = sol.t

interpq = QuadraticInterpolation(U,T)
DataInterpolations.derivative(interpq,1)

DataInterpolations.derivative(interpq,10)

interpl = LagrangeInterpolation(U,T)
DataInterpolations.derivative(interpl,1)

DataInterpolations.derivative(interpl,10)


text:
module Avalon

export
    grad,
    # initialization
    init_constant!,
    init_kaiming_normal!,
    # conv
    conv1d,
    conv2d,
    conv3d,
    maxpool2d,
    # activations
    logistic,
    sigmoid,
    softplus,
    softsign,
    logsigmoid,
    relu,
    leakyrelu,
    elu,
    softmax,
    logsoftmax,
    # losses
    nllloss,
    crossentropyloss,
    mseloss,
    NLLLoss,
    CrossEntropyLoss,
    MSELoss,
    # general layers
    Linear,
    Sequential,
    # conv layers
    Conv1d,
    Conv2d,
    Conv3d,
    # batchnorm
    BatchNorm2d,
    # RNN layers
    RNN,
    LSTM,
    GRU,
    init_hidden,
    # optim
    update!,
    SGD,
    RMSprop,
    Adam,
    # training
    fit!,
    trainmode!,
    testmode!,
    # device API (reexport from Yota)
    best_available_device
    # metrics (not exported by default)

include("core.jl")

end # module


text:
using Random
using MLJ
using DataFrames

"""
Split a given Dataframe into 2 sub-dataframe allowing
for training/testing dataset split
"""
function train_test_split(df::DF; split=0.8, seed=42, shuffled=true) where {DF <: Union{DataFrame, SubDataFrame}}
  Random.seed!(seed)
  (nr, nc) = size(df)
  nrp = round(Int, nr * split)
  row_ixes = shuffled ? shuffle(1:nr) : collect(1:nr)

  df_train = view(df[row_ixes, :], 1:nrp, 1:nc)
  df_test = view(df[row_ixes, :], nrp+1:nr, 1:nc)

  (df_train, df_test)
end

"""
From a Dataframe to matrix
from week2
"""
function get_data(df::DF, features, output) where {DF <: Union{DataFrame, SubDataFrame}}
  df[:, :constant] .= 1.0 # df.constant = fill(1.0, size(df, 1))
  features = [:constant, features...]
  X_matrix = convert(Matrix, select(df, features)) # to get a matrix
  y = df[!, output]                                # => to get a vector
  return (X_matrix, y)
end

"""
Assume feature_matrix is a matrix containing the features as columns
and weights is a corresponding array

from week2
"""
function predict_output(X::Matrix{T}, weights::Vector{T}) where {T <: Real}
    X * weights
end

"""
Calculate RSS given weights vector

depends on predict_output
from week2

y::SubArray{Float64,1, Array{Float64,1}, Tuple{UnitRange{Int64}}, true}
"""
function calc_rss(X::Matrix{T}, y, weights::Vector{T})  where {T <: Real}
  preds = predict_output(X, weights)
  sum((preds .- y) .^ 2)
end

"""
Calculate RSS given a (MLJ) machine which is a wrapper combining a model
and its data.

depends on MLJ predict()

mach::MLJBase.Machine{MLJLinearModels.RidgeRegressor,true}
"""
function get_rss(mach::Machine, X::DF, y::Vector{T}) where {DF <: Union{DataFrame, SubDataFrame}, T <: Real}
  ŷ = predict(mach, X)     # First get the predictions
  diff = y .- ŷ            # Then compute the residuals/errors
  rss = sum(diff .* diff)  # Then square and add them up
  return rss
end


text:
# ENV["PYTHON"] = /usr/bin/python3
# Pkg.build("PyCall")

using PyCall
@pyimport hpipm_python as hp
@pyimport numpy as np

qp_data = hp.hpipm_data()

A = np.array([1, 0, 1, 1])
B = np.array([0, 1])
b = np.array([0, 0])

Q = np.array([1, 0, 0, 1])
S = np.array([0, 0])
R = np.array([1])
q = np.array([1, 1])
r = np.array([0])

qp_data[:A] = PyVector([A, A, A, A, A])
qp_data[:B] = PyVector([B, B, B, B, B])
qp_data[:b] = PyVector([b, b, b, b, b])
qp_data[:Q] = PyVector([Q, Q, Q, Q, Q, Q])
qp_data[:S] = PyVector([S, S, S, S, S, S])
qp_data[:R] = PyVector([R, R, R, R, R, R])
qp_data[:q] = PyVector([q, q, q, q, q, q])
qp_data[:r] = PyVector([r, r, r, r, r, r])

x0 = np.array([1, 1])

qp_data[:d_lb] = PyVector([x0])
qp_data[:d_ub] = PyVector([x0])

qp_data[:idxb] = PyVector([np.array([1, 2])])

qp_dims = hp.hpipm_dims()

qp_dims[:nx]   = np.array([2, 2, 2, 2, 2, 2])
qp_dims[:nu]   = np.array([1, 1, 1, 1, 1, 0])
qp_dims[:nb]   = np.array([2, 0, 0, 0, 0, 0])
qp_dims[:nbx]  = np.array([2, 0, 0, 0, 0, 0])
qp_dims[:nbu]  = np.array([0, 0, 0, 0, 0, 0])
qp_dims[:ng]   = np.array([0, 0, 0, 0, 0, 0])
qp_dims[:ns]   = np.array([0, 0, 0, 0, 0, 0])
qp_dims[:nsbx] = np.array([0, 0, 0, 0, 0, 0])
qp_dims[:nsbu] = np.array([0, 0, 0, 0, 0, 0])
qp_dims[:nsg]  = np.array([0, 0, 0, 0, 0, 0])
qp_dims[:N]    = 5

# set up solver
solver = hp.hpipm_solver(qp_dims, qp_data)

# solve qp
return_flag = solver[:solve]()

println("HPIPM returned with flag $(return_flag)\n")

if return_flag == 0
    println("-> QP solved Solution:\n")
    solver[:print_sol]()
else
    println("-> Solver failed.\n")
end


text:
@testset "Mask" begin
    @testset "unwrap" begin
        img = testpattern()
        m = Augmentor.Mask(img)
        @test Augmentor.unwrap(m) === img
    end
end

@testset "shouldapply" begin
    @testset "type arguments" begin
        ops = [Rotate90(), ColorJitter()]
        whats = [Augmentor.Mask(testpattern())]
        for op in ops, what in whats
            @test Augmentor.shouldapply(op, what) == Augmentor.shouldapply(typeof(op), typeof(what))
        end
    end

    @testset "affine operations" begin
        ops = [FlipX, FlipY, NoOp, Rotate, Rotate180, Rotate270, Rotate90,
               Scale, ShearX, ShearY]
        applicable = [Augmentor.Mask]
        notapplicable = []
        for op in ops
            for what in applicable
                @test Augmentor.shouldapply(op, what) == Val(true)
            end
            for what in notapplicable
                @test Augmentor.shouldapply(op, what) == Val(false)
            end
        end
    end

    @testset "color operations" begin
        ops = [ColorJitter, GaussianBlur]
        applicable = []
        notapplicable = [Augmentor.Mask]
        for op in ops
            for what in applicable
                @test Augmentor.shouldapply(op, what) == Val(true)
            end
            for what in notapplicable
                @test Augmentor.shouldapply(op, what) == Val(false)
            end
        end
    end
end


text:
using LinearAlgebra
using Random, Distributions
using Plots
Plots.scalefontsizes(1.5)

# Global Terms
Random.seed!(1234)

# Logistic Regression

function f(β)
    #=
    Objective Function
    x: n x d
    y: d x 1
    β: d x 1
    λ: penalty
    =#
    result =  sum((X * β - y) .^ 2) / 2 / n + λ *norm(β,1)
    return result
end

function ∇f(β)
    #=
    This is already noisy gradient.
    β: 100 x 1
    σ: noise
    =#
    return X' * (X * β - y) ./n .+ λ + rand(Normal(0,σ^2), (d,1))
end

function ACSA(ν)
    α(t) = 2/(t+1)
    γ(t) = 4 * 2 * L / (ν * t * (t+1))

    xt = xtag = zeros(d)
    f_history = zeros(η+1)
    f_history[1] = f(xt)

    for t in 1:η
        xp1 = (1 - α(t)) * (μ + γ(t)) / (γ(t) + (1-α(t)^2) * μ) * xtag
        xp2 = α(t)*((1-α(t)) * μ +γ(t))/(γ(t) + (1-α(t)^2) * μ) * xt
        xtmd = xp1 + xp2

        p1 = (∇f(xtmd) * (μ-1) * α(t) .- α(t) * λ .+ ((1-α(t))*μ + γ(t))*∇f(xt)) / (μ + γ(t))
        p2 = (p1 .- λ) .*n + X' * y
        xt = inv(X' * X) * p2

        xtag = α(t) * xt + (1- α(t)) * xtag
        f_history[t+1] = f(xt)
    end

    return f_history
end

function ACSA_modified(ν,t_empirical)

    α(p) = 2/(p+1)
    γ(p) = 4 * 2 * L / (ν * p * (p+1))
    xt = xtag = zeros(d)
    f_history = zeros(η+1)
    f_history[1] = f(xt)

    for t in 1:t_empirical
        xp1 = (1 - α(t)) * (μ + γ(t)) / (γ(t) + (1-α(t)^2) * μ) * xtag
        xp2 = α(t)*((1-α(t)) * μ +γ(t))/(γ(t) + (1-α(t)^2) * μ) * xt
        xtmd = xp1 + xp2

        p1 = (∇f(xtmd) * (μ-1) * α(t) .- α(t) * λ .+ ((1-α(t))*μ + γ(t))*∇f(xt)) / (μ + γ(t))
        p2 = (p1 .- λ) .*n + X' * y
        xt = inv(X' * X) * p2

        xtag = α(t) * xt + (1- α(t)) * xtag
        f_history[t+1] = f(xt)

        if t % 100 == 0
            @show t
        end
    end
    γ_1(p) = 4 * L / (ν * p * (p+1)) # + (2*10)/(ν+sqrt(p))
    for t in t_empirical+1:η
        xp1 = (1 - α(t)) * (μ + γ_1(t)) / (γ_1(t) + (1-α(t)^2) * μ) * xtag
        xp2 = α(t)*((1-α(t)) * μ +γ_1(t))/(γ_1(t) + (1-α(t)^2) * μ) * xt
        xtmd = xp1 + xp2

        p1 = (∇f(xtmd) * (μ-1) * α(t) .- α(t) * λ .+ ((1-α(t))*μ + γ_1(t))*∇f(xt)) / (μ + γ_1(t))
        p2 = (p1 .- λ) .*n + X' * y
        xt = inv(X' * X) * p2

        xtag = α(t) * xt + (1- α(t)) * xtag
        f_history[t+1] = f(xt)

        if t % 100 == 0
            @show t
        end
    end

    return f_history
end

## Data setup
global n = 2000   # number of datapoints
global d = 100    # number of dimensions
global λ = 0.01   # Penalty for beta term
global ρ = 0
global X = rand(Normal(0,1), (n,d))
β = rand(Normal(0, 1), (d, 1))
global βstart = β
global y = X * β  .+ rand(Normal(0,1),(n,1))

# Experiments
global K = 10     # number of iterations for MASG
global η = 2000   # Number of iterations for other methods
global L = maximum(eigvals(X' * X))
# global L = norm(X,2)^2/n
global μ = 0

global σ = 1e-1

ν = norm(X,2) / n
f_ACSA = ACSA(ν)
f_ACSA_modified_300 = ACSA_modified(ν,300)
# f_ACSA_modified_500 = ACSA_modified(ν,500)
#f_ACSA_modified_5000 = ACSA_modified(ν,5000)
f_ACSA_modified_0 = ACSA_modified(ν,0)
# f_ACSA_modified_500 = ACSA_modified(ν,50)

plot(f_ACSA, label="ACSA",w=2, color=:red,xscale=:log,yscale=:log,legend=:bottomleft)
plot!(f_ACSA_modified_300, label="ACSA+300",w=2, color=:blue)
plot!(f_ACSA_modified_0, label="ACSA+0",w=2, color=:orange)

savefig("./figures/ACSA_new_1.png")

plot!(f_ACSA_modified_500, label="ACSA+500",w=2, color=:green)
plot!(f_ACSA_modified_5000, label="ACSA+1000",w=2, color=:purple)
# plot!(f_ACSA_modified_500, label="ACSA+500",w=2, color=:orange)

γ(p) = 4 * 2 * L / (ν * p * (p+1))
γ_1(p) = 4 * L / (ν * p * (p+1))#  + (2*10)/(ν+sqrt(p))

gridline = [1:2000;]
plot(γ.(gridline),xscale=:log,yscale=:log)
plot!(γ_1.(gridline))

max_stars_repo_name: chenyangzhu/noisy-gradients

max_issues_repo_name: chenyangzhu/noisy-gradients

max_forks_repo_name: chenyangzhu/noisy-gradients


text:
@testset "table" begin
    # drop database "test_db" if exists
    drop_db(DBNAME)

    # create db
    create_database(DBNAME)
    @test DBNAME in from_sql("SELECT datname FROM pg_database;")[!, 1]

    # enable uuid utils in "test_db"
    @test enable_uuid(DBNAME) == 0

    # modify dbconfig
    dbconfig = copy(current_dbconfig())
    dbconfig[:dbname] = DBNAME

    # create "SqueezedStatesData" table including uuid column into "test_db"
    create_table(string(SqueezedStatesData), gen_table_schema(SqueezedStatesData), dbconfig=dbconfig)

    column_names = [
        "id",
        "r", "theta",
        "dim", "rho",
        "n_points", "bhd",
        "w_range", "w",
    ]

    column_names_from_sql = from_sql(
        """
        SELECT column_name
            FROM information_schema.columns
            WHERE table_schema='public' AND table_name='$(string(SqueezedStatesData))';
        """,
        dbconfig=dbconfig
    )[!, 1]

    @test all(lowercase(col) in column_names_from_sql for col in column_names)

    # create "SqueezedThermalStatesData" table including uuid column into "test_db"
    create_table(string(SqueezedThermalStatesData), gen_table_schema(SqueezedThermalStatesData), dbconfig=dbconfig)

    column_names = [
        "id",
        "r", "theta", "nbar",
        "dim", "rho",
        "n_points", "bhd",
        "w_range", "w",
    ]

    column_names_from_sql = from_sql(
        """
        SELECT column_name
            FROM information_schema.columns
            WHERE table_schema='public' AND table_name='$(string(SqueezedThermalStatesData))';
        """,
        dbconfig=dbconfig
    )[!, 1]

    @test all(lowercase(col) in column_names_from_sql for col in column_names)
end


text:
# Read status of the power supply

# command not working reliably ???

export psp_status!, psp_status

const READ_STATUS = [0x4c, 0x0d]

"""
Power supply status

**fields**

- `voltage`     -- Output voltage
- `current`     -- Output current
- `power`       -- Output power
- `vlimit`      -- Voltage limit
- `ilimit`      -- Current limit
- `plimit`      -- Power limit
- `isrelayon`   -- True if output relay is on
- `isoverheat`  -- True if power supply temperature is too high
- `isknobfine`  -- True if knob is in fine mode
- `isknobunlock`-- True if knob is unlocked
- `isremote`    -- True if power supply is in remote control mode
- `islock`      -- True if front panel controls are locked

"""
type Status
  "Output voltage"
  voltage :: Float64
  "Output current"
  current :: Float64
  "Output power"
  power   :: Float64
  "Voltage limit"
  vlimit  :: Float64
  "Current limit"
  ilimit  :: Float64
  "Power limit"
  plimit  :: Float64
  "True if output relay is on"
  isrelayon  :: Bool
  "True if power supply temperature is too high"
  isoverheat :: Bool
  "True if knob is in fine mode"
  isknobfine :: Bool
  "True if knob is unlocked"
  isknobunlock :: Bool
  "True if power supply is in remote control mode"
  isremote   :: Bool
  "True if front panel controls are locked"
  islock     :: Bool
end

Status() = Status(0.0,0.0,0.0,0.0,0.0,0.0,false,false,false,false,false,false)

function Base.show(io::IO, s::Status)
  println(io,"Status")
  println(io,"voltage = ",s.voltage)
  println(io,"curent  = ",s.current)
  println(io,"power   = ",s.power)
  println(io,"vlimit  = ",s.vlimit)
  println(io,"ilimit  = ",s.ilimit)
  println(io,"plimit  = ",s.plimit)
  println(io,"isrelayon    = ",s.isrelayon)
  println(io,"isoverheat   = ",s.isoverheat)
  println(io,"isknobfine   = ",s.isknobfine)
  println(io,"isknobunlock = ",s.isknobunlock)
  println(io,"isremote     = ",s.isremote)
  println(io,"islock       = ",s.islock)
end

"""
    psp_status!(io, status)

Read status of the power supply and update `status`.
"""
function psp_status!(io_psp::IO, s::Status)
    write(io_psp, READ_STATUS)  # obtain all status from power supply
    status = bytestring(read(io_psp,UInt8, 40))
    if  status[1] != 'V' ||
        status[4] != '.' ||
        status[7] != 'A' ||
        status[9] != '.' ||
        status[13] != 'W' ||
        status[17] != '.' ||
        status[19] != 'U' ||
        status[22] != 'I' ||
        status[24] != '.' ||
        status[27] != 'P' ||
        status[31] != 'F' ||
        status[38:40] != "\r\r\n"
          error("PSP Invalid Response: ",status)
    end
    s.voltage = parse(Float64,status[2:6])
    s.current = parse(Float64,status[8:12])
    s.power = parse(Float64,status[14:18])
    s.vlimit = parse(Float64,status[20:21])
    s.ilimit = parse(Float64,status[23:26])
    s.plimit = parse(Float64,status[28:30])
    s.isrelayon = status[32] == '1'
    s.isoverheat = status[33] == '1'
    s.isknobfine = status[34] == '1'
    s.isknobunlock = status[35] == '1'
    s.isremote = status[36] == '1'
    s.islock = status[37] == '1'
    return s
end

"""
    psp_status(io)

Read status of the power supply and return status.
"""
function psp_status(io_psp::IO)
  s = Status()
  psp_status!(io_psp,s)
end


text:
##################################################
## CasadiSymbolicObject types have field x::PyCall.PyObject

## Symbol class for controlling dispatch
abstract type CasadiSymbolicObject <: Real end

struct SX <: CasadiSymbolicObject
    __pyobject__::PyCall.PyObject
end
struct MX <: CasadiSymbolicObject
    __pyobject__::PyCall.PyObject
end

Base.convert(::Type{C}, x::PyCall.PyObject) where C <: CasadiSymbolicObject = C(x)

##################################################

## important override
## this allows most things to flow though PyCall
PyCall.PyObject(x::CasadiSymbolicObject) = x.__pyobject__

##################################################

## text/plain
jprint(x::CasadiSymbolicObject) = PyCall.pycall(pybuiltin("str"), String, PyObject(x))
Base.show(io::IO, s::CasadiSymbolicObject) = print(io, jprint(s))

function Base.getproperty(o::C, s::Symbol) where C <: CasadiSymbolicObject
    if s in fieldnames(C)
        getfield(o, s)
    else
        getproperty(PyCall.PyObject(o), s)
    end
end


text:
using DifferentialEquations, DiffEqBiological, Latexify, Plots
fmt = :svg
pyplot(fmt=fmt)
rn = @reaction_network begin
    hillr(D₂,α,K,n), ∅ --> m₁
    hillr(D₁,α,K,n), ∅ --> m₂
    (δ,γ), m₁ ↔ ∅
    (δ,γ), m₂ ↔ ∅
    β, m₁ --> m₁ + P₁
    β, m₂ --> m₂ + P₂
    μ, P₁ --> ∅
    μ, P₂ --> ∅
    (k₊,k₋), 2P₁ ↔ D₁
    (k₊,k₋), 2P₂ ↔ D₂
    (k₊,k₋), P₁+P₂ ↔ T
end α K n δ γ β μ k₊ k₋;

latexify(rn; env=:chemical)

x = latexify(rn; env=:chemical, starred=true, mathjax=true);
display("text/latex", "$x");

species(rn)

params(rn)

substratesymstoich(rn, 11)

substratesymstoich.(rn, 1:numreactions(rn))

netstoich.(rn, 1:numreactions(rn))

rxtospecies_depgraph(rn)

species(rn)[[3,4,7]]

speciestorx_depgraph(rn)[1]

findall(depset -> in(:m₁, depset), dependents.(rn, 1:numreactions(rn)))

rxtorx_depgraph(rn)

rnmin = @min_reaction_network begin
    (δ,γ), m₁ ↔ ∅
    (δ,γ), m₂ ↔ ∅
    β, m₁ --> m₁ + P₁
    β, m₂ --> m₂ + P₂
    μ, P₁ --> ∅
    μ, P₂ --> ∅
end δ γ β μ;

addspecies!(rnmin, :D₁)
addspecies!(rnmin, :D₂)
addspecies!(rnmin, :T)

addparam!(rnmin, :α)
addparam!(rnmin, :K)
addparam!(rnmin, :n)
addparam!(rnmin, :k₊)
addparam!(rnmin, :k₋)

addreaction!(rnmin, :(hillr(D₁,α,K,n)), :(∅ --> m₂))
addreaction!(rnmin, :((k₊,k₋)), :(2P₂ ↔ D₂))
addreaction!(rnmin, :k₊, :(2P₁ --> D₁))
addreaction!(rnmin, :k₋, :(D₁ --> 2P₁))

# signature is addreaction!(rnmin, paramexpr, substratestoich, productstoich)
addreaction!(rnmin, :(hillr(D₂,α,K,n)), (), (:m₁ => 1,))
addreaction!(rnmin, :k₊, (:P₁=>1, :P₂=>1), (:T=>1,))
addreaction!(rnmin, :k₋, (:T=>1,), (:P₁=>1, :P₂=>1))

setdiff(species(rn), species(rnmin))

setdiff(params(rn), params(rnmin))

rxidx = numreactions(rn)
setdiff(substrates(rn, rxidx), substrates(rnmin, rxidx))

setdiff(products(rn, rxidx), products(rnmin, rxidx))

rateexpr(rn, rxidx) == rateexpr(rnmin, rxidx)

addodes!(rnmin)

odeexprs(rnmin)

latexify(rnmin)

x = latexify(rnmin, starred=true);
display("text/latex", "$x");

latexify(jacobianexprs(rnmin))

x = latexify(jacobianexprs(rnmin), starred=true);
display("text/latex", "$x");

N = 64
h = 1 / N

rn = @empty_reaction_network

for i = 1:N
    addspecies!(rn, Symbol(:u, i))
end

addparam!(rn, :β)

for i = 1:N
    (i < N) && addreaction!(rn, :β, (Symbol(:u,i)=>1,), (Symbol(:u,i+1)=>1,))
    (i > 1) && addreaction!(rn, :β, (Symbol(:u,i)=>1,), (Symbol(:u,i-1)=>1,))
end

addodes!(rn)

u₀ = zeros(N)
u₀[div(N,2)] = 10000
p = [1/(h*h)]
tspan = (0.,.01)
oprob = ODEProblem(rn, u₀, tspan, p)

sol = solve(oprob, Rodas5())
times = [0., .0001, .001, .01]
plt = plot()
for time in times
    plot!(plt, 1:N, sol(time), fmt=fmt, xlabel="i", ylabel="uᵢ", label=string("t = ", time), lw=3)
end
plot(plt, ylims=(0.,10000.))

addjumps!(rn, build_regular_jumps=false, minimal_jumps=true)

# make the initial condition integer valued
u₀ = zeros(Int, N)
u₀[div(N,2)] = 10000

# setup and solve the problem
dprob = DiscreteProblem(rn, u₀, tspan, p)
jprob = JumpProblem(dprob, DirectCR(), rn, save_positions=(false,false))
jsol = solve(jprob, SSAStepper(), saveat=times)

times = [0., .0001, .001, .01]
plts = []
for i = 1:4
    b = bar(1:N, jsol[i], legend=false, fmt=fmt, xlabel="i", ylabel="uᵢ", title=string("t = ", times[i]))
    plot!(b,sol(times[i]))
    push!(plts,b)
end
plot(plts...)

using DiffEqTutorials
DiffEqTutorials.tutorial_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file], remove_homedir=true)


text:
using JLD2;
using Distributed;
using Printf;
using IterTools;
using Distributions
@everywhere include("runit.jl");
@everywhere include("../utilities/thresholds.jl");
include("../utilities/experiment_helpers.jl");

dist = Gaussian();
L = 0.01; # Lipschitz constant
arms = [0,96,97,98, 99, 100, 101, 102,103,104]
μ = [ 1.06,0.99, 0.99, 0.99, 0.99, 1, 0.99,0.99,0.99,0.99];

pep = LipschitzBestArm(dist, arms, μ, L);
# methods to be compared
srs = [FWSampling(), Menard(CTracking, 1/oracle(pep, μ)[1]), TrackAndStop(DTracking)];

δs =  (0.1,0.01); # confidence
βs = GK16.(δs);
repeats = 100;
seed = 1234;

println("arms=$arms, μ=$μ, repeats=$repeats");
# compute
@time data = pmap(
    ((sr, i),) -> runit(seed + i, sr, μ, pep, βs),
    Iterators.product(srs, 1:repeats),
);
dump_stats(pep, μ, δs, βs, srs, data, repeats);
# save
@save isempty(ARGS) ? "BAI2.dat" : ARGS[1] dist μ pep srs data δs βs repeats seed
# visualise by loading viz_bai.jl


text:
module DebuggingUtilities

export @showln, @showfl, test_showline, time_showline

"""
DebuggingUtilities contains a few tools that may help debug julia code. The
exported tools are:

- `@showln`: a macro for displaying variables and corresponding function, file, and line number information
- `@showfl`: a crude, but faster, version of `@showln`
- `test_showline`: a function that displays progress as it executes a file
- `time_showline`: a function that displays execution time for each expression in a file
"""
DebuggingUtilities

## @showln

# look up an instruction pointer

function print_btinfo(io, frm)
    print(io, "in ", frm.func, " at ", frm.file, ":", frm.line)
end
function show_backtrace1(io, bt)
    st = stacktrace(bt)
    for frm in st
        funcname = frm.func
        if funcname != :backtrace && funcname != Symbol("macro expansion")
            print_btinfo(io, frm)
            break
        end
    end
end

mutable struct FlushedIO <: IO
    io
end

function Base.println(io::FlushedIO, args...)
    println(io.io, args...)
    flush(io.io)
end
Base.getindex(io::FlushedIO) = io.io
Base.setindex!(io::FlushedIO, newio) = io.io = newio

const showlnio = FlushedIO(stdout)

"""
`@showln x` prints "x = val", where `val` is the value of `x`, along
with information about the function, file, and line number at which
this statement was executed. For example:

    function foo()
        x = 5
        @showln x
        x = 7
        @showln x
        nothing
    end

might produce output like

            x = 5
            (in foo at ./error.jl:26 at /tmp/showln_test.jl:52)
            x = 7
            (in foo at ./error.jl:26 at /tmp/showln_test.jl:54)

This macro causes a backtrace to be taken, and looking up the
corresponding code information is relatively expensive, so using
`@showln` can have a substantial performance cost.

The indentation of the line is proportional to the length of the
backtrace, and consequently is an indication of recursion depth.

Line numbers are not typically correct on julia-0.4.
"""
macro showln(exs...)
    blk = showexprs(exs)
    blk = quote
        local bt = backtrace()
        local indent = length(bt)  # to mark recursion
        $blk
        print(showlnio[], " "^indent*"(")
        show_backtrace1(showlnio[], bt)
        println(showlnio[], ")")
    end
    if !isempty(exs)
        push!(blk.args, :value)
    end
    blk
end

function showexprs(exs)
    blk = Expr(:block)
    for ex in exs
        push!(blk.args, :(println(showlnio[], " "^indent, sprint(Base.show_unquoted,$(Expr(:quote, ex)),indent)*" = ", repr(begin value=$(esc(ex)) end))))
    end
    blk
end

"""
`@showfl(@__FILE__, @__LINE__, expressions...)` is similar to
`@showln`, but has much less overhead (and is uglier to use).
"""
macro showfl(fl, ln, exs...)
    blk = showexprs(exs)
    blk = quote
        local indent = 0
        $blk
        println(showlnio[], "(at file ", $fl, ", line ", $ln, ')')
    end
    if !isempty(exs)
        push!(blk.args, :value)
    end
    blk
end

"""
`test_showline(filename)` is equivalent to `include(filename)`, except
that it also displays the expression and file-offset (in characters) for
each expression it executes. This can be useful for debugging errors,
especially those that cause a segfault.
"""
function test_showline(filename)
    str = read(filename, String)
    Core.eval(Main, Meta.parse("using Test"))
    idx = 1
    while idx < length(str)
        ex, idx = Meta.parse(str, idx)
        try
            println(showlnio, idx, ": ", ex)
        catch
            println(showlnio, "failed to print line starting at file-offset ", idx)
        end
        Core.eval(Main,ex)
    end
    println(showlnio, "done")
end

"""
`time_showline(filename)` is equivalent to `include(filename)`, except
that it also analyzes the time expended on each expression within the
file. Once finished, it displays the file-offset (in characters),
elapsed time, and expression in order of increasing duration.  This
can help you identify bottlenecks in execution.

This is less useful now that julia has package precompilation, but can
still be handy on occasion.
"""
function time_showline(filename)
    str = read(filename, String)
    idx = 1
    exprs = Any[]
    t = Float64[]
    pos = Int[]
    while idx < length(str)
        ex, idx = Meta.parse(str, idx)
        push!(exprs, ex)
        push!(t, @elapsed Core.eval(Main,ex))
        push!(pos, idx)
    end
    perm = sortperm(t)
    for p in perm
        print(showlnio[], pos[p], " ", t[p], ": ")
        str = string(exprs[p])
        println(showlnio[], split(str, '\n')[1])
    end
    println(showlnio[], "done")
end

function __init__()
    showlnio[] = stdout
end

end # module


text:
"""
Module JustSayIt

Enables offline, low latency, highly accurate speech to command translation and is usable as software or API.

# General overview and examples
https://github.com/omlins/JustSayIt.jl

# Software usage
```
> julia
julia> using JustSayIt
julia> just_say_it()
```
Type `?just_say_it` to learn about customization keywords.

# Application Programming Interface (API)

#### Macros
- [`@voiceargs`](@ref)

#### Functions
- [`just_say_it`](@ref)
- [`init_jsi`](@ref)
- [`finalize_jsi`](@ref)
- [`is_next`](@ref)
- [`are_next`](@ref)
- [`pause_recording`](@ref)
- [`restart_recording`](@ref)

#### Submodules
- [`Help`](@ref)
- [`Keyboard`](@ref)
- [`Mouse`](@ref)
- [`Email`](@ref)
- [`Internet`](@ref)

To see a description of a function, macro or module type `?<functionname>`, `?<macroname>` (including the `@`) or `?<modulename>`, respectively.
"""
module JustSayIt

## Alphabetical include of submodules, except commmand-submodules (below)
include("Exceptions.jl")
using .Exceptions

## Include of shared constant parameters, types and syntax sugar and voiceargs macro
include("shared.jl")
include("voiceargs.jl")

## Alphabetical include of files
include("finalize_jsi.jl")
include("next_token.jl")
include("init_jsi.jl")
include("recorder.jl")

## Include of command-submodules for peripherics control and help
include("Commands/Keyboard.jl")
include("Commands/Mouse.jl")
include("Commands/Help.jl")

## Alphabetical include of command-submodules (must be at end as needs to import from JustSayIt, .e.g. next_recognition, next_partial_recognition)
include("Commands/Email.jl")
include("Commands/Internet.jl")

## Include of main application (must be at end as needs to import potentially anything available in JustSayIt, in particular the Commands submodule).
include("just_say_it.jl")

## Exports (need to be after include of submodules if re-exports from them)
export just_say_it
export @voiceargs
export init_jsi, finalize_jsi, is_next, are_next, pause_recording, restart_recording
export Keyboard, Mouse, Help, Email, Internet

end # module JustSayIt


text:
using QSimulator

function speed_test(repeats)
	T_f = 1000.0
	sys = sim_setup(4, T_f)
	dt = 1.0

	# run once for JIT
	unitary_propagator(sys, dt, 0.0, T_f, parallelize=false)
	unitary_propagator(sys, dt, 0.0, T_f, parallelize=true)

	# then time each
	tserial = @elapsed for ct = 1:repeats
		unitary_propagator(sys, dt, 0.0, T_f, parallelize=false)
	end

	tparallel = @elapsed for ct = 1:repeats
		unitary_propagator(sys, dt, 0.0, T_f, parallelize=true)
	end

	println("Serial execution in $(tserial/repeats) seconds")
	println("Parallel execution in $(tparallel/repeats) seconds")

	@test isapprox(abs(tserial/repeats - tparallel/repeats), 0.0, atol = 1e-1)
end

function sim_setup(dimension, numTimeSteps)
	α = -0.350
	J = 0.002
	Ω = 0.010
	Q1 = Duffing("q1", 0.0, α, dimension)
	Q2 = Duffing("q2", 1.0, α, dimension)

	drive = MicrowaveControl("CR", 1.0, timeStep=1.0)
	sys = Q1 + Q2 + FlipFlop(Q1, Q2, J) + RotatingSemiClassicalDipole(Field(drive), Q2, Ω)

	seq = ones(round(Int, numTimeSteps))
	load_sequence!(drive, seq)

	return sys
end

@testset "speed test" begin
    speed_test(5)
end


text:
# %% Import NaMaster for the power spectrum analysis on a partial sky
# Reference: Alonso et al.,  MNRAS, 484, 4127 (2019), https://github.com/LSSTDESC/NaMaster
nmt = pyimport("pymaster")
# Read in and apodize mask (C^2 apodization, 10 degrees)
m = readMapFromFITS(maskfile, 1, Float64)
mask = nmt.mask_apodization(m.pixels, 10.0, apotype = "C2")
# Create binning scheme. We will use Δℓ multipoles per bandpower.
b = nmt.NmtBin.from_nside_linear(nside, Δℓ, is_Dell = true)
# Array with effective multipole per bandpower
ell_eff = b.get_effective_ells()
nbands = b.get_n_bands()
# The function defined below will compute the power spectrum between two
# NmtFields f_a and f_b, using the coupling matrix stored in the
# NmtWorkspace wsp and subtracting the deprojection bias clb.
# Note that the most expensive operations in the MASTER algorithm are
# the computation of the coupling matrix and the deprojection bias. Since
# these two objects are precomputed, this function should be pretty fast!
function compute_master(f_a, f_b, wsp)
    # Compute the power spectrum (a la anafast) of the masked fields
    # Note that we only use n_iter=0 here to speed up the computation,
    # but the default value of 3 is recommended in general.
    cl_coupled = nmt.compute_coupled_cell(f_a, f_b)
    # Decouple power spectrum into bandpowers inverting the coupling matrix
    cl_decoupled = wsp.decouple_cell(cl_coupled)
    return cl_decoupled
end


text:
export Freq

"""
    Freq(min, max, [ndec])
    Freq(<keyword arguments>)

`Freq` objects `show` a FastHenry .freq command.

Keyword arguments are `min`, `max`, `ndec`.
"""
struct Freq <: Element
  min :: Float64
  max :: Float64
  ndec :: Float64
  function Freq(min, max, ndec)
    if min>max
      throw(ArgumentError("max must be >= min"))
    end
    if min<0 || max<0 || ndec<0
      throw(ArgumentError("all parameters must be >0"))
    end
    if isnan(max)
      throw(ArgumentError("must specify max frequency"))
    end
    return new(min,max,ndec)
  end
end

Freq(min,max) = Freq(min,max,NaN)
Freq(;min=0, max=NaN, ndec=NaN) = Freq(min,max,ndec)


text:
# Julia wrapper for header: /home/astukalov/.julia/v0.5/Sundials/deps/usr/include/kinsol/kinsol.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function KINCreate()
    ccall((:KINCreate,libsundials_kinsol),KINMemPtr,())
end

function __KINSetErrHandlerFn(kinmem::KINMemPtr,ehfun::KINErrHandlerFn,eh_data::Ptr{Void})
    ccall((:KINSetErrHandlerFn,libsundials_kinsol),Cint,(KINMemPtr,KINErrHandlerFn,Ptr{Void}),kinmem,ehfun,eh_data)
end

function KINSetErrHandlerFn(kinmem,ehfun,eh_data)
    __KINSetErrHandlerFn(convert(KINMemPtr,kinmem),ehfun,pointer(eh_data))
end

function __KINSetErrFile(kinmem::KINMemPtr,errfp::Ptr{FILE})
    ccall((:KINSetErrFile,libsundials_kinsol),Cint,(KINMemPtr,Ptr{FILE}),kinmem,errfp)
end

function KINSetErrFile(kinmem,errfp)
    __KINSetErrFile(convert(KINMemPtr,kinmem),errfp)
end

function __KINSetInfoHandlerFn(kinmem::KINMemPtr,ihfun::KINInfoHandlerFn,ih_data::Ptr{Void})
    ccall((:KINSetInfoHandlerFn,libsundials_kinsol),Cint,(KINMemPtr,KINInfoHandlerFn,Ptr{Void}),kinmem,ihfun,ih_data)
end

function KINSetInfoHandlerFn(kinmem,ihfun,ih_data)
    __KINSetInfoHandlerFn(convert(KINMemPtr,kinmem),ihfun,pointer(ih_data))
end

function __KINSetInfoFile(kinmem::KINMemPtr,infofp::Ptr{FILE})
    ccall((:KINSetInfoFile,libsundials_kinsol),Cint,(KINMemPtr,Ptr{FILE}),kinmem,infofp)
end

function KINSetInfoFile(kinmem,infofp)
    __KINSetInfoFile(convert(KINMemPtr,kinmem),infofp)
end

function __KINSetUserData(kinmem::KINMemPtr,user_data::Any)
    ccall((:KINSetUserData,libsundials_kinsol),Cint,(KINMemPtr,Any),kinmem,user_data)
end

function KINSetUserData(kinmem,user_data)
    __KINSetUserData(convert(KINMemPtr,kinmem),user_data)
end

function __KINSetPrintLevel(kinmemm::KINMemPtr,printfl::Cint)
    ccall((:KINSetPrintLevel,libsundials_kinsol),Cint,(KINMemPtr,Cint),kinmemm,printfl)
end

function KINSetPrintLevel(kinmemm,printfl)
    __KINSetPrintLevel(convert(KINMemPtr,kinmemm),convert(Cint,printfl))
end

function __KINSetNumMaxIters(kinmem::KINMemPtr,mxiter::Clong)
    ccall((:KINSetNumMaxIters,libsundials_kinsol),Cint,(KINMemPtr,Clong),kinmem,mxiter)
end

function KINSetNumMaxIters(kinmem,mxiter)
    __KINSetNumMaxIters(convert(KINMemPtr,kinmem),convert(Clong,mxiter))
end

function __KINSetNoInitSetup(kinmem::KINMemPtr,noInitSetup::Cint)
    ccall((:KINSetNoInitSetup,libsundials_kinsol),Cint,(KINMemPtr,Cint),kinmem,noInitSetup)
end

function KINSetNoInitSetup(kinmem,noInitSetup)
    __KINSetNoInitSetup(convert(KINMemPtr,kinmem),convert(Cint,noInitSetup))
end

function __KINSetNoResMon(kinmem::KINMemPtr,noNNIResMon::Cint)
    ccall((:KINSetNoResMon,libsundials_kinsol),Cint,(KINMemPtr,Cint),kinmem,noNNIResMon)
end

function KINSetNoResMon(kinmem,noNNIResMon)
    __KINSetNoResMon(convert(KINMemPtr,kinmem),convert(Cint,noNNIResMon))
end

function __KINSetMaxSetupCalls(kinmem::KINMemPtr,msbset::Clong)
    ccall((:KINSetMaxSetupCalls,libsundials_kinsol),Cint,(KINMemPtr,Clong),kinmem,msbset)
end

function KINSetMaxSetupCalls(kinmem,msbset)
    __KINSetMaxSetupCalls(convert(KINMemPtr,kinmem),convert(Clong,msbset))
end

function __KINSetMaxSubSetupCalls(kinmem::KINMemPtr,msbsetsub::Clong)
    ccall((:KINSetMaxSubSetupCalls,libsundials_kinsol),Cint,(KINMemPtr,Clong),kinmem,msbsetsub)
end

function KINSetMaxSubSetupCalls(kinmem,msbsetsub)
    __KINSetMaxSubSetupCalls(convert(KINMemPtr,kinmem),convert(Clong,msbsetsub))
end

function __KINSetEtaForm(kinmem::KINMemPtr,etachoice::Cint)
    ccall((:KINSetEtaForm,libsundials_kinsol),Cint,(KINMemPtr,Cint),kinmem,etachoice)
end

function KINSetEtaForm(kinmem,etachoice)
    __KINSetEtaForm(convert(KINMemPtr,kinmem),convert(Cint,etachoice))
end

function __KINSetEtaConstValue(kinmem::KINMemPtr,eta::realtype)
    ccall((:KINSetEtaConstValue,libsundials_kinsol),Cint,(KINMemPtr,realtype),kinmem,eta)
end

function KINSetEtaConstValue(kinmem,eta)
    __KINSetEtaConstValue(convert(KINMemPtr,kinmem),eta)
end

function __KINSetEtaParams(kinmem::KINMemPtr,egamma::realtype,ealpha::realtype)
    ccall((:KINSetEtaParams,libsundials_kinsol),Cint,(KINMemPtr,realtype,realtype),kinmem,egamma,ealpha)
end

function KINSetEtaParams(kinmem,egamma,ealpha)
    __KINSetEtaParams(convert(KINMemPtr,kinmem),egamma,ealpha)
end

function __KINSetResMonParams(kinmem::KINMemPtr,omegamin::realtype,omegamax::realtype)
    ccall((:KINSetResMonParams,libsundials_kinsol),Cint,(KINMemPtr,realtype,realtype),kinmem,omegamin,omegamax)
end

function KINSetResMonParams(kinmem,omegamin,omegamax)
    __KINSetResMonParams(convert(KINMemPtr,kinmem),omegamin,omegamax)
end

function __KINSetResMonConstValue(kinmem::KINMemPtr,omegaconst::realtype)
    ccall((:KINSetResMonConstValue,libsundials_kinsol),Cint,(KINMemPtr,realtype),kinmem,omegaconst)
end

function KINSetResMonConstValue(kinmem,omegaconst)
    __KINSetResMonConstValue(convert(KINMemPtr,kinmem),omegaconst)
end

function __KINSetNoMinEps(kinmem::KINMemPtr,noMinEps::Cint)
    ccall((:KINSetNoMinEps,libsundials_kinsol),Cint,(KINMemPtr,Cint),kinmem,noMinEps)
end

function KINSetNoMinEps(kinmem,noMinEps)
    __KINSetNoMinEps(convert(KINMemPtr,kinmem),convert(Cint,noMinEps))
end

function __KINSetMaxNewtonStep(kinmem::KINMemPtr,mxnewtstep::realtype)
    ccall((:KINSetMaxNewtonStep,libsundials_kinsol),Cint,(KINMemPtr,realtype),kinmem,mxnewtstep)
end

function KINSetMaxNewtonStep(kinmem,mxnewtstep)
    __KINSetMaxNewtonStep(convert(KINMemPtr,kinmem),mxnewtstep)
end

function __KINSetMaxBetaFails(kinmem::KINMemPtr,mxnbcf::Clong)
    ccall((:KINSetMaxBetaFails,libsundials_kinsol),Cint,(KINMemPtr,Clong),kinmem,mxnbcf)
end

function KINSetMaxBetaFails(kinmem,mxnbcf)
    __KINSetMaxBetaFails(convert(KINMemPtr,kinmem),convert(Clong,mxnbcf))
end

function __KINSetRelErrFunc(kinmem::KINMemPtr,relfunc::realtype)
    ccall((:KINSetRelErrFunc,libsundials_kinsol),Cint,(KINMemPtr,realtype),kinmem,relfunc)
end

function KINSetRelErrFunc(kinmem,relfunc)
    __KINSetRelErrFunc(convert(KINMemPtr,kinmem),relfunc)
end

function __KINSetFuncNormTol(kinmem::KINMemPtr,fnormtol::realtype)
    ccall((:KINSetFuncNormTol,libsundials_kinsol),Cint,(KINMemPtr,realtype),kinmem,fnormtol)
end

function KINSetFuncNormTol(kinmem,fnormtol)
    __KINSetFuncNormTol(convert(KINMemPtr,kinmem),fnormtol)
end

function __KINSetScaledStepTol(kinmem::KINMemPtr,scsteptol::realtype)
    ccall((:KINSetScaledStepTol,libsundials_kinsol),Cint,(KINMemPtr,realtype),kinmem,scsteptol)
end

function KINSetScaledStepTol(kinmem,scsteptol)
    __KINSetScaledStepTol(convert(KINMemPtr,kinmem),scsteptol)
end

function __KINSetConstraints(kinmem::KINMemPtr,constraints::N_Vector)
    ccall((:KINSetConstraints,libsundials_kinsol),Cint,(KINMemPtr,N_Vector),kinmem,constraints)
end

function KINSetConstraints(kinmem,constraints)
    __constraints = convert(NVector,constraints)
    __KINSetConstraints(convert(KINMemPtr,kinmem),convert(N_Vector,__constraints))
end

function __KINSetSysFunc(kinmem::KINMemPtr,func::KINSysFn)
    ccall((:KINSetSysFunc,libsundials_kinsol),Cint,(KINMemPtr,KINSysFn),kinmem,func)
end

function KINSetSysFunc(kinmem,func)
    __KINSetSysFunc(convert(KINMemPtr,kinmem),KINSysFn_wrapper(func))
end

function __KINInit(kinmem::KINMemPtr,func::KINSysFn,tmpl::N_Vector)
    ccall((:KINInit,libsundials_kinsol),Cint,(KINMemPtr,KINSysFn,N_Vector),kinmem,func,tmpl)
end

function KINInit(kinmem,func,tmpl)
    __tmpl = convert(NVector,tmpl)
    __KINInit(convert(KINMemPtr,kinmem),KINSysFn_wrapper(func),convert(N_Vector,__tmpl))
end

function __KINSol(kinmem::KINMemPtr,uu::N_Vector,strategy::Cint,u_scale::N_Vector,f_scale::N_Vector)
    ccall((:KINSol,libsundials_kinsol),Cint,(KINMemPtr,N_Vector,Cint,N_Vector,N_Vector),kinmem,uu,strategy,u_scale,f_scale)
end

function KINSol(kinmem,uu,strategy,u_scale,f_scale)
    __uu = convert(NVector,uu)
    __u_scale = convert(NVector,u_scale)
    __f_scale = convert(NVector,f_scale)
    __KINSol(convert(KINMemPtr,kinmem),convert(N_Vector,__uu),convert(Cint,strategy),convert(N_Vector,__u_scale),convert(N_Vector,__f_scale))
end

function __KINGetWorkSpace(kinmem::KINMemPtr,lenrw::Ptr{Clong},leniw::Ptr{Clong})
    ccall((:KINGetWorkSpace,libsundials_kinsol),Cint,(KINMemPtr,Ptr{Clong},Ptr{Clong}),kinmem,lenrw,leniw)
end

function KINGetWorkSpace(kinmem,lenrw,leniw)
    __KINGetWorkSpace(convert(KINMemPtr,kinmem),pointer(lenrw),pointer(leniw))
end

function __KINGetNumNonlinSolvIters(kinmem::KINMemPtr,nniters::Ptr{Clong})
    ccall((:KINGetNumNonlinSolvIters,libsundials_kinsol),Cint,(KINMemPtr,Ptr{Clong}),kinmem,nniters)
end

function KINGetNumNonlinSolvIters(kinmem,nniters)
    __KINGetNumNonlinSolvIters(convert(KINMemPtr,kinmem),pointer(nniters))
end

function __KINGetNumFuncEvals(kinmem::KINMemPtr,nfevals::Ptr{Clong})
    ccall((:KINGetNumFuncEvals,libsundials_kinsol),Cint,(KINMemPtr,Ptr{Clong}),kinmem,nfevals)
end

function KINGetNumFuncEvals(kinmem,nfevals)
    __KINGetNumFuncEvals(convert(KINMemPtr,kinmem),pointer(nfevals))
end

function __KINGetNumBetaCondFails(kinmem::KINMemPtr,nbcfails::Ptr{Clong})
    ccall((:KINGetNumBetaCondFails,libsundials_kinsol),Cint,(KINMemPtr,Ptr{Clong}),kinmem,nbcfails)
end

function KINGetNumBetaCondFails(kinmem,nbcfails)
    __KINGetNumBetaCondFails(convert(KINMemPtr,kinmem),pointer(nbcfails))
end

function __KINGetNumBacktrackOps(kinmem::KINMemPtr,nbacktr::Ptr{Clong})
    ccall((:KINGetNumBacktrackOps,libsundials_kinsol),Cint,(KINMemPtr,Ptr{Clong}),kinmem,nbacktr)
end

function KINGetNumBacktrackOps(kinmem,nbacktr)
    __KINGetNumBacktrackOps(convert(KINMemPtr,kinmem),pointer(nbacktr))
end

function __KINGetFuncNorm(kinmem::KINMemPtr,fnorm::Ptr{realtype})
    ccall((:KINGetFuncNorm,libsundials_kinsol),Cint,(KINMemPtr,Ptr{realtype}),kinmem,fnorm)
end

function KINGetFuncNorm(kinmem,fnorm)
    __KINGetFuncNorm(convert(KINMemPtr,kinmem),pointer(fnorm))
end

function __KINGetStepLength(kinmem::KINMemPtr,steplength::Ptr{realtype})
    ccall((:KINGetStepLength,libsundials_kinsol),Cint,(KINMemPtr,Ptr{realtype}),kinmem,steplength)
end

function KINGetStepLength(kinmem,steplength)
    __KINGetStepLength(convert(KINMemPtr,kinmem),pointer(steplength))
end

function __KINGetReturnFlagName(flag::Clong)
    ccall((:KINGetReturnFlagName,libsundials_kinsol),Ptr{UInt8},(Clong,),flag)
end

function KINGetReturnFlagName(flag)
    __KINGetReturnFlagName(convert(Clong,flag))
end

function KINFree(kinmem::Ref{KINMemPtr})
    ccall((:KINFree,libsundials_kinsol),Void,(Ref{KINMemPtr},),kinmem)
end
# Julia wrapper for header: /home/astukalov/.julia/v0.5/Sundials/deps/usr/include/kinsol/kinsol_direct.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function __KINDlsSetDenseJacFn(kinmem::KINMemPtr,jac::KINDlsDenseJacFn)
    ccall((:KINDlsSetDenseJacFn,libsundials_kinsol),Cint,(KINMemPtr,KINDlsDenseJacFn),kinmem,jac)
end

function KINDlsSetDenseJacFn(kinmem,jac)
    __KINDlsSetDenseJacFn(convert(KINMemPtr,kinmem),jac)
end

function __KINDlsSetBandJacFn(kinmem::KINMemPtr,jac::KINDlsBandJacFn)
    ccall((:KINDlsSetBandJacFn,libsundials_kinsol),Cint,(KINMemPtr,KINDlsBandJacFn),kinmem,jac)
end

function KINDlsSetBandJacFn(kinmem,jac)
    __KINDlsSetBandJacFn(convert(KINMemPtr,kinmem),jac)
end

function __KINDlsGetWorkSpace(kinmem::KINMemPtr,lenrwB::Ptr{Clong},leniwB::Ptr{Clong})
    ccall((:KINDlsGetWorkSpace,libsundials_kinsol),Cint,(KINMemPtr,Ptr{Clong},Ptr{Clong}),kinmem,lenrwB,leniwB)
end

function KINDlsGetWorkSpace(kinmem,lenrwB,leniwB)
    __KINDlsGetWorkSpace(convert(KINMemPtr,kinmem),pointer(lenrwB),pointer(leniwB))
end

function __KINDlsGetNumJacEvals(kinmem::KINMemPtr,njevalsB::Ptr{Clong})
    ccall((:KINDlsGetNumJacEvals,libsundials_kinsol),Cint,(KINMemPtr,Ptr{Clong}),kinmem,njevalsB)
end

function KINDlsGetNumJacEvals(kinmem,njevalsB)
    __KINDlsGetNumJacEvals(convert(KINMemPtr,kinmem),pointer(njevalsB))
end

function __KINDlsGetNumFuncEvals(kinmem::KINMemPtr,nfevalsB::Ptr{Clong})
    ccall((:KINDlsGetNumFuncEvals,libsundials_kinsol),Cint,(KINMemPtr,Ptr{Clong}),kinmem,nfevalsB)
end

function KINDlsGetNumFuncEvals(kinmem,nfevalsB)
    __KINDlsGetNumFuncEvals(convert(KINMemPtr,kinmem),pointer(nfevalsB))
end

function __KINDlsGetLastFlag(kinmem::KINMemPtr,flag::Ptr{Clong})
    ccall((:KINDlsGetLastFlag,libsundials_kinsol),Cint,(KINMemPtr,Ptr{Clong}),kinmem,flag)
end

function KINDlsGetLastFlag(kinmem,flag)
    __KINDlsGetLastFlag(convert(KINMemPtr,kinmem),pointer(flag))
end

function __KINDlsGetReturnFlagName(flag::Clong)
    ccall((:KINDlsGetReturnFlagName,libsundials_kinsol),Ptr{UInt8},(Clong,),flag)
end

function KINDlsGetReturnFlagName(flag)
    __KINDlsGetReturnFlagName(convert(Clong,flag))
end
# Julia wrapper for header: /home/astukalov/.julia/v0.5/Sundials/deps/usr/include/kinsol/kinsol_spils.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function __KINSpilsSetMaxRestarts(kinmem::KINMemPtr,maxrs::Cint)
    ccall((:KINSpilsSetMaxRestarts,libsundials_kinsol),Cint,(KINMemPtr,Cint),kinmem,maxrs)
end

function KINSpilsSetMaxRestarts(kinmem,maxrs)
    __KINSpilsSetMaxRestarts(convert(KINMemPtr,kinmem),convert(Cint,maxrs))
end

function __KINSpilsSetPreconditioner(kinmem::KINMemPtr,pset::KINSpilsPrecSetupFn,psolve::KINSpilsPrecSolveFn)
    ccall((:KINSpilsSetPreconditioner,libsundials_kinsol),Cint,(KINMemPtr,KINSpilsPrecSetupFn,KINSpilsPrecSolveFn),kinmem,pset,psolve)
end

function KINSpilsSetPreconditioner(kinmem,pset,psolve)
    __KINSpilsSetPreconditioner(convert(KINMemPtr,kinmem),pset,psolve)
end

function __KINSpilsSetJacTimesVecFn(kinmem::KINMemPtr,jtv::KINSpilsJacTimesVecFn)
    ccall((:KINSpilsSetJacTimesVecFn,libsundials_kinsol),Cint,(KINMemPtr,KINSpilsJacTimesVecFn),kinmem,jtv)
end

function KINSpilsSetJacTimesVecFn(kinmem,jtv)
    __KINSpilsSetJacTimesVecFn(convert(KINMemPtr,kinmem),jtv)
end

function __KINSpilsGetWorkSpace(kinmem::KINMemPtr,lenrwSG::Ptr{Clong},leniwSG::Ptr{Clong})
    ccall((:KINSpilsGetWorkSpace,libsundials_kinsol),Cint,(KINMemPtr,Ptr{Clong},Ptr{Clong}),kinmem,lenrwSG,leniwSG)
end

function KINSpilsGetWorkSpace(kinmem,lenrwSG,leniwSG)
    __KINSpilsGetWorkSpace(convert(KINMemPtr,kinmem),pointer(lenrwSG),pointer(leniwSG))
end

function __KINSpilsGetNumPrecEvals(kinmem::KINMemPtr,npevals::Ptr{Clong})
    ccall((:KINSpilsGetNumPrecEvals,libsundials_kinsol),Cint,(KINMemPtr,Ptr{Clong}),kinmem,npevals)
end

function KINSpilsGetNumPrecEvals(kinmem,npevals)
    __KINSpilsGetNumPrecEvals(convert(KINMemPtr,kinmem),pointer(npevals))
end

function __KINSpilsGetNumPrecSolves(kinmem::KINMemPtr,npsolves::Ptr{Clong})
    ccall((:KINSpilsGetNumPrecSolves,libsundials_kinsol),Cint,(KINMemPtr,Ptr{Clong}),kinmem,npsolves)
end

function KINSpilsGetNumPrecSolves(kinmem,npsolves)
    __KINSpilsGetNumPrecSolves(convert(KINMemPtr,kinmem),pointer(npsolves))
end

function __KINSpilsGetNumLinIters(kinmem::KINMemPtr,nliters::Ptr{Clong})
    ccall((:KINSpilsGetNumLinIters,libsundials_kinsol),Cint,(KINMemPtr,Ptr{Clong}),kinmem,nliters)
end

function KINSpilsGetNumLinIters(kinmem,nliters)
    __KINSpilsGetNumLinIters(convert(KINMemPtr,kinmem),pointer(nliters))
end

function __KINSpilsGetNumConvFails(kinmem::KINMemPtr,nlcfails::Ptr{Clong})
    ccall((:KINSpilsGetNumConvFails,libsundials_kinsol),Cint,(KINMemPtr,Ptr{Clong}),kinmem,nlcfails)
end

function KINSpilsGetNumConvFails(kinmem,nlcfails)
    __KINSpilsGetNumConvFails(convert(KINMemPtr,kinmem),pointer(nlcfails))
end

function __KINSpilsGetNumJtimesEvals(kinmem::KINMemPtr,njvevals::Ptr{Clong})
    ccall((:KINSpilsGetNumJtimesEvals,libsundials_kinsol),Cint,(KINMemPtr,Ptr{Clong}),kinmem,njvevals)
end

function KINSpilsGetNumJtimesEvals(kinmem,njvevals)
    __KINSpilsGetNumJtimesEvals(convert(KINMemPtr,kinmem),pointer(njvevals))
end

function __KINSpilsGetNumFuncEvals(kinmem::KINMemPtr,nfevalsS::Ptr{Clong})
    ccall((:KINSpilsGetNumFuncEvals,libsundials_kinsol),Cint,(KINMemPtr,Ptr{Clong}),kinmem,nfevalsS)
end

function KINSpilsGetNumFuncEvals(kinmem,nfevalsS)
    __KINSpilsGetNumFuncEvals(convert(KINMemPtr,kinmem),pointer(nfevalsS))
end

function __KINSpilsGetLastFlag(kinmem::KINMemPtr,flag::Ptr{Clong})
    ccall((:KINSpilsGetLastFlag,libsundials_kinsol),Cint,(KINMemPtr,Ptr{Clong}),kinmem,flag)
end

function KINSpilsGetLastFlag(kinmem,flag)
    __KINSpilsGetLastFlag(convert(KINMemPtr,kinmem),pointer(flag))
end

function __KINSpilsGetReturnFlagName(flag::Clong)
    ccall((:KINSpilsGetReturnFlagName,libsundials_kinsol),Ptr{UInt8},(Clong,),flag)
end

function KINSpilsGetReturnFlagName(flag)
    __KINSpilsGetReturnFlagName(convert(Clong,flag))
end
# Julia wrapper for header: /home/astukalov/.julia/v0.5/Sundials/deps/usr/include/kinsol/kinsol_band.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function __KINBand(kinmem::KINMemPtr,N::Clong,mupper::Clong,mlower::Clong)
    ccall((:KINBand,libsundials_kinsol),Cint,(KINMemPtr,Clong,Clong,Clong),kinmem,N,mupper,mlower)
end

function KINBand(kinmem,N,mupper,mlower)
    __KINBand(convert(KINMemPtr,kinmem),convert(Clong,N),convert(Clong,mupper),convert(Clong,mlower))
end
# Julia wrapper for header: /home/astukalov/.julia/v0.5/Sundials/deps/usr/include/kinsol/kinsol_bbdpre.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function __KINBBDPrecInit(kinmem::KINMemPtr,Nlocal::Clong,mudq::Clong,mldq::Clong,mukeep::Clong,mlkeep::Clong,dq_rel_uu::realtype,gloc::KINLocalFn,gcomm::KINCommFn)
    ccall((:KINBBDPrecInit,libsundials_kinsol),Cint,(KINMemPtr,Clong,Clong,Clong,Clong,Clong,realtype,KINLocalFn,KINCommFn),kinmem,Nlocal,mudq,mldq,mukeep,mlkeep,dq_rel_uu,gloc,gcomm)
end

function KINBBDPrecInit(kinmem,Nlocal,mudq,mldq,mukeep,mlkeep,dq_rel_uu,gloc,gcomm)
    __KINBBDPrecInit(convert(KINMemPtr,kinmem),convert(Clong,Nlocal),convert(Clong,mudq),convert(Clong,mldq),convert(Clong,mukeep),convert(Clong,mlkeep),dq_rel_uu,gloc,gcomm)
end

function __KINBBDPrecGetWorkSpace(kinmem::KINMemPtr,lenrwBBDP::Ptr{Clong},leniwBBDP::Ptr{Clong})
    ccall((:KINBBDPrecGetWorkSpace,libsundials_kinsol),Cint,(KINMemPtr,Ptr{Clong},Ptr{Clong}),kinmem,lenrwBBDP,leniwBBDP)
end

function KINBBDPrecGetWorkSpace(kinmem,lenrwBBDP,leniwBBDP)
    __KINBBDPrecGetWorkSpace(convert(KINMemPtr,kinmem),pointer(lenrwBBDP),pointer(leniwBBDP))
end

function __KINBBDPrecGetNumGfnEvals(kinmem::KINMemPtr,ngevalsBBDP::Ptr{Clong})
    ccall((:KINBBDPrecGetNumGfnEvals,libsundials_kinsol),Cint,(KINMemPtr,Ptr{Clong}),kinmem,ngevalsBBDP)
end

function KINBBDPrecGetNumGfnEvals(kinmem,ngevalsBBDP)
    __KINBBDPrecGetNumGfnEvals(convert(KINMemPtr,kinmem),pointer(ngevalsBBDP))
end
# Julia wrapper for header: /home/astukalov/.julia/v0.5/Sundials/deps/usr/include/kinsol/kinsol_dense.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function __KINDense(kinmem::KINMemPtr,N::Clong)
    ccall((:KINDense,libsundials_kinsol),Cint,(KINMemPtr,Clong),kinmem,N)
end

function KINDense(kinmem,N)
    __KINDense(convert(KINMemPtr,kinmem),convert(Clong,N))
end
# Julia wrapper for header: /home/astukalov/.julia/v0.5/Sundials/deps/usr/include/kinsol/kinsol_impl.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

# Julia wrapper for header: /home/astukalov/.julia/v0.5/Sundials/deps/usr/include/kinsol/kinsol_spbcgs.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function __KINSpbcg(kinmem::KINMemPtr,maxl::Cint)
    ccall((:KINSpbcg,libsundials_kinsol),Cint,(KINMemPtr,Cint),kinmem,maxl)
end

function KINSpbcg(kinmem,maxl)
    __KINSpbcg(convert(KINMemPtr,kinmem),convert(Cint,maxl))
end
# Julia wrapper for header: /home/astukalov/.julia/v0.5/Sundials/deps/usr/include/kinsol/kinsol_spgmr.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function __KINSpgmr(kinmem::KINMemPtr,maxl::Cint)
    ccall((:KINSpgmr,libsundials_kinsol),Cint,(KINMemPtr,Cint),kinmem,maxl)
end

function KINSpgmr(kinmem,maxl)
    __KINSpgmr(convert(KINMemPtr,kinmem),convert(Cint,maxl))
end
# Julia wrapper for header: /home/astukalov/.julia/v0.5/Sundials/deps/usr/include/kinsol/kinsol_sptfqmr.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function __KINSptfqmr(kinmem::KINMemPtr,maxl::Cint)
    ccall((:KINSptfqmr,libsundials_kinsol),Cint,(KINMemPtr,Cint),kinmem,maxl)
end

function KINSptfqmr(kinmem,maxl)
    __KINSptfqmr(convert(KINMemPtr,kinmem),convert(Cint,maxl))
end


text:
# Autogenerated wrapper script for libgeotiff_jll for x86_64-w64-mingw32
export applygeo, geotifcp, libgeotiff, listgeo, makegeo

using PROJ_jll
using Libtiff_jll
using LibCURL_jll
JLLWrappers.@generate_wrapper_header("libgeotiff")
JLLWrappers.@declare_library_product(libgeotiff, "libgeotiff.dll")
JLLWrappers.@declare_executable_product(applygeo)
JLLWrappers.@declare_executable_product(geotifcp)
JLLWrappers.@declare_executable_product(listgeo)
JLLWrappers.@declare_executable_product(makegeo)
function __init__()
    JLLWrappers.@generate_init_header(PROJ_jll, Libtiff_jll, LibCURL_jll)
    JLLWrappers.@init_library_product(
        libgeotiff,
        "bin\\libgeotiff.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        applygeo,
        "bin\\applygeo.exe",
    )

    JLLWrappers.@init_executable_product(
        geotifcp,
        "bin\\geotifcp.exe",
    )

    JLLWrappers.@init_executable_product(
        listgeo,
        "bin\\listgeo.exe",
    )

    JLLWrappers.@init_executable_product(
        makegeo,
        "bin\\makegeo.exe",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Test

@testset "OCReractTest" begin
    include("test_tesseract.jl")
end


text:
module CastepBin
using FortranFiles
CastepFortranFile(io) = FortranFile(io, convert="big-endian")
const IntF = Int32

"""
Macro for compact representation of how data should read from the binary file

For example, the following definition:
```julia
@readtag f outdict tags begin

    CELL_VERSION_NUMBER
    VERSION_NUMBER::Float64
    skip
    X::(Float64, :B)

end
```
expands into to:
```julia
if "CELL_VERSION_NUMBER" in keys(tags)
    gototag("CELL_VERSION_NUMBER", f, tags)
    output[:VERSION_NUMBER] = read(f, Float64)
    skiptag(f)
    output[:X] = read(f, outdict[:B])
end

Note that tags that are not present will simply be skipped.
```

"""
macro readtag(f, outdict, tags, expr)
    @assert expr.head == :block "Expect a block definition"
    eout = []
    thistag = Pair(nothing, [])

    function flush_if_block!(eout, thistag)
        if !isnothing(thistag.first)
            push!(eout, Expr(:if,
            quote
                $(thistag.first) in keys($tags)
            end,
            Expr(:block, thistag.second...))
            )
        end
    end

    for line in expr.args
        isa(line, LineNumberNode) && continue
        # This defines a new tags
        # TODO - add condition to check if the tag exists or not
        # Do nothing is the tags is not there!!!
        # Symbol as tag name
        if isa(line, Symbol)
            if line == :skip || line == :SKIP
                push!(thistag.second, :(skiptag($f)))
                continue
            else
                line = String(line)
            end
        end
        # String - treat as the name of the required tag
        if isa(line, String)
            # Create the condition - only parse if the tag exists
            if !isnothing(thistag.first)
                flush_if_block!(eout, thistag)
            end
            thistag = Pair(line, [])
            push!(thistag.second, quote
                gototag($line, $f, $tags)
            end)
            continue
        end
        # This is a read line
        if line.head == :(::)
            var = line.args[1]
            type = line.args[2]
            # Parse the existing parameters
            # Check if any parameters refer to exiting ones
            # X::(Float64, :A)
            # Expands to read(f, (Float64, outdict[:A]))
            if isa(type, Expr) && type.head == :tuple
                type_args = []
                # Construct the type field
                for arg in  type.args
                    if isa(arg, QuoteNode)
                        push!(type_args, quote
                            $outdict[$arg]
                        end
                        )
                    else
                        push!(type_args, arg)
                    end
                end
                type = Expr(:tuple, type_args...)
            end
            push!(thistag.second, quote
                $outdict[$(Meta.quot(var))] = read($f, $type)
            end)
        end
    end

    # Flust out the last entry
    if !isnothing(thistag.first)
        flush_if_block!(eout, thistag)
    end

    esc(Expr(:block, eout...))
end

"""
Scan the binary file for tags
"""
function locate_tags(f::FortranFile;max_length=256)
    REC_LENGTH=4
    out = Dict{String, Int}()
    while true
        # All patterns has been found? just break the loop
        local rec
        try
            rec = FortranFiles.Record(f)
        catch y
            if isa(y, EOFError)
                break
            end
        end
        rec.subreclen <= max_length && record_tag!(out, rec;offset=REC_LENGTH)
        close(rec)
    end
    out
end

function locate_tags(name::AbstractString)
    open(name) do f
        locate_tags(CastepFortranFile(f))
    end
end

"""
    record_tag!(out::Dict{String,Int}, rec;offset::Int)

Record the positions of eligible records
"""
function record_tag!(out::Dict{String,Int}, rec;offset::Int)
    loc = position(rec.io)
    fs = read(rec, FString{rec.subreclen})  # Read as string
    any(x -> x < 0, fs.data) && return nothing
    str = strip(convert(String, fs))
    m = match(r"[A-Z_0-9]*$", str)
    isnothing(m) && return
    if str in keys(out)
        str = str * "_SECOND"
    end
    out[str] = loc - offset
end

"""
    read_cell!(output, f::FortranFile, tags::Dict{String,Int})

Read in the unit cell information. Note that the data is taken from the *second* set which stores
the *current* cell, e.g. consistent with the wavefunction data stored.
"""
function read_cell!(output, f::FortranFile, tags::Dict{String,Int})

    @readtag f output tags begin
        BEGIN_UNIT_CELL_SECOND
        skip
        REAL_LATTICE::(Float64, 3, 3)
        skip
        RECIP_LATTICE::(Float64, 3, 3)
        skip
        VOLUME::Float64
        skip
        NUM_SPECIES::IntF
        skip
        NUM_IONS::IntF
        skip
        MAX_IONS_IN_SPECIES::IntF

        CELL_VERSION_NUMBER_SECOND
        VERSION_NUMBER::Float64

        "CELL%NUM_IONS_IN_SPECIES_SECOND"

        NUM_IONS_IN_SPECIES::(IntF, :NUM_SPECIES)

        "CELL%IONIC_POSITIONS_SECOND"

        IONIC_POSITIONS::(Float64, 3, :MAX_IONS_IN_SPECIES, :NUM_SPECIES)

        "CELL%IONIC_VELOCITIES_SECOND"

        IONIC_VELOCITIES::(Float64, 3, :MAX_IONS_IN_SPECIES, :NUM_SPECIES)

        "CELL%SPECIES_SYMBOL_SECOND"

        SPECIES_SYMBOL::(FString{8}, :NUM_SPECIES)
    end

    @readtag f output tags begin
        NKPTS_SECOND
        NKPTS::IntF

        KPOINTS_SECOND
        KPOINTS::(Float64, 3, :NKPTS)

        KPOINT_WEIGHTS_SECOND
        KPOINT_WEIGHTS::(Float64, :NKPTS)
    end

    output
end

"""
    read_parameters_electronic(output, f::FortranFile, tags::Dict{String, Int})

Read the data store after `BEGIN_ELECTRONIC` header. These are internal parameters used by
CASTEP.
There are duplicated entries in the same section as some of the field are written as different
types for backward compatibility.
"""
function read_parameters_electronic!(output, f::FortranFile, tags::Dict{String, Int})
    @readtag f output tags begin
        BEGIN_ELECTRONIC
        skip
        skip
        skip
        skip
        PARAM_NBANDS::IntF
        PARAM_ELEC_TEMP::Float64
        skip
        skip
        PARAM_SPIN_POLARIZED::IntF
        PARAM_ELECTRONIC_MINIMIZER::FString{10}
        PARAM_NELECTRONS::Float64
        PARAM_NUP::Float64
        PARAM_NDOWN::Float64
        PARAM_SPIN::Float64
        PARAM_CHARGE::Float64
    end
end

function read_grid_properties!(output, f::FortranFile, tags::Dict{String, Int}, with_wavefunction=false)

    @readtag f output tags begin
        END_CELL_GLOBAL_SECOND
        FOUND_GROUND_STATE_WAVEFUNCTION::IntF
        FOUND_GROUND_STATE_DENSITY::IntF
        TOTAL_ENERGY::Float64
        FERMI_ENERGY::Float64
    end
    nbands, nspins = read(f, IntF, IntF)
    output[:NBANDS] = nbands
    output[:NSPINS] = nspins
    if with_wavefunction
        read_wavefunction_complex!(output, f, tags)
    end
    read_eigenvalue_and_occ!(output, f, tags)
    skiptag(f)
    ngxf, ngyf, ngzf = read(f, IntF, IntF, IntF)
    output[:NGX_FINE] = ngxf
    output[:NGY_FINE] = ngyf
    output[:NGZ_FINE] = ngzf
    read_charge_density!(output, f, tags)
end

"""
Read eigenvalues and occupations
"""
function read_eigenvalue_and_occ!(output, f, tags)
    nkpts::Int = output[:NKPTS]
    nspin::Int = output[:NSPINS]
    nbands::Int = output[:NBANDS]

    kpoints = zeros(3, nkpts)
    occ = zeros(nbands, nkpts, nspin)
    eignvalues = zeros(nbands, nkpts, nspin)
    for ik in 1:nkpts
        kpoints[:, ik] = read(f, (Float64, 3))
        for is in 1:nspin
            occ[:, ik, is] = read(f, (Float64, nbands))
            eignvalues[:, ik, is] = read(f, (Float64, nbands))
        end
    end
    output[:OCCUPATIONS] = occ
    output[:EIGENVALUES] = eignvalues
    output[:KPOINTS_OF_EIGENVALUES] = kpoints
end

"""
Read the charge density section from the binary file

NOTE: only support linear spin-polarisation for now
"""
function read_charge_density!(output, f, tags)
    ngxf::Int = output[:NGX_FINE]
    ngyf::Int = output[:NGY_FINE]
    ngzf::Int = output[:NGZ_FINE]
    nspin::Int = output[:NSPINS]
    # NOTE - assume not NCM
    spin_density = zeros(ComplexF64, ngxf, ngyf, ngzf)
    charge_density = zeros(ComplexF64, ngxf, ngyf, ngzf)

    for _ in 1:ngxf
        for _ in 1:ngyf
            if nspin == 2
                nx, ny, zcol, spincol = read(f, IntF, IntF, (ComplexF64, ngzf), (ComplexF64, ngzf))
                charge_density[nx, ny, :] = zcol
                spin_density[nx, ny, :] = spincol
            else
                nx, ny, zcol = read(f, IntF, IntF, (ComplexF64, ngzf))
                charge_density[nx, ny, :] = zcol
            end
        end
    end
    if nspin == 2
        output[:SPIN_DENSITY] = spin_density
    end
    output[:CHARGE_DENSITY] = charge_density
end

function read_wavefunction_complex!(output, f, tags)
    header = read(f, FString{4})
    ngx, ngy, ngz = read(f, IntF, IntF, IntF)
    coeff_size_1, spinorcomps, nbands_max, nkpts, nspins = read(f, IntF, IntF, IntF, IntF, IntF)
    coeffs = zeros(ComplexF64, coeff_size_1, spinorcomps, nbands_max, nkpts, nspins)
    nwaves_at_kp = zeros(Int, nkpts)
    kpts = zeros(3, nkpts)
    pw_grid_coord = zeros(IntF, 3, coeff_size_1, nkpts)

    for is in 1:nspins
        for ik in 1:nkpts
            kpts[:, ik], nwaves = read(f, (Float64, 3), IntF)
            nwaves_at_kp[ik] = nwaves

            # Grid coordinates
            coords_x = read(f, (IntF, nwaves))
            coords_y = read(f, (IntF, nwaves))
            coords_z = read(f, (IntF, nwaves))
            pw_grid_coord[1, 1:nwaves, ik] = coords_x
            pw_grid_coord[2, 1:nwaves, ik] = coords_y
            pw_grid_coord[3, 1:nwaves, ik] = coords_z

            for ib in 1:nbands_max
                for ispinor in 1:spinorcomps
                    coeff = read(f, (ComplexF64, nwaves))
                    coeffs[1:nwaves, ispinor, ib, ik, is] = coeff
                end
            end
        end
    end
    data = (
        ngx=ngx,
        ngy=ngy,
        ngz=ngz,
        pw_grid_coord=pw_grid_coord,
        coeffs=coeffs,
        kpts=kpts,
        nwaves_at_kp=nwaves_at_kp,
    )
    output[:WAVEFUNCTION_DATA] = data
end

"""
Write the wavefunction coefficients back to the file
"""
function write_wavefunction_complex!(output, f, tags, coeffs)
    header = read(f, FString{4})
    ngx, ngy, ngz = read(f, IntF, IntF, IntF)
    coeff_size_1, spinorcomps, nbands_max, nkpts, nspins = read(f, IntF, IntF, IntF, IntF, IntF)
    nwaves_at_kp = zeros(Int, nkpts)
    kpts = zeros(3, nkpts)
    pw_grid_coord = zeros(IntF, 3, coeff_size_1, nkpts)
    for is in 1:nspins
        for ik in 1:nkpts
            kpts[:, ik], nwaves = read(f, (Float64, 3), IntF)
            nwaves_at_kp[ik] = nwaves

            # Grid coordinates
            coords_x = read(f, (IntF, nwaves))
            coords_y = read(f, (IntF, nwaves))
            coords_z = read(f, (IntF, nwaves))
            pw_grid_coord[1, 1:nwaves, ik] = coords_x
            pw_grid_coord[2, 1:nwaves, ik] = coords_y
            pw_grid_coord[3, 1:nwaves, ik] = coords_z

            for ib in 1:nbands_max
                for ispinor in 1:spinorcomps
                    write(f, coeffs[1:nwaves, ispinor, ib, ik, is])
                end
            end
        end
    end
    data = (
        ngx=ngx,
        ngy=ngy,
        ngz=ngz,
        pw_grid_coord=pw_grid_coord,
        coeffs=coeffs,
        kpts=kpts,
        nwaves_at_kp=nwaves_at_kp,
    )
    output[:WAVEFUNCTION_DATA] = data
end

"""
Read the forces
"""
function read_forces!(output, f::FortranFile, tags::Dict{String,Int})
    @readtag f output tags begin
        FORCES
        FORCES::(Float64, 3, :MAX_IONS_IN_SPECIES, :NUM_SPECIES)
    end
    output
end

function read_stress!(output, f::FortranFile, tags::Dict{String,Int})
    @readtag f output tags begin
        STRESS
        STRESS::(Float64, 6)
        STRAIN::(Float64, 3,3)
    end
end

"""
Seek to a tag and be read to read in the data following it
"""
function gototag(tag, f::FortranFile, tags::Dict{String,Int})
    seek(f.io, tags[string(tag)])
    rec = FortranFiles.Record(f)
    close(rec)
end

"""Skip to the next tag"""
function skiptag(f)
    rec = FortranFiles.Record(f)
    close(rec)
end

"""
READ RAW data from CASTEP binary file

Expect everyting to be in the atomic units
Atoms are indexed by (species, index in specie)
"""
function read_castep_check(fname)
    output = Dict{Symbol, Any}()
    open(fname) do fhandle
        f = CastepFortranFile(fhandle)
        tags = locate_tags(f)
        read_parameters_electronic!(output, f, tags)
        read_cell!(output, f, tags)
        read_forces!(output, f, tags)
        read_stress!(output, f, tags)
        has_wave = !("CASTEP_BIN" in keys(tags))
        read_grid_properties!(output, f, tags, has_wave)
    end
    output
end

"""
Read the grid properties but update the plane wave coefficients as given
"""
function update_wavefunction!(output, f::FortranFile, tags::Dict{String, Int}, coeffs)

    @readtag f output tags begin
        END_CELL_GLOBAL_SECOND
        FOUND_GROUND_STATE_WAVEFUNCTION::IntF
        FOUND_GROUND_STATE_DENSITY::IntF
        TOTAL_ENERGY::Float64
        FERMI_ENERGY::Float64
    end
    nbands, nspins = read(f, IntF, IntF)
    output[:NBANDS] = nbands
    output[:NSPINS] = nspins
    # Write the wave function
    write_wavefunction_complex!(output, f, tags, coeffs)
end

"""
    update_wavefunction(fname, coeffs::AbstractArray{T, 5}) where {T}

Update the wavefunction (coefficients) stored in the check file.
The coefficients array passed should be modified based on that was originally
read from the file.

* **NOTE** this function will overwrite target file.
"""
function update_wavefunction!(fname, coeffs::AbstractArray{T, 5}) where {T}
    output = Dict{Symbol, Any}()
    open(fname, "r+") do fhandle
        f = CastepFortranFile(fhandle)
        tags = locate_tags(f)
        has_wave = !("CASTEP_BIN" in keys(tags))
        @assert has_wave "The file does not contain any wavefunction data"
        update_wavefunction!(output, f, tags, coeffs)
    end
end

"""
    update_wavefunction(src, dst, coeffs::AbstractArray{T, 5}) where {T}

Update the wavefunction (coefficients) stored in the check file.
The coefficients array passed should be modified based on that was originally
read from the file.
"""
function update_wavefunction(src, dst, coeffs::AbstractArray{T, 5}) where {T}
    cp(src, dst)
    update_wavefunction!(dst, coeffs)
end

precompile(read_castep_check, (String, ))
export read_castep_check, update_wavefunction, update_wavefunction!

end

using .CastepBin
import .CastepBin


text:
### A Pluto.jl notebook ###
# v0.15.1

using Markdown
using InteractiveUtils

# ╔═╡ 97539bb6-c0b0-11eb-2d83-316749d179f1
begin
	first_setup = false

	using Pkg
	Pkg.activate(".")
	if first_setup
		Pkg.add([
			"Pluto", "PlutoUI", "PlutoSliderServer",
			"Flux",
			"MLJ", "MLJFlux",
			"DifferentialEquations", "DiffEqFlux",
			"MLDatasets", "MLBase",
			"Plots", "PlotlyJS",
			"ProgressMeter",
			"Statistics",
		])
		Pkg.update()
		Pkg.build()
		Pkg.precompile()
	else
		Pkg.instantiate()
	end

	using
		Pluto, PlutoUI, PlutoSliderServer,
		Flux,
		MLJ, MLJFlux,
		DifferentialEquations, DiffEqFlux,
		MLDatasets, MLBase,
		Plots, PlotlyJS,
		ProgressMeter,
		Statistics
end

# ╔═╡ baab9408-fe06-436d-8158-6ec55e1b49cf
begin
	# MLJ.default_resource(CUDALibs())
	# MLJ.default_resource(CPUProcesses())
	# MLJ.default_resource(CPUThreads())
	MLJ.default_resource(CPU1())

	plotlyjs()
end

# ╔═╡ 5bafc48d-eb03-4625-8d87-17170133d999
begin
	just_test = true
	n_sam = 100

	use_gpu = false
end

# ╔═╡ e593f5a3-9c57-42e7-9ba5-fcf4f0c4bb76
begin
	train_x, train_y = MNIST.traindata()
	test_x,  test_y  = MNIST.testdata()

	if just_test
		train_x = coerce(train_x, GrayImage)[1:n_sam]
		train_y = coerce(train_y, Multiclass)[1:n_sam]
		test_x = coerce(test_x, GrayImage)[1:n_sam]
		test_y = coerce(test_y, Multiclass)[1:n_sam]
	else
		train_x = coerce(train_x, GrayImage)
		train_y = coerce(train_y, Multiclass)
		test_x = coerce(test_x, GrayImage)
		test_y = coerce(test_y, Multiclass)
	end

	@assert scitype(train_x) <: AbstractVector{<:GrayImage}
	@assert scitype(train_y) <: AbstractVector{<:Multiclass}
	@assert scitype(test_x) <: AbstractVector{<:GrayImage}
	@assert scitype(test_y) <: AbstractVector{<:Multiclass}
end

# ╔═╡ 3d5dac10-a6a9-4e13-b6c7-85c86d5f61ac
begin
	mutable struct ImgShortDense <: MLJFlux.Builder
		n_hidden_l::Int
		dropout::Float64
		σ
	end
	ImgShortDense(; n_hidden_l=1, dropout=1/8, σ=tanh) = ImgShortDense(n_hidden_l, dropout, σ)
	function MLJFlux.build(builder::ImgShortDense, rng, n, m, c)
		n = prod(n)

		ni_p = floor(Int, log2(n))
		no_p = ceil(Int, log2(m))
		n_middles = [2 ^ i for i in no_p:ni_p]

		imgsize_afterconvs = (n, 1)

		layers = [
		    # x -> float(cat(x..., dims=4)),
			Flux.flatten,
		    BatchNorm(prod(imgsize_afterconvs), builder.σ),
            Dropout(builder.dropout),
			Flux.Dense(prod(imgsize_afterconvs), n_middles[builder.n_hidden_l], builder.σ),
			BatchNorm(n_middles[builder.n_hidden_l], builder.σ),
            Dropout(builder.dropout),
		]
		for i in reverse(2:builder.n_hidden_l)
		    layers = vcat(layers, [
		        Flux.Dense(n_middles[i], n_middles[i-1], builder.σ),
			    BatchNorm(n_middles[i-1], builder.σ),
                Dropout(builder.dropout),
		    ])
	    end
	    layers = vcat(layers, [
	        Flux.Dense(n_middles[1], m, builder.σ),
		    BatchNorm(m, builder.σ),
            Dropout(builder.dropout),
	    ])

		return Flux.Chain(layers...) |> f64
	end
end

# ╔═╡ b48567da-38c0-4d01-b866-08413433d115
begin
	mutable struct ImgShortLite <: MLJFlux.Builder
		σ
	end
	ImgShortLite(; σ=tanh) = ImgShortLite(σ)
	function MLJFlux.build(builder::ImgShortLite, rng, n, m, c)
		n = prod(n)
		imgsize_afterconvs = (13, 13, 16)

		return Flux.Chain(
			# x -> float(cat(x..., dims=4)),
			Conv((3, 3), c => 16, builder.σ), # (26, 26, 16, N)
			MaxPool((2, 2)), # (13, 13, 16, N)
			Flux.flatten,
			Flux.Dense(prod(imgsize_afterconvs), m, builder.σ),
		) |> f64
	end
end

# ╔═╡ 85bb6d5a-5117-4b24-beb1-c1b00f0e9462
begin
	mutable struct ImgShortConv <: MLJFlux.Builder
		dropout::Float64
		σ
	end
	ImgShortConv(; dropout=1/8, σ=tanh) = ImgShortConv(dropout, σ)
	function MLJFlux.build(builder::ImgShortConv, rng, n, m, c)
		return Flux.Chain(
			# x -> float(cat(x..., dims=4)),
			Conv((3, 3), c => 16, builder.σ), # (26, 26, 16, N)
			MaxPool((2, 2)), # (13, 13, 16, N)
			Conv((3, 3), 16 => 32, builder.σ), # (11, 11, 32, N)
			MaxPool((2, 2)), # (5, 5, 32, N)
			Conv((3, 3), 32 => 32, builder.σ), # (3, 3, 32, N)
			MaxPool((2, 2)), # (1, 1, 32, N)
			Flux.flatten,
			BatchNorm(32, builder.σ),
			Dropout(builder.dropout),
			Dense(32, 64, builder.σ),
			BatchNorm(64, builder.σ),
			Dropout(builder.dropout),
			Dense(64, 32, builder.σ),
			BatchNorm(32, builder.σ),
			Dropout(builder.dropout),
			NeuralODE(Chain(
				# BatchNorm(32, builder.σ),
				Dropout(builder.dropout),
				Dense(32, 16, builder.σ),
				# BatchNorm(16, builder.σ),
				Dropout(builder.dropout),
				Dense(16, 16, builder.σ),
				# BatchNorm(16, builder.σ),
				Dropout(builder.dropout),
				Dense(16, 32, builder.σ),
				# BatchNorm(32, builder.σ),
				Dropout(builder.dropout),
			) |> f64, 0.0:1.0, Tsit5()),
			x -> first(x.u),
		    BatchNorm(32, builder.σ),
            Dropout(builder.dropout),
			Flux.Dense(32, 16, builder.σ),
			BatchNorm(16, builder.σ),
            Dropout(builder.dropout),
			Flux.Dense(16, m, builder.σ),
		    BatchNorm(m, builder.σ),
            Dropout(builder.dropout),
		) |> f64
	end
end

# ╔═╡ 984cbe16-be12-4fdf-8345-65b601b2c4cd
begin
	itd = ImgShortDense()
	mitd = MLJFlux.build(itd, nothing, (28, 28), 10, 1)
	mitd(float(cat(train_x[1:10]..., dims=4)))
end

# ╔═╡ 4a988da3-1a6c-4bb0-a15c-779e2f2b5cee
begin
	itc = ImgShortConv()
	mitc = MLJFlux.build(itc, nothing, (28, 28), 10, 1)
	mitc(float(cat(train_x[1:10]..., dims=4)))
end

# ╔═╡ 03dd99ff-edf2-4c66-ac15-a4dff9b2c7c6
begin
	itl = ImgShortLite()
	mitl = MLJFlux.build(itc, nothing, (28, 28), 10, 1)
	mitl(float(cat(train_x[1:10]..., dims=4)))
end

# ╔═╡ 58ffa080-cf63-43d9-bc40-d9fedb085755
const ImgShort = ImgShortConv

# ╔═╡ d0e5c847-6bd5-484d-9e62-24deb3a8f0a5
begin
	losses = []
	training_losses = []

	add_loss(loss) = push!(losses, loss)
	add_training_loss(losses) = push!(training_losses, losses[end])

	f_mdl = IteratedModel(
		model=ImageClassifier(
			builder=ImgShort(),
			finaliser=softmax,
			optimiser=ADAM(),
			loss=Flux.crossentropy,
			epochs=2^4,
			batch_size=2^12,
			acceleration=(use_gpu ? CUDALibs() : CPU1()),
		),
		controls=[
			Step(),
			TimeLimit(),
			NumberLimit(),
			NumberSinceBest(),
			Threshold(),
			GL(),
			PQ(),
			Patience(),

			# NotANumber(),
			# InvalidValue(),

			# Info(),
			# Warn(Bool),
			# Error(Bool),
			Callback(),
			WithNumberDo(),
			WithIterationsDo(),
			WithLossDo(),
			WithTrainingLossesDo(),

			WithLossDo(add_loss),
			WithTrainingLossesDo(add_training_loss),
			# Save(joinpath(@__DIR__, "mnist_machine.jlso")),
		],
		resampling=Holdout(),
		measure=cross_entropy,
		retrain=true,
	)
	# f_ranges = [
	# 	range(f_mdl, :(model.builder.n_hidden_l), lower=1, upper=6, scale=:linear),
	# 	range(f_mdl, :(model.builder.dropout), lower=0, upper=1, scale=:linear),
	# ]
	f_mach1 = machine(
		f_mdl,
		train_x,
		train_y,
	)
	# f_mach2 = machine(
	# 	TunedModel(
	# 		model=f_mdl,
	# 		tuning=Grid(
	# 			resolution=3,
	# 		),
	# 		resampling=Holdout(),
	# 		measure=cross_entropy,
	# 		range=f_ranges,
	# 		acceleration=CPU1(),
	# 		acceleration_resampling=CPU1(),
	# 	),
	# 	train_x,
	# 	train_y,
	# )
end

# ╔═╡ a1f4cd40-654b-4586-8bf0-3729b8415d3b
f_mach1

# ╔═╡ 08910ec3-1495-4e94-bdee-785c6148f1ee
fit!(f_mach1)

# ╔═╡ 016eb179-7add-4797-be11-a6011a7ec57a
report(f_mach1)

# ╔═╡ 0dd8e8a5-afe1-492f-a4f5-fb37ea15564e
f_mach1

# ╔═╡ 48b1b6c9-9beb-49fd-94d4-831712bccdb3
predicted_labels = MLJ.predict(f_mach1)

# ╔═╡ d54771fc-2ece-44b6-bfd7-138ebf60bb56
cross_entropy(predicted_labels, train_y) |> mean

# ╔═╡ 3d6d9bf8-c1f9-4089-bc72-301331ed749a
evaluate!(f_mach1, measure=cross_entropy)

# ╔═╡ 9c847976-1d10-4d8c-83dd-a9cc6600febc
# begin
# 	r = report(f_mach1).plotting
# 	fr = hcat(r.parameter_values, -1 .* r.measurements)

# 	Statistics.cor(fr)
# end

# ╔═╡ 9b60e4c1-c29a-49d6-846c-b2c7a942b9ca
# begin
# 	curve = report(f_mach1).plotting
# 	Plots.scatter(
# 		(x -> x[1]).(eachrow(curve.parameter_values)),
# 		(x -> x[2]).(eachrow(curve.parameter_values)),
# 		curve.measurements,
#      	xlab=curve.parameter_names[1],
# 		ylab=curve.parameter_names[2],
#      	xscale=curve.parameter_scales[1],
# 		yscale=curve.parameter_scales[2],
# 		zlab="CV estimate of error",
# 		# st=:surface,
# 	)
# end

# ╔═╡ a55dd90f-a2d1-48e5-97dd-d9afc91aa38a
size(training_losses), size(losses)

# ╔═╡ 4a9a00f0-6a97-436c-a483-77db50a123de
begin
	Plots.plot(
		losses,
		title="Cross Entropy",
		xlab = "epoch",
		label="out-of-sample",
	)
	Plots.plot!(
		training_losses,
		label="training",
	)
end

# ╔═╡ Cell order:
# ╠═97539bb6-c0b0-11eb-2d83-316749d179f1
# ╠═baab9408-fe06-436d-8158-6ec55e1b49cf
# ╠═5bafc48d-eb03-4625-8d87-17170133d999
# ╠═e593f5a3-9c57-42e7-9ba5-fcf4f0c4bb76
# ╠═3d5dac10-a6a9-4e13-b6c7-85c86d5f61ac
# ╠═b48567da-38c0-4d01-b866-08413433d115
# ╠═85bb6d5a-5117-4b24-beb1-c1b00f0e9462
# ╠═984cbe16-be12-4fdf-8345-65b601b2c4cd
# ╠═4a988da3-1a6c-4bb0-a15c-779e2f2b5cee
# ╠═03dd99ff-edf2-4c66-ac15-a4dff9b2c7c6
# ╠═58ffa080-cf63-43d9-bc40-d9fedb085755
# ╠═d0e5c847-6bd5-484d-9e62-24deb3a8f0a5
# ╠═a1f4cd40-654b-4586-8bf0-3729b8415d3b
# ╠═08910ec3-1495-4e94-bdee-785c6148f1ee
# ╠═016eb179-7add-4797-be11-a6011a7ec57a
# ╠═0dd8e8a5-afe1-492f-a4f5-fb37ea15564e
# ╠═48b1b6c9-9beb-49fd-94d4-831712bccdb3
# ╠═d54771fc-2ece-44b6-bfd7-138ebf60bb56
# ╠═3d6d9bf8-c1f9-4089-bc72-301331ed749a
# ╠═9c847976-1d10-4d8c-83dd-a9cc6600febc
# ╠═9b60e4c1-c29a-49d6-846c-b2c7a942b9ca
# ╠═a55dd90f-a2d1-48e5-97dd-d9afc91aa38a
# ╠═4a9a00f0-6a97-436c-a483-77db50a123de

max_stars_repo_name: prbzrg/mnist-nn

max_issues_repo_name: prbzrg/mnist-nn

max_forks_repo_name: prbzrg/mnist-nn


text:
# ---
# title: 1578. Minimum Deletion Cost to Avoid Repeating Letters
# id: problem1578
# author: Indigo
# date: 2021-01-27
# difficulty: Medium
# categories: Greedy
# link: <https://leetcode.com/problems/minimum-deletion-cost-to-avoid-repeating-letters/description/>
# hidden: true
# ---
#
# Given a string `s` and an array of integers `cost` where `cost[i]` is the cost
# of deleting the `ith` character in `s`.
#
# Return the minimum cost of deletions such that there are no two identical
# letters next to each other.
#
# Notice that you will delete the chosen characters at the same time, in other
# words, after deleting a character, the costs of deleting other characters will
# not change.
#
#
#
# **Example 1:**
#
#
#
#     Input: s = "abaac", cost = [1,2,3,4,5]
#     Output: 3
#     Explanation: Delete the letter "a" with cost 3 to get "abac" (String without two identical letters next to each other).
#
#
# **Example 2:**
#
#
#
#     Input: s = "abc", cost = [1,2,3]
#     Output: 0
#     Explanation: You don't need to delete any character because there are no identical letters next to each other.
#
#
# **Example 3:**
#
#
#
#     Input: s = "aabaa", cost = [1,2,3,4,1]
#     Output: 2
#     Explanation: Delete the first and the last character, getting the string ("aba").
#
#
#
#
# **Constraints:**
#
#   * `s.length == cost.length`
#   * `1 <= s.length, cost.length <= 10^5`
#   * `1 <= cost[i] <= 10^4`
#   * `s` contains only lowercase English letters.
#
#
## @lc code=start
using LeetCode

function min_cost_1578(s::String, cost::Vector{Int})
    if length(s) == 0
        return 0
    end
    s *= "0"
    lst_char = s[1]
    len, res = 1, 0
    for i in 2:length(s)
        if s[i] == lst_char
            len += 1
        else
            res += sum(cost[i-len:i-1]) - maximum(cost[i-len:i-1])
            len, lst_char = 1, s[i]
        end
    end
    res
end
## @lc code=end


text:
using CSV
using DataFrames

using Base.Filesystem: path_separator
using AWS: @service

@service S3
sync(bucket::AbstractString, dest::AbstractString) = _sync(bucket, dest, nothing)
function _sync(bucket, dest, token)
    params = Dict{String, String}()
    token === nothing || push!(params, "continuation-token" => token)
    resp = S3.list_objects_v2(bucket, params)
    for item in resp["Contents"]
        key = item["Key"]
        bytes = S3.get_object(bucket, key)
        path = joinpath(dest, replace(key, "/" => path_separator))
        mkpath(dirname(path))
        write(path, bytes)
    end
    if resp["IsTruncated"] == "true"
        _sync(bucket, dest, resp["NextContinuationToken"])
    end
end

sync("openneuro.org/ds002680", "download")

max_stars_repo_name: TheCedarPrince/NeuriViz

max_issues_repo_name: TheCedarPrince/NeuriViz

max_forks_repo_name: TheCedarPrince/NeuriViz