text:
struct ManualString <: AbstractString
    ptr::Manual{UInt8}
    len::Int64 # in bytes
end

# creation

function Base.unsafe_copy!(ps::ManualString, string::Union{ManualString, String})
    @assert ps.len >= string.len
    unsafe_copy!(convert(Ptr{UInt8}, ps.ptr.ptr), pointer(string), string.len)
end

function Base.String(ps::ManualString)
    unsafe_string(convert(Ptr{UInt8}, ps.ptr.ptr), ps.len)
end

function unsafe_manual_string(ptr::Ptr{UInt8}, len::Int64)
    ManualString(Manual{UInt8}(convert(Ptr{Void}, ptr)), len)
end

# fields

function Base.pointer(ps::ManualString)
    convert(Ptr{UInt8}, ps.ptr.ptr)
end

function Base.pointer(ps::ManualString, i::Integer)
    pointer(ps) + i - 1
end

# string interface - this is largely copied from Base and will almost certainly break when we move to 0.7

Base.sizeof(s::ManualString) = s.len

@inline function Base.codeunit(s::ManualString, i::Integer)
    @boundscheck if (i < 1) | (i > s.len)
        throw(BoundsError(s,i))
    end
    unsafe_load(pointer(s),i)
end

Base.write(io::IO, s::ManualString) = unsafe_write(io, pointer(s), reinterpret(UInt, s.len))

function Base.cmp(a::ManualString, b::ManualString)
    c = ccall(:memcmp, Int32, (Ptr{UInt8}, Ptr{UInt8}, UInt),
              a, b, min(a.len,b.len))
    return c < 0 ? -1 : c > 0 ? +1 : cmp(a.len,b.len)
end

function Base.:(==)(a::ManualString, b::ManualString)
    a.len == b.len && 0 == ccall(:memcmp, Int32, (Ptr{UInt8}, Ptr{UInt8}, UInt), a, b, a.len)
end

function Base.prevind(s::ManualString, i::Integer)
    j = Int(i)
    e = s.len
    if j > e
        return endof(s)
    end
    j -= 1
    @inbounds while j > 0 && Base.is_valid_continuation(codeunit(s,j))
        j -= 1
    end
    j
end

function Base.nextind(s::ManualString, i::Integer)
    j = Int(i)
    if j < 1
        return 1
    end
    e = s.len
    j += 1
    @inbounds while j <= e && Base.is_valid_continuation(codeunit(s,j))
        j += 1
    end
    j
end

Base.byte_string_classify(s::ManualString) =
    ccall(:u8_isvalid, Int32, (Ptr{UInt8}, Int), s, s.len)

Base.isvalid(::Type{ManualString}, s::ManualString) = byte_string_classify(s) != 0
Base.isvalid(s::ManualString) = isvalid(ManualString, s)

function Base.endof(s::ManualString)
    p = pointer(s)
    i = s.len
    while i > 0 && Base.is_valid_continuation(unsafe_load(p,i))
        i -= 1
    end
    i
end

function Base.length(s::ManualString)
    p = pointer(s)
    cnum = 0
    for i = 1:s.len
        cnum += !Base.is_valid_continuation(unsafe_load(p,i))
    end
    cnum
end

Base.done(s::ManualString, state) = state > s.len

@inline function Base.next(s::ManualString, i::Int)
    @boundscheck if (i < 1) | (i > s.len)
        throw(BoundsError(s,i))
    end
    p = pointer(s)
    b = unsafe_load(p, i)
    if b < 0x80
        return Char(b), i + 1
    end
    return Base.slow_utf8_next(p, b, i, s.len)
end

function Base.reverseind(s::ManualString, i::Integer)
    j = s.len + 1 - i
    p = pointer(s)
    while Base.is_valid_continuation(unsafe_load(p,j))
        j -= 1
    end
    return j
end

Base.isvalid(s::ManualString, i::Integer) =
    (1 <= i <= s.len) && !Base.is_valid_continuation(unsafe_load(pointer(s),i))

function Base.getindex(s::ManualString, r::UnitRange{Int})
    isempty(r) && return ""
    i, j = first(r), last(r)
    l = s.len
    if i < 1 || i > l
        throw(BoundsError(s, i))
    end
    @inbounds si = Base.codeunit(s, i)
    if Base.is_valid_continuation(si)
        throw(UnicodeError(UTF_ERR_INVALID_INDEX, i, si))
    end
    if j > l
        throw(BoundsError())
    end
    j = nextind(s,j)-1
    unsafe_manual_string(pointer(s,i), j-i+1)
end

function Base.search(s::ManualString, c::Char, i::Integer = 1)
    if i < 1 || i > sizeof(s)
        i == sizeof(s) + 1 && return 0
        throw(BoundsError(s, i))
    end
    if Base.is_valid_continuation(Base.codeunit(s,i))
        throw(UnicodeError(UTF_ERR_INVALID_INDEX, i, Base.codeunit(s,i)))
    end
    c < Char(0x80) && return search(s, c%UInt8, i)
    while true
        i = search(s, Base.first_utf8_byte(c), i)
        (i==0 || s[i] == c) && return i
        i = next(s,i)[2]
    end
end

function Base.search(a::ManualString, b::Union{Int8,UInt8}, i::Integer = 1)
    if i < 1
        throw(BoundsError(a, i))
    end
    n = sizeof(a)
    if i > n
        return i == n+1 ? 0 : throw(BoundsError(a, i))
    end
    p = pointer(a)
    q = ccall(:memchr, Ptr{UInt8}, (Ptr{UInt8}, Int32, Csize_t), p+i-1, b, n-i+1)
    q == C_NULL ? 0 : Int(q-p+1)
end

function Base.rsearch(s::ManualString, c::Char, i::Integer = s.len)
    c < Char(0x80) && return rsearch(s, c%UInt8, i)
    b = Base.first_utf8_byte(c)
    while true
        i = rsearch(s, b, i)
        (i==0 || s[i] == c) && return i
        i = prevind(s,i)
    end
end

function Base.rsearch(a::ManualString, b::Union{Int8,UInt8}, i::Integer = s.len)
    if i < 1
        return i == 0 ? 0 : throw(BoundsError(a, i))
    end
    n = sizeof(a)
    if i > n
        return i == n+1 ? 0 : throw(BoundsError(a, i))
    end
    p = pointer(a)
    q = ccall(:memrchr, Ptr{UInt8}, (Ptr{UInt8}, Int32, Csize_t), p, b, i)
    q == C_NULL ? 0 : Int(q-p+1)
end

function Base.string(a::ManualString...)
    if length(a) == 1
        return String(a[1]::ManualString)
    end
    n = 0
    for str in a
        n += str.len
    end
    out = Base._string_n(n)
    offs = 1
    for str in a
        unsafe_copy!(pointer(out,offs), pointer(str), str.len)
        offs += str.len
    end
    return out
end


text:
function vis(KS::KineticSolver,psi::Array{Float64,3})

	rho = computeMoments(KS,psi)
	#println("Final mass = $(sum(rho[:]))." )
	nicewrite(string(KS.outputfolderprefix,"data/rhofinal.txt"),rho)

	if KS.ictype == "linesource"
		fig, ax = subplots(figsize=(10.5, 8), dpi=100)
		minVal = 0.0; maxVal = 0.5;
	#	title("$(KS.Q.quadtype), nq=$(KS.nquadpoints), cs=$(KS.convolutionmagnitude), rot=$(KS.rotationmagnitude)", fontsize=20)
		pcolormesh(rho, cmap="RdBu_r",vmin=minVal, vmax=maxVal)
		colorbar()
		PyPlot.savefig(string(KS.outputfolderprefix,"rho.png"))

		rhoana = readdlm("../vis/exactLineSource.txt",',')
		fig, ax = subplots(figsize=(10.5, 8), dpi=100)
		plot(rhoana[:,1],rhoana[:,2])
		rhong  = rho[3:end-2,3:end-2] # remove ghosts
		cuthori = rhong[Int64(KS.ny/2),:]
		x = range(-1.5,stop=1.5,length=KS.nx-4)
		plot(x,cuthori)
		cutverti = rhong[:,Int64(KS.nx/2)]
		plot(x,cutverti)
		cutdiag = [rhong[i,i] for i=1:KS.nx-4]
		plot(sqrt(2)*x,cutdiag)
		cutdiag = [rhong[KS.ny-4+1-i,i] for i=1:KS.nx-4]
		plot(sqrt(2)*x,cutdiag)
		title("$(KS.Q.quadtype), nq=$(KS.nquadpoints), convolution=$(KS.convolutionmagnitude) ", fontsize=30)
		PyPlot.savefig(string(KS.outputfolderprefix,"rhocut.png"))

	elseif KS.ictype == "checkerboard"
		fig, ax = subplots(figsize=(10.5, 8), dpi=100)
		minVal = -7.0; maxVal = 0.0;
		title("$(KS.Q.quadtype), nq=$(KS.nquadpoints), cs=$(KS.convolutionmagnitude)", fontsize=20)
		pcolormesh(log10.(abs.(rho)), cmap="RdBu_r",vmin=minVal, vmax=maxVal)
		colorbar()
		PyPlot.savefig(string(KS.outputfolderprefix,"rho_log10.png"))

	elseif KS.ictype == "beam"

		fig, ax = subplots(figsize=(10.5, 8), dpi=100)
		minVal = -7.0; maxVal = 0.0;
		title("$(KS.Q.quadtype), nq=$(KS.nquadpoints), cs=$(KS.convolutionmagnitude)", fontsize=20)
		pcolormesh(log10.(abs.(rho[3:end-2,3:end-2])), cmap="RdBu_r")#,vmin=minVal, vmax=maxVal)
		colorbar()
		PyPlot.savefig(string(KS.outputfolderprefix,"rho_log10.png"))

		fig, ax = subplots(figsize=(10.5, 8), dpi=100)
		minVal = -7.0; maxVal = 0.0;
		title("$(KS.Q.quadtype), nq=$(KS.nquadpoints), cs=$(KS.convolutionmagnitude)", fontsize=20)
		pcolormesh(((rho[3:end-2,3:end-2])), cmap="RdBu_r")#,vmin=minVal, vmax=maxVal)
		colorbar()
		PyPlot.savefig(string(KS.outputfolderprefix,"rho.png"))

	end

	if KS.convolutionmagnitude>0
		nicewrite(string(KS.outputfolderprefix,"data/convolutionmat.txt"),Matrix(KS.ScatteringKernelConvolution))
		fig, ax = subplots(figsize=(10.5, 8), dpi=100)
		title("Convolutionmatrix", fontsize=20)
		pcolormesh((Matrix(KS.ScatteringKernelConvolution)), cmap="RdBu_r")
		colorbar()
		PyPlot.savefig(string(KS.outputfolderprefix,"convmat.png"))
	end
end

function visMoments(KS::KineticSolver,phi,name::String)

	rho = phi[1,:,:]
	#println("Final mass = $(sum(rho[:]))." )
	#nicewrite(string(KS.outputfolderprefix,"data/rhofinalMoments",name,".txt"),rho)

	if KS.ictype == "linesource"
		fig, ax = subplots(figsize=(10.5, 8), dpi=100)
		minVal = 0.0; maxVal = 0.5;
	#	title("$(KS.Q.quadtype), nq=$(KS.nquadpoints), cs=$(KS.convolutionmagnitude), rot=$(KS.rotationmagnitude)", fontsize=20)
		pcolormesh(rho, cmap="RdBu_r",vmin=minVal, vmax=maxVal)
		colorbar()
		PyPlot.savefig(string(KS.outputfolderprefix,"rho",name,".png"))

		rhoana = readdlm("../vis/exactLineSource.txt",',')
		fig, ax = subplots(figsize=(10.5, 8), dpi=100)
		plot(rhoana[:,1],rhoana[:,2])
		rhong  = rho[3:end-2,3:end-2] # remove ghosts
		cuthori = rhong[Int64(KS.ny/2),:]
		x = range(-1.5,stop=1.5,length=KS.nx-4)
		plot(x,cuthori)
		cutverti = rhong[:,Int64(KS.nx/2)]
		plot(x,cutverti)
		cutdiag = [rhong[i,i] for i=1:KS.nx-4]
		plot(sqrt(2)*x,cutdiag)
		cutdiag = [rhong[KS.ny-4+1-i,i] for i=1:KS.nx-4]
		plot(sqrt(2)*x,cutdiag)
		title("$(KS.Q.quadtype), nq=$(KS.nquadpoints), convolution=$(KS.convolutionmagnitude) ", fontsize=30)
		PyPlot.savefig(string(KS.outputfolderprefix,"rhocut",name,".png"))

	elseif KS.ictype == "checkerboard"
		fig, ax = subplots(figsize=(10.5, 8), dpi=100)
		minVal = -7.0; maxVal = 0.0;
		title("$(KS.Q.quadtype), nq=$(KS.nquadpoints), cs=$(KS.convolutionmagnitude)", fontsize=20)
		pcolormesh(log10.(abs.(rho)), cmap="RdBu_r",vmin=minVal, vmax=maxVal)
		colorbar()
		PyPlot.savefig(string(KS.outputfolderprefix,"rho_log10.png"))

	elseif KS.ictype == "beam"

		fig, ax = subplots(figsize=(10.5, 8), dpi=100)
		minVal = -7.0; maxVal = 0.0;
		title("$(KS.Q.quadtype), nq=$(KS.nquadpoints), cs=$(KS.convolutionmagnitude)", fontsize=20)
		pcolormesh(log10.(abs.(rho[3:end-2,3:end-2])), cmap="RdBu_r")#,vmin=minVal, vmax=maxVal)
		colorbar()
		PyPlot.savefig(string(KS.outputfolderprefix,"rho_log10.png"))

		fig, ax = subplots(figsize=(10.5, 8), dpi=100)
		minVal = -7.0; maxVal = 0.0;
		title("$(KS.Q.quadtype), nq=$(KS.nquadpoints), cs=$(KS.convolutionmagnitude)", fontsize=20)
		pcolormesh(((rho[3:end-2,3:end-2])), cmap="RdBu_r")#,vmin=minVal, vmax=maxVal)
		colorbar()
		PyPlot.savefig(string(KS.outputfolderprefix,"rho.png"))

	end

	if KS.convolutionmagnitude>0
		nicewrite(string(KS.outputfolderprefix,"data/convolutionmat.txt"),Matrix(KS.ScatteringKernelConvolution))
		fig, ax = subplots(figsize=(10.5, 8), dpi=100)
		title("Convolutionmatrix", fontsize=20)
		pcolormesh((Matrix(KS.ScatteringKernelConvolution)), cmap="RdBu_r")
		colorbar()
		PyPlot.savefig(string(KS.outputfolderprefix,"convmat.png"))
	end
end

function summarizeAsPlot(KS::KineticSolver,psi::Array{Float64,3})
	figure()
	subplot(221)
	pcolormesh(KS.SigmaS)
	subplot(222)
	pcolormesh(KS.SigmaT-KS.SigmaS)
	subplot(223)
	pcolormesh(computeMoments(KS,KS.Source))
	subplot(224)
	pcolormesh(computeMoments(KS,KS.psi0))
	PyPlot.savefig(string(KS.outputfolderprefix,"testcase.png"))

end

ext: jl

max_stars_repo_name: camminady/SN

max_issues_repo_name: camminady/SN

max_forks_repo_name: camminady/SN


text:
struct ChainDataFrame{NT<:NamedTuple}
    name::String
    nt::NT
    n_rows::Int
    n_cols::Int
    digits::Int

    function ChainDataFrame(name::String, nt::NamedTuple; digits::Int=4)
        lengths = length(first(nt))
        all(x -> length(x) == lengths, nt) || error("Lengths must be equal.")

        return new{typeof(nt)}(name, nt, lengths, length(nt), digits)
    end
end

ChainDataFrame(nt::NamedTuple) = ChainDataFrame("", nt)

Base.size(c::ChainDataFrame) = (c.n_rows, c.n_cols)
Base.names(c::ChainDataFrame) = collect(keys(c.nt))

function Base.show(io::IO, c::ChainDataFrame)
    println(io, c.name)

    # Preallocations
    n_cols = c.n_cols
    n_rows = c.n_rows
    column_names = collect(keys(c.nt))
    digits = c.digits
    f_format_str = "%.$(digits)f"
    i_format_str = "%n"

    # Create printed array.
    arr = Array{String, 2}(undef, n_rows + 1, n_cols)

    # Add headers
    arr[1, :] = string.(column_names)

    # Add values to array, accumulate string lengths.
    lengths = length.(arr[1,:])
    for i in 1:n_cols
        k = column_names[i]
        values = c.nt[k]
        etype = eltype(values)

        if etype <: Real
            arr[2:end, i] = sprintf1.(f_format_str, values)
        elseif etype <: Integer
            arr[2:end, i] = sprintf1.(i_format_str, values)
        else
            arr[2:end, i] = string.(values)
        end
        lengths[i] = max(maximum(length.(arr[2:end, i])), lengths[i])
    end

    # Do it array style.
    bufs = [IOBuffer() for _ in 1:(n_rows+2)]

    for i in 1:n_cols
        for j in 1:(n_rows+2)
            # Print the headers
            if j == 1
                print(bufs[j],  "  ", lpad(arr[j,i], lengths[i]))
            end

            # Print sep row.
            if j == 2
                print(bufs[j], "  ", repeat("─", lengths[i]))
            end

            # Print values.
            if j > 2
                bufs[j]
                print(bufs[j], "  ", lpad(arr[j-1,i], lengths[i]))
            end
        end
    end

    for j in 1:length(bufs)
        s = String(take!(bufs[j]))
        println(io, s)
    end
end

Base.isequal(cs1::Vector{ChainDataFrame}, cs2::Vector{ChainDataFrame}) = isequal.(cs1, cs2)
Base.isequal(c1::ChainDataFrame, c2::ChainDataFrame) = isequal(c1, c2)

function Base.show(io::IO, cs::Vector{C}) where C<:ChainDataFrame
    println(io, summary(cs))
    for i in cs
        println(io)
        show(io, i)
    end
end

# Allows overriding of `display`
function Base.show(io::IO, ::MIME"text/plain", cs::Vector{ChainDataFrame})
    show(io, cs)
end

# Index functions
function Base.getindex(c::ChainDataFrame, s::Union{Colon, Integer, UnitRange}, g::Union{Colon, Integer, UnitRange})
    convert(Array, getindex(c, c.nt[:parameters][s], collect(keys(c.nt))[g]))
end

function Base.getindex(c::ChainDataFrame, s::Union{Symbol, Vector{Symbol}}, g::Colon)
    getindex(c, c.nt[:parameters], collect(keys(c.nt)))
end

function Base.getindex(c::ChainDataFrame, s::Union{Symbol, Vector{Symbol}})
    getindex(c, s, collect(keys(c.nt)))
end

function Base.getindex(c::ChainDataFrame, s::Union{Colon, Integer, UnitRange}, ks::Union{Symbol, Vector{Symbol}})
    getindex(c, c.nt[:parameters][s], ks)
end

function Base.getindex(
    c::ChainDataFrame,
    s::Union{String, Vector{String}, Symbol, Vector{Symbol}},
    ks::Union{Symbol, Vector{Symbol}}
)
    s = s isa AbstractArray ? s : [s]
    ks = ks isa AbstractArray ? ks : [ks]
    ind = indexin(Symbol.(s), Symbol.(c.nt[:parameters]))

    not_found = map(x -> x === nothing, ind)

    any(not_found) && error("Cannot find parameters $(s[not_found]) in chain")

    # If there are multiple columns, return a new CDF.
    if length(ks) > 1
        if !(:parameters in ks)
            ks = vcat(:parameters, ks)
        end
        nt = NamedTuple{tuple(ks...)}(tuple([c.nt[k][ind] for k in ks]...))
        return ChainDataFrame(c.name, nt, digits=c.digits)
    else
        # Otherwise, return a vector if there's multiple parameters
        # or just a scalar if there's one parameter.
        if length(s) == 1
            return c.nt[ks[1]][ind][1]
        else
            return c.nt[ks[1]][ind]
        end
    end
end

function Base.lastindex(c::ChainDataFrame, i::Integer)
    if i == 1
        return c.n_rows
    elseif i ==2
        return c.n_cols
    else
        error("No such dimension")
    end
end

function Base.convert(::Type{Array}, c::C) where C<:ChainDataFrame
    T = promote_eltype_namedtuple_tail(c.nt)
    arr = Array{T, 2}(undef, c.n_rows, c.n_cols - 1)

    for (i, k) in enumerate(Iterators.drop(keys(c.nt), 1))
        arr[:, i] = c.nt[k]
    end

    return arr
end

Base.convert(::Type{Array{ChainDataFrame,1}}, cs::Array{ChainDataFrame,1}) = cs
function Base.convert(::Type{Array}, cs::Array{C,1}) where C<:ChainDataFrame
    return mapreduce((x, y) -> cat(x, y; dims = Val(3)), cs) do c
        reshape(convert(Array, c), Val(3))
    end
end

"""

# Summarize a Chains object formatted as a ChainsDataFrame

Summarize method for a Chains object.

### Method
```julia
  summarize(
    chn::Chains,
    funs...;
    sections::Vector{Symbol}=[:parameters],
    func_names=[],
    etype=:bm
  )
```

### Required arguments
```julia
* `chn` : Chains object to convert to a ChainsDataFrame-formatted summary
```

### Optional arguments
```julia
* `funs...` : zero or more vector functions, e.g. mean, std, etc.
* `sections = [:parameters]` : Sections from the Chains object to be included
* `etype = :bm` : Default for df_mcse
```

### Examples
```julia
* `summarize(chns)` : Complete chain summary
* `summarize(chns[[:parm1, :parm2]])` : Chain summary of selected parameters
* `summarize(chns, sections=[:parameters])`  : Chain summary of :parameters section
* `summarize(chns, sections=[:parameters, :internals])` : Chain summary for multiple sections
```

"""
function summarize(chains::Chains, funs...;
        sections::Union{Symbol, Vector{Symbol}}=Symbol[:parameters],
        func_names::AbstractVector{Symbol} = Symbol[],
        append_chains::Bool=true,
        showall::Bool=false,
        name::String="",
        additional_df=nothing,
        digits::Int=4,
        sorted::Bool=false)
    # Check that we actually have :parameters.
    showall = showall || !in(:parameters, keys(chains.name_map))

    # If we weren't given any functions, fall back on summary stats.
    if length(funs) == 0
        return summarystats(chains,
            sections=sections,
            showall=showall)
    end

    # Generate a chain to work on.
    chn = Chains(chains, sections, sorted=sorted)

    # Obtain names of parameters.
    names_of_params = names(chn)

    # If no function names were given, make a new list.
    fnames = isempty(func_names) ? collect(nameof.(funs)) : func_names

    # Obtain the additional named tuple.
    additional_nt = additional_df === nothing ? NamedTuple() : additional_df.nt

    if append_chains
        # Evaluate the functions.
        data = to_matrix(chn)
        fvals = [[f(data[:, i]) for i in axes(data, 2)] for f in funs]

        # Build the ChainDataFrame.
        nt = merge((; parameters = names_of_params, zip(fnames, fvals)...), additional_nt)
        df = ChainDataFrame(name, nt; digits=digits)

        return df
    else
        # Evaluate the functions.
        data = to_vector_of_matrices(chn)
        vector_of_fvals = [[[f(x[:, i]) for i in axes(x, 2)] for f in funs] for x in data]

        # Build the ChainDataFrames.
        vector_of_nt = [
            merge((; parameters = names_of_params, zip(fnames, fvals)...), additional_nt)
            for fvals in vector_of_fvals
        ]
        vector_of_df = [
            ChainDataFrame(name * " (Chain $i)", nt; digits=digits)
            for (i, nt) in enumerate(vector_of_nt)
        ]

        return vector_of_df
    end
end


text:
@params struct PartialPolynomialFit{T}
    data::Tuple{Matrix{T}, Vector{T}}
    order::Ref{Int}
    coeffs::AbstractVector{T}
    oldest::Ref{Int}
end
function Base.convert(::Type{PartialPolynomialFit}, n::Int)
    PartialPolynomialFit{Float64}(n)
end
function Base.convert(::Type{PartialPolynomialFit{T}}, n::Int) where T
    PartialPolynomialFit{T}((zeros(T, n+1, n+1), zeros(T, n+1)), Ref(-1),
        zeros(T, n+1), Ref(1))
end
Base.eltype(::PartialPolynomialFit{T}) where T = T

@params struct Derivative
    f
end
# nth order derivatives of polynomial function f using f'''(x) notation
@inline Base.transpose(p::Union{PartialPolynomialFit, Derivative}) = Derivative(p)
@inline (poly::PartialPolynomialFit)(x) = poly(x, 0)
@inline (poly::PartialPolynomialFit)(x, n) = nthderivativeat(poly, n, x)
@inline (der::Derivative)(x) = der.f(x, 1)
@inline (der::Derivative)(x, n) = der.f(x, n+1)
function nthderivativeat(poly::PartialPolynomialFit{T}, n, x) where T
    o = zero(T)
    for j in (n+1):order(poly)+1
        p = n == 0 ? one(T) : prod(i->(j-i), 1:n)
        o += p*poly.coeffs[j]*x^(j-1-n)
    end
    return o
end
Base.eltype(d::Derivative) = eltype(d.f)

maxorder(p::PartialPolynomialFit) = length(p.coeffs)-1
order(p::PartialPolynomialFit) = p.order[]
increaseorder!(p::PartialPolynomialFit) = p.order[]+= 1
reset!(p::PartialPolynomialFit) = p.order[] = -1
oldest(p::PartialPolynomialFit) = p.oldest[]
newest(p::PartialPolynomialFit) = p.oldest[] == 1 ? order(p) + 1 : p.oldest[] - 1

previous(p::PartialPolynomialFit, i) = i == 1 ? order(p) + 1 : i - 1
next(p::PartialPolynomialFit, i) = i == order(p) + 1 ? 1 : i + 1

function newvalues!(p::PartialPolynomialFit, X, V)
    for i in 1:length(X)
        newvalue!(p, X[i], V[i])
    end
    oldest(p)
end
function newvalue!(p::PartialPolynomialFit, x, v)
    n = incrementoldest!(p)
    if order(p) < maxorder(p)
        increaseorder!(p)
    end
    newvalue!(p, x, v, n)
end
function newvalue!(p::PartialPolynomialFit, x, v, i)
    for j in 1:maxorder(p)+1
        p.data[1][i,j] = x^(j-1)
    end
    p.data[2][i] = v
    oldest(p)
end

function newderivatives!(p::PartialPolynomialFit, X, V)
    for i in 1:length(X)
        newderivative!(p, X[i], V[i])
    end
    oldest(p)
end
function newderivative!(p::PartialPolynomialFit, x, v)
    n = incrementoldest!(p)
    if order(p) < maxorder(p)
        increaseorder!(p)
    end
    newderivative!(p, x, v, n)
end
function newderivative!(p::PartialPolynomialFit, x, v, i)
    for j in 1:maxorder(p)+1
        p.data[1][i,j] = (j-1)x^(j-2)
    end
    p.data[2][i] = v
    oldest(p)
end

function incrementoldest!(p::PartialPolynomialFit)
    o = p.oldest[]
    if o == maxorder(p)+1
        p.oldest[] = 1
    else
        p.oldest[] = o + 1
    end
    o
end

function solve!(p::PartialPolynomialFit)
    n = order(p)+1
    p.coeffs[1:n] .= p.data[1][1:n, 1:n] \ p.data[2][1:n]
    @view(p.coeffs[1:n])
end

function _filter(root_ranges)
    unique_root_ranges = filter!(x->(x.status == :unique), root_ranges)
    return [(r.interval.lo + r.interval.lo)/2 for r in unique_root_ranges]
end

function _filter(optimiser_ranges::Tuple)
    return [(r.lo + r.lo)/2 for r in optimiser_ranges[2]]
end

function roots(poly::Union{TP, Derivative{TP}}, min=T(0), max=T(15)) where {T, TP<:PartialPolynomialFit{T}}
    root_ranges = IntervalRootFinding.roots(x->poly(x), min..max)
    return _filter(root_ranges)
end

@params struct PolynomialQuotient
    num
    denom
end
Base.:/(p1::Union{PartialPolynomialFit, Derivative}, p2::Union{PartialPolynomialFit, Derivative}) = PolynomialQuotient(p1, p2)
Base.eltype(p::PolynomialQuotient) = promote_type(eltype(p.num), eltype(p.denom))

@params struct PolynomialEquality
    LHS
    RHS
end
Base.:(==)(p::Union{PolynomialQuotient, PartialPolynomialFit, Derivative}, n) = PolynomialEquality(p, n)
Base.eltype(p::PolynomialEquality) = promote_type(eltype(p.LHS), typeof(p.RHS))

function roots(poly::PolynomialEquality{TP}, min=eltype(poly)(0), max=eltype(poly)(15)) where {TP<:Union{PartialPolynomialFit, Derivative{<:PartialPolynomialFit}}}
    root_ranges = IntervalRootFinding.roots(x->(poly.LHS(x) - poly.RHS), min..max)
    return _filter(root_ranges)
end

function roots(poly::PolynomialEquality{<:PolynomialQuotient}, min=eltype(poly)(0), max=eltype(poly)(15))
    root_ranges = IntervalRootFinding.roots(x->(poly.LHS.num(x) - poly.RHS * poly.LHS.denom(x)), min..max)
    return _filter(root_ranges)
end

# TO-DO: derivative of a quotient using the quotient rule


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    Base.precompile(Tuple{Core.kwftype(typeof(Tullio.Type)),Any,Type{Tullio.DotDict}})
    Base.precompile(Tuple{Type{Base.Broadcast.Broadcasted{Base.Broadcast.Style{Tuple}, Axes, F, Args} where Args<:Tuple where F where Axes},typeof(Tullio.dollarstrip),Tuple{Tuple{Expr, Expr}}})
    Base.precompile(Tuple{typeof(Base.Broadcast.instantiate),Base.Broadcast.Broadcasted{Base.Broadcast.Style{Tuple}, Nothing, typeof(Tullio.dollarstrip), Tuple{Tuple{Expr, Expr}}}})
    Base.precompile(Tuple{typeof(Base.setindex_widen_up_to),Vector{typeof(Tullio.subranges)},Expr,Int64})
    Base.precompile(Tuple{typeof(Tullio._trymatch),Any,Any})
    Base.precompile(Tuple{typeof(Tullio._trymatch),Expr,Union{Val{:ref}, Val{:curly}}})
    Base.precompile(Tuple{typeof(Tullio._tullio),Any,Vararg{Any, N} where N})
    Base.precompile(Tuple{typeof(Tullio._tullio),Any})
    Base.precompile(Tuple{typeof(Tullio.action_functions),Any})
    Base.precompile(Tuple{typeof(Tullio.arrayplusepsilon),Symbol,Any,Any})
    Base.precompile(Tuple{typeof(Tullio.callcost),Symbol,Tullio.DotDict})
    Base.precompile(Tuple{typeof(Tullio.csewalk),Any,Any,Any,Any})
    Base.precompile(Tuple{typeof(Tullio.csewalk),Expr,Dict,Vector{T} where T,Set})
    Base.precompile(Tuple{typeof(Tullio.illegal),Any,Any})
    Base.precompile(Tuple{typeof(Tullio.illegal),Expr,Any})
    Base.precompile(Tuple{typeof(Tullio.index_ranges),Any})
    Base.precompile(Tuple{typeof(Tullio.make_many_actors),Any,Any,Any,Vector{T} where T,Any,Vector{T} where T,Any,Any,Any,Any})
    Base.precompile(Tuple{typeof(Tullio.make_many_actors),Any,Any,Any,Vector{T} where T,Any,Vector{T} where T,Any,Any,Any})
    Base.precompile(Tuple{typeof(Tullio.output_array),Any})
    Base.precompile(Tuple{typeof(Tullio.padmodclamp_ind),Any,Any,Any})
    Base.precompile(Tuple{typeof(Tullio.padmodclamp_ind),Expr,Any,Any})
    Base.precompile(Tuple{typeof(Tullio.padmodclamp_pair),Any,Any,Any})
    Base.precompile(Tuple{typeof(Tullio.padmodclamp_replace),Expr,Any,Any})
    Base.precompile(Tuple{typeof(Tullio.parse_input),Any,Any})
    Base.precompile(Tuple{typeof(Tullio.parse_ranges),Any,Any})
    Base.precompile(Tuple{typeof(Tullio.postwalk),Any,Any})
    Base.precompile(Tuple{typeof(Tullio.range_expr_walk),Any,Expr,Any})
    Base.precompile(Tuple{typeof(Tullio.resolveintersect),Any,Any,Any})
    Base.precompile(Tuple{typeof(Tullio.resolvestrict),Any,Any,Any})
    Base.precompile(Tuple{typeof(Tullio.saveconstraints),Any,Any,Any,Any})
    Base.precompile(Tuple{typeof(Tullio.simplitimes),Number,Any})
    Base.precompile(Tuple{typeof(Tullio.tidyleftraw),Any,Any})
    Base.precompile(Tuple{typeof(Tullio.walk),Any,Any,Any})
    Base.precompile(Tuple{typeof(Tullio.walk),Expr,Any,Any})
    Base.precompile(Tuple{typeof(copy),Base.Broadcast.Broadcasted{Base.Broadcast.Style{Tuple}, Nothing, typeof(Tullio.dollarstrip), Tuple{Tuple{Expr, Expr}}}})
    isdefined(Tullio, Symbol("#113#114")) && Base.precompile(Tuple{getfield(Tullio, Symbol("#113#114")),Any})
    isdefined(Tullio, Symbol("#139#141")) && Base.precompile(Tuple{getfield(Tullio, Symbol("#139#141")),Any})
    isdefined(Tullio, Symbol("#140#142")) && Base.precompile(Tuple{getfield(Tullio, Symbol("#140#142")),Any})
    isdefined(Tullio, Symbol("#23#24")) && Base.precompile(Tuple{getfield(Tullio, Symbol("#23#24")),Any})
    isdefined(Tullio, Symbol("#25#26")) && Base.precompile(Tuple{getfield(Tullio, Symbol("#25#26")),Any})
    isdefined(Tullio, Symbol("#31#34")) && Base.precompile(Tuple{getfield(Tullio, Symbol("#31#34")),Any})
    isdefined(Tullio, Symbol("#37#38")) && Base.precompile(Tuple{getfield(Tullio, Symbol("#37#38")),Any})
    isdefined(Tullio, Symbol("#39#40")) && Base.precompile(Tuple{getfield(Tullio, Symbol("#39#40")),Any})
    isdefined(Tullio, Symbol("#59#62")) && Base.precompile(Tuple{getfield(Tullio, Symbol("#59#62")),Any})
    isdefined(Tullio, Symbol("#61#64")) && Base.precompile(Tuple{getfield(Tullio, Symbol("#61#64")),Any,Any})
    isdefined(Tullio, Symbol("#65#67")) && Base.precompile(Tuple{getfield(Tullio, Symbol("#65#67")),Any})
    isdefined(Tullio, Symbol("#66#68")) && Base.precompile(Tuple{getfield(Tullio, Symbol("#66#68")),Any})
    isdefined(Tullio, Symbol("#75#77")) && Base.precompile(Tuple{getfield(Tullio, Symbol("#75#77")),Any})
    isdefined(Tullio, Symbol("#76#78")) && Base.precompile(Tuple{getfield(Tullio, Symbol("#76#78")),Any})
    isdefined(Tullio, Symbol("#79#84")) && Base.precompile(Tuple{getfield(Tullio, Symbol("#79#84")),Any})
    isdefined(Tullio, Symbol("#80#85")) && Base.precompile(Tuple{getfield(Tullio, Symbol("#80#85")),Any})
    isdefined(Tullio, Symbol("#81#86")) && Base.precompile(Tuple{getfield(Tullio, Symbol("#81#86")),Any})
    isdefined(Tullio, Symbol("#82#87")) && Base.precompile(Tuple{getfield(Tullio, Symbol("#82#87")),Any,Any})
end


text:
####################################
#
# loss tr(Σx⋅θ⋅Σy⋅θ)/2 - tr(Θ⋅E_ab)    ---> computes inverse of (Σy ⊗ Σx)
#
####################################
struct CDInverseKroneckerLoss{T<:AbstractFloat, S} <: CoordinateDifferentiableFunction
  Σx::Symmetric{T, S}
  Σy::Symmetric{T, S}
  A::Matrix{T}    # stores Σx⋅Θ⋅Σy
  a::Int
  b::Int
  p::Int64
end

function CDInverseKroneckerLoss(Σx::Symmetric{T,S}, Σy::Symmetric{T,S}, a::Int, b::Int) where {T<:AbstractFloat} where S
  (issymmetric(Σx) && issymmetric(Σy)) || throw(DimensionMismatch())
  (p = size(Σx, 1)) == size(Σy, 1) || throw(DimensionMismatch())
  CDInverseKroneckerLoss{T,S}(Σx, Σy, zeros(T, p, p), a, b, p)
end

CoordinateDescent.numCoordinates(f::CDInverseKroneckerLoss) = f.p*f.p
function CoordinateDescent.initialize!(f::CDInverseKroneckerLoss, x::SparseIterate)
# compute residuals for the loss

  Σx = f.Σx
  Σy = f.Σy
  A = f.A
  p = f.p

  for ac=1:p, ar=1:p
      @inbounds A[ar,ac] = A_mul_X_mul_B_rc(Σx, x, Σy, ar, ac)
  end

  nothing
end

function CoordinateDescent.gradient(
  f::CDInverseKroneckerLoss{T},
  x::SparseIterate{T},
  j::Int64) where {T <: AbstractFloat}

  Σx = f.Σx
  Σy = f.Σy
  A = f.A
  # ri, ci = ind2sub(Σx, j)
  ri, ci = Tuple(CartesianIndices(Σx)[j])
  @inbounds v = A[ri,ci]
  return ri == f.a && ci == f.b ? v - 1. : v
end

function CoordinateDescent.descendCoordinate!(
  f::CDInverseKroneckerLoss{T},
  g::ProxL1{T},
  x::SparseIterate{T},
  j::Int64) where {T <: AbstractFloat}

  Σx = f.Σx
  Σy = f.Σy
  A = f.A
  p = size(Σx, 1)

  # ri, ci = ind2sub(Σx, j)
  ri, ci = Tuple(CartesianIndices(Σx)[j])

  a = zero(T)
  b = zero(T)
  @inbounds a  = Σx[ri,ri] * Σy[ci,ci]
  @inbounds b = A[ri,ci]
  b = ri == f.a && ci == f.b ? b - 1. : b

  @inbounds oldVal = x[j]
  a = one(T) / a
  @inbounds x[j] -= b * a
  newVal = cdprox!(g, x, j, a)
  h = newVal - oldVal

  # update internals
  for ac=1:p, ar=1:p
    @inbounds A[ar, ac] += h * Σx[ar, ri] * Σy[ci, ac]
  end
  h
end

############################################
#
#   asymmetric version of invHessian
#
############################################

function _computeVarElem(
    q, r,
    f::CDInverseKroneckerLoss,
    SxΘ, ΘSy,
    X, Y, row, col)

    nx, p = size(X)
    ny    = size(Y, 1)

    A = f.A
    a = f.a
    b = f.b

    δ = row == a && col == b ? - 1. : 0.

    fill!(q, 0.)
    fill!(r, 0.)

    # compute t_ab
    @inbounds t_ab = A[row, col] + δ

    # compute qk
    for k=1:nx
        v = 0.
        for l=1:p
            @inbounds v += X[k, l] * ΘSy[l, col]
        end
        @inbounds q[k] = X[k, row] * v + δ
    end

    # compute rk
    for k=1:ny
        v = 0.
        for l=1:p
            @inbounds v += SxΘ[row, l] * Y[k, l]
        end
        @inbounds r[k] = v * Y[k, col] + δ
    end

    σ1 = sum(abs2, q) / (nx - 1) - nx / (nx - 1) * t_ab^2
    σ2 = sum(abs2, r) / (ny - 1) - ny / (ny - 1) * t_ab^2

    return (σ1/nx + σ2/ny)

end

function _computeVar!(ω, f::CDInverseKroneckerLoss, X, Y, Θ)
    nx, p = size(X)
    ny    = size(Y, 1)

    q = zeros(nx)
    r = zeros(ny)

    SxΘ = _mul(f.Σx, Θ)
    ΘSy = _mul(Θ, f.Σy)

    for col=1:p
        for row=1:p
            linInd = (col-1)*p + row
            ω[linInd] = sqrt(_computeVarElem(q, r, f, SxΘ, ΘSy, X, Y, row, col))
        end
    end

end

function invAsymHessianEstimation(Sx::Symmetric, Sy::Symmetric, ri, ci, X, Y, λ, options=CDOptions())
    nx, p = size(X)
    ny = size(Y, 1)

    f = CDInverseKroneckerLoss(Sx, Sy, ri, ci)
    x = SparseIterate(p*p)
    x[(ci-1)*p + ri] = 1.
    CoordinateDescent.initialize!(f, x)

    ##################
    #
    #  first stage
    #
    ##################
    # compute initial variance
    ω = Array{eltype(Sx)}(undef, length(x))
    _computeVar!(ω, f, X, Y, x)
    ω[(ci-1)*p + ri] = 0.

    # compute initial estimate
    g = ProxL1(λ, ω)
    coordinateDescent!(x, f, g, options)

    ##################
    #
    #  second stage
    #
    ##################

    # recompute variance
    opt1 = CDOptions(;
        maxIter=options.maxIter,
        optTol=options.optTol,
        randomize=options.randomize,
        warmStart=true,
        numSteps=options.numSteps)

    _computeVar!(ω, f, X, Y, x)
    ω[(ci-1)*p + ri] = 0.

    # recompute estimate
    g = ProxL1(λ, ω)
    coordinateDescent!(x, f, g, opt1)

    return x
end


text:
using Random

abstract type NNModule end

"""
    forward(m, x)

Forward `x` through NNModule `m`.
"""
forward(m::NNModule, x) = throw("unimplemented")

"""
    params(m)

Return all parameters of `m`.
"""
params(m::NNModule) = throw("unimplemented")

function to!(m::NNModule, device::Device)
    for p ∈ params(m)
        to!(p, device)
    end
end

struct LinearLayer <: NNModule
    W::Tensor
    b::Union{Tensor,Nothing}
    bias::Bool

    function LinearLayer(input_dim::Integer, output_dim::Integer; bias::Bool=true)
        a = sqrt(1 / input_dim)
        new(
            Tensor(rand(input_dim, output_dim) * 2a .- a),
            bias ? Tensor(rand(1, output_dim) * 2a .- a) : nothing,
            bias,
        )
    end
end

forward(m::LinearLayer, x) = m.bias ? x * m.W .+ m.b : x * m.W
params(m::LinearLayer) = m.bias ? [m.W, m.b] : [m.W]

struct Activation <: NNModule
    activation::Function
end

forward(m::Activation, x) = m.activation(x)
params(m::Activation) = []

SigmoidActivation() = Activation(logistic)

struct ModuleList <: NNModule
    modules::Array{NNModule}
end

function forward(m::ModuleList, x)
    for mod ∈ m.modules
        x = forward(mod, x)
    end
    x
end

params(m::ModuleList) = vcat([params(mod) for mod ∈ m.modules]...)

max_stars_repo_name: taraspiotr/DeepJulia

max_issues_repo_name: taraspiotr/DeepJulia

max_forks_repo_name: taraspiotr/DeepJulia


text:
# AUTO GENERATED FILE - DO NOT EDIT

export 'feffery'_antdcascader

"""
    'feffery'_antdcascader(;kwargs...)

An AntdCascader component.

Keyword arguments:
- `id` (String; optional)
- `bordered` (Bool; optional)
- `changeOnSelect` (Bool; optional)
- `className` (String; optional)
- `defaultValue` (Array of Strings; optional)
- `disabled` (Bool; optional)
- `expandTrigger` (String; optional)
- `loading_state` (optional): . loading_state has the following type: lists containing elements 'is_loading', 'prop_name', 'component_name'.
Those elements have the following types:
  - `is_loading` (Bool; optional): Determines if the component is loading or not
  - `prop_name` (String; optional): Holds which property is loading
  - `component_name` (String; optional): Holds the name of the component that is loading
- `options` (optional)
- `placeholder` (String; optional)
- `popupPlacement` (String; optional)
- `size` (String; optional)
- `style` (Dict; optional)
- `value` (Array of Strings; optional)
"""
function 'feffery'_antdcascader(; kwargs...)
        available_props = Symbol[:id, :bordered, :changeOnSelect, :className, :defaultValue, :disabled, :expandTrigger, :loading_state, :options, :placeholder, :popupPlacement, :size, :style, :value]
        wild_props = Symbol[]
        return Component("'feffery'_antdcascader", "AntdCascader", "feffery_antd_components", available_props, wild_props; kwargs...)
end

max_stars_repo_name: Beddow/feffery-antd-components

max_issues_repo_name: Beddow/feffery-antd-components

max_forks_repo_name: Beddow/feffery-antd-components


text:
__precompile__(false)

module ScriptManager

mutable struct Script
  name::String
  initalized::Bool
  inputs::Dict{Symbol, Any}
  events::Dict{Symbol, Function}
  init::Union{Nothing, Function}
  update::Union{Nothing, Function}
  render::Union{Nothing, Function}
  clean::Union{Nothing, Function}

  Script(name::String) = new(name,false,Dict(),Dict(),nothing,nothing,nothing,nothing)
end

main_script = nothing
reload_script = nothing
setReload(func::Function) = global reload_script = func

function reload(script::Union{Nothing, Script}=main_script)
  if script == nothing return end
  reload_script(script)
  global main_script = script
  GC.gc()
  script.init()
end

getFunction(func::Union{Nothing, Function}, args...) = () -> begin if func != nothing eval(:($func($args...))) end; end

getEventCall(script::Script, key::Symbol) = haskey(script.events, key) ? script.events[key] : nothing

callEvent(script::Script, key::Symbol, args...) = getFunction(getEventCall(script, key), args...)()
callEvent(key::Symbol, args...) = if main_script!= nothing getFunction(getEventCall(main_script, key), args...)() end

end #ScriptManager

max_stars_repo_name: Gilga/JuliaOpenGL

max_issues_repo_name: Gilga/JuliaOpenGL

max_forks_repo_name: Gilga/JuliaOpenGL


text:
abstract type Shape{T} end

mutable struct Mesh{T} <: Shape{T}
    bodyids::Vector{Int64}
    m::T
    J::SMatrix{3,3,T,9}

    xoff::SVector{3,T}
    qoff::Quaternion{T}

    path::String
    color::RGBA

    function Mesh(path::String, m::T, J::AbstractMatrix{T};color = RGBA(0.75, 0.75, 0.75), xoff::AbstractVector{T} = zeros(3), qoff::Quaternion{T} = Quaternion{T}()) where T
        bodyids = Int64[]
        new{T}(bodyids, m, J, xoff, qoff, path, color)
    end
end

mutable struct Box{T} <: Shape{T}
    bodyids::Vector{Int64}
    m::T
    J::SMatrix{3,3,T,9}

    xoff::SVector{3,T}
    qoff::Quaternion{T}

    xyz::SVector{3,T}
    color::RGBA

    function Box(x, y, z, m::T;color = RGBA(0.75, 0.75, 0.75), xoff::AbstractVector{T} = zeros(3), qoff::Quaternion{T} = Quaternion{T}()) where T
        J = 1 / 12 * m * diagm([y^2 + z^2;x^2 + z^2;x^2 + y^2])
        bodyids = Int64[]

        new{T}(bodyids, m, J, xoff, qoff, [x;y;z], color)
    end
end

mutable struct Cylinder{T} <: Shape{T}
    bodyids::Vector{Int64}
    m::T
    J::SMatrix{3,3,T,9}

    xoff::SVector{3,T}
    qoff::Quaternion{T}

    rh::SVector{2,T}
    color::RGBA

    # by default the cylinder points in the z direction
    # TODO other direction
    function Cylinder(r, h, m::T;color = RGBA(0.75, 0.75, 0.75), xoff::AbstractVector{T} = zeros(3), qoff::Quaternion{T} = Quaternion{T}()) where T
        J = 1 / 2 * m * diagm([r^2 + 1 / 6 * h^2;r^2 + 1 / 6 * h^2;r^2])
        bodyids = Int64[]

        new{T}(bodyids, m, J, xoff, qoff, [r;h], color)
    end
end

mutable struct Sphere{T} <: Shape{T}
    bodyids::Vector{Int64}
    m::T
    J::SMatrix{3,3,T,9}

    xoff::SVector{3,T}
    qoff::Quaternion{T}

    r::T
    color::RGBA

    function Sphere(r, m::T;color = RGBA(0.75, 0.75, 0.75), xoff::AbstractVector{T} = zeros(3), qoff::Quaternion{T} = Quaternion{T}()) where T
        J = 2 / 5 * m * diagm([r^2 for i = 1:3])
        bodyids = Int64[]

        new{T}(bodyids, m, J, xoff, qoff, r, color)
    end
end


text:
module ChainRulesCore
using Base.Broadcast: broadcasted, Broadcasted, broadcastable, materialize, materialize!
using LinearAlgebra
using SparseArrays: SparseVector, SparseMatrixCSC
using Compat: hasfield

export frule, rrule  # core function
# rule configurations
export RuleConfig, HasReverseMode, NoReverseMode, HasForwardsMode, NoForwardsMode
export frule_via_ad, rrule_via_ad
# definition helper macros
export @non_differentiable, @scalar_rule, @thunk, @not_implemented
export canonicalize, extern, unthunk  # differential operations
export add!!  # gradient accumulation operations
# differentials
export Tangent, NoTangent, InplaceableThunk, Thunk, ZeroTangent, AbstractZero, AbstractThunk

include("compat.jl")
include("debug_mode.jl")

include("differentials/abstract_differential.jl")
include("differentials/abstract_zero.jl")
include("differentials/thunks.jl")
include("differentials/composite.jl")
include("differentials/notimplemented.jl")

include("differential_arithmetic.jl")
include("accumulation.jl")

include("config.jl")
include("rules.jl")
include("rule_definition_tools.jl")

include("deprecated.jl")

end # module


text:
module CGT

using CurrencyAmounts
using DataFrames
using Dates
using DelimitedFiles
using FixedPointDecimals
using HTTP
using XLSX

Fixed = FixedDecimal{Int, 2}

@currencies USD, EUR

const Tax = Fixed(0.33)
const Allowance = 1270EUR

struct ExchangeRates{From <: Currency, To <: Currency}
	cache::Dict{Date, ExchangeRate{<: Number, To, From}}
end
ExchangeRates(::From, ::To) where {From <: Currency, To <: Currency} = ExchangeRates{From, To}(Dict())

function fetch(date, from, to)
	url = "https://sdw-wsrest.ecb.europa.eu/service/data/EXR/D.$from.$to.SP00.A?detail=dataonly&startPeriod=$date&endPeriod=$date"
	(data, header) = readdlm(HTTP.get(url, ["Accept" => "text/csv"]).body, ','; header=true)
	Fixed(data[findfirst(isequal("OBS_VALUE"), header)])
end

function Base.getindex(rates::ExchangeRates{From, To}, date) where {From, To}
	if !haskey(rates.cache, date)
		from = From()
		to = To()
		rates.cache[date] = fetch(date, from, to) * from / to
	end
	rates.cache[date]
end

struct Sales{C <: Currency}
	df::DataFrame
end
Sales(c::Currency, df::DataFrame) = Sales{typeof(c)}(df)

function Base.convert(::Type{Sales{To}}, s::Sales{From}, rates::ExchangeRates{From, To}) where {To <: Currency, From <: Currency}
	exchange(date::Date, c::CurrencyAmount{<: Number, From}) =  convert(To(), c, rates[date])
	exchange(::Date, c::Any) = c
	transform(s.df, AsTable(:) => ByRow(r -> (map(c -> exchange(r.Date, c), r)..., Rate=rates[r.Date])) => AsTable) |> Sales{To}
end
Base.convert(t::Type{Sales{To}}, s::Sales{From}) where {To <: Currency, From <: Currency} = convert(t, s, ExchangeRates(From(), To()))

function load(file)
	fmt = DateFormat("m/d/y")
	to_usd(s) = Fixed(s)*USD

	df = DataFrame(XLSX.gettable(XLSX.readxlsx(file)[1])...)
	filter!("Record Type" => isequal("Sell"), df)
	select!(df, "Date Sold" => ByRow(d -> Date(d, fmt)) => "Date", "Plan Type" => "Type", "Total Proceeds" => ByRow(to_usd) => "Proceeds", "Adjusted Gain/Loss" => ByRow(to_usd) => "Gains")
	sort!(df, :Date)

	Sales(USD, df)
end

function period(c::Currency, df::AbstractDataFrame; details=false)
	gain = sum(df.Gains)
	cols = propertynames(df)

	if details
		show(df; allrows=true, allcols=true, summary=false, eltypes=false, show_row_number=false, vlines=:all, newline_at_end=true, hlines=[:begin, :header, :end], filters_col=((_, i) -> cols[i] != :Late,))
	end
	println("Gain/Loss: ", gain)
	println("Tax: ", max(0c, gain) * Tax)

	gain
end

function by_period(sales::Sales{C}; details=false) where {C <: Currency}
	if !(:Late in propertynames(sales.df))
		transform!(sales.df, :Date => ByRow(d -> Dates.month(d) == 12) => :Late)
	end

	year = Dates.year(sales.df[1, :Date])
	total_gain = 0 * C()
	for (key, sdf) in pairs(groupby(sales.df, :Late; sort=true))
		println("Period ", key.Late ? "$year-12-01 - $year-12-31" : "$year-01-01 - $year-11-30")
		total_gain += period(C(), sdf; details)
		println("")
	end

	total_gain
end

function compute(sales::Sales{Currency{:EUR}}; details=false)
	total_gain = by_period(sales; details)
	net_gain = total_gain >= Allowance ? total_gain - Allowance : total_gain
	tax = max(net_gain, 0EUR) * Tax

	println("Total gain/loss: ", total_gain)
	println("Net gain/loss excluding the allowance: ", net_gain)
	println("Total tax: ", tax)

	tax
end

end


text:
using Constructs
using Documenter

DocMeta.setdocmeta!(Constructs, :DocTestSetup, :(using Constructs); recursive=true)

makedocs(;
    modules=[Constructs],
    authors="Yong-an Lu <miroox@outlook.com>",
    repo="https://github.com/miRoox/Constructs.jl/blob/{commit}{path}#{line}",
    sitename="Constructs.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://miRoox.github.io/Constructs.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/miRoox/Constructs.jl",
)


text:
"""
    module MICS

MICS (Mixture of Independently Collected Samples) is a method for performing free-energy
computation and property reweighting using samples collected at multiple states. Its main
distinguishing feature is the possibility of employing all data from autocorrelated samples,
without sub-sampling.
"""
module MICS

using DataFrames
using DataFramesMeta
using Base.LinAlg.BLAS

export State,
       Mixture,
       freeEnergies,
       reweighting,
       covariance,
       multimap,
       @prop

include("aux.jl")

"""
    struct State

An equilibrium state aimed to be part of a mixture of independently collected samples (MICS)
"""
struct State
  sample::DataFrame          # properties of sampled configurations
  potential::Function        # reduced potential function
  n::Int                     # number of configurations
  b::Int                     # block size
  neff::Int                  # effective sample size
end

"""
# Constructor

    State( sample, potential[, autocorr] )

## Arguments

* `sample::DataFrames.DataFrame`: a data frame whose rows represent configurations sampled
     according to a given probability distribution and whose columns contain a number of
     properties evaluated for such configurations.
* `potential::Function`: the reduced potential that defines the equilibrium state. This
     function might for instance receive `x` and return the result of an element-wise
     calculation involving `x[:a]`, `x[:b]`, etc, with `a`, `b`, etc being names of
     properties in `sample`.
* `autocorr::Function=potential`: a function similar to `potential`, but whose result is an
     autocorrelated property to be used for determining the effective sample size.

## Note

Formally, functions `potential` and `autocorr` must receive `x::DataFrames.DataFrame` and
return `y::T`, where `T<:AbstractArray{Float64,1}` with `length(y) == nrow(x)`.
"""
function State( sample, potential, autocorr=potential )
  n = nrow(sample)
  b = round(Int,sqrt(n))
  y = multimap( [autocorr], sample )
  ym = mean(y,1)
  neff = n*covariance(y,ym,1)[1]/covariance(y,ym,b)[1]
  isfinite(neff) || error( "unable to determine effective sample size" )
  State( sample, potential, n, b, round(Int,neff) )
end

"""
    struct Mixture

A mixture of independently collected samples (MICS)
"""
struct Mixture
  title::String
  state::Vector{State}
  names::Vector{Symbol}
  m::Int                       # number of states
  n::Vector{Int}               # sample size at each state
  π::Vector{Float64}           # mixture composition
  f::Vector{Float64}           # free energy of each state
  P::Vector{Matrix{Float64}}   # probabilities of each configuration at all states
  pm::Vector{Matrix{Float64}}  # mean probabilities at all states
  u0::Vector{Matrix{Float64}}  # reduced energy of each configuration at the mixture state
  B0⁺::Matrix{Float64}
  Θ::Matrix{Float64}
end

"""
# Constructor

    Mixture( states; <keyword arguments> )

## Arguments

* `states::Vector{MICS.State}`:
* `title::String=\"Untitled\"`:
* `verbose::Bool=false`:
* `tol::Float64=1.0E-8`:
"""
function Mixture( states::Vector{State}; title::String = "Untitled",
                  verbose::Bool = false, tol::Float64 = 1.0E-8 )

  verbose && aux.info( "Setting up MICS case: ", title )

  m = length(states)
  verbose && aux.info( "Number of states: ", m )
  m == 0 && error( "state set is empty" )

  properties = names(states[1].sample)
  verbose && aux.info( "Properties: ", aux.str(properties) )
  all([all(names(states[i].sample) .== properties) for i=1:m]) ||
    error( "inconsistent data" )

  n = [nrow(states[i].sample) for i=1:m]
  verbose && aux.info( "Sample sizes: ", aux.str(n) )

  neff = [states[i].neff for i=1:m]
  verbose && aux.info( "Effective sample sizes: ", aux.str(neff) )

  π = neff/sum(neff)
  verbose && aux.info( "Mixture composition: ", π )

  potentials = [states[i].potential for i=1:m]
  u = [multimap(potentials, states[i].sample) for i=1:m]

  P = Vector{Matrix{Float64}}(m)
  u0 = Vector{Matrix{Float64}}(m)
  for i = 1:m
    P[i] = Matrix{Float64}(n[i],m)
    u0[i] = Matrix{Float64}(n[i],1)
  end

  f = aux.overlapSampling( u )
  verbose && aux.info( "Initial free-energy guess:", f )

  iter = 1
  aux.compute!( u0, P, π, f, u )
  p0 = [mean(P[j][:,i]) for i=1:m, j=1:m]*π
  B0 = Symmetric(diagm(p0) - sum(syrk('U', 'T', π[i]/n[i], P[i]) for i=1:m))
  Δf = B0[2:m,2:m]\(π - p0)[2:m]
  while any(abs.(Δf) .> tol)
    iter += 1
    f[2:m] += Δf
    aux.compute!( u0, P, π, f, u )
    p0 = [mean(P[j][:,i]) for i=1:m, j=1:m]*π
    B0 = Symmetric(diagm(p0) - sum(syrk('U', 'T', π[i]/n[i], P[i]) for i=1:m))
    Δf = B0[2:m,2:m]\(π - p0)[2:m]
  end
  verbose && aux.info( "Free energies after $(iter) iterations:", f )

  (D, V) = eig(B0)
  D⁺ = diagm(map(x-> abs(x) < tol ? 0.0 : 1.0/x, D))
  B0⁺ = Symmetric(V*D⁺*V')
  pm = [mean(P[i],1) for i=1:m]
  Σ0 = sum(π[i]^2*covariance(P[i], pm[i], states[i].b) for i=1:m)
  Θ = Symmetric(B0⁺*Σ0*B0⁺)
  verbose && aux.info( "Free-energy covariance matrix:", full(Θ) )

  Mixture( title, states, properties, m, n, π, f, P, pm, u0, B0⁺, Θ )
end

"""
    freeEnergies( mixture )

Returns a data frame containing the relative free energies of the sampled states of a
`mixture`, as well as their standard errors.
"""
function freeEnergies( mixture::Mixture )
  δf = sqrt.([mixture.Θ[i,i] - 2*mixture.Θ[i,1] + mixture.Θ[1,1] for i=1:mixture.m])
  return DataFrame( f = mixture.f, δf = δf )
end

"""
    reweighting( mixture, functions, potential, parameter )

Performs reweighting of the properties computed by `functions` from the mixture to the
states determined by the provided `potential` with all `parameter` values.
"""
function reweighting( mixture::Mixture,
                      properties::Union{Function,Vector{T}} where T<:Function,
                      potential::Function,
                      parameter::AbstractArray )

  state = mixture.state
  m = mixture.m
  π = mixture.π
  P = mixture.P
  pm = mixture.pm

  # Compute z so that z[i] contains a matrix of all properties evaluated for sample i:
  z = [hcat(ones(state[i].n), multimap( properties, state[i].sample )) for i=1:m]

  np = length(parameter)
  y0 = Vector{Vector{Float64}}(np)
  Ξ = Vector{Matrix{Float64}}(np)

  # Carry out reweighting with every parameter value:
  for k = 1:np
    u = x->potential(x,parameter[k])
    Δu = [mixture.u0[i] - multimap( u, state[i].sample ) for i=1:m]
    maxΔu = maximum([maximum(Δu[i]) for i=1:m])
    y = [exp.(Δu[i] - maxΔu) .* z[i] for i=1:m]
    ym = [mean(y[i],1) for i=1:m]
    Σy0y0 = sum(π[i]^2*covariance(y[i], ym[i], state[i].b) for i=1:m)
    Σp0y0 = sum(π[i]^2*covariance(P[i], pm[i], y[i], ym[i], state[i].b) for i=1:m)
    Z0 = -sum(gemm('T', 'N', π[i]/state[i].n, P[i], y[i]) for i=1:m)
    A = Z0'*mixture.B0⁺*Σp0y0
    y0[k] = vec(sum(π.*ym))
    Ξ[k] = Σy0y0 + A + A' + Z0'*mixture.Θ*Z0
  end

end

"""
    covariance( y, ym, [ z, zm,] b )

Computes either the covariance matrix of the columns of matrix `y` among themselves or the
cross-covariance matrix between the columns of matrix `y` with those of matrix `z`. The
method of Overlap Batch Mean (OBM) is employed with blocks of size `b`.
"""
function covariance( y::Matrix{T}, ym::Matrix{T}, b::Integer ) where T<:AbstractFloat
  S = aux.SumOfDeviationsPerBlock( y, ym, b )
  nmb = size(y,1) - b
  return Symmetric(syrk('U', 'T', 1.0/(b*nmb*(nmb+1)), S))
end

function covariance( y::Matrix{T}, ym::Matrix{T},
                     z::Matrix{T}, zm::Matrix{T}, b::Integer ) where T<:AbstractFloat
  Sy = aux.SumOfDeviationsPerBlock( y, ym, b )
  Sz = aux.SumOfDeviationsPerBlock( z, zm, b )
  nmb = size(y,1) - b
  return gemm('T', 'N', 1.0/(b*nmb*(nmb+1)), Sy, Sz)
end

"""
    multimap( functions, frame )

Applies an array of `functions` to a data `frame` and returns a matrix whose number of rows
is the same as in `frame` and the number of columns is equal to the length of `functions`.

# Note
Each function of the array might for instance receive `x` and return the result of an
element-wise calculation involving `x[:a]`, `x[:b]`, etc, with `a`, `b`, etc being names of
properties in `frame`.
"""
function multimap( functions::Array{T} where T <: Function, frame::DataFrame )
  m = length(functions)
  f = Matrix{Float64}(nrow(frame),m)
  for i in eachindex(functions)
    f[:,i] = functions[i]( frame )
  end
  return f
end

multimap( f::Function, frame::DataFrame ) = multimap( [f], frame )

"""
    @prop( expr )
"""
macro prop( expr )
  esc(DataFramesMeta.with_anonymous(expr))
end

end


text:
"""
This module demonstrates how to create your first window using GLFW.
"""
module FirstWindow

using GLFW
using ModernGL

"""
  main()

Application entrypoint for a bare minimum demo which shows how to create a
window using GLFW.
"""
function main()
  # Specify minimum versions
  GLFW.WindowHint(GLFW.CONTEXT_VERSION_MAJOR, 4)
  GLFW.WindowHint(GLFW.CONTEXT_VERSION_MINOR, 6)
  GLFW.WindowHint(GLFW.OPENGL_PROFILE, GLFW.OPENGL_CORE_PROFILE)
  GLFW.WindowHint(GLFW.OPENGL_FORWARD_COMPAT, GL_TRUE)

  # Create the window
  window = GLFW.CreateWindow(800, 600, "Context creation")
  GLFW.ShowWindow(window)
  GLFW.MakeContextCurrent(window)

  GLFW.SetInputMode(window, GLFW.STICKY_KEYS, GL_TRUE)

  while !GLFW.WindowShouldClose(window)
      GLFW.SwapBuffers(window)
      GLFW.PollEvents()
      if GLFW.GetKey(window, GLFW.KEY_ESCAPE) == GLFW.PRESS
          GLFW.SetWindowShouldClose(window, true)
      end
  end

  GLFW.DestroyWindow(window)
end

end # Module FirstWindow


text:
using Base.Test
using Einsum

## Test that vars in Main aren't overwritten by einsum
let
  i = -1
  y = randn(10)
  @einsum x[i] := y[i]
  @test i == -1
end

## Test that B is overwritten by := operator
let
  B = randn(10,10)
  A = randn(5,10)
  @einsum B[i,j] := A[i,j] # this should run without a problem
  @test size(B) == size(A)
end

## CP decomposition test case ##
let

  # preallocated test case
  A = zeros(5,6,7);
  B = zeros(5,6,7);
  X = randn(5,2);
  Y = randn(6,2);
  Z = randn(7,2);

  @einsum A[i,j,k] = X[i,r]*Y[j,r]*Z[k,r]
  @einsimd B[i,j,k] = X[i,r]*Y[j,r]*Z[k,r]

  for i = 1:5
    for j = 1:6
      for k = 1:7
        s = 0.0
        for r = 1:2
          s += X[i,r]*Y[j,r]*Z[k,r]
        end
        @test isapprox(A[i,j,k],s)
        @test isapprox(B[i,j,k],s)
      end
    end
  end

  # without preallocation
  @einsum A2[i,j,k] := X[i,r]*Y[j,r]*Z[k,r]
  @test isapprox(A,A2)

end

# Interesting test case, can throw an error that
# local vars are declared twice.
let
  A = zeros(5,6,7);
  X = randn(5,2);
  Y = randn(6,2);
  Z = randn(7,2);
  if true
    @einsum A[i,j,k] = X[i,r]*Y[j,r]*Z[k,r]
  else
    @einsum A[i,j,k] = X[i,r]*Y[j,r]*Z[k,r]
  end
end

# At one point this threw an error because the lhs
# had no indices/arguments
let
  x = randn(10)
  y = randn(10)
  @einsum k := x[i]*y[i]
  @test isapprox(k,dot(x,y))
end

# Elementwise multiplication (this should create nested loops with no
# no summation.)
let
  x = randn(10)
  y = randn(10)
  @einsum k[i] := x[i]*y[i]
  @test isapprox(k,x.*y)
end

# Transpose a block matrix
let
  z = Any[ rand(2,2) for i=1:2, j=1:2]
  @einsum t[i,j] := transpose(z[j,i])
  @test isapprox(z[1,1], t[1,1]')
  @test isapprox(z[2,2], t[2,2]')
  @test isapprox(z[1,2], t[2,1]')
  @test isapprox(z[2,1], t[1,2]')
end

# Mapping functions
let
  A = randn(10,10)
  @einsum B[i,j] := exp(A[i,j])
  @test isapprox(exp(A),B)
end

# Example from numpy
let
  A = reshape(collect(1:25),5,5)
  @einsum B[i] := A[i,i]
  @test all(B .== [1,7,13,19,25])
end

# TODO: consider adding support for this:
# @einsum A[i,j] = A[i,j] + 50

## Test in-place operations ##
let
  A = randn(5,6,7);
  B = randn(5,6,7);
  A1 = copy(A);
  B1 = copy(B);

  X = randn(5,2);
  Y = randn(6,2);
  Z = randn(7,2);

  @einsum A[i,j,k] += X[i,r]*Y[j,r]*Z[k,r]
  @einsimd B[i,j,k] += X[i,r]*Y[j,r]*Z[k,r]

  for i = 1:5
      for j = 1:6
          for k = 1:7
              s = 0.0
              for r = 1:2
                  s += X[i,r]*Y[j,r]*Z[k,r]
              end
              @test isapprox(A[i,j,k],A1[i,j,k]+s)
              @test isapprox(B[i,j,k],B1[i,j,k]+s)
          end
      end
  end

  x = randn(10)
  y = randn(10)
  k0 = randn()
  k = k0
  @einsum k += x[i]*y[i]
  @test isapprox(k,k0+dot(x,y))

  # test multiplication

  A1[:] = A[:]
  B1[:] = B[:]

  @einsum A[i,j,k] *= X[i,r]*Y[j,r]*Z[k,r]
  @einsimd B[i,j,k] *= X[i,r]*Y[j,r]*Z[k,r]

  for i = 1:5
      for j = 1:6
          for k = 1:7
              s = 0.0
              for r = 1:2
                  s += X[i,r]*Y[j,r]*Z[k,r]
              end
              @test isapprox(A[i,j,k],A1[i,j,k]*s)
              @test isapprox(B[i,j,k],B1[i,j,k]*s)
          end
      end
  end

  x = randn(10)
  y = randn(10)
  k0 = randn()
  k = k0
  @einsum k *= x[i]*y[i]
  @test isapprox(k,k0*dot(x,y))
end

# Test offsets
let
  X = randn(10)

  # without preallocation
  @einsum A[i] := X[i+5]
  @test size(A) == (5,)
  @test all(A .== X[6:end])

  # with preallocation
  B = zeros(10)
  @einsum B[i] = X[i+5]
  @test size(B) == (10,)
  @test all(B[1:5] .== X[6:end])
end

# Test symbolic offsets
let
  offset = 5
  X = randn(10)

  # without preallocation
  @einsum A[i] := X[i+:offset]
  @test size(A) == (5,)
  @test all(A .== X[6:end])

  # with preallocation
  B = zeros(10)
  @einsum B[i] = X[i+:offset]
  @test size(B) == (10,)
  @test all(B[1:5] .== X[6:end])
end

# Test adding/subtracting constants
let
  k = 5
  X = randn(10)

  # without preallocation
  @einsum A[i] := X[i] + k
  @einsum B[i] := X[i] - k
  @test isapprox(A,X+k)
  @test isapprox(B,X-k)

  # with preallocation
  C,D = zeros(10),zeros(10)
  @einsum C[i] = X[i] + k
  @einsum D[i] = X[i] - k
  @test isapprox(C,X+k)
  @test isapprox(D,X-k)
end

# Test multiplying/dividing constants
let
  k = 5
  X = randn(10)

  # without preallocation
  @einsum A[i] := X[i]*k
  @einsum B[i] := X[i]/k
  @test isapprox(A,X.*k)
  @test isapprox(B,X./k)

  # with preallocation
  C,D = zeros(10),zeros(10)
  @einsum C[i] = X[i]*k
  @einsum D[i] = X[i]/k
  @test isapprox(C,X.*k)
  @test isapprox(D,X./k)
end

# Test indexing with a constant
let
  A = randn(10,2)
  j = 2
  @einsum B[i] := A[i,:j]
  @test all(B .== A[:,j])
  @einsum C[i] := A[i,1]
  @test all(C .== A[:,1])

  D = zeros(10,3)
  @einsum D[i,1] = A[i,:j]
  @test isapprox(D[:,1],A[:,j])
  @einsum D[i,:j] = A[i,:j]
  @test isapprox(D[:,j],A[:,j])
end

# Better type inference on allocating arrays
let
  B1 = ones(Int,5)
  B2 = ones(Float32,5)
  B3 = ones(5)
  C = randn(5)
  @einsum A1[i,j] := B1[i]*C[j]
  @einsum A2[i,j] := B2[i]*C[j]
  @einsum A3[i,j] := B3[i]*C[j]

  @test eltype(A1) == Float64
  @test eltype(A2) == Float64
  @test eltype(A3) == Float64
  @test isapprox(A1,A3)
  @test isapprox(A2,A3)
end


text:
module EntityManager

export Entity

using CoreExtended
using MatrixMath

type Entity
	id::Symbol
	enabled::Bool

  position::Vec3f
  rotation::Vec3f
	scaling::Vec3f

	Entity(id=:_) = new(id,true,zeros(Vec3f),zeros(Vec3f),ones(Vec3f))
end

# preset for manager: includes code here
presetManager(EntityManager,Entity)

end #EntityManager

max_stars_repo_name: Gilga/JGE

max_issues_repo_name: Gilga/JGE

max_forks_repo_name: Gilga/JGE


text:
HOME_DIR = dirname(@__DIR__)
PACKAGE_DIR = dirname(HOME_DIR)
pushfirst!(LOAD_PATH, PACKAGE_DIR)
using MPCHelper

using ControlSystems
using LinearAlgebra
using LazySets

using Test

using Plots

## Initialization ---------------------------------------------------------------
# Linear 2D System
A = [1 -2; 0 1]
B= Matrix([1 1]')
R = 0.1*I
Q = Matrix(I,size(A))/10
K = dlqr(A,B,Q,R)
A_stab = A - B*K; eigvals(A_stab)
P = dlyap(A_stab, Q)
l_x = [3,3]
l_u = 1
z = [0.0 ;0.0]; v = 0.0

# # Linear 3D System
# A = [0.8 -2 3; 0 0.5 0; 0 0 0.3]
# B= [1 0; 0 1; 0 0.4]
# R = 0.1*Matrix(I,size(B,2), size(B,2))
# Q = 0.1*Matrix(I,size(A))
# K = dlqr(A,B,Q,R)
# A_stab = A - B*K; eigvals(A_stab)
# P = dlyap(A_stab, Q)
# l_x = [3,3,3]
# l_u = [2,2]
# z = [0.0 ;0.0; 0.0]; v = [0.0, 0.0]

# Main Part
α = 2
E = get_Ellipsoid(P, α)
l(x,u) = dot(l_x,x) + dot(l_u,u)
γ = l_x + K'*l_u
ωstar = ω_opt(γ, P, α)
hΩ = hΩ_opt(γ, P, α)
# Scale γ -> vec_γ, such that vec_γ'*γ the same cost as ω*γ
vec_γ = γ*hΩ/norm(γ)^2
@test ωstar'*γ ≈ vec_γ'*γ
# lmax(z,v) for Ellispoid defined by P, α
lmax(z,v) = lmax_mask(z,v, P, α, l_x, l_u, K)
# lmax(z,v) for Ellispoid defined by Ellipsoid E (Type LazySet)
lmax_lazy(z,v) = lmax_mask(z,v, E, l_x, l_u, K)
lmax_opt = lmax(z,v)
lmax2_opt = lmax_lazy(z,v)
@test hΩ ≈ lmax_opt
x_opt = z .+ ωstar
u_opt = vec(v .+ K*ωstar)
@test l(x_opt, u_opt) ≈ lmax_opt

## Plotting (only for 2D systems)
grid_res = -5:0.2:5
plot(E,1e-3, alpha=0.5, lab="Omega", size=(500,500))
plot_vec!(ωstar, line=(3,:blue), lab="w_opt")
plot_vec!(vec_γ, line=(2,:black, :dash), lab="lx + K'lu")
# plot_vec!([1,1], line=(2,:black, :dash), lab="lx + K'lu")
plot_lin_contour!(l_x, grid_res; color=:bluesreds, line=(1,:dash), colorbar=false)
plot_lin_contour!(γ, grid_res; color=:bluesreds, line=(1))#, contour_labels = true)
l_opt(x) = (-γ[1]*x + hΩ) / γ[2] # γ'x = hα
plot!(grid_res, l_opt.(grid_res), line=(:red), lab="lopt")
xlims!(-5,5)
ylims!(-5,5)

## PLotting Surfaces
grid_res = 0:0.1:10
surface(grid_res,grid_res, (x,y) -> lmax([x,y],0) )
surface!(grid_res,grid_res, (x,y) -> l([x,y],0))

## Solve OP for lmax and compare to analytical solution
using JuMP
using Ipopt
import LazySets.Approximations.overapproximate
n = size(A,2)

function lmax_optimized(z,v)
    HPolyrep_epsilon = overapproximate(E, 1e-5) # approximate ellipse by polytope
    model = Model(with_optimizer(Ipopt.Optimizer))
    JuMP.register(model, :l, 2, l, autodiff=true )
    @variable(model, ω[1:n])
    constrain_variable_by_set2!(model, ω, HPolyrep_epsilon)
    @objective(model, Max, l(z+ω, K*ω .+ v))
    # print(model)
    optimize!(model)
    opt_val = objective_value(model)
    opt_ω = value.(ω)
    @show opt_ω, opt_val
    return opt_val
end

lmax_optimized(z,v) - lmax(z,v)
@test abs(lmax_optimized(z,v) - lmax(z,v)) <= 1e-4


text:
## Packages and settings.
include("../src/StencilBuilder.jl")

using BenchmarkTools
using LoopVectorization
using .StencilBuilder

float_type = Float32

## CPU dynamics kernel.
function dynamics_w_kernel!(
    wt, u, v, w, s,
    visc,
    dxi, dyi, dzi, dzhi,
    is, ie, js, je, ks, ke)

    alpha = 9.81/300;

    @fast3d begin
        @fd (wt, u, v, w, s) wt += (
            - gradx(interpz(u) * interpx(w)) + visc * (gradx(gradx(w)))
            - grady(interpz(v) * interpy(w)) + visc * (grady(grady(w)))
            - gradz(interpz(w) * interpz(w)) + visc * (gradz(gradz(w)))
            + alpha*interpz(s) )
    end
end

## User input
itot = 384; jtot = 384; ktot = 384;
igc = 1; jgc = 1; kgc = 1;

icells = itot + 2igc; jcells = jtot + 2jgc; kcells = ktot + 2kgc
is = igc+1; ie = igc+itot; js = jgc+1; je = jgc+jtot; ks = kgc+1; ke = kgc+ktot

wt = rand(float_type, icells, jcells, kcells)
u = rand(float_type, icells, jcells, kcells)
v = rand(float_type, icells, jcells, kcells)
w = rand(float_type, icells, jcells, kcells)
s = rand(float_type, icells, jcells, kcells)
dxi = rand(float_type)
dyi = rand(float_type)
dzi = rand(float_type, kcells)
dzhi = rand(float_type, kcells)
visc = convert(float_type, 1)

## Run CPU kernel
@btime dynamics_w_kernel!(
    $wt, $u, $v, $w, $s,
    $visc,
    $dxi, $dyi, $dzi, $dzhi,
    $is, $ie, $js, $je, $ks, $ke)


text:
function RunDebugIDE(in)
  global interp = in
  tocall!(in)
  debugmode(true)
  stepto(in)
  wait(cond)
  finish!(in)
  in.retval
end

# Manual Entry

function step(args...)
  global interp
  RunDebugIDE(enter_call_expr(nothing, :($(args...)())))
end

macro step(ex)
  @capture(ex, f_(args__)) || error("Syntax: @step f(...)")
  :(step($(esc(f)), $(map(esc, args)...)))
end

# Breakpoints

using Gallium
import Gallium: stackwalk, process_lowlevel_conditionals, bps_at_location,
  Location, LocalSession, CStackFrame, JuliaStackFrame, Hooking, FileLineSource

import Base.Meta: quot

function breakpoint(args...)
  Gallium.breakpoint(args...)
  return
end

Gallium.breakpoint_hit(hook, RC) = _breakpoint_hit(hook, RC)

_breakpoint_hit = function (hook, RC)
  if !process_lowlevel_conditionals(Location(LocalSession(), hook.addr), RC)
    return
  end
  stack = stackwalk(RC; fromhook = true)[1]
  stacktop = pop!(stack)
  linfo = stacktop.linfo
  fT = linfo.def.sig.parameters[1]
  def_linfo = linfo.def.lambda_template
  argnames = [Symbol("#target_self#");def_linfo.slotnames[2:def_linfo.nargs]]
  spectypes = [fT;linfo.specTypes.parameters[2:end]...]
  bps = bps_at_location[Location(LocalSession(),hook.addr)]
  target_line = minimum(map(bps) do bp
      idx = findfirst(s->isa(s, FileLineSource), bp.sources)
      idx != 0 ? bp.sources[idx].line : def_linfo.def.line
  end)
  conditions = reduce(vcat,map(bp->bp.conditions, bps))
  thunk = Expr(:->,Expr(:tuple,map(x->Expr(:(::),x[1],x[2]),zip(argnames,spectypes))...),Expr(:block,
    :(linfo = $(quot(linfo))),
    :((loctree, code) = ASTInterpreter.reparse_meth(linfo)),
    :(__env = ASTInterpreter.prepare_locals(linfo.def.lambda_template)),
    :(copy!(__env.sparams, linfo.sparam_vals)),
    [ :(__env.locals[$i] = Nullable{Any}($(argnames[i]))) for i = 1:length(argnames) ]...,
    :(interp = ASTInterpreter.enter(linfo,__env,
      $(collect(filter(x->!isa(x,CStackFrame),stack)));
        loctree = loctree, code = code)),
    (target_line != linfo.def.line ?
      :(ASTInterpreter.advance_to_line(interp, $target_line)) :
      :(nothing)),
    # :(tty_state = Gallium.suspend_other_tasks()),
    :((isempty($conditions) ||
      any(c->Gallium.matches_condition(interp,c),$conditions)) &&
      Atom.Debugger.RunDebugIDE(interp)),
    # :(Gallium.restore_other_tasks(tty_state)),
    :(ASTInterpreter.finish!(interp)),
    :(return interp.retval::$(linfo.rettype))))
  f = eval(thunk)
  t = Tuple{spectypes...}
  faddr = Hooking.get_function_addr(f, t)
  Hooking.Deopt(faddr)
end


text:
import Plots
import OrderedCollections
using Test

import Dates
import JSON

import NCDatasets
const NC = NCDatasets

import StatsBase
import Dierckx
import PrettyTables
import Printf

import ArtifactWrappers
const AW = ArtifactWrappers

ENV["GKSwstype"] = "nul"

# Get PyCLES_output dataset folder:
#! format: off
PyCLES_output_dataset = AW.ArtifactWrapper(
    @__DIR__,
    isempty(get(ENV, "CI", "")),
    "PyCLES_output",
    AW.ArtifactFile[
    AW.ArtifactFile(url = "https://caltech.box.com/shared/static/johlutwhohvr66wn38cdo7a6rluvz708.nc", filename = "Rico.nc",),
    AW.ArtifactFile(url = "https://caltech.box.com/shared/static/zraeiftuzlgmykzhppqwrym2upqsiwyb.nc", filename = "Gabls.nc",),
    AW.ArtifactFile(url = "https://caltech.box.com/shared/static/toyvhbwmow3nz5bfa145m5fmcb2qbfuz.nc", filename = "DYCOMS_RF01.nc",),
    AW.ArtifactFile(url = "https://caltech.box.com/shared/static/ivo4751camlph6u3k68ftmb1dl4z7uox.nc", filename = "TRMM_LBA.nc",),
    AW.ArtifactFile(url = "https://caltech.box.com/shared/static/4osqp0jpt4cny8fq2ukimgfnyi787vsy.nc", filename = "ARM_SGP.nc",),
    AW.ArtifactFile(url = "https://caltech.box.com/shared/static/jci8l11qetlioab4cxf5myr1r492prk6.nc", filename = "Bomex.nc",),
    AW.ArtifactFile(url = "https://caltech.box.com/shared/static/pzuu6ii99by2s356ij69v5cb615200jq.nc", filename = "Soares.nc",),
    AW.ArtifactFile(url = "https://caltech.box.com/shared/static/7upt639siyc2umon8gs6qsjiqavof5cq.nc", filename = "Nieuwstadt.nc",),
    ],
)
PyCLES_output_dataset_path = AW.get_data_folder(PyCLES_output_dataset)

SCAMPy_output_dataset = AW.ArtifactWrapper(
    @__DIR__,
    isempty(get(ENV, "CI", "")),
    "SCAMPy_output",
    AW.ArtifactFile[
    AW.ArtifactFile(url = "https://caltech.box.com/shared/static/1dzpydqiagjvzfpyv9lbic3atvca93hl.nc", filename = "Rico.nc",),
    AW.ArtifactFile(url = "https://caltech.box.com/shared/static/agkycoum93bd6xjduyaeo3oqg6asoru5.nc", filename = "GABLS.nc",),
    AW.ArtifactFile(url = "https://caltech.box.com/shared/static/fqpq1q74uxfh1e8018hwhqogw1htn2lq.nc", filename = "DYCOMS_RF01.nc",),
    AW.ArtifactFile(url = "https://caltech.box.com/shared/static/wevi0rqiwo6sgkqdhcddr72u5ylt0tqp.nc", filename = "TRMM_LBA.nc",),
    AW.ArtifactFile(url = "https://caltech.box.com/shared/static/fis6n0g9x9lts70m0zmve5ullqnw0pzq.nc", filename = "ARM_SGP.nc",),
    AW.ArtifactFile(url = "https://caltech.box.com/shared/static/t6qq6plt2oxcmmy40r1szgokahykqggp.nc", filename = "Bomex.nc",),
    AW.ArtifactFile(url = "https://caltech.box.com/shared/static/wp8k4m7ta1hs0c6e4j2fpsp3kj05wdip.nc", filename = "Soares.nc",),
    AW.ArtifactFile(url = "https://caltech.box.com/shared/static/vbuxzg85scwy9mg0ziiuzy6fimo0e389.nc", filename = "Nieuwstadt.nc",),
    AW.ArtifactFile(url = "https://caltech.box.com/shared/static/72t6fr1gq10tg3jjputtp35nfzex0o4k.nc", filename = "DryBubble.nc",),
    AW.ArtifactFile(url = "https://caltech.box.com/shared/static/7axeussneeg8g3k0ndvagsn0pkmbij3e.nc", filename = "life_cycle_Tan2018.nc",),
    AW.ArtifactFile(url = "https://caltech.box.com/shared/static/r6t7dk6g35bmbvc86h006yusb6rd8vkw.nc", filename = "SP.nc",),
    ],
)
SCAMPy_output_dataset_path = AW.get_data_folder(SCAMPy_output_dataset)
#! format: on

function get_time(ds, var)
    if haskey(ds, var)
        return ds[var][:]
    elseif haskey(ds.group["timeseries"], var)
        return ds.group["profiles"][var][:]
    end
    error("No key for $var found in the nc file.")
end

append_dict(::Nothing, sym::Symbol, dict = Dict()) = dict
function append_dict(filename::AbstractString, sym::Symbol, dict = Dict())
    if isfile(filename)
        dict[sym] = NC.Dataset(filename, "r")
    end
    return dict
end

function find_latest_dataset_folder(; dir = pwd())
    matching_paths = String[]
    for file in readdir(dir)
        !ispath(joinpath(dir, file)) && continue
        push!(matching_paths, joinpath(dir, file))
    end
    isempty(matching_paths) && return ""
    # sort by timestamp
    sorted_paths = sort(matching_paths; by = f -> Dates.unix2datetime(stat(f).mtime))
    return pop!(sorted_paths)
end

function compute_mse_wrapper(
    case_name,
    best_mse,
    ds_tc_filename,
    ds_tc_main_filename = nothing;
    ds_les_filename = joinpath(PyCLES_output_dataset_path, "$case_name.nc"),
    ds_scm_filename = joinpath(SCAMPy_output_dataset_path, "$case_name.nc"),
    plot_dir = joinpath(dirname(ds_tc_filename), "comparison"),
    kwargs...,
)
    # Note: cluster_data_prefix is also defined in utils/move_output.jl
    if haskey(ENV, "BUILDKITE_COMMIT")
        cluster_data_prefix = "/central/scratch/climaci/turbulenceconvection-main"
        path = find_latest_dataset_folder(; dir = cluster_data_prefix)

        # TODO: make this more robust in case folder/file changes
        folder_name = joinpath("Output.$case_name.01", "stats")
        ds_tc_main_filename = joinpath(path, folder_name, "Stats.$case_name.nc")
        @info "TC.jl main dataset: $ds_tc_main_filename"
    end
    # Note that we're using a closure over
    #  - PyCLES_output_dataset_path
    #  - SCAMPy_output_dataset_path
    all_args = (case_name, best_mse, plot_dir)

    ds_dict = append_dict(ds_les_filename, :ds_pycles)
    ds_dict = append_dict(ds_scm_filename, :ds_scampy, ds_dict)
    ds_dict = append_dict(ds_tc_filename, :ds_tc, ds_dict)
    ds_dict = append_dict(ds_tc_main_filename, :ds_tc_main, ds_dict)

    computed_mse_success = false

    local computed_mse
    try
        computed_mse = compute_mse(all_args...; ds_dict, kwargs...)
        computed_mse_success = true
    finally
        for ds in values(ds_dict)
            close(ds)
        end
    end
    if computed_mse_success
        return computed_mse
    else
        error("MSE computation unsuccessful.")
    end
end

function compute_mse(case_name, best_mse, plot_dir; ds_dict, plot_comparison = true, group_figs = true, t_start, t_stop)

    TC = TurbulenceConvection
    ds_scampy = haskey(ds_dict, :ds_scampy) ? ds_dict[:ds_scampy] : nothing
    ds_tc_main = haskey(ds_dict, :ds_tc_main) ? ds_dict[:ds_tc_main] : nothing
    ds_pycles = haskey(ds_dict, :ds_pycles) ? ds_dict[:ds_pycles] : nothing
    ds_tc = haskey(ds_dict, :ds_tc) ? ds_dict[:ds_tc] : nothing
    # Make pycles optional
    have_pycles_ds = ds_pycles ≠ nothing
    have_scampy_ds = ds_scampy ≠ nothing
    have_tc_main = ds_tc_main ≠ nothing
    if !have_pycles_ds
        ds_pycles = ds_tc
        @warn "No PyCLES data, using TC.jl instead"
    end
    if !have_scampy_ds
        ds_scampy = ds_tc
        @warn "No SCAMPy data, using TC.jl instead"
    end
    if !have_tc_main
        ds_tc_main = ds_tc
        @warn "No TC.jl main data, using TC.jl instead"
    end
    mse = OrderedCollections.OrderedDict()
    time_tcc = get_time(ds_tc, "t")
    time_les = get_time(ds_pycles, "t")
    time_tcm = get_time(ds_tc_main, "t")
    time_scm = get_time(ds_scampy, "t")
    n_time_points = length(time_tcc)

    mkpath(plot_dir)
    # Ensure domain matches:
    z_les = Array(TC.get_nc_data(ds_pycles, "zc"))
    z_tcc_c = Array(TC.get_nc_data(ds_tc, "zc"))
    z_tcc_f = Array(TC.get_nc_data(ds_tc, "zf"))
    z_tcm_c = Array(TC.get_nc_data(ds_tc_main, "zc"))
    z_tcm_f = Array(TC.get_nc_data(ds_tc_main, "zf"))
    z_scm = Array(TC.get_nc_data(ds_scampy, "zc"))
    n_grid_points = length(z_tcc_c)
    @info "z extrema (les,scm,tcm,tcc): $(extrema(z_les)), $(extrema(z_scm)), $(extrema(z_tcm_c)), $(extrema(z_tcc_c))"
    @info "n-grid points (les,scm,tcm,tcc): $(length(z_les)), $(length(z_scm)), $(length(z_tcm_c)), $(length(z_tcc_c))"

    @info "time extrema (les,scm,tcm,tcc): $(extrema(time_les)), $(extrema(time_scm)), $(extrema(time_tcm)), $(extrema(time_tcc))"
    @info "n-time points (les,scm,tcm,tcc): $(length(time_les)), $(length(time_scm)), $(length(time_tcm)), $(length(time_tcc))"

    # Find the nearest matching final time:
    t_cmp = min(time_tcc[end], time_tcm[end], time_les[end], time_scm[end], t_stop)
    @info "time compared: $t_cmp"

    if length(time_tcc) == 1
        msg = "IO was collected for a single point in time."
        msg *= "This script (compute_mse) requires comparing "
        msg *= "space-time fields with multiple timesteps."
        error(msg)
    end

    # Accidentally running a short simulation
    # could improve MSE. So, let's test that
    # we run for at least t_compare. We should
    # increase this as we can reach higher CFL.
    # @test t_cmp >= t_compare

    # Ensure z_tcc and fields are consistent lengths:
    @test length(z_tcc_c) == n_grid_points
    tcc_variables = []
    computed_mse = []
    table_best_mse = []
    mse_reductions = []
    data_scales_scm = []
    data_scales_les = []
    data_scales_tcc = []
    data_scales_tcm = []
    pycles_weight = []
    plots_dict = Dict()
    plots_dict["profiles"] = Dict()
    plots_dict["contours"] = Dict()
    plot_attr = Dict()
    plot_attr[true] = Dict()
    plot_attr[false] = Dict()

    local fig_height

    for tc_var in keys(best_mse)
        data_les_arr = TC.get_nc_data(ds_pycles, tc_var)
        data_tcm_arr = TC.get_nc_data(ds_tc_main, tc_var)
        data_tcc_arr = TC.get_nc_data(ds_tc, tc_var)
        data_scm_arr = TC.get_nc_data(ds_scampy, tc_var)
        # Only compare fields that exist in the nc files
        missing_les_var = data_les_arr == nothing
        missing_tcm_var = data_tcm_arr == nothing
        missing_tcc_var = data_tcc_arr == nothing
        missing_scm_var = data_scm_arr == nothing

        coord_name = first(NC.dimnames(data_tcc_arr))
        coord_name == "zf" || coord_name == "zc" || error("Bad coord_name")

        if coord_name == "zf"
            z_tcc = z_tcc_f
        elseif coord_name == "zc"
            z_tcc = z_tcc_c
        end
        if have_tc_main
            if coord_name == "zf"
                z_tcm = z_tcm_f
            elseif coord_name == "zc"
                z_tcm = z_tcm_c
            end
        else
            z_tcm = z_tcc
        end
        if !have_scampy_ds
            z_scm = z_tcc
        end
        if !have_pycles_ds
            z_les = z_tcc
        end

        if missing_les_var
            @warn "les missing data for variable $tc_var, filling with zeros"
            data_les_arr = zeros(length(z_les), length(time_les))
            warn_msg_les = " Warning: missing data"
        else
            warn_msg_les = ""
        end
        if missing_tcm_var
            @warn "tcm missing data for variable $tc_var, filling with zeros"
            data_tcm_arr = zeros(length(z_tcm), length(time_tcm))
            warn_msg_tcm = " Warning: missing data"
        else
            warn_msg_tcm = ""
        end
        if missing_tcc_var
            @warn "tcc missing data for variable $tc_var, filling with zeros"
            data_tcc_arr = zeros(length(z_tcc), length(time_tcc))
            warn_msg_tcc = " Warning: missing data"
        else
            warn_msg_tcc = ""
        end
        if missing_scm_var
            @warn "scm missing data for variable $tc_var, filling with zeros"
            data_scm_arr = zeros(length(z_scm), length(time_scm))
            warn_msg_scm = " Warning: missing data"
        else
            warn_msg_scm = ""
        end

        # TC and TC main are the same locally
        if haskey(ENV, "BUILDKITE_COMMIT")
            new_variable = all((missing_les_var, missing_tcm_var, missing_scm_var))
        else
            new_variable = all((missing_les_var, missing_scm_var))
        end

        data_les_arr = Array(data_les_arr)'
        data_tcm_arr = Array(data_tcm_arr)'
        data_tcc_arr = Array(data_tcc_arr)'
        data_scm_arr = Array(data_scm_arr)'
        push!(tcc_variables, tc_var)

        @info "Assembling plots for $tc_var"

        @debug "     Data sizes (les,scm,tcc,tcm): $(size(data_les_arr)), $(size(data_scm_arr)), $(size(data_tcc_arr)), $(size(data_tcm_arr))"
        # Scale the data for comparison
        push!(pycles_weight, "1")

        # Interpolate data
        steady_data = size(data_les_arr, 1) == 1
        if steady_data
            data_les_cont = Dierckx.Spline1D(z_les, vec(data_les_arr); k = 1)
            data_tcc_cont = Dierckx.Spline1D(z_tcc, vec(data_tcc_arr); k = 1)
            data_tcm_cont = Dierckx.Spline1D(z_tcm, vec(data_tcm_arr); k = 1)
            data_scm_cont = Dierckx.Spline1D(z_scm, vec(data_scm_arr); k = 1)
            data_les_cont_mapped = map(z -> data_les_cont(z), z_tcc)
            data_tcm_cont_mapped = map(z -> data_tcm_cont(z), z_tcc)
            data_tcc_cont_mapped = map(z -> data_tcc_cont(z), z_tcc)
            data_scm_cont_mapped = map(z -> data_scm_cont(z), z_tcc)
        else # unsteady data
            data_les_cont = Dierckx.Spline2D(time_les, z_les, data_les_arr; kx = 1, ky = 1)
            data_tcm_cont = Dierckx.Spline2D(time_tcm, z_tcm, data_tcm_arr; kx = 1, ky = 1)
            data_tcc_cont = Dierckx.Spline2D(time_tcc, z_tcc, data_tcc_arr; kx = 1, ky = 1)
            data_scm_cont = Dierckx.Spline2D(time_scm, z_scm, data_scm_arr; kx = 1, ky = 1)
            R = range(t_start, t_cmp; length = 50)
            data_les_cont_mapped = map(z_tcc) do z
                StatsBase.mean(map(t -> data_les_cont(t, z), R))
            end
            data_tcm_cont_mapped = map(z_tcc) do z
                StatsBase.mean(map(t -> data_tcm_cont(t, z), R))
            end
            data_tcc_cont_mapped = map(z_tcc) do z
                StatsBase.mean(map(t -> data_tcc_cont(t, z), R))
            end
            data_scm_cont_mapped = map(z_tcc) do z
                StatsBase.mean(map(t -> data_scm_cont(t, z), R))
            end

        end

        # Plot comparison
        if plot_comparison
            p = Plots.plot()
            if have_pycles_ds
                Plots.plot!(
                    data_les_cont_mapped,
                    z_tcc ./ 10^3,
                    title = tc_var,
                    ylabel = "z [km]",
                    label = "PyCLES$warn_msg_les",
                )
            end
            Plots.plot!(
                data_scm_cont_mapped,
                z_tcc ./ 10^3,
                title = tc_var,
                ylabel = "z [km]",
                label = "SCAMPy$warn_msg_scm",
            )
            if have_tc_main
                Plots.plot!(
                    data_tcm_cont_mapped,
                    z_tcc ./ 10^3,
                    title = tc_var,
                    ylabel = "z [km]",
                    label = "TC.jl (main)$warn_msg_tcm",
                )
            end
            plots_dict["profiles"][tc_var] = Plots.plot!(
                data_tcc_cont_mapped,
                z_tcc ./ 10^3,
                title = tc_var,
                ylabel = "z [km]",
                label = "TC.jl$warn_msg_tcc",
            )

            clims_min = min(minimum(data_scm_arr), minimum(data_tcc_arr))
            clims_max = max(maximum(data_scm_arr), maximum(data_tcc_arr))
            clims = (clims_min, clims_max)

            width_to_height_ratio = have_tc_main ? 15 / 10 : 15 / 10
            fig_height = 1500

            plot_attr[true]["contour_1_kwargs"] = (
                bottom_margin = 0 * Plots.PlotMeasures.px,
                left_margin = 40 * Plots.PlotMeasures.px,
                right_margin = 0 * Plots.PlotMeasures.px,
                top_margin = 0 * Plots.PlotMeasures.px,
                xticks = false,
                colorbar = false,
                size = (width_to_height_ratio * fig_height, fig_height),
                ylabel = "height (km)",
            )
            plot_attr[true]["contour_2_kwargs"] = (
                bottom_margin = 0 * Plots.PlotMeasures.px,
                left_margin = 0 * Plots.PlotMeasures.px,
                right_margin = 0 * Plots.PlotMeasures.px,
                top_margin = 0 * Plots.PlotMeasures.px,
                xticks = false,
                yticks = false,
                colorbar = false,
                size = (width_to_height_ratio * fig_height, fig_height),
            )
            plot_attr[true]["contour_3_kwargs"] = (
                bottom_margin = 0 * Plots.PlotMeasures.px,
                left_margin = 0 * Plots.PlotMeasures.px,
                right_margin = -20 * Plots.PlotMeasures.px,
                top_margin = 0 * Plots.PlotMeasures.px,
                xticks = false,
                yticks = false,
                colorbar = true,
                size = (width_to_height_ratio * fig_height, fig_height), # extra space for colorbar
            )
            plot_attr[false]["contour_1_kwargs"] = (ylabel = "height (km)",)
            plot_attr[false]["contour_2_kwargs"] = (ylabel = "height (km)",)
            plot_attr[false]["contour_3_kwargs"] = (ylabel = "height (km)",)

            clim_kwarg = clims_min == clims_max ? () : (; clims = clims)

            p1 = Plots.contourf(
                time_scm ./ 3600,
                z_scm ./ 10^3,
                data_scm_arr';
                c = :viridis,
                clim_kwarg...,
                title = "$tc_var (SCAMPy)$warn_msg_scm",
                plot_attr[group_figs]["contour_1_kwargs"]...,
            )
            if have_tc_main
                p2 = Plots.contourf(
                    time_tcm ./ 3600,
                    z_tcm ./ 10^3,
                    data_tcm_arr';
                    c = :viridis,
                    clim_kwarg...,
                    title = "$tc_var (TC.jl main)$warn_msg_tcm",
                    plot_attr[group_figs]["contour_2_kwargs"]...,
                )
            end
            p3 = Plots.contourf(
                time_tcc ./ 3600,
                z_tcc ./ 10^3,
                data_tcc_arr';
                c = :viridis,
                clim_kwarg...,
                title = have_tc_main ? "$tc_var (TC.jl PR)$warn_msg_tcc" : "$tc_var (TC.jl)$warn_msg_tcc",
                plot_attr[group_figs]["contour_3_kwargs"]...,
            )
            if group_figs
                if have_tc_main
                    widths = [0.29, 0.31]
                    push!(widths, 1 - sum(widths))
                    plots_dict["contours"][tc_var] = Plots.plot(p1, p2, p3; layout = Plots.grid(1, 3, widths = widths))
                else
                    widths = [0.45]
                    push!(widths, 1 - sum(widths))
                    plots_dict["contours"][tc_var] = Plots.plot(p1, p3; layout = Plots.grid(1, 2, widths = widths))
                end
            else
                if have_tc_main
                    plots_dict["contours"][tc_var] = Plots.plot(p1, p2, p3; layout = (3, 1))
                else
                    plots_dict["contours"][tc_var] = Plots.plot(p1, p3; layout = (2, 1))
                end
            end
        end

        # Compute data scale
        data_scale_tcm = sum(abs.(data_tcm_arr)) / length(data_tcm_arr)
        data_scale_tcc = sum(abs.(data_tcc_arr)) / length(data_tcc_arr)
        data_scale_scm = sum(abs.(data_scm_arr)) / length(data_scm_arr)
        data_scale_les = sum(abs.(data_les_arr)) / length(data_les_arr)
        push!(data_scales_tcm, data_scale_tcm) # TODO: should we add this to the mse table?
        push!(data_scales_tcc, data_scale_tcc)
        push!(data_scales_scm, data_scale_scm)
        push!(data_scales_les, data_scale_les) # TODO: should we add this to the mse table?

        # Compute mean squared error (mse)
        if have_pycles_ds # LES takes first precedence
            mse_single_var = sum((data_les_cont_mapped .- data_tcc_cont_mapped) .^ 2)
            data_scale_used = data_scale_les
            scaled_mse = mse_single_var / data_scale_used^2 # Normalize by data scale
        elseif have_tc_main # TC.jl main takes second precedence
            mse_single_var = sum((data_tcm_cont_mapped .- data_tcc_cont_mapped) .^ 2)
            data_scale_used = data_scale_tcm
            scaled_mse = mse_single_var / data_scale_used^2 # Normalize by data scale
        elseif have_scampy_ds # SCAMPy takes third precedence
            mse_single_var = sum((data_scm_cont_mapped .- data_tcc_cont_mapped) .^ 2)
            data_scale_used = data_scale_scm
            scaled_mse = mse_single_var / data_scale_used^2 # Normalize by data scale
        else
            error("No dataset to compute MSE")
        end

        haskey(best_mse, tc_var) || error("No key found in best_mse for variable $tc_var")

        if new_variable || best_mse[tc_var] isa String
            mse[tc_var] = "NA"
            push!(mse_reductions, "NA")
            push!(table_best_mse, "NA")
            push!(computed_mse, "NA")
        else
            new_variable && @warn "$tc_var must be added to the MSE tables."
            mse[tc_var] = scaled_mse
            push!(mse_reductions, (best_mse[tc_var] - mse[tc_var]) / best_mse[tc_var] * 100)
            push!(table_best_mse, best_mse[tc_var])
            push!(computed_mse, mse[tc_var])
        end
    end

    save_plots(plot_dir, plots_dict; group_figs, have_tc_main, fig_height, case_name)

    # Tabulate output
    header = (
        ["Variable", "Weight", "Data scale", "Data scale", "MSE", "MSE", "MSE"],
        ["TC.jl (EDMF)", "PyCLES", "tcc", "scm", "Computed", "Best", "Reduction (%)"],
    )
    table_data = hcat(
        tcc_variables,
        pycles_weight,
        data_scales_tcc,
        data_scales_scm,
        computed_mse,
        table_best_mse,
        mse_reductions,
    )

    @info Printf.@sprintf("case_name comparison: %s at time t=%s\n", case_name, t_cmp)
    hl_worsened_mse = PrettyTables.Highlighter(
        (data, i, j) -> !sufficient_mse(data[i, 5], data[i, 6]) && j == 5,
        PrettyTables.crayon"red bold",
    )
    hl_worsened_mse_reduction = PrettyTables.Highlighter(
        (data, i, j) -> !sufficient_mse(data[i, 5], data[i, 6]) && j == 7,
        PrettyTables.crayon"red bold",
    )
    hl_improved_mse = PrettyTables.Highlighter(
        (data, i, j) -> sufficient_mse(data[i, 5], data[i, 6]) && j == 7,
        PrettyTables.crayon"green bold",
    )
    PrettyTables.pretty_table(
        table_data;
        header,
        formatters = PrettyTables.ft_printf("%.16e", 5:6),
        header_crayon = PrettyTables.crayon"yellow bold",
        subheader_crayon = PrettyTables.crayon"green bold",
        highlighters = (hl_worsened_mse, hl_improved_mse, hl_worsened_mse_reduction),
        crop = :none,
    )

    return mse
end

function save_plots(plot_dir, plots_dict; group_figs = true, have_tc_main, fig_height, case_name)
    vars_to_skip = [
        "thetal_mean",
        "updraft_thetal",
        "u_mean",
        "v_mean",
        "qt_mean",
        "updraft_qt",
        "temperature_mean",
        "total_flux_h",
        "total_flux_qt",
        "massflux",
    ]
    all_contours = plots_dict["contours"]
    filter!(p -> !occursin("var", p.first), all_contours) # skip higher order moment contours
    filter!(p -> !any(map(x -> occursin(x, p.first), vars_to_skip)), all_contours)
    if group_figs
        all_contours = collect(values(all_contours))
        Plots.plot!(all_contours[end], xticks = true, xlabel = "Time (hr)", bottom_margin = 20 * Plots.PlotMeasures.px)

        n_plots = length(all_contours)
        n_cols = 1
        n_rows = ceil(Int, n_plots / n_cols)
        @info "     Saving $(joinpath(plot_dir, "contours.png"))"
        title = Plots.plot(
            title = case_name,
            grid = false,
            showaxis = false,
            xticks = false,
            yticks = false,
            bottom_margin = -20 * Plots.PlotMeasures.px,
        )
        layout = Plots.@layout [a{0.01h}; Plots.grid(n_rows, n_cols)]

        Plots.plot(
            title,
            all_contours...;
            layout = layout,
            framestyle = :box,
            margin = 20 * Plots.PlotMeasures.px,
            dpi = 200,
        )
        Plots.savefig(joinpath(plot_dir, "contours.png"))

        width_to_height_ratio = 15 / 10
        all_profiles = collect(values(plots_dict["profiles"]))
        n_plots = length(all_profiles)
        n_cols = 5
        n_rows = ceil(Int, n_plots / n_cols)
        @info "     Saving $(joinpath(plot_dir, "profiles.png"))"
        n_empty_plots = n_rows * n_cols - n_plots
        all_profiles =
            [all_profiles..., ntuple(i -> Plots.plot(; framestyle = :none, xticks = false), n_empty_plots)...]
        for (k, I) in enumerate(CartesianIndices((1:n_cols, 1:n_rows)))
            k > length(all_profiles) && continue
            j, i = (Tuple(I)...,)
            if j == 1 # left plots
                Plots.plot!(
                    all_profiles[k],
                    left_margin = 50 * Plots.PlotMeasures.px,
                    bottom_margin = -2 * Plots.PlotMeasures.px,
                    right_margin = -2 * Plots.PlotMeasures.px,
                    top_margin = -2 * Plots.PlotMeasures.px,
                )
            elseif j == n_cols # right plots
                Plots.plot!(
                    all_profiles[k],
                    yticks = false,
                    ylabel = "",
                    bottom_margin = -2 * Plots.PlotMeasures.px,
                    left_margin = -2 * Plots.PlotMeasures.px,
                    right_margin = -2 * Plots.PlotMeasures.px,
                    top_margin = -2 * Plots.PlotMeasures.px,
                )
            else # middle plots
                Plots.plot!(
                    all_profiles[k],
                    yticks = false,
                    ylabel = "",
                    bottom_margin = -2 * Plots.PlotMeasures.px,
                    left_margin = -2 * Plots.PlotMeasures.px,
                    right_margin = -2 * Plots.PlotMeasures.px,
                    top_margin = -2 * Plots.PlotMeasures.px,
                )
            end
        end
        title = Plots.plot(
            title = case_name,
            grid = false,
            showaxis = false,
            xticks = false,
            yticks = false,
            bottom_margin = -20 * Plots.PlotMeasures.px,
        )
        layout = Plots.@layout [a{0.01h}; Plots.grid(n_rows, n_cols)]
        Plots.plot(
            title,
            all_profiles...;
            layout = layout,
            size = (width_to_height_ratio * fig_height, fig_height),
            framestyle = :box,
            margin = 20 * Plots.PlotMeasures.px,
            dpi = 200,
        )
        Plots.savefig(joinpath(plot_dir, "profiles.png"))
    else
        for tc_var in keys(all_contours)
            @info "     Saving $(joinpath(plot_dir, "contours_$tc_var.png"))"
            Plots.plot(all_contours[tc_var])
            Plots.savefig(joinpath(plot_dir, "contours_$tc_var.png"))
        end

        all_profiles = plots_dict["profiles"]
        for tc_var in keys(all_profiles)
            @info "     Saving $(joinpath(plot_dir, "profiles_$tc_var.png"))"
            Plots.plot(all_profiles[tc_var])
            Plots.savefig(joinpath(plot_dir, "profiles_$tc_var.png"))
        end
    end
end

sufficient_mse(computed_mse::Real, best_mse::Real) = computed_mse <= best_mse + sqrt(eps())
# New variables report string MSEs (NA), so those should always be sufficient.
sufficient_mse(computed_mse::String, best_mse::String) = true
sufficient_mse(computed_mse::Real, best_mse::String) = true
sufficient_mse(computed_mse::String, best_mse::Real) = true

function test_mse(computed_mse, best_mse, key)
    mse_not_regressed = sufficient_mse(computed_mse[key], best_mse[key])
    mse_not_regressed || @show key
    @test mse_not_regressed
end


text:
export Func, default_funcs, AND, OR, XOR, NAND, NOR, NOT, ZERO, ONE

type Func
    func::Function
    arity::Integer
    name::AbstractString
end

Func(f::Function, a::Integer) = Func(f, a, string(f))

const AND = Func(&, 2)
const OR = Func(|, 2)
const XOR = Func($, 2)
const NAND = Func((x,y)->~(x & y),2,"^&")
const NOR = Func((x,y)->~(x | y),2,"^|")
const NOT = Func(~, 1)
const ZERO = Func(() -> convert(BitString, 0), 0, "0")
const ONE = Func(() -> ~convert(BitString, 0), 0, "1")

function default_funcs()
    return [AND, OR, XOR, NOT, ZERO, ONE]
end


text:
struct HsAtom end
struct HsContext end
struct HsJIT end
struct NativeFunctionObj end
struct NativeFunctionSignature
    arg::Cstring
    res::Cstring
    _ccall::Cstring  # can't name this field `ccall` as that is a built-in in julia
end

struct CRectArray
    data::Ptr{Nothing}
    shape_ptr::Ptr{Int64}
    strides_ptr::Ptr{Int64}
end

"""
    TaggedUnion{Tuple{A, B, ...}}

Represents a tagged union over types `A`, `B` etc.
Must have a first field `tag::UInt64`
and a second field `payload` which must be some isbits type that can hold the largest
element of the union (Probably a `NTuple{UInt8}` or similar).
Julia doens't directly support Unions in the mappting to-from C
so we store the data as arbitary bits then force it to be reinterpretted based on the tag
"""
abstract type TaggedUnion{T<:Tuple} end
function bust_union(x::TaggedUnion{U}) where U
    T = U.parameters[Int(x.tag) + 1]
    return bust_union(force_reinterpret(T, x.payload))
end
bust_union(x) = x  # not a union, leave it alone

function force_reinterpret(::Type{T}, raw) where T
    isbits(raw) || throw(ArgumentError("Can only reinterpret from a isbits type"))
    isbitstype(T) || throw(ArgumentError("Can only reinterpret into a isbits type"))
    return unsafe_load(Ptr{T}(Base.pointer_from_objref(Ref(raw))))
end

struct CLit <: TaggedUnion{Tuple{Int64, Int32, Int8, Float64, Float32}}
    tag::UInt64
    payload::NTuple{8, UInt8}  # actually the Union, needs to be big enough to hold largest which is an Int64
end

struct CAtom <: TaggedUnion{Tuple{CLit, CRectArray}}
    tag::UInt64
    payload::NTuple{3, UInt64}  # actually the Union, needs to be big enough to hold largest which is the CRectArray
end
function CAtom(atm::Ptr{HsAtom})
    result = Ref{CAtom}()
    success = to_CAtom(atm, result)
    iszero(success) && throw_from_dex()
    return result[]
end

"""
    juliaize(x)

Get the corresponding Julia object from some output of Dex.
"""
juliaize(x::CAtom) = bust_union(x)
juliaize(x::Ptr{HsAtom}) = juliaize(CAtom(x))

max_stars_repo_name: AlexDikelsky/dex-lang

max_issues_repo_name: AlexDikelsky/dex-lang

max_forks_repo_name: AlexDikelsky/dex-lang


text:
mutable struct BaseEnumeration{T} <: AbstractStandard{T}
    A::Array{T, 2}
    b::Array{T, 1}
    c::Array{T, 1}
    xstar::Array{T, 1}
    vstar::T
    status::AbstractStatus
    function BaseEnumeration(A::Array{T, 2}, b::Array{T, 1}, c::Array{T, 1}) where T
        m,n = size(A)
        @assert !(typeof(T) <: Integer)  "Type $T cannot be a subtype of Integer"
        @assert length(b) == m "dimension of A and b mismatch, size(A) = ($m, $n), length(b) = $(length(b)) != $m"
        @assert length(c) == n "dimension of A and c mismatch, size(A) = ($m, $n), length(c) = $(length(c)) != $n"
        @assert rank(A) == m "A is not a full rank Matrix"
        be = new{T}(A, b, c, [],  T(Inf), Unknown())
        return be
    end
end

function solve!(be::BaseEnumeration{T}; verbose::Bool = false) where T
    m, n = size(be.A)
    vopt = T(Inf)
    sol = zeros(T, n)
    allComb = [c for c in combinations(1:n, m)]
    for comb in allComb
        verbose && println(comb)
        B = be.A[:, comb]
        x_B = B\be.b
        v_B = dot(be.c[comb], x_B)
        if rank(B) == m && all(x_B .>= -10*eps()) && (v_B < vopt)
            vopt = v_B
            sol[:] .= zero(T)
            sol[comb] = x_B
        end
    end
    if vopt ==  T(Inf)
        be.status = Infeasible()
    else
        be.status = Optimal()
        be.vstar = vopt
        be.xstar = sol
    end
    be.status
end


text:
module Chemostat

import MetNets
import MetNets: lb!, ub!, rxns, MetNet
import MetEP
import MetLP

include("SteadyState/SteadyState.jl")

end # module

max_stars_repo_name: josePereiro/Chemostat

max_issues_repo_name: josePereiro/Chemostat

max_forks_repo_name: josePereiro/Chemostat