text:
using Base: UInt8
using faust_jll

include("uiglue.jl")

export llvm_dsp, llvm_dsp_factory,
    createCDSPFactoryFromString,
    writeCDSPFactoryToIR,
    createCDSPInstance,
    getNumInputsCDSPInstance,
    getNumOutputsCDSPInstance,
    buildUserInterfaceCDSPInstance,
    getSampleRateCDSPInstance,
    initCDSPInstance,
    instanceInitCDSPInstance,
    instanceConstantsCDSPInstance,
    instanceResetUserInterfaceCDSPInstance,
    instanceClearCDSPInstance,
    cloneCDSPInstance,
    metadataCDSPInstance,
    computeCDSPInstance,
    deleteCDSPFactory,
    deleteCDSPInstance,
    freeCMemory

mutable struct llvm_dsp_factory
end

mutable struct llvm_dsp
end

# /**
# * Create a Faust DSP factory from a DSP source code as a string. Note that the library keeps an internal cache of all
# * allocated factories so that the compilation of the same DSP code (that is same source code and
# * same set of 'normalized' compilations options) will return the same (reference counted) factory pointer. You will
# * have to explicitly use deleteCDSPFactory to properly decrement reference counter when the factory is no more needed.
# *
# * @param name_app - the name of the Faust program
# * @param dsp_content - the Faust program as a string
# * @param argc - the number of parameters in argv array
# * @param argv - the array of parameters (Warning : aux files generation options will be filtered (-svg, ...) --> use generateAuxFiles)
# * @param target - the LLVM machine target: like 'i386-apple-macosx10.6.0:opteron',
# *                 using an empty string takes the current machine settings,
# *                 and i386-apple-macosx10.6.0:generic kind of syntax for a generic processor
# * @param error_msg - the error string to be filled, has to be 4096 characters long
# * @param opt_level - LLVM IR to IR optimization level (from -1 to 4, -1 means 'maximum possible value'
# * since the maximum value may change with new LLVM versions)
# *
# * @return a valid DSP factory on success, otherwise a null pointer.
# */
# llvm_dsp_factory* createCDSPFactoryFromString(const char* name_app,
#                                              const char* dsp_content,
#                                              int argc, const char* argv[],
#                                              const char* target,
#                                              char* error_msg,
#                                              int opt_level);
function createCDSPFactoryFromString(
    name_app,
    dsp_content,
    argv,
    target,
    opt_level)
    argc = length(argv)
    error_msg = Vector{UInt8}(undef, 4096)
    output_ptr = cd(() -> ccall(
        (:createCDSPFactoryFromString, libfaust),
        Ptr{llvm_dsp_factory},
        (Cstring, Cstring, Cint, Ptr{Cstring}, Cstring, Ptr{UInt8}, Cint),
        name_app, dsp_content, argc, argv, target, error_msg, opt_level,
    ), dirname(stdfaust_lib))
    if output_ptr == C_NULL
        error = GC.@preserve error_msg unsafe_string(pointer(error_msg))
        throw(ErrorException("Could not initialize C DSP factory: $error"))
    end

    return output_ptr
end

# /**
# * Delete a Faust DSP factory, that is decrements it's reference counter, possibly really deleting the internal pointer.
# * Possibly also delete DSP pointers associated with this factory, if they were not explicitly deleted with deleteCDSPInstance.
# * Beware : all kept factories and DSP pointers (in local variables...) thus become invalid.
# *
# * @param factory - the DSP factory to be deleted.
# *
# * @return true if the factory internal pointer was really deleted, and false if only 'decremented'.
# */
# bool deleteCDSPFactory(llvm_dsp_factory* factory);
function deleteCDSPFactory(factory)
    ret = ccall(
        (:deleteCDSPFactory, libfaust),
        Cuchar,
        (Ptr{llvm_dsp_factory},),
        factory
    )
    return Bool(ret)
end

function startMTDSPFactories()
    return Bool(ccall(
        (:startMTDSPFactories, libfaust),
        Cuchar,
        (),
    ))
end

function stopMTDSPFactories()
    ccall(
        (:stopMTDSPFactories, libfaust),
        Cuchar,
        (),
    )
end

# /**
# * Write a Faust DSP factory into a LLVM IR (textual) string.
# *
# * @param factory - the DSP factory
# *
# * @return the LLVM IR (textual) as a string (to be deleted by the caller using freeCMemory).
# */
# char* writeCDSPFactoryToIR(llvm_dsp_factory* factory);
function writeCDSPFactoryToIR(factory)
    output_str = ccall(
        (:writeCDSPFactoryToIR, libfaust),
        Cstring,
        (Ptr{llvm_dsp_factory},),
        factory
    )
    if output_str == C_NULL
        throw(ErrorException("Could not write DSP factory IR"))
    end
    copied = GC.@preserve output_str unsafe_string(output_str)
    freeCMemory(output_str)
    return copied
end

# int getNumInputsCDSPInstance(llvm_dsp* dsp);
function getNumInputsCDSPInstance(dsp)
    ret = ccall(
        (:getNumInputsCDSPInstance, libfaust),
        Cint,
        (Ptr{llvm_dsp},),
        dsp
    )
    return ret
end

# int getNumOutputsCDSPInstance(llvm_dsp* dsp);
function getNumOutputsCDSPInstance(dsp)
    return ccall(
        (:getNumOutputsCDSPInstance, libfaust),
        Cint,
        (Ptr{llvm_dsp},),
        dsp
    )
end

# void buildUserInterfaceCDSPInstance(llvm_dsp* dsp, UIGlue* interface);
function buildUserInterfaceCDSPInstance(dsp, interface)
    return ccall(
        (:buildUserInterfaceCDSPInstance, libfaust),
        Cvoid,
        (Ptr{llvm_dsp}, Ref{UIGlue}),
        dsp, interface
    )
end

# /**
#  * The free function to be used on memory returned by getCDSPMachineTarget, getCName, getCSHAKey,
#  * getCDSPCode, getCLibraryList, getAllCDSPFactories, writeCDSPFactoryToBitcode,
#  * writeCDSPFactoryToIR, writeCDSPFactoryToMachine,expandCDSPFromString and expandCDSPFromFile.
#  *
#  * This is MANDATORY on Windows when otherwise all nasty runtime version related crashes can occur.
#  *
#  * @param ptr - the pointer to be deleted.
#  */
function freeCMemory(p)
    ccall(
        (:freeCMemory, libfaust),
        Cvoid,
        (Cstring,),
        p
    )
end

# int getSampleRateCDSPInstance(llvm_dsp* dsp);
function getSampleRateCDSPInstance(dsp)
    return ccall(
        (:getSampleRateCDSPInstance, libfaust),
        Cint,
        (Ptr{llvm_dsp},),
        dsp
    )
end

# void initCDSPInstance(llvm_dsp* dsp, int sample_rate);
function initCDSPInstance(dsp, sample_rate)
    ccall(
        (:initCDSPInstance, libfaust),
        Cvoid,
        (Ptr{llvm_dsp}, Cint),
        dsp, sample_rate
    )
end

# void instanceInitCDSPInstance(llvm_dsp* dsp, int sample_rate);
function instanceInitCDSPInstance(dsp, sample_rate)
    ccall(
        (:instanceInitCDSPInstance, libfaust),
        Cvoid,
        (Ptr{llvm_dsp}, Cint),
        dsp, sample_rate
    )
end

# void instanceConstantsCDSPInstance(llvm_dsp* dsp, int sample_rate);
function instanceConstantsCDSPInstance(dsp, sample_rate)
    ccall(
        (:instanceConstantsCDSPInstance, libfaust),
        Cvoid,
        (Ptr{llvm_dsp}, Cint),
        dsp, sample_rate
    )
end

# void instanceResetUserInterfaceCDSPInstance(llvm_dsp* dsp);
function instanceResetUserInterfaceCDSPInstance(dsp)
    ccall(
        (:instanceResetUserInterfaceCDSPInstance, libfaust),
        Cvoid,
        (Ptr{llvm_dsp},),
        dsp
    )
end

# void instanceClearCDSPInstance(llvm_dsp* dsp);
function instanceClearCDSPInstance(dsp)
    ccall(
        (:instanceClearCDSPInstance, libfaust),
        Cvoid,
        (Ptr{llvm_dsp},),
        dsp
    )
end

# llvm_dsp* cloneCDSPInstance(llvm_dsp* dsp);
function cloneCDSPInstance(dsp)
    return ccall(
        (:cloneCDSPInstance, libfaust),
        Ptr{llvm_dsp},
        (Ptr{llvm_dsp},),
        dsp
    )
end

mutable struct MetaGlue
end

# void metadataCDSPInstance(llvm_dsp* dsp, MetaGlue* meta);
function metadataCDSPInstance(dsp, meta)
    ccall(
        (:metadataCDSPInstance, libfaust),
        Cvoid,
        (Ptr{llvm_dsp}, Ptr{MetaGlue}),
        dsp, meta,
    )
end

# void computeCDSPInstance(llvm_dsp* dsp, int count, FAUSTFLOAT** input, FAUSTFLOAT** output);
function computeCDSPInstance(dsp, count, input = nothing, output = nothing)
    if isnothing(input)
        inputChannels = getNumInputsCDSPInstance(dsp)
        input = zeros(Float32, count, inputChannels)
    end
    if isnothing(output)
        outputChannels = getNumOutputsCDSPInstance(dsp)
        output = zeros(Float32, count, outputChannels)
    end
    inputRef = [pointer(input, i) for i=1:size(input, 1):length(input)]
    outputRef = [pointer(output, i) for i=1:size(output, 1):length(output)]

    ccall(
        (:computeCDSPInstance, libfaust),
        Cvoid,
        (Ptr{llvm_dsp}, Cint, Ptr{Ptr{Float32}}, Ptr{Ptr{Float32}}),
        dsp, count, inputRef, outputRef
    )
    return output
end

# /**
# * Create a Faust DSP instance.
# *
# * @param factory - the Faust DSP factory
# *
# * @return the DSP instance on success, otherwise a null pointer.
# */
# llvm_dsp* createCDSPInstance(llvm_dsp_factory* factory);
function createCDSPInstance(factory)
    output_ptr = ccall(
        (:createCDSPInstance, libfaust),
        Ptr{llvm_dsp},
        (Ptr{llvm_dsp_factory},),
        factory
    )
    if output_ptr == C_NULL
        throw(ErrorException("Could not initialize C DSP instance"))
    end
    return output_ptr
end

# /**
#  * Delete a Faust DSP instance.
#  *
#  * @param dsp - the DSP instance to be deleted.
#  */
# void deleteCDSPInstance(llvm_dsp* dsp);
function deleteCDSPInstance(dsp)
    ccall(
        (:deleteCDSPInstance, libfaust),
        Cvoid,
        (Ptr{llvm_dsp},),
        dsp
    )
end


text:
######### StanVariational Bernoulli example  ###########

using StanVariational

bernoulli_model = "
data {
  int<lower=1> N;
  int<lower=0,upper=1> y[N];
}
parameters {
  real<lower=0,upper=1> theta;
}
model {
  theta ~ beta(1,1);
  y ~ bernoulli(theta);
}
";

bernoulli_data = Dict("N" => 10, "y" => [0, 1, 0, 1, 0, 0, 0, 0, 0, 1])

# Keep tmpdir across multiple runs to prevent re-compilation
tmpdir = joinpath(@__DIR__, "tmp")

stanmodel = VariationalModel(
  "bernoulli", bernoulli_model; tmpdir = tmpdir)

(sample_file, log_file) = stan_variational(stanmodel; data=bernoulli_data)

if sample_file !== Nothing

  (chns, cnames) = read_variational(stanmodel)

  # Show the same output in DataFrame format
  sdf = read_summary(stanmodel)
  display(sdf)
  println()

  # Retrieve mean value of theta from the summary
  sdf[:theta, :mean]

end

ext: jl


text:
123.4 + 234.5     # This is a comment, ignored by Julia

7 / 3             # Integer division returns floating point

7 ÷ 3             # Truncated integer division (type \div and press tab)

7 % 3             # Remainder

2^7 - 4           # All integer expression which returns an integer

3 * (4.2 - 9.7)   # Use parentheses to control the order of operations

# We will later learn why the previous result is not exact

5.6e7 ^ 0.1       # Scientific notation, means 5.6⋅10⁷ (not e⁷)

1 + 2
34 / 3

println(1 + 2.1)
println("The result is ", -123 / 1000000)
34 / 3; # The output will not be displayed

typeof(1.23)

typeof(-745)

typeof(1e20)

println(typemin(Int64))
println(typemax(Int64))
println(typemin(Int128))
println(typemax(Int128))

1000^7

Int128(1000)^7


text:
using FdeSolver
using Plots

## inputs
tSpan = [0, 10]     # [intial time, final time]
x0 = [1, 1]             # intial value
β = 2            # order of the derivative
par = [16.0, 4.0] # [spring constant for a mass on a spring, inertial mass]
h = 0.01

## Equation
function F(t, n, β, x, par)

      K = par[1]
      m = par[2]

      - K ./ m .* x[n]

end

## Numerical solution
t, Yapp = FDEsolver(F, tSpan, x0, β, nothing, par, h = h)

#plot
plot(t, Yapp, linewidth = 5, title = "Simple harmonic motion (order=2)",
     xaxis = "Time (t)", yaxis = "x(t)", label = "Approximation")
a = x0[1] .* map(cos, sqrt(par[1] / par[2]) .* t) .+ x0[2] ./ sqrt(par[1] / par[2]) .* map(sin, sqrt(par[1] / par[2]) .* t)
plot!(t, a, lw = 3, ls = :dash, label = "Exact solution")


text:
module MixedLayerThermoclineDynamics

using
  DocStringExtensions,
  OffsetArrays

export
  AbstractGrid,
  AbstractTopology,
  Periodic,
  Bounded,
  Grid1D,
  Grid2D,

  AbstractLocation,
  Face,
  Center,
  Centre,

  AbstractField,
  Field1D,
  Field2D,
  Field,
  fill_halos!,
  𝐼x!, 𝐼y!,
  ∂x!, ∂y!

import Base: show

""" Abstract supertype for topology of grids. """
abstract type AbstractTopology end

""" Type for location at the cell centres. """
struct Periodic <: AbstractTopology end

""" Type for location at the cell centres. """
struct Bounded <: AbstractTopology end

""" Abstract supertype for grids. """
abstract type AbstractGrid end

""" Abstract supertype for location of variables. """
abstract type AbstractLocation end

""" Abstract supertype for fields. """
abstract type AbstractField end

include("grids.jl")
include("fields.jl")

end # module


text:
module MyModule
function complexfunc(dims::Int, iters::Int)::Float64
    x::Array{Float64, 2} = ones(Float64, (dims, dims))
    for i in 1:iters
        x .+= 1.0
    end

    z::Float64 = 0.0
    for j in 1:dims
        z += x[j, j]
    end

    z
end

function kwfunc(dims::Int, iters::Int; kw::Float64=3.0)::Float64
    x::Array{Float64, 2} = ones(Float64, (dims, dims))
    for i in 1:iters
        x .+= 1.0
    end

    z::Float64 = 0.0
    for j in 1:dims
        z += x[j, j]
    end

    z + kw
end

function subfunc()::Float64
    throw(ErrorException("As expected"))
end

function throwingfunc()::Float64
    subfunc()
end
end

max_stars_repo_name: marhkb/jlrs

max_issues_repo_name: marhkb/jlrs

max_forks_repo_name: marhkb/jlrs


text:
module ITensorTTN

using ITensors
using ITensorInfiniteMPS
using Reexport

@reexport using ITensorInfiniteMPS

using ITensorInfiniteMPS: celltags, CelledVector, tag_starting_with

export InfiniteTTN, degree

struct InfiniteTTN
  data::CelledVector{ITensor}
end

function degree(ψ::InfiniteTTN)
  zs = order.(ψ.data[Cell(1)])
  z1 = first(zs)
  @assert all(==(z1), zs)
  # Remove one for site index
  return z1 - 1
end

childtagprefix() = "d="
childtags(n::Integer) = TagSet(string(childtagprefix(), n))
childtags(::Nothing) = ts""

generationtagprefix() = "g="
function generationtags(n1n2::Pair{<:Integer,<:Integer})
  n1, n2 = n1n2
  return TagSet(string(generationtagprefix(), n1, "↔", n2))
end
generationtags(::Nothing) = ts""

function ttntags(; c=nothing, d=nothing, g=nothing)
  return reduce(addtags, [celltags(c), childtags(d), generationtags(g)])
end

function init_space(::Type{Index{Int}}, linkdims::Integer)
  return linkdims
end

function init_space(::Type{ITensors.QNIndex}, linkdims::Integer)
  return (QN() => linkdims)
end

function InfiniteTTN(s::Vector{IndexT}, st; z::Integer, linkdims) where {IndexT<:Index}
  N = length(s)
  s∞ = infsiteinds(s)
  product_state = [state(s∞[n], st(n)) for n in 1:N]
  ttn = Vector{ITensor}(undef, N)
  link_space = init_space(IndexT, linkdims)
  l_parents = Vector{IndexT}(undef, N)
  l_children = Vector{Vector{IndexT}}(undef, N)
  n = 1
  if N == 1
    l_parents[n] = Index(link_space; tags=ttntags(; c=0, d=1, g=(mod1(n - 1, N) => n)))
    l_child = translatecell(l_parents[n], 1)
    l_children[n] = [replacetags(l_child, childtags(1) => childtags(d)) for d in 1:(z - 1)]
    ttn[n] = randomITensor(dag(l_parents[n]), l_children[n]) * product_state[n]
    return InfiniteTTN(ttn)
  end
  l_parents[n] = Index(link_space; tags=ttntags(; c=0, d=1, g=(mod1(n - 1, N) => n)))
  l_child = Index(link_space; tags=ttntags(; c=1, g=(n => mod1(n + 1, N))))
  l_children[n] = [addtags(l_child, childtags(d)) for d in 1:(z - 1)]
  ttn[n] = randomITensor(dag(l_parents[n]), l_children[n]) * product_state[n]
  for n in 2:(N - 1)
    l_parents[n] = first(l_children[n - 1])
    l_child = Index(link_space; tags=ttntags(; c=1, g=(n => mod1(n + 1, N))))
    l_children[n] = [addtags(l_child, childtags(d)) for d in 1:(z - 1)]
    ttn[n] = randomITensor(dag(l_parents[n]), l_children[n]) * product_state[n]
  end
  n = N
  l_parents[n] = first(l_children[n - 1])
  l_child = translatecell(l_parents[1], 1)
  l_children[n] = [replacetags(l_child, childtags(1) => childtags(d)) for d in 1:(z - 1)]
  ttn[n] = randomITensor(dag(l_parents[n]), l_children[n]) * product_state[n]
  return InfiniteTTN(ttn)
end

getchildtag(ts::TagSet) = tag_starting_with(ts, childtagprefix())
function getchild(ts::TagSet)
  childtag = getchildtag(ts)
  return parse(Int, childtag[(length(childtagprefix()) + 1):end])
end
getchildtags(ts::TagSet) = childtags(getchild(ts))
getchildtags(i::Index) = getchildtags(tags(i))

# Get the tensor at generation `g` and branch `b`
function Base.getindex(ψ::InfiniteTTN, g::Integer, b::Integer)
  # TODO: don't allow `g < 0`, make `g == 0` a special case
  @assert 1 ≤ b ≤ (degree(ψ) - 1)
  # Tensor at generation `g` on the first branch
  ψg1 = ψ.data[g]
  b == 1 && return ψg1
  l_parent_1 = commonind(ψ.data[g - 1], ψg1)
  l_parent_b = replacetags(l_parent_1, getchildtags(l_parent_1) => childtags(b))
  return replaceinds(ψg1, l_parent_1 => l_parent_b)
end

function ITensorInfiniteMPS.TransferMatrix(
  ψ::InfiniteTTN, g::Integer, d::Pair{<:Integer,<:Integer}; dir="in"
)
  @assert dir == "in"
  d1, d2 = d
  A = ψ[g, d2]
  l_in = commonind(A, ψ[g + 1, d1])
  l_out = commonind(A, ψ[g - 1, 1])
  Adag = prime(dag(A); inds=(l_in, l_out))
  return ITensorMap(
    [A, Adag]; input_inds=(dag(l_in), l_in'), output_inds=(l_out, dag(l_out)')
  )
end

end


text:
type AFile
    should_reload::Bool
    mtime::Float64
    deps::Vector{UTF8String}
end

type Dep
    should_reload::Bool
    name::UTF8String
end
suppress_warnings = false
const files = Dict{UTF8String,AFile}()
const options = @compat Dict{Symbol,Any}(:constants=>false, :strip_types=>true, :smart_types=>true, :verbose_level=>:warn, :state=>:on)
# verbose_level = :warn
# state = :on


text:
### FIXME: This function must be in ZipFile.jl:
function ZipFile.Writer(thunk :: Function, args...)
    io = ZipFile.Writer(args...)
    try
        thunk(io)
    finally
        close(io)
    end
end

### FIXME: This function must be in ZipFile.jl:
function ZipFile.Reader(thunk :: Function, args...)
    io = ZipFile.Reader(args...)
    try
        thunk(io)
    finally
        close(io)
    end
end

"""
    datafile(name)

Return the name of a file inside the archive with slices.
"""
datafile(name :: AbstractString) = joinpath("data", name)

"""
    slice_prefix(name)

Given the name of the archive `name`, return the name prefix for
slices.
"""
slice_prefix(name :: AbstractString) =
    splitext(splitpath(name)[end])[begin]

function find_file(root :: ZipFile.Reader, name :: AbstractString)
    files = root.files

    pos = findfirst(files) do file
        file.name == name
    end

    return files[pos]
end

"""
    stich(x, y)

Concatenate arrays `x` and `y` along the third axis
"""
stich(x :: AbstractArray{T}, y :: AbstractArray{T}) where T =
    cat(x, y; dims = 3)


text:
using Match

input = joinpath(@__DIR__, "input")
lines = readlines(input)

lines = map(line -> collect(line), lines)

function evalue(line)
    seen = []
    while length(line) > 0
        curr = popfirst!(line)
        @match curr begin
            ')' => pop!(seen) != '(' && return 3
            ']' => pop!(seen) != '[' && return 57
            '}' => pop!(seen) != '{' && return 1197
            '>' => pop!(seen) != '<' && return 25137
            otherwise => push!(seen, otherwise)
        end
    end
    return 0
end

function remvalue(line)
    seen = []
    while length(line) > 0
        curr = popfirst!(line)
        @match curr begin
            ')' => pop!(seen) != '(' && return nothing
            ']' => pop!(seen) != '[' && return nothing
            '}' => pop!(seen) != '{' && return nothing
            '>' => pop!(seen) != '<' && return nothing
            otherwise => push!(seen, otherwise)
        end
    end
    return seen
end

function score(line)
    acc = 0
    while length(line) > 0
        acc *= 5
        curr = pop!(line)
        if curr == '('
            acc += 1
        elseif curr == '['
            acc += 2
        elseif curr == '{'
            acc += 3
        elseif curr == '<'
            acc += 4
        end
    end
    return acc
end

p1 = sum(map(evalue, deepcopy(lines)))

scores = sort(map(score, filter(!isnothing, map(remvalue, deepcopy(lines)))))
midpoint = div(length(scores), 2) + 1
p2 = scores[midpoint]

println("-----------------------------------------------------------------------")
println("syntax scoring -- part one :: $p1")
println("syntax scoring -- part two :: $p2")
println("-----------------------------------------------------------------------")

@assert(p1 == 243939)
@assert(p2 == 2421222841)


text:
module KeywordCalls

using NestedTuples
using MLStyle: @match
using GeneralizedGenerated

export kwcall, @kwcall

"""
    kwcallperm(f, keys::Tuple{Symbol})

Compute the permutation required to get `keys` to the ordering declared by `@kwcall`
"""
function kwcallperm(f, keys)
    sortedkeys = Tuple(sort(collect(keys)))

    # The permutation to get from sorted order to the preferred order
    π = baseperm[(f, sortedkeys)]

    # The permutation to get from the call order to the sorted order
    σ = sortperm(collect(keys))

    # Composing the permutations
    return σ[π]
end

# `baseperm[(f, sortedargs::Tuple{Symbol})]` gives the permutation
# from sorted arguments to the ordering declared with @kwcall
const baseperm = Dict()

"""
    @kwcall f(b,a,d)

Declares that any call `f(::NamedTuple{N})` with `sort(N) == (:a,:b,:d)`
should be dispatched to the method already defined on `f(::NamedTuple{(:b,:a,:d)})`
"""
macro kwcall(call)
    esc(_kwcall(call))
end

function _kwcall(call)
    @match call begin
        :($f($(args...))) => begin
            π = invperm(sortperm(collect(args)))
            sargs = Tuple(sort(args))
            targs = Tuple(args)
            quote
                KeywordCalls.baseperm[($f, $sargs)] = $π

                $f(nt::NamedTuple) = kwcall($f, nt)

                $f(; kwargs...) = $f(kwargs.data)

                # $f($(args...)) = $f(NamedTuple{$targs}($(args...)))
            end
        end
        _ => @error "`@kwcall` declaration must be of the form `@kwcall f(b,a,d)`"
    end
end

"""
    kwcall(f, ::NamedTuple)

Dispatch to the permuted `f(::NamedTuple)` call declared using `@kwcall`
"""
@gg function kwcall(::Type{F}, nt::NamedTuple{N}) where {F,N}
    π = Tuple(kwcallperm(F, N))
    Nπ = Tuple((N[p] for p in π))
    quote
        v = values(nt)
        valind(n) = @inbounds v[n]
        $F(NamedTuple{$Nπ}(Tuple(valind.($π))))
    end
end

@gg function kwcall(::F, nt::NamedTuple{N}) where {F<:Function,N}
    f = F.instance
    π = Tuple(kwcallperm(f, N))
    Nπ = Tuple((N[p] for p in π))
    quote
        v = values(nt)
        valind(n) = @inbounds v[n]
        $f(NamedTuple{$Nπ}(Tuple(valind.($π))))
    end
end

end # module


text:
using Test

@testset "register" begin
    include("register.jl")
end

@testset "symengine" begin
    include("symengine/backend.jl")
end


text:
# This file should contain site-specific commands to be executed on Julia startup
# Users may store their own personal commands in the user home directory `homedir()`, in a file named `.juliarc.jl`

using DelimitedFiles

max_stars_repo_name: YingboMa/julia

max_issues_repo_name: JuliaLangEs/julieta

max_forks_repo_name: JuliaLangEs/julieta


text:
#========== tuples ==========#

struct TupleVector{T,DT,L} <: AbstractVector{T}
    data::DT
    TupleVector(tup::DT) where {DT <: Tuple} =
        new{mapreduce(typeof, Base.promote_typejoin, tup), DT, length(tup)}(tup)
end

Base.size(v::TupleVector{T,DT,L}) where {T,DT,L} = (L,)
Base.@propagate_inbounds Base.getindex(v::TupleVector, i::Integer) = getindex(v.data, i)

Core.Tuple(v::TupleVector) = v.data
Base.iterate(v::TupleVector) = iterate(v.data)
Base.iterate(v::TupleVector, state) = iterate(v.data, state)

transmute(tup::Tuple, perm) = transmute(TupleVector(tup), perm)
transmute(tup::Tuple, ::Val{perm}) where {perm} = transmute(TupleVector(tup), Val(perm))

transmutedims(tup::Tuple, perm=(2,1)) = collect(transmute(tup, perm))

function Base.showarg(io::IO, v::TupleVector{T,DT}, toplevel) where {T,DT}
    if all(==(T), DT.parameters)
        toplevel || print(io, "::")
        print(io, "TupleVector{", T, "}")
    else
        print(io, "TupleVector(::", DT, ")")
        toplevel && print(io, " with eltype ", T)
    end
end

#========== dropdims ==========#

# Especially when dropping a trivial dimension, we don't want to produce
# reshape(TransmutedDimsArray(::Array ...

Base.dropdims(A::TransmutedDimsArray; dims) = _dropdims(A, dims...)

_dropdims(A) = A
function _dropdims(A::TransmutedDimsArray{T,N,P}, d::Int, dims...) where {T,N,P}
    if P[d]==0
        perm = ntuple(n -> n<d ? P[n] : P[n+1], N-1)
        newdims = map(n -> n<d ? n : n-1, dims)
        _dropdims(transmute(A.parent, perm), newdims...)
    else
        perm = ntuple(N-1) do n
            Pn = n<d ? P[n] : P[n+1]
            Pn<d ? Pn : Pn-1
        end
        newdims = map(n -> n<d ? n : n-1, dims)
        _dropdims(transmute(dropdims(A.parent; dims=P[d]), perm), newdims...)
    end
end

#=

@btime        dropdims($(transmute(rand(3,3), (2,0,1))), dims=2)
@code_warntype dropdims(transmute(rand(3,3), (2,0,1)), dims=2)

@btime        dropdims($(transmute(rand(3,1,3), (3,2,1))), dims=2)
@code_warntype dropdims(transmute(rand(3,1,3), (3,2,1)), dims=2)

=#

#========== reshape ==========#

function Base.vec(A::TransmutedDimsArray{T,N,P}) where {T,N,P}
    if increasing_or_zero(P)  # the case which transmute() will avoid creating
        vec(A.parent)
    else
        reshape(A, length(A))
    end
end

#========== transpose, etc ==========#

Base.transpose(A::TransmutedDimsArray{<:Number, 1}) = transmute(A, Val((2,1)))
Base.transpose(A::TransmutedDimsArray{<:Number, 2}) = transmute(A, Val((2,1)))
Base.adjoint(A::TransmutedDimsArray{<:Real, 1}) = transmute(A, Val((2,1)))
Base.adjoint(A::TransmutedDimsArray{<:Real, 2}) = transmute(A, Val((2,1)))

Base.PermutedDimsArray(A::TransmutedDimsArray, perm) = transmute(A, perm)

#========== reductions ==========#

# Same strategy as in https://github.com/JuliaLang/julia/pull/39513

function Base.mapreducedim!(f, op, B::AbstractArray, A::TransmutedDimsArray{T,N,P,Q}) where {T,N,P,Q}
    if unique_or_zero(Val(P))
        # any dense transmutation
        Base.mapreducedim!(f, op, transmute(B, Q), parent(A))  # using Val(Q) changes nothing
    else
        # default next step
        Base._mapreducedim!(f, op, B, A)
    end
    B
end

if VERSION > v"1.6-"
    Base._mapreduce_dim(f, op, init::Base._InitialValue, A::TransmutedDimsArray, dims::Colon) =
        _mapreduce_scalar(f, op, init, A, dims)
else
    Base._mapreduce_dim(f, op, init::NamedTuple{()}, A::TransmutedDimsArray, dims::Colon) =
        _mapreduce_scalar(f, op, init, A, dims)
end

@inline function _mapreduce_scalar(f, op, init, A::TransmutedDimsArray{T,N,P}, dims::Colon) where {T,N,P}
    if dims === Colon() && f === identity && op === Base.add_sum
        # safe & easy
        Base._mapreduce_dim(f, op, init, parent(A), dims)
    elseif unique_or_zero(Val(P))
        # any dense transmutation
        Base._mapreduce_dim(f, op, init, parent(A), dims)
    elseif op === Base.add_sum && iszero(f(zero(T)))
        # like sum(::Diagonal)
        Base._mapreduce_dim(f, op, init, parent(A), dims)
    else
        # default next step
        Base._mapreduce(f, op, IndexStyle(A), A)
    end
end

#========== copyto! ==========#

function Base.copyto!(dst::AbstractArray, src::TransmutedDimsArray)
    if axes(dst) == axes(src)
        copy!(dst, src)
    elseif length(dst) == length(src)
        copy!(reshape(dst, axes(src)), src)  # could save a reshape when increasing_or_zero(P)
    elseif length(dst) < length(src)
        throw(BoundsError(dst, lastindex(src)))
    else
        throw(BoundsError(src, lastindex(dst)))
    end
    dst
end

# @propagate_inbounds
function Base.copy!(dst::AbstractArray, src::TransmutedDimsArray{T,N,P,Q}) where {T,N,P,Q}
    @boundscheck axes(dst) == axes(src) || throw(ArgumentError("arrays must have the same axes for copy! (consider using copyto!"))
    if increasing_or_zero(P)  # just a reshape
        copyto!(dst, parent(src))
    else
        if unique_or_zero(P)
            _densecopy_permuted!(dst, parent(src), Val(P))
            # this is happy to reshape... should it be limited to
        else
            fill!(dst, zero(T))  # Diagonal-like
            _copy_into!(dst, parent(src), Val(P))
        end
    end
    dst
end

# For Arrays, this dispatches to use Strided.jl version. Second best:
@generated function _densecopy_permuted!(dst::DenseArray, src::AbstractArray, val::Val{P}) where {P}
    Pminus = filter(!=(0), collect(P))
    if 0 in P
        SB = [:(axes(src,$p)) for p in Pminus]
        Bex = :(reshape(dst, ($(SB...),)))
    else
        Bex = :dst
    end
    if sort(Pminus) == 1:ndims(src)
        Aex = :src
        perm = Tuple(Pminus)
    else
        SA = [:(axes(src,$d)) for d in 1:ndims(src) if d in Pminus]
        Aex = :(reshape(src, ($(SA...),)))
        perm = Tuple(sortperm(Pminus))
    end
    :(permutedims!($Bex, $Aex, $perm); nothing)
end

# Fallback option:
_densecopy_permuted!(dst::AbstractArray, src::AbstractArray, val::Val) =
    _copy_into!(dst, src, val)

function _copy_into!(dst::AbstractArray, parent::AbstractArray, ::Val{P}) where {P}
    @inbounds @simd for I in CartesianIndices(parent)
        J = CartesianIndex(map(p -> p==0 ? 1 : I[p], P))
        dst[J] = parent[I]
    end
    nothing
end

#========== view ==========#

function Base.view(A::TransmutedDimsArray{T,N,P,Q}, inds::Vararg{Union{Int,Colon},N}) where {T,N,P,Q}
    if _is_simple(inds, P)
        parent_inds = genperm_zero(inds, Q, missing)
        view(parent(A), parent_inds...)
    else
        view(A, Base.to_indices(A, inds)...)
    end
end

# Only allow one colon, and P there is not zero
@inline function _is_simple(inds::Tuple, P::NTuple{N,Int}) where {N}
    sum(map(i -> Int(i isa Colon), inds)) == 1 || return false
    n = sum(ntuple(d -> inds[d] isa Colon ? d : 0, N))
    return P[n] != 0
end
# @btime TransmuteDims._is_simple((1,2,:,3), (1,0,2,3)) # 0.041 ns

# @inline function _is_simple(inds, P::Tuple)
#     count(i -> i isa Colon, inds) == 1 || return false
#     n = findfirst(i -> i isa Colon, inds)
#     return P[n] != 0
# end
# @btime _is_simple((1,2,:,3), (1,0,2,3)) # 58.365 ns

#========== the end. ==========#


text:
export get_prenex_normal_form

function handle_negation(formula)
    result = @match formula begin
    f::AFormula => EFormula(f.variable, handle_negation(f.formula))
    f::EFormula => AFormula(f.variable, handle_negation(f.formula))
    _ => Negation(formula)
    end
    return result
end

function handle_junction(formula1, formula2, junction)
    if typeof(formula1) in (AFormula, EFormula)
        return typeof(formula1)(
            formula1.variable,
            handle_junction(formula1.formula, formula2, junction),
        )
    elseif typeof(formula2) in (AFormula, EFormula)
        return typeof(formula2)(
            formula2.variable,
            handle_junction(formula1, formula2.formula, junction),
        )
    else
        return junction(formula1, formula2)
    end
end

"""
    get_prenex_normal_form(formula)

Return the prenex normal form of `formula`.

The returned formula is logically equivalent to `formula`.
"""
function get_prenex_normal_form(formula)
    prepared_formula = get_renamed_quantifiers_form(formula)
    @match prepared_formula begin
    f::Negation => handle_negation(
        get_prenex_normal_form(f.formula),
    )
    f::Conjunction || f::Disjunction => handle_junction(
        get_prenex_normal_form(f.formula1),
        get_prenex_normal_form(f.formula2),
        typeof(f),
    )
    _ => return formula
    end
end


text:
@testset "statesyms ODE" begin
    @named m = samplesystem()
    @test statesyms(m) == [:x, :RHS]
    @test parsyms(m) == [ :τ]
end

@testset "symbol" begin
    @named m = samplesystem()
    @test symbol(m.x) == :m₊x # calls Num which calls Par
    @test symbol(m.τ) == :m₊τ
end;

@testset "strip_namespace" begin
    @test strip_namespace(:x) == :x
    @test strip_namespace(:m₊x) == :x
    @test strip_namespace(:s₊m₊x) == :x
    @test strip_namespace("s.m.x") == "x"
end;

@testset "embed with name different than m" begin
    @named m2 = samplesystem()
    @named sys = embed_system(m2)
    prob = ODEProblem(sys, [m2.x => 0.0], (0.0,10.0), [m2.τ => 3.0])
    sol = solve(prob)
    @test first(sol[m2.x]) == 0.0
    #plot(sol, vars=[m2.x,m2.RHS])
end

max_stars_repo_name: bgctw/MTKHelpers

max_issues_repo_name: bgctw/MTKHelpers

max_forks_repo_name: bgctw/MTKHelpers


text:
type ZeroElementException <: Exception end

#macro debug_stmt(x)
#  x
#end

macro debug_stmt(x)
end


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

JLLWrappers.@generate_main_file_header("Goalign")
JLLWrappers.@generate_main_file("Goalign", UUID("ed1357da-d9e6-57e3-b342-6130e0ca2d5e"))
end  # module Goalign_jll


text:
# ---
# title: 263. Ugly Number
# id: problem263
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Math
# link: <https://leetcode.com/problems/ugly-number/description/>
# hidden: true
# ---
#
# Write a program to check whether a given number is an ugly number.
#
# Ugly numbers are **positive numbers** whose prime factors only include `2, 3,
# 5`.
#
# **Example 1:**
#
#
#
#     Input: 6
#     Output: true
#     Explanation: 6 = 2 × 3
#
# **Example 2:**
#
#
#
#     Input: 8
#     Output: true
#     Explanation: 8 = 2 × 2 × 2
#
#
# **Example 3:**
#
#
#
#     Input: 14
#     Output: false
#     Explanation:14 is not ugly since it includes another prime factor 7.
#
#
# **Note:**
#
#   1. `1` is typically treated as an ugly number.
#   2. Input is within the 32-bit signed integer range: [−231,  231 − 1].
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
#=

Plots reservation wage against the job separation rate

=#

grid_size = 25
alpha_vals = linspace(0.05, 0.5, grid_size)
w_bar_vals = similar(alpha_vals)

mcm = McCallModel()

for (i, alpha) in enumerate(alpha_vals)
    mcm.alpha = alpha
    w_bar = compute_reservation_wage(mcm)
    w_bar_vals[i] = w_bar
end

plot(alpha_vals,
    w_bar_vals,
    lw=2,
    alpha=0.7,
    xlabel="job separation rate",
    ylabel="reservation wage",
    label=L"$\bar w$ as a function of $\alpha$")


text:
include("..\\libs\\FunctionsRobot.jl")

function putmarkersOpposite!(r::Robot, sides::Vector{HorizonSide})
    arr = deepcopy(sides)
    x, y = countingCoordinate(arr)
    sides = [HorizonSide(i) for i in 3:-1:0] #reversed enumerated HorizonSide

    for _ in 1:2
        for coord in (y, x)
            side = pop!(sides)
            moves!(r, side, coord)
            putmarker!(r)
            n = movesAndCounting!(r, side)
            if coord == x
                x = n
            else
                y = n
            end
        end
    end
end

function naprotiv!(r::Robot)
    sides = moveAndReturnDirections!(r)
    putmarkersOpposite!(r, sides)
    moveToStartplace!(r)
    moveToBeginplace!(r, sides)
end


text:
using Distributions
using LinearAlgebra

export rand_constrained_Wish
function rand_constrained_Wish(Psi0, nu, h)
    """
    Psi0 is some pos-def matrix that complies with restrictions imposed by h
    dm is the dimension
    L is the Upper Triangular cholesky factor for Psi0
    A is some random matrix where we impose restrictions to meet our needs
    """

    # On real data, force 0s and 1s where necessary
    # It's like we don't even sample this which is ok I think
    dm = size(Psi0, 1)
    zeroidx = findall(h .== 0)
    nonzeroidx = setdiff(1:1:dm, zeroidx)

    if size(zeroidx,1) == dm
        return Matrix{Float64}(I, dm, dm) * nu
    elseif size(zeroidx,1) == 0
        subh = h
        subPsi0 = Psi0
    else
        subh = h[nonzeroidx]
        subPsi0 = @inbounds Psi0[nonzeroidx, nonzeroidx]
    end

    subdm = size(subPsi0, 1)

    U = cholesky(Hermitian(subPsi0)).U
    A = zeros((subdm, subdm))

    # Sample the diagonal of A
    for i in 1:subdm
        A[i,i] = sqrt(rand(Chisq(nu - i + 1)))
    end

    # Construct the pairs of dimensions of A
    # in the order in which we need to simulate them
    npairs = binomial(dm,2)
    dimpairs = Matrix{Int64}(undef, npairs, 2)
    counter = 1
    for k = 1:1:(dm-1)
        for i=1:1:k
            #global counter
            @inbounds dimpairs[counter,:] = [i, k+1]
            counter = counter + 1
        end
    end

    if size(zeroidx, 1) > 0
        rmpairs = findall(mapslices(x -> size(intersect(x, zeroidx), 1), dimpairs; dims = 2)[:,1] .!= 0)
        dimpairs = dimpairs[setdiff(1:1:npairs, rmpairs), :]
        npairs = size(dimpairs,1)
    end

    # Simulate the IW matrix, starting from the top left corner
    # and working our way down
    for idx = 1:1:npairs
        @inbounds m = dimpairs[idx,1]
        @inbounds p = dimpairs[idx,2]

        mstar = findall(nonzeroidx .== m)[1]
        pstar = findall(nonzeroidx .== p)[1]

        # Compute bounds on A[m,n]
        if mstar == 1
            @inbounds premult = U[mstar,mstar] * A[mstar,mstar]
            term4 = 0.0
            for i in mstar:(pstar-1)
                #global term4
                @inbounds term4 += U[i,pstar] * A[i,mstar]
            end

            if @inbounds sign(h[m] * h[p]) == -1
                lb = -Inf
                ub = @inbounds ( -term4 ) / U[pstar, pstar]#(0 - outterm4) / denom
            else
                lb = @inbounds ( -term4 ) / U[pstar, pstar]#(0 - outterm4) / denom
                ub = Inf
            end
            #llb = min(lb, ub)
            #uub = max(lb, ub)
            #@inbounds A[pstar,mstar] = rand(truncated(Normal(), llb, uub))
            @inbounds A[pstar,mstar] = rand(truncated(Normal(), lb, ub))
        else
            term1 = 0.0
            for k in 1:(mstar-1)
                outerterm = 0.0
                for j in 1:k
                    innerterm = 0.0
                    for i in j:pstar
                        @inbounds innerterm += (A[i,j] * U[i,pstar])
                    end
                    @inbounds outerterm += (innerterm * A[k,j])
                end
                #global term1
                @inbounds term1 += (U[k,mstar] * outerterm)
            end

            outerterm2 = 0.0
            for j in 1:(mstar-1)
                innerterm2 = 0.0
                for i in j:pstar
                    @inbounds innerterm2 += (A[i,j] * U[i,pstar])
                end
                #global outerterm2
                @inbounds outerterm2 += (innerterm2 * A[mstar,j])
            end
            @inbounds term2 = U[mstar,mstar] * outerterm2

            term3 = 0.0
            for i in mstar:(pstar-1)
                #global term3
                @inbounds term3 += (A[i,mstar] * U[i,pstar])
            end
            @inbounds term3 *= (A[mstar,mstar] * U[mstar,mstar])

            if @inbounds sign(h[m] * h[p]) == -1
                lb = -Inf
                ub = @inbounds ((-term1 - term2 - term3) / (U[pstar,pstar] * U[mstar,mstar] * A[mstar,mstar]) )
            else
                lb = @inbounds ((-term1 - term2 - term3) / (U[pstar,pstar] * U[mstar,mstar] * A[mstar,mstar]) )
                ub = Inf
            end
            #llb = min(lb, ub)
            #uub = max(lb, ub)
            #@inbounds A[pstar,mstar] = rand(truncated(Normal(), llb, uub))
            @inbounds A[pstar,mstar] = rand(truncated(Normal(), lb, ub))
        end
    end

    if size(zeroidx,1) == 0
        return U' * A * A' * U
    else
        matmult = (U' * A * A' * U)
        out = Matrix{Float64}(I, dm, dm)
        out[nonzeroidx, nonzeroidx] = matmult

        [out[zz, zz] = nu for zz in zeroidx]
        return out
    end
end

export rand_constrained_MVN
function rand_constrained_MVN(Sigma, mu0, h)
    dm = size(Sigma, 1)
    zeroidx = findall(h .== 0)
    nonzeroidx = setdiff(1:1:dm, zeroidx)

    if size(zeroidx,1) == dm
        return zeros(dm)
    elseif size(zeroidx,1) == 0
        subh = h
        subSigma = Sigma
        submu0 = mu0
    else
        subh = @inbounds h[nonzeroidx]
        subSigma = @inbounds Sigma[nonzeroidx, nonzeroidx]
        submu0 = @inbounds mu0[nonzeroidx]
    end

    subdm = size(nonzeroidx, 1)
    A = cholesky(Hermitian(subSigma)).L

    z = Array{Float64,1}(undef, subdm)
    for i in 1:1:subdm
        bound = -mu0[i]
        if i == 1
            bound /= A[i,i]
        else
            bound = @inbounds (bound - A[i,1:1:(i-1)]' * z[1:1:(i-1)]) / A[i,i]
        end

        if subh[i] == 1
            z[i] = rand(truncated(Normal(), bound, Inf))
        else
            z[i] = rand(truncated(Normal(), -Inf, bound))
        end
        bound = nothing
    end

    subout = A * z .+ submu0
    out = zeros(dm)
    out[nonzeroidx] = subout
    out
end


text:
struct CorrelatedOU <: TimeSeriesModel{2,Float64}
    σ::Float64
    θ::Float64
    ρ::Float64
    σ²::Float64
    θ²::Float64
    function CorrelatedOU(σ,θ,ρ)
        σ > 0 || throw(ArgumentError("OU process requires 0 < σ."))
        θ > 0 || throw(ArgumentError("OU process requires 0 < θ."))
        1>ρ>0 || throw(ArgumentError("OU process requires 0 < ρ < 1."))
        new(σ,θ,ρ,σ^2,θ^2)
    end
function CorrelatedOU(x::AbstractVector{Float64})
        @boundscheck checkparameterlength(x,CorrelatedOU)
        @inbounds CorrelatedOU(x[1], x[2], x[3])
    end
end

npars(::Type{CorrelatedOU}) = 3
lowerbounds(::Type{CorrelatedOU}) = [0,0,0]
upperbounds(::Type{CorrelatedOU}) = [Inf,Inf,1.0]

@propagate_inbounds function add_sdf!(out, m::CorrelatedOU, ω)
    @boundscheck checkbounds(out,1:3)
    @inbounds begin
        s = m.σ²/(π*(m.θ²+ω^2))
        out[1] += s
        out[2] += s*m.ρ
        out[3] += s
    end
end
@propagate_inbounds function acv!(out, m::CorrelatedOU, τ::Number)
    @boundscheck checkbounds(out,1:3)
    @inbounds begin
        a = exp(-m.θ*abs(τ)) * m.σ² / m.θ
        out[1] = a
        out[2] = a*m.ρ
        out[3] = a
    end
end

@propagate_inbounds function grad_add_sdf!(out, m::CorrelatedOU, ω)
    @boundscheck checkbounds(out,1:3,1:3)
    @inbounds begin
        σ, θ, ρ = m.σ, m.θ, m.ρ
        σ_part = 2σ / (π*(m.θ²+ω^2))
        # ∂σ
        out[1,1] += σ_part
        out[2,1] += σ_part * ρ
        out[3,1] += σ_part
        # ∂θ
        θ_part = -π*θ/2*(σ_part)^2
        out[1,2] += θ_part
        out[2,2] += θ_part * ρ
        out[3,2] += θ_part
        # ∂ρ
        # out[1,3] += 0
        out[2,3] += σ_part * σ / 2
        # out[3,3] += 0
    end
    nothing
end

@propagate_inbounds function hess_add_sdf!(out, m::CorrelatedOU, ω)
    @boundscheck checkbounds(out,1:3,1:5)
    @inbounds begin
        σ, θ, ρ = m.σ, m.θ, m.ρ
        σ_part = 2 / (π*(m.θ²+ω^2))
        θσ_part = -π*θ*σ*(σ_part)^2
        θ_part = -π*σ_part*σ*(σ_part*σ/2 + θ*θσ_part)

        # ∂σ²
        out[1,1] += σ_part
        out[2,1] += σ_part * ρ
        out[3,1] += σ_part

        # ∂σθ
        out[1,2] += θσ_part
        out[2,2] += θσ_part * ρ
        out[3,2] += θσ_part

        #∂σρ
        # out[1,3] += 0
        out[2,3] += σ_part * σ
        # out[3,3] += 0

        # ∂θ²
        out[1,4] += θ_part
        out[2,4] += θ_part * ρ
        out[3,4] += θ_part

        #∂θρ
        # out[1,5] += 0
        out[2,5] += σ/2 * θσ_part
        # out[3,5] += 0

        #∂ρ² all zero
    end

    nothing
end

@propagate_inbounds function grad_acv!(out, m::CorrelatedOU, τ::Number)
    @boundscheck checkbounds(out,1:3,1:3)
    @inbounds begin
        σ, θ, ρ = m.σ, m.θ, m.ρ
        absτ = abs(τ)
        σ_part = 2exp(-θ*absτ)*σ/θ
        # ∂σ
        out[1,1] = σ_part
        out[2,1] = σ_part * ρ
        out[3,1] = σ_part
        # ∂θ
        θ_part = -σ_part*σ*(1/θ+absτ)/2
        out[1,2] = θ_part
        out[2,2] = θ_part * ρ
        out[3,2] = θ_part
        # ∂ρ
        # out[1,3] = 0
        out[2,3] = σ_part * σ/2
        # out[3,3] = 0
    end
    nothing
end

@propagate_inbounds function hess_acv!(out, m::CorrelatedOU, τ::Number)
    @boundscheck checkbounds(out,1:3,1:5)
    @inbounds begin
        σ, θ, ρ = m.σ, m.θ, m.ρ
        absτ = abs(τ)
        σ_part = 2exp(-θ*absτ)/θ
        θσ_part = -σ_part*σ*(1/θ+absτ)
        θ_part = σ_part*(m.σ²)*(1/(m.θ²) + absτ/θ + (absτ^2)/2)

        # ∂σ²
        out[1,1] = σ_part
        out[2,1] = σ_part * ρ
        out[3,1] = σ_part

        # ∂σθ
        out[1,2] = θσ_part
        out[2,2] = θσ_part * ρ
        out[3,2] = θσ_part

        #∂σρ
        # out[1,3] = 0
        out[2,3] = σ_part * σ
        # out[3,3] = 0

        # ∂θ²
        out[1,4] = θ_part
        out[2,4] = θ_part * ρ
        out[3,4] = θ_part

        #∂θρ
        # out[1,5] = 0
        out[2,5] = -σ_part*m.σ²*(1/θ+absτ)/2
        # out[3,5] = 0

        #∂ρ² all zero
    end
    nothing
end


text:
using Test

@testset "Core" begin
  include("MarkovCategories.jl")
end


text:
module SingleLayerQG

export
  Problem,
  set_q!,
  updatevars!,

  energy,
  kinetic_energy,
  potential_energy,
  energy_dissipation,
  energy_work,
  energy_drag,
  enstrophy,
  enstrophy_dissipation,
  enstrophy_work,
  enstrophy_drag

using
  CUDA,
  Reexport,
  DocStringExtensions

@reexport using FourierFlows

using LinearAlgebra: mul!, ldiv!
using FourierFlows: parsevalsum, parsevalsum2

nothingfunction(args...) = nothing

"""
    Problem(dev::Device=CPU();
                      nx = 256,
                      ny = nx,
                      Lx = 2π,
                      Ly = Lx,
                       β = 0.0,
      deformation_radius = Inf,
                     eta = nothing,
                       ν = 0.0,
                      nν = 1,
                       μ = 0.0,
                      dt = 0.01,
                 stepper = "RK4",
                   calcF = nothingfunction,
              stochastic = false,
        aliased_fraction = 1/3,
                       T = Float64)

Construct a single-layer quasi-geostrophic `problem` on device `dev`.

Keyword arguments
=================
    - `dev`: (required) `CPU()` or `GPU()`; computer architecture used to time-step `problem`.
    - `nx`: Number of grid points in ``x``-domain.
    - `ny`: Number of grid points in ``y``-domain.
    - `Lx`: Extent of the ``x``-domain.
    - `Ly`: Extent of the ``y``-domain.
    - `β`: Planetary vorticity ``y``-gradient.
    - `deformation_radius`: Rossby radius of deformation; set `Inf` for purely barotropic.
    - `eta`: Topographic potential vorticity.
    - `ν`: Small-scale (hyper)-viscosity coefficient.
    - `nν`: (Hyper)-viscosity order, `nν```≥ 1``".
    - `μ`: Linear drag coefficient.
    - `dt`: Time-step.
    - `stepper`: The extent of the ``y``-domain.
    - `calcF`: Function that calculates the Fourier transform of the forcing, ``F̂``.
    - `stochastic`: `true` or `false`; boolean denoting whether `calcF` is temporally stochastic.
    - `aliased_fraction`: the fraction of high-wavenubers that are zero-ed out by `dealias!()`.
    - `T`: `Float32` or `Float64`; floating point type used for `problem` data.
"""
function Problem(dev::Device=CPU();
  # Numerical parameters
                  nx = 256,
                  ny = nx,
                  Lx = 2π,
                  Ly = Lx,
  # Physical parameters
                   β = 0.0,
  deformation_radius = Inf,
                 eta = nothing,
  # Drag and (hyper-)viscosity
                   ν = 0.0,
                  nν = 1,
                   μ = 0.0,
  # Timestepper and equation options
                  dt = 0.01,
             stepper = "RK4",
               calcF = nothingfunction,
          stochastic = false,
  # Float type and dealiasing
    aliased_fraction = 1/3,
                   T = Float64)

  # the grid
  grid = TwoDGrid(dev, nx, Lx, ny, Ly; aliased_fraction=aliased_fraction, T=T)
  x, y = gridpoints(grid)

  # topographic PV
  eta === nothing && (eta = zeros(dev, T, (nx, ny)))

  params = deformation_radius == Inf ? BarotropicQGParams(grid, T(β), eta, T(μ), T(ν), nν, calcF) : EquivalentBarotropicQGParams(grid, T(β), T(deformation_radius), eta, T(μ), T(ν), nν, calcF)

  vars = calcF == nothingfunction ? DecayingVars(dev, grid) : (stochastic ? StochasticForcedVars(dev, grid) : ForcedVars(dev, grid))

  equation = Equation(params, grid)

  return FourierFlows.Problem(equation, stepper, dt, grid, vars, params, dev)
end

# ----------
# Parameters
# ----------

abstract type SingleLayerQGParams <: AbstractParams end

"""
    Params{T, Aphys, Atrans, ℓ}(β, deformation_radius, eta, etah, μ, ν, nν, calcF!)

A struct containing the parameters for the SingleLayerQG problem. Included are:

$(TYPEDFIELDS)
"""
struct Params{T, Aphys, Atrans, ℓ} <: SingleLayerQGParams
    "planetary vorticity y-gradient"
                   β :: T
    "Rossby radius of deformation"
  deformation_radius :: ℓ
    "topographic potential vorticity"
                 eta :: Aphys
    "Fourier transform of topographic potential vorticity"
                etah :: Atrans
    "linear drag coefficient"
                   μ :: T
    "small-scale (hyper)-viscosity coefficient"
                   ν :: T
    "(hyper)-viscosity order, `nν```≥ 1``"
                  nν :: Int
    "function that calculates the Fourier transform of the forcing, ``F̂``"
              calcF! :: Function
end

const BarotropicQGParams = Params{<:AbstractFloat, <:AbstractArray, <:AbstractArray, Nothing}
const EquivalentBarotropicQGParams = Params{<:AbstractFloat, <:AbstractArray, <:AbstractArray, <:AbstractFloat}

"""
    EquivalentBarotropicQGParams(grid::TwoDGrid, β, deformation_radius, eta, μ, ν, nν::Int, calcF

Constructor for EquivalentBarotropicQGParams (finite Rossby radius of deformation).
"""
function EquivalentBarotropicQGParams(grid::AbstractGrid{T, A}, β, deformation_radius, eta, μ, ν, nν::Int, calcF) where {T, A}
  eta_grid = typeof(eta) <: AbstractArray ? eta : FourierFlows.on_grid(eta, grid)
  eta_gridh = rfft(eta_grid)
  return Params(β, deformation_radius, eta_grid, eta_gridh, μ, ν, nν, calcF)
end

"""
    BarotropicQGParams(grid::TwoDGrid, β, eta, μ, ν, nν::Int, calcF

Constructor for BarotropicQGParams (infinite Rossby radius of deformation).
"""
BarotropicQGParams(grid::AbstractGrid, β, eta, μ, ν, nν::Int, calcF) =
    EquivalentBarotropicQGParams(grid, β, nothing, eta, μ, ν, nν, calcF)

# ---------
# Equations
# ---------

"""
    Equation(params::BarotropicQGParams, grid)

Return the `equation` for a barotropic QG problem with `params` and `grid`. Linear operator
``L`` includes bottom drag ``μ``, (hyper)-viscosity of order ``n_ν`` with coefficient ``ν``
and the ``β`` term:

```math
L = - μ - ν |𝐤|^{2 n_ν} + i β k_x / |𝐤|² .
```
Nonlinear term is computed via `calcN!` function.
"""
function Equation(params::BarotropicQGParams, grid::AbstractGrid)
  L = @. - params.μ - params.ν * grid.Krsq^params.nν + im * params.β * grid.kr * grid.invKrsq
  CUDA.@allowscalar L[1, 1] = 0

  return FourierFlows.Equation(L, calcN!, grid)
end

"""
    Equation(params::EquivalentBarotropicQGParams, grid)

Return the `equation` for an equivalent-barotropic QG problem with `params` and `grid`.
Linear operator ``L`` includes bottom drag ``μ``, (hyper)-viscosity of order ``n_ν`` with
coefficient ``ν`` and the ``β`` term:

```math
L = -μ - ν |𝐤|^{2 n_ν} + i β k_x / (|𝐤|² + 1/ℓ²) .
```
Nonlinear term is computed via `calcN!` function.
"""
function Equation(params::EquivalentBarotropicQGParams, grid::AbstractGrid)
  L = @. - params.μ - params.ν * grid.Krsq^params.nν + im * params.β * grid.kr / (grid.Krsq + 1 / params.deformation_radius^2)
  CUDA.@allowscalar L[1, 1] = 0

  return FourierFlows.Equation(L, calcN!, grid)
end

# ----
# Vars
# ----

abstract type SingleLayerQGVars <: AbstractVars end

"""
    Vars{Aphys, Atrans, F, P}(q, ψ, u, v, qh, , ψh, uh, vh, Fh, prevsol)

The variables for SingleLayer QG:

$(FIELDS)
"""
struct Vars{Aphys, Atrans, F, P} <: SingleLayerQGVars
    "relative vorticity (+ vortex stretching)"
        q :: Aphys
    "streamfunction"
        ψ :: Aphys
    "x-component of velocity"
        u :: Aphys
    "y-component of velocity"
        v :: Aphys
    "Fourier transform of relative vorticity (+ vortex stretching)"
       qh :: Atrans
    "Fourier transform of streamfunction"
       ψh :: Atrans
    "Fourier transform of x-component of velocity"
       uh :: Atrans
    "Fourier transform of y-component of velocity"
       vh :: Atrans
    "Fourier transform of forcing"
       Fh :: F
    "`sol` at previous time-step"
  prevsol :: P
end

const ForcedVars = Vars{<:AbstractArray, <:AbstractArray, <:AbstractArray, Nothing}
const StochasticForcedVars = Vars{<:AbstractArray, <:AbstractArray, <:AbstractArray, <:AbstractArray}

"""
    DecayingVars(dev, grid)

Return the `vars` for unforced single-layer QG problem on device `dev` and with `grid`
"""
function DecayingVars(dev::Dev, grid::AbstractGrid) where Dev
  T = eltype(grid)

  @devzeros Dev T (grid.nx, grid.ny) q u v ψ
  @devzeros Dev Complex{T} (grid.nkr, grid.nl) qh uh vh ψh

  Vars(q, ψ, u, v, qh, ψh, uh, vh, nothing, nothing)
end

"""
    ForcedVars(dev, grid)

Return the `vars` for forced single-layer QG problem on device dev and with `grid`.
"""
function ForcedVars(dev::Dev, grid::AbstractGrid) where Dev
  T = eltype(grid)

  @devzeros Dev T (grid.nx, grid.ny) q u v ψ
  @devzeros Dev Complex{T} (grid.nkr, grid.nl) qh uh vh ψh Fh

  return Vars(q, ψ, u, v, qh, ψh, uh, vh, Fh, nothing)
end

"""
    StochasticForcedVars(dev, grid)

Return the vars for stochastically forced barotropic QG problem on device dev and with `grid`.
"""
function StochasticForcedVars(dev::Dev, grid::AbstractGrid) where Dev
  T = eltype(grid)
  @devzeros Dev T (grid.nx, grid.ny) q u v ψ
  @devzeros Dev Complex{T} (grid.nkr, grid.nl) qh uh vh ψh Fh prevsol

  return Vars(q, ψ, u, v, qh, ψh, uh, vh, Fh, prevsol)
end

# -------
# Solvers
# -------

"""
    calcN_advection!(N, sol, t, clock, vars, params, grid)

Calculate the Fourier transform of the advection term, ``- 𝖩(ψ, q+η)`` in conservative
form, i.e., ``- ∂_x[(∂_y ψ)(q+η)] - ∂_y[(∂_x ψ)(q+η)]`` and store it in `N`:

```math
N = - \\widehat{𝖩(ψ, q+η)} = - i k_x \\widehat{u (q+η)} - i k_y \\widehat{v (q+η)} .
```
"""
function calcN_advection!(N, sol, t, clock, vars, params, grid)
  @. vars.qh = sol
  streamfunctionfrompv!(vars.ψh, vars.qh, params, grid)
  @. vars.uh = -im * grid.l  * vars.ψh
  @. vars.vh =  im * grid.kr * vars.ψh

  ldiv!(vars.q, grid.rfftplan, vars.qh)
  ldiv!(vars.u, grid.rfftplan, vars.uh)
  ldiv!(vars.v, grid.rfftplan, vars.vh)

  uq_plus_η = vars.u                                            # use vars.u as scratch variable
  @. uq_plus_η *= vars.q + params.eta                           # u*(q+η)
  vq_plus_η = vars.v                                            # use vars.v as scratch variable
  @. vq_plus_η *= vars.q + params.eta                           # v*(q+η)

  uq_plus_ηh = vars.uh                                          # use vars.uh as scratch variable
  mul!(uq_plus_ηh, grid.rfftplan, uq_plus_η)                    # \hat{u*(q+η)}
  vq_plus_ηh = vars.vh                                          # use vars.vh as scratch variable
  mul!(vq_plus_ηh, grid.rfftplan, vq_plus_η)                    # \hat{v*(q+η)}

  @. N = -im * grid.kr * uq_plus_ηh - im * grid.l * vq_plus_ηh  # -∂[u*(q+η)]/∂x -∂[v*(q+η)]/∂y
  return nothing
end

"""
    calcN!(N, sol, t, clock, vars, params, grid)

Calculate the nonlinear term, that is the advection term and the forcing,

```math
N = - \\widehat{𝖩(ψ, q+η)} + F̂ .
```
"""
function calcN!(N, sol, t, clock, vars, params, grid)
  dealias!(sol, grid)

  calcN_advection!(N, sol, t, clock, vars, params, grid)

  addforcing!(N, sol, t, clock, vars, params, grid)

  return nothing
end

"""
    addforcing!(N, sol, t, clock, vars, params, grid)

When the problem includes forcing, calculate the forcing term ``F̂`` and add it to the
nonlinear term ``N``.
"""
addforcing!(N, sol, t, clock, vars::Vars, params, grid) = nothing

function addforcing!(N, sol, t, clock, vars::ForcedVars, params, grid)
  params.calcF!(vars.Fh, sol, t, clock, vars, params, grid)
  @. N += vars.Fh

  return nothing
end

function addforcing!(N, sol, t, clock, vars::StochasticForcedVars, params, grid)
  if t == clock.t # not a substep
    @. vars.prevsol = sol # sol at previous time-step is needed to compute budgets for stochastic forcing
    params.calcF!(vars.Fh, sol, t, clock, vars, params, grid)
  end

  @. N += vars.Fh

  return nothing
end

"""
    streamfunctionfrompv!(ψh, qh, params, grid)

Invert the Fourier transform of PV `qh` to obtain the Fourier transform of the streamfunction `ψh`.
"""
function streamfunctionfrompv!(ψh, qh, params::BarotropicQGParams, grid)
  @. ψh =  - grid.invKrsq * qh
end

function streamfunctionfrompv!(ψh, qh, params::EquivalentBarotropicQGParams, grid)
  @. ψh  = - qh / (grid.Krsq + 1 / params.deformation_radius^2)
end

# ----------------
# Helper functions
# ----------------

"""
    updatevars!(sol, vars, params, grid)

Update the variables in `vars` with the solution in `sol`.
"""
function updatevars!(sol, vars, params, grid)
  dealias!(sol, grid)

  @. vars.qh = sol
  streamfunctionfrompv!(vars.ψh, vars.qh, params, grid)
  @. vars.uh = -im * grid.l  * vars.ψh
  @. vars.vh =  im * grid.kr * vars.ψh

  ldiv!(vars.q, grid.rfftplan, deepcopy(vars.qh))
  ldiv!(vars.ψ, grid.rfftplan, deepcopy(vars.ψh))
  ldiv!(vars.u, grid.rfftplan, deepcopy(vars.uh))
  ldiv!(vars.v, grid.rfftplan, deepcopy(vars.vh))

  return nothing
end

updatevars!(prob) = updatevars!(prob.sol, prob.vars, prob.params, prob.grid)

"""
    set_q!(prob, q)

Set the solution of problem, `prob.sol` as the transform of ``q`` and update variables `prob.vars`.
"""
function set_q!(prob, q)
  sol, vars, params, grid = prob.sol, prob.vars, prob.params, prob.grid

  mul!(vars.qh, grid.rfftplan, q)
  @. sol = vars.qh

  updatevars!(sol, vars, params, grid)

  return nothing
end

"""
    kinetic_energy(prob)
    kinetic_energy(sol, grid, vars, params)

Return the domain-averaged kinetic energy of the fluid. Since ``u² + v² = |{\\bf ∇} ψ|²``, the
domain-averaged kinetic energy is
```math
\\int \\frac1{2} |{\\bf ∇} ψ|² \\frac{𝖽x 𝖽y}{L_x L_y} = \\sum_{𝐤} \\frac1{2} |𝐤|² |ψ̂|² .
```
"""
function kinetic_energy(sol, vars, params, grid)
  streamfunctionfrompv!(vars.ψh, sol, params, grid)
  @. vars.uh = sqrt.(grid.Krsq) * vars.ψh      # vars.uh is a dummy variable

  return parsevalsum2(vars.uh , grid) / (2 * grid.Lx * grid.Ly)
end

kinetic_energy(prob) = kinetic_energy(prob.sol, prob.vars, prob.params, prob.grid)

"""
    potential_energy(prob)
    potential_energy(sol, grid, vars, params)

Return the domain-averaged potential energy of the fluid,
```math
\\int \\frac1{2} \\frac{ψ²}{ℓ²} \\frac{𝖽x 𝖽y}{L_x L_y} = \\sum_{𝐤} \\frac1{2} \\frac{|ψ̂|²}{ℓ²} .
```
"""
function potential_energy(sol, vars, params::EquivalentBarotropicQGParams, grid)
  streamfunctionfrompv!(vars.ψh, sol, params, grid)
  return 1 / params.deformation_radius^2 * parsevalsum2(vars.ψh, grid) / (2 * grid.Lx * grid.Ly)
end

@inline potential_energy(sol, vars, params::BarotropicQGParams, grid) = 0

@inline potential_energy(prob) = potential_energy(prob.sol, prob.vars, prob.params, prob.grid)

"""
    energy(prob)
    energy(sol, grid, vars, params)

Return the domain-averaged total energy of the fluid, that is, the kinetic energy for a
pure barotropic flow or the sum of kinetic and potential energies for an equivalent barotropic
flow.
"""
@inline energy(prob) = energy(prob.sol, prob.vars, prob.params, prob.grid)

@inline energy(sol, vars, params::BarotropicQGParams, grid) = kinetic_energy(sol, vars, params, grid)

@inline energy(sol, vars, params::EquivalentBarotropicQGParams, grid) = kinetic_energy(sol, vars, params, grid) + potential_energy(sol, vars, params, grid)

"""
    enstrophy(prob)
    enstrophy(sol, vars, params, grid)

Return the domain-averaged enstrophy
```math
\\int \\frac1{2} (q + η)² \\frac{𝖽x 𝖽y}{L_x L_y} = \\sum_{𝐤} \\frac1{2} |q̂ + η̂|² .
```
"""
function enstrophy(sol, vars, params, grid)
  @. vars.qh = sol
  return parsevalsum2(vars.qh + params.etah, grid) / (2 * grid.Lx * grid.Ly)
end

@inline enstrophy(prob) = enstrophy(prob.sol, prob.vars, prob.params, prob.grid)

"""
    energy_dissipation(prob)
    energy_dissipation(sol, vars, params, grid)

Return the domain-averaged energy dissipation rate. nν must be >= 1.
"""
@inline function energy_dissipation(sol, vars, params::BarotropicQGParams, grid)
  energy_dissipationh = vars.uh # use vars.uh as scratch variable

  @. energy_dissipationh = params.ν * grid.Krsq^(params.nν-1) * abs2(sol)
  CUDA.@allowscalar energy_dissipationh[1, 1] = 0
  return 1 / (grid.Lx * grid.Ly) * parsevalsum(energy_dissipationh, grid)
end

energy_dissipation(sol, vars, params::EquivalentBarotropicQGParams, grid) = error("not implemented for finite deformation radius")

@inline energy_dissipation(prob) = energy_dissipation(prob.sol, prob.vars, prob.params, prob.grid)

"""
    enstrophy_dissipation(prob)
    enstrophy_dissipation(sol, vars, params, grid)

Return the domain-averaged enstrophy dissipation rate. nν must be >= 1.
"""
@inline function enstrophy_dissipation(sol, vars, params::BarotropicQGParams, grid)
  enstrophy_dissipationh = vars.uh # use vars.uh as scratch variable

  @. enstrophy_dissipationh = params.ν * grid.Krsq^params.nν * abs2(sol)
  CUDA.@allowscalar enstrophy_dissipationh[1, 1] = 0
  return 1 / (grid.Lx * grid.Ly) * parsevalsum(enstrophy_dissipationh, grid)
end

@inline enstrophy_dissipation(sol, vars, params::EquivalentBarotropicQGParams, grid) = error("not implemented for finite deformation radius")

@inline enstrophy_dissipation(prob) = enstrophy_dissipation(prob.sol, prob.vars, prob.params, prob.grid)

"""
    energy_work(prob)
    energy_work(sol, vars, params, grid)

Return the domain-averaged rate of work of energy by the forcing `Fh`.
"""
@inline function energy_work(sol, vars::ForcedVars, params::BarotropicQGParams, grid)
  energy_workh = vars.uh # use vars.uh as scratch variable

  @. energy_workh = grid.invKrsq * sol * conj(vars.Fh)
  return 1 / (grid.Lx * grid.Ly) * parsevalsum(energy_workh, grid)
end

@inline function energy_work(sol, vars::StochasticForcedVars, params::BarotropicQGParams, grid)
  energy_workh = vars.uh # use vars.uh as scratch variable

  @. energy_workh = grid.invKrsq * (vars.prevsol + sol)/2 * conj(vars.Fh)
  return 1 / (grid.Lx * grid.Ly) * parsevalsum(vars.uh, grid)
end

@inline energy_work(sol, vars, params::EquivalentBarotropicQGParams, grid) = error("not implemented for finite deformation radius")

@inline energy_work(prob) = energy_work(prob.sol, prob.vars, prob.params, prob.grid)

"""
    enstrophy_work(prob)
    enstrophy_work(sol, vars, params, grid)

Return the domain-averaged rate of work of enstrophy by the forcing `Fh`.
"""
@inline function enstrophy_work(sol, vars::ForcedVars, params::BarotropicQGParams, grid)
  enstrophy_workh = vars.uh # use vars.uh as scratch variable

  @. enstrophy_workh = sol * conj(vars.Fh)
  return 1 / (grid.Lx * grid.Ly) * parsevalsum(enstrophy_workh, grid)
end

@inline function enstrophy_work(sol, vars::StochasticForcedVars, params::BarotropicQGParams, grid)
  enstrophy_workh = vars.uh # use vars.uh as scratch variable

  @. enstrophy_workh = (vars.prevsol + sol) / 2 * conj(vars.Fh)
  return 1 / (grid.Lx * grid.Ly) * parsevalsum(enstrophy_workh, grid)
end

@inline enstrophy_work(sol, vars, params::EquivalentBarotropicQGParams, grid) = error("not implemented for finite deformation radius")

@inline enstrophy_work(prob) = enstrophy_work(prob.sol, prob.vars, prob.params, prob.grid)

"""
    energy_drag(prob)
    energy_drag(sol, vars, params, grid)

Return the extraction of domain-averaged energy by drag μ.
"""
@inline function energy_drag(sol, vars, params::BarotropicQGParams, grid)
  energy_dragh = vars.uh # use vars.uh as scratch variable

  @. energy_dragh = params.μ * grid.invKrsq * abs2(sol)
  CUDA.@allowscalar energy_dragh[1, 1] = 0
  return  1 / (grid.Lx * grid.Ly) * parsevalsum(energy_dragh, grid)
end

@inline energy_drag(sol, vars, params::EquivalentBarotropicQGParams, grid) = error("not implemented for finite deformation radius")

@inline energy_drag(prob) = energy_drag(prob.sol, prob.vars, prob.params, prob.grid)

"""
    enstrophy_drag(prob)
    enstrophy_drag(sol, vars, params, grid)

Return the extraction of domain-averaged enstrophy by drag/hypodrag μ.
"""
@inline function enstrophy_drag(sol, vars, params::BarotropicQGParams, grid)
  enstrophy_dragh = vars.uh # use vars.uh as scratch variable

  @. enstrophy_dragh = params.μ * abs2(sol)
  CUDA.@allowscalar enstrophy_dragh[1, 1] = 0
  return 1 / (grid.Lx * grid.Ly) * parsevalsum(enstrophy_dragh, grid)
end

@inline enstrophy_drag(sol, vars, params::EquivalentBarotropicQGParams, grid) = error("not implemented for finite deformation radius")

@inline enstrophy_drag(prob) = enstrophy_drag(prob.sol, prob.vars, prob.params, prob.grid)

end # module


text:
using Clang.Generators
using LLVM_full_jll

cd(@__DIR__)

llvm_include_dir = joinpath(LLVM_full_jll.artifact_dir, "include") |> normpath
include_dir = joinpath(ENV["LIBCLANGEX_INSTALL_PREFIX"], "include") |> normpath

options = load_options(joinpath(@__DIR__, "generator.toml"))

args = Generators.get_default_args()
push!(args, "-isystem$llvm_include_dir")
push!(args, "-I$include_dir")

headers = detect_headers(include_dir, args)

@add_def time_t
@add_def LLVMModuleRef
@add_def LLVMOpaqueModule
@add_def LLVMOpaqueContext
@add_def LLVMContextRef
@add_def LLVMMemoryBufferRef
@add_def LLVMGenericValueRef

ctx = create_context(headers, args, options)

build!(ctx)


text:
#############################################################################
# transpose.jl
# Returns the transpose of a matrix
# All expressions and atoms are subtpyes of AbstractExpr.
# Please read expressions.jl first.
#############################################################################

import Base.transpose, Base.adjoint

struct TransposeAtom <: AbstractExpr
    head::Symbol
    id_hash::UInt64
    children::Tuple{AbstractExpr}
    size::Tuple{Int,Int}

    function TransposeAtom(x::AbstractExpr)
        children = (x,)
        return new(:transpose, hash(children), children, (x.size[2], x.size[1]))
    end
end

function sign(x::TransposeAtom)
    return sign(x.children[1])
end

function monotonicity(x::TransposeAtom)
    return (Nondecreasing(),)
end

function curvature(x::TransposeAtom)
    return ConstVexity()
end

function evaluate(x::TransposeAtom)
    return transpose(evaluate(x.children[1]))
end

# Since everything is vectorized, we simply need to multiply x by a permutation
# matrix such that coeff * vectorized(x) - vectorized(x') = 0
function conic_form!(x::TransposeAtom, unique_conic_forms::UniqueConicForms)
    if !has_conic_form(unique_conic_forms, x)
        objective = conic_form!(x.children[1], unique_conic_forms)

        sz = length(x)

        num_rows = x.size[1]
        num_cols = x.size[2]

        I = Array{Int}(undef, sz)
        J = Array{Int}(undef, sz)

        k = 1
        for r in 1:num_rows
            for c in 1:num_cols
                I[k] = (c - 1) * num_rows + r
                J[k] = (r - 1) * num_cols + c
                k += 1
            end
        end

        transpose_matrix = sparse(I, J, 1.0)

        objective = transpose_matrix * objective
        cache_conic_form!(unique_conic_forms, x, objective)
    end
    return get_conic_form(unique_conic_forms, x)
end

transpose(x::AbstractExpr) = TransposeAtom(x)

struct AdjointAtom <: AbstractExpr
    head::Symbol
    id_hash::UInt64
    children::Tuple{AbstractExpr}
    size::Tuple{Int,Int}

    function AdjointAtom(x::AbstractExpr)
        children = (x,)
        return new(:adjoint, hash(children), children, (x.size[2], x.size[1]))
    end
end

function sign(x::AdjointAtom)
    return sign(x.children[1])
end

function monotonicity(x::AdjointAtom)
    return (Nondecreasing(),)
end

function curvature(x::AdjointAtom)
    return ConstVexity()
end

function evaluate(x::AdjointAtom)
    return evaluate(x.children[1])'
end

# Since everything is vectorized, we simply need to multiply x by a permutation
# matrix such that coeff * vectorized(x) - vectorized(x') = 0
function conic_form!(x::AdjointAtom, unique_conic_forms::UniqueConicForms)
    if !has_conic_form(unique_conic_forms, x)
        objective = conic_form!(x.children[1], unique_conic_forms)

        sz = length(x)

        num_rows = x.size[1]
        num_cols = x.size[2]

        I = Array{Int}(undef, sz)
        J = Array{Int}(undef, sz)

        k = 1
        for r in 1:num_rows
            for c in 1:num_cols
                I[k] = (c - 1) * num_rows + r
                J[k] = (r - 1) * num_cols + c
                k += 1
            end
        end

        transpose_matrix = sparse(I, J, 1.0)
        objective = transpose_matrix * objective

        for var in keys(objective)
            x1 = conj(objective[var][1])
            x2 = conj(objective[var][2])
            objective[var] = (x1, x2)
        end
        cache_conic_form!(unique_conic_forms, x, objective)
    end
    return get_conic_form(unique_conic_forms, x)
end

adjoint(x::AbstractExpr) = AdjointAtom(x)
adjoint(x::Constant) = Constant(x.value')


text:
module Bombe

    using ParetoSmooth
    using Statistics
    using StructArrays
    using Tullio
    using VectorizationBase

    include("Loo.jl")
    using .LeaveOneOut
    export psis_loo
end

max_stars_repo_name: ParadaCarleton/Bombe

max_issues_repo_name: ParadaCarleton/Bombe

max_forks_repo_name: ParadaCarleton/Bombe


text:
using OrdinaryDiffEq, ModelingToolkit, CardiovascularModels

@variables t
D = Differential(t)

#left Ventricle
@named l_ventricle = Ventricle(Ees = 200e6,
    Vd = 0., V0 = 0., λ = 33e3, P0 = 10.)

@named card_driver = automaticDriver()

@named aortic_valve = Valve(R = 6e6)

# Aorta
@named aorta = ElasticVessel(Ees = 200e6, Vd = 100e-6)

@named systemic_resistance = Resistor(R = 50e6)

eqs_driver = [
    l_ventricle.drv ~ card_driver.contraction
]

@named vein = ElasticVessel(Ees = 5e6, Vd = 500e-6)
@named mitral_valve = Valve(R = 10e6)

eqs_con = serial_connect(vein,
    mitral_valve,
    l_ventricle,
    aortic_valve,
    aorta,
    systemic_resistance,
    vein)

eqs_comb = [eqs_driver; eqs_con]

# ODE parameters
volume_start = [l_ventricle.V => 1e-4,
    aorta.V => 1.5e-4,
    vein.V => 2e-3]
time_span = (0.0, 10.0)

@named hemo_sys = ODESystem(eqs_comb)

@named connected = compose(hemo_sys, l_ventricle, card_driver,
                            aortic_valve, aorta, systemic_resistance, vein, mitral_valve)

problem = ODEProblem(structural_simplify(connected), volume_start, time_span, [])

sol = solve(problem, Tsit5(), dtmax = 0.01, reltol = 1e-6)

# Plots.plot(sol)


text:
include("Parser.jl")

struct StaticAnalysisException <: Exception
  msg::String
end

default_int_value = -1
default_bool_value = false
default_string_value = ""

@enum MPLType begin
  MInt
  MBool
  MString
end

struct SValue
  type::MPLType
  value::Union{Int,Bool,String}
end

function SValue(tc::TokenClass)
  if tc == kw_int return SValue(MInt) end
  if tc == kw_bool return SValue(MBool) end
  if tc == kw_string return SValue(MString) end
end

function SValue(type::MPLType)
  if type == MInt return SValue(type, default_int_value) end
  if type == MBool return SValue(type, default_bool_value) end
  if type == MString return SValue(type, default_string_value) end
end

(*)(left::SValue, right::SValue) = SValue(MInt, left.value * right.value)

function (+)(left::SValue, right::SValue)
  if left.type == MInt
    return SValue(MInt, left.value + right.value)
  elseif left.type == MString
    return SValue(MString, left.value * right.value)
  end
end

(-)(left::SValue, right::SValue) = SValue(MInt, left.value - right.value)
(÷)(left::SValue, right::SValue) = SValue(MInt, left.value ÷ right.value)
(==)(left::SValue, right::SValue) = SValue(MBool, left.value == right.value)
(<)(left::SValue, right::SValue) = SValue(left.type, left.value < right.value)
(&)(left::SValue, right::SValue) = SValue(left.type, left.value & right.value)
(!)(operand::SValue) = SValue(MBool, !operand.value)

operator_to_function = Dict(
  times => *,
  plus => +,
  minus => -,
  divide => ÷,
  equals => ==,
  log_and => &,
  log_not => !,
  less_than => <
)

unary_result_types = Dict(
  (MBool, !) => MBool
)

binary_result_types = Dict(
  (MBool, &) => MBool,
  (MInt, *) => MInt,
  (MInt, +) => MInt,
  (MInt, -) => MInt,
  (MInt, ÷) => MInt,
  (MInt, <) => MBool,
  (MBool, <) => MBool,
  (MString, <) => MBool,
  (MInt, ==) => MBool,
  (MBool, ==) => MBool,
  (MString, ==) => MBool,
  (MInt, ÷) => MInt,
  (MString, +) => MString,
)

"""
A SymbolTable contains a Dict where the keys are variable names and
the values are (<variable value>, <variable value is locked>) pairs.
Iteration variables are locked during the execution of the for loop;
all other variables are unlocked.
"""
mutable struct SymbolTable
  table::Dict{String,Tuple{SValue, Bool}}
end

SymbolTable() = SymbolTable(Dict{String,Tuple{SValue, Bool}}())

function addOrUpdateSymbol(st::SymbolTable, symbol::String, sv::SValue, line=0::Int)
  if hasSymbol(st, symbol) && st.table[symbol][2]
    throw(StaticAnalysisException(
      "Attempting to update iteration variable $(symbol) on line $(line)."
    ))
  end
  push!(st.table, symbol => (sv, false))
end

function hasSymbol(st::SymbolTable, symbol::String)
  return symbol ∈ keys(st.table)
end

function getValue(st::SymbolTable, symbol::String)
  return st.table[symbol][1]
end

function setIterationVariable(st::SymbolTable, symbol::String, value::Int)
  push!(st.table, symbol => (SValue(MInt, value), true))
end

# Unlock an iteration variable upon for loop completion.
function releaseIterationVariable(st::SymbolTable, symbol::String)
  value = st.table[symbol][1]
  push!(st.table, symbol => (value, false))
end

# The entry point for the static analyzer.
function staticAnalysis(AST::Statements)
  st = SymbolTable()
  for stmt in AST.statements
    staticAnalysis(stmt, st)
  end
end

function staticAnalysis(r::Read, st::SymbolTable)
  varName = r.variable.lexeme
  if !hasSymbol(st, varName)
    throw(StaticAnalysisException(
      "Attempting to read into an unassigned variable on line $(r.line)."
    ))
  end
end

function staticAnalysis(a::Assert, st::SymbolTable)
  if typeCheck(a.argument, st) != MBool
    println("Asserting a non-boolean value on line $(a.line).")
  end
end

function staticAnalysis(d::Declaration, st::SymbolTable)
  varName = d.variable.lexeme
  if hasSymbol(st, varName)
    throw(StaticAnalysisException(
      "Attempting to declare existing variable $(varName) on line $(d.line)."
    ))
  else
    addOrUpdateSymbol(st, varName, SValue(d.type.class), d.line)
  end
end

function staticAnalysis(d::DecAssignment, st::SymbolTable)
  varName = d.variable.lexeme
  if hasSymbol(st, varName)
    throw(StaticAnalysisException(
      "Attempting to declare existing variable $(varName) on line $(d.line)."
    ))
  else
    addOrUpdateSymbol(st, varName, SValue(typeCheck(d.value, st)), d.line)
  end
end

function staticAnalysis(f::For, st::SymbolTable)
  varName = f.variable.lexeme
  if (hasSymbol(st, varName) && getValue(st, varName).type != MInt)
    throw(StaticAnalysisException(
      "Iteration variable on line $(f.line) has non-integer type."
    ))
  end
  if (typeCheck(f.from, st) != MInt)
    throw(StaticAnalysisException(
      "Range start value on line $(f.line) has non-integer type."
    ))
  end
  if (typeCheck(f.to, st) != MInt)
    throw(StaticAnalysisException(
      "Range end value on line $(f.line) has non-integer type."
    ))
  end
  setIterationVariable(st, varName, 1)
  for stmt in f.body.statements
    staticAnalysis(stmt, st)
  end
  releaseIterationVariable(st, varName)
end

function staticAnalysis(p::Print, st::SymbolTable)
  typeCheck(p.argument, st)
end

function staticAnalysis(a::Assignment, st::SymbolTable)
  varName = a.variable.lexeme
  if !hasSymbol(st, varName)
    throw(StaticAnalysisException(
      "Attempting to assign a value to an undeclared variable on line $(a.line)."
    ))
  end
  variableType = getValue(st, varName).type
  valueType = typeCheck(a.value, st)
  if variableType != valueType
    throw(StaticAnalysisException(
      """Attempting to assign a value of type $(valueType)
      to a variable of type $(variableType) on line $(a.line)."""
    ))
  end
  addOrUpdateSymbol(st, varName, SValue(valueType), a.line)
end

function typeCheck(l::Literal, st::SymbolTable)
  class = l.token.class
  class == int_literal && return MInt
  class ∈ [kw_true, kw_false] && return MBool
  class == string_literal && return MString
end

function typeCheck(u::UnaryOperation, st::SymbolTable)
  operandType = typeCheck(u.operand, st)
  operatorFunction = operator_to_function[u.operator.class]
  return unary_result_types[(operandType, operatorFunction)]
end

function typeCheck(b::BinaryOperation, st::SymbolTable)
  leftType = typeCheck(b.leftOperand, st)
  rightType = typeCheck(b.rightOperand, st)
  if leftType != rightType
    throw(StaticAnalysisException(
      "Operand type mismatch on line $(b.line)"
    ))
  end
  operatorFunction = operator_to_function[b.operator.class]
  if (leftType, operatorFunction) ∈ keys(binary_result_types)
    return binary_result_types[(leftType, operatorFunction)]
  else
    throw(StaticAnalysisException(
      "Operator '$(b.operator.lexeme)' on line $(b.line) is not valid for provided argument types."
    ))
  end
end

function typeCheck(v::Var, st::SymbolTable)
  varName = v.variable.lexeme
  if !hasSymbol(st, varName)
    throw(StaticAnalysisException(
      "Variable '$(varName)' on line $(v.line) has not been declared."
    ))
  end
  return getValue(st, varName).type
end

max_stars_repo_name: juhamyllari/mini-pl-julia

max_issues_repo_name: juhamyllari/mini-pl-julia

max_forks_repo_name: juhamyllari/mini-pl-julia


text:
""" FS::CholeskySpecialShifted{T,MT}

Special Cholesky type to allow for "removal" of rows and columns from a Cholesky factorization.
Factorizes the shifted (positive definite A+ϵI) and used iterative refinement to solve system.

`F0 = cholesky(A+ϵI)`
`F = CholeskySpecialShifted(F0, A, ϵ)`
or
`F = CholeskySpecialShifted(A, shift=...)`

Will destroy F0!

"Removing" row and column generates the factorization corresponding to
`A[:,i] = 0; A[i,:] = 0; A[i,i] = 1`

`deleterowcol!(F::CholeskySpecialShifted, i)`
    Delete row and column i

`addrowcol!(F::CholeskySpecialShifted, i)`
    Reverse the deletion row and column i

Supports `ldiv!`

Requires roughly n^2 + 3n extra memory
Update of row k requires O(n^2) computations
Solving is equivalent to a normal cholesky

"""
struct CholeskySpecialShifted{T,MT} <: AbstractCholeskySpecial{T,MT}
    F::Cholesky{T,MT}
    idx::Set{Int}               # Indices of deleted rows and column
    tmp::Vector{T}              # tmp storage in deleterowcol
    tmp2::Vector{T}             # tmp storage in deleterowcol
    tmp3::Vector{T}             # tmp storage in ldiv2!
    tmp4::Vector{T}             # tmp storage in ldiv2!
    tmp5::Vector{T}             # tmp storage in ldiv2!
    M::MT
    shift::T
end

# TODO We should try to make some deafult choise here
function CholeskySpecialShifted(F::Cholesky{T,MT}, M, shift) where {T,MT}
    if F.uplo != 'U'
        error("Only implemented for U type")
    end
    n = size(F,1)
    CholeskySpecialShifted{T,MT}(F, Set{Int}(),                 # F and idx
        fill(zero(T), n), fill(zero(T), n), fill(zero(T), n),   # tmp, tmp2, tmp3
        fill(zero(T), n), fill(zero(T), n),                     # tmp4, tmp5
        M, shift)                                               # M (unshifted), shift
end

function CholeskySpecialShifted(M::AbstractMatrix{T}, shift=sqrt(sqrt(eps(T))*eps(T))) where T
    F = cholesky(M + I*shift)
    return CholeskySpecialShifted(F, M, shift)
end

# Helper function norm(a-b)
function normdiff(a,b)
    n = length(a)
    @assert n == length(b)
    s = zero(eltype(a))
    @inbounds @simd for i in eachindex(a)
        s += (a[i]-b[i])^2
    end
    return sqrt(s)
end

"""
b, projection = ldiv2!(F::CholeskySpecialShifted{T,MT}, b::AbstractVector{T}; x0=zero(T))

Solve `Mx=b` where F is Cholesky factorization if M with some rows and columns set to identity

The solution `x` is the projection onto `Mx=b` from `x0`.
`x0` can be a Vector or Number (treated as `fill(x0,n)`).

If no solution exists, finds projection of `b` onto space `Mx=0`.

`projection=false` if soulution `Mx=b` found, and
`projection=true` if not.

Overwrites and returns `b`.

"""
function ldiv2!(F::CholeskySpecialShifted{T,MT}, b::AbstractVector{T}; x0=zero(T), data=nothing) where {T,MT}
    done = false
    nleft = 2 # Number of iterations to run when done
    projection = false
    rk = F.tmp
    xk = F.tmp2
    rkold = F.tmp3
    xkold = F.tmp4
    xkoldold = F.tmp5

    rnorm = typemax(norm(zero(T)/one(T))) # Inf of right type
    xnorm = rnorm   # norm(x_k)
    drnorm = rnorm  # norm(r_k - r_(k-1))
    dxnorm = rnorm  # norm(x_k - x_(k-1))
    ddxnorm = rnorm # norm(x_k -2x_(k+1) + x_(k-2))

    # Initial value
    xk .= x0
    # Set rk (xkdiffold) to Inf
    rk .= rnorm

    for j in F.idx  # M should be M with some cols identity
        xk[j] = b[j] # This might ruin projection?
    end
    # Iterative refinement
    for i = 1:100
        for j in F.idx  # let xk be zero where M shouldn't act
            xk[j] = zero(T)
        end
        mul!(rk, F.M, xk)   # rk = M*xk
        rk .= b .- rk       # rk = b - M*xk
        for j in F.idx  # M should be M with some cols identity
            rk[j] = zero(T) # zero error here, This might ruing projection?
        end
        rnorm = norm(rk)
        drnorm = normdiff(rk, rkold)
        # Now save rk
        rkold .= rk
        #push!(data.rk, copy(rk))
        ldiv!(F.F, rk)      # rk = (M+ϵI)⁻¹(b - M*xk)
        for j in F.idx  # M should be M with some cols identity
            rk[j] = zero(T) # zero error here, This might ruing projection?
        end
        xk .+= rk           # xk .= xk + (M+ϵI)⁻¹(b - M*xk)

        #push!(data.xk, copy(xk))

        xnorm = norm(xk)
        dxnorm = norm(rk)
        # We can now use xkoldold as temp storage for ddx
        xkoldold .= rk .- xkold .+ xkoldold # ddx
        #rk .= rk .- xkold .+ xkoldold # ddx
        ddxnorm = norm(xkoldold)
        # Now save xk
        xkoldold .= xkold
        xkold .= xk

        # DEBUG && println("rnorm: $rnorm")
        # DEBUG && println("xnorm: $xnorm")
        # DEBUG && println("drnorm: $drnorm")
        # DEBUG && println("dxnorm: $(dxnorm)")
        # println("ddxnorm: $(ddxnorm)")
        # println("dxnorm/xnorm: $(dxnorm/xnorm)")
        # println("ddxnorm/xnorm: $(ddxnorm/xnorm)")

        # This should re robust
        # TODO 1e-11 for robustness?
        if rnorm < 1e-10 && drnorm < 1e-10 # Break if we found Mx=b solution

            done = true
            projection = false
            if nleft < 1
                DEBUG && println("Projection: $projection i: $i")
                break
            end
            nleft -= 1
        end
        # Maybe not as robust, so only use one case
        if ddxnorm/xnorm < 1e-7
            if dxnorm/xnorm < 1e-2 # Let the robust case catch this
                # done = true
                # projection = false
                # break
            else
                done = true
                projection = true
                if nleft < 1
                    DEBUG && println("Projection: $projection i: $i")
                    break
                end
                nleft -= 1
            end
        end
    end
    for j in F.idx  # M should be M with some cols identity
        xk[j] = b[j] # This might ruin projection?
    end

    if !done # Allow larger erros
        DEBUG && println("No convergence")
        if rnorm < 1e-10 # Assume solution
            if dxnorm > 1e-10 # We don't want this
                error("ldiv! did not converge to a solution, case 1, rnorm: $rnorm, dxnorm: $(dxnorm)")
            end
            projection = false
            DEBUG && println("Projection: $projection after")
        elseif ddxnorm/xnorm < 1e-5 && dxnorm/xnorm > 1e-3
            projection = true
            DEBUG && println("Projection: $projection after")
        else
            # Try to catch case where it doesnt really matter
            #@warn "Testing direction we are not sure about"
            #projection = true
            error("ldiv! did not converge to a solution, case 2, rnorm: $rnorm, dxnorm: $(dxnorm), dxnorm/xnorm: $(dxnorm/xnorm), ddxnorm/xnorm: $(ddxnorm/xnorm)")
        end
    end
    DEBUG && println("Projection: $projection at end test")
    if projection
        b .= rk.*F.shift
    else
        b .= xk
    end

    return b, projection
end

"""

ldiv!(F::CholeskySpecialShifted{T,MT}, b::AbstractVector{T}; x0=zero(T))

Solve Mx=b where F is Cholesky factorization if M with some rows and columns set to identity

"""
function LinearAlgebra.ldiv!(F::CholeskySpecialShifted{T,MT}, b::AbstractVector{T}; x0=zero(T)) where {T,MT}
    _, projection = ldiv2!(F, b, x0=x0)
    if projection
        @error "ldiv! did not converge to a solution"
        error(SingularException)
    end
    return b
end

Base.size(F::CholeskySpecialShifted) = size(F.F)
Base.size(F::CholeskySpecialShifted, i::Int) = size(F.F)[i]

function Base.getproperty(F::CholeskySpecialShifted{T,MT}, s::Symbol) where {T,MT}
    if s == :U
        return F.F.U
    else
        getfield(F,s)
    end
end

""" `deleterowcol!(F::CholeskySpecialShifted{T}, i)`
    Given F=U*U', update factorization corresponding to setting
    F[i, :] = 0, F[:, i] = 0
    F[i,i] = 1.
"""
function deleterowcol!(F::CholeskySpecialShifted{T,MT}, i) where {T,MT}
    if i in F.idx
        error("Can not delete row and column that is already deleted")
    end
    n = size(F.F,1)
    # Create vector for rank update
    F.tmp[1:i] .= zero(T)
    F.tmp[(i+1):n] .= F.F.U[i, (i+1):n]
    # Record that  this row/column is "removed"
    push!(F.idx, i)
    # Do rank update
    lowrankupdate!(F.F, F.tmp)
    # Zero out
    for i in F.idx
        F.F.U[1:(i-1),i] .= zero(T)
        F.F.U[i,(i+1):end] .= zero(T)
        F.F.U[i,i] = one(T)
    end
    return
end

"""  Reverse deleterowcol!
"""
function addrowcol!(F::CholeskySpecialShifted{T,MT}, i) where {T,MT}
    if !(i in F.idx)
        error("Can not add row and column that has not been deleted")
    end
    n = size(F.F,1)
    b = F.tmp
    b .= F.M[:,i]
    b[i] += F.shift # Add shift from [i,i]
    b[i:end] .= zero(T)
    for j in F.idx
        #This column should be zero in M
        b[j] = zero(T)
    end
    ldiv!(F.F.U', b) # b is now [S12; junk]
    S12 = view(b, 1:(i-1))
    # Scalar result
    S22 = sqrt((F.M[i,i] + F.shift) - S12'S12)

    # Get row from M (under assumption that some rows/cols are identity)
    M23 = view(F.tmp2, (i+1):n)
    M23 .= F.M[i,(i+1):end] # No shift here
    for j in F.idx
        if j > i
            #This column should be zero in M
            F.tmp2[j] = zero(T)
            # Equivalent to
            # M23[j-i] = zero(T)
        end
    end
    # TODO see if view on factor matrix directly is faster
    # Now set the new row to
    # S23 = (M23 - F.U[1:(i-1),(i+1):end]'S12)/S22 # Transposed
    S23 = view(F.F.U, i, (i+1):n) # This is where we write
    S13 = view(F.F.U, 1:(i-1), (i+1):n)
    mul!(S23, S13', S12)
    S23 .= (M23 .- S23)./S22

    # Set the diagonal
    F.F.U[i,i] = S22
    # And the column
    F.F.U[1:(i-1),i] .= S12

    # No longer need S12, so reuse tmp for rank update
    F.tmp[1:i] .= zero(T)
    F.tmp[(i+1):n] .= S23

    #Update the triangular part S33
    lowrankdowndate!(F.F, F.tmp)

    # Rememeber that i is added back
    pop!(F.idx, i)
    return
end

max_stars_repo_name: mfalt/QPDAS

max_issues_repo_name: mfalt/QPDAS

max_forks_repo_name: mfalt/QPDAS


text:
export
    Packed_Normal,
    Packed_MvNormal,
    Packed_AliasingScalarSampler,
    Packed_BallTreeDensity

mutable struct Packed_Normal
  mean::Float64
  std::Float64
  distType::String
end

mutable struct Packed_MvNormal
  mean::Vector{Float64}
  cov::Vector{Float64}
  distType::String
end

mutable struct Packed_AliasingScalarSampler
  samples::Vector{Float64}
  weights::Vector{Float64}
  quantile::Union{Float64, Nothing}
  distType::String # AliasingScalarSampler
end

mutable struct Packed_BallTreeDensity
  dim::Int  # point dimension
  kernelType::String
  bandwidth::Vector{Float64}
  weights::Vector{Float64}
  points::Vector{Float64} # Stack points as a 1D vector.
end


text:
# bello() = runi(`brishz.dash redo2 2 bell-greencase`, wait=false)
bello() = runi(`brishz.dash awaysh bello`, wait=false)
export bello

bellj() = runi(`brishz.dash awaysh bellj`, wait=false)
export bellj

function bella(host=nothing)
    if host === nothing ||
        (host == :mbp && mbpP())
        runi(`brishz.dash awaysh bella`, wait=false)
    end
end
export bella

okj() = runi(`brishz.dash awaysh okj`, wait=false)
export okj

function firstbell()
    if ! @isdefined firstLoad
        bellj()
    else
        bello()
    end
    global firstLoad = true
end
export firstbell


text:
struct BSONIndexKey
    name_p::Ptr{UInt8}
    name_len::Int32
    parent_offset::Int32
end

@inline BSONIndexKey(name::AbstractString, parent_offset::Integer) = BSONIndexKey(
    pointer(name),
    sizeof(name) % Int32,
    parent_offset % Int32
)

@inline BSONIndexKey(name::Symbol, parent_offset::Integer) = BSONIndexKey(
    Base.unsafe_convert(Ptr{UInt8}, name),
    ccall(:strlen, Csize_t, (Cstring,), Base.unsafe_convert(Ptr{UInt8}, name)) % Int32,
    parent_offset % Int32
)

@inline function Base.isequal(x::BSONIndexKey, y::BSONIndexKey)
    x.name_len == y.name_len &&
        x.parent_offset == y.parent_offset &&
        ccall(:memcmp, Cint, (Ptr{UInt8}, Ptr{UInt8}, Csize_t), x.name_p, y.name_p, x.name_len % Csize_t) == 0
end

struct BSONIndexValue
    offset::Int32
    type::UInt8
end

mutable struct BSONIndex
    entries::Vector{Tuple{BSONIndexKey, Int32, BSONIndexValue}}
    version::Int32
    size_mask::Int32
    include_arrays::Bool

    function BSONIndex(size::Integer; include_arrays::Bool = false)
        pow2size = nextpow(2, size)
        new(
            fill((BSONIndexKey(Ptr{UInt8}(0), 0, 0), 0, BSONIndexValue(0, 0)), pow2size),
            0,
            pow2size - 1,
            include_arrays
        )
    end
end

function build_index_(index::BSONIndex, reader::BSONReader, ::Val{include_arrays}) where include_arrays
    let parent_offset = reader.offset % Int32
        f = @inline x -> begin
            name = x.first
            field_reader = x.second
            key = BSONIndexKey(name.ptr, name.len % Int32, parent_offset)
            value = BSONIndexValue(field_reader.offset % Int32, field_reader.type)
            index[key] = value
            if field_reader.type == BSON_TYPE_DOCUMENT || (include_arrays && field_reader.type == BSON_TYPE_ARRAY)
                build_index_(index, field_reader, Val{include_arrays}())
            end
            nothing
        end
        foreach(f, Map(identity), reader)
    end
    nothing
end

@inline function index_(index::BSONIndex, key::BSONIndexKey)
    nh = fnv1a(UInt32, key.name_p, key.name_len)
    oh = reinterpret(UInt32, key.parent_offset) * 0x9e3779b1
    h = nh ⊻ (oh + 0x9e3779b9 + (nh << 6) + (nh >> 2))
    (h & index.size_mask) % Int + 1
end

@inline function Base.setindex!(index::BSONIndex, reader::BSONReader)
    index.version += Int32(1)
    if index.include_arrays
        build_index_(index, reader, Val{true}())
    else
        build_index_(index, reader, Val{false}())
    end
    nothing
end

@inline function Base.setindex!(index::BSONIndex, value::BSONIndexValue, key::BSONIndexKey)
    @inbounds index.entries[index_(index, key)] = (key, index.version, value)
    nothing
end

@inline function Base.getindex(index::BSONIndex, key::BSONIndexKey)
    @inbounds e_key, e_version, e_value = index.entries[index_(index, key)]
    (e_version == index.version && isequal(e_key, key)) ? e_value : nothing
end


text:
@testset "▶ /layout                     " begin
    # subplot(abc)
    f = Figure()
    subplot(222)
    @test length(f.axes) == 4
    @test gca() == f.axes[2]
    @test_throws ArgumentError subplot(10010)

    # subplot(a, b, c)
    @test_throws ArgumentError subplot(10, 1, 1)
    @test_throws ArgumentError subplot(1, 10, 1)
    @test_throws ArgumentError subplot(2, 2, 10)
    @test_throws ArgumentError subplot(2, 3, 5)

    # layout
    clf()
    @test isempty(f.axes)
    W, H = gcf().size
    layout(f, [0.1 0.1 0.3 0.3;
               0.1 0.5 0.3 0.3])
    @test length(f.axes) == 2
    @test f.axes[1].origin == (0.1W, 0.1H)
    @test f.axes[2].origin == (0.1W, 0.5H)
    @test f.axes[1].size == (0.3W, 0.3H)
    @test f.axes[2].size == (0.3W, 0.3H)

    @test_throws ArgumentError layout(f, [0.1 0.1 0.3 0.3 0.0; 0.1 0.5 0.3 0.3 0.0])
    @test_throws ArgumentError layout(f, [0.1 5.1 0.3 0.3 0.0; 0.1 0.5 0.3 0.3 0.0])
end


text:
using ThreePhotons
using Compat

function check_git_status()
    status = readstring(`git ls-files --modified`)
    @assert length(status) == 0 "Please commit all changes before running server scripts."

    return readstring(`git rev-parse HEAD`)
end

function jobname(directory, number::Int64=0)
    main = replace(directory, "/", "_")
    suffix = number > 0 ? "_$(number)" : ""
    return "$main$suffix"
end

function environment_path(local_path::String)
    return "$(ENV["DETERMINATION_DATA"])/$(ENV_root)/$(local_path)"
end

if contains(readstring(`hostname`), "owl")
    include("environment_owl.jl")
elseif contains(readstring(`hostname`), "hydra")
    include("environment_hydra.jl")
elseif contains(readstring(`hostname`), "gwdu103")
    include("environment_gwdg.jl")
else
    include("environment_local.jl")
end

include("run_generate_histograms.jl")
include("run_optimal.jl")
include("run_parallel_determination.jl")


text:
# vim: set et ts=2 sw=2:
@with_kw mutable struct Some
  pos=0; txt=""; w=1; n=0;
  _all=[]; max=it.some.max; stale=false end

@with_kw mutable struct Sym
  pos=0; txt=""; w=1; n=0;
  seen=Dict(); mode=no; most=0 end

#---- Create -----------------------------------------------
function col(;txt="",pos=0, c=it.char)
  x= c.less in txt||c.more in txt||c.num in txt ? Some : Sym
  w= it.char.less in txt ? -1 : 1
  x(txt=txt, pos=pos, w=w) end

#---- Update -----------------------------------------------
inc!(i,x)=
  x==it.char.skip ? x : begin i.n += 1; inc1!(i,x); x end

function inc1!(i::Sym, x)
  new = i.seen[x] = 1 + get(i.seen,x,0)
  if new > i.most
    i.mode, i.most = x,new end end

function inc1!(i::Some, x::Number)
  m = length(i._all)
  if m < i.max
    i.stale=true
    push!(i._all,x)
  elseif rand() < m/i.n
    i.stale=true
    i._all[ int(m*rand())+1 ]=x end end

#---- Query ------------------------------------------------
norm!(i::Some,x, a=all(i), skip=x==it.char.skip) =
  skip ? x : max(0,min(1, (x-a[1]) / (a[end] - a[1]+1E-31)))

mid(i::Some;lo=no,hi=no) = per(i,p=.5,lo=lo,hi=hi)

sd(i::Some;lo=no,hi=no)  = (per(i,p=.9,lo=lo,hi=hi) -
                            per(i,p=.1,lo=lo,hi=hi)) / 2.564

function per(i::Some;p=.5,lo=no,hi=no, lst=all(i))
  hi = hi==no ? length(lst) : hi
  lo = lo==no ? 1           : lo
  lst[ int(lo + p*(hi - lo +1)) ] end

function all(i::Some)
  if i.stale i._all=sort(i._all) end
  i.stale=false
  i._all end

max_stars_repo_name: timm/jut

max_issues_repo_name: timm/jut

max_forks_repo_name: timm/jut


text:
function varsassign(mystr::String, varname::String, n::Int, flg::Int)

if flg==0
    for i = 1:n
        str1 = varname * string(i);
        str2 = varname * "[" * string(i-1) * "]";
        mystr = mystr * "\t\tT " * str1 * " = " * str2 * ";\n";
    end
else
    for i = 1:n
        str1 = varname * string(i);
        str2 = varname * "[" * string(i-1) * "*ng+i" * "]";
        mystr = mystr * "\t\tT " * str1 * " = " * str2 * ";\n";
    end
end

return mystr;

end

max_stars_repo_name: glwagner/Exasim

max_issues_repo_name: glwagner/Exasim

max_forks_repo_name: glwagner/Exasim


text:
struct TSIterator{T, S}
    source::S
end

IteratorInterfaceExtensions.isiterable(x::Temporal.TS) = true
TableTraits.isiterabletable(x::Temporal.TS) = true

function IteratorInterfaceExtensions.getiterator(ta::S) where {S<:Temporal.TS}
    etype = eltype(ta.values)

    T = NamedTuple{(:Index, ta.fields...), Tuple{eltype(ta.index), fill(etype, length(ta.fields))...}}

    return TSIterator{T,S}(ta)
end

function Base.length(iter::TSIterator)
    return length(iter.source.index)
end

function Base.eltype(iter::TSIterator{T,TS}) where {T,TS}
    return T
end

Base.eltype(::Type{TSIterator{T,TS}}) where {T,TS} = T

@generated function Base.iterate(iter::TSIterator{T,TS}, state=1) where {T,TS}
    constructor_call = Expr(:call, :($T))

    push!(constructor_call.args, Expr(:tuple))

    # Add timestamp column
    push!(constructor_call.args[2].args, :(iter.source.index[i]))

    for i in 1:fieldcount(T)-1
        push!(constructor_call.args[2].args, :(iter.source.values[i,$i]))
    end

    quote
        if state>length(iter.source.index)
            return nothing
        else
            i = state
            a = $constructor_call
            return a, state+1
        end
    end
end

# Sink

function Temporal.TS(x; index_column::Symbol=:Index)
    TableTraits.isiterabletable(x) || error("Cannot create a TS from something that is not a table.")

    iter = getiterator(x)

    et = eltype(iter)

    if fieldcount(et)<2
        error("Need at least two columns")
    end

    names = fieldnames(et)

    timestep_col_index = findfirst(isequal(index_column), names)

    if timestep_col_index===nothing
        error("No timestamp column found.")
    end

    timestep_col_index = something(timestep_col_index)

    col_types = [fieldtype(et, i) for i=1:fieldcount(et)]

    data_columns = collect(Iterators.filter(i->i[2][1]!=index_column, enumerate(zip(names, col_types))))

    orig_data_type = data_columns[1][2][2]

    data_type = orig_data_type <: DataValues.DataValue ? orig_data_type.parameters[1] : orig_data_type

    orig_timestep_type = col_types[timestep_col_index]

    timestep_type = orig_timestep_type <: DataValues.DataValue ? orig_timestep_type.parameters[1] : orig_timestep_type

    if any(i->i[2][2]!=orig_data_type, data_columns)
        error("All data columns need to be of the same type.")
    end

    t_column = Vector{timestep_type}(undef,0)
    d_array = Vector{Vector{data_type}}(undef,0)
    for i in data_columns
        push!(d_array, Vector{data_type}(undef,0))
    end

    for v in iter
        if orig_timestep_type <: DataValue
            push!(t_column, get(v[timestep_col_index]))
        else
            push!(t_column, v[timestep_col_index])
        end

        if orig_data_type <: DataValue
            for (i,c) in enumerate(data_columns)
                push!(d_array[i],get(v[c[1]]))
            end
        else
            for (i,c) in enumerate(data_columns)
                push!(d_array[i],v[c[1]])
            end
        end
    end

    d_array = hcat(d_array...)

    ta = Temporal.TS(d_array, t_column,[i[2][1] for i in data_columns])
    return ta
end