text:
export matches, recursivematches, replace, find_matches, find_matches!, evaluate, associative_ordering, isequal_associative

function reconstruct_replacement(replacement_expr::AbstractExpression, replacements::Dict{LogicalSymbol, AbstractExpression})
    if istree(replacement_expr)
        rr = reconstruct_replacement.(arguments(replacement_expr), Iterators.repeat([replacements], length(arguments(replacement_expr))))
        return LogicalExpression(Vector{AbstractExpression}(rr), operation(replacement_expr))
    end

    return replacements[replacement_expr]
end

function Base.replace(expr::AbstractExpression, rule::Pair{T, U}) where {T <: AbstractExpression, U <: AbstractExpression}
    replacements = Dict{LogicalSymbol, AbstractExpression}()
    if !find_matches!(expr, rule.first, replacements)
        return expr
    end

    reconstruct_replacement(rule.second, replacements)
end

function recursivereplace(expr::AbstractExpression, rule::Pair{T, U}) where {T <: AbstractExpression, U <: AbstractExpression}
    replaced_expr = replace(expr, rule)
    if replaced_expr == expr && istree(replaced_expr)
        new_args = Vector{AbstractExpression}(recursivereplace.(arguments(expr), Iterators.repeat([rule], length(arguments(expr)))))
        if any(new_args .!= arguments(expr))
            return LogicalExpression(new_args, operation(expr))
        end
    end

    return replaced_expr
end

matches(expr::AbstractExpression, pattern::AbstractExpression) = find_matches!(expr, pattern, Dict{LogicalSymbol, AbstractExpression}())

function find_matches(expr::AbstractExpression, pattern::AbstractExpression)
    matches = Dict{LogicalSymbol, AbstractExpression}()
    find_matches!(expr, pattern, matches)

    matches
end
function find_matches!(expr::AbstractExpression, pattern::AbstractExpression, matches::Dict{LogicalSymbol, AbstractExpression})
    if istree(pattern)
        if !istree(expr) || operation(pattern) != operation(expr)
            return false
        end

        return all(find_matches!.(arguments(expr), arguments(pattern), Iterators.repeat([matches], length(arguments(expr)))))
    elseif pattern isa ExtensionalSet
        if !(expr isa ExtensionalSet) || cardinality(expr) != cardinality(pattern)
            return false
        end

        # return all(find_matches!.())
    end

    if haskey(matches, pattern)
        return isequal(matches[pattern], expr)
    else
        matches[pattern] = expr
    end
    true
end

function recursivematches(expr::AbstractExpression, pattern::AbstractExpression)
    if matches(expr, pattern)
        return true
    end
    return any(recursivematches.(arguments(expr), flat_repeat(pattern, length(arguments(expr)))))
end

substitute(sym::LogicalSymbol, symbol_values::Dict{LogicalSymbol, AbstractExpression}) = haskey(symbol_values, sym) ? symbol_values[sym] : sym
function substitute(expr::LogicalExpression, symbol_values::Dict{LogicalSymbol, AbstractExpression})
    return LogicalExpression(Vector{AbstractExpression}(substitute.(arguments(expr), flat_repeat(symbol_values, length(arguments(expr))))), operation(expr))
end

# simplification schemes
function repeated_chain_simplify(expr::AbstractExpression, rules::Vector{Pair{AbstractExpression, AbstractExpression}})
    modified_expression = expr

    i = 1
    while i <= length(rules)
        rule = rules[i]
        replacement_expression = replace(modified_expression, rule)
        i += 1
        if modified_expression != replacement_expression
            modified_expression = replacement_expression
            i = 1
        end
    end

    modified_expression
end

evaluate(sym::LogicalSymbol, values::Dict{LogicalSymbol, Bool}) = values[sym]
evaluate(expr::LogicalExpression, values::Dict{LogicalSymbol, Bool}=Dict{LogicalSymbol, Bool}()) = operation(expr).bool_fn((evaluate.(arguments(expr), Iterators.repeat([values], length(arguments(expr)))))...)

# Association
_isleft(x) = x == 1
_isright(x) = x == 2
_opposite_direction(x) = _isleft(x) ? 2 : 1
function associative_step(expr::AbstractExpression, leaf::AbstractExpression, direction)
    expr_copy = deepcopy(expr)
    associative_step!(expr_copy, leaf, direction)
end
function associative_step!(expr::AbstractExpression, leaf::AbstractExpression, direction)
    leaf_path = locate_expression(expr, leaf)

    descension_expr = nothing
    for i ∈ (length(leaf_path)):-1:1
        if leaf_path[i] != direction
            subpath = leaf_path[1:(i-1)]
            descension_expr = select_by_path(expr, subpath)
            break
        end
    end

    if isnothing(descension_expr)
        return expr
    end

    downwards_direction = _opposite_direction(direction)
    parent_node = descension_expr
    current_node = arguments(descension_expr)[direction]
    while istree(current_node)
        parent_node = current_node
        current_node = arguments(current_node)[downwards_direction]
    end

    associated_expr = LogicalExpression(_isleft(direction) ? AbstractExpression[current_node, leaf] : AbstractExpression[leaf, current_node], operation(expr))

    other_expr = select_by_path(expr, [leaf_path[1:end-1]..., _opposite_direction(leaf_path[end])])
    if !isequal(current_node, other_expr)

        t = select_by_path(expr, leaf_path[1:end-2])
        set_argument(parent_node, downwards_direction, associated_expr)
        if length(leaf_path[1:end-2]) == 0
            return arguments(expr)[direction]
        end
        set_argument(t, downwards_direction, other_expr)
    end
    @info expr

    expr
end


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

using Markdown
using InteractiveUtils

# ╔═╡ c65cb242-1f99-11ec-3f99-4fcf53839d59
using Pkg; Pkg.develop(path=".."); Pkg.activate("..")

# ╔═╡ 15c5ae4b-1910-4728-ae27-c3f8a601746e
begin
	using SqState
	using Plots
	using QuantumStateBase
	using JLD2
end

# ╔═╡ f51d1bc5-5ffc-4f6d-a5b1-cee378907dd1
readdir(joinpath(SqState.training_data_path(), "sqth"))

# ╔═╡ ec885764-37f6-4aa0-a399-6e0b20312ede
jldopen(joinpath(SqState.training_data_path(), "sqth", "sqth_2021-09-28T19_18_55.338.jld2"))["𝛒s"]

# ╔═╡ d65ed2f8-9d3a-4916-9874-a6c33676bb60
sum(isnothing.(jldopen(joinpath(SqState.training_data_path(), "sqth", "sqth_2021-09-28T19_19_33.958.jld2"))["𝛒s"]))

# ╔═╡ f4c98102-6b12-4f04-b46a-d02ef0c72882
σs = jldopen(joinpath(SqState.training_data_path(), "sqth", "sqth_2021-09-28T19_18_55.338.jld2"))["σs"]

# ╔═╡ c97e24d0-3c31-4624-8c88-6f182c0c161b
scatter(σs[:, 7])

# ╔═╡ Cell order:
# ╠═c65cb242-1f99-11ec-3f99-4fcf53839d59
# ╠═15c5ae4b-1910-4728-ae27-c3f8a601746e
# ╠═f51d1bc5-5ffc-4f6d-a5b1-cee378907dd1
# ╠═ec885764-37f6-4aa0-a399-6e0b20312ede
# ╠═d65ed2f8-9d3a-4916-9874-a6c33676bb60
# ╠═f4c98102-6b12-4f04-b46a-d02ef0c72882
# ╠═c97e24d0-3c31-4624-8c88-6f182c0c161b

ext: jl


text:
# Tracing charged particle in a static EM field.
#
# Hongyang Zhou, hyzhou@umich.edu

using TestParticle
using Meshes
using OrdinaryDiffEq
using Plots

## Initialize grid and field

x = range(-10, 10, length=15)
y = range(-10, 10, length=20)
z = range(-10, 10, length=25)
B = fill(0.0, 3, length(x), length(y), length(z)) # [T]
E = fill(0.0, 3, length(x), length(y), length(z)) # [V/m]

B[3,:,:,:] .= 1e-11
E[3,:,:,:] .= 5e-13

Δx = x[2] - x[1]
Δy = y[2] - y[1]
Δz = z[2] - z[1]

grid = CartesianGrid((length(x)-1, length(y)-1, length(z)-1),
   (x[1], y[1], z[1]),
   (Δx, Δy, Δz))

## Initialize particles

isAnalytic = false
trajectories = 1

x0 = [0.0, 0.0, 0.0] # initial position, [m]
u0 = [1.0, 0.0, 0.0] # initial velocity, [m/s]
stateinit = [x0..., u0...]

param_electron = prepare(grid, E, B, species=Electron)
tspan_electron = (0.0, 15.0)

param_proton = prepare(grid, E, B, species=Proton)
tspan_proton = (0.0, 10.0)

## Solve for the trajectories

prob_e = ODEProblem(trace_numeric!, stateinit, tspan_electron, param_electron)
prob_p = ODEProblem(trace_numeric!, stateinit, tspan_proton, param_proton)

sol_e = solve(prob_e, Tsit5(); save_idxs=[1,2,3])
sol_p = solve(prob_p, Tsit5(); save_idxs=[1,2,3])

## Visualization

xyz = plot(sol_e, vars=(1,2,3), xlabel="x", ylabel="y", zlabel="z",
   aspect_ratio=:equal, label="electron", lw=2)

plot!(sol_p, vars=(1,2,3), label="proton", lw=2, camera = (70,50))


text:
const MetaDiGraph = MetaGraph{<:Any,<:Any,<:DiGraph}

SimpleDiGraph(g::MetaDiGraph) = g.graph

is_directed(::Type{<:MetaDiGraph}) = true

function arrange(g::MetaDiGraph, label_1, label_2, code...)
    label_1, label_2
end


text:
module BrainFuck
using SparseArrays
export @bf_str, compile, interpret

"""

    bracketindices(tokens)

Returns a Dict of the indices of square brackets and their corresponding matching bracket.
"""
function bracketindices(tokens)
    stack = []
    bracketpairs = Dict{Int,Int}()
    for (index, token) in Iterators.enumerate(tokens)
        if token == '['
            push!(stack, index)
        elseif token == ']'
            if isempty(stack)
                throw(ArgumentError("Unmatched ']' at position $(index)."))
            end
            otherindex = pop!(stack)
            push!(bracketpairs, index => otherindex, otherindex => index)
        end
    end
    if !isempty(stack)
        throw(ArgumentError(
            "Unmatched '[' at position$(length(stack) > 1 ? "s" : "") $(join(stack, ", ", " and "))."
        ))
    end
    bracketpairs
end

"""

    clean(code::AbstractString)

Strips out all non-Brainfuck characters from the given string.
"""
clean(code::AbstractString) = replace(code, r"[^\-+,.<>[\]]" => "")

"""

    interpret(code; kwargs...)

Interpret the given string as Brainfuck `code`.

## Keyword arguments

- `memorysize::Integer` How many bytes to allocate for memory, default is `30000`.
- `input::Union{AbstractString,IO}` The string or buffer which the `,` command reads from.
    Defaults to `stdin`.
- `output::Union{IO,Type{String}}` The buffer which the `.` command write to. Returns the
    output as a string if `String` is passed as to `output`. Defaults to `stdout`.
"""
function interpret(code;
    memorysize::Integer = 30000,
    input::Union{AbstractString,IO} = stdin,
    output::Union{IO,Type{String}} = stdout
)
    tokens = clean(code)
    brackets = bracketindices(tokens)
    memory = spzeros(UInt8, memorysize)
    pointer = 1
    i = 1
    if input isa AbstractString
        input = IOBuffer(input)
    end
    returnstring = output isa Type{String}
    if returnstring
        output = IOBuffer()
    end
    while i <= length(tokens)
        token = tokens[i]
        if token == '-'
            memory[pointer] -= 0x01
        elseif token == '+'
            memory[pointer] += 0x01
        elseif token == '<'
            pointer = max(1, pointer - 1)
        elseif token == '>'
            pointer = min(memorysize, pointer + 1)
        elseif token == ','
            if !eof(input)
                skipchars(c -> c == 0x0d && peek(input) == 0x0a, input)
                memory[pointer] = read(input, UInt8)
            end
        elseif token == '.'
            print(output, Char(memory[pointer]))
        elseif token == '['
            if iszero(memory[pointer])
                i = brackets[i]
            end
        elseif token == ']'
            if !iszero(memory[pointer])
                i = brackets[i]
            end
        end
        i += 1
    end
    if returnstring
        String(take!(output))
    end
end

"""

    compile(code)

Compiles the BrainFuck code to a Julia function. The returned function has the same
keywords arguments as `interpret(code)`.
"""
compile(code) = eval(expr(code))

"""

    bf"+++++++++++++++++++++."

Defines a BrainFuck program.

## Example

```jldoctest
julia> helloworld = bf"\""
    ++++++++[>++++[>++>+++>+++>+<<<<-]>+>+>->>+[<]<-]>>.>
    ---.+++++++..+++.>>.<-.<.+++.------.--------.>>+.>++.
"\"";

julia> helloworld()
Hello World!

julia> helloworld(output = String)
"Hello World!\\n"

```
"""
macro bf_str(code)
    :(compile($code))
end

"""
Compiles the BrainFuck code to a Julia expr which will `eval` into a function.
"""
function expr(code)
    stack = [quote
        if input isa AbstractString
            input = IOBuffer(input)
        end
        returnstring = output isa Type{String}
        if returnstring
            output = IOBuffer()
        end
        pointer = 1
        memory = spzeros(UInt8, memorysize)
    end]
    tokens = Iterators.Stateful(clean(code))
    for token in tokens
        if occursin(token, "+-")
            i = token == '+' ? 0x01 : -0x01
            while occursin(peek(tokens), "+-")
                i += (popfirst!(tokens) == '+' ? 0x01 : -0x01)
            end
            push!(stack[end].args, :(memory[pointer] += $i))
        elseif occursin(token, "<>")
            i = token == '>' ? 1 : -1
            while occursin(peek(tokens), "<>")
                i += (popfirst!(tokens) == '>' ? 1 : -1)
            end
            push!(stack[end].args, :(pointer = clamp(pointer + $i, 1:memorysize)))
        elseif token == ','
            push!(stack[end].args, :(
                if !eof(input)
                    skipchars(c -> c == 0x0d && peek(input) == 0x0a, input)
                    memory[pointer] = read(input, UInt8)
                end
            ))
        elseif token == '.'
            push!(stack[end].args, :(print(output, Char(memory[pointer]))))
        elseif token == '['
            loop = Expr(:while, :(!iszero(memory[pointer])), Expr(:block))
            push!(stack[end].args, loop)
            push!(stack, loop.args[end])
        elseif token == ']'
            if length(stack) == 1
                throw(ArgumentError("Unmatched ']'."))
            end
            pop!(stack)
        end
    end
    if length(stack) != 1
        throw(ArgumentError("Unmatched '['."))
    end
    push!(stack[end].args, :(
        if returnstring
            String(take!(output))
        end
    ))
    Expr(
        :function,
        :((;
            memorysize::Integer = 30000,
            input::Union{AbstractString,IO} = stdin,
            output::Union{IO,Type{String}} = stdout
        )),
        only(stack)
    )
end

end # module

max_stars_repo_name: stellartux/Esolangs

max_issues_repo_name: stellartux/Esolangs

max_forks_repo_name: stellartux/Esolangs


text:
function find_memory_type(physical_device, type_flag, properties::Vk.MemoryPropertyFlag)
    mem_props = get_physical_device_memory_properties(physical_device)
    indices = findall(x -> properties in x.property_flags, mem_props.memory_types[1:mem_props.memory_type_count]) .- 1
    if isempty(indices)
        error("Could not find memory with properties $properties")
    else
        ind = findfirst(i -> type_flag & (1 << i) ≠ 0, indices)
        if isnothing(ind)
            error("Could not find memory with type $type_flag")
        else
            indices[ind], mem_props[indices[ind]]
        end
    end
end

"""
Memory domains:
- `MEMORY_DOMAIN_HOST` is host-visible memory, ideal for uploads to the GPU. It is preferably coherent and non-cached.
- `MEMORY_DOMAIN_HOST_CACHED` is host-visible, cached memory, ideal for readbacks from the GPU. It is preferably coherent.
- `MEMORY_DOMAIN_DEVICE` is device-local. It may be visible (integrated GPUs).
"""
@enum MemoryDomain::Int8 begin
    MEMORY_DOMAIN_HOST
    MEMORY_DOMAIN_HOST_CACHED
    MEMORY_DOMAIN_DEVICE
end

function score(domain::MemoryDomain, properties)
    @match domain begin
        &MEMORY_DOMAIN_HOST => 10 * MEMORY_PROPERTY_HOST_VISIBLE_BIT in properties + MEMORY_PROPERTY_HOST_COHERENT_BIT in properties - MEMORY_PROPERTY_HOST_CACHED_BIT in propetries
        &MEMORY_DOMAIN_HOST_CACHED => 10 * MEMORY_PROPERTY_HOST_VISIBLE_BIT | MEMORY_PROPERTY_HOST_CACHED_BIT in properties + MEMORY_PROPERTY_HOST_COHERENT_BIT in properties
        &MEMORY_DOMAIN_DEVICE => MEMORY_PROPERTY_DEVICE_LOCAL_BIT in properties
    end
end

function find_memory_type(f, physical_device, type_flag)
    mem_props = get_physical_device_memory_properties(physical_device)
    memtypes = mem_props.memory_types[1:mem_props.memory_type_count]
    candidate_indices = findall(i -> type_flag & (1 << i) ≠ 0, 0:length(memtypes) - 1)
    index, mem_prop = findmax(i -> f(memtypes[i + 1].property_flags), candidate_indices)
    index, mem_prop
end

find_memory_type(physical_device, type_flag, domain::MemoryDomain) = find_memory_type(Base.Fix1(score, domain), physical_device, type_flag)

mutable struct LinearAllocator
    memory::DeviceMemory
    size::Int
    properties::MemoryPropertyFlag
    last_offset::Int
end

device(allocator::LinearAllocator) = allocator.memory.device

function LinearAllocator(device::Device, size, properties)
    memory = DeviceMemory(device, size, properties)
    la = LinearAllocator(memory, size, properties, 0)
    map_memory(la)
    finalizer(unmap_memory, la)
end

function Vulkan.bind_buffer_memory(buffer, allocator::LinearAllocator, size, alignment)::ResultTypes.Result{Result,VulkanError}
    offset = alignment * cld(allocator.last_offset, alignment)
    ret = bind_buffer_memory(device(allocator), buffer, allocator.memory, offset)
    if !iserror(ret)
        allocator.last_offset = offset + size
    end
    ret
end

function Vulkan.map_memory(allocator::LinearAllocator, offset::Integer, size::Integer)
    if MEMORY_PROPERTY_HOST_COHERENT_BIT ∉ allocator.properties
        invalidate_mapped_memory_ranges(device(allocator), [MappedMemoryRange(allocator.memory, offset, size)])
    end
    map_memory(device(allocator), allocator.memory, offset, size)
end

function Vulkan.unmap_memory(allocator::LinearAllocator)
    if MEMORY_PROPERTY_HOST_COHERENT_BIT ∉ allocator.properties
        flush_mapped_memory_ranges(device(allocator), [MappedMemoryRange(allocator.memory, 0, allocator.size)])
    end
    unmap_memory(device(allocator), allocator.memory)
end

function reset!(allocator::LinearAllocator)
    allocator.last_offset = 0
end

buffer_size(data::DenseVector{T}) where {T} = sizeof(T) * length(data)
buffer_size(data) = sizeof(data)

# struct MemoryManager
#     allocators::Dictionary{MemoryDomain,}
# end

function Vulkan.DeviceMemory(device, memory_requirements::MemoryRequirements, properties)
    i, _ = find_memory_type(device.physical_device, memory_requirements.memory_type_bits, properties)
    DeviceMemory(device, memory_requirements.size, i)
end

"""
Upload data to the specified memory.

!!! warning
    The `memory` must be host coherent and host visible, otherwise the operation will fail.
"""
function upload_data(memory::DeviceMemory, data::DenseArray{T}; offset=0) where {T}
    memptr = unwrap(map_memory(memory.device, memory, offset, buffer_size(data)))
    GC.@preserve data unsafe_copyto!(Ptr{T}(memptr), pointer(data), length(data))
    unwrap(unmap_memory(memory.device, memory))
end

upload_data(resource::Allocated, data; offset=0) = upload_data(memory(resource), data; offset)

"""
Download data from the specified memory to an `Array`.

If `copy` if set to true, then all the data mapped from `memory` will be copied. If not, care should be taken to preserve the `memory` mapped and valid as long as the returned data is in use.
"""
function download_data(::Type{<:DenseArray{T}}, memory::DeviceMemory, dims; offset = 0, copy = true, unmap = true) where {T}
    size = sizeof(T) * prod(dims)
    memptr = unwrap(map_memory(memory.device, memory, offset, size))
    data = unsafe_wrap(Array, convert(Ptr{T}, memptr), dims; own = false)
    if unmap
        unwrap(unmap_memory(memory.device, memory))
    end
    if copy
        deepcopy(data)
    else
        data
    end
end

"""
Allocate a `DeviceMemory` object with the specified properties and bind it to the `buffer` using memory requirements from `get_buffer_memory_requirements`.
"""
function Vulkan.DeviceMemory(buffer::Buffer, properties::MemoryPropertyFlag) where {T}
    device = buffer.device
    memreqs = get_buffer_memory_requirements(device, buffer)
    memory = DeviceMemory(device, memreqs, properties)
    unwrap(bind_buffer_memory(device, buffer, memory, 0))
    memory
end

"""
Allocate a host visible and coherent `DeviceMemory` object, bind it to the `buffer` using memory requirements from `get_buffer_memory_requirements` and upload `data` to it.
"""
function Vulkan.DeviceMemory(buffer::Buffer, data::DenseArray{T}) where {T}
    memory = DeviceMemory(buffer, MEMORY_PROPERTY_HOST_VISIBLE_BIT | MEMORY_PROPERTY_HOST_COHERENT_BIT)
    upload_data(memory, data)
    memory
end

"""
Allocate a `DeviceMemory` object and bind it to the `image` using memory requirements from `get_image_memory_requirements`.
"""
function Vulkan.DeviceMemory(image::Image, properties::MemoryPropertyFlag)
    memreqs = get_image_memory_requirements(image.device, image)
    memory = DeviceMemory(image.device, memreqs, properties)
    unwrap(bind_image_memory(image.device, image, memory, 0))
    memory
end


text:
module SkipLists

#===========================
Includes
===========================#

include("core.jl")

include("node.jl")
include("node_nonconcurrent.jl")

include("list.jl")
include("list_nonconcurrent.jl")

#===========================
Exports
===========================#

export height
export SkipList, SkipListSet

end # module


text:
using HorizonSideRobots

function sherlock(r)
    kek=3::Int
    lol=1::Int
    while isborder(r,HorizonSide(0))
        for a in 1:lol
            move!(r,HorizonSide(kek % 4))
            if !isborder(r, HorizonSide(0))
            break
            end
        end
    lol += 2
    kek += 2
    end
end

sherlock(r)

max_stars_repo_name: Cheburkova/mirea-progs

max_issues_repo_name: Cheburkova/mirea-progs

max_forks_repo_name: Cheburkova/mirea-progs


text:
"""
# Description

A logical node that represents a group-by operation. It has two fields:

1. `source`: A node describing a relation that will be divided into groups.
2. `expressions`: A tuple of `Expression` objects that define the groups.
"""
struct GroupBy <: LogicalNode
    source::Any
    expressions::Tuple
end

function Base.print(io::IO, node::GroupBy)
    println(io, "GroupBy LogicalNode")
    println(io, type(node.source))
    for e in node.expressions
        println(io, e)
    end
end


text:
function plmdcasym(filename::AbstractString;
                decimation::Bool=false,
                fracmax::Float64 = 0.3,
                fracdec::Float64 = 0.1,
                blockdecimate::Bool=true, # decimate on a per-block base (all J[:,:,i,j] = 0)
                remove_dups::Bool = true,
                min_separation::Int = 1,
                max_gap_fraction::Real = 0.9,
                theta = :auto,
                lambdaJ::Real=0.01,
                lambdaH::Real=0.01,
                gaugecol::Int=-1,
                epsconv::Real=1.0e-5,
                maxit::Int=1000,
                verbose::Bool=true,
                method::Symbol=:LD_LBFGS)

    gaugecol >= 1 && println("Warning: gaugecol not implemented. Proceeding ...")

    W,Z,N,M,q = ReadFasta(filename,max_gap_fraction, theta, remove_dups)

    plmalg = PlmAlg(method,verbose, epsconv ,maxit, nothing)
    plmvar = PlmVar(N,M,q,q*q,gaugecol,lambdaJ,lambdaH,Z,W)

    if decimation == false
        Jmat, pseudolike = MinimizePLSym(plmalg,plmvar)
    else
        if blockdecimate
            decvar = DecVar{1}(fracdec, fracmax, blockdecimate, ones(Bool, binomial(N,2)))
        else
            decvar = DecVar{1}(fracdec, fracmax, blockdecimate, ones(Bool, binomial(N,2)*q*q+N*q))
        end
        Jmat, pseudolike = DecimateSym!(plmvar, plmalg, decvar)
    end
    score, Jtens, htens = ComputeScoreSym(Jmat, plmvar, min_separation)
    return output = PlmOut{3}(pseudolike, Jtens, htens, score)
end

function MinimizePLSym(alg::PlmAlg, var::PlmVar)

    N  = var.N
    q  = var.q
    q2 = var.q2

    Nc2 = binomial(N,2)
    LL  = Nc2 * q2  + N * q

    x0 = zeros(Float64, LL)

    opt = Opt(alg.method, length(x0))

    ftol_abs!(opt, alg.epsconv)
    maxeval!(opt, alg.maxit)
    min_objective!(opt, (x,g)->optimfunwrapper(x,g,var))
    elapstime = @elapsed  (minf, minx, ret) = optimize(opt, x0)
    @printf("pl = %.4f\t time = %.4f\t exit status = ", minf, elapstime)
    println(ret)

    return minx, minf
end

function ComputeScoreSym(Jvec::Array{Float64,1}, var::PlmVar, min_separation::Int)

    LL = length(Jvec)
    N=var.N
    q=var.q
    Nc2 = binomial(N,2)

    Jtens=reshape(Jvec[1:LL-N*q],q,q,Nc2)
    htens=fill(0.0,q,N)
    J = zeros(q,q,Nc2)

    htens=reshape(Jvec[LL-N*q + 1:end],q,N)

    for l=1:Nc2
        J[:,:,l] = Jtens[:,:,l] - repmat(mean(Jtens[:,:,l],1),q,1)-repmat(mean(Jtens[:,:,l],2),1,q) .+ mean(Jtens[:,:,l])
    end

    FN = zeros(Float64, N,N)
    l = 1
    for i=1:N-1
        for j=i+1:N
            FN[i,j] = vecnorm(J[:,:,l],2)
            FN[j,i] = FN[i,j]
            l+=1
        end
    end
    FN = GaussDCA.correct_APC(FN)
    score = GaussDCA.compute_ranking(FN,min_separation)
    return score, Jtens,htens
end

function PLsiteAndGradSym!(vecJ::Array{Float64,1}, grad::Array{Float64,1}, plmvar::PlmVar)

    LL = length(vecJ)
    q2 = plmvar.q2
    q = plmvar.q
    N = plmvar.N
    M = plmvar.M
    Z = plmvar.Z
    W = plmvar.W
    lambdaJ = plmvar.lambdaJ
    lambdaH = plmvar.lambdaH

    for i=1:LL-N*q
        grad[i] = 2.0 * vecJ[i] * lambdaJ
    end
    for i=(LL-N*q + 1):LL
        grad[i] = 4.0 * vecJ[i] * lambdaH
    end
    pseudolike = 0
    for a = 1:M
        pseudolike += ComputePatternPLSym!(grad, vecJ, Z[:,a], W[a], N, q, q2)
    end

    pseudolike += L2norm_sym(vecJ, plmvar)
    return pseudolike
end

function ComputePatternPLSym!(grad::Array{Float64,1}, vecJ::Array{Float64,1}, Z::Array{Int,1}, Wa::Float64, N::Int, q::Int, q2::Int)

    vecene = zeros(Float64,q)
    expvecenesunorm = zeros(Float64,q)
    pseudolike = 0
    offset = mygetindex(N-1, N, q, q, N, q, q2)
    @inbounds for site=1:N    # site < i
        fillvecenesym!(vecene, vecJ, Z, site, q,N)
        norm = sumexp(vecene)
        expvecenesunorm = exp.(vecene .- log(norm))
        pseudolike -= Wa * ( vecene[Z[site]] - log(norm) )
	for i = 1:(site-1)
            for s = 1:q
                grad[ mygetindex(i, site, Z[i], s, N, q, q2) ] += 0.5 * Wa * expvecenesunorm[s]
            end
            grad[ mygetindex(i, site , Z[i], Z[site],  N,q,q2)] -= 0.5 * Wa
        end
	for i = (site+1):N
            for s = 1:q
                grad[ mygetindex(site, i , s,  Z[i], N,q,q2) ] += 0.5 * Wa * expvecenesunorm[s]
            end
            grad[ mygetindex(site, i , Z[site], Z[i], N,q,q2)] -= 0.5* Wa
        end
        @simd for s = 1:q
            grad[ offset + s ] += Wa *  expvecenesunorm[s]
        end
	grad[ offset + Z[site] ] -= Wa
 	offset += q
    end
    return pseudolike
end

function fillvecenesym!(vecene::Array{Float64,1}, vecJ::Array{Float64,1}, Z::Array{Int64,1}, site::Int, q::Int ,N::Int)
    q2 = q*q
#    Z = sdata(sZ)

    @inbounds begin
        for l = 1:q
            offset::Int = 0
            scra::Float64 = 0.0

            for i=1:1:site-1
                scra += vecJ[ mygetindex(i, site, Z[i], l,  N, q, q2)]
            end
    	    for i = site+1:N
                scra += vecJ[ mygetindex(site, i, l, Z[i], N, q, q2)]
            end # End sum_i \neq site J
            scra *= 0.5
            offset = mygetindex(N-1, N, q, q, N, q, q2)  + ( site - 1) * q  # last J element + (site-1)*q
#            println(length(vecJ), " offset ", offset)
            scra += vecJ[offset + l] # sum H
            vecene[l] = scra
        end
    end
end

function L2norm_sym(vec::Array{Float64,1}, var::PlmVar)

    q = var.q
    N = var.N
    lambdaJ = var.lambdaJ
    lambdaH = var.lambdaH

    LL = length(vec)

    mysum1 = 0.0
    @inbounds @simd for i=1:(LL-N*q)
        mysum1 += vec[i] * vec[i]
    end
    mysum1 *= lambdaJ

    mysum2 = 0.0
    @inbounds @simd for i=(LL - N*q + 1):LL
        mysum2 += vec[i] * vec[i]
    end
    mysum2 *= 2lambdaH

    return mysum1+mysum2
end

function mygetindex( i::Int, j::Int, coli::Int, colj::Int, N::Int, q::Int, q2::Int)
    offset_i = ( (i-1) * N  - ( (i * ( i -1 ) ) >> 1 ) ) * q2 # (i-1) N q2 + i (i-1) q2 / 2
    offset_j = (j - i - 1 ) * q2
    return offset_i + offset_j + coli + q * (colj - 1)
end

nothing

max_stars_repo_name: YoshitakaMo/PlmDCA

max_issues_repo_name: YoshitakaMo/PlmDCA

max_forks_repo_name: YoshitakaMo/PlmDCA


text:
# function for computing and applying the mass matrix and its inverse

@doc """
  This function calculates the inverse mass matrix and returns it.
  Because we use SBP operators, the mass matrix is diagonal, so it is stored
  in a vector.  mesh.dofs is used to put the components of the inverse
  mass matrix in the same place as the corresponding values in eqn.res_vec

  Arguments:
    mesh: AbstractMesh
    sbp: SBP operator
    eqn: an implementation of AbstractSolutionData. Does not have to be fully initialized.

  Outputs:
    Minv: vector containing inverse mass matrix

"""->
# used by AbstractSolutionData Constructor
# mid level functions
function calcMassMatrixInverse(mesh::AbstractMesh{Tmsh},
                               sbp::AbstractOperator,
                               eqn::AbstractSolutionData{Tsol, Tres}) where {Tmsh,  Tsol, Tres}
# calculate the inverse mass matrix so it can be applied to the entire solution vector
# mass matrix is diagonal, stores in vector eqn.Minv

  Minv = zeros(Tmsh, mesh.numDof)

  for i=1:mesh.numEl
    for j=1:sbp.numnodes
      for k=1:mesh.numDofPerNode
        dofnum_k = mesh.dofs[k,j,i]
        # multiplication is faster than division, so do the divisions here
        # and then multiply solution vector times Minv
        Minv[dofnum_k] += (sbp.w[j]/mesh.jac[j,i])
      end
    end
  end

  for i=1:mesh.numDof
    Minv[i] = 1/Minv[i]
  end

  return Minv

end     # end of calcMassMatrixInverse function

@doc """
  This function calculate the mass matrix and returns it.
  Beause w are using SBP operators, the mass matrix is diagonal, so it is
  stored in a vector.

  Arguments:
    mesh: AbstractMesh
    sbp: SBP operator
    eqn: an implementation of AbstractSolutionData. Does not have to be fully initialized.

  Outputs:
    M: vector containing mass matrix

"""->
function calcMassMatrix(mesh::AbstractMesh{Tmsh},
                        sbp::AbstractOperator,
                        eqn::AbstractSolutionData{Tsol, Tres}) where {Tmsh,  Tsol, Tres}
# calculate the (diagonal) mass matrix as a vector
# return the vector M

  M = zeros(Tmsh, mesh.numDof)

  for i=1:mesh.numEl
    for j=1:sbp.numnodes
      for k=1:mesh.numDofPerNode
        dofnum_k = mesh.dofs[k,j,i]
        # multiplication is faster than division, so do the divions here
        # and then multiply solution vector times M
        M[dofnum_k] += (sbp.w[j]/mesh.jac[j,i])
      end
    end
  end

  return M

end     # end of calcMassMatrix function

@doc """

  This function calculates the mass matrix and returns it, in a 3D array format suitable for
  application to eqn.res.

  Arguments:
    mesh: AbstractMesh
    sbp: SBP operator
    eqn: an implementation of AbstractSolutionData. Does not have to be fully initialized.

  Outputs:
    M: vector containing mass matrix

"""->
# calcMassMatrixInverse3D:
#   calculates the inverse mass matrix, returning it as a 3D array suitable for application to eqn.res
function calcMassMatrixInverse3D(mesh::AbstractMesh{Tmsh},
                               sbp::AbstractOperator,
                               eqn::AbstractSolutionData{Tsol, Tres}) where {Tmsh,  Tsol, Tres}

  Minv3D = zeros(Tmsh, mesh.numDofPerNode, sbp.numnodes, mesh.numEl)

  for i=1:mesh.numEl
    for j=1:sbp.numnodes
      for k=1:mesh.numDofPerNode
        dofnum_k = mesh.dofs[k,j,i]
        # multiplication is faster than division, so do the divisions here
        # and then multiply solution vector times Minv
        Minv3D[k, j, i] += (sbp.w[j]/mesh.jac[j,i])
      end
    end
  end

  for i=1:mesh.numEl
    for j=1:sbp.numnodes
      for k=1:mesh.numDofPerNode
        Minv3D[k, j, i] = 1/Minv3D[k, j, i]
      end
    end
  end

  return Minv3D

end

@doc """
  This function multiplies eqn.res_vec (the residual in vector form  by eqn.Minv,
  the diagonal mass matrix.  This is a very fast function because all values
  are precomputed and stored linearly in memory.

  This is a mid level function, and does the correct thing regardless of the
  dimension of the equation.

  Aliasing restrictions: none
"""->
# mid level function (although it doesn't really need to Tdim)
function applyMassMatrixInverse(eqn::AbstractSolutionData{Tsol, Tres},
                                res_vec::AbstractVector{Tsol}) where {Tsol, Tres}
  # apply the inverse mass matrix stored eqn to res_vec

  ndof = length(res_vec)
  for i=1:ndof
    res_vec[i] *= eqn.Minv[i]
  end

  return nothing
end


text:
"""
    execute_sts_monthly_api_call(user_model::PSUserModel, stock_symbol::String;
        data_type::Symbol = :json, outputsize::Symbol = :compact, logger::Union{Nothing,AbstractLogger} = nothing) -> PSResult

Stuff will go here. Awesome stuff, the most beautiful stuff ever.
"""
function execute_sts_monthly_api_call(user_model::PSUserModel, stock_symbol::String;
    data_type::Symbol = :json, outputsize::Symbol = :compact, logger::Union{Nothing,AbstractLogger} = nothing)::PSResult

    try

        # get key's and id -
        api_key = user_model.alphavantage_api_key

        # use alpha_vantage_api to access data change function to monthly in this case
        # change this line for different time periods
        url = "$(alphavantage_api_url_string)?function=TIME_SERIES_MONTHLy&symbol=$(stock_symbol)&apikey=$(api_key)&datatype=$(string(data_type))&outputsize=$(string(outputsize))"
        api_call_raw_data = http_get_call_with_url(url) |> checkresult

        # call to logger if we have one -
        if (isnothing(logger) == false)
            log_api_call(logger, user_model, url)
        end

        # call a different method depending upon data type -
        if (data_type == :json)

            # process json
            data_series_key = "Monthly Time Series"
            return process_raw_json_api_data_sts(api_call_raw_data, data_series_key)
        elseif (data_type == :csv)

            # process csv results -
            return process_raw_csv_api_data(api_call_raw_data)
        else

            # formulate the error message -
            error_message = "$(data_type) is not supported. Supported values are {:json,:csv}"
            return PSResult{PSError}(PSError(error_message))
        end

    catch error
        return PSResult(error)
    end
end

"""
    execute_sts_adjusted_monthly_api_call(user_model::PSUserModel, stock_symbol::String;
        data_type::Symbol = :json, outputsize::Symbol = :compact, logger::Union{Nothing,AbstractLogger} = nothing) -> PSResult

Stuff will go here. Awesome stuff, the most beautiful stuff ever.
"""
function execute_sts_adjusted_monthly_api_call(user_model::PSUserModel, stock_symbol::String;
    data_type::Symbol = :json, outputsize::Symbol = :compact, logger::Union{Nothing,AbstractLogger} = nothing)::PSResult

    try

        # get key's and id -
        api_key = user_model.alphavantage_api_key

        # use alpha_vantage_api to access data change function to monthly in this case
        # change this line for different time periods
        url = "$(alphavantage_api_url_string)?function=TIME_SERIES_MONTHLY_ADJUSTED&symbol=$(stock_symbol)&apikey=$(api_key)&datatype=$(string(data_type))&outputsize=$(string(outputsize))"
        api_call_raw_data = http_get_call_with_url(url) |> checkresult

        # call to logger
        if (isnothing(logger) == false)
            log_api_call(logger, user_model, url)
        end

        if (data_type == :json)

            # process json
            # this string line should match the .json of the type we are calling
            data_series_key = "Monthly Adjusted Time Series"
            return process_raw_json_api_data_sts_adjusted(api_call_raw_data, data_series_key)

        elseif (data_type == :csv)
            return process_raw_csv_api_data(api_call_raw_data)
        else
            # formulate the error message -
            error_message = "$(data_type) is not supported. Supported values are {:json,:csv}"
            return PSResult{PSError}(PSError(error_message))
        end

    catch error
        return PSResult(error)
    end
end


text:
using ACE
using Printf, Test, LinearAlgebra, StaticArrays, BenchmarkTools
using ACE: evaluate, evaluate_d, evaluate_ed,
      Rn1pBasis, Ylm1pBasis,
      PositionState, Product1pBasis, getlabel, get_spec,
      State, DState, rand_vec3, rand_radial, rand_sphere, Scal1pBasis,
      valtype, gradtype, acquire_B!, acquire_dB!,
      discrete_jacobi, release!
using Random: shuffle
using ACEbase.Testing: dirfdtest, fdtest, print_tf, test_fio, println_slim

##

@info "Build a 1p basis from scratch"

maxdeg = 20
maxL = 10
r0 = 1.0
rcut = 3.0
maxorder = 3
Bsel = SimpleSparseBasis(maxorder, maxdeg)

trans = PolyTransform(1, r0)   # r -> x = 1/r^2
J = discrete_jacobi(maxdeg; pcut=2, xcut = rcut, pin = 0, xin = 0.0, trans=trans)
Rn = Rn1pBasis(J, trans)
Ylm = Ylm1pBasis(maxL)
Pk = Scal1pBasis(:u, nothing, :k, J)
A_nlm = Product1pBasis( (Rn, Ylm) )
ACE.init1pspec!(A_nlm, Bsel)
A_nlmk = Product1pBasis( (Rn, Ylm, Pk) )
ACE.init1pspec!(A_nlmk, Bsel)

nX = 10
Xs = [ State(rr = rand_vec3(Rn), u = randn() ) for _=1:nX ]
X = Xs[1]
cfg = ACEConfig(Xs)

# A = evaluate(A_nlm, Xs)

##

@info("Evaluate Pk, Ylm, Rn WITHOUT release")
@btime B = evaluate($Pk, $X)
@btime B = evaluate($Ylm, $X)
@btime B = evaluate($Rn, $X)

@info("Evaluate Pk, Ylm, Rn WITH release")
@btime begin B = evaluate($Pk, $X); release!(B);  end
@btime begin B = evaluate($Ylm, $X); release!(B); end
@btime begin B = evaluate($Rn, $X); release!(B);  end

@info("Gradient Pk, Ylm, Rn WITHOUT release")
@btime begin B, dB = evaluate_ed($Pk, $X);  release!(B); release!(dB); end
@btime begin B, dB = evaluate_ed($Ylm, $X); release!(B); release!(dB); end
@btime begin B, dB = evaluate_ed($Rn, $X);  release!(B); release!(dB); end

##

@btime begin B, dB = evaluate_ed($(Ylm.basis), $(X.rr));
             release!(B); release!(dB); end

S = ACE.SphericalHarmonics.cart2spher(X.rr)
SH = Ylm.basis
P, dP = ACE.SphericalHarmonics._evaluate_ed(SH.alp, S)
Y = ACE.acquire!(SH.B_pool, length(SH), valtype(SH, X.rr))
dY = ACE.acquire!(SH.dB_pool, length(SH), gradtype(SH, X.rr))
_maxL = ACE.SphericalHarmonics.maxL(SH)

@btime begin P, dP = ACE.SphericalHarmonics._evaluate_ed($(SH.alp), $S);
         release!(P); release!(dP); end

@btime ACE.SphericalHarmonics.cYlm_ed!($Y, $dY, $_maxL, $S, $P, $dP)

@btime ACE.SphericalHarmonics.evaluate_ed!($Y, $dY, $SH, $(X.rr))

ACE.SphericalHarmonics.evaluate_ed!(Y, dY, SH, X.rr)

##

function runn(N, f, args...)
   for n = 1:N
      Y, dY = f(args...)
      release!(Y)
      release!(dY)
   end
end

##

runn(10, evaluate_ed, Ylm.basis, X.rr)

@profview runn(100_000, evaluate_ed, Ylm.basis, X.rr)

@btime runn(10, $evaluate_ed, $(Ylm.basis), $(X.rr))


text:
# module for calculating standard deviation parameters
# robot rotates at same position
# including random noise and bias
# this is iterated for 100 times

module MotionTestRotBias
  using DataFrames, Statistics

  include(joinpath(split(@__FILE__, "src")[1], "src/model/world/world.jl"))
  include(joinpath(split(@__FILE__, "src")[1], "src/model/robot/differential_wheeled_robot/differential_wheeled_robot.jl"))

  function main(time_interval=0.1)
    w = World(-5.0, 5.0, -5.0, 5.0)

    rs = []
    init_pose = [0.0, 0.0, 0.0]
    rotate = Agent(0.0, 0.1)
    for i in 1:100
      r = DifferentialWheeledRobot(init_pose, 0.2, "black",
                                   rotate, time_interval,
                                   noise_per_meter=5,
                                   noise_std=pi/60,
                                   bias_rate_stds=[0.1, 0.1])
      push!(rs, r)
    end

    for t in 0:time_interval:40
      draw(w)
      annotate!(-3.5, 4.5, "t = $(time_interval)", "black")
      for r in rs
        draw!(r)
      end
    end

    poses = DataFrame(range=[sqrt(r.pose[1]^2 + r.pose[2]^2) for r in rs],
                      theta=[r.pose[3] for r in rs])
    theta_var = var(poses.theta, corrected=false)
    theta_mean = mean(poses.theta)
    println("Variance of range and theta: ", poses)
    println("Variance of theta: ", theta_var)
    println("Mean of theta: ", theta_mean)
    println("Standard deviation per deg of theta: ", sqrt(theta_var/theta_mean))
  end
end

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy


text:
using BenchmarkTools
using Printf
using Transducers

#### Types

struct Group
    records::Array{Set{Char}, 1}
end

function questionsWithYesAnswers(group::Group)
    reduce(union, group.records; init = Set{Char}())
end

function questionsEveryoneAnsweredYesTo(group::Group)
    reduce(intersect, group.records; init = Set{Char}(collect('a':'z')))
end

#### Record source (iterator)

struct RecordSource
    io::IOStream
end

function Base.iterate(iter::RecordSource, state = ())
    records = Array{Set{Char}, 1}()
    while !eof(iter.io)
        line = readline(iter.io)
        line == "" && !isempty(records) && break
        if line != ""
            push!(records, Set{Char}(collect(line)))
        end
    end
    !isempty(records) ? (Group(records), ()) : nothing
end

Base.IteratorSize(::Type{RecordSource}) = Base.SizeUnknown()

Base.IteratorEltype(::Type{RecordSource}) = Base.HasEltype()

Base.eltype(::Type{RecordSource}) = Group

#### Factory methods

function recordSourceFactory(filename)
    factory() = RecordSource(open(filename))
    return factory
end

####

exampleRecordSource = recordSourceFactory("src/advent/day6/example.txt")
inputRecordSource = recordSourceFactory("src/advent/day6/input.txt")

#### Part 1

function part1(recordSource)
    recordSource |> Map(group -> length(questionsWithYesAnswers(group))) |> sum
end

@assert part1(exampleRecordSource()) == 11

@printf("part1 = %d\n", part1(inputRecordSource()))

#collectedRecords = collect(inputRecordSource())
#@btime part1(collectedRecords);

#### Part 2

function part2(recordSource)
    recordSource |> Map(group -> length(questionsEveryoneAnsweredYesTo(group))) |> sum
end

@printf("part2 = %d\n", part2(inputRecordSource()))

#collectedRecords = collect(inputRecordSource())
#@btime part2(collectedRecords);


text:
using Random, Plots, StatsPlots, Distributions, DataFrames, Printf, SpecialFunctions
# Define student t distribution with loc and scale parameter
# https://juliastats.org/Distributions.jl/v0.21/univariate/#Distributions.LocationScale
# https://discourse.julialang.org/t/scale-and-location-to-tdist-in-distributions/74360
# p(x|\nu, \mu, \sigma^2) = \frac{\Gamma(\frac{\nu+1}{2})}{\Gamma(\frac{\nu}{2})\sqrt{\pi\nu\sigma^2}}\left[1 + \frac{(x-\mu)^2}{\nu\sigma^2}\right]^{-\frac{\nu+1}{2}}
"""
t-distribution probability density function from definition.
"""
function tdist_from_def(x, ν, μ, σ)
    return (gamma((ν+1)/2) / (gamma(ν/2) * sqrt(pi*ν*σ^2)) * (1 + (x-μ)^2/(ν*σ^2))^(-(ν+1)/2))
end

"""
t-distribution probability density function usinf Distributions.LocationScale
"""
function tdist_from_LocationScale(ν, μ, σ)
    return LocationScale(μ, σ, TDist(ν))
end

# tdist_from_def(5, 7, 10, 0.7)
# pdf(tdist_from_LocationScale(7, 10, 0.7), 5)

# Define calculation of posterior
"""
Conjugate Priors for Normal Dist parameters (μ|σ^2 ~ Normal, σ^2 ~ InverseGamma)

    Inputs
    ------
    data: data (0 or 1)
    α0: shape parameter of prior (Gamma dist)
    β0: inversed scale parameter of prior (Gamma dist)
    prob: credible interval probability (0 < prob < 1)

    Outputs
    -------
    results: posterior summary stats table (DataFrames)
    α_star: shape parameter of posterior (Gamma dist)
    β: inversed scale parameter (1/θ) of posterior (Gamma dist)

"""
function gaussian_stats(data, μ0, n0, ν0, λ0, prob)
    # posterior
    n = size(data)[1]
    mean_data = mean(data)
    ssd_data = (n-1) * var(data) # var is unbiased variance scaled by n-1
    n_star = n + n0
    μ_star = (n * mean_data + n0 * μ0) / n_star
    ν_star = n + ν0
    λ_star = ssd_data + (n * n0 / n_star) * (μ0 - mean_data)^2 + λ0
    τ_star = sqrt(λ_star / (ν_star * n_star))

    posterior_σ_sq = InverseGamma(ν_star/2, λ_star/2)
    marginal_posterior_μ = LocationScale(μ_star, τ_star, TDist(ν_star))

    sd_μ = std(marginal_posterior_μ)
    ci_μ_upper = quantile(marginal_posterior_μ, prob - (1-prob)/2)
    ci_μ_lower = quantile(marginal_posterior_μ, (1-prob)/2)
    ci_μ = [ci_μ_lower, ci_μ_upper]
    mean_σ_sq = mean(posterior_σ_sq)
    mode_σ_sq = mode(posterior_σ_sq)
    median_σ_sq = median(posterior_σ_sq)
    sd_σ_sq = std(posterior_σ_sq)
    ci_σ_sq_upper = quantile(posterior_σ_sq, prob - (1-prob)/2)
    ci_σ_sq_lower = quantile(posterior_σ_sq, (1-prob)/2)
    ci_σ_sq = [ci_σ_sq_lower, ci_σ_sq_upper]

    stats_μ = reshape([μ_star, μ_star, μ_star, sd_μ, ci_μ], 1, 5)
    stats_σ_sq = reshape([mean_σ_sq, median_σ_sq, mode_σ_sq, sd_σ_sq, ci_σ_sq], 1, 5)

    stats_string = ["Mean", "Median", "Mode", "SD", "CI"]
    results = DataFrame([stats_μ; stats_σ_sq], stats_string)
    return results, marginal_posterior_μ, posterior_σ_sq
end

# Generate moc data
μ = 1.0
σ = 2.0
n = 50
Random.seed!(99)
data = rand(Normal(μ, σ), n)
# histogram(data, bins=20)

# Set prior
μ0 = 0.0
n0 = 0.2
ν0 = 5.0
λ0 = 7.0
τ0 = sqrt(λ0 / (ν0 * n0))

prior_σ_sq = InverseGamma(ν0/2, λ0/2)
marginal_prior_μ = LocationScale(μ0, τ0, TDist(ν0))

# Calculate posterior statistics
prob = 0.95
results, marginal_posterior_μ, posterior_σ_sq = gaussian_stats(data, μ0, n0, ν0, λ0, prob)
print(results)

# Visualize Posterior
l = @layout [grid(1, 2)]

p1 = plot(
    marginal_prior_μ, xlims=(-6, 6), ylims=(0, 1.55), xlabel="μ", ylabel="Probability Density",
    label=@sprintf("prior: T(ν=%.1f, μ=%.1f, σ^2=%.1f)", params(params(marginal_prior_μ)[3])[1], params(marginal_prior_μ)[1], params(marginal_prior_μ)[2]), linestyle=:dash)
plot!(marginal_posterior_μ, label=@sprintf("posterior: T(ν=%.1f, μ=%.1f, σ^2=%.1f)", params(params(marginal_posterior_μ)[3])[1], params(marginal_posterior_μ)[1], params(marginal_posterior_μ)[2]), linestyle=:solid)

p2 = plot(
    prior_σ_sq, xlims=(0, 10), ylims=(0, 0.65), xlabel="σ^2", ylabel="Probability Density",
    label=@sprintf("prior: InverseGamma(α=%.1f, θ=%.1f)", params(prior_σ_sq)[1], params(prior_σ_sq)[2]), linestyle=:dash)
plot!(posterior_σ_sq, label=@sprintf("posterior: InverseGamma(α=%.1f, θ=%.1f)", params(posterior_σ_sq)[1], params(posterior_σ_sq)[2]), linestyle=:solid)

plot(p1, p2, laypout=l)
# μの推定ヘタクソじゃね?


text:
const CONFIG = MOIT.Config(Int)

const OPTIMIZER = CPLEXCP.Optimizer()
MOI.set(OPTIMIZER, MOI.Silent(), true)
const BRIDGED_OPTIMIZER = MOI.Bridges.full_bridge_optimizer(OPTIMIZER, Float64)

COIT.runtests(
    BRIDGED_OPTIMIZER,
    CONFIG
)


text:
const MD_TOGGLE = true

using MatrixDepot
using Test
using LinearAlgebra
using SparseArrays

import MatrixDepot: DataError

macro inc(a)
    :(@testset $a begin let n, p, i; include($a) end end)
end

include("clean.jl")
data_save = save_target(data_dir)
user_save = save_target(user_dir)

# that will download the index files if necessary and initialize internal data
MatrixDepot.init()
@testset "MatrixDepot generator tests" begin

    @inc("test_magic.jl")
    @inc("test_cauchy.jl")
    @inc("test_circul.jl")
    @inc("test_hadamard.jl")
    @inc("test_hilb.jl")
    @inc("test_dingdong.jl")
    @inc("test_frank.jl")
    @inc("test_invhilb.jl")
    @inc("test_forsythe.jl")
    @inc("test_grcar.jl")
    @inc("test_triw.jl")
    @inc("test_moler.jl")
    @inc("test_pascal.jl")
    @inc("test_kahan.jl")
    @inc("test_pei.jl")
    @inc("test_vand.jl")
    @inc("test_invol.jl")
    @inc("test_chebspec.jl")
    @inc("test_lotkin.jl")
    @inc("test_clement.jl")
    @inc("test_fiedler.jl")
    @inc("test_minij.jl")
    @inc("test_binomial.jl")
    @inc("test_tridiag.jl")
    @inc("test_lehmer.jl")
    @inc("test_parter.jl")
    @inc("test_chow.jl")
    @inc("test_randcorr.jl")
    @inc("test_poisson.jl")
    @inc("test_neumann.jl")
    @inc("test_rosser.jl")
    @inc("test_sampling.jl")
    @inc("test_wilkinson.jl")
    @inc("test_rando.jl")
    @inc("test_randsvd.jl")
    @inc("test_rohess.jl")
    @inc("test_kms.jl")
    @inc("test_wathen.jl")
    @inc("test_oscillate.jl")
    @inc("test_toeplitz.jl")
    @inc("test_hankel.jl")
    @inc("test_golub.jl")
    @inc("test_companion.jl")
    @inc("test_erdrey.jl")
    @inc("test_gilbert.jl")
    @inc("test_smallworld.jl")
end

@testset "MatrixDepot remote matrix tests" begin
    tests = [
            "regu",
            "include_generator",
            "download",
            "common",
            "number",
            "property",
            ]

    try
        @testset "$t" for t in tests
            tp = joinpath(@__DIR__(), "$(t).jl")
            println("running $(tp) ...")
            include(tp)
            println("finished $(tp)")
        end

        xdir = MatrixDepot.DATA_DIR
        xtmpdir = string(xdir, ".tmp")
        println("mv $xdir $xtmpdir")
        mv(xdir, xtmpdir, force=true)
        MatrixDepot.toggle_remote()
        MatrixDepot.init(ignoredb=true)
        MatrixDepot.update()
        @test mdopen(sp(1)) != nothing
        mv(xtmpdir, xdir, force=true)
    finally
        revert_target(user_save, user_dir)
        revert_target(data_save, data_dir)
    end
end

println("Success in all tests.")


text:
using DifferentialEquations
using GLMakie

#=
 # 1. simplest case

  ̇x = - λ ̇x

=#

testeq(u, p, t) = @. - p.k * ( p.a - u ) * u

u0  = [0.1]
tspan = (0.0, 60.0)
p = (; k = 1.0, a = 1.6)

prob = ODEProblem(testeq, u0, tspan, p)
sol = solve(prob) #, reltol = 1e-6, abstol = 1e-8)

begin
    f = Figure()

    ax = Axis(f[1,1], title = "Numerical solution and interpolation", xlabel=L"x(t)", ylabel=L"t")
    scatter!(sol.t, sol[1,:])
    lines!(sol.t, sol[1,:])

    sol_int = sol(LinRange(tspan..., 100))
    lines!(sol_int.t, sol_int[1,:])

    mkpath("plots")
    save("plots/00_first_order.png", f)
    f
end

max_stars_repo_name: SteffenPL/Julia-for-mathematicians

max_issues_repo_name: SteffenPL/Julia-for-mathematicians

max_forks_repo_name: SteffenPL/Julia-for-mathematicians


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

function build_harfbuzz(ARGS, name::String)

    icu = name == "HarfBuzz_ICU"

    version = v"2.8.1"

    # Collection of sources required to build Harfbuzz
    sources = [
        ArchiveSource("https://github.com/harfbuzz/harfbuzz/releases/download/$(version)/harfbuzz-$(version).tar.xz",
                      "4124f663ec4bf4e294d9cf230668370b4249a48ff34deaf0f06e8fc82d891300"),
    ]

    # Bash recipe for building across all platforms
    # Side note: this is a great use-case for
    # https://github.com/JuliaPackaging/BinaryBuilder.jl/issues/778
    script = "ICU=$(icu)\n" * raw"""
cd $WORKSPACE/srcdir/harfbuzz-*/
mkdir build && cd build
meson .. \
    --cross-file="${MESON_TARGET_TOOLCHAIN}" \
    --buildtype=release \
    -Dcairo=enabled \
    -Dfreetype=enabled \
    -Dglib=enabled \
    -Dgobject=enabled \
    -Dgraphite=enabled \
    -Dintrospection=disabled \
    -Ddocs=disabled \
    -Dtests=disabled \
    -Dicu=auto \
    -Dicu_builtin=false \
    -Dcoretext=enabled
ninja -j${nproc}
if [[ "${ICU}" == true ]]; then
    # Manually install only ICU-related files
    cp src/libharfbuzz-icu*${dlext}* ${libdir}/.
    cp meson-private/harfbuzz-icu.pc ${prefix}/lib/pkgconfig/.
    cp ../src/hb-icu.h ${includedir}/harfbuzz/.
else
    ninja install
fi
"""

    # These are the platforms we will build for by default, unless further
    # platforms are passed in on the command line
    platforms = filter!(p -> arch(p) != "armv6l", supported_platforms(; experimental=true))

    # The products that we will ensure are always built
    products = if icu
        [
            LibraryProduct("libharfbuzz-icu", :libharfbuzz_icu),
        ]
    else
        [
            LibraryProduct("libharfbuzz", :libharfbuzz),
            LibraryProduct("libharfbuzz-subset", :libharfbuzz_subset),
            LibraryProduct("libharfbuzz-gobject", :libharfbuzz_gobject),
        ]
    end

    # Dependencies that must be installed before this package can be built
    dependencies = [
        Dependency("Cairo_jll"),
        Dependency("Fontconfig_jll"),
        Dependency("FreeType2_jll"),
        Dependency("Glib_jll"; compat="2.68.1"),
        Dependency("Graphite2_jll"),
        Dependency("Libffi_jll"; compat="~3.2.2"),
        BuildDependency("Xorg_xorgproto_jll"),
    ]

    if icu
        append!(dependencies, [
            Dependency("HarfBuzz_jll"; compat="$(version)"),
            Dependency("ICU_jll"; compat="69.1.0"),
        ])
    end

    # Build the tarballs, and possibly a `build.jl` as well.
    build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version=v"5", julia_compat="1.6")
end

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
## Integrator Dispatches

# Can get rid of an allocation here with a function
# get_tmp_arr(integrator.cache) which gives a pointer to some
# cache array which can be modified.

@inline function ode_addsteps!{calcVal,calcVal2,calcVal3}(integrator,f=integrator.f,always_calc_begin::Type{Val{calcVal}} = Val{false},allow_calc_end::Type{Val{calcVal2}} = Val{true},force_calc_end::Type{Val{calcVal3}} = Val{false})
  if !(typeof(integrator.cache) <: CompositeCache)
    ode_addsteps!(integrator.k,integrator.tprev,integrator.uprev,integrator.u,integrator.dt,f,integrator.cache,always_calc_begin,allow_calc_end,force_calc_end)
  else
    ode_addsteps!(integrator.k,integrator.tprev,integrator.uprev,integrator.u,integrator.dt,f,integrator.cache.caches[integrator.cache.current],always_calc_begin,allow_calc_end,force_calc_end)
  end
end

@inline function ode_interpolant(Θ,integrator::DEIntegrator,idxs,deriv)
  ode_addsteps!(integrator)
  if !(typeof(integrator.cache) <: CompositeCache)
    val = ode_interpolant(Θ,integrator.dt,integrator.uprev,integrator.u,integrator.k,integrator.cache,idxs,deriv)
  else
    val = ode_interpolant(Θ,integrator.dt,integrator.uprev,integrator.u,integrator.k,integrator.cache.caches[integrator.cache.current],idxs,deriv)
  end
  val
end

@inline function ode_interpolant!(val,Θ,integrator::DEIntegrator,idxs,deriv)
  ode_addsteps!(integrator)
  if !(typeof(integrator.cache) <: CompositeCache)
    ode_interpolant!(val,Θ,integrator.dt,integrator.uprev,integrator.u,integrator.k,integrator.cache,idxs,deriv)
  else
    ode_interpolant!(val,Θ,integrator.dt,integrator.uprev,integrator.u,integrator.k,integrator.cache.caches[integrator.cache.current],idxs,deriv)
  end
end

@inline function current_interpolant(t::Number,integrator::DEIntegrator,idxs,deriv)
  Θ = (t-integrator.tprev)/integrator.dt
  ode_interpolant(Θ,integrator,idxs,deriv)
end

@inline function current_interpolant(t,integrator::DEIntegrator,idxs,deriv)
  Θ = (t.-integrator.tprev)./integrator.dt
  [ode_interpolant(ϕ,integrator,idxs,deriv) for ϕ in Θ]
end

@inline function current_interpolant!(val,t::Number,integrator::DEIntegrator,idxs,deriv)
  Θ = (t-integrator.tprev)/integrator.dt
  ode_interpolant!(val,Θ,integrator,idxs,deriv)
end

@inline function current_interpolant!(val,t,integrator::DEIntegrator,idxs,deriv)
  Θ = (t.-integrator.tprev)./integrator.dt
  [ode_interpolant!(val,ϕ,integrator,idxs,deriv) for ϕ in Θ]
end

@inline function current_extrapolant(t::Number,integrator::DEIntegrator,idxs=size(integrator.uprev),deriv=Val{0})
  Θ = (t-integrator.tprev)/(integrator.t-integrator.tprev)
  ode_extrapolant(Θ,integrator,idxs,deriv)
end

@inline function current_extrapolant!(val,t::Number,integrator::DEIntegrator,idxs=size(integrator.uprev),deriv=Val{0})
  Θ = (t-integrator.tprev)/(integrator.t-integrator.tprev)
  ode_extrapolant!(val,Θ,integrator,idxs,deriv)
end

@inline function current_extrapolant(t::AbstractArray,integrator::DEIntegrator,idxs=size(integrator.uprev),deriv=Val{0})
  Θ = (t.-integrator.tprev)./(integrator.t-integrator.tprev)
  [ode_extrapolant(ϕ,integrator,idxs,deriv) for ϕ in Θ]
end

@inline function current_extrapolant!(val,t,integrator::DEIntegrator,idxs=size(integrator.uprev),deriv=Val{0})
  Θ = (t.-integrator.tprev)./(integrator.t-integrator.tprev)
  [ode_extrapolant!(val,ϕ,integrator,idxs,deriv) for ϕ in Θ]
end

@inline function ode_extrapolant!(val,Θ,integrator::DEIntegrator,idxs,deriv)
  ode_addsteps!(integrator)
  if !(typeof(integrator.cache) <: CompositeCache)
    ode_interpolant!(val,Θ,integrator.t-integrator.tprev,integrator.uprev2,integrator.uprev,integrator.k,integrator.cache,idxs,deriv)
  else
    ode_interpolant!(val,Θ,integrator.t-integrator.tprev,integrator.uprev2,integrator.uprev,integrator.k,integrator.cache.caches[integrator.cache.current],idxs,deriv)
  end
end

@inline function ode_extrapolant(Θ,integrator::DEIntegrator,idxs,deriv)
  ode_addsteps!(integrator)
  if !(typeof(integrator.cache) <: CompositeCache)
    ode_interpolant(Θ,integrator.t-integrator.tprev,integrator.uprev2,integrator.uprev,integrator.k,integrator.cache,idxs,deriv)
  else
    ode_interpolant(Θ,integrator.t-integrator.tprev,integrator.uprev2,integrator.uprev,integrator.k,integrator.cache.caches[integrator.cache.current],idxs,deriv)
  end
end

##

"""
ode_interpolation(tvals,ts,timeseries,ks)

Get the value at tvals where the solution is known at the
times ts (sorted), with values timeseries and derivatives ks
"""
@inline function ode_interpolation(tvals,id,idxs,deriv)
  @unpack ts,timeseries,ks,f,notsaveat_idxs,cache = id
  tdir = sign(ts[end]-ts[1])
  idx = sortperm(tvals,rev=tdir<0)
  i = 2 # Start the search thinking it's between ts[1] and ts[2]
  tvals[idx[end]] > ts[end] && error("Solution interpolation cannot extrapolate past the final timepoint. Either solve on a longer timespan or use the local extrapolation from the integrator interface.")
  if idxs == nothing
    if (eltype(timeseries) <: AbstractArray) && !(eltype(timeseries) <: Array)
      vals = Vector{Vector{eltype(first(timeseries))}}(length(tvals))
    else
      vals = Vector{eltype(timeseries)}(length(tvals))
    end
  elseif typeof(idxs) <: Number
    vals = Vector{eltype(first(timeseries))}(length(tvals))
  else
    vals = Vector{Vector{eltype(first(timeseries))}}(length(tvals))
  end
  @inbounds for j in idx
    t = tvals[j]
    i = searchsortedfirst(@view(ts[@view(notsaveat_idxs[i:end])]),t,rev=tdir<0)+i-1 # It's in the interval ts[i-1] to ts[i]
    if ts[notsaveat_idxs[i]] == t
      if idxs == nothing
        vals[j] = timeseries[notsaveat_idxs[i]]
      else
        vals[j] = timeseries[notsaveat_idxs[i]][idxs]
      end
    elseif ts[notsaveat_idxs[i-1]] == t # Can happen if it's the first value!
      if idxs == nothing
        vals[j] = timeseries[notsaveat_idxs[i-1]]
      else
        vals[j] = timeseries[notsaveat_idxs[i-1]][idxs]
      end
    else
      dt = ts[notsaveat_idxs[i]] - ts[notsaveat_idxs[i-1]]
      Θ = (t-ts[notsaveat_idxs[i-1]])/dt
      if idxs == nothing && eltype(timeseries) <: AbstractArray
        idxs_internal = size(timeseries[notsaveat_idxs[i-1]])
      else
        idxs_internal = idxs
      end
      if typeof(cache) <: CompositeCache
        ode_addsteps!(ks[i],ts[notsaveat_idxs[i-1]],timeseries[notsaveat_idxs[i-1]],timeseries[notsaveat_idxs[i]],dt,f,cache.caches[id.alg_choice[notsaveat_idxs[i-1]]]) # update the kcurrent
        vals[j] = ode_interpolant(Θ,dt,timeseries[notsaveat_idxs[i-1]],timeseries[notsaveat_idxs[i]],ks[i],cache.caches[id.alg_choice[notsaveat_idxs[i-1]]],idxs_internal,deriv)
      elseif typeof(cache) <: (DiscreteCache) || typeof(cache) <: DiscreteConstantCache
        vals[j] = ode_interpolant(Θ,dt,timeseries[notsaveat_idxs[i-1]],timeseries[notsaveat_idxs[i]],0,cache,idxs_internal,deriv)
      else
        ode_addsteps!(ks[i],ts[notsaveat_idxs[i-1]],timeseries[notsaveat_idxs[i-1]],timeseries[notsaveat_idxs[i]],dt,f,cache) # update the kcurrent
        vals[j] = ode_interpolant(Θ,dt,timeseries[notsaveat_idxs[i-1]],timeseries[notsaveat_idxs[i]],ks[i],cache,idxs_internal,deriv)
      end
    end
  end
  vals
end

"""
ode_interpolation(tvals,ts,timeseries,ks)

Get the value at tvals where the solution is known at the
times ts (sorted), with values timeseries and derivatives ks
"""
@inline function ode_interpolation!(vals,tvals,id,idxs,deriv)
  @unpack ts,timeseries,ks,f,notsaveat_idxs,cache = id
  tdir = sign(ts[end]-ts[1])
  idx = sortperm(tvals,rev=tdir<0)
  i = 2 # Start the search thinking it's between ts[1] and ts[2]
  tvals[idx[end]] > ts[end] && error("Solution interpolation cannot extrapolate past the final timepoint. Either solve on a longer timespan or use the local extrapolation from the integrator interface.")
  @inbounds for j in idx
    t = tvals[j]
    i = searchsortedfirst(@view(ts[@view(notsaveat_idxs[i:end])]),t,rev=tdir<0)+i-1 # It's in the interval ts[i-1] to ts[i]
    if ts[notsaveat_idxs[i]] == t
      if idxs == nothing
        vals[j] = timeseries[notsaveat_idxs[i]]
      else
        vals[j] = timeseries[notsaveat_idxs[i]][idxs]
      end
    elseif ts[notsaveat_idxs[i-1]] == t # Can happen if it's the first value!
      if idxs == nothing
        vals[j] = timeseries[notsaveat_idxs[i-1]]
      else
        vals[j] = timeseries[notsaveat_idxs[i-1]][idxs]
      end
    else
      dt = ts[notsaveat_idxs[i]] - ts[notsaveat_idxs[i-1]]
      Θ = (t-ts[notsaveat_idxs[i-1]])/dt
      if idxs == nothing && eltype(vals) <: AbstractArray
        idxs_internal = eachindex(vals[j])
      else
        idxs_internal = idxs
      end
      if typeof(cache) <: CompositeCache
        ode_addsteps!(ks[i],ts[notsaveat_idxs[i-1]],timeseries[notsaveat_idxs[i-1]],timeseries[notsaveat_idxs[i]],dt,f,cache.caches[id.alg_choice[notsaveat_idxs[i-1]]]) # update the kcurrent
        if eltype(timeseries) <: AbstractArray
          ode_interpolant!(vals[j],Θ,dt,timeseries[notsaveat_idxs[i-1]],timeseries[notsaveat_idxs[i]],ks[i],cache.caches[id.alg_choice[notsaveat_idxs[i-1]]],idxs_internal,deriv)
        else
          vals[j] = ode_interpolant(Θ,dt,timeseries[notsaveat_idxs[i-1]],timeseries[notsaveat_idxs[i]],ks[i],cache.caches[id.alg_choice[notsaveat_idxs[i-1]]],idxs_internal,deriv)
        end
      elseif typeof(cache) <: (DiscreteCache) || typeof(cache) <: DiscreteConstantCache
        if eltype(timeseries) <: AbstractArray
          ode_interpolant!(vals[j],Θ,dt,timeseries[notsaveat_idxs[i-1]],timeseries[notsaveat_idxs[i]],0,cache,idxs_internal,deriv)
        else
          vals[j] = ode_interpolant(Θ,dt,timeseries[notsaveat_idxs[i-1]],timeseries[notsaveat_idxs[i]],0,cache,idxs_internal,deriv)
        end
      else
        ode_addsteps!(ks[i],ts[notsaveat_idxs[i-1]],timeseries[notsaveat_idxs[i-1]],timeseries[notsaveat_idxs[i]],dt,f,cache) # update the kcurrent
        if eltype(timeseries) <: AbstractArray
          ode_interpolant!(vals[j],Θ,dt,timeseries[notsaveat_idxs[i-1]],timeseries[notsaveat_idxs[i]],ks[i],cache,idxs_internal,deriv)
        else
          vals[j] = ode_interpolant(Θ,dt,timeseries[notsaveat_idxs[i-1]],timeseries[notsaveat_idxs[i]],ks[i],cache,idxs_internal,deriv)
        end
      end
    end
  end
end

"""
ode_interpolation(tval::Number,ts,timeseries,ks)

Get the value at tval where the solution is known at the
times ts (sorted), with values timeseries and derivatives ks
"""
@inline function ode_interpolation(tval::Number,id,idxs,deriv)
  @unpack ts,timeseries,ks,f,notsaveat_idxs,cache = id
  tval > ts[end] && error("Solution interpolation cannot extrapolate past the final timepoint. Either solve on a longer timespan or use the local extrapolation from the integrator interface.")
  tdir = sign(ts[end]-ts[1])
  @inbounds i = searchsortedfirst(@view(ts[notsaveat_idxs]),tval,rev=tdir<0) # It's in the interval ts[i-1] to ts[i]
  @inbounds if ts[notsaveat_idxs[i]] == tval
    if idxs == nothing
      val = timeseries[notsaveat_idxs[i]]
    else
      val = timeseries[notsaveat_idxs[i]][idxs]
    end
  elseif ts[notsaveat_idxs[i-1]] == tval # Can happen if it's the first value!
    if idxs == nothing
      val = timeseries[notsaveat_idxs[i-1]]
    else
      val = timeseries[notsaveat_idxs[i-1]][idxs]
    end
  else
    dt = ts[notsaveat_idxs[i]] - ts[notsaveat_idxs[i-1]]
    Θ = (tval-ts[notsaveat_idxs[i-1]])/dt
    if idxs == nothing && eltype(timeseries) <: AbstractArray
      idxs_internal = size(timeseries[notsaveat_idxs[i-1]])
    else
      idxs_internal = idxs
    end
    if typeof(cache) <: CompositeCache
      ode_addsteps!(ks[i],ts[notsaveat_idxs[i-1]],timeseries[notsaveat_idxs[i-1]],timeseries[notsaveat_idxs[i]],dt,f,cache.caches[id.alg_choice[notsaveat_idxs[i-1]]]) # update the kcurrent
      val = ode_interpolant(Θ,dt,timeseries[notsaveat_idxs[i-1]],timeseries[notsaveat_idxs[i]],ks[i],cache.caches[id.alg_choice[notsaveat_idxs[i-1]]],idxs_internal,deriv)
    elseif typeof(cache) <: (DiscreteCache) || typeof(cache) <: DiscreteConstantCache
      val = ode_interpolant(Θ,dt,timeseries[notsaveat_idxs[i-1]],timeseries[notsaveat_idxs[i]],0,cache,idxs_internal,deriv)
    else
      ode_addsteps!(ks[i],ts[notsaveat_idxs[i-1]],timeseries[notsaveat_idxs[i-1]],timeseries[notsaveat_idxs[i]],dt,f,cache) # update the kcurrent
      val = ode_interpolant(Θ,dt,timeseries[notsaveat_idxs[i-1]],timeseries[notsaveat_idxs[i]],ks[i],cache,idxs_internal,deriv)
    end
  end
  val
end

"""
ode_interpolation!(out,tval::Number,ts,timeseries,ks)

Get the value at tval where the solution is known at the
times ts (sorted), with values timeseries and derivatives ks
"""
@inline function ode_interpolation!(out,tval::Number,id,idxs,deriv)
  @unpack ts,timeseries,ks,f,notsaveat_idxs,cache = id
  tval > ts[end] && error("Solution interpolation cannot extrapolate past the final timepoint. Either solve on a longer timespan or use the local extrapolation from the integrator interface.")
  tdir = sign(ts[end]-ts[1])
  @inbounds i = searchsortedfirst(@view(ts[notsaveat_idxs]),tval,rev=tdir<0) # It's in the interval ts[i-1] to ts[i]
  @inbounds if ts[notsaveat_idxs[i]] == tval
    if idxs == nothing
      copy!(out,timeseries[notsaveat_idxs[i]])
    else
      copy!(out,timeseries[notsaveat_idxs[i]][idxs])
    end
  elseif ts[notsaveat_idxs[i-1]] == tval # Can happen if it's the first value!
    if idxs == nothing
      copy!(out,timeseries[notsaveat_idxs[i-1]])
    else
      copy!(out,timeseries[notsaveat_idxs[i-1]][idxs])
    end
  else
    dt = ts[notsaveat_idxs[i]] - ts[notsaveat_idxs[i-1]]
    Θ = (tval-ts[notsaveat_idxs[i-1]])/dt
    if idxs == nothing
      idxs_internal = eachindex(out)
    else
      idxs_internal = idxs
    end
    if typeof(cache) <: CompositeCache
      ode_addsteps!(ks[i],ts[notsaveat_idxs[i-1]],timeseries[notsaveat_idxs[i-1]],timeseries[notsaveat_idxs[i]],dt,f,cache.caches[id.alg_choice[notsaveat_idxs[i-1]]]) # update the kcurrent
      ode_interpolant!(out,Θ,dt,timeseries[notsaveat_idxs[i-1]],timeseries[notsaveat_idxs[i]],ks[i],cache.caches[id.alg_choice[notsaveat_idxs[i-1]]],idxs_internal,deriv)
    elseif typeof(cache) <: (DiscreteCache) || typeof(cache) <: DiscreteConstantCache
      ode_interpolant!(out,Θ,dt,timeseries[notsaveat_idxs[i-1]],timeseries[notsaveat_idxs[i]],0,cache,idxs_internal,deriv)
    else
      ode_addsteps!(ks[i],ts[notsaveat_idxs[i-1]],timeseries[notsaveat_idxs[i-1]],timeseries[notsaveat_idxs[i]],dt,f,cache) # update the kcurrent
      ode_interpolant!(out,Θ,dt,timeseries[notsaveat_idxs[i-1]],timeseries[notsaveat_idxs[i]],ks[i],cache,idxs_internal,deriv)
    end
  end
end

"""
By default, simpledense
"""
function ode_addsteps!{calcVal,calcVal2,calcVal3}(k,t,uprev,u,dt,f,cache,always_calc_begin::Type{Val{calcVal}} = Val{false},allow_calc_end::Type{Val{calcVal2}} = Val{true},force_calc_end::Type{Val{calcVal3}} = Val{false})
  if length(k)<2 || calcVal
    if !(typeof(uprev)<:AbstractArray)
      copyat_or_push!(k,1,f(t,uprev))
      copyat_or_push!(k,2,f(t+dt,u))
    else
      rtmp = similar(cache.fsalfirst)
      f(t,uprev,rtmp)
      copyat_or_push!(k,1,rtmp)
      f(t+dt,u,rtmp)
      copyat_or_push!(k,2,rtmp,Val{false})
    end
  end
  nothing
end

@inline function ode_interpolant(Θ,dt,y₀,y₁,k,cache::OrdinaryDiffEqMutableCache,idxs,T::Type{Val{0}})
  if typeof(idxs) <: Tuple
    out = similar(y₀,idxs)
    idxs_internal=eachindex(y₀)
  else
    out = similar(y₀,length(idxs))
    idxs_internal=idxs
  end
  ode_interpolant!(out,Θ,dt,y₀,y₁,k,cache,idxs_internal,T)
  if typeof(idxs) <: Number
    return first(out)
  else
    return out
  end
end

"""
Hairer Norsett Wanner Solving Ordinary Differential Euations I - Nonstiff Problems Page 190

Herimte Interpolation, chosen if no other dispatch for ode_interpolant
"""
@inline function ode_interpolant(Θ,dt,y₀,y₁,k,cache,idxs,T::Type{Val{0}}) # Default interpolant is Hermite
  if typeof(y₀) <: AbstractArray
    if typeof(idxs) <: Tuple
      out = similar(y₀,idxs)
    else
      out = similar(y₀,length(idxs))
    end
    @inbounds for (j,i) in enumerate(idxs)
      out[j] = (1-Θ)*y₀[i]+Θ*y₁[i]+Θ*(Θ-1)*((1-2Θ)*(y₁[i]-y₀[i])+(Θ-1)*dt*k[1][i] + Θ*dt*k[2][i])
    end
  else
    out = (1-Θ)*y₀+Θ*y₁+Θ*(Θ-1)*((1-2Θ)*(y₁-y₀)+(Θ-1)*dt*k[1] + Θ*dt*k[2])
  end
  out
end

"""
Hairer Norsett Wanner Solving Ordinary Differential Euations I - Nonstiff Problems Page 190

Herimte Interpolation, chosen if no other dispatch for ode_interpolant
"""
@inline function ode_interpolant!(out,Θ,dt,y₀,y₁,k,cache,idxs,T::Type{Val{0}}) # Default interpolant is Hermite
  @inbounds for (j,i) in enumerate(idxs)
    out[j] = (1-Θ)*y₀[i]+Θ*y₁[i]+Θ*(Θ-1)*((1-2Θ)*(y₁[i]-y₀[i])+(Θ-1)*dt*k[1][i] + Θ*dt*k[2][i])
  end
end


text:
# TODO: Support `import` for APIs marked as overloadbale
"""
    PublicAPI.@strict using Module: name₁, name₂, …, nameₙ
    PublicAPI.@strict import Module

Enable strict import; i.e., fail on using non-public API.

The simple form `import Module` create a dummy object named `Module` that acts
like the original `Module` but forbids access to the internal names.

# Extended help

Limitation: Currently, `PublicAPI.@strict` with the simple form `import Module`
creates a dummy local module and a global constant named `Module` is bind to it.
Thus, unlike `import Module`, the expression `PublicAPI.@strict import Module`
cannot be evaluated more than once inside a module.  This is an implementation
detail that may be fixed in the future if we find a better implementation.
"""
macro strict(import_statement::Expr)
    @gensym source_module
    if is_explicit_using(import_statement)  # using M: a.b.c
        import_expr = import_statement
        qualified_module_path = import_statement.args[1].args[1]  # `M` in `using M: ...`

        imported_names = map(import_statement.args[1].args[2:end]) do ex
            @assert Meta.isexpr(ex, :., 1)   # what else?
            ex.args[1]::Symbol
        end
    elseif is_simple_import(import_statement)
        fullpath = import_statement.args[1].args
        shim = gensym("PublicAPI_wrapper_" * join(fullpath, "."))
        lastname = fullpath[end]::Symbol
        import_expr = quote
            baremodule $shim
            $(Expr(:import, import_statement.args[1]))
            for api in $PublicAPI.of($lastname)
                name = $nameof(api)
                $Base.@eval $shim const $(Expr(:$, :name)) =
                    $lastname.$(Expr(:$, :name))
            end
            for name in $names($lastname)
                $Base.@eval $shim export $(Expr(:$, :name))
            end
            end
            const $lastname = $shim
        end
        @assert import_expr.head === :block
        import_expr = Expr(:toplevel, import_expr.args...)
        qualified_module_path = import_statement.args[1]  # `M` in `using M``
        imported_names = Symbol[]  # no need to check
    else
        error("Unsupported syntax: ", import_statement)
    end

    if VERSION ≥ v"1.6"
        import_source_module =
            Expr(:import, Expr(:as, qualified_module_path, source_module))
    else
        @gensym scratch
        import_source_module = quote
            module $scratch
            $(Expr(:import, qualified_module_path))
            const $source_module = $(qualified_module_path.args[end])
            end
            using .$scratch: $source_module
        end
        @assert import_source_module.head === :block
        import_source_module = Expr(:toplevel, import_source_module.args...)
    end

    quote
        $import_expr # e.g., using M: a, b, c
        $import_source_module # import M as $source_module
        $strict_check($source_module, $(QuoteNode(imported_names)))
        nothing
    end |> esc
end

"""
    is_simple_import(ex::Expr) :: Bool

Check if `ex` is of the form `import Module`.
"""
is_simple_import(ex::Expr) =
    ex.head === :import &&
    length(ex.args) == 1 &&
    Meta.isexpr(ex.args[1], :., 1) &&
    all(x -> x isa Symbol, ex.args[1].args)

"""
    is_explicit_using(ex::Expr) :: Bool

Check if `ex` is of the form `using Module: name₁, name₂, …, nameₙ`.
"""
is_explicit_using(ex::Expr) =
    ex.head === :using &&
    length(ex.args) == 1 &&
    Meta.isexpr(ex.args[1], :(:)) &&
    all((x -> Meta.isexpr(x, :., 1)), ex.args[1].args)

function strict_check(source_module::Module, imported_names::Vector{Symbol})
    internals = setdiff(imported_names, nameof.(PublicAPI.of(source_module)))
    isempty(internals) && return
    error("Non-public API imported from $source_module: ", join(internals, ", "))
end


text:
module CPortGraphDynam

using Catlab.WiringDiagrams
using Catlab.WiringDiagrams.CPortGraphs
using  Catlab.CategoricalAlgebra
using  Catlab.CategoricalAlgebra.FinSets
using Catlab.Graphs
import Catlab.Graphs: Graph
using Catlab.Graphics.GraphvizGraphs: to_graphviz

import Catlab.WiringDiagrams: oapply

import ..UWDDynam: nstates, nports, eval_dynamics, euler_approx
import ..DWDDynam: AbstractMachine, ContinuousMachine, DiscreteMachine,
ninputs, noutputs
import Catlab.CategoricalAlgebra.CSets: migrate!

using Base.Iterators
import Base: show, eltype

ContinuousMachine{T}(nports::Int, nstates::Int, d::Function, r::Function) where T =
    ContinuousMachine{T}(nports, nstates, nports, d, r)
DiscreteMachine{T}(nports::Int, nstates::Int, d::Function, r::Function) where T =
    DiscreteMachine{T}(nports, nstates, nports, d, r)

migrate!(g::Graph, pg::OpenCPortGraph) = migrate!(g, migrate!(CPortGraph(), pg))

draw(g::Graph) = to_graphviz(g, prog="neato", edge_labels=true, node_labels=true)
draw(pg::OpenCPortGraph) = draw(migrate!(Graph(), pg))

concat(xs::Vector) = (collect ∘ Iterators.flatten)(xs)

nports(d::OpenCPortGraph, b::Int) = incident(d, b, :box) |> length
nports(d::OpenCPortGraph, b) = map(length, incident(d, b, :box))
nports(d::OpenCPortGraph, b::Colon) = map(length, incident(d, :, :box))

colimitmap!(f::Function, output, C::Colimit, input) = begin
    for (i,x) in enumerate(input)
        y = f(i, x)
        I = legs(C)[i](1:length(y))
        # length(I) == length(y) || error("colimitmap! attempting to fill $(length(I)) slots with $(length(y)) values")
        output[I] .= y
    end
    return output
end

@inline fillreadouts!(y, d, xs, Ports, statefun) = colimitmap!(y, Ports, xs) do i,x
    return x.readout(statefun(i))
end

@inline fillstates!(y, d, xs, States, statefun, inputfun, p, t) = colimitmap!(y, States, xs) do i, x
    return x.dynamics(statefun(i), inputfun(i), p, t)
end

fillreadins!(readins, d, readouts) = begin
    for b in parts(d, :Box)
        ports = incident(d, b, :box)
        for p in ports
            ws = incident(d, p, :tgt)
            qs = d[ws, :src]
            readins[p] += sum(readouts[qs])
        end
    end
    return readins
end

"""    oapply(d::OpenCPortGraph, ms::Vector)

Implements the operad algebras for directed composition of dynamical systems given a
composition pattern (implemented by an open circular port graph `d`)
and primitive systems (implemented by a collection of
machines `ms`).

Each box of the composition pattern `d` is filled by a machine with the
appropriate type signature. Returns the composite machine.
"""
function oapply(d::OpenCPortGraph, xs::Vector{Machine}) where {T, Machine<:AbstractMachine{T}}
    x -> FinSet(x.nstates)
    S = coproduct((FinSet∘nstates).(xs))
    Inputs = coproduct((FinSet∘ninputs).(xs))
    Ports = coproduct([FinSet.(nports(d, b)) for b in parts(d, :Box)])
    state(u::Vector, b::Int) = u[legs(S)[b](1:xs[b].nstates)]
    readouts = zeros(T, length(apex(Ports)))
    readins  = zeros(T, length(apex(Ports)))
    ϕ = zeros(T, length(apex(S)))
    v = (u::AbstractVector, x::AbstractVector, p, t::Real) -> begin
        # length(p) == length(d[:, :con]) || error("Expected $(length(d[:, :con])) parameters, have $(length(p))")
        statefun(b) = state(u,b)
        inputfun(b) = readins[incident(d, b, :box)]
        fillreadouts!(readouts, d, xs, Ports, statefun)
        # communicate readouts to the ports at the other end of the wires, external connections directly fill ports
        readins .= 0
        fillreadins!(readins, d, readouts)
        readins[d[:, :con]] .+= x
        fillstates!(ϕ, d, xs, S, statefun, inputfun, p, t)
        return ϕ
    end
    function readout(u::AbstractVector)
        statefun(b) = state(u,b)
        fillreadouts!(readouts, d, xs, Ports, statefun)
        return readouts[d[:, :con]]
    end
    return Machine( nparts(d, :OuterPort), apex(S).set, v, readout)
end

"""    oapply(d::OpenCPortGraph, m::AbstractMachine)

A version of `oapply` where each box of `d` is filled with the machine `m`.
"""
function oapply(d::OpenCPortGraph, x::AbstractMachine)
    oapply(d, collect(repeated(x, nparts(d, :Box))))
end

"""    barbell(n::Int)

Constructs an open CPG with two boxes each with `n` ports. The ``i``th ports on each box are connected.
"""
barbell(k::Int) = begin
  g = OpenCPortGraph()
  add_parts!(g, :Box, 2)
  add_parts!(g, :Port, 2k; box=[fill(1,k); fill(2,k)])
  add_parts!(g, :Wire, k; src=1:k, tgt=k+1:2k)
  add_parts!(g, :Wire, k; tgt=1:k, src=k+1:2k)
  return g
end

meshpath(n::Int) = begin
    gt = @acset OpenCPortGraph begin
        Box = 1
        Port = 3
        Wire = 0
        OuterPort = 2
        box= ones(Int, 3)
        con= [3,2]
    end
    gm = @acset OpenCPortGraph begin
        Box = 1
        Port = 4
        Wire = 0
        OuterPort = 2
        box= ones(Int, 4)
        con= [4,2]
    end
    subs = [gt]
    for i in 2:n-1
        push!(subs, gm)
    end
    push!(subs, gt)
    X = coproduct(subs)
    for i in 1:n-1
        xi = subs[i]
        xj = subs[i+1]
        p = legs(X)[i][:Port](nparts(xi, :Port)-1)
        q = legs(X)[i+1][:Port](1)
        add_parts!(apex(X), :Wire, 2, src=[p,q], tgt=[q,p])
    end
    c₁ = apex(X)[1:2:nparts(apex(X),:OuterPort) ,:con]
    c₂ = apex(X)[2:2:nparts(apex(X),:OuterPort) ,:con]
    apex(X)[:,:con] = vcat(c₁,c₂)
    return X
end

function gridpath(n::Int, width::Int)
    node = @acset OpenCPortGraph begin
        Box = 1
        Port = 0
        Wire = 0
        box = 1
    end
    add_parts!(node, :Port, 2width, box=1)
    X = coproduct(collect(repeated(node, n)))
    L = legs(X)
    A = apex(X)
    for i in 1:n-1
        for j in 1:width
            s = L[i][:Port](j)
            t = L[i+1][:Port](j+width)
            add_part!(A, :Wire, src=s, tgt=t)
            add_part!(A, :Wire, src=t, tgt=s)
        end
    end
    upstream = L[1][:Port](width+1:2width)
    add_parts!(A, :OuterPort, width, con=upstream)
    downstream = L[end][:Port](1:width)
    add_parts!(A, :OuterPort, width, con=downstream)
    return X
end

grid(n::Int, m::Int) = ocompose(apex(gridpath(n,m)), collect(repeated(apex(meshpath(m)), n)))

end


text:
using DataFrames
using Random
using XLSX
using StructArrays
using StatsBase
using CSV
using CPLEX

include("../src/JOH.jl")
include("../MDMKP/MDMKP.jl")

""" create a variety of SSIT methods. Accept a parameter to multiply each time
limit by. """
function make_SSIT_methods(m=60; n_threads=6)
    [
        JOH.Matheur.SSIT.make_SSIT_method(
			[.001, .005, .01, .02, .05],
			[m*5, m*5, m*5, m*5, m*5],
			"even time", n_threads)
        JOH.Matheur.SSIT.make_SSIT_method(
			[.001, .001, .001, .001, .001],
			[m*5,m*5,m*5,m*5,m*5],
			"one tolerance", n_threads)
        JOH.Matheur.SSIT.make_SSIT_method(
			[.001, .001, .005, .005, .005],
			[m*5, m*5, m*5, m*5, m*5],
			"tight tolerances", n_threads)
        JOH.Matheur.SSIT.make_SSIT_method(
			[.001, .005, .01, .02, .05],
			[m*2, m*4, m*5, m*6, m*8],
            "increasing time", n_threads)
        JOH.Matheur.SSIT.make_SSIT_method(
			[.001, .005, .01, .02, .05],
			[m*8, m*6, m*5, m*4, m*2],
            "decreasing time", n_threads)
    ]
end

struct MethodProblemResult
	method_name
	problem_id
	lowest_gap
	highest_reached_tolerance
	total_time
	last_phase_time
	cplex_obj
	true_obj
	infeasibility
	n_phases
end

MethodProblemResult(method, problem, solution, last_row) = MethodProblemResult(
	method.name,
	problem.id.id,
	last_row[!, :gap],
	last_row[!, :rtol],
	last_row[!, :elapsed_time],
	last_row[!, :solve_time],
	last_row[!, :objective],
	solution.objective,
	solution.infeasibility,
	last_row[!, :index])

mutable struct ExperimentResults
	problem_ids::Vector{MDMKP.Problem_ID}
	SSIT_phases::Vector{DataFrame}
	method_problem_results::Vector{MethodProblemResult}
end

ExperimentResults() = ExperimentResults([], [], [])

function flatten_ssit(df::DataFrame, tolerances)
	times = []
	gaps = []

	n_rows = length(df[!, 1])

	for i in 1:n_rows
		push!(times, df[i, :][:elapsed_time])
		push!(gaps, df[i, :][:gap])
	end

	flat_df = DataFrame()
	for i in 1:n_rows
		flat_df[!, Symbol("phase $i time")] = [times[i]]
		flat_df[!, Symbol("phase $i gap")] = [gaps[i]]
	end
	flat_df[!, :termination] = [last(df)[:term_stat]]
	flat_df[!, :lowest_gap] = [last(df)[:gap]]
	flat_df
end

function include_aux_data(df::DataFrame, method, problem_id)
	df[!, :method] .= method.name
	for field in fieldnames(typeof(problem_id))
		val = getfield(problem_id, field)
		df[!, Symbol("problem_$(field)")] = [val]
	end
	df
end

function include_sol_data(df, ssit_phases, solution_constructor, problem)
	lp = last(ssit_phases, 1)
	solution = solution_constructor(lp[!, :bitarr][1], problem)
	df[!, :objective] = [solution.objective]
	df[!, :infeasibility] = [solution.infeasibility]
	df[!, :solution] = [solution.value]
	df
end

function summarize_ssit(ssit_phases::DataFrame, method, problem, solution)
	df = include_aux_data(flatten_ssit(ssit_phases, method.tolerances),
		method, problem.id)
	df = include_sol_data(df, ssit_phases, solution, problem)
	df
end

function generate_comparison_data(
		methods::Vector{JOH.Matheur.SSIT.SSIT_method},
		problem_groups::Vector{Vector{T}},
		solution; results_dir="results") where T <: JOH.Problem

	results = []
	index = 1
	for (method, problem_group) in zip(methods, problem_groups)
		for problem in problem_group
			model = MDMKP.create_MIPS_model(problem, CPLEX.Optimizer)
			ssit_phases = JOH.Matheur.evaluate(model, method, id=problem.id)
			result_df = summarize_ssit(ssit_phases, method, problem, solution)
			CSV.write("$(results_dir)/$(index).csv", result_df)
			index += 1

			push!(results, result_df)
		end
	end
	results
end

function ba_rep(ba)
	join([b == 1 ? "1" : "0" for b in ba], "")
end

function split_problems(problems; n_groups=5, datasets=1:9, tightnesses=[.25, .5, .75],
			cases=1:6)::Vector{Vector{JOH.Problem}}
	groups = [[] for _ in 1:n_groups]
	for dataset in datasets
		for tightness in tightnesses
			for case in cases
				random_assignment_order = shuffle(collect(1:n_groups))
				subset = filter(x->x.id.tightness==tightness
					&& x.id.dataset==dataset
					&& x.id.case==case, problems)
				j = 1
				for i in random_assignment_order
					push!(groups[i], subset[j])
					j += 1
				end
			end
		end
	end
	groups
end

function record_dataset(datasets, methods, problems, solution, results_dir="results")
	problems = filter(x->x.id.dataset in datasets, problems)
	problem_groups = split_problems(problems, datasets=datasets)

	generate_comparison_data(methods, problem_groups, solution, results_dir=results_dir)
end

all_problems = MDMKP.load_folder()
ssit_methods = make_SSIT_methods()

res_dir = "results/Gurobi_test"
mkpath(res_dir)
data = record_dataset(1:9, ssit_methods, all_problems, MDMKP.MDMKP_Sol, res_dir)

max_stars_repo_name: Dibillilia/JuliaOptHeuristics

max_issues_repo_name: Dibillilia/JuliaOptHeuristics

max_forks_repo_name: Dibillilia/JuliaOptHeuristics


text:
module ElPhDynamics

using LinearAlgebra
using FFTW
using TOML

#######################
## INCLUDING MODULES ##
#######################

include("Utilities.jl")

include("UnitCells.jl")

include("Lattices.jl")

include("Checkerboard.jl")

include("IterativeSolvers.jl")

include("TimeFreqFFTs.jl")

include("Models.jl")

include("KPMPreconditioners.jl")

include("InitializePhonons.jl")

include("PhononAction.jl")

include("FourierAcceleration.jl")

include("LangevinDynamics.jl")

include("HMC.jl")

include("SpecialUpdates.jl")

include("GreensFunctions.jl")

include("SimulationParams.jl")

include("Measurements.jl")

include("MuFinder.jl")

include("RunSimulation.jl")

include("SimulationSummary.jl")

include("ProcessInputFile.jl")

####################################
## DEFINING HIGHET LEVEL FUNCTION ##
##     TO RUN A SIMULATION        ##
####################################

using ..RunSimulation: run_simulation!
using ..ProcessInputFile: process_input_file, process_checkpoint, initialize_model
using ..Models: read_phonons!
using ..MuFinder: MuTuner, update_μ!
using ..SimulationSummary: write_simulation_summary!

export simulate, load_model

"""
Highest level function used to run a QMC simulation.
To run a simulation execute the following command:
`julia -O3 -e "using ElPhDynamics; simulate(ARGS)" -- input.toml`
"""
function simulate(args)

    # set the number of threads used by BLAS and FFTW libraries
    BLAS.set_num_threads(1)
    FFTW.set_num_threads(1)

    ######################
    ## READ CONFIG FILE ##
    ######################

    # getting iput filename
    config_file = args[1]

    # read input file
    input = TOML.parsefile(config_file)

    # construct data folder name
    filepath   = input["simulation"]["filepath"]
    foldername = input["simulation"]["foldername"]
    if length(args)==1
        datafolder = name_datafolder(filepath,foldername)
    else
        id = parse(Int,args[2])
        datafolder = name_datafolder(filepath,foldername,id)
    end
    input["simulation"]["datafolder"] = datafolder

    ###########################
    ## INITIALIZE SIMULATION ##
    ###########################

    if isdir(datafolder) # resume previous simulation

        # extract state from checkpoint
        (model, Gr, μ_tuner, sim_params, simulation_dynamics, burnin_dynamics,
         burnin_reflect_update, sim_reflect_update, burnin_swap_update, sim_swap_update,
         fourier_accelerator, preconditioner, container, burnin_start, sim_start, sim_stats) = process_checkpoint(input)

    else # start new simulation

        # initialize new simulation
        (model, Gr, μ_tuner, sim_params, simulation_dynamics, burnin_dynamics,
         burnin_reflect_update, sim_reflect_update, burnin_swap_update, sim_swap_update,
         fourier_accelerator, preconditioner, container, burnin_start, sim_start, sim_stats) = process_input_file(config_file,input)
    end

    ####################
    ## RUN SIMULATION ##
    ####################

    if sim_start < sim_params.nsteps
        sim_stats = run_simulation!(model, Gr, μ_tuner, sim_params, simulation_dynamics, burnin_dynamics,
                                    sim_reflect_update, burnin_reflect_update, sim_swap_update, burnin_swap_update,
                                    fourier_accelerator, container, preconditioner, sim_stats, burnin_start, sim_start)
    end

    #################################
    ## SUMARIZE SIMULATION RESULTS ##
    #################################

    if sim_start < sim_params.nsteps
        write_simulation_summary!(datafolder)
    end

    return nothing
end

"""
Pass this function a directory name generated by a simulation and it returns a HolsteinModel object
with the phonons intialized to the final phonon configuration sampled in the simulation.
"""
function load_model(dir::String)

    files  = readdir(dir)
    config = findall(f -> endswith(f, r"\.toml|\.TOML"), files)
    phonon = findall(f -> endswith(f, "_config.out"), files)
    @assert length(config) == length(phonon) == 1
    config_file = joinpath(dir, files[config[1]])
    phonon_file = joinpath(dir, files[phonon[1]])
    input = TOML.parsefile(config_file)
    input["simulation"]["datafolder"] = dir
    model = initialize_model(input)
    read_phonons!(model, phonon_file)

    return model
end

###########################################
## PRIVATE METHODS THAT ARE NOT EXPORTED ##
###########################################

"""
Name the data folder.
"""
function name_datafolder(filepath::String,foldername::String,id::Int=0)

    # data folder name without id at end
    key = joinpath(filepath,foldername)

    if id>0 # if ID given then attach to end of key
        datafolder = "$key-$id"
    else # determine the ID
        id = 1
        while true
            datafolder = "$key-$id"
            if !isdir(datafolder)
                break
            else
                id += 1
            end
        end
    end

    return datafolder
end

end

max_stars_repo_name: cohensbw/Langevin

max_issues_repo_name: cohensbw/ElPhDynamics

max_forks_repo_name: cohensbw/ElPhDynamics


text:
@with_kw struct Options
  nbest :: Int = 5
  besttol :: Float64 = 1e-4
  maxtrials :: Int = 100
  fine :: Int = 100
  p0_range :: Vector{Float64} = [ 0., 0. ]
  debug :: Bool = false
end
export Options


text:
include("Granular_Module.jl")
include("calculate_Pressure.jl")
include("Min_Energy.jl")
include("Deform_Config.jl")
include("Hertz.jl")
using JuLIP
using Hertz
using Minimise_Energy
using Deform
using Pressure
# =====================================
# ----------- Initializaiton ----------
# =====================================

tic()

Grain_Num = 1024
α = 2.0
ϕ_Initial = 0.839
Lx = 1.0
Ly = 1.0
ratio = 1.4

# read the initial infomation from the initial_config.txt file

Initial_Info = readdlm("initial_config.txt")
x = Initial_Info[:, 1]'
y = Initial_Info[:, 2]'
z = zeros(Grain_Num)'

grain_type = Initial_Info[:, 3] + 1
ρ = sqrt(2.0*ϕ_Initial*Lx*Ly/Grain_Num/π/(1+ratio^2))
M = ρ*(0.4*grain_type + 0.6)

X = [x; y; z]
X = vecs(X)

Config_Initial = gran.setup_cell(Grain_Num, ϕ_Initial, α)
Config_Initial.X .= X
Config_Initial.Z = grain_type
Config_Initial.M = M

while true
    deformType = "grow"
    deformValue = 1.0
    Config_Initial = Deform.deform_Config(Config_Initial, deformType, deformValue)
    V_Initial = get_data(Config_Initial, "Energy")
    if V_Initial < 1.00001 * eps()
        break
    end
end

time0 = toc()

# ===============================================================================
# ---------- Find the Jamming State Roughly by Letting the Grains Grow ----------
# ===============================================================================
# setup the growth parameter for fetching the rough jamming configuration

tic()

ϕ_Target = 0.843  # this is a volume fraction which is already known bigger than the critical volume fraction
N_Grow = 50      # this is the growth step
ϕ_Increment = (ϕ_Target - ϕ_Initial)/N_Grow
 # the growth history
Config_Grow_Hist = Array{AbstractAtoms, 1}(N_Grow)
V_Grow_Hist = zeros(N_Grow, 1)
p_Grow_Hist = zeros(N_Grow, 1)

deformType = "phi"
Config_Temp1= deepcopy(Config_Initial)
ϕ_Temp = get_data(Config_Temp1, "Phi")

 # setup the signal for jamming state
Jamming_Signal = 0
Jamming_Index = 0                              # if you need to use a variable out of loop/control body
Config_RoughJamming = deepcopy(Config_Initial) # you should declare it firstly

for i = 1:N_Grow
    deformValue = ϕ_Temp + ϕ_Increment
    Config_Temp2 = Deform.deform_Config(Config_Temp1, deformType, deformValue)
    Config_Grow_Hist[i] = deepcopy(Config_Temp2)

    V_Grow_Hist[i] = get_data(Config_Temp2, "Energy")
    p_Grow_Hist[i] = trace(get_data(Config_Temp2, "Pressure"))/get_data(Config_Temp2, "Dimension")

    if (Jamming_Signal == 0) & (V_Grow_Hist[i] > 5.0e-16) & (p_Grow_Hist[i] > 1e-9)
        Config_RoughJamming = deepcopy(Config_Grow_Hist[i])
        Jamming_Index = i
        Jamming_Signal = 1
    end

    if V_Grow_Hist[i] < 5.0e-16
        Jamming_Signal = 0
    end

    Config_Temp1 = deepcopy(Config_Temp2)
    ϕ_Temp = get_data(Config_Temp1, "Phi")

end

# ==============================================================================
# ---------- Find the Jamming State Accurately Using Bisection Method ----------
# ==============================================================================

ϕ_Jamming = get_data(Config_RoughJamming, "Phi")
Config_Jamming = deepcopy(Config_RoughJamming)

ϕ_Temp = get_data(Config_Grow_Hist[Jamming_Index - 1], "Phi")

deformType = "phi"

# the bisection method history
V_Mid_Hist = Array{Float64, 1}(0)
p_Mid_Hist = Array{Float64, 1}(0)
ϕ_Mid_Hist = Array{Float64, 1}(0)
Config_Mid_Hist = Array{AbstractAtoms, 1}(0)

while (ϕ_Jamming - ϕ_Temp) > 5.0e-8
    ϕ_Mid = (ϕ_Jamming + ϕ_Temp)/2
    deformValue = ϕ_Mid
    Config_Mid = deform_Config(Config_Jamming, deformType, deformValue)
    V_Mid = get_data(Config_Mid, "Energy")
    p_Mid = trace(get_data(Config_Mid, "Pressure"))/get_data(Config_Mid, "Dimension")
    push!(V_Mid_Hist, V_Mid)
    push!(p_Mid_Hist, p_Mid)
    push!(ϕ_Mid_Hist, ϕ_Mid)
    push!(Config_Mid_Hist, Config_Mid)

    if V_Mid > 5.0e-16
        ϕ_Jamming = ϕ_Mid
        Config_Jamming = deepcopy(Config_Mid)
    else
        ϕ_Temp = ϕ_Mid
    end

end
time1 = toc()

# ==============================================================================
# ---------- Test the Scaling Law (p vs ϕ) ----------
# Pick some ϕ greater than ϕ_Jamming and calculate the pressure correspondingly
# ==============================================================================

tic()
N_ϕ = 51

ϕ_Test1 = ϕ_Jamming + 10.^linspace(-6, -2, N_ϕ)
ϕ_Test1_Hist = Array{Float64, 1}(0)

Config_Test1 = Array{AbstractAtoms, 1}(N_ϕ) # This is the first main test so ...
Config_Test1_Hist = Array{AbstractAtoms, 1}(0)

Pressure_Test1 = Array{Float64, 1}(N_ϕ)
Pressure_Test1_Hist = Array{Float64, 1}(0)

V_Test1 = Array{Float64, 1}(N_ϕ)
V_Test1_Hist = Array{Float64, 1}(0)

deformType = "phi"

Config_Temp = deepcopy(Config_Jamming)
for i = 1:N_ϕ
    deformValue = ϕ_Test1[i]
    ϕ_diff = deformValue - get_data(Config_Temp, "Phi")
    ϕ_tol = 1.0e-4
    if ϕ_diff > ϕ_tol
        k = ceil(ϕ_diff/ϕ_tol)
        ϕ_Increment = ϕ_diff/k
        for j = 1:k
            deformValue = get_data(Config_Temp, "Phi") + ϕ_Increment
            Config_Test1[i] = deform_Config(Config_Temp, deformType, deformValue)
            push!(ϕ_Test1_Hist, deformValue)
            push!(V_Test1_Hist, get_data(Config_Temp, "Energy"))
            push!(Pressure_Test1_Hist, trace(get_data(Config_Temp, "Pressure"))/get_data(Config_Temp, "Dimension"))
            push!(Config_Test1_Hist, Config_Temp)
            Config_Temp = deepcopy(Config_Test1[i])
        end
    else
        Config_Test1[i] = deform_Config(Config_Temp, deformType, deformValue)
        push!(ϕ_Test1_Hist, deformValue)
        push!(V_Test1_Hist, get_data(Config_Temp, "Energy"))
        push!(Pressure_Test1_Hist, trace(get_data(Config_Temp, "Pressure"))/get_data(Config_Temp, "Dimension"))
        push!(Config_Test1_Hist, Config_Temp)
        Config_Temp = deepcopy(Config_Test1[i])
    end

    Pressure_Test1[i] = trace(get_data(Config_Temp, "Pressure"))/get_data(Config_Temp, "Dimension")
    V_Test1[i] = get_data(Config_Temp, "Energy")

end
time3 = toc()