text:
module ClearStacktrace

const MODULECOLORS = [:light_blue, :light_yellow, :light_red, :light_green, :light_magenta, :light_cyan, :blue, :yellow, :red, :green, :magenta, :cyan]

const EXPAND_BASE_PATHS = Ref(true)
const CONTRACT_USER_DIR = Ref(true)
const REPLACE_BACKSLASHES = Ref(true)

function expandbasepath(str)

    basefileregex = if Sys.iswindows()
        r"^\.\\\w+\.jl$"
    else
        r"^\./\w+\.jl$"
    end

    if !isnothing(match(basefileregex, str))
        sourcestring = Base.find_source_file(str[3:end]) # cut off ./
    else
        str
    end
end

function replaceuserpath(str)
    str1 = replace(str, homedir() => "~")
    # seems to be necessary for some paths with small letter drive c:// etc
    replace(str1, lowercasefirst(homedir()) => "~")
end

function replacebackslashes(str)
    replace(str, raw"\\" => '/')
end

getline(frame) = frame.line
function getfile(frame)
    file = string(frame.file)
    if EXPAND_BASE_PATHS[]
        file = expandbasepath(file)
    end
    if CONTRACT_USER_DIR[]
        file = replaceuserpath(file)
    end
    if REPLACE_BACKSLASHES[]
        file = replacebackslashes(file)
    end

    file
end
getfunc(frame) = string(frame.func)
getmodule(frame) = try; string(frame.linfo.def.module) catch; "" end

getsigtypes(frame) = try;  frame.linfo.specTypes.parameters[2:end] catch; "" end

function convert_trace(trace)
    files = getfile.(trace)
    lines = getline.(trace)

    methodss = map(trace) do t
        try
            t.linfo.def
        catch
            nothing
        end
    end

    arguments = map(methodss) do m
        if isnothing(m)
            []
        else
            tv, decls, file, line = Base.arg_decl_parts(m)
            decls[2:end]
        end
    end

    funcs = getfunc.(trace)
    moduls = getmodule.(trace)
    sigtypes = getsigtypes.(trace)
    inlineds = getfield.(trace, :inlined)

    # replace empty modules if there is another frame from the same file
    for (i_this, mo) in enumerate(moduls)
        if mo == ""
            for i_other in 1:length(trace)
                if files[i_this] == files[i_other] && moduls[i_other] != ""
                    moduls[i_this] = moduls[i_other]
                end
            end
        end
    end

    (files = files, lines = lines, funcs = funcs,
        moduls = moduls, sigtypes = sigtypes, inlineds = inlineds, arguments = arguments)
end

function printtrace(io::IO, converted_stacktrace)

    files, lines, funcs, moduls, sigtypes, inlineds, arguments = converted_stacktrace

    n = length(files)
    ndigits = length(digits(n))
    length_numstr = ndigits + 2

    uniquemodules = setdiff(unique(moduls), [""])
    modulecolors = Dict(u => c for (u, c) in
        Iterators.zip(uniquemodules, Iterators.cycle(MODULECOLORS)))

    for (i, (func, inlined, modul, file, line, stypes, args)) in enumerate(
            zip(funcs, inlineds, moduls, files, lines, sigtypes, arguments))

        modulecolor = get(modulecolors, modul, :default)
        print_frame(io, i, func, inlined, modul, file, line, stypes, args, length_numstr, modulecolor)
        println(io)
        println(io)
    end
end

function print_frame(io, i, func, inlined, modul, file, line, stypes,
        args, length_numstr, modulecolor)

    # frame number
    print(io, lpad("[" * string(i) * "]", length_numstr))
    print(io, " ")

    # function name
    printstyled(io, func, bold = true)

    # type signature
    printstyled(io, "(", color = :light_black)

    i = 1
    for (stype, (varname, vartype)) in zip(stypes, args)
        if i > 1
            printstyled(io, ", ", color = :light_black)
        end
        printstyled(io, string(varname), color = :light_black, bold = true)
        printstyled(io, "::")
        printstyled(io, string(stype), color = :light_black)
        i += 1
    end

    printstyled(io, ")", color = :light_black)

    println(io)

    # @
    printstyled(io, " " ^ (length_numstr - 1) * "@ ", color = :light_black)

    # module
    if !isempty(modul)
        printstyled(io, modul, color = modulecolor)
        print(io, " ")
    end

    # filepath
    pathparts = splitpath(file)
    for p in pathparts[1:end-1]
        printstyled(io, p * "/", color = :light_black)
    end

    # filename, separator, line
    # bright black (90) and underlined (4)
    print(io, "\033[90;4m$(pathparts[end] * ":" * string(line))\033[0m")

    # inlined
    printstyled(io, inlined ? " [inlined]" : "", color = :light_black)
end

@warn "Overloading Base.show_backtrace(io::IO, t::Vector) with custom version"
function Base.show_backtrace(io::IO, t::Vector)

    ### this part is copied from the original function
    resize!(Base.LAST_SHOWN_LINE_INFOS, 0)
    filtered = Base.process_backtrace(t)
    isempty(filtered) && return

    if length(filtered) == 1 && StackTraces.is_top_level_frame(filtered[1][1])
        f = filtered[1][1]
        if f.line == 0 && f.file == Symbol("")
            # don't show a single top-level frame with no location info
            return
        end
    end
    ###

    if length(filtered) > Base.BIG_STACKTRACE_SIZE
        Base.show_reduced_backtrace(IOContext(io, :backtrace => true), filtered, true)
        return
    end

    println(io, "\nStacktrace:")

    # process_backtrace returns a Tuple{Frame, Int}
    frames = first.(filtered)

    converted_stacktrace = convert_trace(frames)

    printtrace(io, converted_stacktrace)
end

end # module


text:
import Base: +, -, *, sum, getindex, circshift
import LinearAlgebra: dot, adjoint

# rules for creating the computation graph

for op in (:+, :-, :*, :dot)
    @eval begin
        $op(a::TrackedValue, b::TrackedValue) = add_call!($op, Any[a, b])
        $op(a::TrackedValue, @nospecialize(b)) = add_call!($op, Any[a, b])
        $op(@nospecialize(a), b::TrackedValue) = add_call!($op, Any[a, b])
    end
end

for op in (:+, :-, :sum)
    @eval begin
        $op(x::TrackedValue) = add_call!($op, Any[x])
    end
end

import Base.Broadcast: broadcasted, materialize
import Base: map

broadcasted(f, x::TrackedValue, y...) = add_call!(broadcasted, Any[f, x, y...])

materialize(x::TrackedValue) = add_call!(materialize, Any[x])

map(f, x::TrackedValue, y...) = add_call!(map, Any[f, x, y...])

getindex(x::TrackedValue, i) = add_call!(getindex, Any[x, i])

adjoint(x::TrackedValue) = add_call!(adjoint, Any[x])

circshift(x::TrackedValue, i) = add_call!(circshift, Any[x, i])

eval_tracked(x, y) = x(y)
eval_tracked(x::TrackedValue, y) = add_call!(eval_tracked, Any[x, y])
(x::TrackedValue)(y) = eval_tracked(x, y)

# rules for differentiation

diff_rule(::typeof(getindex), ::Val{1}, x, i) = y -> y[i]
diff_rule(::typeof(eval_tracked), ::Val{1}, x, i) = y -> y(i)

diff_rule(::typeof(adjoint), ::Val{1}, x) = y -> adjoint(y)
diff_rule(::typeof(circshift), ::Val{1}, x, i) = y -> circshift(y, i)

diff_rule(::typeof(+), ::Val{1}, x1, x2) = identity
diff_rule(::typeof(+), ::Val{2}, x1, x2) = identity

diff_rule(::typeof(-), ::Val{1}, x1, x2) = identity
diff_rule(::typeof(-), ::Val{2}, x1, x2) = y -> -y

diff_rule(::typeof(*), ::Val{1}, x1, x2) = y -> y * x2
diff_rule(::typeof(*), ::Val{2}, x1, x2) = y -> x1 * y

diff_rule(::typeof(dot), ::Val{1}, x1, x2) = y -> dot(y, x2)
diff_rule(::typeof(dot), ::Val{2}, x1, x2) = y -> dot(x1, y)

diff_rule(::typeof(sum), ::Val{1}, x) = y -> sum(y)

function partial_derivative(f, x::Tuple, ::Val{i}) where i
    return ForwardDiff.derivative(t -> f(Base.setindex(x, t, i)...), x[i])
end

function diff_rule(::typeof(broadcasted), ::Val{i}, f, x...) where i
    s = Base.broadcasted((x_...,) -> partial_derivative(f, x_, Val(i-1)), x...)
    s = Base.materialize(s)

    return y -> s .* y
end

function diff_rule(::typeof(map), ::Val{i}, f, x...) where i
    i_minus_1 = Val(i-1)
    s = map((x_...,) -> partial_derivative(f, x_, i_minus_1), x...)

    return y -> map(*, s, y)
end

diff_rule(::typeof(broadcasted), ::Val{2}, ::typeof(*), x1, x2) = y -> y .* x2
diff_rule(::typeof(broadcasted), ::Val{3}, ::typeof(*), x1, x2) = y -> x1 .* y

diff_rule(::typeof(materialize), ::Val{1}, x) = identity

ext: jl


text:
#=

# Minimal working compute example

The amount of control offered by Vulkan is not a very welcome property for
users who just want to run a simple shader to compute something quickly, and
the effort required for the "first good run" is often quite deterrent. To ease
the struggle, this tutorial gives precisely the small "bootstrap" piece of code
that should allow you to quickly run a compute shader on actual data. In short,
we walk through the following steps:

- Opening a device and finding good queues families and memory types
- Allocating memory and buffers
- Compiling a shader program and filling up the structures necessary to run it:
  - specialization constants
  - push constants
  - descriptor sets and layouts
- Making a command buffer and submitting it to the queue, efficiently running
  the shader

## Initialization

=#

using SwiftShader_jll # hide
using Vulkan
@set_driver :SwiftShader # hide

instance = Instance(["VK_LAYER_KHRONOS_validation"], [])

# Take the first available physical device (you might check that it is an
# actual GPU, using [`get_physical_device_properties`](@ref)).
physical_device = first(unwrap(enumerate_physical_devices(instance)))

# At this point, we need to choose a queue family index to use. For this
# example, have a look at `vulkaninfo` command and pick the good queue manually
# from the list of `VkQueueFamilyProperties` -- you want one that has
# `QUEUE_COMPUTE` in the flags. In a production environment, you would use
# [`get_physical_device_queue_family_properties`](@ref) to find a good index.
qfam_idx = 0

# Create a device object and make a queue for our purposes.
device = Device(physical_device, [DeviceQueueCreateInfo(qfam_idx, [1.0])], [], [])

# ## Allocating the memory
#
# Same way, you need to find a good memory type. Again, you can find a good one
# using `vulkaninfo` or with [`get_physical_device_memory_properties`](@ref).
# For compute, you want something that is both at the device (contains
# `MEMORY_PROPERTY_DEVICE_LOCAL_BIT`) and visible from the host
# (`..._HOST_VISIBLE_BIT`).
memorytype_idx = 0

# Let's create some data. We will work with 100 flimsy floats.
data_items = 100
mem_size = sizeof(Float32) * data_items

# Allocate the memory of the correct type
mem = DeviceMemory(device, mem_size, memorytype_idx)

# Make a buffer that will be used to access the memory, and bind it to the
# memory. (Memory allocations may be quite demanding, it is therefore often
# better to allocate a single big chunk of memory, and create multiple buffers
# that view it as smaller arrays.)
buffer = Buffer(
    device,
    mem_size,
    BUFFER_USAGE_STORAGE_BUFFER_BIT,
    SHARING_MODE_EXCLUSIVE,
    [qfam_idx],
)

bind_buffer_memory(device, buffer, mem, 0)

# ## Uploading the data to the device
#
# First, map the memory and get a pointer to it.
memptr = unwrap(map_memory(device, mem, 0, mem_size))

# Here we make Julia to look at the mapped data as a vector of `Float32`s, so
# that we can access it easily:
data = unsafe_wrap(Vector{Float32}, convert(Ptr{Float32}, memptr), data_items, own = false);

# For now, let's just zero out all the data, and *flush* the changes to make
# sure the device can see the updated data. This is the simplest way to move
# array data to the device.
data .= 0
unwrap(flush_mapped_memory_ranges(device, [MappedMemoryRange(mem, 0, mem_size)]))
# The flushing is not required if you have verified that the memory is
# host-coherent (i.e., has `MEMORY_PROPERTY_HOST_COHERENT_BIT`).

# Eventually, you may need to allocate memory types that are not visible from
# host, because these provide better capacity and speed on the discrete GPUs.
# At that point, you may need to use the transfer queue and memory transfer
# commands to get the data from host-visible to GPU-local memory, using e.g.
# [`cmd_copy_buffer`](@ref).

# ## Compiling the shader
#
# Now we need to make a shader program. We will use `glslangValidator` packaged
# in a JLL to compile a GLSL program from a string into a spir-v bytecode,
# which is later passed to the GPU drivers.

shader_code = """
#version 430

layout(local_size_x_id = 0) in;
layout(local_size_x = 1, local_size_y = 1, local_size_z = 1) in;

layout(constant_id = 0) const uint blocksize = 1; // manual way to capture the specialization constants

layout(push_constant) uniform Params
{
    float val;
    uint n;
} params;

layout(std430, binding=0) buffer databuf
{
    float data[];
};

void
main()
{
    uint i = gl_GlobalInvocationID.x;
    if(i < params.n) data[i] = params.val * i;
}
"""

# Push constants are small packs of variables that are used to quickly send
# configuration data to the shader runs. Make sure that this structure
# corresponds to what is declared in the shader.
struct ShaderPushConsts
    val::Float32
    n::UInt32
end

# Specialization constants are similar to push constants, but less dynamic: You
# can change them before compiling the shader for the pipeline, but not
# dynamically. This may have performance benefits for "very static" values,
# such as block sizes.
struct ShaderSpecConsts
    local_size_x::UInt32
end

# Let's now compile the shader to SPIR-V with `glslang`. We can use the artifact `glslang_jll` which provides the binary through the [Artifact system](https://pkgdocs.julialang.org/v1/artifacts/).
using glslang_jll: glslangValidator
glslang = glslangValidator(identity)
shader_bcode = mktempdir() do dir
    inpath = joinpath(dir, "shader.comp")
    outpath = joinpath(dir, "shader.spv")
    open(f -> write(f, shader_code), inpath, "w")
    status = run(`$glslang -V -S comp -o $outpath $inpath`)
    @assert status.exitcode == 0
    reinterpret(UInt32, read(outpath))
end

# We can now make a shader module with the compiled code:
shader = ShaderModule(device, sizeof(UInt32) * length(shader_bcode), shader_bcode)

# ## Assembling the pipeline
#
# Descriptor set layout describes how many resources of what kind are going to
# be used by the shader. In this case, we only use a single buffer:
dsl = DescriptorSetLayout(
    device,
    [
        DescriptorSetLayoutBinding(
            0,
            DESCRIPTOR_TYPE_STORAGE_BUFFER,
            SHADER_STAGE_COMPUTE_BIT;
            descriptor_count = 1,
        ),
    ],
)

# Pipeline layout describes the descriptor set together with the location of
# push constants:
pl = PipelineLayout(
    device,
    [dsl],
    [PushConstantRange(SHADER_STAGE_COMPUTE_BIT, 0, sizeof(ShaderPushConsts))],
)

# Shader compilation can use "specialization constants" that get propagated
# (and optimized) into the shader code. We use them to make the shader
# workgroup size "dynamic" in the sense that the size (32) is not hardcoded in
# GLSL, but instead taken from here.
const_local_size_x = 32
spec_consts = [ShaderSpecConsts(const_local_size_x)]

# Next, we create a pipeline that can run the shader code with the specified layout:
pipeline_info = ComputePipelineCreateInfo(
    PipelineShaderStageCreateInfo(
        SHADER_STAGE_COMPUTE_BIT,
        shader,
        "main", # this needs to match the function name in the shader
        specialization_info = SpecializationInfo(
            [SpecializationMapEntry(0, 0, 4)],
            UInt64(4),
            Ptr{Nothing}(pointer(spec_consts)),
        ),
    ),
    pl,
    -1,
)
ps, _ = unwrap(create_compute_pipelines(device, [pipeline_info]))
p = first(ps)

# Now make a descriptor pool to allocate the buffer descriptors from (not a big
# one, just 1 descriptor set with 1 descriptor in total), ...
dpool = DescriptorPool(device, 1, [DescriptorPoolSize(DESCRIPTOR_TYPE_STORAGE_BUFFER, 1)])

# ... allocate the descriptors for our layout, ...
dsets = unwrap(allocate_descriptor_sets(device, DescriptorSetAllocateInfo(dpool, [dsl])))
dset = first(dsets)

# ... and make the descriptors point to the right buffers.
update_descriptor_sets(
    device,
    [
        WriteDescriptorSet(
            dset,
            0,
            0,
            DESCRIPTOR_TYPE_STORAGE_BUFFER,
            [],
            [DescriptorBufferInfo(buffer, 0, WHOLE_SIZE)],
            [],
        ),
    ],
    [],
)

# ## Executing the shader
#
# Let's create a command pool in the right queue family, and take a command
# buffer out of that.
cmdpool = CommandPool(device, qfam_idx)
cbufs = unwrap(
    allocate_command_buffers(
        device,
        CommandBufferAllocateInfo(cmdpool, COMMAND_BUFFER_LEVEL_PRIMARY, 1),
    ),
)
cbuf = first(cbufs)

# Now that we have a command buffer, we can fill it with commands that cause
# the kernel to be run. Basically, we bind and fill everything, and then
# dispatch a sufficient amount of invocations of the shader to span over the
# array.
begin_command_buffer(
    cbuf,
    CommandBufferBeginInfo(flags = COMMAND_BUFFER_USAGE_ONE_TIME_SUBMIT_BIT),
)

cmd_bind_pipeline(cbuf, PIPELINE_BIND_POINT_COMPUTE, p)

const_buf = [ShaderPushConsts(1.234, data_items)]
cmd_push_constants(
    cbuf,
    pl,
    SHADER_STAGE_COMPUTE_BIT,
    0,
    sizeof(ShaderPushConsts),
    Ptr{Nothing}(pointer(const_buf)),
)

cmd_bind_descriptor_sets(cbuf, PIPELINE_BIND_POINT_COMPUTE, pl, 0, [dset], [])

cmd_dispatch(cbuf, div(data_items, const_local_size_x, RoundUp), 1, 1)

end_command_buffer(cbuf)

# Finally, find a handle to the compute queue and send the command to execute
# the shader!
compute_q = get_device_queue(device, qfam_idx, 0)
unwrap(queue_submit(compute_q, [SubmitInfo([], [], [cbuf], [])]))

# ## Getting the data
#
# After submitting the queue, the data is being crunched in the background. To
# get the resulting data, we need to wait for completion and invalidate the
# mapped memory (so that whatever data updates that happened on the GPU get
# transferred to the mapped range visible for the host).
#
# While [`queue_wait_idle`](@ref) will wait for computations to be carried out,
# we need to make sure that the required data is kept alive during queue
# operations. In non-global scopes, such as functions, the compiler may skip
# the allocation of unused variables or garbage-collect objects that the
# runtime thinks are no longer used. If garbage-collected, objects will call
# their finalizers which imply the destruction of the Vulkan objects
# (via `vkDestroy...`). In this particular case, the runtime is not aware
# that for example the pipeline and buffer objects are still used and that
# there's a dependency with these variables until the command returns, so we
# tell it manually.
GC.@preserve buff dsl pl p const_buf begin
    unwrap(queue_wait_idle(compute_q))
end

# Free the command buffers and the descriptor sets. These are the only handles that are not cleaned up automatically (see [Automatic finalization](@ref)).

unwrap(free_command_buffers(device, cmdpool, cbufs))
unwrap(free_descriptor_sets(device, dpool, dsets))

# Just as with flushing, the invalidation is only required for memory that is
# not host-coherent. You may skip this step if you check that the memory has
# the host-coherent property flag.
unwrap(invalidate_mapped_memory_ranges(device, [MappedMemoryRange(mem, 0, mem_size)]))

# Finally, let's have a look at the data created by your compute shader!
data # WHOA

max_stars_repo_name: JuliaGPU/VulkanAbstraction

max_issues_repo_name: JuliaGPU/VulkanAbstraction

max_forks_repo_name: JuliaGPU/VulkanAbstraction


text:
"""
    NeoHook(; E, ν)

Hyperelastic material
#Arguments
- `λ::Float64`: Lamé parameter
- `μ::Float64`: Lamé parameter (shear modulus)
"""

struct NeoHook <: AbstractMaterial
    λ::Float64
    μ::Float64
end

struct NeoHookState <: AbstractMaterialState
end

function initial_material_state(::NeoHook)
    return NeoHookState()
end

function NeoHook(; λ::T, μ::T) where T

    return NeoHook(λ, μ)
end

strainmeasure(::NeoHook) = RightCauchyGreen()

function ψ(mp::NeoHook, C::SymmetricTensor{2,3})
    J = sqrt(det(C))
    I = tr(C)
    return mp.μ/2 * (I-3) - mp.μ*log(J) + mp.λ/2 * log(J)^2
end

function material_response(mp::NeoHook, C::SymmetricTensor{2,3}, state::NeoHookState = NeoHookState(),
                           Δt=nothing; cache=nothing, options=nothing)
    invC = inv(C)
    J = sqrt(det(C))
    S = mp.μ*(one(SymmetricTensor{2,3}) - inv(C)) + mp.λ*log(J)*inv(C)
    ∂S∂C = mp.λ*(invC⊗invC) + 2*(mp.μ-mp.λ*log(J))*otimesu(invC,invC)
    ∂S∂C = symmetric(∂S∂C)
    return S, ∂S∂C, NeoHookState()
end


text:
"""
    init_input_layer(res_size::Int, in_size::Int, sigma::Float64)

Return a weighted input layer matrix, with random non-zero elements drawn from \$ [-\\text{sigma}, \\text{sigma}] \$, as described in [1].

[1] Lu, Zhixin, et al. "Reservoir observers: Model-free inference of unmeasured variables in chaotic systems." Chaos: An Interdisciplinary Journal of Nonlinear Science 27.4 (2017): 041102.
"""
function init_input_layer(res_size::Int,
        in_size::Int,
        sigma::Float64)

    W_in = zeros(Float64, res_size, in_size)
    q = floor(Int, res_size/in_size) #need to fix the reservoir input size. Check the constructor
    for i=1:in_size
        W_in[(i-1)*q+1 : (i)*q, i] = (2*sigma).*(rand(Float64, 1, q).-0.5)
    end
    return W_in

end

"""
    init_dense_input_layer(res_size::Int, in_size::Int, sigma::Float64)

Return a fully connected input layer matrix, with random non-zero elements drawn from \$ [-sigma, sigma] \$.
"""
function init_dense_input_layer(res_size::Int,
        in_size::Int,
        sigma::Float64)

    W_in = rand(Float64, res_size, in_size)
    W_in = 2.0 .*(W_in.-0.5)
    W_in = sigma .*W_in
    return W_in
end

"""
    init_sparse_input_layer(res_size::Int, in_size::Int, sigma::Float64, sparsity::Float64)

Return a sparsely connected input layer matrix, with random non-zero elements drawn from \$ [-sigma, sigma] \$ and given sparsity.
"""
function init_sparse_input_layer(res_size::Int,
        in_size::Int,
        sigma::Float64,
        sparsity::Float64)

    W_in = Matrix(sprand(Float64, res_size, in_size, sparsity))
    W_in = 2.0 .*(W_in.-0.5)
    replace!(W_in, -1.0=>0.0)
    W_in = sigma .*W_in
    return W_in
end

#from "minimum complexity echo state network" Rodan
"""
    min_complex_input(res_size::Int, in_size::Int, weight::Float64)

Return a fully connected input layer matrix with the same weights and sign drawn from a Bernoulli distribution, as described in [1].

[1] Rodan, Ali, and Peter Tino. "Minimum complexity echo state network." IEEE transactions on neural networks 22.1 (2010): 131-144.
"""
function min_complex_input(res_size::Int,
        in_size::Int,
        weight::Float64)

    W_in = Array{Float64}(undef, res_size, in_size)
    for i=1:res_size
        for j=1:in_size
            if rand(Bernoulli()) == true
                W_in[i, j] = weight
            else
                W_in[i, j] = -weight
            end
        end
    end
    return W_in
end

#from "minimum complexity echo state network" Rodan
#and "simple deterministically constructed cycle reservoirs with regular jumps" by Rodan and Tino

"""
    irrational_sign_input(res_size::Int, in_size::Int , weight::Float64 [, start::Int, irrational::Irrational])

Return a fully connected input layer matrix with the same weights and sign decided by the values of an irrational number, as described in [1] and [2].

[1] Rodan, Ali, and Peter Tino. "Minimum complexity echo state network." IEEE transactions on neural networks 22.1 (2010): 131-144.
[2] Rodan, Ali, and Peter Tiňo. "Simple deterministically constructed cycle reservoirs with regular jumps." Neural computation 24.7 (2012): 1822-1852.
"""
function irrational_sign_input(res_size::Int,
        in_size::Int,
        weight::Float64;
        start::Int = 1,
        irrational::Irrational = pi)

    setprecision(BigFloat, Int(ceil(log2(10)*(res_size*in_size+start+1))))
    ir_string = string(BigFloat(irrational)) |> collect
    deleteat!(ir_string, findall(x->x=='.', ir_string))
    ir_array = Array{Int}(undef, length(ir_string))
    W_in = Array{Float64}(undef, res_size, in_size)

    for i =1:length(ir_string)
        ir_array[i] = parse(Int, ir_string[i])
    end

    counter = start

    for i=1:res_size
        for j=1:in_size
            if ir_array[counter] < 5
                W_in[i, j] = -weight
            else
                W_in[i, j] = weight
            end
            counter += 1
        end
    end
    return W_in
end

"""
physics_informed_input(res_size::Int, in_size::Int, sigma::Float64, γ::Float64)

Return a weighted input layer matrix, with random non-zero elements drawn from \$ [-\\text{sigma}, \\text{sigma}] \$, where some γ
of reservoir nodes are connected exclusively to the raw inputs, and the rest to the outputs of the prior knowledge model , as described in [1].

[1] Jaideep Pathak et al. "Hybrid Forecasting of Chaotic Processes: Using Machine Learning in Conjunction with a Knowledge-Based Model" (2018)
"""
function physics_informed_input(res_size::Int,
        in_size::Int,
        sigma::Float64,
        γ::Float64,
        model_in_size::Int)

    state_size = in_size - model_in_size
    if state_size <= 0
        throw(DimensionMismatch("in_size must be greater than model_in_size"))
    end
    W_in = zeros(Float64, res_size, in_size)
    #Vector used to find res nodes not yet connected
    zero_connections = zeros(in_size)
    #Num of res nodes allotted for raw states
    num_for_state = floor(Int, res_size*γ)
    #Num of res nodes allotted for prior model input
    num_for_model = floor(Int, (res_size*(1-γ)))
    for i in 1:num_for_state
        #find res nodes with no connections
        idxs = findall(Bool[zero_connections == W_in[i,:] for i=1:size(W_in,1)])
        random_row_idx = idxs[rand(1:end)]
        random_clm_idx = range(1, state_size, step = 1)[rand(1:end)]
        W_in[random_row_idx,random_clm_idx] = rand(Uniform(-sigma, sigma))
    end

    for i in 1:num_for_model
        idxs = findall(Bool[zero_connections == W_in[i,:] for i=1:size(W_in,1)])
        random_row_idx = idxs[rand(1:end)]
        random_clm_idx = range(state_size+1, in_size, step = 1)[rand(1:end)]
        W_in[random_row_idx,random_clm_idx] = rand(Uniform(-sigma, sigma))
    end
    return W_in
end


text:
#=
day10:
- Julia version: 1.7.0
- Author: Paul.Mealor
- Date: 2021-12-10
=#

using AoC, Test

SCORES = Dict(')'=>3, ']'=>57, '}'=>1197, '>'=>25137)
CLOSERS = Dict('('=>')', '['=>']', '{'=>'}', '<'=>'>')
AC_SCORES = Dict(')'=>1, ']'=>2, '}'=>3, '>'=>4)

function syntaxCheck(line)
    expectedClosers = []
    for c in line
        closer = get(CLOSERS, c, nothing)
        if isnothing(closer)
            if isempty(expectedClosers)
                return (nothing, expectedClosers)
            end
            expectedC = popfirst!(expectedClosers)
            if expectedC != c
                return (c, [])
            end
        else
            pushfirst!(expectedClosers, CLOSERS[c])
        end
    end
    return (nothing, expectedClosers)
end

incorrectChar(line) = syntaxCheck(line)[1]
@test incorrectChar("{>") == '>'
@test incorrectChar("{()>") == '>'
@test incorrectChar("{<)") == ')'
@test incorrectChar("{()}") == nothing
@test incorrectChar("{()") == nothing

part1(lines) = incorrectChar.(lines) |> cc -> sum(isnothing(c) ? 0 : SCORES[c] for c in cc)

@test part1(exampleLines(10,1)) == 26397

autocompleteChars(line) = syntaxCheck(line)[2]

function acScore(chars)
    score = 0
    for c in chars
        score *= 5
        score += AC_SCORES[c]
    end
    return score
end
@test acScore(['}','}',']',']',')','}',')',']']) == 288957

middleItem(A) = A[(length(A)+1) ÷ 2]
@test middleItem(['a','b','c']) == 'b'

function part2(lines)
    lines = filter(l->isnothing(incorrectChar(l)), lines)
    allChars = autocompleteChars.(lines)
    scores = acScore.(allChars)
    sort!(scores)
    return middleItem(scores)
end

@test part2(exampleLines(10,1)) == 288957

lines(10) |> ll -> @time part1(ll) |> show
lines(10) |> ll -> @time part2(ll) |> show


text:
data = [parse(Int,x) for x in readline("data/Dec8_data.txt")]

# Question 1
r = reshape(data, (25,6,:))
l = findmin(vec(mapslices(x->count(==(0),x),r, dims = 1:2)))[2]
count(==(1), r[:,:,l]) * count(==(2), r[:,:,l])

# Question 2
merg(a1, a2) = a1 == 2 ? a2 : a1
res = foldl((x,y) -> merg.(x,y), eachslice(r, dims = 3))
using Plots
heatmap(res', aspect_ratio = 1, yflip = true)


text:
export softmax,∂softmax

function softmax(rᵢ) #for prediction
  Z=sum(exp(rᵢ),1)
  rₒ=exp(rᵢ)./Z
end
function ∂softmax(rᵢ,rₒ)
  ∂ij=Diagm(rₒ) - rₒ⊗rₒ  # ∂rₒᵢ/∂rᵢⱼ
end


text:
export AbstractAgent, ABM, AgentBasedModel

#######################################################################################
# %% Fundamental type definitions
#######################################################################################
"""
All agents must be a mutable subtype of `AbstractAgent`.
Your agent type **must have** the `id` field as first field.
Depending on the space structure there might be a `pos` field of appropriate type
and a `vel` field of appropriate type.

Your agent type may have other additional fields relevant to your system,
for example variable quantities like "status" or other "counters".

## Examples
Imagine agents who have extra properties `weight, happy`. For a [`GraphSpace`](@ref)
we would define them like
```julia
mutable struct ExampleAgent <: AbstractAgent
    id::Int
    pos::Int
    weight::Float64
    happy::Bool
end
```
while for e.g. a [`ContinuousSpace`](@ref) we would use
```julia
mutable struct ExampleAgent <: AbstractAgent
    id::Int
    pos::NTuple{2, Float64}
    vel::NTuple{2, Float64}
    weight::Float64
    happy::Bool
end
```
where `vel` is optional, useful if you want to use [`move_agent!`](@ref) in continuous
space.
"""
abstract type AbstractAgent end

abstract type AbstractSpace end
SpaceType=Union{Nothing, AbstractSpace}

abstract type DiscreteSpace <: AbstractSpace end

# This is a collection of valid position types, sometimes used for ambiguity resolution
ValidPos = Union{Int, NTuple{N, Int}, NTuple{M, <:AbstractFloat}} where {N, M}

struct AgentBasedModel{A<:AbstractAgent, S<:SpaceType, F, P}
    agents::Dict{Int,A}
    space::S
    scheduler::F
    properties::P
    maxid::Base.RefValue{Int64}
end

const ABM = AgentBasedModel

agenttype(::ABM{A}) where {A} = A
spacetype(::ABM{A, S}) where {A, S} = S

"""
    union_types(U)
Return a set of types within a `Union`. Preserves order.
"""
union_types(x::Union) = union_types(x.a, x.b)
union_types(a::Union, b::Type) = (union_types(a)..., b)
union_types(a::Type, b::Type) = (a, b)
union_types(x::Type) = (x,)
# For completness
union_types(a::Type, b::Union) = (a, union_types(b)...)

"""
    AgentBasedModel(AgentType [, space]; scheduler, properties) → model
Create an agent based model from the given agent type and `space`.
You can provide an agent _instance_ instead of type, and the type will be deduced.
 `ABM` is equivalent with `AgentBasedModel`.

The agents are stored in a dictionary that maps unique ids (integers)
to agents. Use `model[id]` to get the agent with the given `id`.

`space` is a subtype of `AbstractSpace`: [`GraphSpace`](@ref), [`GridSpace`](@ref) or
[`ContinuousSpace`](@ref).
If it is ommited then all agents are virtually in one position and have no spatial structure.

**Note:** Spaces are mutable objects and are not designed to be shared between models.
Create a fresh instance of a space with the same properties if you need to do this.

`properties = nothing` is additional model-level properties (typically a dictionary)
that can be accessed as `model.properties`. However, if `properties` is a dictionary with
key type `Symbol`, or of it is a struct, then the syntax
`model.name` is short hand for `model.properties[:name]` (or `model.properties.name`
for structs).
This syntax can't be used for `name` being `agents, space, scheduler, properties`,
which are the fields of `AgentBasedModel`.

`scheduler = fastest` decides the order with which agents are activated
(see e.g. [`by_id`](@ref) and the scheduler API).

Type tests for `AgentType` are done, and by default
warnings are thrown when appropriate. Use keyword `warn=false` to supress that.
"""
function AgentBasedModel(
        ::Type{A}, space::S = nothing;
        scheduler::F = fastest, properties::P = nothing, warn = true
        ) where {A<:AbstractAgent, S<:SpaceType, F, P}
    agent_validator(A, space, warn)

    agents = Dict{Int, A}()
    return ABM{A, S, F, P}(agents, space, scheduler, properties, Ref(0))
end

function AgentBasedModel(agent::AbstractAgent, args...; kwargs...)
    return ABM(typeof(agent), args...; kwargs...)
end

#######################################################################################
# %% Model accessing api
#######################################################################################
export random_agent, nagents, allagents, allids, nextid

"""
    model[id]
    getindex(model::ABM, id::Integer)

Return an agent given its ID.
"""
Base.getindex(m::ABM, id::Integer) = m.agents[id]

"""
    model[id] = agent
    setindex!(model::ABM, agent::AbstractAgent, id::Int)

Add an `agent` to the `model` at a given index: `id`.
Note this method will return an error if the `id` requested is not equal to `agent.id`.
**Internal method**, use [`add_agents!`](@ref) instead to actually add an agent.
"""
function Base.setindex!(m::ABM, a::AbstractAgent, id::Int)
    a.id ≠ id && throw(ArgumentError("You are adding an agent to an ID not equal with the agent's ID!"))
    m.agents[id] = a
    m.maxid[] < id && (m.maxid[] += 1)
    return a
end

"""
    nextid(model::ABM) → id
Return a valid `id` for creating a new agent with it.
"""
nextid(model::ABM) = model.maxid[] + 1

"""
    model.prop
    getproperty(model::ABM, prop::Symbol)

Return a property from the current `model`, assuming the model `properties` are either
a dictionary with key type `Symbol` or a Julia struct.
For example, if a model has the set of properties `Dict(:weight => 5, :current => false)`,
retrieving these values can be obtained via `model.weight`.

The property names `:agents, :space, :scheduler, :properties, :maxid` are internals
and **should not be accessed by the user**.
"""
function Base.getproperty(m::ABM{A, S, F, P}, s::Symbol) where {A, S, F, P}
    if s === :agents
        return getfield(m, :agents)
    elseif s === :space
        return getfield(m, :space)
    elseif s === :scheduler
        return getfield(m, :scheduler)
    elseif s === :properties
        return getfield(m, :properties)
    elseif s === :maxid
        return getfield(m, :maxid)
    elseif P <: Dict
        return getindex(getfield(m, :properties), s)
    else # properties is assumed to be a struct
        return getproperty(getfield(m, :properties), s)
    end
end

function Base.setproperty!(m::ABM{A, S, F, P}, s::Symbol, x) where {A, S, F, P}
    properties = getfield(m, :properties)
    if properties ≠ nothing && haskey(properties, s)
        properties[s] = x
    else
        throw(ErrorException("Cannot set $(s) in this manner. Please use the `AgentBasedModel` constructor."))
    end
end

"""
    random_agent(model)
Return a random agent from the model.
"""
random_agent(model) = model[rand(keys(model.agents))]

"""
    nagents(model::ABM)
Return the number of agents in the `model`.
"""
nagents(model::ABM) = length(model.agents)

"""
    allagents(model)
Return an iterator over all agents of the model.
"""
allagents(model) = values(model.agents)

"""
    allids(model)
Return an iterator over all agent IDs of the model.
"""
allids(model) = keys(model.agents)

#######################################################################################
# %% Model construction validation
#######################################################################################
"""
    agent_validator(agent, space)
Validate the user supplied agent (subtype of `AbstractAgent`).
Checks for mutability and existence and correct types for fields depending on `SpaceType`.
"""
function agent_validator(::Type{A}, space::S, warn::Bool) where {A<:AbstractAgent, S<:SpaceType}
    # Check A for required properties & fields
    if isconcretetype(A)
        do_checks(A, space, warn)
    else
        warn && @warn "AgentType is not concrete. If your agent is parametrically typed, you're probably seeing this warning because you gave `Agent` instead of `Agent{Float64}` (for example) to this function. You can also create an instance of your agent and pass it to this function. If you want to use `Union` types for mixed agent models, you can silence this warning."
        for type in union_types(A)
            do_checks(type, space, warn)
        end
    end
end

"""
    do_checks(agent, space)
Helper function for `agent_validator`.
"""
function do_checks(::Type{A}, space::S, warn::Bool) where {A<:AbstractAgent, S<:SpaceType}
    if warn
        isbitstype(A) && @warn "AgentType should be mutable. Try adding the `mutable` keyword infront of `struct` in your agent definition."
    end
    (any(isequal(:id), fieldnames(A)) && fieldnames(A)[1] == :id) || throw(ArgumentError("First field of Agent struct must be `id` (it should be of type `Int`)."))
    fieldtype(A, :id) <: Integer || throw(ArgumentError("`id` field in Agent struct must be of type `Int`."))
    if space != nothing
        (any(isequal(:pos), fieldnames(A)) && fieldnames(A)[2] == :pos) || throw(ArgumentError("Second field of Agent struct must be `pos` when using a space."))
        # Check `pos` field in A has the correct type
        pos_type = fieldtype(A, :pos)
        space_type = typeof(space)
        if space_type <: GraphSpace && !(pos_type <: Integer)
            throw(ArgumentError("`pos` field in Agent struct must be of type `Int` when using GraphSpace."))
        elseif space_type <: GridSpace && !(pos_type <: NTuple{D, Integer} where {D})
            throw(ArgumentError("`pos` field in Agent struct must be of type `NTuple{Int}` when using GridSpace."))
        elseif space_type <: ContinuousSpace
            if !(pos_type <: NTuple{D, <:AbstractFloat} where {D})
                throw(ArgumentError("`pos` field in Agent struct must be of type `NTuple{<:AbstractFloat}` when using ContinuousSpace."))
            end
            if warn && any(isequal(:vel), fieldnames(A)) && !(fieldtype(A, :vel) <: NTuple{D, <:AbstractFloat} where {D})
                @warn "`vel` field in Agent struct should be of type `NTuple{<:AbstractFloat}` when using ContinuousSpace."
            end
        end
    end
end

#######################################################################################
# %% Pretty printing
#######################################################################################
function Base.show(io::IO, abm::ABM{A}) where {A}
    n = isconcretetype(A) ? nameof(A) : string(A)
    s = "AgentBasedModel with $(nagents(abm)) agents of type $(n)"
    if abm.space == nothing
        s*= "\n no space"
    else
        s*= "\n space: $(sprint(show, abm.space))"
    end
    s*= "\n scheduler: $(schedulername(abm.scheduler))"
    print(io, s)
    if abm.properties ≠ nothing
        print(io, "\n properties: ", abm.properties)
    end
end
schedulername(x::Union{Function, DataType}) = nameof(x)
schedulername(x) = string(x)


text:
module JungleHelperSpiderBoss

using ..Ahorn, Maple

@mapdef Entity "JungleHelper/SpiderBoss" SpiderBoss(x::Integer, y::Integer, color::String="Blue", sprite::String="", webSprite::String="", flag::String="")

const bossColors = String["Blue", "Purple", "Red"]

const bossSprites = Dict{String, String}(
	"Blue" => "JungleHelper/SpiderBoss/spider_b_00",
	"Purple" => "JungleHelper/SpiderBoss/spider_p_00",
	"Red" => "JungleHelper/SpiderBoss/spider_r_00"
)

const placements = Ahorn.PlacementDict(
	"Spider Boss ($(color)) (Jungle Helper)" => Ahorn.EntityPlacement(
		SpiderBoss,
		"point",
		Dict{String, Any}(
			"color" => color
		)
	) for color in bossColors
)

Ahorn.editingOptions(entity::SpiderBoss) = Dict{String, Any}(
	"color" => bossColors
)

function Ahorn.selection(entity::SpiderBoss)
    x, y = Ahorn.position(entity)
    return Ahorn.Rectangle(x - 9, y - 9, 17, 17)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::SpiderBoss, room::Maple.Room)
	color = get(entity.data, "color", "Blue")
	Ahorn.drawSprite(ctx, bossSprites[color], 0, 0)
end

end

max_stars_repo_name: Jackal-Celeste/JungleHelper

max_issues_repo_name: Jackal-Celeste/JungleHelper

max_forks_repo_name: Jackal-Celeste/JungleHelper


text:
module NN

include("backprop.jl")
include("errors.jl")
include("activation.jl")

using LinearAlgebra

struct perceptron
  activation::Function
end

struct layer
  layer::Array{perceptron,1}
end

struct NN
  NN::Array{layer,1}
end

perceptron_compute(x::perceptron,W::Array{Float64,1},I::Array{Float64,2}) = x.activation.(I * W)
layer_compute(l::layer,W::Array{Float64,2},I::Array{Float64,2}) = [pcom(i,W[:,j],I) for (i,j) in zip(l.layer,1:size(W,2)]

NeuralNetwork(x::Array{Tuple{Function,Int64},1}) = NN([layer([perceptron(x[i][1]) for k in 1:(x[i][2]+1)]) for i in 1:length(x)])

end


text:
using JavaCall
using JProxies
using Test

#const pxyptr = JavaCall.pxyptr
@testset "proxy_initialization" begin
    @test JProxies.init(["-Djava.class.path=$(@__DIR__)"])==nothing
end

@testset "proxy_array_list" begin
    JAL = @jimport java.util.ArrayList
    @test JProxy(JavaCall.jnew(Symbol("java.lang.Integer"), (jint,), 3)).toString() == "3"
    @test JProxy(@jimport(java.lang.Integer)).valueOf(3) == 3
    a = JProxy(JAL())
    @test a.size() == 0
    a.add("one")
    @test a.size() == 1
    @test a.toString() == "[one]"
    removed = a.remove(0)
    @test typeof(removed) == String
    @test removed == "one"
    a.add(1)
    @test a.get(0) == 1
    @test a.toString() == "[1]"
    b = JProxy(JAL())
    b.addAll(a)
    @test a.toString() == b.toString()
    a.add("two")
    @test collect(a) == [1, "two"]
end

@testset "proxy_test_class" begin
    T = @jimport(Test)
    t = JProxy(T(()))
    t.integerField = 3
    @test(t.integerField == 3)
    t.stringField = "hello"
    @test(t.stringField == "hello")
    @test(t.toString() == "Test(3, hello)")
    t.objectField = t
    @test(t.objectField == t)
    @test(t.objectField.stringField == "hello")
    @test(t.objectField.getInt() == 3)
    @test(t.objectField.getString() == "hello")
end

@testset "proxy_meta" begin
    @test(JProxy(@jimport(java.lang.Integer)).MAX_VALUE == 2147483647)
    @test(JProxy(@jimport(java.lang.Long)).MAX_VALUE == 9223372036854775807)
    @test(JProxy(@jimport(java.lang.Double)).MAX_VALUE == 1.7976931348623157e308)
    @test("class java.lang.Object" == JProxy(JavaCall.metaclass("java.lang.Class")).forName("java.lang.Object").toString())
    @test("class java.io.PrintStream" == JProxy(JavaObject{Symbol("java.lang.System")}).out.getClass().toString())
    @test("static class java.util.Arrays" == string(JProxy(@jimport(java.util.Arrays))))
end

@testset "proxy_array" begin
    s,ptr=JavaCall.convert_arg(JProxy{Array{Int,1}, false}, [1,2]) # convert Julia array to an unwrapped java array
    p=JProxy(ptr) # wrap it
    @test(length(p) == 2)
    @test(p[1] == 1)
    @test(p[2] == 2)
    @test(collect(p) == [1, 2])
    p[1] = 3
    @test(p[1] == 3)
end


text:
#!/usr/bin/env julia
using Pkg;

# Add support for various MIMI models
Pkg.Registry.add(RegistrySpec(url = "https://github.com/mimiframework/MimiRegistry.git"));

Pkg.add.([
    Pkg.PackageSpec(;name="MimiFUND", version="3.8.6"),
    Pkg.PackageSpec(;name="MimiIWG", version="1.1.0"),
    Pkg.PackageSpec(;name="MimiDICE2010", version="1.0.1"),
    Pkg.PackageSpec(;name="MimiDICE2013", version="1.0.2"),
    Pkg.PackageSpec(;name="MimiPAGE2009", version="3.1.0"),
    Pkg.PackageSpec(;name="MimiRICE2010", version="3.0.3"),
    Pkg.PackageSpec(;name="Query", version="1.0.0"),
    Pkg.PackageSpec(;name="VegaLite", version="2.6.0"),
    Pkg.PackageSpec(;name="CSVFiles", version="1.0.1"),
    Pkg.PackageSpec(;name="Distributions", version="0.23.11"),
    Pkg.PackageSpec(;name="DataFrames", version="0.21.8"),
    Pkg.PackageSpec(;name="Plots", version="1.24.3"),
    Pkg.PackageSpec(;name="Images", version="0.24.1"),
    Pkg.PackageSpec(;name="PyPlot", version="2.10.0"),
    Pkg.PackageSpec(;name="PyCall", version="1.92.5"),
    Pkg.PackageSpec(;name="Roots", version="1.3.11"),
    Pkg.PackageSpec(;name="Polynomials", version="2.0.18"),
    Pkg.PackageSpec(;name="DifferentialEquations", version="6.18.0"),
    Pkg.PackageSpec(;name="Mimi", version="1.4.0"),
    Pkg.PackageSpec(;name="Optim", version="1.5.0")
]);

# Precompiling installed packages
Pkg.precompile();

max_stars_repo_name: petersudmant/datahub

max_issues_repo_name: petersudmant/datahub

max_forks_repo_name: petersudmant/datahub


text:
"""
    sumr(;binfile::AbstractString="", outdir::AbstractString=".", pseudocount::Number=1.0)
Extract the summary information of data matrix.
`binfile` is a Julia Binary file generated by `csv2bin` function.
`outdir` is specified the directory you want to save the result.
`pseudocount` is specified to avoid NaN by log10(0) and used when `Row_LogMeans.csv` <log10(mean+pseudocount) value of each Row> is generated.
Output Files
---------
- `Col_Means.csv` : Mean in each column.
- `Col_LogMeans.csv` : Mean of Log10(count+pseudocount) in each column.
- `Col_FTTMeans.csv` : Mean of sqrt(count) + sqrt(count+1) in ech column.
- `Col_Vars.csv` : Variance in each column.
- `Col_LogVars.csv` : Variance of Log10(count+pseudocount) in each column.
- `Col_FTTVars.csv` : Variance of sqrt(count) + sqrt(count+1) in ech column.
- `Col_SDs.csv` : Standard deviation in each column.
- `Col_LogSDs.csv` : Standard deviation of Log10(count+pseudocount) in each column.
- `Col_FTTSDs.csv` : Standard deviation of sqrt(count) + sqrt(count+1) in ech column.
- `Col_CV2s.csv` : Coefficient of Variation in each column.
- `Col_NoZeros.csv` : Number of zero-elements in each column.
- `Col_NoCounts.csv` : Sum of counts in each column.
- `Row_Means.csv` : Mean in each row.
- `Row_LogMeans.csv` : Mean of Log10(count+pseudocount) in each row.
- `Row_FTTMeans.csv` : Mean of sqrt(count) + sqrt(count+1) in ech row.
- `Row_Vars.csv` : Variance in each row.
- `Row_LogVars.csv` : Variance of Log10(count+pseudocount) in each row.
- `Row_FTTVars.csv` :  Variance of sqrt(count) + sqrt(count+1) in ech row.
- `Row_SDs.csv` : Standard deviation in each row.
- `Row_LogSDs.csv` : Standard deviation of Log10(count+pseudocount) in each row.
- `Row_FTTSDs.csv` : Standard deviation of sqrt(count) + sqrt(count+1) in ech row.
- `Row_CV2s.csv` : Coefficient of Variation in each row.
- `Row_NoZeros.csv` : Number of zero-elements in each row.
- `Row_NoCounts.csv` : Sum of counts in each row.
"""
function sumr(;binfile::AbstractString="", outdir::AbstractString=".", pseudocount::Number=1.0)
    # 1 / 2 : Column-wise statistics
    println("1 / 2 : Column-wise statistics are calculated...")
    Col_Means, Col_LogMeans, Col_FTTMeans, Col_Vars, Col_LogVars, Col_FTTVars, Col_SDs, Col_LogSDs, Col_FTTSDs, Col_CV2s, Col_NoZeros, Col_NoCounts = colstats(binfile, pseudocount)

    # 2 / 2 : Row-wise statistics
    println("2 / 2 : Row-wise statistics are calculated...")
    Row_Means, Row_LogMeans, Row_FTTMeans, Row_Vars, Row_LogVars, Row_FTTVars, Row_SDs, Row_LogSDs, Row_FTTSDs,  Row_CV2s, Row_NoZeros, Row_NoCounts = rowstats(binfile, pseudocount)

    # Save
    # Col-wise
    writecsv(joinpath(outdir, "Col_Means.csv"), Col_Means)
    writecsv(joinpath(outdir, "Col_LogMeans.csv"), Col_LogMeans)
    writecsv(joinpath(outdir, "Col_FTTMeans.csv"), Col_FTTMeans)
    writecsv(joinpath(outdir, "Col_Vars.csv"), Col_Vars)
    writecsv(joinpath(outdir, "Col_LogVars.csv"), Col_LogVars)
    writecsv(joinpath(outdir, "Col_FTTVars.csv"), Col_FTTVars)
    writecsv(joinpath(outdir, "Col_SDs.csv"), Col_SDs)
    writecsv(joinpath(outdir, "Col_LogSDs.csv"), Col_LogSDs)
    writecsv(joinpath(outdir, "Col_FTTSDs.csv"), Col_FTTSDs)
    writecsv(joinpath(outdir, "Col_CV2s.csv"), Col_CV2s)
    writecsv(joinpath(outdir, "Col_NoZeros.csv"), Col_NoZeros)
    writecsv(joinpath(outdir, "Col_NoCounts.csv"), Col_NoCounts)

    # Row-wise
    writecsv(joinpath(outdir, "Row_Means.csv"), Row_Means)
    writecsv(joinpath(outdir, "Row_LogMeans.csv"), Row_LogMeans)
    writecsv(joinpath(outdir, "Row_FTTMeans.csv"), Row_FTTMeans)
    writecsv(joinpath(outdir, "Row_Vars.csv"), Row_Vars)
    writecsv(joinpath(outdir, "Row_LogVars.csv"), Row_LogVars)
    writecsv(joinpath(outdir, "Row_FTTVars.csv"), Row_FTTVars)
    writecsv(joinpath(outdir, "Row_SDs.csv"), Row_SDs)
    writecsv(joinpath(outdir, "Row_LogSDs.csv"), Row_LogSDs)
    writecsv(joinpath(outdir, "Row_FTTSDs.csv"), Row_FTTSDs)
    writecsv(joinpath(outdir, "Row_CV2s.csv"), Row_CV2s)
    writecsv(joinpath(outdir, "Row_NoZeros.csv"), Row_NoZeros)
    writecsv(joinpath(outdir, "Row_NoCounts.csv"), Row_NoCounts)
end

# Column-wise statistics
function colstats(binfile::AbstractString, pseudocount::Number)
    N, M = nm([binfile])
    N = N[]
    M = M[]
    tmpN = zeros(UInt32, 1)
    tmpM = zeros(UInt32, 1)
    x = zeros(UInt32, M)
    m = zeros(M)
    lm = zeros(M)
    fttm = zeros(M)
    v = zeros(M)
    lv = zeros(M)
    fttv = zeros(M)
    sd = zeros(M)
    lsd = zeros(M)
    fttsd = zeros(M)
    c = zeros(M)
    nz = zeros(M)
    nc = zeros(M)
    # Mean Calculator
    progress = Progress(N)
    open(binfile) do file
        stream = ZstdDecompressorStream(file)
        read!(stream, tmpN)
        read!(stream, tmpM)
        for n = 1:N
            # Data Import
            read!(stream, x)
            # Update
            m = m .+ (x .- m) / n
            lm = lm .+ (log10.(x .+ pseudocount) .- lm) / n
            fttm = fttm .+ (sqrt.(x) .+ sqrt.(x .+ 1) .- lm) / n
            # Non-zero counts
            nc = nc .+ x
            tmpx = x
            for mm = 1:M
                if tmpx[mm] != 0
                    tmpx[mm] = 1
                end
            end
            nz = nz .+ tmpx
            # Progress Bar
            next!(progress)
        end
        close(stream)
    end
    # Sum of squared deviation Calculator
    progress = Progress(N)
    open(binfile) do file
        stream = ZstdDecompressorStream(file)
        read!(stream, tmpN)
        read!(stream, tmpM)
        for n = 1:N
            # Data Import
            read!(stream, x)
            # Update
            tmpv = x .- m
            tmplv = log10.(x .+ pseudocount) .- lm
            tmpfttv = sqrt.(x) .+ sqrt.(x .+ 1) .- fttm
            v = v .+ tmpv .* tmpv
            lv = lv .+ tmplv .* tmplv
            fttv = fttv .+ tmpfttv .* tmpfttv
            # Progress Bar
            next!(progress)
        end
        close(stream)
    end
    # Variance
    v = v ./ N
    lv = lv ./ N
    fttv = fttv ./ N
    # Var → SD
    sd = sqrt.(v)
    lsd = sqrt.(lv)
    fttsd = sqrt.(fttv)
    # CV2
    c = v ./ m .* m
    return m, lm, fttm, v, lv, fttv, sd, lsd, fttsd, c, nz, nc
end

# Row-wise statistics
function rowstats(binfile::AbstractString, pseudocount::Number)
    N, M = nm([binfile])
    N = N[]
    M = M[]
    tmpN = zeros(UInt32, 1)
    tmpM = zeros(UInt32, 1)
    x = zeros(UInt32, M)
    m = zeros(N)
    lm = zeros(N)
    fttm = zeros(N)
    v = zeros(N)
    lv = zeros(N)
    fttv = zeros(N)
    sd = zeros(N)
    lsd = zeros(N)
    fttsd = zeros(N)
    c = zeros(N)
    nz = zeros(UInt32, N)
    nc = zeros(UInt32, N)
    progress = Progress(N)
    open(binfile) do file
        stream = ZstdDecompressorStream(file)
        read!(stream, tmpN)
        read!(stream, tmpM)
        for n = 1:N
            # Data Import
            read!(stream, x)
            # Update
            m[n] = mean(x)
            lm[n] = mean(log10.(x .+ pseudocount))
            fttm[n] = mean(sqrt.(x) .+ sqrt.(x .+ 1))
            v[n] = var(x)
            lv[n] = var(log10.(x .+ pseudocount))
            fttv[n] = var(sqrt.(x) .+ sqrt.(x .+ 1))
            sd[n] = var(x)^(1/2)
            lsd[n] = var(log10.(x .+ pseudocount))^(1/2)
            fttsd[n] = var(sqrt.(x) .+ sqrt.(x .+ 1))^(1/2)
            c[n] = v[n] ./ (m[n] .* m[n])
            for mm = 1:M
                if x[mm] != 0
                    nz[n] = nz[n] + 1
                end
            end
            nc[n] = sum.(x)[1]
            # Progress Bar
            next!(progress)
        end
        close(stream)
    end
    return m, lm, fttm, v, lv, fttv, sd, lsd, fttsd, c, nz, nc
end


text:
using StochasticDiffEq, Test, Random, DiffEqNoiseProcess
Random.seed!(100)

struct NoIndexArray{T, N} <: AbstractArray{T, N}
    x::Array{T, N}
end
Base.size(x::NoIndexArray) = size(x.x)
Base.axes(x::NoIndexArray) = axes(x.x)
Base.similar(x::NoIndexArray, dims::Union{Integer, AbstractUnitRange}...) = NoIndexArray(similar(x.x, dims...))
Base.copyto!(x::NoIndexArray, y::NoIndexArray) = NoIndexArray(copyto!(x.x, y.x))
Base.copy(x::NoIndexArray) = NoIndexArray(copy(x.x))
Base.fill!(x::NoIndexArray, y) = NoIndexArray(fill!(x.x, y))
Base.mapreduce(f, op, x::NoIndexArray; kwargs...) = mapreduce(f, op, x.x; kwargs...)
Base.any(f::Function, x::NoIndexArray; kwargs...) = any(f, x.x; kwargs...)
Base.all(f::Function, x::NoIndexArray; kwargs...) = all(f, x.x; kwargs...)
Base.:(==)(x::NoIndexArray, y::NoIndexArray) = x.x == y.x

struct NoIndexStyle{N} <: Broadcast.AbstractArrayStyle{N} end
NoIndexStyle(::Val{N}) where N = NoIndexStyle{N}()
NoIndexStyle{M}(::Val{N}) where {N,M} = NoIndexStyle{N}()
Base.BroadcastStyle(::Type{<:NoIndexArray{T,N}}) where {T,N} = NoIndexStyle{N}()
Base.similar(bc::Base.Broadcast.Broadcasted{NoIndexStyle{N}}, ::Type{ElType}) where {N, ElType} = NoIndexArray(similar(Array{ElType, N}, axes(bc)))
Base.Broadcast._broadcast_getindex(x::NoIndexArray, i) = x.x[i]
Base.Broadcast.extrude(x::NoIndexArray) = x

@inline function Base.copyto!(dest::NoIndexArray, bc::Base.Broadcast.Broadcasted{<:NoIndexStyle})
    axes(dest) == axes(bc) || throwdm(axes(dest), axes(bc))
    bc′ = Base.Broadcast.preprocess(dest, bc)
    dest′ = dest.x
    @simd for I in eachindex(bc′)
        @inbounds dest′[I] = bc′[I]
    end
    return dest
end
Base.show_vector(io::IO, x::NoIndexArray) = Base.show_vector(io, x.x)

Base.show(io::IO, x::NoIndexArray) = (print(io, "NoIndexArray");show(io, x.x))
function Base.show(io::IO, ::MIME"text/plain", x::NoIndexArray)
    println(io, Base.summary(x), ":")
    Base.print_array(io, x.x)
end

prob = SDEProblem((du, u, p, t)->copyto!(du, u),(du, u, p, t)->copyto!(du, u),NoIndexArray(ones(10, 10)),(0.0,10.0))
algs = [SOSRI(), SOSRA()]
DiffEqNoiseProcess.wiener_randn!(rng::AbstractRNG,rand_vec::NoIndexArray) = randn!(rng,rand_vec.x)

for alg in algs
  sol = @test_nowarn solve(prob, alg)
  @test_nowarn sol(0.1)
  @test_nowarn sol(similar(prob.u0), 0.1)
end


text:
struct DefaultInit <: DiffEqBase.DAEInitializationAlgorithm end
struct NoInit <: DiffEqBase.DAEInitializationAlgorithm end

struct ShampineCollocationInit{T} <: DiffEqBase.DAEInitializationAlgorithm
  initdt::T
end
ShampineCollocationInit() = ShampineCollocationInit(nothing)

struct BrownFullBasicInit{T} <: DiffEqBase.DAEInitializationAlgorithm
  abstol::T
end
BrownFullBasicInit() = BrownFullBasicInit(1e-10)

## Notes

#=
differential_vars = [any(!iszero,x) for x in eachcol(M)]

A column should be zero for an algebraic variable, since that means that the
derivative term doesn't show up in any equations (i.e. is an algebraic variable).
The rows are not necessarily non-zero, for example a flux condition between two
differential variables. But if it's a condition that doesn't involve the algebraic
variable, then the system is not Index 1!

=#

## Expansion

function DiffEqBase.initialize_dae!(integrator::ODEIntegrator,
                  initializealg = integrator.initializealg)
  _initialize_dae!(integrator, integrator.sol.prob,
           initializealg,
           Val(DiffEqBase.isinplace(integrator.sol.prob)))
end

## Default algorithms

function _initialize_dae!(integrator, prob::ODEProblem,
             alg::DefaultInit, x::Val{true})
  _initialize_dae!(integrator, prob,
          BrownFullBasicInit(), x)
end

function _initialize_dae!(integrator, prob::ODEProblem,
             alg::DefaultInit, x::Val{false})
  _initialize_dae!(integrator, prob,
          BrownFullBasicInit(), x)
end

function _initialize_dae!(integrator, prob::DAEProblem,
             alg::DefaultInit, x::Val{false})
  if prob.differential_vars === nothing
    _initialize_dae!(integrator, prob,
            ShampineCollocationInit(), x)
  else
    _initialize_dae!(integrator, prob,
            BrownFullBasicInit(), x)
  end
end

function _initialize_dae!(integrator, prob::DAEProblem,
             alg::DefaultInit, x::Val{true})
  if prob.differential_vars === nothing
    _initialize_dae!(integrator, prob,
            ShampineCollocationInit(), x)
  else
    _initialize_dae!(integrator, prob,
            BrownFullBasicInit(), x)
  end
end

## NoInit

function _initialize_dae!(integrator, prob::ODEProblem,
             alg::NoInit, x::Val{true})
end

function _initialize_dae!(integrator, prob::ODEProblem,
             alg::NoInit, x::Val{false})
end

function _initialize_dae!(integrator, prob::DAEProblem,
             alg::NoInit, x::Val{false})
end

function _initialize_dae!(integrator, prob::DAEProblem,
             alg::NoInit, x::Val{true})
end

## ShampineCollocationInit

#=
The method:

du = (u-u0)/h
Solve for `u`

=#

function _initialize_dae!(integrator, prob::ODEProblem, alg::ShampineCollocationInit, ::Val{true})
  @unpack p, t, f = integrator
  M = integrator.f.mass_matrix
  dtmax = integrator.opts.dtmax
  tmp = first(get_tmp_cache(integrator))
  u0 = integrator.u

  dt = if alg.initdt === nothing
    integrator.dt != 0 ? min(integrator.dt/5, dtmax) : 1//1000 # Haven't implemented norm reduction
  else
    alg.initdt
  end

  algebraic_vars = [all(iszero,x) for x in eachcol(M)]
  algebraic_eqs  = [all(iszero,x) for x in eachrow(M)]
  (iszero(algebraic_vars) || iszero(algebraic_eqs)) && return
  update_coefficients!(M, u0, p, t)
  f(tmp, u0, p, t)
  tmp .= ArrayInterface.restructure(tmp,algebraic_eqs .* vec(tmp))

  integrator.opts.internalnorm(tmp, t) <= integrator.opts.abstol && return

  if isdefined(integrator.cache, :nlsolver)
    # backward Euler
    nlsolver = integrator.cache.nlsolver
    oldγ, oldc, oldmethod, olddt = nlsolver.γ, nlsolver.c, nlsolver.method, integrator.dt
    nlsolver.tmp .= integrator.uprev
    nlsolver.γ, nlsolver.c = 1, 1
    nlsolver.method = DIRK
    integrator.dt = dt
    z = nlsolve!(nlsolver, integrator, integrator.cache)
    nlsolver.γ, nlsolver.c, nlsolver.method, integrator.dt = oldγ, oldc, oldmethod, olddt
    # TODO: failure handling
    nlsolvefail(nlsolver) && @warn "ShampineCollocationInit DAE initialization algorithm failed with dt=$dt. Try to adjust initdt like `ShampineCollocationInit(initdt)`."
    @.. integrator.u = integrator.uprev + z
  else
    isad = alg_autodiff(integrator.alg)
    _tmp = isad ? PreallocationTools.dualcache(tmp, ForwardDiff.pickchunksize(length(tmp))) : tmp
    nlequation! = @closure (out,u) -> begin
      update_coefficients!(M,u,p,t)
      #M * (u-u0)/dt - f(u,p,t)
      tmp = isad ? PreallocationTools.get_tmp(_tmp, u) : _tmp
      @. tmp = (u - u0)/dt
      mul!(vec(out),M,vec(tmp))
      f(tmp,u,p,t)
      out .-= tmp
      nothing
    end
    r = nlsolve(nlequation!, u0, autodiff=isad ? :forward : :central, method = :newton)
    integrator.u .= r.zero
  end
  recursivecopy!(integrator.uprev,integrator.u)
  if alg_extrapolates(integrator.alg)
    recursivecopy!(integrator.uprev2,integrator.uprev)
  end
  return nothing
end

function _initialize_dae!(integrator, prob::ODEProblem, alg::ShampineCollocationInit, ::Val{false})
  @unpack p, t, f = integrator
  u0 = integrator.u
  M = integrator.f.mass_matrix
  dtmax = integrator.opts.dtmax

  dt = if alg.initdt === nothing
    integrator.dt != 0 ? min(integrator.dt/5, dtmax) : 1//1000 # Haven't implemented norm reduction
  else
    alg.initdt
  end

  algebraic_vars = [all(iszero,x) for x in eachcol(M)]
  algebraic_eqs  = [all(iszero,x) for x in eachrow(M)]
  (iszero(algebraic_vars) || iszero(algebraic_eqs)) && return
  update_coefficients!(M,u0,p,t)
  du = f(u0,p,t)
  resid = @view vec(du)[algebraic_eqs]

  integrator.opts.internalnorm(resid,t) <= integrator.opts.abstol && return

  if isdefined(integrator.cache, :nlsolver)
    # backward Euler
    nlsolver = integrator.cache.nlsolver
    oldγ, oldc, oldmethod, olddt = nlsolver.γ, nlsolver.c, nlsolver.method, integrator.dt
    nlsolver.tmp .= integrator.uprev
    nlsolver.γ, nlsolver.c = 1, 1
    nlsolver.method = DIRK
    integrator.dt = dt
    z = nlsolve!(nlsolver, integrator, integrator.cache)
    nlsolver.γ, nlsolver.c, nlsolver.method, integrator.dt = oldγ, oldc, oldmethod, olddt
    # TODO: failure handling
    nlsolvefail(nlsolver) && @warn "ShampineCollocationInit DAE initialization algorithm failed with dt=$dt. Try to adjust initdt like `ShampineCollocationInit(initdt)`."
    @.. integrator.u = integrator.uprev + z
  else
    nlequation_oop = @closure u -> begin
      update_coefficients!(M,u,p,t)
      M * (u-u0)/dt - f(u,p,t)
    end

    nlequation! = @closure (out,u) -> out .= nlequation_oop(u)

    integrator.u = nlsolve(nlequation!, u0).zero
  end
  integrator.uprev = integrator.u
  if alg_extrapolates(integrator.alg)
    integrator.uprev2 = integrator.uprev
  end
  return
end

function _initialize_dae!(integrator, prob::DAEProblem,
             alg::ShampineCollocationInit, ::Val{true})

  @unpack p, t, f = integrator
  u0 = integrator.u

  dtmax = integrator.opts.dtmax
  tmp = get_tmp_cache(integrator)[1]
  resid = get_tmp_cache(integrator)[2]

  dt = t != 0 ? min(t/1000,dtmax) : dtmax # Haven't implemented norm reduction

   nlequation! = @closure (out,u) -> begin
     #M * (u-u0)/dt - f(u,p,t)
    @. tmp = (u - u0)/dt
    f(out,tmp,u,p,t)
    nothing
   end

  nlequation!(tmp,u0)
  f(resid,integrator.du,u0,p,t)
   integrator.opts.internalnorm(resid,t) <= integrator.opts.abstol && return

  integrator.u .= nlsolve(nlequation!, u0).zero
  recursivecopy!(integrator.uprev,integrator.u)
  if alg_extrapolates(integrator.alg)
    recursivecopy!(integrator.uprev2,integrator.uprev)
  end
  return
end

function _initialize_dae!(integrator, prob::DAEProblem,
             alg::ShampineCollocationInit, ::Val{false})

  @unpack p, t, f = integrator
  u0 = integrator.u
  dtmax = integrator.opts.dtmax

  dt = t != 0 ? min(t/1000,dtmax/10) : dtmax # Haven't implemented norm reduction

  nlequation_oop = u -> begin
    f((u-u0)/dt,u,p,t)
  end

  nlequation! = @closure (out,u) -> out .= nlequation_oop(u)

  resid = f(integrator.du,u0,p,t)
  integrator.opts.internalnorm(resid,t) <= integrator.opts.abstol && return

  integrator.u = nlsolve(nlequation!, u0).zero
  integrator.uprev = integrator.u
  if alg_extrapolates(integrator.alg)
    integrator.uprev2 = integrator.uprev
  end
end

## BrownFullBasic

#=
The method:

Keep differential variables constant
Solve for the algebraic variables

=#

function _initialize_dae!(integrator, prob::ODEProblem, alg::BrownFullBasicInit, ::Val{true})
  @unpack p, t, f = integrator
  u = integrator.u
  M = integrator.f.mass_matrix
  update_coefficients!(M,u,p,t)
  algebraic_vars = [all(iszero,x) for x in eachcol(M)]
  algebraic_eqs  = [all(iszero,x) for x in eachrow(M)]
  (iszero(algebraic_vars) || iszero(algebraic_eqs)) && return
  tmp = get_tmp_cache(integrator)[1]

  f(tmp,u,p,t)

  tmp .= ArrayInterface.restructure(tmp,algebraic_eqs .* vec(tmp))

  integrator.opts.internalnorm(tmp,t) <= alg.abstol && return
  alg_u = @view u[algebraic_vars]

  isad = alg_autodiff(integrator.alg)
  if isad
    chunk = ForwardDiff.pickchunksize(count(algebraic_vars))
    _tmp = PreallocationTools.dualcache(tmp, chunk)
    _du_tmp = PreallocationTools.dualcache(tmp, chunk)
  else
    _tmp, _du_tmp = tmp, similar(tmp)
  end
  nlequation = @closure (out, x) -> begin
    uu = isad ? DiffEqBase.get_tmp(_tmp, x) : _tmp
    du_tmp = isad ? DiffEqBase.get_tmp(_du_tmp, x) : _du_tmp
    copyto!(uu, integrator.u)
    alg_uu = @view uu[algebraic_vars]
    alg_uu .= x
    f(du_tmp, uu, p, t)
    out .= @view du_tmp[algebraic_eqs]
    return nothing
  end

  r = nlsolve(nlequation, u[algebraic_vars], autodiff=isad ? :forward : :central, method=:newton)
  alg_u .= r.zero

  recursivecopy!(integrator.uprev,integrator.u)
  if alg_extrapolates(integrator.alg)
    recursivecopy!(integrator.uprev2,integrator.uprev)
  end

  return
end

function _initialize_dae!(integrator, prob::ODEProblem,
             alg::BrownFullBasicInit, ::Val{false})
  @unpack p, t, f = integrator

  u0 = integrator.u
  M = integrator.f.mass_matrix
  update_coefficients!(M,u0,p,t)
  algebraic_vars = [all(iszero,x) for x in eachcol(M)]
  algebraic_eqs  = [all(iszero,x) for x in eachrow(M)]
  (iszero(algebraic_vars) || iszero(algebraic_eqs)) && return

  du = f(u0,p,t)
  resid = @view vec(du)[algebraic_eqs]

  integrator.opts.internalnorm(resid,t) <= alg.abstol && return

  if u0 isa Number
    # This doesn't fix static arrays!
    u = [u0]
  else
    u = u0
  end

  alg_u = @view u[algebraic_vars]

  nlequation = @closure (out,x) -> begin
    alg_u .= x
    du = f(u,p,t)
    out .= @view du[algebraic_eqs]
  end
  r = nlsolve(nlequation, u0[algebraic_vars])
  alg_u .= r.zero

  if u0 isa Number
    # This doesn't fix static arrays!
    integrator.u = first(u)
  else
    integrator.u = u
  end

  integrator.uprev = integrator.u
  if alg_extrapolates(integrator.alg)
    integrator.uprev2 = integrator.uprev
  end

  return
end

function _initialize_dae!(integrator, prob::DAEProblem,
             alg::BrownFullBasicInit, ::Val{true})
  @unpack p, t, f = integrator
  differential_vars = prob.differential_vars
  u = integrator.u
  du = integrator.du

  tmp = get_tmp_cache(integrator)[1]
  f(tmp, du, u, p, t)

  if integrator.opts.internalnorm(tmp,t) <= alg.abstol
    return
  elseif differential_vars === nothing
    error("differential_vars must be set for DAE initialization to occur. Either set consistent initial conditions, differential_vars, or use a different initialization algorithm.")
  end

  nlequation = @closure (out, x) -> begin
    @. du = ifelse(differential_vars,x,du)
    @. u  = ifelse(differential_vars,u,x)
    f(out, du, u, p, t)
  end

  r = nlsolve(nlequation, ifelse.(differential_vars,du,u))

  @. du = ifelse(differential_vars,r.zero,du)
  @. u  = ifelse(differential_vars,u,r.zero)

  recursivecopy!(integrator.uprev,integrator.u)
  if alg_extrapolates(integrator.alg)
    recursivecopy!(integrator.uprev2,integrator.uprev)
  end

  return
end

function _initialize_dae!(integrator, prob::DAEProblem,
              alg::BrownFullBasicInit, ::Val{false})
  @unpack p, t, f = integrator
  differential_vars = prob.differential_vars

  if integrator.opts.internalnorm(f(integrator.du, integrator.u, p, t),t) <= alg.abstol
    return
  elseif differential_vars === nothing
    error("differential_vars must be set for DAE initialization to occur. Either set consistent initial conditions, differential_vars, or use a different initialization algorithm.")
  end

  if integrator.u isa Number && integrator.du isa Number
    # This doesn't fix static arrays!
    u = [integrator.u]
    du = [integrator.du]
  else
    u = integrator.u
    du = integrator.du
  end

  nlequation = @closure (out,x) -> begin
    @. du = ifelse(differential_vars,x,du)
    @. u  = ifelse(differential_vars,u,x)
    out .= f(du, u, p, t)
  end

  r = nlsolve(nlequation, ifelse.(differential_vars,du,u))

  @. du = ifelse(differential_vars,r.zero,du)
  @. u  = ifelse(differential_vars,u,r.zero)

  if integrator.u isa Number && integrator.du isa Number
    # This doesn't fix static arrays!
    integrator.u = first(u)
    integrator.du = first(du)
  else
    integrator.u = u
    integrator.du = du
  end

  integrator.uprev = integrator.u
  if alg_extrapolates(integrator.alg)
    integrator.uprev2 = integrator.uprev
  end

  return
end


text:
push!(LOAD_PATH, "../src/", "../test/src/")

DOCUMENTER_DEBUG=true

using Documenter, PartialSvdStoch

makedocs(
    format = Documenter.HTML(prettyurls = false),
    sitename = "PartialSvdStoch",
    pages = Any[
        "Introduction" => "INTRO.md",
        "PartialSvdStoch.jl " => "index.md",
        "Tests" => "Test.md"
    ]
)


text:
# Parts of this code was written by @jpfairbanks.

# Parallel Breadth-first search / traversal using a frontier based parallelized
# approach.

#################################################
#
# Parallel frontier based Breadth-first search approach
#
#################################################

using Base.Threads

import Base: push!, shift!, isempty, getindex

export bfs_tree, LevelSynchronousBFS

mutable struct LevelSynchronousBFS <: AbstractGraphVisitAlgorithm end

"""
    ThreadQueue

A thread safe queue implementation for using as the queue for BFS.
"""

immutable ThreadQueue{T, N<:Integer}
    data::Vector{T}
    head::Atomic{N} #Index of the head
    tail::Atomic{N} #Index of the tail
end

function ThreadQueue(T::Type, maxlength::N) where N <: Integer
    q = ThreadQueue(Vector{T}(maxlength), Atomic{N}(1), Atomic{N}(1))
    return q
end

function push!{T, N}(q::ThreadQueue{T, N}, val::T)
    # TODO: check that head > tail
    offset = atomic_add!(q.tail, one(N))
    q.data[offset] = val
    return offset
end

function shift!{T, N}(q::ThreadQueue{T, N})
    # TODO: check that head < tail
    offset = atomic_add!(q.head, one(N))
    return q.data[offset]
end

function isempty{T, N}(q::ThreadQueue{T, N})
    return ( q.head[] == q.tail[] ) && q.head != one(N)
end

function getindex{T}(q::ThreadQueue{T}, iter)
    return q.data[iter]
end

# Traverses the vertices in the queue and adds newly found successors to the queue.
function bfskernel{T <: Integer}(
        alg::LevelSynchronousBFS,
        next::ThreadQueue, # Thread safe queue to add vertices to
        g::AbstractGraph, # The graph
        parents::Array{Atomic{T}}, # Parents array
        level::Array{T} # Vertices in the current frontier
    )
    @threads for src in level
        vertexneighbors = neighbors(g, src) # Get the neighbors of the vertex
        for vertex in vertexneighbors
            # Atomically check and set parent value if not set yet.
            parent = atomic_cas!(parents[vertex], zero(T), src)
            if parent==0
                push!(next, vertex) #Push onto queue if newly found
            end
        end
    end
end

"""
    bfs_tree!(LevelSynchronousBFS(), g, src, parents)

Provide a parallel breadth-first traversal of the graph `g` starting with source vertex `s`,
and return a parents array. The returned array is an Array of `Atomic` integers.

### Implementation Notes
This function uses `@threads` for parallelism which depends on the `JULIA_NUM_THREADS`
environment variable to decide the number of threads to use. Refer `@threads` documentation
for more details.
"""
function bfs_tree!(
        alg::LevelSynchronousBFS,
        next::ThreadQueue, # Thread safe queue to add vertices to
        g::AbstractGraph, # The graph
        source::T, # Source vertex
        parents::Array{Atomic{T}} # Parents array
    ) where T<:Integer
    parents[source][]=source # Set source to source
    push!(next, source) # Add source to the queue
    while !isempty(next)
        level = next[next.head[]:next.tail[]-1] # Get vertices in the frontier
        next.head[] = next.tail[] # reset the queue
        bfskernel(alg, next, g, parents, level) # Find new frontier
    end
    return parents
end

"""
    bfs_tree(LevelSynchronousBFS(), g, s, nv)

Provide a parallel breadth-first traversal of the graph `g` starting with source vertex `s`,
and return a directed acyclic graph of vertices in the order they were discovered
using a frontier based parallel approach.

### Implementation Notes
This function uses `@threads` for parallelism which depends on the `JULIA_NUM_THREADS`
environment variable to decide the number of threads to use. Refer `@threads` documentation
for more details.
This function is a high level wrapper around [`bfs_tree!`](@ref); use that function for more performance.
"""
function bfs_tree{T <: Integer}(alg::LevelSynchronousBFS, g::AbstractGraph, source::T, nv::T)
    next = ThreadQueue(T, nv) # Initialize threadqueue
    parents = [Atomic{T}(0) for i=1:nv] # Create parents array
    bfs_tree!(alg, next, g, source, parents)
    tree([i[] for i in parents])
end

function bfs_tree{T <: Integer}(alg::LevelSynchronousBFS, g::AbstractGraph, source::T)
    nvg = nv(g)
    bfs_tree(alg, g, source, nvg)
end

max_stars_repo_name: AsileBcd/JohnsonCycles

max_issues_repo_name: AsileBcd/JohnsonCycles

max_forks_repo_name: AsileBcd/JohnsonCycles


text:
using CUDA
using Adapt
using StructArrays
using Images
#using ImageView
using Random
using Test
include("../util/Matrix_utils.jl")
include("../util/Enum_Collection.jl")
include("../ui/Window.jl")
include("../ui/Button.jl")
include("../ui/Drawables.jl")

struct Jadx_Environment{A,B,C,D}
    all_windows::A
    all_buttons::B
    windows::C
    frame_buffer::D
end
Adapt.@adapt_structure Jadx_Environment

function init_components()
end

function init_app_close_button(number_individuals)
    matrix_unclicked = CuArray(get_array_of_image("close_window_button_large_unclicked.png"))
    return Button(380,1,size(matrix_unclicked,3),size(matrix_unclicked,2),2,matrix_unclicked,matrix_unclicked,CUDA.fill(false,(2,number_individuals)),true,placeholder)
end

function init_dropdown_button_datei(position,number_individuals)
    matrix_unclicked = CuArray(get_array_of_image("drpdwn_datei_unclicked.png"))
    matrix_clicked = CuArray(get_array_of_image("drpdwn_datei_clicked.png"))
    return Button(position,12,size(matrix_unclicked,3),size(matrix_unclicked,2),2,matrix_clicked,matrix_unclicked,CUDA.fill(false,(2,number_individuals)),true,placeholder)
end
function init_dropdown_button_anzeigen(position,number_individuals)
    matrix_unclicked = CuArray(get_array_of_image("drpdwn_anzeigen_unclicked.png"))
    matrix_clicked = CuArray(get_array_of_image("drpdwn_anzeigen.png"))
    return Button(position,12,size(matrix_unclicked,3),size(matrix_unclicked,2),2,matrix_clicked,matrix_unclicked,CUDA.fill(false,(2,number_individuals)),true,placeholder)
end
function init_dropdown_button_navigation(position,number_individuals)
    matrix_unclicked = CuArray(get_array_of_image("drpdwn_navigation_unclicked.png"))
    matrix_clicked = CuArray(get_array_of_image("drpdwn_navigation_clicked.png"))
    return Button(position,12,size(matrix_unclicked,3),size(matrix_unclicked,2),2,matrix_clicked,matrix_unclicked,CUDA.fill(false,(2,number_individuals)),true,placeholder)
end
function init_dropdown_button_tools(position,number_individuals)
    matrix_unclicked = CuArray(get_array_of_image("drpdwn_tools_unclicked.png"))
    matrix_clicked = CuArray(get_array_of_image("drpdwn_tools_clicked.png"))
    return Button(position,12,size(matrix_unclicked,3),size(matrix_unclicked,2),2,matrix_clicked,matrix_unclicked,CUDA.fill(false,(2,number_individuals)),true,placeholder)
end
function init_dropdown_button_hilfe(position,number_individuals)
    matrix_unclicked = CuArray(get_array_of_image("drpdwn_hilfe_unclicked.png"))
    matrix_clicked = CuArray(get_array_of_image("drpdwn_hilfe_clicked.png"))
    return Button(position,12,size(matrix_unclicked,3),size(matrix_unclicked,2),2,matrix_clicked,matrix_unclicked,CUDA.fill(false,(2,number_individuals)),true,placeholder)
end

function init_main_window_small_buttons(number_individuals)
    position=1
    matrix_unclicked = CuArray(get_array_of_image("small_button_1_unclicked.png"))
    matrix_clicked = CuArray(get_array_of_image("small_button_1_clicked.png"))
    small_button_1 = Button(position,22,size(matrix_unclicked,3),size(matrix_unclicked,2),2,matrix_clicked,matrix_unclicked,CUDA.fill(false,(2,number_individuals)),true,placeholder)
    position += 7
    matrix_unclicked = CuArray(get_array_of_image("small_button_2_unclicked.png"))
    matrix_clicked = CuArray(get_array_of_image("small_button_2_clicked.png"))
    small_button_2 = Button(position,22,size(matrix_unclicked,3),size(matrix_unclicked,2),2,matrix_clicked,matrix_unclicked,CUDA.fill(false,(2,number_individuals)),true,placeholder)
    position += 7
    matrix_unclicked = CuArray(get_array_of_image("small_button_3_unclicked.png"))
    matrix_clicked = CuArray(get_array_of_image("small_button_3_clicked.png"))
    small_button_3 = Button(position,22,size(matrix_unclicked,3),size(matrix_unclicked,2),2,matrix_clicked,matrix_unclicked,CUDA.fill(false,(2,number_individuals)),true,placeholder)
    position += 7
    matrix_unclicked = CuArray(get_array_of_image("small_button_4_unclicked.png"))
    matrix_clicked = CuArray(get_array_of_image("small_button_4_clicked.png"))
    small_button_4 = Button(position,22,size(matrix_unclicked,3),size(matrix_unclicked,2),2,matrix_clicked,matrix_unclicked,CUDA.fill(false,(2,number_individuals)),true,placeholder)
    position += 7
    matrix_unclicked = CuArray(get_array_of_image("small_button_5_unclicked.png"))
    matrix_clicked = CuArray(get_array_of_image("small_button_5_clicked.png"))
    small_button_5 = Button(position,22,size(matrix_unclicked,3),size(matrix_unclicked,2),2,matrix_clicked,matrix_unclicked,CUDA.fill(false,(2,number_individuals)),true,placeholder)
    position += 7
    matrix_unclicked = CuArray(get_array_of_image("small_button_6_unclicked.png"))
    matrix_clicked = CuArray(get_array_of_image("small_button_6_clicked.png"))
    small_button_6 = Button(position,22,size(matrix_unclicked,3),size(matrix_unclicked,2),2,matrix_clicked,matrix_unclicked,CUDA.fill(false,(2,number_individuals)),true,placeholder)
    position += 7
    matrix_unclicked = CuArray(get_array_of_image("small_button_7_unclicked.png"))
    matrix_clicked = CuArray(get_array_of_image("small_button_7_clicked.png"))
    small_button_7 = Button(position,22,size(matrix_unclicked,3),size(matrix_unclicked,2),2,matrix_clicked,matrix_unclicked,CUDA.fill(false,(2,number_individuals)),true,placeholder)
    position += 7
    matrix_unclicked = CuArray(get_array_of_image("small_button_8_unclicked.png"))
    matrix_clicked = CuArray(get_array_of_image("small_button_8_clicked.png"))
    small_button_8 = Button(position,22,size(matrix_unclicked,3),size(matrix_unclicked,2),2,matrix_clicked,matrix_unclicked,CUDA.fill(false,(2,number_individuals)),true,placeholder)
    position += 7
    matrix_unclicked = CuArray(get_array_of_image("small_button_9_unclicked.png"))
    matrix_clicked = CuArray(get_array_of_image("small_button_9_clicked.png"))
    small_button_9 = Button(position,22,size(matrix_unclicked,3),size(matrix_unclicked,2),2,matrix_clicked,matrix_unclicked,CUDA.fill(false,(2,number_individuals)),true,placeholder)
    position += 7
    matrix_unclicked = CuArray(get_array_of_image("small_button_10_unclicked.png"))
    matrix_clicked = CuArray(get_array_of_image("small_button_10_clicked.png"))
    small_button_10 = Button(position,22,size(matrix_unclicked,3),size(matrix_unclicked,2),2,matrix_clicked,matrix_unclicked,CUDA.fill(false,(2,number_individuals)),true,placeholder)
    position += 7
    matrix_unclicked = CuArray(get_array_of_image("small_button_11_unclicked.png"))
    matrix_clicked = CuArray(get_array_of_image("small_button_11_clicked.png"))
    small_button_11 = Button(position,22,size(matrix_unclicked,3),size(matrix_unclicked,2),2,matrix_clicked,matrix_unclicked,CUDA.fill(false,(2,number_individuals)),true,placeholder)
    position += 7
    matrix_unclicked = CuArray(get_array_of_image("small_button_12_unclicked.png"))
    matrix_clicked = CuArray(get_array_of_image("small_button_12_clicked.png"))
    small_button_12 = Button(position,22,size(matrix_unclicked,3),size(matrix_unclicked,2),2,matrix_clicked,matrix_unclicked,CUDA.fill(false,(2,number_individuals)),true,placeholder)
    position += 7
    matrix_unclicked = CuArray(get_array_of_image("small_button_13_unclicked.png"))
    matrix_clicked = CuArray(get_array_of_image("small_button_13_clicked.png"))
    small_button_13 = Button(position,22,size(matrix_unclicked,3),size(matrix_unclicked,2),2,matrix_clicked,matrix_unclicked,CUDA.fill(false,(2,number_individuals)),true,placeholder)

    return begin
        small_button_1,
        small_button_2,
        small_button_3,
        small_button_4,
        small_button_5,
        small_button_6,
        small_button_7,
        small_button_8,
        small_button_9,
        small_button_10,
        small_button_11,
        small_button_12,
        small_button_13
    end
end

function init_all_buttons(number_individuals)
    close_button = init_app_close_button(number_individuals)
    position = 1
    button_datei = init_dropdown_button_datei(position,number_individuals)
    position += button_datei.width
    button_anzeigen = init_dropdown_button_anzeigen(position,number_individuals)
    position += button_anzeigen.width
    button_navigation = init_dropdown_button_navigation(position,number_individuals)
    position += button_navigation.width
    button_tools = init_dropdown_button_tools(position,number_individuals)
    position += button_tools.width
    button_hilfe = init_dropdown_button_hilfe(position,number_individuals)
    begin
        small_button_1,
        small_button_2,
        small_button_3,
        small_button_4,
        small_button_5,
        small_button_6,
        small_button_7,
        small_button_8,
        small_button_9,
        small_button_10,
        small_button_11,
        small_button_12,
        small_button_13 = init_main_window_small_buttons(number_individuals)
    end
    return All_Buttons(
        close_button,
        button_datei,
        button_anzeigen,
        button_navigation,
        button_tools,
        button_hilfe,
        small_button_1,
        small_button_2,
        small_button_3,
        small_button_4,
        small_button_5,
        small_button_6,
        small_button_7,
        small_button_8,
        small_button_9,
        small_button_10,
        small_button_11,
        small_button_12,
        small_button_13
        )
end

function init_main_window(number_individuals)
    windows = CuArray{Window_Names}(undef,(5,number_individuals))
    buttons = Array{Button_Names}(undef,19)
    buttons[1] = app_close_button
    buttons[2] = dropdown_button_datei
    buttons[3] = dropdown_button_anzeigen
    buttons[4] = dropdown_button_navigation
    buttons[5] = dropdown_button_tools
    buttons[6] = dropdown_button_hilfe
    buttons[7] = small_button_1
    buttons[8] = small_button_2
    buttons[9] = small_button_3
    buttons[10] = small_button_4
    buttons[11] = small_button_5
    buttons[12] = small_button_6
    buttons[13] = small_button_7
    buttons[14] = small_button_8
    buttons[15] = small_button_9
    buttons[16] = small_button_10
    buttons[17] = small_button_11
    buttons[18] = small_button_12
    buttons[19] = small_button_13
    matrix_self = CuArray(get_array_of_image("main_window.png"))
    current_matrix = CuArray{N0f8}(undef,(size(matrix_self,1),size(matrix_self,2),size(matrix_self,3),number_individuals))

    return Window(matrix_self,current_matrix,1,1,size(matrix_self,3),size(matrix_self,2),windows,CuArray(buttons),true,false)
end

function init_all_windows(number_individuals)
    main = init_main_window(number_individuals)
    return All_Windows(main)
end

function initialize()
    number_individuals = 100

    all_buttons = init_all_buttons(number_individuals)
    #init All Buttons

    all_windows = init_all_windows(number_individuals)
    #init all Windows
    windows = CuArray{Window_Names}(undef,(5,number_individuals))
    #frame_buffer = similar(all_windows.main_window.matrix_self)
    frame_buffer = CuArray{N0f8}(undef,(size(all_windows.main_window.matrix_self,1),size(all_windows.main_window.matrix_self,2),size(all_windows.main_window.matrix_self,3),number_individuals))
    return Jadx_Environment(all_windows,all_buttons,windows,frame_buffer)
end
#display(env)

function reset(env::Jadx_Environment)
    for i in 1:size(env.windows,1)
        env.windows[i,blockIdx().x] = no_window
    end
    env.windows[1,blockIdx().x] = main_window
    sync_threads()
    current_window = get_window(env.windows[1,blockIdx().x],env.all_windows)

    draw_self(current_window,env.all_windows,env.all_buttons)

    sync_threads()

    kernel_blit_image_inplace(threadIdx().x,blockIdx().x,env.frame_buffer,current_window.current_matrix,current_window.x_Coord,current_window.y_Coord)

end

function step(env::Jadx_Environment,action,input)
    reward = 0
    number_of_windows_to_be_removed  = 0
    should_restack = false
    click_x = clamp(ceil(action[1] * 400),0,400)
    click_y = clamp(ceil(action[2] * 268),0,268)
    for i in 1:size(env.windows,1)
        if env.windows[i,blockIdx().x] ≠ no_window
            #index = i - windows_to_be_removed
            #last_clicked_index = size(env.windows,1) - index
            current_window = get_window(env.windows[i,blockIdx().x],env.all_windows)

            reward,window_includes_point,coords_x,coords_y = click(click_x,click_y,current_window,env.all_buttons)
            sync_threads()
            if window_includes_point
                if !should_restack
                    draw_self(current_window,env.all_windows,env.all_buttons)
                    sync_threads()
                    kernel_blit_image_inplace(threadIdx().x,blockIdx().x,env.frame_buffer,current_window.current_matrix,coords_x,coords_y)
                end
                break

            else
                if current_window.modal
                    break
                elseif current_window.autoclose

                end
            end
        end
    end
    if should_restack
        #stack_windows
    end
    should_restack = false
    sync_threads()
    grayscale(input,env.frame_buffer)
    sync_threads()
    return reward
end

function kernel3(env::Jadx_Environment,input)
    action = @cuDynamicSharedMem(Float32,2)
    fitness_total = 0
    for round in 1:1

        reset(env)
        for i in 1:1
            if threadIdx().x == 1
                action[1] = rand(Float32)
                action[2] = rand(Float32)
            end
            sync_threads()
            reward = step(env,action,input)
        end
    end
    return

end
env = initialize()
width = env.all_windows.main_window.width
height = env.all_windows.main_window.height
input = CUDA.fill(0.0f0,(268,400,100))
@cuda threads=400 blocks=100 shmem=sizeof(Float32)*2 kernel3(env,input)
CUDA.synchronize()

#display(env.windows)
#imshow(colorview(Gray,Array(input[:,:,1])))

#imshow(colorview(RGB,Array(env.frame_buffer[:,:,:,1])))

#imshow(colorview(RGB,Array(env.frame_buffer[:,:,:,2])))
#imshow(colorview(RGB,Array(env.frame_buffer[:,:,:,3])))
#imshow(colorview(RGB,Array(env.frame_buffer[:,:,:,4])))
#imshow(colorview(RGB,Array(env.frame_buffer[:,:,:,5])))
#imshow(colorview(RGB,Array(env.frame_buffer[:,:,:,6])))
#imshow(colorview(RGB,Array(env.frame_buffer[:,:,:,7])))
#imshow(colorview(RGB,Array(env.frame_buffer[:,:,:,8])))
#imshow(colorview(RGB,Array(env.frame_buffer[:,:,:,9])))
#imshow(colorview(RGB,Array(env.frame_buffer[:,:,:,10])))

function episode(number_time_steps,number_rounds,environment_cfg::Jadx_Environment)
    action = @cuDynamicSharedMem(Float32,2)
    fitness_total = 0
    #brain_initialize
    for round in 1:number_rounds
        fitness_current = 0
        #reset brain
        #env_initialize()
        obs = reset(environment_cfg)
        sync_threads()
        for timestep in 1:number_time_steps
            #action = brain_step(obs)
            if threadIdx().x == 1
                action[1] = rand(Float32)
                action[2] = rand(Float32)
                @cuprintln()
            end
            sync_threads()
            obs,reward = env_step(environment_cfg,action)
        end

    end
end

#@cuda threads=400 blocks=5 shmem=sizeof(Float32)*2 episode(1000,1,env)


text:
module RoME

using Reexport

@reexport using IncrementalInference
@reexport using TransformUtils
@reexport using Distributions
@reexport using KernelDensityEstimate
@reexport using ApproxManifoldProducts

using
  Distributed,
  LinearAlgebra,
  Statistics,
  Graphs,  # TODO determine how many parts still require Graphs still directly
  Rotations,
  CoordinateTransformations,
  JLD2,
  ProgressMeter,
  DocStringExtensions

import Base: +, \, convert
import TransformUtils: ⊖, ⊕, convert, compare, ominus, veeQuaternion
import IncrementalInference: convert, getSample, reshapeVec2Mat, extractdistribution, DFG

# const AMP = ApproxManifoldProducts

export
  IIF,  # Aliases for various packages
  KDE,
  TU,
  AMP,
  # initfg,
  # RoME specific functions
  measureMeanDist,
  predictBodyBR,
  getLastPose,
  getLastPose2D,
  odomKDE,
  initFactorGraph!,
  addOdoFG!,
  addposeFG!,
  newLandm!,
  addBRFG!,
  addMMBRFG!,
  addAutoLandmBR!,
  projNewLandm!,
  malahanobisBR,
  veePose3,
  veePose,
  \,
  RangeAzimuthElevation,

  # helper functions
  get2DSamples,
  # getAll2D,
  get2DSampleMeans,
  getAll2DMeans,
  getAll2DPoses,
  get2DPoseSamples,
  get2DPoseMeans,
  getKDE,
  getVertKDE,
  get2DPoseMax,
  getAll2DLandmarks,
  get2DLandmSamples,
  get2DLandmMeans,
  get2DLandmMax,

  # helper functions
  getLastLandm2D,
  getLastPose2D,
  getNextLbl,

  # RobotUtils
  getRangeKDEMax2D,

  # some transform functions
  cart2pol,
  pol2cart,

  # Feature tracking code
  Feature,
  initTrackersFrom,
  propAllTrackers!,
  measUpdateTrackers!,
  assocMeasWFeats!,

  lsrBR,

  # Didson model
  evalPotential,
  LinearRangeBearingElevation,
  project!,
  project,
  backprojectRandomized!,
  residual!,
  residualLRBE!,
  reuseLBRA,
  ominus,
  ominus!,
  +,
  evalPotential,
  getSample!,
  getSample,
  # obsolete
  WrapParam,
  WrapParamArray,

  # Didson convenience function
  addLinearArrayConstraint,

  # camera model -- TODO --separate out
  CameraIntrinsic,
  CameraExtrinsic,
  CameraModelFull,
  project!,
  project,
  backprojectRandomized!,
  # keep
  cameraResidual!,

  # Point2D
  Point2,
  Point2Point2,
  PackedPoint2Point2,
  Point2Point2WorldBearing,
  PackedPoint2Point2WorldBearing,
  Point2Point2Range,
  PackedPoint2Point2Range,
  PriorPoint2,
  PackedPriorPoint2,
  Pose2Point2BearingRange,
  Pose2Point2BearingRangeMH,
  PackedPose2Point2BearingRange,
  PackedPose2Point2BearingRangeMH,
  Pose2Point2Bearing,
  PackedPose2Point2Bearing,
  Pose2Point2Range,
  PackedPose2Point2Range,
  PriorPoint2,
  PackedPriorPoint2,
  # Point2D with null hypotheses
  PriorPoint2DensityNH,
  PackedPriorPoint2DensityNH,

  # Velocity in Point2 types
  DynPoint2,
  DynPoint2VelocityPrior,
  DynPoint2DynPoint2,
  VelPoint2VelPoint2,
  Point2Point2Velocity,
  PackedDynPoint2VelocityPrior,
  PackedVelPoint2VelPoint2,

  # likely to be deprecated
  solveLandm,
  solvePose2,
  solveSetSeps,
  addPose2Pose2!,

  # acoustics
  Pose2Point2BearingRangeDensity,
  PackedPose2Point2BearingRangeDensity,
  Pose2Point2RangeDensity,
  PackedPose2Point2RangeDensity,

  # Pose2D
  Pose2,
  PriorPose2,
  PackedPriorPose2,
  PartialPriorYawPose2,
  PackedPartialPriorYawPose2,
  Pose2Pose2,
  PackedPose2Pose2,
  # velocity in Pose2
  DynPose2,
  DynPose2VelocityPrior,
  PackedDynPose2VelocityPrior,
  VelPose2VelPose2,
  PackedVelPose2VelPose2,
  DynPose2Pose2,
  PackedDynPose2Pose2,
  # Will be deprecated
  addPose2Pose2,

  # Polar types
  Polar,
  PolarPolar,
  PriorPolar,

  # MultipleFeatures2D constraint functions
  MultipleFeatures2D,
  getUvecScaleFeature2D,
  getUvecScaleBaseline2D,

  # Pose3, Three dimensional
  Pose3,
  Point3,
  PriorPoint3,
  PriorPose3,
  PackedPriorPose3,
  Pose3Pose3,
  PackedPose3Pose3,
  projectParticles,
  ⊕,
  Pose3Pose3NH,
  PackedPose3Pose3NH,

  # partial Pose3
  PriorPose3ZRP,
  Pose3Pose3XYYaw,
  # Various utilities
  passTypeThrough,

  # SLAM specific functions
  SLAMWrapper,

  # FG Analysis tools
  rangeErrMaxPoint2,
  rangeCompAllPoses,
  rangeCompAllPoses,

  # new robot navigation functionality
  triggerPose,
  GenericInSituSystem,
  InSituSystem,
  makeInSituSys,
  makeGenericInSituSys,
  advOdoByRules,
  poseTrigAndAdd!,
  poseTrigAndAdd!,
  processTreeTrackersUpdates!,
  addSoftEqualityPoint2D,
  vectoarr2,
  basicFactorGraphExample,
  predictVariableByFactor,

  # jld required Features Type
  LaserFeatures,

  # Deprecated
  PartialPriorRollPitchZ,
  PackedPartialPriorRollPitchZ,
  PartialPose3XYYaw,
  PackedPartialPose3XYYaw,
  PartialPose3XYYawNH,
  PackedPartialPose3XYYawNH
  # Point2DPoint2DRange,
  # PackedPoint2DPoint2DRange,
  # PackedPoint2DPoint2DRange,
  # PackedPriorPoint2D,
  # Pose2DPoint2DRangeDensity, # to be deprecated
  # Point2DPoint2D, # deprecated
  # Point2DPoint2DRange, # deprecated
  # PriorPoint2D, # deprecated
  # Pose2DPoint2DBearingRange, # begin deprecated
  # Pose2DPoint2DBearing, # deprecated
  # PriorPoint2D, # deprecated
  # PackedPriorPoint2D # deprecated`

# doesnt seem to work
# @info "Setting IncrementalInference de-serialization namespace RoME"
# setSerializationNamespace!("RoME" => RoME)
# @info "done..."

  # # solve with isam in pytslam
  # doISAMSolve,
  # drawCompPosesLandm,
  #
  # # Victoria Park data specific
  # addLandmarksFactoGraph!,
  # appendFactorGraph!,
  # doBatchRun,
  # rotateFeatsToWorld

include("SpecialDefinitions.jl")

include("BayesTracker.jl")

include("SensorModels.jl")
include("CameraModel.jl")

# 2D
include("variables/Point2D.jl")
include("variables/Pose2D.jl")
include("variables/DynPoint2D.jl")
include("variables/DynPose2D.jl")

# 3D
include("variables/Point3D.jl")
include("variables/Pose3D.jl")

include("factors/Point2D.jl")
include("factors/Polar.jl")
include("factors/Pose2D.jl")
include("factors/Bearing2D.jl")
include("factors/Range2D.jl")
include("factors/BearingRange2D.jl")
include("factors/DynPoint2D.jl")
include("factors/DynPose2D.jl")
include("factors/Point3D.jl")
include("factors/Pose3Pose3.jl")
include("factors/PartialPose3.jl")
include("factors/MultipleFeaturesConstraint.jl")
include("factors/InertialPose3.jl")

include("Slam.jl")

include("RobotUtils.jl")

include("SimulationUtils.jl")

include("FactorGraphAnalysisTools.jl")

include("RobotDataTypes.jl") #WheeledRobotUtils
include("NavigationSystem.jl")

include("Deprecated.jl")

# include("dev/ISAMRemoteSolve.jl")

end


text:
using MonteCarlo, Distributions, PyPlot, DataFrames, JLD

Tdist = Normal(MonteCarlo.IsingTc, .64)
n_Ts = 2^8 # 2^8
Ts = sort!(rand(Tdist, n_Ts))
Ts = Ts[Ts.>=1.2]
Ts = Ts[Ts.<=3.8]
therm = 10^4 # 10^4
sweeps = 10^3 # 10^3

df = DataFrame(L=Int[], T=Float64[], M=Float64[], χ=Float64[], E=Float64[], C_V=Float64[])

for L in 2 .^ [3, 4, 5, 6] #[3, 4, 5, 6]
	println("L = ", L)
	for (i, T) in enumerate(Ts)
		println("\t T = ", T)
		beta = 1/T
		model = IsingModel(dims=2, L=L)
		mc = MC(model, beta=beta)
		# mc.p.global_moves = true # enable Wolff cluster
		# mc.p.global_rate = 1
		run!(mc, sweeps=sweeps, thermalization=therm, verbose=false)
		meas = MonteCarlo.measurements(mc)[:ME]
		push!(df, [L, T, mean(meas[:Magn].m), mean(meas[:Magn].chi), mean(meas[:Energy].e), mean(meas[:Energy].C)])
	end
	flush(stdout)
end

sort!(df, [:L, :T])
@save "ising2d.jld" df

grps = groupby(df, :L)
# plot results individually
for g in grps
	L = g[:L][1]
	fig, ax = subplots(2,2, figsize=(12,8))
	ax[1].plot(g[:T], g[:E], "o", color="darkred", markeredgecolor="black")
	ax[1].set_ylabel("Energy")
	ax[1].set_xlabel("Temperature")

	ax[2].plot(g[:T], g[:C_V], "o", color="darkred", markeredgecolor="black")
	ax[2].set_ylabel("Specific heat")
	ax[2].set_xlabel("Temperature")

	ax[3].plot(g[:T], g[:M], "o", color="C0", markeredgecolor="black")
	ax[3].axvline(x=MonteCarlo.IsingTc, color="black", linewidth=2.0, label="\$ T_c \$")
	# ax[3].legend(loc="best")
	ax[3].set_ylabel("Magnetization")
	ax[3].set_xlabel("Temperature")

	ax[4].plot(g[:T], g[:χ], "o", color="C0", markeredgecolor="black")
	ax[4].axvline(x=MonteCarlo.IsingTc, color="black", linewidth=2.0, label="\$ T_c \$")
	# ax[4].legend(loc="best")
	ax[4].set_ylabel("Susceptibility χ")
	ax[4].set_xlabel("Temperature")
	tight_layout()
	savefig("ising2d_L_$(L).pdf")
end

# plot results together
fig, ax = subplots(2,2, figsize=(12,8))

for g in grps
	L = g[:L][1]
	ax[1].plot(g[:T], g[:E], "o", markeredgecolor="black", label="L=$L")
	ax[2].plot(g[:T], g[:C_V], "o", markeredgecolor="black", label="L=$L")
	ax[3].plot(g[:T], g[:M], "o", markeredgecolor="black", label="L=$L")
	ax[4].plot(g[:T], g[:χ], "o", markeredgecolor="black", label="L=$L")
end
ax[1].legend(loc="best")
ax[1].set_ylabel("Energy")
ax[1].set_xlabel("Temperature")

ax[2].set_ylabel("Specific heat")
ax[2].set_xlabel("Temperature")
ax[2].axvline(x=MonteCarlo.IsingTc, color="black", linestyle="dashed", label="\$ T_c \$")
ax[2].legend(loc="best")

ax[3].set_ylabel("Magnetization")
ax[3].set_xlabel("Temperature")
# ax[3].axvline(x=MonteCarlo.IsingTc, linewidth=2.0, color="black", label="\$ T_c \$")
x = range(1.2, stop=MonteCarlo.IsingTc, length=100)
y = (1 .- sinh.(2.0 ./ (x)).^(-4)).^(1/8)
ax[3].plot(x,y, "k--", label="exact")
ax[3].plot(range(MonteCarlo.IsingTc, stop=3.8, length=100), zeros(100), "k--")
ax[3].legend(loc="best")

ax[4].set_ylabel("Susceptibility χ")
ax[4].set_xlabel("Temperature")
ax[4].axvline(x=MonteCarlo.IsingTc, color="black", linestyle="dashed", label="\$ T_c \$")
ax[4].legend(loc="best")
tight_layout()
savefig("ising2d.pdf")
savefig("ising2d.svg")


text:
#! /usr/bin/julia

# Rosetta Code, Stem-and-leaf plot

function stemleaf{T<:Real}(a::Array{T,1}, leafsize=1)
    ls = 10^int(log10(leafsize))
    (stem, leaf) = divrem(sort(int(a/ls)), 10)
    leaf[sign(stem) .== -1] *= -1
    negzero = leaf .< 0
    if any(negzero)
        leaf[negzero] *= -1
        nz = @sprintf "%10s | " "-0"
        nz *= join(map(string, leaf[negzero]), " ")
        nz *= "\n"
        stem = stem[!negzero]
        leaf = leaf[!negzero]
    else
        nz = ""
    end
    slp = ""
    for i in stem[1]:stem[end]
        i != 0 || (slp *= nz)
        slp *= @sprintf "%10d | " i
        slp *= join(map(string, leaf[stem .== i]), " ")
        slp *= "\n"
    end
    slp *= " Leaf Unit = " * string(convert(T, ls)) * "\n"
    return slp
end

println("Using the Task's Test Data")
test = """12 127 28 42 39 113 42 18 44 118 44 37 113 124 37 48 127 36 29
   31 125 139 131 115 105 132 104 123 35 113 122 42 117 119 58 109 23 105
   63 27 44 105 99 41 128 121 116 125 32 61 37 127 29 113 121 58 114 126
   53 114 96 25 109 7 31 141 46 13 27 43 117 116 27 7 68 40 31 115 124 42
   128 52 71 118 117 38 27 106 33 117 116 111 40 119 47 105 57 122 109
   124 115 43 120 43 27 27 18 28 48 125 107 114 34 133 45 120 30 127 31
   116 146"""
test = map(parseint, split(test, r"\s"))
println(stemleaf(test))

println("Test with Reals and Negative Zero Stem")
test = [-23.678758, -12.45, -3.4, 4.43, 5.5, 5.678, 16.87, 24.7, 56.8]
println(stemleaf(test))

println("Test with Leaf Size Scaling")
test = int(500*randn(20))
println("Using:  ", test)
println(stemleaf(test, 10))

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code