text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct ServiceAssociationLinkPropertiesFormat2 <: SwaggerModel
    linkedResourceType::Any # spec type: Union{ Nothing, String } # spec name: linkedResourceType
    link::Any # spec type: Union{ Nothing, String } # spec name: link
    provisioningState::Any # spec type: Union{ Nothing, ProvisioningState } # spec name: provisioningState
    allowDelete::Any # spec type: Union{ Nothing, Bool } # spec name: allowDelete
    locations::Any # spec type: Union{ Nothing, Vector{String} } # spec name: locations

    function ServiceAssociationLinkPropertiesFormat2(;linkedResourceType=nothing, link=nothing, provisioningState=nothing, allowDelete=nothing, locations=nothing)
        o = new()
        validate_property(ServiceAssociationLinkPropertiesFormat2, Symbol("linkedResourceType"), linkedResourceType)
        setfield!(o, Symbol("linkedResourceType"), linkedResourceType)
        validate_property(ServiceAssociationLinkPropertiesFormat2, Symbol("link"), link)
        setfield!(o, Symbol("link"), link)
        validate_property(ServiceAssociationLinkPropertiesFormat2, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        validate_property(ServiceAssociationLinkPropertiesFormat2, Symbol("allowDelete"), allowDelete)
        setfield!(o, Symbol("allowDelete"), allowDelete)
        validate_property(ServiceAssociationLinkPropertiesFormat2, Symbol("locations"), locations)
        setfield!(o, Symbol("locations"), locations)
        o
    end
end # type ServiceAssociationLinkPropertiesFormat2

const _property_map_ServiceAssociationLinkPropertiesFormat2 = Dict{Symbol,Symbol}(Symbol("linkedResourceType")=>Symbol("linkedResourceType"), Symbol("link")=>Symbol("link"), Symbol("provisioningState")=>Symbol("provisioningState"), Symbol("allowDelete")=>Symbol("allowDelete"), Symbol("locations")=>Symbol("locations"))
const _property_types_ServiceAssociationLinkPropertiesFormat2 = Dict{Symbol,String}(Symbol("linkedResourceType")=>"String", Symbol("link")=>"String", Symbol("provisioningState")=>"ProvisioningState", Symbol("allowDelete")=>"Bool", Symbol("locations")=>"Vector{String}")
Base.propertynames(::Type{ ServiceAssociationLinkPropertiesFormat2 }) = collect(keys(_property_map_ServiceAssociationLinkPropertiesFormat2))
Swagger.property_type(::Type{ ServiceAssociationLinkPropertiesFormat2 }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ServiceAssociationLinkPropertiesFormat2[name]))}
Swagger.field_name(::Type{ ServiceAssociationLinkPropertiesFormat2 }, property_name::Symbol) =  _property_map_ServiceAssociationLinkPropertiesFormat2[property_name]

function check_required(o::ServiceAssociationLinkPropertiesFormat2)
    true
end

function validate_property(::Type{ ServiceAssociationLinkPropertiesFormat2 }, name::Symbol, val)
end


text:
#####
##### Dependencies
#####

import Base.display

#####
##### Structs
#####

@doc raw"""
"""
struct KernelData{K<:AbstractCorrelation}
    kern::K # Kernel Correlation function
    Ku::Matrix{Float64} # Correlation at locations U
    Kux::DataMatrix{Float64} # Cross correlation of observation locations
    rKx::IncrementalCholesky{Float64} # Reduced kernel system at locations X
end

#####
##### Construction
#####

function kernel_data!(kern, Ku, Kux, rKx, U)
    correlation!(Ku, kern, U)
    return KernelData(kern, Ku, Kux, rKx)
end

function valid_update(kd::KernelData, k)
    return valid_update(kd.Kux, k) && valid_update(kd.rKx, k)
end

function view_next(kd::KernelData, sd)
    k = size(sd.X, 2) - size(kd.Kux, 2)
    valid_update(kd, k) || throw(ErrorException("Size cannot exceed capacity"))
    Kxx′, Kx′ = view_next(kd.rKx, k)
    return (Kux′ = view_next(kd.Kux, k), Kxx′ = Kxx′, Kx = Kx′)
end

function compute_next!(kd::KernelData, sd)
    n = size(kd.Kux, 2)
    k = size(sd.X, 2) - n
    @views X, X′ = sd.X[:, 1:n], sd.X[:, n+1:end]
    @views W, W′ = sd.W[:, 1:n], sd.W[:, n+1:end]
    Kux′, Kxx′, Kx′ = view_next(kd, sd)

    cross_correlation!(Kux′, kd.kern, sd.U, X′)
    cross_correlation!(Kxx′, kd.kern, X, X′)
    correlation!(Kx′, kd.kern, X′)

    reduced_update!(Kxx′, Kx′, kd.Ku, kd.Kux, Kux′, W, W′)

    compute_next!(kd.Kux, k)
    compute_next!(kd.rKx, k)
    return nothing
end

function remove_last!(kd::KernelData, sd)
    k = size(kd.Kux, 2) - size(sd.X, 2)
    valid_update(kd, -k) || throw(ErrorException("Size cannot be negative"))
    remove_last!(kd.Kux, k)
    remove_last!(kd.rKx, k)
    return nothing
end

function reset!(kd::KernelData)
    remove_last!(kd.Kux, size(kd.Kux, 2))
    remove_last!(kd.rKx, size(kd.rKx, 2))
    return nothing
end

#####
##### Functions
#####

function Base.display(kd::KernelData)
    print("\n    Correlation Kernel:\n")
    display(kd.kern)
    println("\n    Correlation of U points:\n")
    display(kd.Ku)
    println("\n    Cross correlation between U and X:\n")
    display(kd.Kux)
    println("\n    Solver for the reduced kernel system:\n")
    display(kd.rKx)
    return nothing
end

function reduced_update!(Kxx′, Kx′, Ku, Kux, Kux′, W, W′)
    tmp = Kux′ .- Ku * W′

    mul!(Kxx′, W', tmp, -1, 1)
    mul!(Kxx′, Kux', W′, -1, 1)
    mul!(Kx′, (W′)', tmp, -1, 1)
    mul!(Kx′, (Kux′)', W′, -1, 1)

    return Kxx′, Kx′
end

function diag_reduced_schur!(Dx′, Kxx′, kd, Kux′, W, W′)
    tmp = Kux′ .- kd.Ku * W′

    mul!(Kxx′, W', tmp, -1, 1)
    mul!(Kxx′, kd.Kux', W′, -1, 1)
    ldiv!(get_chol(kd.rKx).L, Kxx′)

    for i in axes(Dx′, 1)
        @views Dx′[i] = 1 - dot(W′[:, i], tmp[:, i]) - dot(Kux′[:, i], W′[:, i])
    end

    r = similar(Dx′, 1, length(Dx′))
    sum!(abs2, r, Kxx′)
    @views Dx′ .-= r[:]
    return Dx′
end

ext: jl

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
using Random
using Test
using FFTW
using CUDA

Random.seed!(0)

FT = Float64
ε = eps(FT)

# Test both even and odd sizes
for N in (4, 5)
    A_cpu = rand(FT, N)
    A_gpu = CuArray(A_cpu)

    B_cpu = FFTW.fft(A_cpu)
    B_gpu = CUDA.CUFFT.fft(A_gpu)

    @test isapprox.(B_cpu, Vector(B_gpu), atol=ε) |> all

    C_cpu = FFTW.ifft(B_cpu)
    C_gpu = CUDA.CUFFT.ifft(B_gpu)

    @test isapprox.(C_cpu, Vector(C_gpu), atol=ε) |> all

    @test isapprox.(A_cpu, C_cpu, atol=ε) |> all
    @test isapprox.(Vector(A_gpu), Vector(C_gpu), atol=ε) |> all
end


text:
"""
$(TYPEDEF)
"""
struct PDEProblem{P,E,S} <: AbstractPDEProblem
  prob::P
  extrapolation::E
  space::S
end


text:
using Spelling
using Base.Test

@test spellcheck("radio") == true
@test spellcheck("aklsdfjlak") == false
@test score("test") > 5000
@test (w->w.match).(spellcheck_sentence("fox house spoon asdkflj horse")) == ["asdkflj"]


text:
# NOTE: Matrix implicitly assumes all agents have same number of actions
mutable struct PerStateMPStats
    agent_action_n::Matrix{Int64} # N X A
    agent_action_q::Matrix{Float64}
    edge_action_n::Matrix{Int64} # |E| X A^2
    edge_action_q::Matrix{Float64}
end

"""
Tracks the specific informations and statistics we need to use Max-Plus to coordinate_action
the joint action in Factored-Value MCTS. Putting parameters here is a little ugly but coordinate_action can't have them since VarEl doesn't use those args.

Fields:
    adjmatgraph::SimpleGraph
        The coordination graph as a LightGraphs SimpleGraph.

    message_iters::Int64
        Number of rounds of message passing.

    message_norm::Bool
        Whether to normalize the messages or not after message passing.

    use_agent_utils::Bool
        Whether to include the per-agent utilities while computing the best agent action (see our paper for details)

    node_exploration::Bool
        Whether to use the per-node UCB style bonus while computing the best agent action (see our paper for details)

    edge_exploration::Bool
        Whether to use the per-edge UCB style bonus after the message passing rounds (see our paper for details). One of this or node_exploration MUST be true for exploration.

    all_states_stats::Dict{AbstractVector{S},PerStateMPStats}
        Maps each joint state in the tree to the per-state statistics.
"""
mutable struct MaxPlusStatistics{S} <: CoordinationStatistics
    adjmatgraph::SimpleGraph
    message_iters::Int64
    message_norm::Bool
    use_agent_utils::Bool
    node_exploration::Bool
    edge_exploration::Bool # NOTE: One of this or node exploration must be true
    all_states_stats::Dict{S,PerStateMPStats}
end

function clear_statistics!(mp_stats::MaxPlusStatistics)
    empty!(mp_stats.all_states_stats)
end

function update_statistics!(mdp::JointMDP{S,A}, tree::FVMCTSTree{S,A,MaxPlusStatistics{S}},
                            s::S, ucb_action::A, q::AbstractFloat) where {S,A}

    update_statistics!(mdp, tree, s, ucb_action, ones(typeof(q), n_agents(mdp)) * q)
end

"""
Take the q-value from the MCTS step and distribute the updates across the per-node and per-edge q-stats as per the formula in our paper.
"""
function update_statistics!(mdp::JointMDP{S,A}, tree::FVMCTSTree{S,A,MaxPlusStatistics{S}},
                            s::S, ucb_action::A, q::AbstractVector{Float64}) where {S,A}

    state_stats = tree.coordination_stats.all_states_stats[s]
    nagents = n_agents(mdp)

    # Update per agent action stats
    for i = 1:nagents
        ac_idx = agent_actionindex(mdp, i, ucb_action[i])
        lock(tree.lock) do
            state_stats.agent_action_n[i, ac_idx] += 1
            state_stats.agent_action_q[i, ac_idx] +=
                (q[i] - state_stats.agent_action_q[i, ac_idx]) / state_stats.agent_action_n[i, ac_idx]
        end
    end

    # Now update per-edge action stats
    for (idx, e) in enumerate(edges(tree.coordination_stats.adjmatgraph))
        # NOTE: Need to be careful about action ordering
        # Being more general to have unequal agent actions
        edge_comp = (e.src,e.dst)
        edge_tup = Tuple(1:length(tree.all_agent_actions[c]) for c in edge_comp)
        edge_ac_idx = LinearIndices(edge_tup)[agent_actionindex(mdp, e.src, ucb_action[e.src]),
                                              agent_actionindex(mdp, e.dst, ucb_action[e.dst])]
        q_edge_value = q[e.src] + q[e.dst]

        lock(tree.lock) do
            state_stats.edge_action_n[idx, edge_ac_idx] += 1
            state_stats.edge_action_q[idx, edge_ac_idx] +=
                (q_edge_value - state_stats.edge_action_q[idx, edge_ac_idx]) / state_stats.edge_action_n[idx, edge_ac_idx]
        end
    end

    lock(tree.lock) do
        tree.coordination_stats.all_states_stats[s] = state_stats
    end

end

function init_statistics!(tree::FVMCTSTree{S,A,MaxPlusStatistics{S}}, planner::FVMCTSPlanner,
                          s::S) where {S,A}

    n_agents = length(s)

    # NOTE: Assuming all agents have the same actions here
    n_all_actions = length(tree.all_agent_actions[1])

    agent_action_n = zeros(Int64, n_agents, n_all_actions)
    agent_action_q = zeros(Float64, n_agents, n_all_actions)

    # Loop over agents and then actions
    # TODO: Need to define init_N and init_Q for single agent
    for i = 1:n_agents
        for (j, ac) in enumerate(tree.all_agent_actions[i])
            agent_action_n[i, j] = init_N(planner.solver.init_N, planner.mdp, s, i, ac)
            agent_action_q[i, j] = init_Q(planner.solver.init_Q, planner.mdp, s, i, ac)
        end
    end

    n_edges = ne(tree.coordination_stats.adjmatgraph)
    edge_action_n = zeros(Int64, n_edges, n_all_actions^2)
    edge_action_q = zeros(Float64, n_edges, n_all_actions^2)

    # Loop over edges and then action_i \times action_j
    for (idx, e) in enumerate(edges(tree.coordination_stats.adjmatgraph))
        edge_comp = (e.src, e.dst)
        n_edge_actions = prod([length(tree.all_agent_actions[c]) for c in edge_comp])
        edge_tup = Tuple(1:length(tree.all_agent_actions[c]) for c in edge_comp)

        for edge_ac_idx = 1:n_edge_actions
            ct_idx = CartesianIndices(edge_tup)[edge_ac_idx]
            edge_action = [tree.all_agent_actions[c] for c in edge_comp]

            edge_action_n[idx, edge_ac_idx] = init_N(planner.solver.init_N, planner.mdp, s, edge_comp, edge_action)
            edge_action_q[idx, edge_ac_idx] = init_Q(planner.solver.init_Q, planner.mdp, s, edge_comp, edge_action)
        end
    end

    state_stats = PerStateMPStats(agent_action_n, agent_action_q, edge_action_n, edge_action_q)
    lock(tree.lock) do
        tree.coordination_stats.all_states_stats[s] = state_stats
    end
end

"""
Runs Max-Plus at the current state using the per-state MaxPlusStatistics to compute the best joint action with either or both of node-wise and edge-wise exploration bonus. Rounds of message passing are followed by per-node maximization.
"""
function coordinate_action(mdp::JointMDP{S,A}, tree::FVMCTSTree{S,A,MaxPlusStatistics{S}}, s::S,
                           exploration_constant::Float64=0.0, node_id::Int64=0) where {S,A}

    state_stats = lock(tree.lock) do
        tree.coordination_stats.all_states_stats[s]
    end
    adjgraphmat = lock(tree.lock) do
        tree.coordination_stats.adjmatgraph
    end
    k = tree.coordination_stats.message_iters
    message_norm = tree.coordination_stats.message_norm

    n_agents = length(s)
    state_agent_actions = [agent_actions(mdp, i, si) for (i, si) in enumerate(s)]
    n_all_actions = length(tree.all_agent_actions[1])
    n_edges = ne(tree.coordination_stats.adjmatgraph)

    # Init forward and backward messages and q0
    fwd_messages = zeros(Float64, n_edges, n_all_actions)
    bwd_messages = zeros(Float64, n_edges, n_all_actions)

    if tree.coordination_stats.use_agent_utils
        q_values = state_stats.agent_action_q / n_agents
    else
        q_values = zeros(size(state_stats.agent_action_q))
    end

    state_total_n = lock(tree.lock) do
        (node_id > 0) ? tree.total_n[node_id] : 1
    end

    # Iterate over passes
    for t = 1:k
        fnormdiff, bnormdiff = perform_message_passing!(fwd_messages, bwd_messages, mdp, tree.all_agent_actions,
                                    adjgraphmat, state_agent_actions, n_edges, q_values, message_norm,
                                    0, state_stats, state_total_n)

        if !tree.coordination_stats.use_agent_utils
            q_values = zeros(size(state_stats.agent_action_q))
        end

        # Update Q value with messages
        for i = 1:n_agents

            # need indices of all edges that agent is involved in
            nbrs = neighbors(tree.coordination_stats.adjmatgraph, i)

            edgelist = collect(edges(adjgraphmat))

            if tree.coordination_stats.use_agent_utils
                @views q_values[i, :] = state_stats.agent_action_q[i, :]/n_agents
            end
            for n in nbrs
                if Edge(i,n) in edgelist # use backward message
                    q_values[i,:] += bwd_messages[findfirst(isequal(Edge(i,n)), edgelist), :]
                elseif Edge(n,i) in edgelist
                    q_values[i,:] += fwd_messages[findfirst(isequal(Edge(n,i)), edgelist), :]
                else
                    @warn "Neither edge found!"
                end
            end
        end

        # If converged, break
        if isapprox(fnormdiff, 0.0) && isapprox(bnormdiff, 0.0)
            break
        end
    end # for t = 1:k

    # If edge exploration flag enabled, do a final exploration bonus
    if tree.coordination_stats.edge_exploration
        perform_message_passing!(fwd_messages, bwd_messages, mdp, tree.all_agent_actions,
                                 adjgraphmat, state_agent_actions, n_edges, q_values, message_norm,
                                 exploration_constant, state_stats, state_total_n)
    end # if edge_exploration

    # Maximize q values for agents
    best_action = Vector{eltype(A)}(undef, n_agents)
    for i = 1:n_agents

        # NOTE: Again can't just iterate over agent actions as it may be a subset
        exp_q_values = zeros(length(state_agent_actions[i]))
        if tree.coordination_stats.node_exploration
            for (idx, ai) in enumerate(state_agent_actions[i])
                ai_idx = agent_actionindex(mdp, i, ai)
                exp_q_values[idx] = q_values[i, ai_idx] + exploration_constant*sqrt((log(state_total_n + 1.0))/(state_stats.agent_action_n[i, ai_idx] + 1.0))
            end
        else
            for (idx, ai) in enumerate(state_agent_actions[i])
                ai_idx = agent_actionindex(mdp, i, ai)
                exp_q_values[idx] = q_values[i, ai_idx]
            end
        end

        # NOTE: Can now look up index in exp_q_values and then again look at state_agent_actions
        _, idx = findmax(exp_q_values)
        best_action[i] = state_agent_actions[i][idx]
    end

    return best_action
end

function perform_message_passing!(fwd_messages::AbstractArray{F,2}, bwd_messages::AbstractArray{F,2},
                                  mdp, all_agent_actions,
                                  adjgraphmat, state_agent_actions, n_edges::Int, q_values, message_norm,
                                  exploration_constant, state_stats, state_total_n) where {F}
    # Iterate over edges
    fwd_messages_old = deepcopy(fwd_messages)
    bwd_messages_old = deepcopy(bwd_messages)
    for (e_idx, e) in enumerate(edges(adjgraphmat))

        i = e.src
        j = e.dst
        edge_tup_indices = LinearIndices(Tuple(1:length(all_agent_actions[c]) for c in (i,j)))

        # forward: maximize sender
        # NOTE: Can't do enumerate as action set might be smaller
        # Need to look up global index of agent action and use that
        # Need to break up vectorized loop
        @inbounds for aj in state_agent_actions[j]
            aj_idx = agent_actionindex(mdp, j, aj)
            fwd_message_vals = zeros(length(state_agent_actions[i]))
            # TODO: Should we use inbounds here again?
            @inbounds for (idx, ai) in enumerate(state_agent_actions[i])
                ai_idx = agent_actionindex(mdp, i, ai)
                fwd_message_vals[idx] = q_values[i, ai_idx] - bwd_messages_old[e_idx, ai_idx] + state_stats.edge_action_q[e_idx, edge_tup_indices[ai_idx, aj_idx]]/n_edges + exploration_constant * sqrt( (log(state_total_n + 1.0)) / (state_stats.edge_action_n[e_idx, edge_tup_indices[ai_idx, aj_idx]] + 1) )
            end
            fwd_messages[e_idx, aj_idx] = maximum(fwd_message_vals)
        end

        @inbounds for ai in state_agent_actions[i]
            ai_idx = agent_actionindex(mdp, i, ai)
            bwd_message_vals = zeros(length(state_agent_actions[j]))
            @inbounds for (idx, aj) in enumerate(state_agent_actions[j])
                aj_idx = agent_actionindex(mdp, j, aj)
                bwd_message_vals[idx] = q_values[j, aj_idx] - fwd_messages_old[e_idx, aj_idx] + state_stats.edge_action_q[e_idx, edge_tup_indices[ai_idx, aj_idx]]/n_edges + exploration_constant * sqrt( (log(state_total_n + 1.0))/ (state_stats.edge_action_n[e_idx, edge_tup_indices[ai_idx, aj_idx]] + 1) )
            end
            bwd_messages[e_idx, ai_idx] = maximum(bwd_message_vals)
        end

        # Normalize messages for better convergence
        if message_norm
            @views fwd_messages[e_idx, :] .-= sum(fwd_messages[e_idx, :])/length(fwd_messages[e_idx, :])
            @views bwd_messages[e_idx, :] .-= sum(bwd_messages[e_idx, :])/length(bwd_messages[e_idx, :])
        end

    end # (idx,edges) in enumerate(edges)

    # Return norm of message difference
    return norm(fwd_messages - fwd_messages_old), norm(bwd_messages - bwd_messages_old)
end


text:
@testset "49.group-anagrams.jl" begin
    function check(input, reference)
        res = group_anagrams(input)
        res = sort(map(sort, res))
        ref = sort(map(sort, res))
        return res == ref
    end
    @test check(["eat","tea","tan","ate","nat","bat"], [["bat"],["nat","tan"],["ate","eat","tea"]])
    @test check([""], [""])
    @test check(["a"], ["a"])
end


text:
using Decimals
using Base.Test

# Addition
@test decimal(0.1) + 0.2 == 0.1 + decimal(0.2) == decimal(0.1) + decimal(0.2) == decimal(0.3)
@test decimal([0.1 0.2]) + [0.3 0.1] == decimal([0.4 0.3])
@test decimal(2147483646) + decimal(1) == decimal(2147483647)
@test Decimal(1,3,-2) + decimal("0.2523410412138103") == Decimal(0,2223410412138103,-16)

# Subtraction
@test decimal(0.3) - 0.1 == 0.3 - decimal(0.1)
@test 0.3 - decimal(0.1) == decimal(0.3) - decimal(0.1)
@test decimal(0.3) - decimal(0.1) == decimal(0.2)
@test decimal([0.3 0.1]) - [0.1 0.5] == decimal([0.2 -0.4])

# Negation
@test -decimal([0.3 0.2]) == [-decimal(0.3) -decimal(0.2)]

# Multiplication
@test decimal(12.21) * decimal(2.12) == Decimal(0,258852,-4)
@test decimal(12.2112543) * decimal(2.121352) == Decimal(0,259043687318136,-13)
@test decimal(0.2) * 0.1 == 0.2 * decimal(0.1)
@test 0.2 * decimal(0.1) == decimal(0.02)
@test decimal(12.34) * 0.1234 == 12.34 * decimal(0.1234)
@test 12.34 * decimal(0.1234) == decimal(1.522756)
@test decimal(0.21084210) * -2 == -2 * decimal(0.21084210)
@test -2 * decimal(0.21084210) == decimal(-0.4216842)
# @test decimal([0.1 0.2]) .* [3 2] == [Decimal(0,3,-1) Decimal(0,4,-1)]
# @test decimal([0.1 0.2]) .* [1.2 2.1] == [Decimal(0,12,-2) Decimal(0,42,-2)]
# @test decimal([0.1,0.2]) * [1.2 2.1] == [Decimal(0,12,-2) Decimal(0,21,-2); Decimal(0,24,-2) Decimal(0,42,-2)]
@test Decimal(0, 2, -1) * 0.0 == decimal(0)
@test decimal([0.3, 0.6]) .* 5 == [decimal(0.3)*5, decimal(0.6)*5]

# Inversion
@test inv(Decimal(0, 1, -1)) == Decimal(0, 1, 1)
@test inv(Decimal(0, 1, 1)) == Decimal(0, 1, -1)
@test inv(Decimal(1, 2, -1)) == Decimal(1, 5, 0)
@test inv(Decimal(1, 5, 0)) == Decimal(1, 2, -1)
@test inv(Decimal(0, 2, -2)) == Decimal(0, 5, 1)
@test inv(Decimal(0, 5, 1)) == Decimal(0, 2, -2)
@test inv(Decimal(1, 4, -1)) == Decimal(1, 25, -1)
@test inv(Decimal(1, 25, -1)) == Decimal(1, 4, -1)
@test inv(d) == map(inv, d)

# Division
@test decimal(0.2) / decimal(0.1) == decimal(2)
@test decimal(0.3) / decimal(0.1) == Decimal(0,3,0)
@test [decimal(0.3)/decimal(0.1), decimal(0.6)/decimal(0.1)] == [decimal(0.3), decimal(0.6)] ./ decimal(0.1)
@test [decimal(0.3)/0.1, decimal(0.6)/0.1] == [decimal(0.3), decimal(0.6)] ./ 0.1


text:
using GlobalDiffEq, OrdinaryDiffEq, LinearAlgebra
using Test

l = 1.0                             # length [m]
m = 1.0                             # mass[m]
g = 9.81                            # gravitational acceleration [m/s²]

function pendulum!(du,u,p,t)
    du[1] = u[2]                    # θ'(t) = ω(t)
    du[2] = -3g/(2l)*sin(u[1]) + 3/(m*l^2)*p(t) # ω'(t) = -3g/(2l) sin θ(t) + 3/(ml^2)M(t)
end

θ₀ = 0.01                           # initial angular deflection [rad]
ω₀ = 0.0                            # initial angular velocity [rad/s]
u₀ = [θ₀, ω₀]                       # initial state vector
tspan = (0.0,10.0)                  # time interval

M = t->0.1sin(t)                    # external torque [Nm]

prob = ODEProblem(pendulum!,u₀,tspan,M)

v0 = solve(prob, Tsit5(), dt=0.1, reltol=1e-12, abstol=0).(1:10)
ve = solve(prob, GlobalRichardson(SSPRK33()), dt=0.2, reltol=1e-12, abstol=0).(1:10)
@test norm(ve - v0) / norm(v0) < 1e-10


text:
@testset "Transition configuration tests" begin
    mpSim = MPsim()
    filePath = joinpath( "base", "baseConfig" )
    fileName = string( filePath, ".xlsx" )
    catFileName = joinpath( "base", "catalogue.xlsx" )

    XLSX.openxlsx( fileName ) do xf
        XLSX.openxlsx( catFileName ) do catXF
            MPIO.readAttributes( mpSim, xf, catXF )
            MPIO.readBaseNodes( mpSim, xf, catXF )
            MPIO.readCompoundNodes( mpSim, xf, catXF )
        end  # XLSX.openXLSX( catFileName ) do xf

        MPIO.readTransitions( mpSim, xf )
    end  # XLSX.openxlsx( fileName ) do xf

    @test length( mpSim.transitionTypeOrder ) == 5
    @test ( length( mpSim.recruitmentByName ) == 1 ) &&
        ( length( mpSim.recruitmentByTarget ) == 3 )
    @test ( length( mpSim.transitionsByName ) == 4 ) &&
        ( length( mpSim.transitionsBySource ) == 6 ) &&
        ( length( mpSim.transitionsByTarget ) == 4 )
end  # @testset "Transition configuration tests"


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    SimplificationMethod

A method for simplifying geometries.
"""
abstract type SimplificationMethod end

"""
    simplify(geometry, method)

Simplify `geometry` with given `method`.
"""
function simplify end

function simplify(multi::Multi, method::SimplificationMethod)
  Multi([simplify(geom, method) for geom in multi])
end

function simplify(polygon::Polygon, method::SimplificationMethod)
  c = [simplify(chain, method) for chain in chains(polygon)]
  PolyArea(c[1], c[2:end])
end

# ----------------
# IMPLEMENTATIONS
# ----------------

include("simplification/douglaspeucker.jl")


text:
module Select

using Nullables

export @select

# =========================================================================================
# Custom concurrency primitives needed to support `@select`
# ------------------

function isready_put(c::Channel, sibling_tasks)
    # TODO: To fix the circular dependency, I think it might be enough to just add a check
    # here that there is at least one ready task that _isn't_ one of our siblings! We can
    # take another argument to this function, which is the list of tasks, and cross-reference it?
    return if Base.isbuffered(c)
        length(c.data) != c.sz_max
    else
        # TODO: No this isn't enough. I need to do it for the _wait_ function, not the wait_put. :'(
        #@info sibling_tasks
        #@info "isready_put:" (!isempty(c.cond_take.waitq))#, collect(c.cond_take.waitq))
        !isempty(c.cond_take.waitq) && any(t->!in(t, sibling_tasks), c.cond_take.waitq)
    end
end

function wait_put(c::Channel, sibling_tasks)
    #isready_put(c, sibling_tasks) && return
    # TODO: Is this sufficiently thread-safe?
    lock(c)
    try
        while !isready_put(c, sibling_tasks)
            Base.check_channel_state(c)
            wait(c.cond_put)  # Can be cancelled while waiting here...
        end
    finally
        unlock(c)
    end
    nothing
end

isready_wait_nosibs(c::Channel, sibling_tasks) = n_avail_nosibs(c, sibling_tasks) > 0
function n_avail_nosibs(c::Channel, sibling_tasks)
    if Base.isbuffered(c)
        length(c.data)
    else
        #@info "isready_wait_nosibs:" (isempty(c.cond_put.waitq), collect(c.cond_put.waitq))
        length(filter(t->0==count(x->x==t, sibling_tasks), collect(c.cond_put.waitq)))
    end
end
function wait_nosibs(c::Channel, sibling_tasks)
    # I don't understand why its okay to access this outside the lock...?
    #isready_wait_nosibs(c, sibling_tasks) && return
    lock(c)
    try
        while !isready_wait_nosibs(c, sibling_tasks)
            Base.check_channel_state(c)
            wait(c.cond_wait)
        end
    finally
        unlock(c)
    end
    nothing
end

wait_select(c::Channel, parent_task, sibling_tasks) = wait_nosibs(c, sibling_tasks)
wait_select(c::Base.GenericCondition, parent_task, sibling_tasks) = wait_from_parent(c, parent_task)
wait_select(x, parent_task, sibling_tasks) = wait(x)

# ---- Conditions ---------
assert_parent_haslock(c::Base.GenericCondition, parent_task) = assert_parent_haslock(c.lock, parent_task)
assert_parent_haslock(l::ReentrantLock, parent_task) =
    (islocked(l) && l.locked_by === parent_task) ? nothing : Base.concurrency_violation()
assert_parent_haslock(l::Base.AlwaysLockedST, parent_task) =
    (islocked(l) && l.ownertid === parent_task) ? nothing : Base.concurrency_violation()

function wait_from_parent(c::Base.GenericCondition, parent_task)
    ct = current_task()
    # Note that the parent task is guaranteed to be blocking on us, so this is okay.
    assert_parent_haslock(c, parent_task)
    push!(c.waitq, ct)
    token = unlockall_from_parent(c.lock, parent_task)
    try
        return wait()
    catch
        ct.queue === nothing || Base.list_deletefirst!(ct.queue, ct)
        rethrow()
    finally
        # Note that now _this task_ gets the lock, so we can execute the remaining body w/ the lock
        Base.relockall(c.lock, token)
        # eww, manually re-assign the parent to own this lock
        c.lock.locked_by = parent_task
    end
end

function unlockall_from_parent(rl::ReentrantLock, parent_task)
    n = rl.reentrancy_cnt
    rl.locked_by === parent_task || error("unlock from wrong thread")
    n == 0 && error("unlock count must match lock count")
    lock(rl.cond_wait)
    rl.reentrancy_cnt = 0
    rl.locked_by = nothing
    if !isempty(rl.cond_wait.waitq)
        try
            notify(rl.cond_wait)
        catch
            unlock(rl.cond_wait)
            rethrow()
        end
    end
    unlock(rl.cond_wait)
    return n
end

# =========================================================================================

## Implementation of 'select' mechanism to block on the disjunction of
## of 'waitable' objects.

@enum SelectClauseKind SelectPut SelectTake SelectDefault

# Represents a single parsed select "clause" of a @select macro call.
# eg, the (channel |> value) part of
# @select if channel |> value
#    println(value)
# ...
# end
struct SelectClause{ChannelT, ValueT}
    kind::SelectClauseKind
    channel::Nullable{ChannelT}
    value::Nullable{ValueT}
end

const select_take_symbol = :|>
const select_put_symbol = :<|

#  A 'structured' select clause is one of the form "channel|>val" or
#  "channel<|val". All other clauses are considered "non-structured", meaning
#  the entire clause is assumed to be an expression that evaluates to a
#  conditional to which "_take!" will be applied.
is_structured_select_clause(clause::Expr) =
    clause.head == :call &&
    length(clause.args) == 3 &&
    clause.args[1] ∈ (select_take_symbol, select_put_symbol)

is_structured_select_clause(clause) = false

function parse_select_clause(clause)
    if is_structured_select_clause(clause)
        if clause.args[1] == select_take_symbol
            SelectClause(SelectTake, Nullable(clause.args[2]), Nullable(clause.args[3]))
        elseif clause.args[1] == select_put_symbol
            SelectClause(SelectPut, Nullable(clause.args[2]), Nullable(clause.args[3]))
        end
    else
        # Assume this is a 'take' clause whose return value isn't wanted.
        # To simplify the rest of the code to not have to deal with this special case,
        # the return value is assigned to a throw-away gensym.
        SelectClause(SelectTake, Nullable(clause), Nullable(gensym()))
    end
end

"""
`@select`
A select expression of the form:
```julia
@select begin
     clause1 => body1
     clause2 => body2
     _       => default_body
    end
end
```
Wait for multiple clauses simultaneously using a pattern matching syntax, taking a different action depending on which clause is available first.
A clause has three possible forms:
1) `event |> value`
If `event` is an `AbstractChannel`, wait for a value to become available in the channel and assign `take!(event)` to `value`.
if `event` is a `Task`, wait for the task to complete and assign `value` the return value of the task.
2) `event |< value`
Only suppored for `AbstractChannel`s. Wait for the channel to capabity to store an element, and then call `put!(event, value)`.
3) `event`
Calls `wait` on `event`, discarding the return value. Usable on any "waitable" events", which include channels, tasks, `Condition` objects, and processes.

If a default branch is provided, `@select` will check arbitrary choose any event which is ready and execute its body, or will execute `default_body` if none of them are.

Otherise, `@select` blocks until at least one event is ready.

For example,

```julia
channel1 = Channel()
channel2 = Channel()
task = @task ...
result = @select begin
    channel1 |> value => begin
            info("Took from channel1")
            value
        end
    channel2 <| :test => info("Put :test into channel2")
    task              => info("task finished")
end
```
"""
macro select(expr)
    clauses = Tuple{SelectClause, Any}[]
    # @select can operate in blocking or nonblocking mode, determined by whether
    # an 'else' clause is present in the @select body (in which case it will be
    # nonblocking).
    mode = :blocking
    for se in expr.args
        # skip line nodes
        isa(se, Expr) || continue
        # grab all the pairs
        if se.head == :call && se.args[1] == :(=>)
            if se.args[2] != :_
                push!(clauses, (parse_select_clause(se.args[2]), se.args[3]))
            else
                # The defaule case (_). If present, the select
                # statement is considered non-blocking and will return this
                # section if none of the other conditions are immediately available.
                push!(clauses, (SelectClause(SelectDefault, Nullable(), Nullable()), se.args[3]))
                mode = :nonblocking
            end
        elseif se.head != :block && se.head != :line
            # if we run into an expression that is not a block. line or pair throw an error
            throw(ErrorException("Selection expressions must be Pairs. Found: $(se.head)"))
        end
    end
    if mode == :nonblocking
        _select_nonblock_macro(clauses)
    else
        _select_block_macro(clauses)
    end
end
# These defintions allow for any condition-like object to be used
# with select.
# @select if x |> value  ... will ultimately insert an expression value=_take!(x).
_take!(c::AbstractChannel) = take!(c)
_take!(x) = fetch(x)
# @select if x <| value .... will ultimately inset value=put!(x), which currently
# is only meanginful for channels and so no underscore varirant is used here.
# These are used with the non-blocking variant of select, which will
# only work with channels and tasks. Arbitrary conditionals can't be supported
# since "wait" is level-triggered.
_isready(c::AbstractChannel) = isready(c)
_isready(t::Task) = istaskdone(t)

_wait_condition(c::AbstractChannel) = c.cond_wait
_wait_condition(x) = x

# helper function to place the default case in the proper position
function set_default_first!(clauses)
    default_pos = findall(clauses) do x
        clause, body = x
        clause.kind == SelectDefault
    end
    l = length(default_pos)
    l == 0 && return # bail out if there is no default case
    l  > 1 && throw(ErrorException("Select takes at most one default case. Found: $l"))
    # swap elements to sure make SelectDefault comes first
    clauses[1], clauses[default_pos[1]] = clauses[default_pos[1]], clauses[1]
    clauses
end

function _select_nonblock_macro(clauses)
    set_default_first!(clauses)
    branches = Expr(:block)
    for (clause, body) in clauses
        branch =
        if clause.kind == SelectPut
            channel_var = gensym("channel")
            channel_assignment_expr = :($channel_var = $(clause.channel|>get|>esc))
            :(if ($channel_assignment_expr; isready_put($channel_var, []))
                put!($channel_var, $(clause.value|>get|>esc))
                $(esc(body))
            end)
        elseif clause.kind == SelectTake
            channel_var = gensym("channel")
            channel_assignment_expr = :($channel_var = $(clause.channel|>get|>esc))
            :(if ($channel_assignment_expr; _isready($channel_var))
                $(clause.value|>get|>esc) = _take!($channel_var)
                $(esc(body))
            end)
        elseif clause.kind == SelectDefault
            :($(esc(body)))
        end

        # the next two lines build an if / elseif chain from the bottom up
        push!(branch.args, branches)
        branches = branch
    end
    :($branches)
end

# The strategy for blocking select statements is to create a set of "rival"
# tasks, one per condition. When a rival "wins" by having its conditional be
# the first available, it sends a special interrupt to its rivals to kill them.
# The interrupt includes the task where control should be resumed
# once the rival has shut itself down.
struct SelectInterrupt <: Exception
    parent::Task
end
# Kill all tasks in "tasks" besides  a given task. Used for killing the rivals
# of the winning waiting task.
function select_kill_rivals(tasks, myidx)
    #@info myidx
    for (taskidx, task) in enumerate(tasks)
        taskidx == myidx && continue
        #@info taskidx, task
        #if task.state == :waiting || task.state == :queued
            # Rival is blocked waiting for its channel; send it a message that it's
            # lost the race.
            Base.schedule(task, SelectInterrupt(current_task()), error=true)
        # TODO: Is this still a legit optimization?:
        # elseif task.state==:queued
        #     # Rival hasn't starting running yet and so hasn't blocked or set up
        #     # a try-catch block to listen for SelectInterrupt.
        #     # Just delete it from the workqueue.
        #     queueidx = findfirst(Base.Workqueue.==task)
        #     deleteat!(Base.Workqueue, queueidx)
        # end
    end
    #@info "done killing"
end
function _select_block_macro(clauses)
    branches = Expr(:block)
    body_branches = Expr(:block)
    clause_lock = gensym("clause_lock")
    lock_assignment_expr = :($clause_lock = Base.ReentrantLock())
    for (i, (clause, body)) in enumerate(clauses)
        channel_var = gensym("channel")
        value_var = clause.value|>get|>esc
        channel_declaration_expr = :(local $channel_var)
        channel_assignment_expr = :($channel_var = $(clause.channel|>get|>esc))
        if clause.kind == SelectPut
            isready_func = isready_put
            wait_for_channel =  :(wait_put($channel_var, tasks))
            mutate_channel =  :(put!($channel_var, $value_var))
            bind_variable = :(nothing)
        elseif clause.kind == SelectTake
            isready_func = _isready
            wait_for_channel =  :(wait_select($channel_var, maintask, tasks))
            mutate_channel =  :(_take!($channel_var))
            bind_variable = :($value_var = branch_val)
        end
        branch = quote
            tasks[$i] = @async begin
                $channel_declaration_expr
                try  # Listen for genuine errors to throw to the main task
                    $channel_assignment_expr

                    # ---- Begin the actual `wait_and_select` algorithm ----
                    # TODO: Is this sufficiently thread-safe?
                    # Listen for SelectInterrupt messages so we can shutdown
                    # if a rival's channel unblocks first.
                    try
                        #@info "Task $($i) about to wait"
                        $wait_for_channel

                        # TODO: Because of this gap, where no locks are held, it's possible
                        # that multiple tasks can be woken-up due to a `put!` or `take!` on
                        # a channel they were waiting for. Only once will proceed in this
                        # @select, but a channel running _outside this macro_ may yet proceed
                        # and cause a problem.. I think this is bad. Fix this (probably) by
                        # returning the lock to unlock from `wait_for_channel`.

                        # NOTE: This is _not a deadock_ because there is a global ordering
                        # to the locks: we _ALWAYS_ wait on the channel before waiting on
                        # the clause_lock. This invariant must not be violated.
                        #@info "Task $($i) about to lock"
                        lock($clause_lock)
                        # We got the lock, so run this task to completion.
                        try
                            #@info "Task $($i): got lock"
                            # This block is atomic, so it _shouldn't_ matter whether we kill
                            # rivals first or mutate_channel first. It only matters if one
                            # case is accidentally synchronizing w/ another case, which
                            # should be specifically prohibited (somehow).
                            # For now, I'm killing rivals first so that at least we'll get
                            # an exception, rather than a deadlock, if we end up waiting on
                            # our rival, sibling cases.
                            #@info "Task $($i): killing rivals"
                            select_kill_rivals(tasks, $i)

                            #@info "Task $($i): mutating"
                            event_val = $mutate_channel
                            #@info "Got event_val: $event_val"
                            put!(winner_ch, ($i, event_val))
                        finally
                            #@info "Task $($i)) unlock"
                            unlock($clause_lock)
                        end
                    catch err
                        if isa(err, SelectInterrupt)
                            #@info "CAUGHT SelectInterrupt: $err"
                            #yieldto(err.parent)  # TODO: is this still a thing we should do?
                            return
                        else
                            rethrow()
                        end
                    end
                catch err
                    Base.throwto(maintask, err)
                end
            end # if
        end # for
        push!(branches.args, branch)

        body_branch = :(if branch_id == $i; $bind_variable; $(esc(body)); end)
        # the next two lines build an if / elseif chain from the bottom up
        push!(body_branch.args, body_branches)
        body_branches = body_branch
    end
    quote
        winner_ch = Channel(1)
        tasks = Array{Task}(undef, $(length(clauses)))
        maintask = current_task()
        $lock_assignment_expr
        $branches # set up competing tasks
        (branch_id, branch_val) = take!(winner_ch) # get the id of the winning task
        $body_branches # execute the winning block in the original lexical context
    end
end
# The following methods are the functional (as opposed to macro) forms of
# the select statement.
function _select_nonblock(clauses)
    for (i, clause) in enumerate(clauses)
        if clause[1] == :put
            if isready_put(clause[2], [])
                return (i, put!(clause[2], clause[3]))
            end
        elseif clause[1] == :take
            if _isready(clause[2])
                return (i, _take!(clause[2]))
            end
        else
            error("Invalid select clause: $clause")
        end
    end
    return (0, nothing)
end
function _select_block(clauses)
    winner_ch = Channel{Tuple{Int, Any}}(1)
    tasks = Array{Task}(undef, length(clauses))
    maintask = current_task()
    for (i, clause) in enumerate(clauses)
        tasks[i] = Threads.@spawn begin
            try
                try
                    if clause[1] == :put
                        wait_put(clause[2], tasks)
                    elseif clause[1] ==  :take
                        wait_select(clause[2], maintask, tasks)
                    end
                catch err
                    if isa(err, SelectInterrupt)
                        yieldto(err.parent)
                        return
                    else
                        rethrow()
                    end
                end
                select_kill_rivals(tasks, i)
                if clause[1] == :put
                    ret = put!(clause[2], clause[3])
                elseif clause[1] == :take
                    ret = _take!(clause[2])
                end
                put!(winner_ch, (i, ret))
            catch err
                Base.throwto(maintask, err)
            end
        end
    end
    take!(winner_ch)
end
"""
`select(clauses[, block=true]) -> (clause_index, clause_value)`

Functional form of the `@select` macro, intended to be used when the set of clauses is dynamic. In general, this method will be less performant than the macro variant.

Clauses are specified as an array of tuples. Each tuple is expected to have 2 or 3 elements, as follows:

1) The clause type (`:take` or `:put`)
2) The waitable object
3) If the clause type is `:put`, the value to insert into the object.

If `block` is `true` (the default), wait for at least one clause to be satisfied and return a tuple whose first elmement is the index of the clause which unblocked first and whose whose second element is the value of the clause (see the manual on `select` for the meaning of clause value).

Otherwise, an arbitrary available clause will be executed, or a return value of `(0, nothing)` will be returned  immediately if no clause is available.
"""
function select(clauses, block=true)
    if block
        _select_block(clauses)
    else
        _select_nonblock(clauses)
    end
end
# package code goes here

end # module


text:
using Makie, GeometryTypes, Colors
using AbstractPlotting: slider!, playbutton
using Observables
using AbstractPlotting: textslider
using FileIO, NIfTI

cd(@__DIR__)
#using Pkg
#Pkg.pkg"add NRRD"
# p = download("http://www.bic.mni.mcgill.ca/~vfonov/icbm/2009/mni_icbm152_nlin_asym_09a_nifti.zip")
p = joinpath(homedir(), "Desktop", "brain.nii.gz")
brain_data = Array{Float32}(niread(p))

r = range(-1, stop = 1, length = size(brain_data, 1))

scene3d = Scene(show_axis = false)

contour!(
    scene3d,
    r, r, r, brain_data,
    alpha = 0.1, colorrange = (1, 93), transparency = true, levels = 5
)
c = scene3d[end]
volume = c[4]
planes = (:xy, :xz, :yz)
sliders = ntuple(3) do i
    s, idx = textslider(1:size(volume[], i), ("x", "y", "z")[i], start = size(volume[], i) ÷ 2)
    plane = planes[i]
    indices = ntuple(3) do j
        planes[j] == plane ? 1 : (:)
    end
    hmap = heatmap!(
        scene3d, r, r, volume[][indices...],
        colorrange = (1, 93),
        interpolate = true
    )[end]
    onany(idx, volume) do idx, vol
        transform!(hmap, (plane, r[idx]))
        indices = ntuple(3) do j
            planes[j] == plane ? idx : (:)
        end
        if checkbounds(Bool, vol, indices...)
            hmap[3][] = view(vol, indices...)
        end
    end
    s
end

b2 = Scene()
button!(b2, "3d/2d"; dimensions = (60, 40)) do clicks
    if iseven(clicks)
        cam3d!(scene3d)
    else
        cam2d!(scene3d)
        update_cam!(scene3d, FRect(-1.25, -1.25, 2.5, 2.5))
    end
end

reverse!(scene3d.plots)
display(hbox(
    scene3d,
    vbox(sliders..., b2)
))


text:
# Julia wrapper for header: /usr/local/cuda/include/cusparse.h

#helper functions
function cusparseCreate()
  handle = Ref{cusparseHandle_t}()
  @check ccall( (:cusparseCreate, libcusparse), cusparseStatus_t, (Ptr{cusparseHandle_t},), handle)
  handle[]
end
function cusparseDestroy(handle)
  @check ccall( (:cusparseDestroy, libcusparse), cusparseStatus_t, (cusparseHandle_t,), handle)
end
function cusparseGetVersion(handle, version)
  @check ccall( (:cusparseGetVersion, libcusparse), cusparseStatus_t, (cusparseHandle_t, Ptr{Cint}), handle, version)
end
function cusparseSetStream(handle, streamId)
  @check ccall( (:cusparseSetStream, libcusparse), cusparseStatus_t, (cusparseHandle_t, CuStream_t), handle, streamId)
end
function cusparseGetStream(handle, streamId)
  @check ccall( (:cusparseGetStream, libcusparse), cusparseStatus_t, (cusparseHandle_t, Ptr{CuStream_t}), handle, streamId)
end
function cusparseGetPointerMode(handle, mode)
  @check ccall( (:cusparseGetPointerMode, libcusparse), cusparseStatus_t, (cusparseHandle_t, Ptr{cusparsePointerMode_t}), handle, mode)
end
function cusparseSetPointerMode(handle, mode)
  @check ccall( (:cusparseSetPointerMode, libcusparse), cusparseStatus_t, (cusparseHandle_t, cusparsePointerMode_t), handle, mode)
end
function cusparseCreateHybMat(hybA)
  @check ccall( (:cusparseCreateHybMat, libcusparse), cusparseStatus_t, (Ptr{cusparseHybMat_t},), hybA)
end
function cusparseDestroyHybMat(hybA)
  @check ccall( (:cusparseDestroyHybMat, libcusparse), cusparseStatus_t, (cusparseHybMat_t,), hybA)
end
function cusparseCreateSolveAnalysisInfo(info)
  @check ccall( (:cusparseCreateSolveAnalysisInfo, libcusparse), cusparseStatus_t, (Ptr{cusparseSolveAnalysisInfo_t},), info)
end
function cusparseDestroySolveAnalysisInfo(info)
  @check ccall( (:cusparseDestroySolveAnalysisInfo, libcusparse), cusparseStatus_t, (cusparseSolveAnalysisInfo_t,), info)
end
function cusparseCreateBsrsm2Info(info)
  @check ccall( (:cusparseCreateBsrsm2Info, libcusparse), cusparseStatus_t, (Ptr{bsrsm2Info_t},), info)
end
function cusparseDestroyBsrsm2Info(info)
  @check ccall( (:cusparseDestroyBsrsm2Info, libcusparse), cusparseStatus_t, (bsrsm2Info_t,), info)
end
function cusparseCreateBsrsv2Info(info)
  @check ccall( (:cusparseCreateBsrsv2Info, libcusparse), cusparseStatus_t, (Ptr{bsrsv2Info_t},), info)
end
function cusparseDestroyBsrsv2Info(info)
  @check ccall( (:cusparseDestroyBsrsv2Info, libcusparse), cusparseStatus_t, (bsrsv2Info_t,), info)
end
function cusparseCreateCsrsv2Info(info)
  @check ccall( (:cusparseCreateCsrsv2Info, libcusparse), cusparseStatus_t, (Ptr{csrsv2Info_t},), info)
end
function cusparseDestroyCsrsv2Info(info)
  @check ccall( (:cusparseDestroyCsrsv2Info, libcusparse), cusparseStatus_t, (csrsv2Info_t,), info)
end
function cusparseCreateCsric02Info(info)
  @check ccall( (:cusparseCreateCsric02Info, libcusparse), cusparseStatus_t, (Ptr{csric02Info_t},), info)
end
function cusparseDestroyCsric02Info(info)
  @check ccall( (:cusparseDestroyCsric02Info, libcusparse), cusparseStatus_t, (csric02Info_t,), info)
end
function cusparseCreateCsrilu02Info(info)
  @check ccall( (:cusparseCreateCsrilu02Info, libcusparse), cusparseStatus_t, (Ptr{csrilu02Info_t},), info)
end
function cusparseDestroyCsrilu02Info(info)
  @check ccall( (:cusparseDestroyCsrilu02Info, libcusparse), cusparseStatus_t, (csrilu02Info_t,), info)
end
function cusparseCreateBsric02Info(info)
  @check ccall( (:cusparseCreateBsric02Info, libcusparse), cusparseStatus_t, (Ptr{bsric02Info_t},), info)
end
function cusparseDestroyBsric02Info(info)
  @check ccall( (:cusparseDestroyBsric02Info, libcusparse), cusparseStatus_t, (bsric02Info_t,), info)
end
function cusparseCreateBsrilu02Info(info)
  @check ccall( (:cusparseCreateBsrilu02Info, libcusparse), cusparseStatus_t, (Ptr{bsrilu02Info_t},), info)
end
function cusparseDestroyBsrilu02Info(info)
  @check ccall( (:cusparseDestroyBsrilu02Info, libcusparse), cusparseStatus_t, (bsrilu02Info_t,), info)
end

function cusparseGetProperty(property::CUDAapi.libraryPropertyType)
  value_ref = Ref{Cint}()
  @check ccall((:cusparseGetProperty, libcusparse),
               cusparseStatus_t,
               (Cint, Ptr{Cint}),
               property, value_ref)
  value_ref[]
end


text:
using ITensors
using ITensorsGrad
using OptimKit
using Random
using Zygote

using ITensors: data

Random.seed!(1243)

N = 20
χ = 20

s = siteinds("S=1/2", N; conserve_qns=true)
os = OpSum()
for n in 1:(N - 1)
  os .+= 0.5, "S+", n, "S-", n + 1
  os .+= 0.5, "S-", n, "S+", n + 1
  os .+= "Sz", n, "Sz", n + 1
end
Hmpo = MPO(os, s)
ψ₀mps = randomMPS(s, n -> isodd(n) ? "↑" : "↓"; linkdims=χ)

H = data(Hmpo)
ψ₀ = data(ψ₀mps)

# The Rayleigh quotient to minimize
function E(H::Vector{ITensor}, ψ::Vector{ITensor})
  N = length(ψ)
  ψdag = dag.(addtags.(ψ, "bra"; tags="Link"))
  ψ′dag = prime.(ψdag)
  e = ITensor(1.0)
  for n in 1:N
    e = e * ψ′dag[n] * H[n] * ψ[n]
  end
  norm = ITensor(1.0)
  for n in 1:N
    norm = norm * ψdag[n] * ψ[n]
  end
  return e[] / norm[]
end

@show sum(flux, ψ₀)
@show E(H, ψ₀)

E(ψ) = E(H, ψ)
∇E(ψ) = E'(ψ)
fg(ψ) = (E(ψ), ∇E(ψ))

linesearch = HagerZhangLineSearch(; c₁=0.1, c₂=0.9, ϵ=1e-6, θ=1 / 2, γ=2 / 3, ρ=5.0)
algorithm = LBFGS(5; maxiter=100, gradtol=1e-8, linesearch=linesearch, verbosity=2)
ψ, fψ, gψ, numfg, normgradhistory = optimize(fg, ψ₀, algorithm)

sweeps = Sweeps(5)
setmaxdim!(sweeps, χ)
fψ_dmrg, ψ_dmrg = dmrg(Hmpo, ψ₀mps, sweeps)

time_Eψ = @elapsed E(ψ)
time_∇Eψ = @elapsed E'(ψ)
@show time_∇Eψ / time_Eψ;
@show fψ_dmrg
@show E(H, ψ)

max_stars_repo_name: mtfishman/ITensorsGrad

max_issues_repo_name: mtfishman/ITensorsGrad

max_forks_repo_name: mtfishman/ITensorsGrad


text:
# adapted from CuArrays (https://github.com/JuliaGPU/CuArrays.jl) wrap.jl
# https://github.com/JuliaGPU/CuArrays.jl/blob/master/res/wrap/wrap.jl

using Clang, Libdl

function wrap(name, headers...; library="lib$name", defines=[])
    include_dirs = dirname.(headers)

    clang_args = String[]
    append!(clang_args, map(dir->"-I$dir", include_dirs))
    unique!(clang_args)
    for define in defines
        if isa(define, Pair)
            append!(clang_args, ["-D", "$(first(define))=$(last(define))"])
        else
            append!(clang_args, ["-D", "$define"])
        end
    end

    output_file = "lib$(name).jl"
    common_file = "lib$(name)_common.jl"

    context = init(;
                   headers = [headers...],
                   output_file = output_file,
                   common_file = common_file,
                   clang_includes = [include_dirs..., CLANG_INCLUDE],
                   clang_args = clang_args,
                   header_wrapped = (root, current)->root == current,
                   header_library = x->library,
                   clang_diagnostics = true,
                  )
    run(context)

    return output_file, common_file
end

#
# Fixing-up
#

using CSTParser, Tokenize

## pass infrastructure

struct Edit{T}
    loc::T
    text::String
end

function pass(x, state, f = (x, state)->nothing)
    f(x, state)
    if x.args isa Vector
        for a in x.args
            pass(a, state, f)
        end
    else
        state.offset += x.fullspan
    end
    state
end

function apply(text, edit::Edit{Int})
    string(text[1:edit.loc], edit.text, text[nextind(text, edit.loc):end])
end
function apply(text, edit::Edit{UnitRange{Int}})
    # println("Rewriting '$(text[edit.loc])' to '$(edit.text)'")
    string(text[1:prevind(text, first(edit.loc))], edit.text, text[nextind(text, last(edit.loc)):end])
end

## rewrite passes

mutable struct State
    offset::Int
    edits::Vector{Edit}
end

# insert `@check` before each `ccall` when it returns a checked type
const checked_types = [
    "rocblas_status",
    "rocfft_status",
    "rocsparse_status",
    "rocsolver_status",
    "rocrand_status",
    "miopenStatus_t"
]
function insert_check(x, state)
    if x isa CSTParser.EXPR && x.typ == CSTParser.Call && x.args[1].val == "ccall"
        # get the ccall return type
        rv = x.args[5]

        if rv.val in checked_types
            push!(state.edits, Edit(state.offset, "@check "))
        end
    end
end

## indenting passes

mutable struct IndentState
    offset::Int
    lines
    edits::Vector{Edit}
end

function get_lines(text)
    lines = Tuple{Int,Int}[]
    pt = Tokens.EMPTY_TOKEN(Tokens.Token)
    for t in CSTParser.Tokenize.tokenize(text)
        if pt.endpos[1] != t.startpos[1]
            if t.kind == Tokens.WHITESPACE
                nl = findfirst("\n", t.val) != nothing
                if !nl
                    push!(lines, (length(t.val), 0))
                else
                end
            else
                push!(lines, (0, 0))
            end
        elseif t.startpos[1] != t.endpos[1] && t.kind == Tokens.TRIPLE_STRING
            nls = findall(x->x == '\n', t.val)
            for nl in nls
                push!(lines, (t.startpos[2] - 1, nl + t.startbyte))
            end
        elseif t.startpos[1] != t.endpos[1] && t.kind == Tokens.WHITESPACE
            push!(lines, (t.endpos[2], t.endbyte - t.endpos[2] + 1))
        end
        pt = t
    end
    lines
end

function wrap_at_comma(x, state, indent, offset, column)
    comma = nothing
    for y in x
        if column + y.fullspan > 92 && comma !== nothing
            column = indent
            push!(state.edits, Edit(comma, ",\n" * " "^column))
            column += offset - comma[1] - 1 # other stuff might have snuck between the comma and the current expr
            comma = nothing
        elseif y.typ == CSTParser.PUNCTUATION && y.kind == Tokens.COMMA
            comma = (offset+1):(offset+y.fullspan)
        end
        offset += y.fullspan
        column += y.fullspan
    end
end

function indent_ccall(x, state)
    if x isa CSTParser.EXPR && x.typ == CSTParser.Call && x.args[1].val == "ccall"
        # figure out how much to indent by looking at where the expr starts
        line = findlast(y -> state.offset >= y[2], state.lines) # index, not the actual number
        line_indent, line_offset = state.lines[line]
        expr_indent = state.offset - line_offset
        indent = expr_indent + sum(x->x.fullspan, x.args[1:2])

        if length(x.args[7]) > 2    # non-empty tuple type
            # break before the tuple type
            offset = state.offset + sum(x->x.fullspan, x.args[1:6])
            push!(state.edits, Edit(offset:offset, "\n" * " "^indent))

            # wrap tuple type
            wrap_at_comma(x.args[7], state, indent+1, offset, indent+1)
        end

        if length(x.args) > 9
            # break before the arguments
            offset = state.offset + sum(x->x.fullspan, x.args[1:8])
            push!(state.edits, Edit(offset:offset, "\n" * " "^indent))

            # wrap arguments
            wrap_at_comma(x.args[9:end], state, indent, offset, indent)
        end
    end
end

function indent_definition(x, state)
    if x isa CSTParser.EXPR && x.typ == CSTParser.FunctionDef
        # figure out how much to indent by looking at where the expr starts
        line = findlast(y -> state.offset >= y[2], state.lines) # index, not the actual number
        line_indent, line_offset = state.lines[line]
        expr_indent = state.offset - line_offset
        indent = expr_indent + x.args[1].fullspan + sum(x->x.fullspan, x.args[2].args[1:2])

        if length(x.args[2]) > 2    # non-empty args
            offset = state.offset + x.args[1].fullspan + sum(x->x.fullspan, x.args[2].args[1:2])
            wrap_at_comma(x.args[2].args[3:end-1], state, indent, offset, indent)
        end
    end
end

function process(name, headers...; kwargs...)
    output_file, common_file = wrap(name, headers...; kwargs...)
    text = read(output_file, String)

    ## rewriting passes

    state = State(0, Edit[])
    ast = CSTParser.parse(text, true)

    state.offset = 0
    pass(ast, state, insert_check)

    # apply
    state.offset = 0
    sort!(state.edits, lt = (a,b) -> first(a.loc) < first(b.loc), rev = true)
    for i = 1:length(state.edits)
        text = apply(text, state.edits[i])
    end

    ## indentation passes

    lines = get_lines(text)
    state = IndentState(0, lines, Edit[])
    ast = CSTParser.parse(text, true)

    state.offset = 0
    pass(ast, state, indent_definition)

    state.offset = 0
    pass(ast, state, indent_ccall)

    # apply
    state.offset = 0
    sort!(state.edits, lt = (a,b) -> first(a.loc) < first(b.loc), rev = true)
    for i = 1:length(state.edits)
        text = apply(text, state.edits[i])
    end

    ## manual patches

    write(output_file, text)

    patchdir = joinpath(@__DIR__, "patches")
    for entry in readdir(patchdir)
        if startswith(entry, name) && endswith(entry, ".patch")
            path = joinpath(patchdir, entry)
            run(`patch -p1 -i $path`)
        end
    end

    ## move to destination

    for src in (output_file, common_file)
        dst = joinpath(dirname(@__DIR__), "src", name[4:end], src)
        cp(src, dst; force=true)
    end

    return
end

# from https://github.com/JuliaGPU/CUDAapi.jl/blob/master/src/discovery.jl
function resolve(path)
    if islink(path)
        dir = dirname(path)
        resolve(joinpath(dir, readlink(path)))
    else
        abspath(path)
    end
end

function find_include_dir(name)
    if Libdl.find_library("lib$(name)") != ""
        path = resolve(Libdl.dlpath("lib$(name)"))
        # installation dir should contain lib folder with library and include folder
        i =  Base.findall(lib, path)[end-1][1] - 1
        include_dir = joinpath(path[1:i], "include")
    else
        include_dir = ""
    end
    return include_dir
end

function main()
    rocm_include = isdir("/opt/rocm/include") ? "/opt/rocm/include" : ""
    function process_if_existing(name, headers...; kwargs...)
        env_var = uppercase(name) * "_INCLUDE"
        include = get(ENV, env_var, rocm_include)
        include = include != "" ? include : find_include_dir(name)
        for header in headers
            if !isfile(joinpath(include, header))
                @error "Cannot wrap $name: could not find $header in $include (override the include directory with the $env_var environment variable)"
                return
            end
        end
        _headers = map(header->resolve(joinpath(include, header)), headers)
        # hack since Clang can not use cstddef header
        headers = Vector{String}()
        for h in _headers
            filename = h[findlast("/", h)[1]+1:end]
            cp(h, joinpath(@__DIR__, filename))
            run(`sed -i 's/cstddef/stddef.h/' $filename`)
            push!(headers, resolve(filename))
        end
        process(name, headers...; kwargs...)
        rm.(headers)
        return
    end

    process_if_existing("rocfft", "rocfft-version.h", "rocfft.h", "rocfft-export.h")
end

if abspath(PROGRAM_FILE) == @__FILE__
    main()
end


text:
@testset "trace" begin
    @testset "basic" begin
        @test isnothing(current_span())
        @test isnothing(span_context())

        @testset "TraceState" begin
            ts = TraceState(
                "1a-2f@foo" => "bar1",
                "1a-_*/2b@foo" => "bar2",
                "foo" => "bar3",
                "foo-_*/bar" => "bar4",
                "^foo-_*/bar" => "bar5",
                "hello" => "世界！",
            )

            @test string(ts) == "1a-2f@foo=bar1,1a-_*/2b@foo=bar2,foo=bar3,foo-_*/bar=bar4"
            @test length(ts) == 4
            @test haskey(ts, :foo)
            @test !haskey(ts, "^foo-_*/bar")
            @test ts[:foo] == "bar3"
        end
    end

    @testset "TracerProvider" begin
        tracer = Tracer()
        s = Span("test", tracer)
        @test span_context(s) === INVALID_SPAN_CONTEXT
        @test is_recording(s) == false

        @testset "non recording behaviors" begin
            s["foo"] = "bar"
            push!(s, API.Event(; name = "test"))
            push!(s, Link(INVALID_SPAN_CONTEXT, StaticAttrs()))
            set_status!(s, SPAN_STATUS_OK)

            @test haskey(s, "foo") == false
            @test length(s.events) == 0
            @test length(s.links) == 0
            @test s.status[].code == SPAN_STATUS_UNSET
        end

        s.end_time[] = nothing # !!! for test only

        @test_throws ErrorException with_span(s) do
            s["foo"] = "bar"
            push!(s, API.Event(; name = "test"))
            push!(s, Link(INVALID_SPAN_CONTEXT, StaticAttrs()))
            throw(ErrorException("!!!"))
        end

        @test is_recording(s) == false
        @test s["foo"] == "bar"
        @test length(s.events) == 2  # note that the ErrorException is also recorded
        @test length(s.links) == 1
        @test s.status[].code == SPAN_STATUS_ERROR

        with_span(Span("foo", tracer)) do
            with_span(Span("bar", tracer)) do
                with_span(Span("baz", tracer)) do
                    println("hello")
                end
                with_span(Span("baz", tracer)) do
                    println("world!")
                end
            end
        end
    end
end


text:
#   This file is part of Reduce.jl. It is licensed under the MIT license
#   Copyright (C) 2017 Michael Reed

ispsl() = haskey(ENV,"PSL")
red() = ispsl() ? `redpsl` : `redcsl -w`
red(path) = ispsl() ? `$path/redpsl` : `$path/redcsl -w`
redsl() = ispsl() ? "psl" : "csl"

ρ = 7

date = Dict(
  0 =>  "2017-05-16",
  1 =>  "2018-01-17",
  2 =>  "2018-04-06",
  3 =>  "2019-04-13",
  4 =>  "2019-09-15",
  5 =>  "2020-03-01",
  6 =>  "2020-06-24",
  7 =>  "2020-10-07")

rsvn = Dict(
  0 =>  "4052",
  1 =>  "4372",
  2 =>  "4567",
  3 =>  "4961",
  4 =>  "5129",
  5 =>  "5286",
  6 =>  "5377",
  7 =>  "5424")

function redsys(dir)
    if !(Sys.iswindows())
        if Sys.islinux()
            red(joinpath(dir,"usr","bin"))
        elseif Sys.isapple()
            red(joinpath(dir,ispsl() ? "psl" : "csl"))
        else
            rsl = red(joinpath(dir,"Reduce-svn$(rsvn[ρ])-src","bin"))
        end
    else
        if ispsl()
            `"$(dir)\psl\bpsl.exe" -td 16000000 -f "$(dir)\red\reduce.img"`
        else
            `"$(dir)\reduce.exe" --nogui`
        end
    end
end

function _spawn(cmd, input=devnull, output=devnull)
    run(pipeline(cmd, stdin=input, stdout=output, stderr=stderr), wait=false)
end


text:
function cross2d(x,y)
    return x[1]*y[2] - x[2]*y[1]
end

function intersect_line(ls1::NTuple{2,T},ls2::NTuple{2,S}) where {T<:AbstractVector,S<:AbstractVector}
    r = ls1[2] - ls1[1]
    s = ls2[2] - ls2[1]
    qmp = ls2[1]-ls1[1]
    d = cross2d(r,s)
    un = cross2d(qmp,r)
    d == 0.0 && (return un == 0.0)
    return (0 < un/d < 1) && (0 < cross2d(qmp, s)/d < 1)
end

function get_angle(p_pro,p_cur,p_pre)
    r = p_pro - p_cur
    p = p_pre - p_cur
    t = atan(cross2d(p,r),dot(p,r))
    return t > 0 ? t : t + 2pi
end

function intersect_hull(ls, hull)
    v = hull.vertices
    np = length(v)
    if np < 3
        return false
    end
    for i = 1:(np-2)
        if intersect_line(ls,(v[i+1],v[i]))
            return true
        end
    end
    return false
end

# In Julia v0.7 and above, indmax cannot be used with generators.
# In a future version, we should be able to replace this function
# with Base.argmax(f, itr)
function indmax_f(f::Function, itr::AbstractVector)
    imax = firstindex(itr)
    vmax = f(itr[imax])
    for i in (firstindex(itr)+1):lastindex(itr)
        v = f(itr[i])
        if v > vmax
            vmax = v
            imax = i
        end
    end
    imax
end

function concave_hull(tree::KDTree, k::Int)
    npoints = length(tree.data)
    k = clamp(k, 3, npoints-1)
    hull = Hull(eltype(tree.data), k)

    #Start at point with largest x value
    i0 = indmax_f(v -> v[1], tree.data)

    ishull = zeros(Bool,npoints)
    ishull[i0] = true

    #Initiate hull
    p0 = copy(tree.data[i0])
    push!(hull.vertices,p0)

    #Allocate point arrays
    p_cur = copy(p0)
    p_pre = p_cur + [1, 0]
    p_pro = copy(p0)
    d01 = zero(eltype(p0))

    nstep = 1
    npick = 1
    while npick < npoints
        if nstep == 3
            npick = npick - 1
            ishull[i0] = false
        end
        kk = clamp(k, 3 , npoints-npick)
        kind, kdist = knn(tree, p_cur, kk, true, i -> ishull[i])
        angles = collect(get_angle(tree.data[i], p_cur, p_pre) for i in kind)
        w = sortperm(angles,rev=true)
        inter = true
        for i in kind[w]
            #p_pre(0) ---------- d02 --------- p_pro(2)
            #p_pre(0) --d01-- p_cur(1) --d12-- p_pro(2)
            p_pro = tree.data[i]
            d12 = norm(p_pro - p_cur)
            d02 = norm(p_pro - p_pre)
            if ~intersect_hull((p_pro,p_cur),hull) || (d01 + d12) == d02
                p_pre = p_cur
                d01 = d12
                p_cur = tree.data[i]
                ishull[i] = true
                inter = false
                break
            end
        end

        if inter
            if k+1 >= npoints
                @info("Unable to construct concave hull")
                return hull
            end
            return concave_hull(tree, k+1)
        end

        if p_cur != p0
            push!(hull.vertices, copy(p_cur))
            npick = npick + 1
        else
            break
        end
        nstep = nstep + 1
    end

    if any(~in_hull(p,hull) for p in tree.data[.~(ishull)])
        if k+1 >= npoints
            @info("Unable to construct concave hull")
            return hull
        end
        return concave_hull(tree, k+1)
    end

    return Hull(hull.vertices, hull.k, true)
end

function concave_hull(points::Vector, k::Int = 3)
    p = unique(points)
    npoints = length(p)
    npoints < 3 && throw(ArgumentError("Number of unique points should be greater than 2"))
    npoints == 3 && return Hull(points,k, true)
    tree = KDTree(hcat(points...),reorder=false)
    return concave_hull(tree, k)
end


text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License. You
# may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable
# law or agreed to in writing, software distributed under the License is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the License for the specific language
# governing permissions and limitations under the License.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

"""
Trajplot jsons into two folders depending on an indicator vector.
"""
module SortPlotsToFolders

export sort_plots_to_folders
export script_APL042017_10K_cassim

using DataFrames
using RLESUtils, PGFPlotUtils
using RLESCAS
include_visualize()

const APL042017_10K_CASSIM_META = Pkg.dir("Datasets/data/APL042017_10K_cassim/_META.csv.gz")

function sort_plots_to_folders{T<:AbstractString}(files::Vector{T},
    indicator::Vector{Bool}, true_folder::AbstractString,
    false_folder::AbstractString; verbose::Bool=true,
    format::Symbol=:TEXPDF)

    @assert length(files) == length(indicator)
    mkpath(true_folder)
    mkpath(false_folder)
    for i = 1:length(files)
        verbose && println("file $i of $(length(files))")
        f = files[i]
        b = indicator[i]
        folder = b ? true_folder : false_folder
        d = trajLoad(f)
        outfileroot = joinpath(folder, split(basename(f), ".")[1])
        trajPlot(outfileroot, d; format=format)
    end
end

function script_APL042017_10K_cassim(basefolder::AbstractString="./";
    format::AbstractString="PDF")

    D = readtable(APL042017_10K_CASSIM_META)
    encounter_ids = convert(Array, D[:encounter_id])
    files = map(encounter_ids) do id
        f = "trajSaveMCTS_ACASX_GM_" * lpad(id, 5, 0) * ".json.gz"
        joinpath(basefolder, f)
    end
    indicator = convert(Array, D[:nmac])
    true_folder = joinpath(basefolder, "nmac")
    false_folder = joinpath(basefolder, "non-nmac")
    sort_plots_to_folders(files, indicator, true_folder, false_folder;
        format=format)
end

end #module

max_stars_repo_name: sisl/RLESCAS

max_issues_repo_name: sisl/RLESCAS

max_forks_repo_name: sisl/RLESCAS


text:
@testset "CP.jl" begin

    include("variables/variables.jl")
    include("core/model.jl")
    include("constraints/constraints.jl")

    include("core/fixPoint.jl")
    include("core/cpModification.jl")

    include("variableselection/variableselection.jl")

    include("valueselection/valueselection.jl")

    include("core/search/search.jl")
    include("core/solver.jl")

    # @testset "solve()" begin
     #    @test (@test_logs (:info, "Solved !") SeaPearl.solve()) == nothing
    # end
end