text:
"""
YarnClientProtocol: Hadoop RPC client for application client to Yarn resource manager protocol.
"""
const YarnClientProtocol = HadoopRpcProtocol{ApplicationClientProtocolServiceBlockingStub}

for fn in (:getClusterMetrics, :getClusterNodes, :getNewApplication, :submitApplication, :forceKillApplication, :getApplicationReport, :getApplicationAttempts)
    @eval begin
       (hadoop.yarn.$fn)(p::YarnClientProtocol, inp) = (hadoop.yarn.$fn)(p.stub, p.controller, inp)
    end
end

"""
YarnAMRMProtocol: Hadoop RPC client for Yarn application master to resource manager protocol.
"""
const YarnAMRMProtocol = HadoopRpcProtocol{ApplicationMasterProtocolServiceBlockingStub}

for fn in (:registerApplicationMaster, :finishApplicationMaster, :allocate)
    @eval begin
        (hadoop.yarn.$fn)(p::YarnAMRMProtocol, inp) = (hadoop.yarn.$fn)(p.stub, p.controller, inp)
    end
end

"""
YarnAMNMProtocol: Hadoop RPC client for Yarn application master to node manager protocol.
"""
const YarnAMNMProtocol = HadoopRpcProtocol{ContainerManagementProtocolServiceBlockingStub}

for fn in (:startContainers, :stopContainers, :getContainerStatuses)
    @eval begin
        (hadoop.yarn.$fn)(p::YarnAMNMProtocol, inp) = (hadoop.yarn.$fn)(p.stub, p.controller, inp)
    end
end

"""
# YarnException
Thrown by Yarn APIs.
"""
mutable struct YarnException <: Exception
    message::String
end

function show(io::IO, serex::SerializedExceptionProto)
    print(io, "Exception: ")
    isfilled(serex, :class_name) && print(io, serex.class_name, ": ")
    isfilled(serex, :message) && print(io, serex.message)
    isfilled(serex, :trace) && print(io, '\n', serex.trace, '\n')

    if isfilled(serex, :cause)
        println(io, "Caused by:")
        show(io, serex.cause)
    end
    nothing
end

function YarnException(cex::ContainerExceptionMapProto)
    cont_id = cex.ccontainer_id
    app_id = cont_id.app_id
    atmpt_id = cont_id.app_attempt_id
    iob = IOBuffer()
    println(iob, "Error in container $(cont_id.id) of app $(app_id.id) attempt $(atmpt_id.attemptId):")
    show(iob, cex.exception)
    YarnException(String(take!(iob)))
end

"""
YarnNode represents a node manager in the yarn cluster and its
communication address, resource state and run state.
"""
mutable struct YarnNode
    host::AbstractString
    port::Int32
    rack::AbstractString
    ncontainers::Int32
    mem::Int32
    cores::Int32
    memused::Int32
    coresused::Int32
    state::Int32
    isrunning::Bool
    report::NodeReportProto

    function YarnNode(node::NodeReportProto)
        host = node.nodeId.host
        port = node.nodeId.port
        rack = node.rackName
        ncontainers = node.numContainers

        state = node.node_state
        isrunning = (state == NodeStateProto.NS_RUNNING)

        if isrunning
            mem = node.capability.memory
            cores = node.capability.virtual_cores

            memused = node.used.memory
            coresused = node.used.virtual_cores
        else
            mem = cores = memused = coresused = 0
        end

        new(host, port, rack, ncontainers, mem, cores, memused, coresused, state, isrunning, node)
    end
end

"""
NODE_STATES: enum value to state map. Used for converting state for display.
"""
const NODE_STATES = [:new, :running, :unhealthy, :decommissioned, :lost, :rebooted]

function show(io::IO, node::YarnNode)
    print(io, "YarnNode: $(node.rack)/$(node.host):$(node.port) $(NODE_STATES[node.state])")
    print(io, node.isrunning ? ", Used mem: $(node.memused)/$(node.mem), cores: $(node.coresused)/$(node.cores)" : "")
    nothing
end

"""
YarnNodes holds node information as visible to the application master.
It also caches connection to node masters. Connection are reused if they are required before a set keepalivesecs time.
"""
mutable struct YarnNodes
    count::Int32
    status::Dict{NodeIdProto,YarnNode}
    conn::Dict{NodeIdProto,Tuple}
    ugi::UserGroupInformation
    keepalivesecs::UInt64
    lck::ReentrantLock

    function YarnNodes(ugi::UserGroupInformation)
        new(0, Dict{NodeIdProto,YarnNode}(), Dict{NodeIdProto,Tuple}(), ugi, YARN_NM_CONN_KEEPALIVE_SECS, ReentrantLock())
    end
end

function show(io::IO, nodes::YarnNodes)
    println(io, "YarnNodes: $(nodes.count) (connected to $(length(nodes.conn)))")
    for n in values(nodes.status)
        show(io, n)
        println(io, "")
    end
    nothing
end

function update(nodes::YarnNodes, arp::AllocateResponseProto)
    isfilled(arp, :num_cluster_nodes) && (nodes.count = arp.num_cluster_nodes)
    if isfilled(arp, :updated_nodes)
        for nrep in arp.updated_nodes
            nodes.status[nrep.nodeId] = YarnNode(nrep)
        end
    end
    nothing
end

function update(nodes::YarnNodes, gcnrp::GetClusterNodesResponseProto)
    isfilled(gcnrp, :nodeReports) || return

    nlist = gcnrp.nodeReports
    nodes.count = length(nlist)
    for nrep in nlist
        nodes.status[nrep.nodeId] = YarnNode(nrep)
    end
    nothing
end

function _new_or_existing_conn(nodes::YarnNodes, nodeid::NodeIdProto)
    lock(nodes.lck)
    try
        t = time()
        (conn,lastusetime,lck) = (nodeid in keys(nodes.conn)) ? nodes.conn[nodeid] : (YarnAMNMProtocol(nodeid.host, nodeid.port, nodes.ugi), t, ReentrantLock())
        if t > (lastusetime + nodes.keepalivesecs)
            try
                disconnect(conn.channel)
            catch
                # ignore exception
            finally
                conn = YarnAMNMProtocol(nodeid.host, nodeid.port, nodes.ugi)
                lastusetime = t
            end
        end
        nodes.conn[nodeid] = (conn, lastusetime, lck)
        return (conn, lastusetime, lck)
    catch ex
        rethrow(ex)
    finally
        unlock(nodes.lck)
    end
end

function get_connection(nodes::YarnNodes, nodeid::NodeIdProto)
    (nodeid in keys(nodes.status)) || throw(YarnException("Unknown Yarn node: $(nodeid.host):$(nodeid.port)"))
    node = nodes.status[nodeid]
    node.isrunning || throw(YarnException("Yarn node $(nodeid.host):$(nodeid.port) is not running"))

    conn, lastusetime, lck = _new_or_existing_conn(nodes, nodeid)
    # lock the connection to mark in use
    lock(lck)
    conn
end

function release_connection(nodes::YarnNodes, nodeid::NodeIdProto, conn::YarnAMNMProtocol, reuse::Bool)
    lock(nodes.lck)
    try
        # disconnect if we don't want to reuse the connection
        if !reuse
            try
                disconnect(conn.channel)
            catch
                # ignore exceptions
            end
        end

        if nodeid in keys(nodes.conn)
            # release the lock if it's an existing connection
            (conn_old, lastusetime, lck) = nodes.conn[nodeid]
            (conn === conn_old) && unlock(lck)
            if reuse
                nodes.conn[nodeid] = (conn, time(), lck)
            else
                delete!(nodes.conn, nodeid)
            end
        else
            # add connection to cache for reuse
            if reuse
                nodes.conn[nodeid] = (conn, time(), ReentrantLock())
            end
        end
    catch ex
        rethrow(ex)
    finally
        unlock(nodes.lck)
    end
    nothing
end

"""
RequestPipeline holds entities while they are requested for from yarn resource manager.
Application master thread extracts pending items and requests them from RM, whereupon they are moved on to the requested state.
"""
mutable struct RequestPipeline{T}
    pending::Vector{T}
    requested::Vector{T}

    function RequestPipeline{T}() where T
        new(T[], T[])
    end
end

pending(pipe::RequestPipeline{T}, item::T) where {T} = push!(pipe.pending, item)
function torequest(pipe::RequestPipeline{T}) where T
    ret = pipe.pending
    if !isempty(ret)
        append!(pipe.requested, ret)
        pipe.pending = T[]
    end
    ret
end
haverequests(pipe::RequestPipeline) = !isempty(pipe.pending)

"""
YarnContainers holds all containers related to the application.
It also holds the allocation and release pipelines that are used by application master for requesting actions from resource manager.
Also schedules callbacks as tasks when containers are allocated or terminated.
"""
mutable struct YarnContainers
    containers::Dict{ContainerIdProto,ContainerProto}
    status::Dict{ContainerIdProto,ContainerStatusProto}
    active::Set{ContainerIdProto}
    busy::Set{ContainerIdProto}

    alloc_pipeline::RequestPipeline{ResourceRequestProto}
    release_pipeline::RequestPipeline{ContainerIdProto}
    ndesired::Int

    on_container_alloc::Union{Nothing,Function}
    on_container_finish::Union{Nothing,Function}

    function YarnContainers()
        new(Dict{ContainerIdProto,ContainerProto}(), Dict{ContainerIdProto,ContainerStatusProto}(), Set{ContainerIdProto}(), Set{ContainerIdProto}(),
            RequestPipeline{ResourceRequestProto}(), RequestPipeline{ContainerIdProto}(), 0, nothing, nothing)
    end
end

function show(io::IO, containers::YarnContainers)
    println(io, "Containers: $(length(containers.active))/$(length(containers.containers)) active, $(length(containers.busy)) in use")
    nothing
end

function callback(containers::YarnContainers, on_container_alloc::Union{Nothing,Function}, on_container_finish::Union{Nothing,Function})
    containers.on_container_alloc = on_container_alloc
    containers.on_container_finish = on_container_finish
    nothing
end

function update(containers::YarnContainers, arp::AllocateResponseProto)
    active = containers.active
    busy = containers.busy
    status = containers.status
    contlist = containers.containers
    cballoc = containers.on_container_alloc
    cbfinish = containers.on_container_finish

    if isfilled(arp, :allocated_containers)
        for cont in arp.allocated_containers
            id = cont.id
            contlist[id] = cont
            push!(active, id)
            @debug("calling callback for alloc")
            (cballoc === nothing) || @async cballoc(id)
        end
    end
    if isfilled(arp, :completed_container_statuses)
        @debug("have completed containers")
        for contst in arp.completed_container_statuses
            id = contst.container_id
            @debug("container $id is finished")
            status[id] = contst
            @debug("id in active: $(id in active)")
            (id in active) && pop!(active, id)
            (id in busy) && pop!(busy, id)
            @debug("calling callback for finish")
            (cbfinish === nothing) || @async cbfinish(id)
        end
    end
    nothing
end

# See Yarn issues related to allocating containers at specific locations:
# https://issues.apache.org/jira/browse/YARN-2027
# https://issues.apache.org/jira/browse/YARN-1412
function request_alloc(containers::YarnContainers, numcontainers::Int;
                    mem::Integer=YARN_CONTAINER_MEM_DEFAULT, cpu::Integer=YARN_CONTAINER_CPU_DEFAULT,
                    loc::AbstractString=YARN_CONTAINER_LOCATION_DEFAULT, priority::Integer=YARN_CONTAINER_PRIORITY_DEFAULT)
    prio = protobuild(PriorityProto, Dict(:priority => priority))
    capability = protobuild(ResourceProto, Dict(:memory => mem, :virtual_cores => cpu))
    req = protobuild(ResourceRequestProto, Dict(:priority => prio,
            :resource_name => loc,
            :num_containers => numcontainers,
            :capability => capability))
    pending(containers.alloc_pipeline, req)
    containers.ndesired += numcontainers
    nothing
end

function request_release(containers::YarnContainers, cids::ContainerIdProto...)
    for cid in cids
        pending(containers.release_pipeline, cid)
    end
    containers.ndesired -= length(cids)
    nothing
end

function set_busy(containers::YarnContainers, cids::ContainerIdProto...)
    busy = containers.busy
    for cid in cids
        push!(busy, cid)
    end
    nothing
end

function set_free(containers::YarnContainers, cids::ContainerIdProto...)
    busy = containers.busy
    for cid in cids
        pop!(busy, cid)
    end
    nothing
end

torequest(containers::YarnContainers) = (torequest(containers.alloc_pipeline), torequest(containers.release_pipeline))
haverequests(containers::YarnContainers) = containers.ndesired != length(containers.active)

# TODO: support local resources
# TODO: support tokens
function launchcontext(;cmd::AbstractString="", env::Dict=Dict(), service_data::Dict=Dict())
    clc = ContainerLaunchContextProto()
    if !isempty(cmd)
        set_field!(clc, :command, AbstractString[cmd])
    end
    if !isempty(env)
        envproto = StringStringMapProto[]
        for (n,v) in env
            (isa(n, AbstractString) && isa(v, AbstractString)) || throw(ArgumentError("non string environment variable specified: $(typeof(n)) => $(typeof(v))"))
            push!(envproto, protobuild(StringStringMapProto, Dict(:key => n, :value => v)))
        end
        set_field!(clc, :environment, envproto)
    end
    if !isempty(service_data)
        svcdataproto = StringBytesMapProto[]
        for (n,v) in service_data
            (isa(n, AbstractString) && isa(v, Vector{UInt8})) || throw(ArgumentError("incompatible service data type specified: $(typeof(n)) => $(typeof(v))"))
            push!(svcdataproto, protobuild(StringBytesMapProto, Dict(:key => n, :value => v)))
        end
        set_field!(clc, :service_data, servicedataproto)
    end
    clc
end


text:
using CompScienceMeshes, BEAST
# Γ = readmesh(joinpath(@__DIR__,"sphere2.in"))
Γ = readmesh(joinpath(dirname(pathof(BEAST)),"../examples/sphere2.in"))
# Γ = meshsphere(1.0, 0.4)
@show length(Γ)

X = raviartthomas(Γ)
Y = buffachristiansen(Γ)

Δt, Nt = 0.6, 200
T0 = timebasisshiftedlagrange(Δt, Nt, 0)
T1 = timebasisshiftedlagrange(Δt, Nt, 1)
iT0 = integrate(T0)
iT1 = integrate(T1)
δ = timebasisdelta(Δt, Nt)

# V = X ⊗ T
# iV = X ⊗ iT
# W = Y ⊗ δ

duration = 20 * Δt
delay = 1.5 * duration
amplitude = 1.0
gaussian = creategaussian(duration, delay, amplitude)
int_gaussian = integrate(gaussian)

direction, polarisation = ẑ, x̂
E = BEAST.planewave(polarisation, direction, gaussian, 1.0)
iE = BEAST.planewave(polarisation, direction, int_gaussian, 1.0)
H = direction × E
iH = direction × iE

DL = TDMaxwell3D.doublelayer(speedoflight=1.0)
I = Identity()
N = NCross()
M = 0.5*N⊗I + 1.0*DL
iM = integrate(M)
NI = N ⊗ I
iNI = integrate(NI)

# error()

@hilbertspace k
@hilbertspace j
# mfie = @discretise (0.5(N⊗I) + 1.0DL)[k,j] == -1.0H[k] k∈W j∈V
mfie = @discretise iM[k,j] == -1.0iH[k] k∈(Y⊗δ) j∈(X⊗T0)

xmfie_int = solve(mfie)

Xmfie, Δω, ω0 = fouriertransform(xmfie_int, Δt, 0.0, 2)
ω = collect(ω0 .+ (0:Nt-1)*Δω)
_, i1 = findmin(abs.(ω .- 1.0))

ω1 = ω[i1]
um = Xmfie[:,i1] / fouriertransform(gaussian)(ω1)

ext: jl


text:
module AmorphousBuildingBlocks

include("atomic_information.jl"), ("construct_blocks.jl")

end # module


text:
"""
    blosc2_context

Context opaque object
"""
mutable struct blosc2_context
end

const FiltersTuple = NTuple{BLOSC2_MAX_FILTERS, UInt8}

"""
    blosc2_cparams

The parameters for creating a context for compression purposes.

In parenthesis it is shown the default value used internally when a 0
(zero) in the fields of the struct is passed to a function.

preffiler must match blosc2_prefilter_fn defenition:
typedef int (*blosc2_prefilter_fn)(blosc2_prefilter_params* params);
"""
struct blosc2_cparams

  compcode::UInt8 #!< The compressor codec.
  compcode_meta::UInt8 #!< The metadata for the compressor codec.
  clevel::UInt8 #!< The compression level (5).
  use_dict::Cint #!< Use dicts or not when compressing (only for ZSTD).
  typesize ::Int32 #!< The type size (8).
  nthreads ::Int16 #!< The number of threads to use internally (1).
  blocksize ::Int32 #!< The requested size of the compressed blocks (0; meaning automatic).
  splitmode ::Int32 #!< Whether the blocks should be split or not.
  schunk ::Ref{Cvoid} #!< The associated schunk, if any (NULL).
  filters ::FiltersTuple #!< The (sequence of) filters.
  filters_meta ::FiltersTuple #!< The metadata for filters.
  preffiler ::Ptr{Cvoid} #!< The prefilter function.
  preparams ::Ptr{Cvoid} #!< The prefilter parameters. (blosc2_prefilter_params)
  udbtune ::Ptr{Cvoid}#!< The user-defined BTune parameters. (blosc2_btune)

  function blosc2_cparams(;compcode = BLOSC_BLOSCLZ,
      compcode_meta = 0,
      clevel =5,
      use_dict =0,
      typesize = 8,
      nthreads = 1,
      blocksize =0,
      splitmode = BLOSC_FORWARD_COMPAT_SPLIT,
      schunk = C_NULL,
      filters = FiltersTuple([0, 0, 0, 0, 0, BLOSC_SHUFFLE]),
      filters_meta = FiltersTuple([0, 0, 0, 0, 0, 0])
    )
        new(
            compcode,
            compcode_meta,
            clevel,
            use_dict,
            typesize,
            nthreads,
            blocksize,
            splitmode,
            schunk,
            filters,
            filters_meta,
            C_NULL,
            C_NULL,
            C_NULL
        )
  end
end

"""
    blosc2_dparams

The parameters for creating a context for decompression purposes.
In parenthesis it is shown the default value used internally when a 0
(zero) in the fields of the struct is passed to a function.
"""
struct blosc2_dparams
  nthreads ::Cint #!< The number of threads to use internally (1).
  schunk ::Ptr{Cvoid} #!< The associated schunk, if any (NULL).
  postfilter ::Ptr{Cvoid} #!< The postfilter function.
  postparams ::Ptr{Cvoid} #!< The postfilter parameters.
  function blosc2_dparams(;
    nthreads = 1,
    schunk = C_NULL
  )
    new(nthreads, schunk, C_NULL, C_NULL)
  end
end

"""
    blosc2_create_cctx(cparams::blosc2_cparams = blosc2_cparams())

Create a context for *_ctx() compression functions.

# Returns

A pointer to the new context. NULL is returned if this fails.
"""
function blosc2_create_cctx(cparams::blosc2_cparams = blosc2_cparams())
    return @ccall lib.blosc2_create_cctx(cparams::blosc2_cparams)::Ptr{blosc2_context}
end

"""
    blosc2_create_dctx(dparams::blosc2_dparams = blosc2_dparams())

Create a context for *_ctx() decompression functions.

# Returns

A pointer to the new context. NULL is returned if this fails.
"""
function blosc2_create_dctx(dparams::blosc2_dparams = blosc2_dparams())
    return @ccall lib.blosc2_create_dctx(dparams::blosc2_dparams)::Ptr{blosc2_context}
end

"""
    blosc2_free_ctx(context::Ptr{blosc2_context})

Free the resources associated with a context.
"""
function blosc2_free_ctx(context::Ptr{blosc2_context})
    @ccall lib.blosc2_free_ctx(context::Ptr{blosc2_context})::Cvoid
end

"""
    blosc2_compress_ctx(ctx::Ptr{blosc2_context}, src, srcsize, dest, destsize)

Context interface to Blosc compression.

This does not require a call to blosc_init and can be called from multithreaded applications without the global lock being used, so allowing Blosc be executed simultaneously in those scenarios.

# Parameters

- `context` – A blosc2_context struct with the different compression params.

- `src` – The buffer containing the data to be compressed.

- `srcsize` – The number of bytes to be compressed from the src buffer.

- `dest` – The buffer where the compressed data will be put.

- `destsize` – The size in bytes of the dest buffer.

Returns
The number of bytes compressed. If src buffer cannot be compressed into destsize, the return value is zero and you should discard the contents of the dest buffer. A negative return value means that an internal error happened. It could happen that context is not meant for compression (which is stated in stderr). Otherwise, please report it back together with the buffer data causing this and compression settings.

"""
function blosc2_compress_ctx(ctx::Ptr{blosc2_context}, src, srcsize, dest, destsize)
    @ccall lib.blosc2_compress_ctx(
            ctx::Ptr{blosc2_context}, src::Ptr{Cvoid}, srcsize::Int32,  dest::Ptr{Cvoid}, destsize::Int32
        )::Cint
end

"""
    blosc2_decompress_ctx(ctx::Ptr{blosc2_context}, src, srcsize, dest, destsize)

Context interface to Blosc decompression.

This does not require a call to blosc_init and can be called from multithreaded applications without the global lock being used,
so allowing Blosc be executed simultaneously in those scenarios.

# Returns
The number of bytes copied to dest or a negative value if some error happens.

"""
function blosc2_decompress_ctx(ctx::Ptr{blosc2_context}, src, srcsize, dest, destsize)
    @ccall lib.blosc2_decompress_ctx(
            ctx::Ptr{blosc2_context}, src::Ptr{Cvoid}, srcsize::Int32,  dest::Ptr{Cvoid}, destsize::Int32
        )::Cint
end

"""
    blosc2_getitem_ctx(ctx::Ptr{blosc2_context}, src, srcsize, start, nitems, dest, destsize)

Context interface counterpart for `Lib.blosc_getitem`.

"""
function blosc2_getitem_ctx(ctx::Ptr{blosc2_context}, src, srcsize, start, nitems, dest, destsize)
    @ccall lib.blosc2_getitem_ctx(
            ctx::Ptr{blosc2_context}, src::Ptr{Cvoid}, srcsize::Int32,  start::Cint, nitems::Cint, dest::Ptr{Cvoid}, destsize::Int32
        )::Cint
end

"""
    blosc2_set_maskout(ctx::Ptr{blosc2_context}, maskout, nblocks)

Set a maskout so as to avoid decompressing specified blocks.

# Remark

The maskout is valid for contexts only meant for decompressing a chunk via `blosc2_decompress_ctx`.
Once a call to `blosc2_decompress_ctx` is done, this mask is reset so that next call to `blosc2_decompress_ctx` will decompress the whole chunk.

# Parameters

- `ctx` – The decompression context to update.

- `maskout` – The boolean mask for the blocks where decompression is to be avoided.

- `nblocks` – The number of blocks in maskout above.

Returns
If success, a 0 values is returned. An error is signaled with a negative int.
"""
function blosc2_set_maskout(ctx::Ptr{blosc2_context}, maskout::DenseArray{Bool}, nblocks)
    @ccall lib.blosc2_set_maskout(ctx::Ptr{blosc2_context}, maskout::Ptr{Bool}, nblocks::Cint)::Cint
end

"""
    blosc2_ctx_get_cparams(ctx::Ptr{blosc2_context})

Create a cparams associated to a context.

# Returns

`blosc2_cparams` if succeeds. Else `nothing`
"""
function blosc2_ctx_get_cparams(ctx::Ptr{blosc2_context})
    cparams = Ref{blosc2_cparams}()

    r = @ccall lib.blosc2_ctx_get_cparams(ctx::Ptr{blosc2_context}, cparams::Ptr{blosc2_cparams})::Cint

    (r < 0 || cparams[] == C_NULL) && return nothing
    return cparams[]
end
"""
    blosc2_ctx_get_dparams(ctx::Ptr{blosc2_context})

Create a dparams associated to a context.

# Returns

`blosc2_cdarams` if succeeds. Else `nothing`
"""
function blosc2_ctx_get_dparams(ctx::Ptr{blosc2_context})
    dparams = Ref{blosc2_dparams}()

    r = @ccall lib.blosc2_ctx_get_dparams(ctx::Ptr{blosc2_context}, dparams::Ptr{blosc2_dparams})::Cint

    (r < 0 || cparams[] == C_NULL) && return nothing
    return dparams[]
end


text:
using RANSACBenchmark
using RANSAC
using Test
using LinearAlgebra
using StaticArrays: SVector

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

@testset "validate benchmark clouds" begin
    include("benchmarkcloud1.jl")
end


text:
function consolidate{T}(a::Array{Set{T},1})
    1 < length(a) || return a
    b = copy(a)
    c = Set{T}[]
    while 1 < length(b)
        x = shift!(b)
        cme = true
        for (i, y) in enumerate(b)
            !isempty(intersect(x, y)) || continue
            cme = false
            b[i] = union(x, y)
            break
        end
        !cme || push!(c, x)
    end
    push!(c, b[1])
    return c
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
struct SimulationFixed <: SimulationSpace

end

function arguments_(neighbors::SimulationFixed, space::Space, platform::String)
    declareVar = Nothing
    declareF = Nothing
    args = Nothing
    argsEval = Nothing
    execInit = Nothing
    execInloop = Nothing
    execAfter = Nothing

    error("Not implemented yet.")

    return declareVar, declareF, args, argsEval, execInit, execInloop, execAfter
end

function loop_(neighbors::SimulationFixed, space::Space, code::Expr)

    error("Not implemented yet.")

    return loop
end

max_stars_repo_name: gatocor/embryogenesisJulia

max_issues_repo_name: gatocor/embryogenesisJulia

max_forks_repo_name: gatocor/embryogenesisJulia


text:
n = rand(1:10)
# no null vector for 1-by-1 chebspec.
if n == 1
    n = 2
end
A = matrixdepot("chebspec", n)
# A has null vector ones(n)
@test_approx_eq_eps A*ones(n) zeros(n) 1e-6
B = matrixdepot("chebspec", n, 1)
# B's eigenvalues have negative real parts.
v = real(eigvals(B))
for i = 1:n
    @test v[i] < 0
end
println("'chebspec' passed test...")


text:
include("Include.jl")

# extra:
using Flux
using Flux: @epochs
using BSON: @load

# load training set -
full_training_data_frame = load_tga_training_data()

# filter -
has_TM_flag = 0
experimental_data_table = filter([:visitid, :TM] => (x, y) -> ((x == 2 || x == 3) && y == has_TM_flag), full_training_data_frame)
input_data = convert.(Float32, Matrix(experimental_data_table[!, 3:14]))
validation_data_array = convert.(Float32, Matrix(experimental_data_table[!, 15:19]))

# what is the size of the system?
(P, D) = size(experimental_data_table)

# How are we going to do this?
δ = 0.90
number_of_samples = 1
indim = 12
outdim = 5
input_range_array = Array{Float32,2}(undef, indim, 2)
for i ∈ 1:indim
    L = (1 - δ) * minimum(input_data[:, i])
    U = (1 + δ) * maximum(input_data[:, i])
    input_range_array[i, 1] = L
    input_range_array[i, 2] = U
end

m = 2
model_name = "deep_coag_model-L$(m)O-TF-TM-$(has_TM_flag)-V2-V3.bson"
model_file_path = joinpath(_PATH_TO_MODELS, "network-1-epoch-12k", model_name)
@load model_file_path deep_coag_model

# baseline -
baseline = sample(deep_coag_model,input_data[1,:])

# create some storage -
tmp_dict = Dict()
sensitivity_dict = Dict()

for i ∈ 1:number_of_samples

    p_input_array = Array{Vector{Float32},1}()
    p_output_array = Array{Vector{Float32},1}()

    for j ∈ 1:indim

        # get the range for this input -
        L = input_range_array[j, 1]
        U = input_range_array[j, 2]

        # create a p + Δp -
        #r = rand()
        #v = L .+ r .* (U - L)

        # create a pertubed input data set -
        perturbed_input = copy(input_data[1,:])
        perturbed_input[j] = 1000.0

        # run the model -
        model_out = sample(deep_coag_model, perturbed_input)

        # grab -
        push!(p_input_array, perturbed_input)
        push!(p_output_array, model_out)
    end

    tmp_dict[i] = p_input_array => p_output_array
end

# get the data for a trial -
for t ∈ 1:number_of_samples

    p = tmp_dict[t].first
    x = tmp_dict[t].second

    # initialize -
    SA = Array{Float32,2}(undef, outdim, (indim-1))

    for i ∈ 1:outdim

        for j ∈ 1:(indim-1)

            # compute Δx -
            Δx = (x[j] .- baseline)[i]

            # compute Δp -
            Δp = (p[j] .- input_data[1,:])[j]

            # sensitivity -
            SF = 1.0
            SA[i,j] = SF*(Δx)/(Δp)
        end
    end

    sensitivity_dict[t] = SA
end

# process individual sensitivity arrays -
master_sens_array = zeros(outdim, indim-1)
for t ∈ 1:number_of_samples

    ASA = abs.(sensitivity_dict[t])

    for i ∈ 1:outdim
        for j ∈ 1:(indim - 1)
            old_value = master_sens_array[i,j]
            master_sens_array[i,j] = old_value + ASA[i,j]
        end
    end
end

max_stars_repo_name: varnerlab/UVM-TopDown-LegacyModel

max_issues_repo_name: varnerlab/UVM-TopDown-LegacyModel

max_forks_repo_name: varnerlab/UVM-TopDown-LegacyModel


text:
struct Timestamp{D, T} <: Dates.AbstractTime
    ts::D
    val::T
end

function Base.show(io::IO, ts::Timestamp)
    print(io, ts.ts, " | ", ts.val)
end

Dates.month(ts::Timestamp) = month(ts.ts)
Dates.year(ts::Timestamp) = year(ts.ts)
Dates.day(ts::Timestamp) = day(ts.ts)

# Should use something like AxisKeys.jl or NamedDims.jl
# single date

# function getindex(tsa::TimestampArray, d::TimeType)
#     counter = Int[]
#     for i in 1:length(tsa)
#         if d == tsa[i].timestamp
#             push!(counter, i)
#         end
#     end
#     return tsa[counter]
# end

# range of dates
#function getindex(tsa::TimestampArray, dates::Union(Vector{Date}, Vector{DateTime}))
#    counter   = Int[]
#    tsa_dates = TimeType[t.timestamp for t in tsa]
#  #  counter = int(zeros(length(dates)))
#    for i in 1:length(dates)
#        if findfirst(tsa_dates, dates[i]) != 0
#        #counter[i] = findfirst(ta.timestamp, dates[i])
#            push!(counter, findfirst(tsa_dates, dates[i]))
#        end
#    end
#    tsa[counter]
#end

# function getindex(tsa::TimestampArray, r::Union(StepRange{Date}, StepRange{DateTime}))
#     tsa[[r;]]
# end


text:
using DEBmicroTrait, Test

n_polymers = 5
n_enzymes  = 6
P          = rand(n_polymers)
E          = rand(n_enzymes)
#
Χ_0        = [10,10,10,10,10]
ρ_P        = [1.35, 1.5, 1.4, 0.4, 1.0]
R_E        = [3.29e-9, 4e-9, 4e-9, 4.4e-9, 5e-9,3.9e-9]
V_E        = [5, 0.24, 4,0.19, 3,4]
MW_E       = [48e3, 65e3, 59.7e3, 70e3, 50e3, 52e3]
D_E        = DEBmicroTrait.aqueous_diffusivity(MW_E)
#
R_P        = 20e-3*ones(n_polymers)
#
α_kin_P    = DEBmicroTrait.binding_sites(ρ_P, R_E)
K_EP_0     = DEBmicroTrait.polymer_affinity(R_E, R_P, V_E, D_E)

#######################################################################
#                     test paper values
# enzyme characteristics
R_E        = [3.29e-9]
V_E        = [5.0]
MW_E       = [48e3]
#D_E        = DEBmicroTrait.aqueous_diffusivity(MW_E)
D_E        = [1e-10]
# Polymer characteristics
ρ_P        = [1.35]
R_P        = [1].*1e-3
#
K_EP_0     = DEBmicroTrait.polymer_affinity(R_E, R_P, V_E, D_E)
@test K_EP_0 ≈ [0.244] atol=0.001
#
k_max_RP   = DEBmicroTrait.max_hydrolyis_rate(R_E, R_P, V_E)
#
R_P        = [1, 10, 20, 30, 40, 50].*1e-3
K_EP_0     = DEBmicroTrait.polymer_affinity(R_E, R_P, V_E, D_E)

V_E        = [42.3, 9.2, 0.41, 34.1, 31.5, 113.5, 128.6, 0.18, 112.2, 0.092, 4.6, 32896, 165, 92.8, 152.13, 69, 0.02, 54, 7590, 16.0, 18.6,
              300.3, 0.39, 0.39]

k_max_RP   = DEBmicroTrait.max_hydrolyis_rate(R_E, R_P, V_E)

kmaxs      = zeros(24,50)
R_P        = collect(range(1,50, length=50))

for i in 1:24
    for j in 1:50
        kmaxs[i,j] = DEBmicroTrait.max_hydrolyis_rate(R_E, [R_P[j]], [V_E[i]])[1]
    end
end

using JLD
save("/Users/glmarschmann/.julia/dev/DEBmicroTrait/files/cellulase_depolymerization.jld", "kmaxs", kmaxs)

using Statistics

median(V_E)

max_stars_repo_name: giannamars/DEBMicroTrait

max_issues_repo_name: giannamars/DEBMicroTrait

max_forks_repo_name: giannamars/DEBMicroTrait


text:
module SetBlasInt

using LinearAlgebra

export setblasint

"""
    setblasint(T, symbols...)

Re-define methods inside the LinearAlgebra.BLAS module to use the specified type of integer.
"""
function setblasint(T, symbs...)
    exprs = setblasint_exprs(T, symbs)
    @eval LinearAlgebra.BLAS begin
        $(exprs...)
    end
    length(exprs)
end

function setblasint_exprs(::Type{BlasInt}, symbs) where {BlasInt<:Union{Int32, Int64}}
    function blasfunc(x)
        BlasInt === Int32 ? Expr(:quote, x) : Expr(:quote, Symbol(x, "64_"))
    end
    exprs = []
    function pushfun(fun, expr)
        if endswith(String(fun), "_")
            fun = Symbol(String(fun)[1:end-1])
        end
        if symbs == (:all,) || fun in symbs
            push!(exprs, expr)
        end
    end

    # The code that appears as arguments to `pushfun` below is taken from Julia v1.7.0
    # https://github.com/JuliaLang/julia/blob/v1.7.0/stdlib/LinearAlgebra/src/blas.jl
    # with BlasInt changed to $BlasInt and @blasfunc($x) changed to $(blasfunc(x)).
    # (This is likely to be incompatible with other versions of Julia.)

    for (gemm, elty) in ((:dgemm_,:Float64),
                        (:sgemm_,:Float32),
                        (:zgemm_,:ComplexF64),
                        (:cgemm_,:ComplexF32))
        pushfun(gemm, quote
            function gemm!(transA::AbstractChar, transB::AbstractChar,
                       alpha::Union{($elty), Bool},
                       A::AbstractVecOrMat{$elty}, B::AbstractVecOrMat{$elty},
                       beta::Union{($elty), Bool},
                       C::AbstractVecOrMat{$elty})
                require_one_based_indexing(A, B, C)
                m = size(A, transA == 'N' ? 1 : 2)
                ka = size(A, transA == 'N' ? 2 : 1)
                kb = size(B, transB == 'N' ? 1 : 2)
                n = size(B, transB == 'N' ? 2 : 1)
                if ka != kb || m != size(C,1) || n != size(C,2)
                    throw(DimensionMismatch("A has size ($m,$ka), B has size ($kb,$n), C has size $(size(C))"))
                end
                chkstride1(A)
                chkstride1(B)
                chkstride1(C)
                ccall(($(blasfunc(gemm)), libblastrampoline), Cvoid,
                    (Ref{UInt8}, Ref{UInt8}, Ref{$BlasInt}, Ref{$BlasInt},
                    Ref{$BlasInt}, Ref{$elty}, Ptr{$elty}, Ref{$BlasInt},
                    Ptr{$elty}, Ref{$BlasInt}, Ref{$elty}, Ptr{$elty},
                    Ref{$BlasInt}, Clong, Clong),
                    transA, transB, m, n,
                    ka, alpha, A, max(1,stride(A,2)),
                    B, max(1,stride(B,2)), beta, C,
                    max(1,stride(C,2)), 1, 1)
                C
            end
        end)
    end

    for (fname, elty) in ((:dsyrk_,:Float64),
                            (:ssyrk_,:Float32),
                            (:zsyrk_,:ComplexF64),
                            (:csyrk_,:ComplexF32))
        pushfun(fname, quote
            function syrk!(uplo::AbstractChar, trans::AbstractChar,
                    alpha::Union{($elty), Bool}, A::AbstractVecOrMat{$elty},
                    beta::Union{($elty), Bool}, C::AbstractMatrix{$elty})
                require_one_based_indexing(A, C)
                n = checksquare(C)
                nn = size(A, trans == 'N' ? 1 : 2)
                if nn != n throw(DimensionMismatch("C has size ($n,$n), corresponding dimension of A is $nn")) end
                k  = size(A, trans == 'N' ? 2 : 1)
                chkstride1(A)
                chkstride1(C)
                ccall(($(blasfunc(fname)), libblastrampoline), Cvoid,
                (Ref{UInt8}, Ref{UInt8}, Ref{$BlasInt}, Ref{$BlasInt},
                    Ref{$elty}, Ptr{$elty}, Ref{$BlasInt}, Ref{$elty},
                    Ptr{$elty}, Ref{$BlasInt}, Clong, Clong),
                uplo, trans, n, k,
                alpha, A, max(1,stride(A,2)), beta,
                C, max(1,stride(C,2)), 1, 1)
                C
            end
        end)
    end

    exprs
end

end # module


text:
using SIAN

println("Setting up the problem")

ode = @ODEmodel(
    s'(t) = mu - bi * s(t) * i(t) - bw * s(t) * w(t) - mu * s(t) + al * r(t),
    i'(t) = bw * s(t) * w(t) + bi * s(t) * i(t) - g * i(t) - mu * i(t),
    w'(t) = dz * (i(t) - w(t)),
    r'(t) = g * i(t) - mu * r(t) - al * r(t),
    y1(t) = k * i(t),
    y2(t) = i(t) + r(t) + s(t)
)

res = identifiability_ode(ode, get_parameters(ode); p=0.99, p_mod=2^29 - 3, nthrds=1)

println(res)

max_stars_repo_name: iliailmer/SIAN-Julia

max_issues_repo_name: iliailmer/SIAN-Julia

max_forks_repo_name: iliailmer/SIAN-Julia


text:
#EasyData: A quick/easy way to save datasets/plots to file.
#-------------------------------------------------------------------------------
module EasyData

import ..CMDimData
using CMDimData.EasyPlot
using MDDatasets
using Colors
import HDF5
import Pkg

import CMDimData.EasyPlot: PlotCollection, Plot, YStrip, GridCartesian
import CMDimData.EasyPlot: Waveform, LineAttributes, GlyphAttributes
import CMDimData.EasyPlot: Extents1D, Axis, FoldedAxis
import CMDimData.EasyPlot: TextAnnotation, HVMarker

include("base.jl")
include("hdf5typed.jl")
include("hdf5datamd.jl")
include("hdf5plots.jl")

#==Initialization
===============================================================================#
function __init__()
	#Needed by @includepkg
end

#==Un-exported public interface to read/write data.
================================================================================
	openwriter(path::String)
	openreader(path::String)
	Base.close(w::EasyDataWriter)
	Base.close(r::EasyDataReader)

	Base.write(w::EasyDataWriter, d::DataMD, name::String)
	Base.write(w::EasyDataWriter, pcoll::PlotCollection, name::String)

	Base.read(::Type{DataMD}, r::EasyDataReader, name::String)
	Base.read(::Type{PlotCollection}, r::EasyDataReader, name::String)
	readdata(r::EasyDataReader, name::String) #Avoids specifying ::Type{DataMD}
	readplot(r::EasyDataReader, name::String) #Avoids specifying ::Type{PlotCollection}

	#High-level open, read/write, close interface for plots:
	writeplot(filepath::String, pcoll::PlotCollection; name::String="_unnamed")
	readplot(filepath::String; name::String="_unnamed")
==#

end #EasyData

#Last line


text:
module TransducersBenchmarksReports

using BenchmarkTools
using PkgBenchmark
using Run
using Transducers

benchpath(pathparts...) =
    joinpath(dirname(dirname(pathof(Transducers))), "benchmark", pathparts...)

loadresult() = PkgBenchmark.readresults(benchpath("result.json"))

function mayberun()
    resultpath = benchpath("result.json")
    if isfile(resultpath)
        @info "result.json already exists at: $resultpath"
    else
        Run.script(benchpath("runbenchmarks.jl"))
    end
    return loadresult()
end

end # module


text:
#### annotation
# P, F, C
# id -> vec of ids

# function load_annot_short()
#
#
#
# end

function load_annot_long(file, ontol ; genecol=:Gene, gocol=:GOID)

    table = CSV.read(file)
    annot = Dict(a => Dict{String, Vector{Int}}() for a in ["P", "F", "C"])
    namedict = Dict(biological_process => "P", molecular_function => "F", cellular_component => "C")
    missingids = 0
    annots = 0
    for (gene, goid) in zip(table[!, genecol], table[!, gocol])
        gid = goid_to_num(goid)
        track = false
        if (gid == 45545) && (gene == "gene6898|chrd")
            @show gene, goid, gid
            track = true
        end
        if !haskey(ontol, gid)
            missingids += 1
            if track
                println("its missing")
                error("")
            end
            continue
        end

        gid = ontol[gid].id ## remove alt ids
        if track
            @show "U", gene, goid, gid
        end
        ns = namedict[ontol[gid].namespace]

        !haskey(annot[ns], gene)  && (annot[ns][gene] = Int[])
        if gid ∉ annot[ns][gene]
            push!(annot[ns][gene], gid)
            annots += 1
        end
    end

    println("Loaded: $(file)\t", annots, " annotations, ", missingids, " not found in ontology.");
    annot

end

reverseannotation(annot) = Dict(o => reverse_ontol_annot(ann) for (o, ann) in annot)

function reverse_ontol_annot(ann)
    r = Dict{Int, Vector{String}}()

    for (gene, ids) in ann
        for id in ids
            r[id] = push!(get(r, id, String[]), gene)
        end
    end
    r
end


text:
using Oceananigans.TimeSteppers: AbstractTimeStepper

struct WickerSkamarockTimeStepper{S, F, I} <: AbstractTimeStepper
      slow_source_terms :: S
      fast_source_terms :: F
    intermediate_fields :: I

    function WickerSkamarockTimeStepper(arch, grid, tracers;
        slow_source_terms = SlowSourceTermFields(arch, grid, tracers),
        fast_source_terms = FastSourceTermFields(arch, grid, tracers),
      intermediate_fields = FastSourceTermFields(arch, grid, tracers))

      return new{typeof(slow_source_terms), typeof(fast_source_terms),
        typeof(intermediate_fields)}(slow_source_terms, fast_source_terms, intermediate_fields)
  end
end

function SlowSourceTermFields(arch, grid, tracernames)
    ρu = XFaceField(arch, grid)
    ρv = YFaceField(arch, grid)
    ρw = ZFaceField(arch, grid)
    tracers = TracerFields(tracernames, arch, grid, ())
    return (ρu = ρu, ρv = ρv, ρw = ρw, tracers = tracers)
end

function FastSourceTermFields(arch, grid, tracernames)
    ρu = XFaceField(arch, grid)
    ρv = YFaceField(arch, grid)
    ρw = ZFaceField(arch, grid)
    tracers = TracerFields(tracernames, arch, grid, ())
    return (ρu = ρu, ρv = ρv, ρw = ρw, tracers = tracers)
end


text:
using RecipesBase

@userplot PsfPlot

@recipe function f(p::PsfPlot)
    model = p.args[1]
    inds = p.args[2:end]
    # if `inds` is a vector/tuple, this effectively unpacks it
    if length(inds) == 1
        inds = inds[1]
    end

    seriestype := :heatmap
    aspect_ratio --> 1
    xlims --> extrema(first(inds))
    ylims --> extrema(last(inds))
    xguide --> "x"
    yguide --> "y"

    arr = map(model, CartesianIndices(inds))

    return inds..., transpose(arr)
end


text:
# -------- TODO!!! --------

# Add tests.

# Add docs.

# Add get, get!, haskey.

# Implement Edge properties.

# Implement weights functionality.

# Rename labels() to keys() ?
# Although if you have both vertex properties and edge properties, it's
# less clear what the keys are and what the values are, since there's actually
# two dictionaries.

# Perhaps use labels() and edge_labels() instead of keys() ?

# Add vertex_properties() and edge_properties() instead of values() ?

# Add show() method.

# Add merge, join, union, etc, methods.

# ---------------------------

struct NoVertexProperties end
struct NoEdgeProperties end

# V and E are any types with getproperty() defined.
struct PropertyGraph{T, G <: AbstractGraph{T}, L, V, E} <: AbstractGraph{T}
    g::G
    vmap::Bijection{L, T}
    vprops::Dict{L, V}    # Make these Union{No*, Dict{T, R}} etc?
    eprops::Dict{L, E}    # This should actually be something like Dict{Edge, EP} or Dict{Tuple{VL,VL}, EP}
end

Base.broadcastable(pg::PropertyGraph) = Ref(pg)

function PropertyGraph{T, G, L, V, E}() where {T, G, L, V, E}
    PropertyGraph(
        G(),
        Bijection{L, T}(),
        Dict{L, V}(),
        Dict{L, E}()
    )
end

function PropertyGraph{T, G, L, V}() where {T, G, L, V}
    PropertyGraph{T, G, L, V, NoEdgeProperties}()
end

function PropertyGraph{T, G, L}() where {T, G, L}
    PropertyGraph{T, G, L, NoVertexProperties}()
end

SimplePropertyGraph{T, L, V, E} = PropertyGraph{T, SimpleGraph{T}, L, V, E}
SimplePropertyDiGraph{T, L, V, E} = PropertyGraph{T, SimpleDiGraph{T}, L, V, E}
const SimpleAlias = Union{SimplePropertyGraph, SimplePropertyDiGraph}

(pg::PropertyGraph)(label) = pg.vmap[label]
(pg::SimpleAlias)(src_label, dst_label) = Edge(pg(src_label), pg(dst_label))

label(pg::PropertyGraph, code) = pg.vmap(code)

labels(pg::PropertyGraph) = domain(pg.vmap)
Base.in(label, pg::PropertyGraph) = (label in labels(pg))

Base.eltype(::PropertyGraph{T}) where {T} = T
LightGraphs.edges(pg::PropertyGraph) = edges(pg.g)
LightGraphs.vertices(pg::PropertyGraph) = vertices(pg.g)
LightGraphs.nv(pg::PropertyGraph) = nv(pg.g)
LightGraphs.ne(pg::PropertyGraph) = ne(pg.g)

LightGraphs.inneighbors(pg::PropertyGraph, code::Integer) = inneighbors(pg.g, code)
LightGraphs.inneighbors(pg::PropertyGraph{T, G, L}, label::L) where {T, G, L} = inneighbors(pg.g, pg(label))

LightGraphs.outneighbors(pg::PropertyGraph, code::Integer) = outneighbors(pg.g, code)
LightGraphs.outneighbors(pg::PropertyGraph{T, G, L}, label::L) where {T, G, L} = outneighbors(pg.g, pg(label))

LightGraphs.has_vertex(pg::PropertyGraph, code::Integer) = has_vertex(pg.g, code)
LightGraphs.has_vertex(pg::PropertyGraph{T, G, L}, label::L) where {T, G, L} = (label in pg)

LightGraphs.has_edge(pg::PropertyGraph, u::Integer, v::Integer) = has_edge(pg.g, u, v)

function LightGraphs.has_edge(pg::PropertyGraph{T, G, L}, u::L, v::L) where {T, G, L}
    if u ∉ pg || v ∉ pg
        false
    else
        has_edge(pg.g, pg(u), pg(v))
    end
end

function LightGraphs.is_directed(::Type{PropertyGraph{T,G,L,V,E}}) where {T,G,L,V,E}
    is_directed(G)
end

function LightGraphs.add_vertex!(pg::SimpleAlias, label)
    if label in pg
        return false
    end

    added = add_vertex!(pg.g)
    if added
        pg.vmap[label] = nv(pg)
    end

    added
end

LightGraphs.add_edge!(pg::SimpleAlias, e::Edge) = add_edge!(pg.g, e)

function LightGraphs.add_edge!(pg::SimpleAlias, src_label, dst_label)
    if  src_label in pg  &&  dst_label in pg
        add_edge!(pg, pg(src_label, dst_label))
    else
        return false
    end
end

Base.setindex!(pg::PropertyGraph, val, label) = ( pg.vprops[label] = val )
Base.getindex(pg::PropertyGraph, label) = pg.vprops[label]

function LightGraphs.rem_vertex!(pg::SimpleAlias, _label)
    code = pg(_label)
    last_code = nv(pg)
    removed = rem_vertex!(pg.g, code)

    if removed
        delete!(pg.vprops, _label)
        last_label = label(pg, last_code)

        delete!(pg.vmap, _label)

        if _label != last_label
            # Bijections.jl requires deleting a key before changing its value.
            delete!(pg.vmap, last_label)
            pg.vmap[last_label] = code
        end
    end

    removed
end

function LightGraphs.rem_edge!(pg::SimpleAlias, src_label, dst_label)
    removed = rem_edge!(pg.g, pg(src_label, dst_label))
    if removed
        # Delete edge properties.
    end
    removed
end

Base.reverse!(pg::SimpleAlias) = reverse!(pg.g)

function Base.copy(pg::PropertyGraph)
    PropertyGraph(
        copy(pg.g),
        Bijection(Dict(pg.vmap)),
        copy(pg.vprops),
        copy(pg.eprops)
    )
end


text:
using LinearAlgebra
using DifferentialEquations
using Distributions
using Random
using Printf
using Plots
using RCall
using CPUTime
using StatsBase
using JLD
# using StatsPlots
include("utils.jl")
include("sir_ode.jl")
include("abc.jl")
include("abcMCMC.jl")
include("abcSMC.jl")
include("BayesianCalibration.jl")

# Bayesian Calibration
algo_list = (ABC=ABC, ABC_MCMC=ABC_MCMC, ABC_SMC=ABC_SMC)

true_prior = (truncated(Normal(2,0.5),0,Inf),truncated(Normal(1,0.25),0,Inf))

algo_parameter_ABC = (prior = true_prior,
                     epsilon = 0.1,eta = identity_mapping, d = compute_norm)

algo_parameter_mcmc = (prior = true_prior,
                        epsilon = 1, eta = identity_mapping,
                        d= compute_norm, proposal = proposal_Normal,
                        proposalRatio = proposalRatio_Normal,sd = (0.5,0.25),
                        thinning = 10, verbose=false,
                        burn_in = 100)

algo_parameters_smc = (prior = true_prior,
                        time_final=5, eps_list =[10, 5, 1, 0.5, 0.1],
                        eta= identity_mapping, d= distanceFunction,
                        kernel=proposal_Normal,
                        kernel_density=proposal_Normal_density,
                        sd=(0.5,0.25), resample_method=systematic_resample,
                        verbose=false)

algo_param_list = (ABC = algo_parameter_ABC,
                    ABC_MCMC = algo_parameter_mcmc,
                    ABC_SMC = algo_parameters_smc)

simulation = BayesianCalibration(Int(2^10),Int(250),0.10,algo_list,algo_param_list,
                                ode_model = solve_ode,
                                initial_state = [50000.0;4.5;350]./50000,
                                time_window=(0,20.0),
                                add_noise = false,
                                true_p_dist=true_prior)

y = solve_ode([50000.0;4.5;350]./50000,(0,20.0),rand.(true_prior))
plot(y[1,:])
plot!(y[2,:])
plot!(y[3,:])

# calibration
show_calibration(simulation)

results_dir = "results/simulation_SIR/"

save(results_dir*"output.jld","output",simulation)

using DataFrames
simulation = load(results_dir*"output.jld")["output"]

function to_dataframe(x)
    tmp = DataFrame(x)
    rename!(tmp,[:ABC,:MCMC,:SMC])
end

CSV.write(results_dir*"calibration.csv",to_dataframe(simulation[2]))
CSV.write(results_dir*"cpu.csv",to_dataframe(simulation[4]))
CSV.write(results_dir*"ess.csv",to_dataframe(simulation[5]))
CSV.write(results_dir*"ess_cpu.csv",to_dataframe(simulation[6]))


text:
abstract type AbstractGroup end
abstract type AbstractAbelianGroup <: AbstractGroup end

abstract type AbstractCharge end
struct Trivial <: AbstractCharge end
abstract type AbstractNonTrivialCharge end

struct Z2Charge <: AbstractCharge
    charge :: Bool
end
"U1Charge for U(1) abelian group"
struct U1Charge <: AbstractCharge
    charge :: Int
end
const U1 = U1Charge

@inline Base.inv(c::U1Charge) = U1Charge(-c.charge)
@inline Base.:*(a::Int, c::U1Charge) = U1Charge(a*c.charge)
@inline Base.:+(c1::U1Charge, c2::U1Charge) = U1Charge(c1.charge + c2.charge)
@inline Base.:+(c::U1Charge, a::Int) = U1Charge(c.charge + a)
@inline Base.:+(c::U1Charge) = c
@inline Base.:-(c1::U1Charge, c2::U1Charge) = U1Charge(c1.charge - c2.charge)
@inline Base.:-(c::U1Charge, a::Int) = U1Charge(c.charge - a)
@inline Base.:-(c::U1Charge) = inv(c)
@inline Base.div(c::U1Charge, a::Int) = div(c.charge, a)
@inline Base.zero(::Type{U1Charge}) = U1Charge(0)
@inline Base.isodd(c::U1Charge) = isodd(c.charge)
@inline Base.isless(c1::U1Charge, c2::U1Charge) = isless(c1.charge, c2.charge)
@inline Base.isequal(c1::U1Charge, c2::U1Charge) = isequal(c1.charge, c2.charge)
@inline Base.isequal(c1::U1Charge, c2::Int) = isequal(c1.charge, c2)
@inline Base.length(::U1Charge) = 1
@inline Base.iterate(c::U1Charge) = (c, nothing)
convert(::Type{U1Charge}, s::Int) = U1Charge(s)
convert(::Type{Int}, c::U1Charge) = c.charge

function Base.show(io::IO, c::AbstractCharge)
    print(io, c.charge)
    nothing
end

max_stars_repo_name: amiragha/mps

max_issues_repo_name: amiragha/mps

max_forks_repo_name: amiragha/mps


text:
import Images
using Color, Base.Test, FixedPointNumbers

# Comparison of each element in arrays with a scalar
approx_equal(ar, v) = all(abs(ar.-v) .< sqrt(eps(v)))
approx_equal(ar::Images.AbstractImage, v) = approx_equal(Images.data(ar), v)

# arithmetic
img = convert(Images.Image, zeros(3,3))
img2 = (img .+ 3)/2
@test all(img2 .== 1.5)
img3 = 2img2
@test all(img3 .== 3)
img3 = copy(img2)
img3[img2 .< 4] = -1
@test all(img3 .== -1)
img = convert(Images.Image, rand(3,4))
A = rand(3,4)
img2 = img .* A
@test all(Images.data(img2) == Images.data(img).*A)
img2 = convert(Images.Image, A)
img2 = img2 .- 0.5
img3 = 2img .* img2
img2 = img ./ A
img2 = (2img).^2
# Same operations with ColorValue images
img = Images.colorim(zeros(Float32,3,4,5))
img2 = (img .+ RGB{Float32}(1,1,1))/2
@test all(img2 .== RGB{Float32}(1,1,1)/2)
img3 = 2img2
@test all(img3 .== RGB{Float32}(1,1,1))
A = fill(2, 4, 5)
@test all(A.*img2 .== fill(RGB{Float32}(1,1,1), 4, 5))
img2 = img2 .- RGB{Float32}(1,1,1)/2
A = rand(Uint8,3,4)
img = reinterpret(Images.Gray{Ufixed8}, Images.grayim(A))
imgm = mean(img)
imgn = img/imgm
@test_approx_eq reinterpret(Float32, Images.data(imgn)) float32(A/mean(A))

# Reductions
let
    A = rand(5,5,3)
    img = Images.colorim(A, "RGB")
    s12 = sum(img, (1,2))
    @test Images.colorspace(s12) == "RGB"
    s3 = sum(img, (3,))
    @test Images.colorspace(s3) == "Unknown"
    A = [NaN, 1, 2, 3]
    @test_approx_eq Images.meanfinite(A, 1) [2]
    A = [NaN 1 2 3;
         NaN 6 5 4]
    @test_approx_eq Images.meanfinite(A, 1) [NaN 3.5 3.5 3.5]
    @test_approx_eq Images.meanfinite(A, 2) [2, 5]'
    @test_approx_eq Images.meanfinite(A, (1,2)) [3.5]
    @test Images.minfinite(A) == 1
    @test Images.maxfinite(A) == 6
    @test Images.maxabsfinite(A) == 6
    A = float32(rand(3,5,5))
    img = Images.colorim(A, "RGB")
    dc = Images.data(Images.meanfinite(img, 1))-reinterpret(RGB{Float32}, mean(A, 2), (1,5))
    @test maximum(map(abs, dc)) < 1e-6
end

# fft and ifft
A = rand(Float32, 3, 5, 6)
img = Images.colorim(A)
imgfft = fft(img)
@test_approx_eq Images.data(imgfft) fft(A, 2:3)
@test Images.colordim(imgfft) == 1
img2 = ifft(imgfft)
@test_approx_eq img2 reinterpret(Float32, img)

# Array padding
let A = [1 2; 3 4]
    @test Images.padindexes(A, 1, 0, 0, "replicate") == [1,2]
    @test Images.padindexes(A, 1, 1, 0, "replicate") == [1,1,2]
    @test Images.padindexes(A, 1, 0, 1, "replicate") == [1,2,2]
    @test Images.padarray(A, (0,0), (0,0), "replicate") == A
    @test Images.padarray(A, (1,2), (2,0), "replicate") == [1 1 1 2; 1 1 1 2; 3 3 3 4; 3 3 3 4; 3 3 3 4]
    @test Images.padindexes(A, 1, 1, 0, "circular") == [2,1,2]
    @test Images.padindexes(A, 1, 0, 1, "circular") == [1,2,1]
    @test Images.padarray(A, [2,1], [0,2], "circular") == [2 1 2 1 2; 4 3 4 3 4; 2 1 2 1 2; 4 3 4 3 4]
    @test Images.padindexes(A, 1, 1, 0, "symmetric") == [1,1,2]
    @test Images.padindexes(A, 1, 0, 1, "symmetric") == [1,2,2]
    @test Images.padarray(A, (1,2), (2,0), "symmetric") == [2 1 1 2; 2 1 1 2; 4 3 3 4; 4 3 3 4; 2 1 1 2]
    @test Images.padarray(A, (1,2), (2,0), "value", -1) == [-1 -1 -1 -1; -1 -1 1 2; -1 -1 3 4; -1 -1 -1 -1; -1 -1 -1 -1]
    A = [1 2 3; 4 5 6]
    @test Images.padindexes(A, 2, 1, 0, "reflect") == [2,1,2,3]
    @test Images.padindexes(A, 2, 0, 1, "reflect") == [1,2,3,2]
    @test Images.padarray(A, (1,2), (2,0), "reflect") == [6 5 4 5 6; 3 2 1 2 3; 6 5 4 5 6; 3 2 1 2 3; 6 5 4 5 6]
    A = [1 2; 3 4]
    @test Images.padarray(A, (1,1)) == [1 1 2 2; 1 1 2 2; 3 3 4 4; 3 3 4 4]
    @test Images.padarray(A, (1,1), "replicate", "both") == [1 1 2 2; 1 1 2 2; 3 3 4 4; 3 3 4 4]
    @test Images.padarray(A, (1,1), "circular", "pre") == [4 3 4; 2 1 2; 4 3 4]
    @test Images.padarray(A, (1,1), "symmetric", "post") == [1 2 2; 3 4 4; 3 4 4]
    A = ["a" "b"; "c" "d"]
    @test Images.padarray(A, (1,1)) == ["a" "a" "b" "b"; "a" "a" "b" "b"; "c" "c" "d" "d"; "c" "c" "d" "d"]
    @test_throws ErrorException Images.padindexes(A, 1, 1, 1, "unknown")
    @test_throws ErrorException Images.padarray(A, (1,1), "unknown")
end

# filtering
EPS = 1e-14
imgcol = Images.colorim(rand(3,5,6))
imgcolf = convert(Images.Image{RGB{Ufixed8}}, imgcol)
for T in (Float64, Int)
    A = zeros(T,3,3); A[2,2] = 1
    kern = rand(3,3)
    @test maximum(abs(Images.imfilter(A, kern) - rot180(kern))) < EPS
    kern = rand(2,3)
    @test maximum(abs(Images.imfilter(A, kern)[1:2,:] - rot180(kern))) < EPS
    kern = rand(3,2)
    @test maximum(abs(Images.imfilter(A, kern)[:,1:2] - rot180(kern))) < EPS
end
kern = zeros(3,3); kern[2,2] = 1
@test maximum(map(abs, imgcol - Images.imfilter(imgcol, kern))) < EPS
@test maximum(map(abs, imgcolf - Images.imfilter(imgcolf, kern))) < EPS
for T in (Float64, Int)
    # Separable kernels
    A = zeros(T,3,3); A[2,2] = 1
    kern = rand(3).*rand(3)'
    @test maximum(abs(Images.imfilter(A, kern) - rot180(kern))) < EPS
    kern = rand(2).*rand(3)'
    @test maximum(abs(Images.imfilter(A, kern)[1:2,:] - rot180(kern))) < EPS
    kern = rand(3).*rand(2)'
    @test maximum(abs(Images.imfilter(A, kern)[:,1:2] - rot180(kern))) < EPS
end
A = zeros(3,3); A[2,2] = 1
kern = rand(3,3)
@test maximum(abs(Images.imfilter_fft(A, kern) - rot180(kern))) < EPS
kern = rand(2,3)
@test maximum(abs(Images.imfilter_fft(A, kern)[1:2,:] - rot180(kern))) < EPS
kern = rand(3,2)
@test maximum(abs(Images.imfilter_fft(A, kern)[:,1:2] - rot180(kern))) < EPS
kern = zeros(3,3); kern[2,2] = 1
@test maximum(map(abs, imgcol - Images.imfilter_fft(imgcol, kern))) < EPS
@test maximum(map(abs, imgcolf - Images.imfilter_fft(imgcolf, kern))) < EPS

@test approx_equal(Images.imfilter(ones(4,4), ones(3,3)), 9.0)
@test approx_equal(Images.imfilter(ones(3,3), ones(3,3)), 9.0)
@test approx_equal(Images.imfilter(ones(3,3), [1 1 1;1 0.0 1;1 1 1]), 8.0)
img = convert(Images.Image, ones(4,4))
@test approx_equal(Images.imfilter(img, ones(3,3)), 9.0)
A = zeros(5,5,3); A[3,3,[1,3]] = 1
@test Images.colordim(A) == 3
kern = rand(3,3)
kernpad = zeros(5,5); kernpad[2:4,2:4] = kern
Af = Images.imfilter(A, kern)

@test_approx_eq Af cat(3, rot180(kernpad), zeros(5,5), rot180(kernpad))
Aimg = permutedims(convert(Images.Image, A), [3,1,2])
@test_approx_eq Images.imfilter(Aimg, kern) permutedims(Af, [3,1,2])
@test approx_equal(Images.imfilter(ones(4,4),ones(1,3),"replicate"), 3.0)

A = zeros(5,5); A[3,3] = 1
kern = rand(3,3)
Af = Images.imfilter(A, kern, "inner")
@test Af == rot180(kern)
Afft = Images.imfilter_fft(A, kern, "inner")
@test_approx_eq Af Afft
h = [0.24,0.87]
@test_approx_eq imfilter(eye(3), h, "inner") imfilter_fft(eye(3), h, "inner")  # issue #204

@test approx_equal(Images.imfilter_gaussian(ones(4,4), [5,5]), 1.0)
A = fill(nan(Float32), 3, 3)
A[1,1] = 1
A[2,1] = 2
A[3,1] = 3
@test Images.imfilter_gaussian(A, [0,0]) === A
@test_approx_eq Images.imfilter_gaussian(A, [0,3]) A
B = copy(A)
B[isfinite(B)] = 2
@test_approx_eq Images.imfilter_gaussian(A, [10^3,0]) B
@test maximum(map(abs, Images.imfilter_gaussian(imgcol, [10^3,10^3]) - mean(imgcol))) < 1e-4
@test maximum(map(abs, Images.imfilter_gaussian(imgcolf, [10^3,10^3]) - mean(imgcolf))) < 1e-4
A = rand(4,5)
img = reinterpret(Gray{Float64}, Images.grayim(A))
imgf = Images.imfilter_gaussian(img, [2,2])
@test_approx_eq reinterpret(Float64, data(imgf)) Images.imfilter_gaussian(A, [2,2])
A = rand(3,4,5)
img = Images.colorim(A)
imgf = Images.imfilter_gaussian(img, [2,2])
@test_approx_eq reinterpret(Float64, data(imgf)) Images.imfilter_gaussian(A, [0,2,2])

A = zeros(Int, 9, 9); A[5, 5] = 1
@test maximum(abs(Images.imfilter_LoG(A, [1,1]) - Images.imlog(1.0))) < EPS

@test Images.imaverage() == fill(1/9, 3, 3)
@test Images.imaverage([3,3]) == fill(1/9, 3, 3)
@test_throws ErrorException Images.imaverage([5])

# restriction
A = reshape(uint16(1:60), 4, 5, 3)
B = Images.restrict(A, (1,2))
Btarget = cat(3, [ 0.96875  4.625   5.96875;
                   2.875   10.5    12.875;
                   1.90625  5.875   6.90625],
                 [ 8.46875  14.625 13.46875;
                  17.875    30.5   27.875;
                   9.40625  15.875 14.40625],
                 [15.96875  24.625 20.96875;
                  32.875    50.5   42.875;
                  16.90625  25.875 21.90625])
@test_approx_eq B Btarget
Argb = reinterpret(RGB, reinterpret(Ufixed16, permutedims(A, (3,1,2))))
B = Images.restrict(Argb)
Bf = permutedims(reinterpret(Float64, B), (2,3,1))
@test_approx_eq_eps Bf Btarget/reinterpret(one(Ufixed16)) 1e-12
Argba = reinterpret(RGBA{Ufixed16}, reinterpret(Ufixed16, A))
B = Images.restrict(Argba)
@test_approx_eq_eps reinterpret(Float64, B) Images.restrict(A, (2,3))/reinterpret(one(Ufixed16)) 1e-12
A = reshape(1:60, 5, 4, 3)
B = Images.restrict(A, (1,2,3))
@test_approx_eq B cat(3, [ 2.6015625  8.71875 6.1171875;
                           4.09375   12.875   8.78125;
                           3.5390625 10.59375 7.0546875],
                         [10.1015625 23.71875 13.6171875;
                          14.09375   32.875   18.78125;
                          11.0390625 25.59375 14.5546875])
Images.restrict(imgcol, (1,2))

# erode/dilate
A = zeros(4,4,3)
A[2,2,1] = 0.8
A[4,4,2] = 0.6
Ae = Images.erode(A)
@test Ae == zeros(size(A))
Ad = Images.dilate(A)
Ar = [0.8 0.8 0.8 0;
      0.8 0.8 0.8 0;
      0.8 0.8 0.8 0;
      0 0 0 0]
Ag = [0 0 0 0;
      0 0 0 0;
      0 0 0.6 0.6;
      0 0 0.6 0.6]
@test Ad == cat(3, Ar, Ag, zeros(4,4))
Ae = Images.erode(Ad)
Ar = [0.8 0.8 0 0;
      0.8 0.8 0 0;
      0 0 0 0;
      0 0 0 0]
Ag = [0 0 0 0;
      0 0 0 0;
      0 0 0 0;
      0 0 0 0.6]
@test Ae == cat(3, Ar, Ag, zeros(4,4))

# opening/closing
A = zeros(4,4,3)
A[2,2,1] = 0.8
A[4,4,2] = 0.6
Ao = Images.opening(A)
@test Ao == zeros(size(A))
A = zeros(10,10)
A[4:7,4:7] = 1
B = copy(A)
A[5,5] = 0
Ac = Images.closing(A)
@test Ac == B

# label_components
A = [true  true  false true;
     true  false true  true]
lbltarget = [1 1 0 2;
             1 0 2 2]
lbltarget1 = [1 2 0 4;
              1 0 3 4]
@test Images.label_components(A) == lbltarget
@test Images.label_components(A, [1]) == lbltarget1
connectivity = [false true  false;
                true  false true;
                false true  false]
@test Images.label_components(A, connectivity) == lbltarget
connectivity = trues(3,3)
lbltarget2 = [1 1 0 1;
              1 0 1 1]
@test Images.label_components(A, connectivity) == lbltarget2

# phantoms

P = [ 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0;
      0.0  0.0  1.0  0.2  0.2  1.0  0.0  0.0;
      0.0  0.0  0.2  0.3  0.3  0.2  0.0  0.0;
      0.0  0.0  0.2  0.0  0.2  0.2  0.0  0.0;
      0.0  0.0  0.2  0.0  0.0  0.2  0.0  0.0;
      0.0  0.0  0.2  0.2  0.2  0.2  0.0  0.0;
      0.0  0.0  1.0  0.2  0.2  1.0  0.0  0.0;
      0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0 ]

Q = Images.shepp_logan(8)
@test norm((P-Q)[:]) < 1e-10

P = [ 0.0  0.0  0.0   0.0   0.0   0.0   0.0  0.0;
      0.0  0.0  2.0   1.02  1.02  2.0   0.0  0.0;
      0.0  0.0  1.02  1.03  1.03  1.02  0.0  0.0;
      0.0  0.0  1.02  1.0   1.02  1.02  0.0  0.0;
      0.0  0.0  1.02  1.0   1.0   1.02  0.0  0.0;
      0.0  0.0  1.02  1.02  1.02  1.02  0.0  0.0;
      0.0  0.0  2.0   1.02  1.02  2.0   0.0  0.0;
      0.0  0.0  0.0   0.0   0.0   0.0   0.0  0.0 ]

Q = Images.shepp_logan(8,highContrast=false)
@test norm((P-Q)[:]) < 1e-10


text:
using OperationsResearchApplications
using Documenter

makedocs(;
    modules=[OperationsResearchApplications],
    authors="Ethan Buck",
    repo="https://github.com/egbuck/OperationsResearchApplications.jl/blob/{commit}{path}#L{line}",
    sitename="OperationsResearchApplications.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://egbuck.github.io/OperationsResearchApplications.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/egbuck/OperationsResearchApplications.jl",
)


text:
# SPDX-License-Identifier: MIT

"""
    struct MatrixFSM{T<:Semiring}
        π        # vector of initial probabilities
        T        # matrix of transition probabilities
        Tᵀ       # transpose of `T`
        C        # matrix mapping state -> pdfindex
        Cᵀ       # tranpose of `C`
        labels   # vector of labels (label can be `String`, `Tuple`, etc.)
    end

FSM where the arcs are stored as a sparse matrix. This type of FSM
is used by the inference algorigthms. In addition to the structure of
the FSM, the `MatrixFSM` store the mapping between states and `pdf-id`
in the matrix `C`.

!!! warning
    In this format, the final state is just a "regular" state.
    Therefore, the dimension of `π` (similarly for `T` and `C`)
    will have one more dimension.

# Constructor

    MatrixFSM(fsm, pdfid_mapping)

"""
struct MatrixFSM{Tv} <: AbstractFSM{Tv}
    π::AbstractSparseVector{Tv}
    T::AbstractSparseMatrix{Tv}
    Tᵀ::AbstractSparseMatrix{Tv}
    C::AbstractSparseMatrix{Tv}
    Cᵀ::AbstractSparseMatrix{Tv}
    labels::Vector
end

function MatrixFSM(fsm::AbstractFSM{Tv}, pdfid_mapping::Dict,
                   keyfn::Function = identity) where Tv
    state2idx = Dict(s => i for (i, s) in enumerate(states(fsm)))
    S = length(fsm)

    π = spzeros(Tv, S+1)
    for s in filter(isinit, states(fsm))
        π[state2idx[s]] = s.initweight
    end

    T = spzeros(Tv, S+1, S+1)
    Tᵀ = spzeros(Tv, S+1, S+1)
    for src in states(fsm)
        for arc in arcs(fsm, src)
            T[state2idx[src], state2idx[arc.dest]] = arc.weight
            Tᵀ[state2idx[arc.dest], state2idx[src]] = arc.weight
        end
    end

    # Special case for the final states.
    for src in filter(isfinal, states(fsm))
        T[state2idx[src], end] = src.finalweight
        Tᵀ[end, state2idx[src]] = src.finalweight
    end
    T[end,end] = one(Tv)
    Tᵀ[end,end] = one(Tv)

    K = maximum(values(pdfid_mapping))
    C = spzeros(Tv, S+1, K+1)
    Cᵀ = spzeros(Tv, K+1, S+1)
    labels = Vector{Any}(undef, S+1)
    for s in states(fsm)
        C[state2idx[s], pdfid_mapping[keyfn(s.label)]] = one(Tv)
        Cᵀ[pdfid_mapping[keyfn(s.label)], state2idx[s]] = one(Tv)
        labels[state2idx[s]] = s.label
    end
    C[end,end] = one(Tv)
    Cᵀ[end,end] = one(Tv)
    labels[end] = nothing

    MatrixFSM(π, T, Tᵀ, C, Cᵀ, labels)
end

function states(fsm::MatrixFSM)
    S = length(fsm.π) - 1
    [State(i, fsm.labels[i], fsm.π[i], fsm.T[i,end]) for i in 1:S]
end

function arcs(fsm::MatrixFSM{Tv}, state) where Tv
    retval = []
    for (i, v) in zip(findnz(fsm.T[state.id, 1:end-1])...)
        dest = State(i, fsm.labels[i], fsm.π[i], fsm.T[i,end])
        push!(retval, Arc{Tv}(dest, v))
    end
    retval
end

"""
    gpu(fsm)

Make a copy of `fsm` on GPU.
"""
function gpu(fsm::MatrixFSM{Tv}) where Tv
    T = CuSparseMatrixCSC(fsm.T)
    Tᵀ = CuSparseMatrixCSC(fsm.Tᵀ)
    C = CuSparseMatrixCSC(fsm.C)
    Cᵀ = CuSparseMatrixCSC(fsm.Cᵀ)
    return MatrixFSM{Tv}(
        CuSparseVector(fsm.π),
        CuSparseMatrixCSR(Tᵀ.colPtr, Tᵀ.rowVal, Tᵀ.nzVal, T.dims),
        CuSparseMatrixCSR(T.colPtr, T.rowVal, T.nzVal, Tᵀ.dims),
        CuSparseMatrixCSR(Cᵀ.colPtr, Cᵀ.rowVal, Cᵀ.nzVal, C.dims),
        CuSparseMatrixCSR(C.colPtr, C.rowVal, C.nzVal, Cᵀ.dims),
        fsm.labels
    )
end

function Base.convert(::Type{MatrixFSM{NT}}, mfsm) where NT <: Semiring
    MatrixFSM{NT}(
        copyto!(similar(mfsm.π, NT), mfsm.π),
        copyto!(similar(mfsm.T, NT), mfsm.T),
        copyto!(similar(mfsm.Tᵀ, NT), mfsm.Tᵀ),
        copyto!(similar(mfsm.C, NT), mfsm.C),
        copyto!(similar(mfsm.Cᵀ, NT), mfsm.Cᵀ),
        mfsm.labels
    )
end
Base.convert(::Type{MatrixFSM{T}}, fsm::MatrixFSM{T}) where T<:Semiring = fsm

struct UnionMatrixFSM{Tv} <: AbstractFSM{Tv}
    ranges::Vector
    π::AbstractSparseVector{Tv}
    T::AbstractSparseMatrix{Tv}
    Tᵀ::AbstractSparseMatrix{Tv}
    C::AbstractSparseMatrix{Tv}
    Cᵀ::AbstractSparseMatrix{Tv}
    labels::Vector
end

function states(fsm::UnionMatrixFSM)
    states = []
    for (ri, r) in enumerate(fsm.ranges)
        for i in r[1]:r[end]-1
            s = State((ri, i), fsm.labels[i], fsm.π[i], fsm.T[i,r[end]])
            push!(states, s)
        end
    end
    states
end

function arcs(fsm::UnionMatrixFSM{Tv}, state) where Tv
    retval = []
    ri, i = state.id
    r = fsm.ranges[ri]
    for (j, v) in zip(findnz(fsm.T[i, r[1]:r[end]-1])...)
        j = r[1]+j-1
        dest = State((ri, j), fsm.labels[j], fsm.π[j], fsm.T[j,r[end]])
        push!(retval, Arc{Tv}(dest, fsm.T[i,j]))
    end
    retval
end

function gpu(fsm::UnionMatrixFSM{Tv}) where Tv
    T = CuSparseMatrixCSC(fsm.T)
    Tᵀ = CuSparseMatrixCSC(fsm.Tᵀ)
    C = CuSparseMatrixCSC(fsm.C)
    Cᵀ = CuSparseMatrixCSC(fsm.Cᵀ)
    return UnionMatrixFSM{Tv}(
        fsm.ranges,
        CuSparseVector(fsm.π),
        CuSparseMatrixCSR(Tᵀ.colPtr, Tᵀ.rowVal, Tᵀ.nzVal, T.dims),
        CuSparseMatrixCSR(T.colPtr, T.rowVal, T.nzVal, Tᵀ.dims),
        CuSparseMatrixCSR(Cᵀ.colPtr, Cᵀ.rowVal, Cᵀ.nzVal, C.dims),
        CuSparseMatrixCSR(C.colPtr, C.rowVal, C.nzVal, Cᵀ.dims),
        fsm.labels
    )
end

function Base.union(fsms::Vararg{MatrixFSM{Tv},N}) where {Tv,N}
    nstates_sofar = 0
    ranges = []
    for fsm in fsms
        S = length(fsm.π)
        push!(ranges, (1:S) .+ nstates_sofar)
        nstates_sofar += S
    end

    UnionMatrixFSM{Tv}(
        ranges,
        vcat(map(x -> x.π, fsms)...),
        blockdiag(map(x -> x.T, fsms)...),
        blockdiag(map(x -> x.Tᵀ, fsms)...),
        blockdiag(map(x -> x.C, fsms)...),
        blockdiag(map(x -> x.Cᵀ, fsms)...),
        vcat(map(x -> x.labels, fsms)...)
    )
end

max_stars_repo_name: MartinKocour/HiddenMarkovModel

max_issues_repo_name: MartinKocour/HiddenMarkovModel

max_forks_repo_name: BUTSpeechFIT/ProbabilisticFSMs


text:
module Process

using Ito.Time
using Ito.TermStructure
using Calendar

#Stochastic processes
abstract StochasticProcess

evolve(sp::StochasticProcess, t0::Real, x0, dt::Real, dw) =
					apply_change(sp, apply_change(sp, x0 , drift_discretize(sp, t0, x0, dt)) , diffusion_discretize(sp, t0,x0,dt) * dw)
apply_change(ps::StochasticProcess, x, dx) = x + dx

covariance(sp::StochasticProcess, t0::Real, x0, dt::Real) = covariance_discretize(sp, t0, x0, dt)

#Euler discretization functions
drift_discretize(sp::StochasticProcess, t0::Real, x0, dt::Real) = drift(sp, t0, x0) .* dt

diffusion_discretize(sp::StochasticProcess, t0::Real, x0, dt::Real) = diffusion(sp, t0, x0) .* sqrt(dt)

function covariance_discretize(sp::StochasticProcess, t0::Real, x0, dt::Real)
	s= diffusion(sp, t0, x0)
	( s * s' ) .* dt
end

type GeometricBrownianMotion <: StochasticProcess
	start::Real
    mu::Real
    sigma::Real
end

diffusion(p::GeometricBrownianMotion, t, x) = p.sigma * x
drift(p::GeometricBrownianMotion, t, x) = p.mu * x
x0(p::GeometricBrownianMotion) = p.start

type GenericBlackScholesProcess
	start::Real
	dts::YieldTermStructure
	yts::YieldTermStructure
	vol::VolatilityTermStructure
end

apply_change(sp::GenericBlackScholesProcess, x, dx) = x + exp(dx)

diffusion(sp::GenericBlackScholesProcess, t, x) = vol(sp.vol, t, x)

function drift(sp::GenericBlackScholesProcess, t, x)
	t1 = t+.0001
    r=forward_rate(yts, t, t1, Continuous, NoFrequency)
    d=forward_rate(dts, t, t1, Continuous, NoFrequency)
    sigma = diffusion(sp, t, x)
    return r - d - 0.5 * sigma * sigma
end

end #module


text:
function main()
    s, t = readline(), readline()
    lₛ, lₜ = length(s), length(t)

    dp = zeros(Int, lₛ+1, lₜ+1)

    for i in 2:(lₛ+1), j in 2:(lₜ+1)
        dp[i, j] = max(dp[i-1, j], dp[i, j-1], dp[i-1, j-1] + (s[i-1]==t[j-1]))
    end

    ans = ""
    i, j = lₛ+1, lₜ+1
    while i ≥ 2 && j ≥ 2
        if s[i-1] == t[j-1]
            ans = "$(s[i-1])$ans"
            i -= 1
            j -= 1
        elseif dp[i, j-1] > dp[i-1, j]
            j -= 1
        else
            i -= 1
        end
    end

    println(ans)
end

main()


text:
include("CommonFuncs.jl")
threshold = 1e7 # σY ≈ 1e8

idx = 0
H_function = spd_Chol_Orth

if length(ARGS)==1
    global idx = parse(Int64, ARGS[1])
elseif length(ARGS)==2
    global idx = parse(Int64, ARGS[1])
    global tid = parse(Int64, ARGS[2])
else
    @warn("No ARGS provided")
end

@show H_function
if H_function==spd_zero_to_H || H_function==spd_Chol_Orth
    global nout = 4
end

if nntype=="ae_scaled"
    global nout = 3
end

if idx == 0
    global config=[20,20,20,nout]
elseif idx == 1
    global config=[20, 20, 20, 20, nout]
elseif idx == 2
    global config=[20, 20, 20, 20, 20, nout]
elseif idx == 3
    global config=[20, 20, 20, 20, 20, 20,nout]
elseif idx == 5
    global config=[nout]
end
printstyled("idx = $idx, config=$config, H_function=$H_function\n", color=:green)

function nn(ε, ε0, σ0) # ε, ε0, σ0 450x3
    local y, z
    global H0
    if nntype=="linear"
        y = ε*H0*stress_scale
        y
    elseif nntype=="ae_scaled"
        x = ε/strain_scale
        if isa(x, Array)
            x = constant(x)
        end
        # y = ae(x, [50,50,50,50,50,50,50,3], nntype)*stress_scale
        # y = ae(x, [50,50,50,50,50,50,50,3], nntype)*stress_scale
        y = ae(x, config, nntype)*stress_scale
    elseif nntype=="stiffmat"
        x = ε/strain_scale
        x = constant(x)

        σ0 = constant(σ0)

        y = ae(x, config, nntype)

        if H_function==spd_H
            z = spd_H(y, H0)
        else
            z = H_function(y)
        end

        σnn = squeeze(tf.matmul(z, tf.reshape((ε-ε0)/strain_scale, (-1,3,1))))

        ####################################
        # σH = (ε-ε0)/strain_scale * H0
        # s = σ0[:,1]^2-σ0[:,1]*σ0[:,2]+σ0[:,2]^2+3*σ0[:,3]^2

        # i = sigmoid(1000*(s-threshold)/1e9)
        # i = [i i i]
        # out = σnn .* i + σH .* (1-i)  + σ0/stress_scale
        ###############################################

        out = σnn + σ0/stress_scale

        out*stress_scale
    elseif nntype=="piecewise"
        x = [ε/strain_scale ε0/strain_scale σ0/stress_scale]
        x = constant(x)
        ε = constant(ε)
        ε0 = constant(ε0)
        σ0 = constant(σ0)

        y = ae(x, config, nntype)

        if H_function==spd_H
            z = spd_H(y, H0)
        else
            z = H_function(y)
        end

	    # op = tf.print(z)
	    # z = bind(z, op)
        # z = sym_H(y)

        σnn = squeeze(tf.matmul(z, tf.reshape((ε-ε0)/strain_scale, (-1,3,1))))
        σH = (ε-ε0)/strain_scale * H0
        s = σ0[:,1]^2-σ0[:,1]*σ0[:,2]+σ0[:,2]^2+3*σ0[:,3]^2

        i = sigmoid(1000*(s-threshold)/1e9)
        i = [i i i]
        out = σnn .* i + σH .* (1-i)  + σ0/stress_scale
        out*stress_scale
    elseif nntype=="doublenn"
        x = [ε/strain_scale ε0/strain_scale σ0/stress_scale]
        x = constant(x)
        ε = constant(ε)
        ε0 = constant(ε0)
        σ0 = constant(σ0)

        y = ae(x, config, nntype)
        i = sigmoid(ae(x, [20,20,1], nntype*"i"))
        # op = tf.print("call spd_H")
        if H_function==spd_H
            z = spd_H(y, H0)
        else
            z = H_function(y)
        end
        # z = bind(z, op)
        # z = sym_H(y)

        σnn = squeeze(tf.matmul(z, tf.reshape((ε-ε0)/strain_scale, (-1,3,1))))
        σH = (ε-ε0)/strain_scale * H0
        s = σ0[:,1]^2-σ0[:,1]*σ0[:,2]+σ0[:,2]^2+3*σ0[:,3]^2

        # i = sigmoid(1000*(s-threshold)/1e9)
        i = [i i i]
        out = σnn .* i + σH .* (1-i)  + σ0/stress_scale
        out*stress_scale
    else
        error("$nntype does not exist")
    end
end

function sigmoid_(z)

    return 1.0 / (1.0 + exp(-z))

end

function nn_helper(ε, ε0, σ0)
    local y1
    if nntype=="linear"
        x = reshape(reshape(ε,1,3)*H0,3,1)
    elseif nntype=="ae_scaled"
        x = reshape(ε/strain_scale,1, 3)
        reshape(nnae_scaled(x)*stress_scale,3,1)
    elseif nntype=="stiffmat"
        ε = ε/strain_scale
        ε0 = ε0/strain_scale
        σ0 = σ0/stress_scale
        x = reshape(ε, 1, 3)
        # y1 = (reshape(ε, 1, 3) - reshape(ε0, 1, 3))*sym_H(nnpiecewise(x))
        if H_function==spd_H
            y1 = (reshape(ε, 1, 3) - reshape(ε0, 1, 3))*spd_H(nnstiffmat(x), H0)
        else
            y1 = (reshape(ε, 1, 3) - reshape(ε0, 1, 3))*H_function(nnstiffmat(x))
        end
        y1 = reshape(y1, 3, 1)
        ################################################
        # y2 = reshape((reshape(ε, 1, 3) - reshape(ε0, 1, 3))*H0, 3,1)
        # s = σ0[1]^2-σ0[1]*σ0[2]+σ0[2]^2+3*σ0[3]^2
        # i = sigmoid_(1000*(s*stress_scale^2 - threshold)/1e9)
        # out = y1 * i + y2 * (1-i)  + reshape(σ0, 3, 1)
        #####################################################

        out = y1  + reshape(σ0, 3, 1)

        out*stress_scale
    elseif nntype=="piecewise"
        ε = ε/strain_scale
        ε0 = ε0/strain_scale
        σ0 = σ0/stress_scale
        x = reshape([ε;ε0;σ0],1, 9)
        # y1 = (reshape(ε, 1, 3) - reshape(ε0, 1, 3))*sym_H(nnpiecewise(x))
        if H_function==spd_H
            y1 = (reshape(ε, 1, 3) - reshape(ε0, 1, 3))*spd_H(nnpiecewise(x), H0)
        else
            y1 = (reshape(ε, 1, 3) - reshape(ε0, 1, 3))*H_function(nnpiecewise(x))
        end
        y1 = reshape(y1, 3, 1)
        y2 = reshape((reshape(ε, 1, 3) - reshape(ε0, 1, 3))*H0, 3,1)
        s = σ0[1]^2-σ0[1]*σ0[2]+σ0[2]^2+3*σ0[3]^2
        i = sigmoid_(1000*(s*stress_scale^2 - threshold)/1e9)
        out = y1 * i + y2 * (1-i)  + reshape(σ0, 3, 1)
        out*stress_scale
    elseif nntype=="doublenn"
        ε = ε/strain_scale
        ε0 = ε0/strain_scale
        σ0 = σ0/stress_scale
        x = reshape([ε;ε0;σ0],1, 9)
        # y1 = (reshape(ε, 1, 3) - reshape(ε0, 1, 3))*sym_H(nnpiecewise(x))
        if H_function==spd_H
            y1 = (reshape(ε, 1, 3) - reshape(ε0, 1, 3))*spd_H(nndoublenn(x), H0)
        else
            y1 = (reshape(ε, 1, 3) - reshape(ε0, 1, 3))*H_function(nndoublenn(x))
        end
        y1 = reshape(y1, 3, 1)
        y2 = reshape((reshape(ε, 1, 3) - reshape(ε0, 1, 3))*H0, 3,1)
        s = σ0[1]^2-σ0[1]*σ0[2]+σ0[2]^2+3*σ0[3]^2
        # i = sigmoid_(1000*(s*stress_scale^2 - threshold)/1e9)
        i = sigmoid_(nndoublenni(x)[1,1])
        out = y1 * i + y2 * (1-i)  + reshape(σ0, 3, 1)
        out*stress_scale
    elseif nntype=="stress"
        ε = ε/strain_scale
        ε0 = ε0/strain_scale
        σ0 = σ0/stress_scale
        x = reshape([ε;ε0;σ0],1, 9)
        out = nnstress(x)
        reshape(out*stress_scale,3,1)
    else
        error("$nntype does not exist")
    end
end

function post_nn(ε, ε0, σ0, Δt)
    f = x -> nn_helper(x, ε0, σ0)
    df = ForwardDiff.jacobian(f, ε)
    return f(ε), df
end


text:
struct UnorderedPair{T}
    a::T
    b::T
end

Base.hash(p::UnorderedPair, h::UInt) = hash(p.a, h) + hash(p.b, h)

function Base.:(==)(x::UnorderedPair, y::UnorderedPair)
    (x.a == y.a && x.b == y.b) || (x.a == y.b && x.b == y.a)
end


text:
module OmniSci

using Thrift, Random, Sockets, Distributed, DataFrames, Dates, DecFP, GeoInterface, LibGEOS, Decimals
import Thrift: process, meta, distribute, ThriftMeta, ThriftMetaAttribs
import Base: show, convert
import DataFrames: DataFrame

#### types and enums ####
export TDatum, TColumn, TExecuteMode

#### functions ####
export connect, disconnect
export get_status, get_hardware_info
export get_tables_meta, get_table_details
export get_users, get_databases
export sql_execute, sql_execute_df, sql_execute_gdf
export render_vega
export get_dashboards, get_dashboard_grantees
export load_table_binary_columnar, load_table_binary_arrow, load_table
export create_table, import_table
export get_roles, get_all_roles_for_user
export DataFrame

#Thrift interface code, mostly generated from Thrift.jl
include("common/common_types.jl")
include("completion_hints/completion_hints_types.jl")
include("serialized_result_set/serialized_result_set_types.jl")
include("extension_functions/extension_functions_types.jl")
include("mapd/mapd_types.jl")
include("mapd/MapD.jl")

#Package code
include("enums.jl") #hand-maintained code for more convenient enums
include("constructors.jl")
include("ipc.jl")
include("client.jl")
include("dataframe_show.jl")

end # module OmniSci


text:
using jInv.Mesh;
using ForwardHelmholtz
using ForwardHelmholtz.PointSourceADR
using Multigrid
using MAT
using EikonalInv

plotting = false;
if plotting
	using jInvVis
	using PyPlot;
	close("all");
end
include("getModels.jl");
function performExperiment(formulation::Int64,n,f,model,coarseSolveType = "MUMPS")
n_tup = tuple(n...);
dim = length(n);

pad    = 25;
if model=="linear"
	domain = [0.0,20.0,0.0,4.0];
	if dim == 3
		domain = [0.0,20.0,0.0,20.0,0.0,4.0];
	end
	m = getLinearModel(1.6,3.5,n);
	m = 1./(m.^2);
elseif model=="Marmousi"
	if dim==3
		error("Marmousi is only 2D");
	end
	domain = [0.0,17.0,0.0,3.5];
	model_dir = "./../../../BenchmarkModels/";
	file = matopen(string(model_dir,"ElasticMarmousi2_single.mat"));
	DICT = read(file); close(file);
	m = DICT["Vp"]';

	m = expandModelNearest(m,[13601;2801],n);
	m = 1./m.^2;
	# ms = copy(m);
	# for k=1:10
		# for j = 2:size(m,2)-1
			# for i = 2:size(m,1)-1
				# @inbounds ms[i,j] = (2*ms[i,j] + (ms[i-1,j-1]+ms[i-1,j]+ms[i-1,j+1]+ms[i,j-1]+ms[i,j+1]+ms[i+1,j-1]+ms[i+1,j]+ms[i,j+1]))/10.0;
			# end
		# end
	# end
	# m = ms;
elseif model=="Overthrust"
	domain = [0.0,20.0,0.0,20.0,0.0,4.65];
	model_dir = "./../../../BenchmarkModels/";
	file = matopen(string(model_dir,"3DOverthrust801801187.mat"));
	DICT = read(file); close(file);
	m = DICT["A"];DICT = 0;
	m = m.*1e-3;
	m = reshape(m,801,801,187);
	m = expandModelNearest(m,[801;801;187],n);
	m = 1./m.^2;

elseif model=="Wedge"
	if dim==3
		error("TBD");
	end
	n_cells = n-1; (m,Minv) = getWedgeModel(n);domain = Minv.domain;
else
	domain = [0.0,1.0,0.0,1.0];
	if dim == 3
		domain = [0.0,1.0,0.0,1.0,0.0,1.0];
	end
	m = ones(n_tup);
end

Minv = getRegularMesh(domain,n-1);
println("RUN: formulation = ",formulation," n = ",n," f = ",f," model = ",model," coarsest = ",coarseSolveType);

## Generating the right hand side
if dim==3
	q = zeros(Complex64,n_tup);
	src = [div(Minv.n[1]+1,2),div(Minv.n[2]+1,2),1];
	q[src[1],src[2],src[3]] = 1./(prod(Minv.h))
else
	q = zeros(Complex64,n_tup);
	src = [div(Minv.n[1]+1,2),1];
	q[src[1],src[2]] = 1./(prod(Minv.h))
end

q = q[:];

w = 2*pi*f

# if model!="const"
	# w = getMaximalFrequency(m,Minv)*0.95;
# end
println("maximal frequency is: f = ",getMaximalFrequency(m,Minv)/(2*pi));
println("sp*omega*h = ",w*maximum(Minv.h)*maximum(sqrt(m)));

if plotting
	figure()
	plotModel(1./sqrt(m),true,Minv,0);
	# return;
end

pad = pad*ones(Int64,3);
maxOmega = getMaximalFrequency(m,Minv);
ABLamp = maxOmega;
NeumannAtFirstDim = true;
gamma = 0.000001*2*pi*ones(Float32,size(m));

# (ADR_long,ADR_short,T,G1,G2) = getHelmholtzADR(Minv, m, w, gamma,NeumannAtFirstDim,pad,ABLamp,src,true);

# H,gamma = GetHelmholtzOperator(Minv, m, w, gamma,NeumannAtFirstDim,pad,ABLamp,true);

shift 		= 0.2;
Shift 		= GetHelmholtzShiftOP(m, w,shift)
levels      = 5;
numCores 	= 4;
maxIter     = 50;
relativeTol = 1e-5;
cycleType   ='K';
relaxType   = "Jac-GMRES";
relaxParam  = 0.8;
relaxPre 	= l->l+1;
relaxPost   = l->l+1;

if formulation == 1 # Pure Helmholtz
	# relaxPre 	= 2;
	# relaxPost   = 2;
	Shift = Shift';
	if eltype(q)==Complex64
		Shift = convert(SparseMatrixCSC{Complex64,Int64},Shift);
	end

	gamma += getABL(Minv,NeumannAtFirstDim,pad,ABLamp);
	######################################
	######################################
	# tic()
	# H_param = HelmholtzParam(Minv,gamma,m,w,NeumannAtFirstDim,true);
	# Hparam_shifted = getShiftedHelmholtzParam(H_param,shift);
	# function restrictParam(mesh_fine,mesh_coarse,param_fine,level)
		# P = getShiftedHelmholtzParam(param_fine,0.0);
		# P.gamma = restrictNodalVariables(P.gamma,mesh_fine.n+1);
		# P.m = restrictNodalVariables(P.m,mesh_fine.n+1);
		# P.Mesh = mesh_coarse;
		# return P;
	# end
	# getOperator = (mesh,Hparam)->GetHelmholtzOperator(mesh,Hparam.m,Hparam.omega,Hparam.gamma,Hparam.NeumannOnTop,Hparam.Sommerfeld);
	# MG = getMGparam(levels,numCores,maxIter,relativeTol,relaxType,relaxParam,
				# relaxPre,relaxPost,cycleType,coarseSolveType,0.5,0.0,"FullWeighting");
	# MGsetup(getMultilevelOperatorConstructor(Hparam_shifted,getOperator,restrictParam),Hparam_shifted.Mesh,MG,eltype(q),1,true);
	# SHT = MG.As[1];
	# set1 = toc();
	# set2 = 0.0;
	#########################################
	#########################################
	println("getting Helmholtz operator")
	tic()
	HT = GetHelmholtzOperator(Minv, m, w, gamma,NeumannAtFirstDim,pad,ABLamp,true)[1]';
	if eltype(q)==Complex64
		HT = convert(SparseMatrixCSC{Complex64,Int64},HT);
	end

	SHT = HT + Shift;HT = 0.0;
	set1 = toc()
	println("setup")
	tic()
	MG = getMGparam(levels,numCores,maxIter,relativeTol,relaxType,relaxParam,
				relaxPre,relaxPost,cycleType,coarseSolveType,0.5,0.0,"FullWeighting");
	MGsetup(SHT,Minv,MG,eltype(q),1,false);
	set2 = toc();
	#########################################
	#########################################

	HT = x->SHT'*x - Shift'*x;
	println("FGMRES-MG")
	tic()
	x = solveGMRES_MG(HT,MG,q,zeros(eltype(q),size(q)),true,5)[1]
	sol1 = toc()
	println(norm(HT(x) - q)/norm(q))

	println("Total setup = ", set1 + set2 );
	println("Total sol = ", sol1 );
elseif formulation == 4 # ADR long up
	println("problem disc.")
	gamma += getABL(Minv,NeumannAtFirstDim,pad,ABLamp);
	#############################################################
	###########################################################
	# tic()
	# T,G,LT = getADRcoefficients(Minv, m,src);
	# set1 = toc();
	# tic()
	# (ADR_longUp,~,~) = getHelmholtzADR(false,Minv, m, w, gamma,T,G,LT,NeumannAtFirstDim,src,true);
	# ADR_longUp = ADR_longUp';
	# ADRparam = getADRparam(Minv,gamma,m,w,T,G,LT,src,NeumannAtFirstDim,true);
	# getOperator = (mesh,ADRparam)->((ADR2,ADR1,~) = getHelmholtzADR(false,ADRparam.Minv, ADRparam.m, ADRparam.w, ADRparam.gamma,
									# ADRparam.T,ADRparam.G,ADRparam.LT,ADRparam.NeumannAtFirstDim,ADRparam.src,ADRparam.Sommerfeld);
									# return 0.25*ADR1 + 0.75*ADR2;);
	# MG = getMGparam(levels,numCores,maxIter,relativeTol,relaxType,relaxParam,
				# relaxPre,relaxPost,cycleType,coarseSolveType,0.5,0.0,"FullWeighting");

	# MGsetup(getMultilevelOperatorConstructor(ADRparam,getOperator,restrictADRParam),ADRparam.Minv,MG,eltype(q),1,true);
	# SHT = MG.As[1];
	# set2 = toc();

	###########################################################
	############################################################
	tic()
	T,G,LT = getADRcoefficients(Minv, m,src);
	set1 = toc();

	# tic()
	(ADR_longUp,ADR_short,~) = getHelmholtzADR(false,Minv, m, w, gamma,T,G,LT,NeumannAtFirstDim,src,true);
	T = 0; G = 0; LT = 0; m = 0; gamma = 0;
	# set1 = toc()
	println("setup")
	tic()
	ADR_longUp = ADR_longUp';
	if eltype(q)==Complex64
		ADR_short = convert(SparseMatrixCSC{Complex64,Int64},0.25*ADR_short') + convert(SparseMatrixCSC{Complex64,Int64},0.75*ADR_longUp);
	else
		ADR_short = 0.25*ADR_short' + 0.75*ADR_longUp;
	end

	MG = getMGparam(levels,numCores,maxIter,relativeTol,relaxType,relaxParam,
				relaxPre,relaxPost,cycleType,coarseSolveType,0.5,0.0,"FullWeighting");
	MGsetup(ADR_short,Minv,MG,eltype(q),1,false);
	set2 = toc()
	######################################################################
	#######################################################################3
	println("FGMRES-MG")
	tic()
	x = solveGMRES_MG(ADR_longUp,MG,q,zeros(eltype(q),size(q)),true,5)[1]
	sol1 = toc()
	println(norm(ADR_longUp'*x - q)/norm(q))
	println("Total setup = ", set1 + set2 );
	println("Total sol = ", sol1 );
elseif formulation == 5 # ADR long shifted as prec to ADR long
	# (ADR_long,ADR_short,T) = getHelmholtzADR(true,Minv, m, w, gamma,NeumannAtFirstDim,pad,ABLamp,src,true);
	# Shift = Shift'
	# M = exp(-1im*w*T);
	# ADR_long = (spdiagm(M)*ADR_long*spdiagm(1./M));
	# SHT = ADR_long' + Shift;ADR_long = 0.0;
	# MG = getMGparam(levels,numCores,maxIter,relativeTol,relaxType,relaxParam,
				# relaxPre,relaxPost,cycleType,coarseSolveType,0.5,0.0,"FullWeighting");
	# MGsetup(SHT,Minv,MG,eltype(q),1,true);
	# HT = x->SHT'*x - Shift'*x;
	# println("FGMRES-MG")
	# x = solveGMRES_MG(HT,MG,q,zeros(eltype(q),size(q)),true,5)[1]
	# println(norm(HT(x) - q)/norm(q))
elseif formulation == 8 # two-phase solution.
	# println("problem disc.")
	gamma += getABL(Minv,NeumannAtFirstDim,pad,ABLamp);
	# #############################################################
	# ###########################################################
	tic()
	T,G,LT = getADRcoefficients(Minv, m,src);
	set1 = toc();
	tic()
	ADRparam = getADRparam(Minv,gamma,m,w,T,G,LT,src,NeumannAtFirstDim,true);
	getOperator = (mesh,ADRparam)->((ADR2,ADR1,~) = getHelmholtzADR(true,ADRparam.Minv, ADRparam.m, ADRparam.w, ADRparam.gamma,
									ADRparam.T,ADRparam.G,ADRparam.LT,ADRparam.NeumannAtFirstDim,ADRparam.src,ADRparam.Sommerfeld);
									return ADR1);
	MG = getMGparam(levels,numCores,1,0.01,relaxType,relaxParam,
				relaxPre,relaxPost,cycleType,coarseSolveType,0.5,0.0,"FullWeighting");
	MGsetup(getMultilevelOperatorConstructor(ADRparam,getOperator,restrictADRParam),ADRparam.Minv,MG,eltype(q),1,true);
	ADRparam = 0;
	ADR_short = MG.As[1];
	set2 = toc();
	# #############################################################
	# ###########################################################
	# println("problem disc.")
	# tic()
	# (ADR_long,ADR_short,T) = getHelmholtzADR(true,Minv, m, w, gamma,NeumannAtFirstDim,pad,ABLamp,src,true);
	# set1 = toc()

	# # Lap = ForwardHelmholtz.getNodalLaplacianMatrix(Minv);
	# # HT = GetHelmholtzOperator(Minv, m, w, gamma,NeumannAtFirstDim,pad,ABLamp,true)[1]';
	# # HT = (spdiagm(1./M)*(HT)*spdiagm(M));
	# println("Setup:")
	# tic()
	# M = exp(-1im*w*T);
	# # ADR_short = (0.25*ADR_short+0.75*ADR_long)';
	# ADR_short = ADR_short';
	# ADR_long = ADR_long';
	# toc()
	# tic()
	# MG = getMGparam(levels,numCores,1,0.01,relaxType,relaxParam,
		# relaxPre,relaxPost,cycleType,coarseSolveType,0.5,0.0,"FullWeighting");
	# MGsetup(ADR_short,Minv,MG,eltype(q),1,false);
	# set2 = toc()
	# #############################################################
	# ###########################################################

	println("FGMRES-MG")
	tic()
	x = solveGMRES_MG(ADR_short,MG,q,zeros(eltype(q),size(q)),true,5)[1]
	sol1 = toc()
	# println(norm(ADR_long'*x - q)/norm(q))

	ADR_short = 0;
	MG.As[1] = MG.As[end];
	MG.As[2] = MG.As[end];
	MG.As[3] = MG.As[end];
	# x = zeros(eltype(q),size(q));

	# #############################################################
	# ###########################################################

	# HT = GetHelmholtzOperator(Minv, m, w, gamma,NeumannAtFirstDim,pad,ABLamp,true)[1]';
	# HT = (spdiagm(M)*(HT)*spdiagm(1./M));
	# Shift = Shift';
	# println("Setup II :")
	# tic()
	# rescaleMatrix(ADR_long,M); # HT = (spdiagm(M)*(ADR_long)*spdiagm(1./M));
	# HT = ADR_long;
	# x = x.*M;
	# relres = norm(HT'*x - q)/norm(q);
	# println("Rescaling the problem, res:", relres);
	# MG.relativeTol = (relativeTol/relres);
	# MG.maxOuterIter = maxIter;
	# SHT = HT + Shift;HT = 0.0;
	# replaceMatrixInHierarchy(MG,SHT);
	# HT = x->SHT'*x - Shift'*x;
	# set3 = toc()
	# #############################################################
	# ###########################################################

	(ADR_long,~,~) = getHelmholtzADR(true,Minv, m, w, gamma,T,G,LT,NeumannAtFirstDim,src,true);
	G = 0; LT = 0; m = 0; gamma = 0;
	ADR_long = ADR_long';
	tic()
	relres = norm(ADR_long'*x - q)/norm(q);
	println("Rescaling the problem, res:", relres);

	M = exp(-1im*w*T);
	if eltype(q)==Complex64
		M = convert(Array{Complex64},M);
	end
	T = 0;

	x = x.*M;
	rescaleMatrix(ADR_long,M);
	Shift = Shift';
	if eltype(q)==Complex64
		Shift = convert(SparseMatrixCSC{Complex64,Int64},Shift);
	end
	###############################################
	###############################################3
	# ADRparam = getADRparam(Minv,gamma,m,w,T,G,LT,src,NeumannAtFirstDim,true,M);
	# getOperator = (mesh,ADRparam)->((ADR2,ADR1,~) = getHelmholtzADR(true,ADRparam.Minv, ADRparam.m, ADRparam.w, ADRparam.gamma,
									# ADRparam.T,ADRparam.G,ADRparam.LT,ADRparam.NeumannAtFirstDim,ADRparam.src,ADRparam.Sommerfeld);
									# H = ADR2 + GetHelmholtzShiftOP(ADRparam.m, ADRparam.w,shift);
									# # M = exp(-1im*w*T);
									# # M = ADRparam.scaling;
									# # rescaleMatrix(H,M);
									# return H);

	# MG = getMGparam(levels,numCores,maxIter,relativeTol/relres,relaxType,relaxParam,
				# relaxPre,relaxPost,cycleType,coarseSolveType,0.5,0.0,"FullWeighting");
	# MGsetup(getMultilevelOperatorConstructor(ADRparam,getOperator,restrictADRParam),ADRparam.Minv,MG,eltype(q),1,true);
	# ADR_long = MG.As[1];
	# HT = x->MG.As[1]'*x - Shift'*x;
	###############################################
	###############################################
	MG.relativeTol = (relativeTol/relres);
	MG.maxOuterIter = maxIter;
	ADR_long_shifted = ADR_long + Shift;ADR_long = 0;
	replaceMatrixInHierarchy(MG,ADR_long_shifted);
	HT = x->MG.As[1]'*x - Shift'*x;

	set3 = toc();

	println("FGMRES-MG")
	tic()
	x = solveGMRES_MG(HT,MG,q,x,true,5)[1]
	sol2 = toc();

	println("Total setup = ", set1 + set2 + set3);
	println("Total sol = ", sol1 + sol2);
	println(norm(HT(x) - q)/norm(q))
end
println("**************************************** THE END ***************************************************");
end

if true
println("*********************************************************************")
println("********************** WARM-UP  ************************************")
println("*********************************************************************")
n = [65;65;33];
performExperiment(1,n,0.5,"linear","GMRES");
performExperiment(8,n,0.5,"linear","GMRES");
performExperiment(4,n,0.5,"linear","GMRES");

performExperiment(1,n,0.5,"Overthrust","GMRES");
performExperiment(8,n,0.5,"Overthrust","GMRES");
performExperiment(4,n,0.5,"Overthrust","GMRES");

# n = [257;257;65];
# performExperiment(1,n,1.5,"linear","GMRES");
# performExperiment(8,n,1.5,"linear","GMRES");
# performExperiment(4,n,1.5,"linear","GMRES");
# println("*********************************************************************")

# performExperiment(1,n,2.0,"linear","GMRES");
# performExperiment(8,n,2.0,"linear","GMRES");
# performExperiment(4,n,2.0,"linear","GMRES");

# n = [385;385;97];
# performExperiment(1,n,2.0,"linear","GMRES");
# performExperiment(8,n,2.0,"linear","GMRES");
# performExperiment(4,n,2.0,"linear","GMRES");
# println("*********************************************************************")

# performExperiment(1,n,3.0,"linear","GMRES");
# performExperiment(8,n,3.0,"linear","GMRES");
# performExperiment(4,n,3.0,"linear","GMRES");
# println("*********************************************************************")

 n = [513;513;129];
 performExperiment(1,n,3.0,"linear","GMRES");
 performExperiment(8,n,3.0,"linear","GMRES");
 performExperiment(4,n,3.0,"linear","GMRES");
 println("*********************************************************************")

 performExperiment(1,n,4.0,"linear","GMRES");
 performExperiment(8,n,4.0,"linear","GMRES");
 performExperiment(4,n,4.0,"linear","GMRES");
 println("*********************************************************************")

# n = [769;769;193];
# performExperiment(1,n,4.0,"linear","GMRES");
# performExperiment(8,n,4.0,"linear","GMRES");
# performExperiment(4,n,4.0,"linear","GMRES");
# println("*********************************************************************")

# performExperiment(1,n,6.0,"linear","GMRES");
# performExperiment(8,n,6.0,"linear","GMRES");
# performExperiment(4,n,6.0,"linear","GMRES");
# println("*********************************************************************")

n = [257;257;65];
performExperiment(1,n,2.0,"Overthrust","GMRES");
performExperiment(8,n,2.0,"Overthrust","GMRES");
performExperiment(4,n,2.0,"Overthrust","GMRES");
println("*********************************************************************")

performExperiment(1,n,3.0,"Overthrust","GMRES");
performExperiment(8,n,3.0,"Overthrust","GMRES");
performExperiment(4,n,3.0,"Overthrust","GMRES");

n = [385;385;97];
performExperiment(1,n,3.0,"Overthrust","GMRES");
performExperiment(8,n,3.0,"Overthrust","GMRES");
performExperiment(4,n,3.0,"Overthrust","GMRES");
println("*********************************************************************")

performExperiment(1,n,4.0,"Overthrust","GMRES");
performExperiment(8,n,4.0,"Overthrust","GMRES");
performExperiment(4,n,4.0,"Overthrust","GMRES");
println("*********************************************************************")

n = [513;513;129];
performExperiment(1,n,4.0,"Overthrust","GMRES");
performExperiment(8,n,4.0,"Overthrust","GMRES");
performExperiment(4,n,4.0,"Overthrust","GMRES");
println("*********************************************************************")

performExperiment(1,n,6.0,"Overthrust","GMRES");
performExperiment(8,n,6.0,"Overthrust","GMRES");
performExperiment(4,n,6.0,"Overthrust","GMRES");
println("*********************************************************************")

end
# 2D results:

if false

println("*********************************************************************")
println("********************** WARM-UP  ************************************")
println("*********************************************************************")

n = [257;129];
performExperiment(1,n,2.0,"linear","GMRES");
performExperiment(8,n,2.0,"linear","GMRES");
performExperiment(4,n,2.0,"linear","GMRES");

println("*********************************************************************")
println("*********************************************************************")
println("*********************************************************************")

n = [769;257];
performExperiment(1,n,3.5,"linear","GMRES");
performExperiment(8,n,3.5,"linear","GMRES");
performExperiment(4,n,3.5,"linear","GMRES");
println("*********************************************************************")
println("*********************************************************************")
println("*********************************************************************")

performExperiment(1,n,5.5,"linear","GMRES");
performExperiment(8,n,5.5,"linear","GMRES");
performExperiment(4,n,5.5,"linear","GMRES");

println("*********************************************************************")
println("*********************************************************************")
println("*********************************************************************")

n = [1025;385];
performExperiment(1,n,5.5,"linear","GMRES");
performExperiment(8,n,5.5,"linear","GMRES");
performExperiment(4,n,5.5,"linear","GMRES");

println("*********************************************************************")
println("*********************************************************************")
println("*********************************************************************")

performExperiment(1,n,7.5,"linear","GMRES");
performExperiment(8,n,7.5,"linear","GMRES");
performExperiment(4,n,7.5,"linear","GMRES");

# println("*********************************************************************")
# println("*********************************************************************")
# println("*********************************************************************")

n = [1537;513];
performExperiment(1,n,7.5,"linear","GMRES");
performExperiment(8,n,7.5,"linear","GMRES");
performExperiment(4,n,7.5,"linear","GMRES");

performExperiment(1,n,11.0,"linear","GMRES");
performExperiment(8,n,11.0,"linear","GMRES");
# performExperiment(9,n,11.0,"linear","GMRES");
performExperiment(4,n,11.0,"linear","GMRES");

println("*********************************************************************")
println("*********************************************************************")
println("*********************************************************************")

n = [769;257];
performExperiment(1,n,3.0,"Marmousi","GMRES");
performExperiment(8,n,3.0,"Marmousi","GMRES");
performExperiment(4,n,3.0,"Marmousi","GMRES");

println("*********************************************************************")
println("*********************************************************************")
println("*********************************************************************")

performExperiment(1,n,4.5,"Marmousi","GMRES");
performExperiment(8,n,4.5,"Marmousi","GMRES");
performExperiment(4,n,4.5,"Marmousi","GMRES");

println("*********************************************************************")
println("*********************************************************************")
println("*********************************************************************")

n = [1025;385];
performExperiment(1,n,4.5,"Marmousi","GMRES");
performExperiment(8,n,4.5,"Marmousi","GMRES");
performExperiment(4,n,4.5,"Marmousi","GMRES");

println("*********************************************************************")
println("*********************************************************************")
println("*********************************************************************")

performExperiment(1,n,6.5,"Marmousi","GMRES");
performExperiment(8,n,6.5,"Marmousi","GMRES");
performExperiment(4,n,6.5,"Marmousi","GMRES");

println("*********************************************************************")
println("*********************************************************************")
println("*********************************************************************")

n = [1537;513];
performExperiment(1,n,6.5,"Marmousi","GMRES");
performExperiment(8,n,6.5,"Marmousi","GMRES");
performExperiment(4,n,6.5,"Marmousi","GMRES");

# println("*********************************************************************")
# println("*********************************************************************")
# println("*********************************************************************")

performExperiment(1,n,9.0,"Marmousi","GMRES");
performExperiment(8,n,9.0,"Marmousi","GMRES");
performExperiment(4,n,9.0,"Marmousi","GMRES");
end