text:
using ImageCore, Colors, FixedPointNumbers, OffsetArrays, Test

if VERSION >= v"1.2.0-DEV.229"
    sumsz(img) = Base.dims2string(size(img)) * ' '
else
    sumsz(img) = ""
end

# N0f8 is shown as either N0f8 or Normed{UInt8, 8}
# RGB is shown as ColorTypes.RGB or RGB
function typestring(::Type{T}) where T
    buf = IOBuffer()
    show(buf, T)
    String(take!(buf))
end
N0f8_str = typestring(N0f8)
N0f16_str = typestring(N0f16)
RGB_str = typestring(RGB)

@testset "show" begin
    rgb32 = rand(RGB{Float32}, 3, 5)
    v = view(rgb32, 2:3, :)
    @test summary(v) == "2×5 view(::Array{RGB{Float32},2}, 2:3, :) with eltype $(RGB_str){Float32}"
    a = channelview(rgb32)
    @test summary(a) == "3×3×5 reinterpret(Float32, ::Array{RGB{Float32},3})"
    num64 = rand(3,5)
    b = colorview(RGB, num64)
    @test summary(b) == "5-element reshape(reinterpret(RGB{Float64}, ::$(typeof(num64))), 5) with eltype $(RGB_str){Float64}"
    rgb8 = rand(RGB{N0f8}, 3, 5)
    c = rawview(channelview(rgb8))
    @test summary(c) == "3×3×5 rawview(reinterpret(N0f8, ::Array{RGB{N0f8},3})) with eltype UInt8"
    @test summary(rgb8) == "3×5 Array{RGB{N0f8},2} with eltype $(RGB_str){$(N0f8_str)}"
    rand8 = rand(UInt8, 3, 5)
    d = normedview(permuteddimsview(rand8, (2,1)))
    @test summary(d) == "5×3 normedview(N0f8, PermutedDimsArray(::$(typeof(rand8)), (2, 1))) with eltype $(N0f8_str)"
    e = permuteddimsview(normedview(rand8), (2,1))
    @test summary(e) == "5×3 PermutedDimsArray(reinterpret(N0f8, ::$(typeof(rand8))), (2, 1)) with eltype $(N0f8_str)"
    rand16 = rand(UInt16, 3, 5)
    f = permuteddimsview(normedview(N0f16, rand16), (2,1))
    @test summary(f) == "5×3 PermutedDimsArray(reinterpret(N0f16, ::$(typeof(rand16))), (2, 1)) with eltype $(N0f16_str)"
    g = channelview(rgb8)
    @test summary(g) == "3×3×5 reinterpret(N0f8, ::Array{RGB{N0f8},3})"
    h = OffsetArray(rgb8, -1:1, -2:2)
    @test summary(h) == "$(sumsz(h))OffsetArray(::Array{RGB{N0f8},2}, -1:1, -2:2) with eltype $(RGB_str){$(N0f8_str)} with indices -1:1×-2:2"
    i = channelview(h)
    @test summary(i) == "$(sumsz(i))reinterpret(N0f8, OffsetArray(::Array{RGB{N0f8},3}, 1:1, -1:1, -2:2)) with indices 1:3×-1:1×-2:2"
    c = channelview(rand(RGB{N0f8}, 2))
    o = OffsetArray(c, -1:1, 0:1)
    @test summary(o) == "$(sumsz(o))OffsetArray(reinterpret(N0f8, ::Array{RGB{N0f8},2}), -1:1, 0:1) with eltype $(N0f8_str) with indices -1:1×0:1"
    # Issue #45
    a = collect(tuple())
    @test summary(a) == "0-element $(typeof(a))"
    b = view(a, :)
    @test summary(b) == "0-element view(::$(typeof(a)), :) with eltype Union{}"
end

nothing


text:
abstract type AbstractTPInterpolation <: SatPureAproximation end

"""
    invTlogPInterpolation(max_length::Int = 256)

performs interpolation of available values of saturation equilibria.

"""
struct invTlogPInterpolation
    max_length::Int
end

function sizehint(cachedmodel::CACHED_SAT_PURE_APROX{invTlogPInterpolation})
    return cachedmodel.sat_pure_aprox.max_length
end

const TPInterpolation = invTlogPInterpolation

TPInterpolation() = TPInterpolation(256) #max length of 200 data points

function x0_sat_pure(model::CACHED_SAT_PURE_APROX{invTlogPInterpolation},T)
    if !haskey(model.cache,:sat_pure)
        return x0_sat_pure(model.model,T)
    end

    T_vec,P_vec,Vl_vec,Vv_vec = model.cache[:sat_pure]
    isone(length(T_vec)) && return x0_sat_pure(model.model,T) #use underlying x0_sat_pure
    Tmin,Tc = first(T_vec),last(T_vec)
    nan = zero(T)/zero(T)
    T > Tc && return [nan,nan] #error, because T>Tc
    T < Tmin && return x0_sat_pure(model.model,T) #use underlying x0_sat_pure
    pos = searchsorted(T_vec,T)

    if length(pos) == 1 #the value is already there
        log10vl = log10(only(Vl_vec[pos]))
        log10vv = log10(only(Vv_vec[pos]))
        return [log10vl,log10vv]
    end
    adyacent_pos = last(pos):first(pos)
    T1,T2 = T_vec[adyacent_pos]
    P1,P2 = P_vec[adyacent_pos]

    #invTlogP: log(Pi)  ≈ A-B/Ti
    #log(P1) - log(P2) = -B(1/T1 - 1/T2)
    #log(P1/P2)/(1/T1 - 1/T2)
    B = -log(P1/P2)/(1/T1 - 1/T2)
    A = log(P1) + B/T1

    P0 = exp(A - B/T)
    vv0 = volume_virial(model.model,P0,T)

    vl1,vl2 = Vl_vec[adyacent_pos]
    #rackett approximation
    #log(V1) = C + D(1-Tr1)^(2/7) = C+D*fT1
    #log(V1) - log(V2) = D(ft1-ft2)
    fT1 = (1-T1/Tc)^(2/7)
    fT2 = (1-T2/Tc)^(2/7)
    D = log(vl1/vl2)/(fT1-fT2)
    C = log(vl1) - D*fT1
    vl0_interp = exp(C+D*(1-T/Tc)^(2/7))
    vl0_lb = x0_volume(model.model,P0,T,SA[1.0],phase=:l)
    #vl00 = 0.9*vl0_interp + 0.1*vl0_lb
    #vl0 = volume_compress(model.model,P0,T,V0=vl00)
    vl0 = vl0_interp
    return [log10(vl0),log10(vv0)]
end

function sat_pure_approx(model::CACHED_SAT_PURE_APROX{invTlogPInterpolation},T)
    if !haskey(model.cache,:sat_pure)
        return sat_pure(model,T)
    end
    T_vec,P_vec,Vl_vec,Vv_vec = model.cache[:sat_pure]
    isone(length(T_vec)) && return sat_pure(model,T) #use underlying x0_sat_pure
    Tmin,Tc = first(T_vec),last(T_vec)
    nan = zero(T)/zero(T)
    T > Tc && return (nan,nan,nan) #error, because T>Tc
    T < Tmin && return sat_pure(model,T) #use underlying x0_sat_pure
    pos = searchsorted(T_vec,T)

    if length(pos) == 1 #the value is already there
        vl = only(Vl_vec[pos])
        vv = only(Vv_vec[pos])
        p0 = only(P_vec[pos])
        return p0,vl,vv
    end
    adyacent_pos = last(pos):first(pos)
    T1,T2 = T_vec[adyacent_pos]
    P1,P2 = P_vec[adyacent_pos]

    #invTlogP: log(Pi)  ≈ A-B/Ti
    #log(P1) - log(P2) = -B(1/T1 - 1/T2)
    #log(P1/P2)/(1/T1 - 1/T2)
    B = -log(P1/P2)/(1/T1 - 1/T2)
    A = log(P1) + B/T1

    P0 = exp(A - B/T)
    vv0 = volume_virial(model.model,P0,T)

    vl1,vl2 = Vl_vec[adyacent_pos]
    #rackett approximation
    #log(V1) = C + D(1-Tr1)^(2/7) = C+D*fT1
    #log(V1) - log(V2) = D(ft1-ft2)
    fT1 = (1-T1/Tc)^(2/7)
    fT2 = (1-T2/Tc)^(2/7)
    D = log(vl1/vl2)/(fT1-fT2)
    C = log(vl1) - D*fT1
    vl0_interp = exp(C+D*(1-T/Tc)^(2/7))
    vl0_lb = x0_volume(model.model,P0,T,SA[1.0],phase=:l)
    vl0 = vl0_interp
    #vl0 = volume_compress(model.model,P0,T,V0=vl00)
    return (P0,vl0,vv0)
end
function sat_pure_cached(::TPInterpolation,model::CachedEoS,T)
    res = sat_pure!(model,T)
end

export sat_pure_approx, TPInterpolation

ext: jl


text:
@compat abstract type AbstractAWSType end

import Base: show

type AWSError
    code::String
    msg::String
end
function AWSError(pd::LightXML.XMLElement)
    code = ""
    msg = ""
    epd = LightXML.find_element(pd, "Error")
    if epd !== nothing
        elem = LightXML.find_element(epd, "Code")
        (elem === nothing) || (code = LightXML.content(elem))
        elem = LightXML.find_element(epd, "Message")
        (elem === nothing) || (msg = LightXML.content(elem))
    end
    AWSError(code, msg)
end
export AWSError

#show(io::IO, o::AWSError) = print("code: $(o.code), msg : $(o.msg), $(o.request_id)")

type AWSResponse
    http_code::Int
    headers::Dict{AbstractString,AbstractString}
    body::Union{String, Void}
    pd::Union{LightXML.XMLElement, Void}
    obj::Any

    AWSResponse() = new(0, Dict{AbstractString, AbstractString}(), "", nothing, nothing)
end
export AWSResponse

const BASIC_TYPES = Union{String, Integer, AbstractFloat, DateTime}
is_basic_type{T<:BASIC_TYPES}(v::T) = true
is_basic_type(v) = false

_ename(name::Symbol) = _ename(string(name))
function _ename(name::String)
    name = lcfirst(name)
    if name in ("type", "return", "end")
        name = "_" * name
    end
    name
end
tbasename{T<:AbstractAWSType}(obj::T) = tbasename(T)
tbasename{T}(::Type{T}) = rsplit(string(T), "."; limit=2)[end]

function aws_execute(env::AWSEnv, action::String, api_name::String, api_ver::String, params::Vector{Tuple}=Vector{Tuple}())
    # prepare the standard params
    push!(params, ("Action", action))
    push!(params, ("AWSAccessKeyId", env.aws_id))
    push!(params, ("Timestamp", get_utc_timestamp()))
    push!(params, ("Version", api_ver))

    amz_headers, data, signed_querystr = canonicalize_and_sign(env, api_name, false, params)

    ep_path = env.ep_path
    endswith(ep_path, '/') && (ep_path = ep_path * "/")
    complete_url = "https://" * ep_host(env, api_name) * ep_path * "?" * signed_querystr
    (env.dbg || env.dry_run) && info("URL:\n$complete_url")

	headers = Dict{String, String}(k=>v for (k,v) in amz_headers)

    apiresp = AWSResponse()
    env.dry_run && (return apiresp)

    # make the request
    resp = Requests.get(complete_url; headers=headers, timeout=env.timeout)

    apiresp.http_code = resp.status
    apiresp.headers = resp.headers
    apiresp.body = String(copy(resp.data))

    env.dbg && info("HTTPCode: ", apiresp.http_code, "\nHeaders: ", apiresp.headers, "\nBody : ", apiresp.body)

    if (resp.status >= 200) && (resp.status <= 299)
        for (n,v) in resp.headers
            if lowercase(n) == "content-type"
                if contains(v, "/xml")
                    apiresp.pd = LightXML.root(LightXML.parse_string(apiresp.body))
                end
                break
            end
        end
    elseif (resp.status >= 400) && (resp.status <= 599)
        if length(apiresp.body) > 0
            apiresp.obj = AWSError(LightXML.root(LightXML.parse_string(apiresp.body)))
        else
            error("HTTP error: $(resp.status)")
        end
    else
        error("HTTP error: $(resp.status), $(apiresp.body)")
    end

    apiresp
end

function _parse_from_xml{T<:AbstractAWSType}(::Type{T}, types::Dict{Symbol,Vector{Pair{Symbol,DataType}}}, pd::Union{LightXML.XMLElement,Void}; member_pfx="")
    typename = Symbol(tbasename(T))
    typeprops = types[typename]
    _parse_from_xml(T, types, typeprops, pd; member_pfx=member_pfx)
end

function _parse_from_xml{T<:AbstractAWSType}(::Type{T}, types::Dict{Symbol,Vector{Pair{Symbol,DataType}}}, typeprops::Vector{Pair{Symbol,DataType}}, pd::Union{LightXML.XMLElement,Void}; member_pfx="")
    kwargs = Dict{Symbol,Any}()
    if pd !== nothing
        pd_type = LightXML.find_element(pd, tbasename(T))
        (pd_type === nothing) && (pd_type = pd)
        for (elemname,elemtype) in typeprops
            ename = string(elemname)
            # handle field names that match julia reserved types
            startswith(ename, "_") && (ename = ename[2:end])
            if ((elemtype <: Dict) || ((elemtype <: Vector) && !(elemtype <: Vector{UInt8}))) && isempty(member_pfx)
                elem = LightXML.get_elements_by_tagname(pd_type, ename)
            else
                elem = LightXML.find_element(pd_type, ename)
            end
            (elem === nothing) && continue
            kwargs[Symbol(_ename(ename))] = _parse_from_xml(elemtype, types, elem; member_pfx=member_pfx)
        end
    end
    T(; kwargs...)
end

_parse_from_xml{T<:BASIC_TYPES}(::Type{T}, types::Dict{Symbol,Vector{Pair{Symbol,DataType}}}, pd::LightXML.XMLElement; member_pfx="") = AWS.safe_parse_as(T, LightXML.content(pd))

function _parse_from_xml{T}(::Type{Vector{T}}, types::Dict{Symbol,Vector{Pair{Symbol,DataType}}}, pd_list::Union{LightXML.XMLElement,Void}; member_pfx="")
    result = Vector{T}()
    if pd_list !== nothing
        for pd in LightXML.child_elements(pd_list)
            push!(result, _parse_from_xml(T, types, pd; member_pfx=member_pfx))
        end
    end
    result
end

function _parse_from_xml{T}(::Type{Vector{T}}, types::Dict{Symbol,Vector{Pair{Symbol,DataType}}}, pd_list::Union{Vector{LightXML.XMLElement},Void}; member_pfx="")
    result = Vector{T}()
    if pd_list !== nothing
        for pd in pd_list
            push!(result, _parse_from_xml(T, types, pd; member_pfx=member_pfx))
        end
    end
    result
end

function _parse_from_xml{Tn,Tv}(::Type{Dict{Tn,Tv}}, types::Dict{Symbol,Vector{Pair{Symbol,DataType}}}, pd_list::Union{Vector{LightXML.XMLElement},Void}; member_pfx="")
    result = Dict{Tn,Tv}()
    for elem in pd_list
        name = _parse_from_xml(Tn, types, LightXML.find_element(elem, "Name"); member_pfx=member_pfx)
        value = _parse_from_xml(Tv, types, LightXML.find_element(elem, "Value"); member_pfx=member_pfx)
        result[name] = value
    end
    result
end

function _prep_params{T<:AbstractAWSType}(obj::T, params::Vector{Tuple}=Vector{Tuple}(), pfx::String=""; member_pfx="")
    for m in fieldnames(T)
        fld_val = getfield(obj, m)
        (fld_val == nothing) && continue

        arg_name = string(m)
        # handle field names that match julia reserved types
        startswith(arg_name, "_") && (arg_name = arg_name[2:end])
        # captitalize the first letter for the argument
        arg_name = ucfirst(arg_name)

        if is_basic_type(fld_val)
            push!(params, (pfx * arg_name, aws_string(fld_val)))
        elseif isa(fld_val, Vector{UInt8})
            push!(params, (pfx * arg_name, aws_string(fld_val)))
        elseif isa(fld_val, Array)
            for (idx, fv) in enumerate(fld_val)
                subarg_name = "$(pfx)$(arg_name).$(member_pfx)$(idx)"
                if is_basic_type(fv)
                    push!(params, (subarg_name, aws_string(fv)))
                else
                    _prep_params(fv, params, subarg_name*".")
                end
            end
        elseif isa(fld_val, Dict)
            for (idx, nv) in enumerate(fld_val)
                subarg_name = "$(pfx)$(arg_name).$(member_pfx)$(idx).Name"
                subarg_value = "$(pfx)$(arg_name).$(member_pfx)$(idx).Value"
                push!(params, (subarg_name, aws_string(nv[1])))
                fv = nv[2]
                if is_basic_type(fv)
                    push!(params, (subarg_value, aws_string(fv)))
                else
                    _prep_params(fv, params, subarg_value*".")
                end
            end
        else # compound type
            _prep_params(fld_val, params, "$(pfx)$(arg_name).")
        end
    end
    params
end

function generate_type(outfile::IO, specname::String, typename::String, typeprops::Vector{Pair{Symbol,DataType}})
    println(outfile, "")
    println(outfile, "type $typename <: AbstractAWSType")
    for (elemname,elemtype) in typeprops
        println(outfile, "    $(_ename(elemname))::Union{Void,$(elemtype)}")
    end

    elemkws = ["$(_ename(elemname))=nothing" for (elemname,elemtype) in typeprops]
    println(outfile, "    function $(typename)(;", join(elemkws, ", "), ")")
    elemparams = [_ename(elemname) for (elemname,elemtype) in typeprops]
    println(outfile, "        new(", join(elemparams, ", "), ")")
    println(outfile, "    end")

    println(outfile, "    $(typename)(pd) = parse_from_xml($(typename), $(specname), pd)")
    println(outfile, "end # $typename")
    println(outfile, "export $typename")
end

function generate_types(filename::String, specname::String, types::Vector{Pair{Symbol,Vector{Pair{Symbol,DataType}}}})
    open(filename, "w") do outfile
        println(outfile, "# generated from $specname via AWS/src/codegen.jl")
        println(outfile, "# do not edit, edit specification file and regenerate instead")
        for (typename,typeprops) in types
            generate_type(outfile, specname, string(typename), typeprops)
        end
    end
end

function generate_apis(filename::String, specname::String, apis::Dict)
    open(filename, "w") do outfile
        println(outfile, "# generated from $specname via AWS/src/codegen.jl")
        println(outfile, "# do not edit, edit specification file and regenerate instead")

        for (apiname,apiprops) in apis
            for (reqresptype,reqrespspec) in apiprops
                generate_type(outfile, specname, "$(apiname)$(reqresptype)", reqrespspec)
            end

            has_req = :Request in keys(apiprops)
            has_resp = :Result in keys(apiprops)
            resptype = specname[1:end-8] * "Response"
            apicall = lowercase(specname[1:end-8]) * "_execute"
            println(outfile, "")
            println(outfile, "function $(apiname)(env::AWSEnv", has_req ? ", msg::$(apiname)Request)" : ")")
            println(outfile, "    resp::$(resptype) = $(apicall)(env, \"$(apiname)\"", has_req ? ", msg)" : ")")
            if has_resp
                println(outfile, "    if (resp.pd !== nothing) && (resp.obj === nothing)")
                println(outfile, "        resp.obj = $(apiname)Result(resp.pd)")
                println(outfile, "    end")
            end
            println(outfile, "    resp")
            println(outfile, "end")
            if has_req
                println(outfile, "$(apiname)(env::AWSEnv; kwargs...) = $(apiname)(env, $(apiname)Request(; kwargs...))")
            end
            println(outfile, "export $(apiname)")
        end
    end
end


text:
# Black pixels
test_ellipse_1 = [
0 0 0 0 0 ;
0 0 0 0 0 ;
0 1 1 1 0 ;
0 0 0 0 0 ;
0 0 0 0 0]

test_ellipse_2 = [
0 0 0 0 0 ;
0 0 0 1 0 ;
0 1 1 0 0 ;
0 0 0 0 0 ;
0 0 0 0 0]

test_ellipse_3 = [
0 0 0 0 0 ;
0 0 0 1 0 ;
0 0 1 0 0 ;
0 1 0 0 0 ;
0 0 0 0 0]

test_ellipse_4 = [
0 0 0 0 0 ;
0 0 1 0 0 ;
0 0 1 0 0 ;
0 1 0 0 0 ;
0 0 0 0 0]

test_ellipse_5 = [
0 0 0 0 0 ;
0 0 1 0 0 ;
0 0 1 0 0 ;
0 0 1 0 0 ;
0 0 0 0 0]

test_ellipse_6 = [
0 0 0 0 0 ;
0 1 0 0 0 ;
0 0 1 0 0 ;
0 0 1 0 0 ;
0 0 0 0 0]

test_ellipse_7 = [
0 0 0 0 0 ;
0 1 0 0 0 ;
0 0 1 0 0 ;
0 0 0 1 0 ;
0 0 0 0 0]

test_ellipse_8 = [
0 0 0 0 0 ;
0 0 0 0 0 ;
0 1 1 0 0 ;
0 0 0 1 0 ;
0 0 0 0 0]

test_ellipse_9 = [
0 0 0 0 0 ;
0 0 0 0 0 ;
0 1 1 1 0 ;
0 0 0 0 0 ;
0 0 0 0 0]

test_ellipse_10 = [
0 0 1 0 0 ;
0 1 1 1 0 ;
1 1 1 1 1 ;
0 1 1 1 0 ;
0 0 1 0 0]

test_ellipse_11 = [
0 0 0 0 0 ;
0 1 1 1 0 ;
1 1 1 1 1 ;
0 1 1 1 0 ;
0 0 0 0 0]

test_ellipse_12 = [
0 0 0 0 0 ;
0 0 1 0 0 ;
1 1 1 1 1 ;
0 0 1 0 0 ;
0 0 0 0 0]

test_ellipse_13 = [
0 0 0 0 0 ;
0 0 1 0 0 ;
0 1 1 1 0 ;
0 0 1 0 0 ;
0 0 0 0 0]

test_ellipse_14 = [
0 0 0 0 0 ;
0 0 0 0 0 ;
0 1 1 1 0 ;
0 0 0 0 0 ;
0 0 0 0 0]

#
# test_ellipse_15 = [
# 0 0 0 0 0 ;
# 0 0 0 0 0 ;
# 0 1 1 0 0 ;
# 0 0 0 1 0 ;
# 0 0 0 0 0]


text:
"""
   Immutable structure that contains recordings. Only pressure fields are sampled.

# Fields

-`nt::Ti`: Number of samples per trace.
-`nr::Ti`: Number of receivers.
-`dt::Tv`: Time sampling rate.
-`irz::Vector{Ti}`: Vertical grid index for receivers.
-`irx::Vector{Ti}`: Horizontal grid index for receivers.
-`spt2rec::Vector{Ti}`:Index mapping of receiver to snapshot.
-`p::Vector{Tv}`: Vector of recodrings of preassure field.
"""
struct Recordings{Ti<:Int64, Tv<:AbstractFloat}
     nt      :: Ti
     nr      :: Ti
     dt      :: Tv
     irz     :: Vector{Ti}  # vertical grid index of receivers
     irx     :: Vector{Ti}  # horizontal grid index of receivers
     spt2rec :: Vector{Ti}  # index mapping of receiver to snapshot
     p       :: Matrix{Tv}  # recordings of pressure field
     vz      :: Matrix{Tv}  # recordings of vertical velocity

end

"""
    R = Recordings(rz, rx, params; <keyword arguments>)

Constructs a Recordings structure.

Please see the extended help with `?? Source`

# Arguments

-`rz`: Receiver coordinate.
-`rx`: Receiver coordinate.
-`params::TdParams`: Parameters for time domain modelling.

# Keyword Arguments:

-`location_flag`: Receiver(s) positioning by `"index"` or `"distance".

# Extended help

      R=Recordings(rz, rx, params; location_flag="index")

The structure `Recordings{Ti<:Int64, Tv<:AbstractFloat}` defines the receiver location by either its position or relative index in the model through the `location_flag`.

"""
function Recordings(rz::Vector, rx::Vector,
         params::TdParams; location_flag="index")

    # number of receivers
    nr  = length(rz)
    nr == length(rx) || error("length of receiver coordinate doesn't match")

    irz = zeros(Int64, nr)
    irx = zeros(Int64, nr)

    # location are provided as index
    if location_flag == "index"
       for i = 1 : nr
           irz[i] = round(Int64, rz[i])
           irx[i] = round(Int64, rx[i])
       end

    # location are given as distance
    elseif location_flag == "distance"
       for i = 1 : nr
           irz[i] = round(Int64, rz[i]/params.dz) + 1
           irx[i] = round(Int64, rx[i]/params.dx) + 1
       end

    else
       error("wrong specification of receiver location")
    end

    # error checking
    for i = 1 : nr
        if irz[i] > params.nz || irx[i] > params.nx || irz[i] < 1 || irx[i] < 1
           error("receiver located outside of modeling area")
        end
    end

    # can't put receivers on the free surface
    if params.free_surface
       for i = 1 : nr
           if irz[i] == 1
              error("can't put receiver at free surface, move it deeper")
           end
       end
    end

    # the auxillary vector mapping snapshot to recordings
    spt2rec = zeros(Int64, nr)
    for i = 1 : nr
        spt2rec[i] = (irx[i]+params.npml-1) * params.Nz + irz[i] + params.ntop
    end

    # initalize empty recordings
    p   = zeros(params.data_format, params.nt, nr)
    vz  = zeros(params.data_format, params.nt, nr)
    rec = Recordings(params.nt, nr, params.dt, irz, irx, spt2rec, p, vz)

    return rec
end

"""
   write recordings into su file
"""
function write_recordings(path::String, rec::Recordings)

    fid = open(path, "w")

    # write number of samples
    write(fid, rec.nt)

    # write number of receivers
    write(fid, rec.nr)

    # write the locations of receivers (index)
    write(fid, rec.irz)
    write(fid, rec.irx)

    # write the index mapping vector
    write(fid, rec.spt2rec)

    # data format code
    if eltype(rec.p) == Float32
       write(fid, 1)             # 1 indicate Float32

    elseif eltype(rec.p) == Float64
       write(fid, 2)             # 2 indicate Float64

    else
       error("non-support data format")
    end

    # write time sampling interval
    write(fid, rec.dt)

    # write the recordings
    write(fid, rec.p)
    write(fid, rec.vz)

    #close file
    close(fid)

    return nothing
end

"""
   read one shot gather
"""
function read_recordings(path::String)

    # open the file
    fid = open(path, "r")

    # number of samples
    nt  = read(fid, Int64)

    # number of traces
    nr  = read(fid, Int64)

    # read the index locations of receivers
    irz = zeros(Int64, nr); read!(fid, irz)
    irx = zeros(Int64, nr); read!(fid, irx)

    # read the index mapping vector
    spt2rec = zeros(Int64, nr); read!(fid, spt2rec)

    # read the data_format code
    code = read(fid, Int64)
    if code == 1
       data_format = Float32
    elseif code == 2
       data_format = Float64
    else
       error("non-support data format")
    end

    # read the time sampling interval
    dt = read(fid, data_format)

    # read the recordings (p and vz)
    p = zeros(data_format, nt * nr)
    read!(fid, p); p = reshape(p, nt, nr);

    vz = zeros(data_format, nt * nr)
    read!(fid, vz); vz = reshape(vz, nt, nr);

    #close file
    close(fid)

    # organize them into a structure
    rec = Recordings(nt, nr, dt, irz, irx, spt2rec, p, vz)

    return rec
end

"""
   sampling one snashot to fill the recordings
"""
function sample_spt2rec!(rec::Recordings, spt::Snapshot, it::Int64)

    # loop over receivers
    for i = 1 : rec.nr
        idx = rec.spt2rec[i]
        rec.p[it,i]  = spt.pz[idx] + spt.px[idx]
        rec.vz[it,i] = spt.vz[idx]
    end

    return nothing
end

"""
   the adjoint of sampling operator, inject recordings to snapshot
"""
function inject_rec2spt!(spt::Snapshot, rec::Recordings, it::Int64)

    # loop over receivers
    for i = 1 : rec.nr
        idx = rec.spt2rec[i]
        spt.pz[idx] = spt.pz[idx] + rec.p[it,i]
        spt.px[idx] = spt.px[idx] + rec.p[it,i]
    end

    return nothing
end

"""
   check whether two shot gathers are equal to each other
"""
function recordings_isequal(rec1::Tr, rec2::Tr) where {Tr<:Recordings}

    rec1.nt      != rec2.nt      && (return false)
    rec1.nr      != rec2.nr      && (return false)
    rec1.dt      != rec2.dt      && (return false)
    rec1.irz     != rec2.irz     && (return false)
    rec1.irx     != rec2.irx     && (return false)
    rec1.spt2rec != rec2.spt2rec && (return false)
    norm(rec1.p-rec2.p)/norm(rec1.p) > 10*eps(eltype(rec1.p)) && (return false)

    return true
end


text:
function OMMParser()
    return OMMParser(())
end

function OMMParser(arg0::DataContext)
    return OMMParser((DataContext,), arg0)
end

function parse(obj::OMMParser, arg0::InputStream)
    return jcall(obj, "parse", OMMFile, (InputStream,), arg0)
end

function parse(obj::OMMParser, arg0::InputStream, arg1::JString)
    return jcall(obj, "parse", OMMFile, (InputStream, JString), arg0, arg1)
end

function parse(obj::OMMParser, arg0::JString)
    return jcall(obj, "parse", OMMFile, (JString,), arg0)
end

function with_conventions(obj::OMMParser, arg0::IERSConventions)
    return jcall(obj, "withConventions", OMMParser, (IERSConventions,), arg0)
end

function with_data_context(obj::OMMParser, arg0::DataContext)
    return jcall(obj, "withDataContext", OMMParser, (DataContext,), arg0)
end

function with_international_designator(obj::OMMParser, arg0::jint, arg1::jint, arg2::JString)
    return jcall(obj, "withInternationalDesignator", OMMParser, (jint, jint, JString), arg0, arg1, arg2)
end

function with_mission_reference_date(obj::OMMParser, arg0::AbsoluteDate)
    return jcall(obj, "withMissionReferenceDate", OMMParser, (AbsoluteDate,), arg0)
end

function with_mu(obj::OMMParser, arg0::jdouble)
    return jcall(obj, "withMu", OMMParser, (jdouble,), arg0)
end

function with_simple_eop(obj::OMMParser, arg0::jboolean)
    return jcall(obj, "withSimpleEOP", OMMParser, (jboolean,), arg0)
end


text:
"""
    Beta(α,β)

The *Beta distribution* has probability density function

```math
f(x; \\alpha, \\beta) = \\frac{1}{B(\\alpha, \\beta)}
 x^{\\alpha - 1} (1 - x)^{\\beta - 1}, \\quad x \\in [0, 1]
```

The Beta distribution is related to the [`Gamma`](@ref) distribution via the
property that if ``X \\sim \\operatorname{Gamma}(\\alpha)`` and ``Y \\sim \\operatorname{Gamma}(\\beta)``
independently, then ``X / (X + Y) \\sim Beta(\\alpha, \\beta)``.

```julia
Beta()        # equivalent to Beta(1, 1)
Beta(a)       # equivalent to Beta(a, a)
Beta(a, b)    # Beta distribution with shape parameters a and b

params(d)     # Get the parameters, i.e. (a, b)
```

External links

* [Beta distribution on Wikipedia](http://en.wikipedia.org/wiki/Beta_distribution)

"""
struct Beta{T<:Real} <: ContinuousUnivariateDistribution
    α::T
    β::T
    Beta{T}(α::T, β::T) where {T} = new{T}(α, β)
end

function Beta(α::T, β::T) where {T<:Real}
    @check_args(Beta, α > zero(α) && β > zero(β))
    return Beta{T}(α, β)
end

function Beta(α::T, β::T, ::NoArgCheck) where {T<:Real}
    return Beta{T}(α, β)
end

Beta(α::Real, β::Real) = Beta(promote(α, β)...)
Beta(α::Integer, β::Integer) = Beta(float(α), float(β))
Beta(α::Real) = Beta(α, α)
Beta() = Beta(1.0, 1.0, NoArgCheck())

@distr_support Beta 0.0 1.0

#### Conversions
function convert(::Type{Beta{T}}, α::Real, β::Real) where T<:Real
    Beta(T(α), T(β))
end
function convert(::Type{Beta{T}}, d::Beta{S}) where {T <: Real, S <: Real}
    Beta(T(d.α), T(d.β), NoArgCheck())
end

#### Parameters

params(d::Beta) = (d.α, d.β)
@inline partype(d::Beta{T}) where {T<:Real} = T

#### Statistics

mean(d::Beta) = ((α, β) = params(d); α / (α + β))

function mode(d::Beta)
    (α, β) = params(d)
    (α > 1 && β > 1) || error("mode is defined only when α > 1 and β > 1.")
    return (α - 1) / (α + β - 2)
end

function mode(d::Beta, ::NoArgCheck)
    (α, β) = params(d)
    return (α - 1) / (α + β - 2)
end

modes(d::Beta) = [mode(d)]

function var(d::Beta)
    (α, β) = params(d)
    s = α + β
    return (α * β) / (abs2(s) * (s + 1))
end

meanlogx(d::Beta) = ((α, β) = params(d); digamma(α) - digamma(α + β))

varlogx(d::Beta) = ((α, β) = params(d); trigamma(α) - trigamma(α + β))
stdlogx(d::Beta) = sqrt(varlogx(d))

function skewness(d::Beta)
    (α, β) = params(d)
    if α == β
        return zero(α)
    else
        s = α + β
        (2(β - α) * sqrt(s + 1)) / ((s + 2) * sqrt(α * β))
    end
end

function kurtosis(d::Beta)
    α, β = params(d)
    s = α + β
    p = α * β
    6(abs2(α - β) * (s + 1) - p * (s + 2)) / (p * (s + 2) * (s + 3))
end

function entropy(d::Beta)
    α, β = params(d)
    s = α + β
    lbeta(α, β) - (α - 1) * digamma(α) - (β - 1) * digamma(β) +
        (s - 2) * digamma(s)
end

#### Evaluation

@_delegate_statsfuns Beta beta α β

gradlogpdf(d::Beta{T}, x::Real) where {T<:Real} =
    ((α, β) = params(d); 0 <= x <= 1 ? (α - 1) / x - (β - 1) / (1 - x) : zero(T))

#### Sampling

struct BetaSampler{T<:Real, S1 <: Sampleable{Univariate,Continuous},
                   S2 <: Sampleable{Univariate,Continuous}} <:
    Sampleable{Univariate,Continuous}
    γ::Bool
    iα::T
    iβ::T
    s1::S1
    s2::S2
end

function sampler(d::Beta{T}) where T
    (α, β) = params(d)
    if (α ≤ 1.0) && (β ≤ 1.0)
        return BetaSampler(false, inv(α), inv(β),
                           sampler(Uniform()), sampler(Uniform()))
    else
        return BetaSampler(true, inv(α), inv(β),
                           sampler(Gamma(α, one(T))),
                           sampler(Gamma(β, one(T))))
    end
end

# From Knuth
function rand(rng::AbstractRNG, s::BetaSampler)
    if s.γ
        g1 = rand(rng, s.s1)
        g2 = rand(rng, s.s2)
        return g1 / (g1 + g2)
    else
        iα = s.iα
        iβ = s.iβ
        while true
            u = rand(rng) # the Uniform sampler just calls rand()
            v = rand(rng)
            x = u^iα
            y = v^iβ
            if x + y ≤ one(x)
                if (x + y > 0)
                    return x / (x + y)
                else
                    logX = log(u) * iα
                    logY = log(v) * iβ
                    logM = logX > logY ? logX : logY
                    logX -= logM
                    logY -= logM
                    return exp(logX - log(exp(logX) + exp(logY)))
                end
            end
        end
    end
end

function rand(rng::AbstractRNG, d::Beta{T}) where T
    (α, β) = params(d)
    if (α ≤ 1.0) && (β ≤ 1.0)
        while true
            u = rand(rng)
            v = rand(rng)
            x = u^inv(α)
            y = v^inv(β)
            if x + y ≤ one(x)
                if (x + y > 0)
                    return x / (x + y)
                else
                    logX = log(u) / α
                    logY = log(v) / β
                    logM = logX > logY ? logX : logY
                    logX -= logM
                    logY -= logM
                    return exp(logX - log(exp(logX) + exp(logY)))
                end
            end
        end
    else
        g1 = rand(rng, Gamma(α, one(T)))
        g2 = rand(rng, Gamma(β, one(T)))
        return g1 / (g1 + g2)
    end
end

#### Fit model

# TODO: add MLE method (should be similar to Dirichlet)

"""
    fit(::Type{<:Beta}, x::AbstractArray{T})

fit a `Beta` distribution
"""
function fit(::Type{<:Beta}, x::AbstractArray{T}) where T<:Real
    x_bar = mean(x)
    v_bar = varm(x, x_bar)
    temp = ((x_bar * (one(T) - x_bar)) / v_bar) - one(T)
    α = x_bar * temp
    β = (one(T) - x_bar) * temp
    return Beta(α, β)
end


text:
#=
test_hello_world:
- Julia version: 1.5.0
- Author: shisa
- Date: 2020-08-05
=#

module TestHelloWorld
    # packages
    using Test

    # external modules
    include("hello_world.jl")

    # methods
    function test()
        @testset "HelloWorld" begin
            @test_nowarn HelloWorld.main()
        end
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    using .TestHelloWorld
    TestHelloWorld.test()
end

max_stars_repo_name: ShisatoYano/JuliaPractice

max_issues_repo_name: ShisatoYano/JuliaPractice

max_forks_repo_name: ShisatoYano/JuliaPractice


text:
export ExponentialDist, nc

"""
    ExponentialDist{FT} <: AbstractSizeDist{FT}

Exponential size distribution defined by two moments: total water content, L̅ and
mean droplet mass, x̅
"""
struct ExponentialDist{FT} <: AbstractSizeDist{FT}
    L̅::FT # Total water content, kg/m3
    x̅::FT # Mean initial droplet mass, kg

    # TODO: Add special constructor to enforce positive args
    function ExponentialDist(;L̅::FT, x̅::FT) where {FT <: Real}
        new{FT}(L̅, x̅)
    end
end

function Base.show(io::IO, d::ExponentialDist{FT}) where {FT <: Real}
    print(io, "ExponentalDist(L̅=$(d.L̅ * 1e3) g/m³, x̅=$(d.x̅) kg)")
end

"""
    nc(dist::ExponentialDist{FT}, x::FT) where {FT <: Real}

Compute the number concentration for droplets of the given mass for a
specific size distribution.
"""
function nc(dist::ExponentialDist{FT}, x::FT) where {FT <: Real}
    return (dist.L̅ / dist.x̅^2) * exp(-x / dist.x̅)
end


text:
# ----------------------------------------
# FUNCTIONS
# ----------------------------------------

function distance_row(left, right; deletion = 1, insertion = 1, substitution = 1)
    function complete(left, right)
        previous = Array{Int64}(undef, length(right) + 1)
        current  = Array{Int64}(undef, length(right) + 1)

        for r ∈ 1:length(right) + 1
            @inbounds previous[r] = (r - 1) * insertion
        end

        for l ∈ 1:length(left)
            @inbounds current[1] = l * deletion

            for r ∈ 1:length(right)
                @inbounds current[r + 1] = min(previous[r + 1] + deletion,
                                               current[r] + insertion,
                                               previous[r] + (left[l] == right[r] ? 0 : substitution))
            end

            ( previous, current ) = ( current, previous )
        end

        return previous
    end

    if length(left) == 0
        return length(right)

    elseif length(right) == 0
        return length(left)

    else
        return complete(left, right)[end]
    end
end

max_stars_repo_name: jlegare/LevenshteinToolkit

max_issues_repo_name: jlegare/LevenshteinToolkit

max_forks_repo_name: jlegare/LevenshteinToolkit


text:
function myAwesomeFunction()
    println("hello there!");
    end;

myAwesomeFunction();

function sum(a::Int, b::Int)::Float16
    summation::Int = a + b;
    return summation;
end

s = sum(1,2);
println(s);

max_stars_repo_name: fahimfarhan/legendary-coding-odyssey

max_issues_repo_name: fahimfarhan/legendary-coding-odyssey

max_forks_repo_name: fahimfarhan/legendary-coding-odyssey


text:
# Hellwarth 1999 PRB - Part IV; T-dep of the Feynman variation parameter
# A Friday afternoon of hacking to try and implement the T-dep electron-phonon coupling from the above PRB
# Which was unusually successful! And more or less reproduced Table III

# one-dimensional numerical integration in Julia using adaptive Gauss-Kronrod quadrature
using QuadGK

# Equation numbers follow above Hellwarth 1999 PRB
# 62b
A(v,w,β)=3/β*( log(v/w) - 1/2*log(2*π*β) - log(sinh(v*β/2)/sinh(w*β/2)))

# 62d
Y(x,v,β)=1/(1-exp(-v*β))*(1+exp(-v*β)-exp(-v*x)-exp(v*(x-β)))
# 62c integrand
f(x,v,w,β)=(exp(β-x)+exp(x))/(w^2*x*(1-x/β)+Y(x,v,β)*(v^2-w^2)/v)^(1/2)
# 62c
B(v,w,β,α) = α*v/(sqrt(π)*(exp(β)-1)) * quadgk(x->f(x,v,w,β),0,β/2)[1]
#62e
C(v,w,β)=3/4*(v^2-w^2)/v * (coth(v*β/2)-2/(v*β))

F(v,w,β,α)=-(A(v,w,β)+B(v,w,β,α)+C(v,w,β)) #(62a)

# Can now evaluate, e.g.
# F(v,w,β,α)=F(7.2,6.5,1.0,1.0)
# BUT - this is just the objective function! Not the optimised parameters.
# Also there's a scary numeric integration (quadgk) buried within...

"Print out F(alpha,beta) for a specific v,w; as a test"
function test_fns()
    @printf("\t\t")
    for α in 1:5
        @printf("α=%d\t\t",α)
    end
    @printf("\n")

    for β in 1:0.25:3.0
        v=w=4
        print("β: $β  \t||")
        for α in 1:5
            @printf("%f\t",F(v,w,β,α))
        end
        println()
    end
end

test_fns()  # OK - very primitive!

"
These are 1D traces along the solution for Alpha=Beta=1 in Helwarth PRB TABLE III,
this was used to correct a transcription error in the above typed-in equations
It was also good to see what F(v,w) looked like as a function of v and w near an optimal solution"
function test_trace()
	v=7.20
	w=6.5
	α=1.0
	β=1.0

	for v=6:0.1:8
    	@printf("%f %f\n",v,F(v,w,β,α))
	end

	@printf("\n")
	v=7.20
	for w=6:0.1:7
    	@printf("%f %f\n",w,F(v,w,β,α))
	end
end
test_trace()

# Angle for the ringside seats, when the fall, don't blame me, Bring on the Major Leagues
using Optim
# Julia package stuffed full of magic, does auto-differentation & etc. etc.

Fopt(x) = F(x[1],x[2],1,1)

function test_Fopt()
    show(Fopt([7.2,6.5]))
# OK! It looks like I can bury the alpha, beta parameters (which we don't optimise), by wrapping our function in a function definition.
    initial=[7.2,6.5]

    show(optimize(Fopt,  initial, LBFGS()))

    show(optimize(Fopt, initial, BFGS(), Optim.Options(autodiff=true)))
end

test_Fopt()

function test_Optim()
	# After a bit of fiddling, I figured out how to add bounds, to stop that 'DomainError',
	# which occurs where the you are evaluating log(-ve Real), i.e. w<0.0 or v<0.0

	initial=[7.2,6.5]

	lower=[0.0,0.0]
	upper=[10.0,10.0]

	@printf("\t\t")
	for α in 1:5
		@printf("α=%d\t\t",α)
	end
	@printf("\n")

	for β in 1:0.25:3.0
		print("β: $β  \t||")
		for α in 1:5
			myf(x) = F(x[1],x[2],β,α)
			solution=optimize(DifferentiableFunction(myf), initial, lower, upper, Fminbox(); optimizer = ConjugateGradient, optimizer_o=Optim.Options(autodiff=true))
			minimum=Optim.minimizer(solution)

			v=minimum[1]
			w=minimum[2]
			#print(solution,"\t")
			@printf("%.2f %.2f\t",v,w)

		end
		println()
	end
end

test_Optim()

function test_Optimisers()
	# So that looks really good! I was super stoked to see how close these values are to TABLE III in Hellwarth
	# However, the solutions all start on (7.20,6.50) so that top-left data point is cheating, whereas the
	# others have some disagreement / noise associated with them
	# I was wondering whether it might be a function of the optimiser, so thought I'd try them all

	initial=[7.1,6.5]
	# Main use of these bounds is stopping v or w going negative, at which you get a NaN error as you are evaluating log(-ve Real)
	lower=[1.0,1.0]
	upper=[10.0,10.0]

	for optimizer in [BFGS, LBFGS, ConjugateGradient] # Newton, GradientDescent, NelderMead - steps outside box & log(-ve)->NaN error
		@printf("\n\t\t##### NOW TRIALING: %s #####\n\n",optimizer)

		@printf("\t\t")
		for α in 1:5
			@printf("α=%d\t\t",α)
		end
		@printf("\n")

		for β in 1:0.25:3.0
			print("β: $β  \t||")
			for α in 1:5
				myf(x) = F(x[1],x[2],β,α)
				res=optimize(DifferentiableFunction(myf), initial, lower, upper, Fminbox(); optimizer = optimizer, optimizer_o=Optim.Options(autodiff=true))
				minimum=Optim.minimizer(res)
				#show(Optim.converged(res)) # All came out as 'true'

				#print(solution,"\t")
				@printf("%.2f %.2f\t",minimum[1],minimum[2])
			end
			println()
		end
	end
end

test_Optimisers()


text:
module PNUTS

include("example.jl")

end # module


text:
# Test the example scripts

module TestSIMPExample
println("PointLoadCantilever Example")
mktempdir() do dir
    cd(dir) do
        include(joinpath(@__DIR__, "../../docs/src/literate/simp.jl"))
    end
end
end

module TestBESOExample
println("BESO Example")
mktempdir() do dir
    cd(dir) do
        include(joinpath(@__DIR__, "../../docs/src/literate/beso.jl"))
    end
end
end

module TestGESOExample
println("GESO Example")
mktempdir() do dir
    cd(dir) do
        include(joinpath(@__DIR__, "../../docs/src/literate/geso.jl"))
    end
end
end


text:
using MPI
using StaticArrays
using CLIMA
using CLIMA.VariableTemplates
using CLIMA.Mesh.Topologies
using CLIMA.Mesh.Grids
using CLIMA.MPIStateArrays
using CLIMA.DGmethods
using CLIMA.DGmethods.NumericalFluxes
using Printf
using LinearAlgebra
using Logging
using GPUifyLoops

@static if haspkg("CuArrays")
  using CUDAdrv
  using CUDAnative
  using CuArrays
  CuArrays.allowscalar(false)
  const ArrayTypes = (CuArray, )
else
  const ArrayTypes = (Array, )
end

import CLIMA.DGmethods: BalanceLaw, vars_aux, vars_state, vars_gradient,
                        vars_diffusive, vars_integrals, integrate_aux!,
                        flux_nondiffusive!, flux_diffusive!, source!, wavespeed,
                        update_aux!, indefinite_stack_integral!,
                        reverse_indefinite_stack_integral!,  boundary_state!,
                        init_aux!, init_state!, init_ode_state,
                        LocalGeometry

struct IntegralTestModel{dim} <: BalanceLaw
end

vars_integrals(::IntegralTestModel, T) = @vars(a::T,b::T)
vars_aux(m::IntegralTestModel,T) = @vars(int::vars_integrals(m,T),
                                         rev_int::vars_integrals(m,T),
                                         coord::SVector{3,T}, a::T, b::T)

vars_state(::IntegralTestModel, T) = @vars()
vars_diffusive(::IntegralTestModel, T) = @vars()

flux_nondiffusive!(::IntegralTestModel, _...) = nothing
flux_diffusive!(::IntegralTestModel, _...) = nothing
source!(::IntegralTestModel, _...) = nothing
boundary_state!(_, ::IntegralTestModel, _...) = nothing
init_state!(::IntegralTestModel, _...) = nothing
wavespeed(::IntegralTestModel,_...) = 1

function init_aux!(::IntegralTestModel{dim}, aux::Vars,
                   g::LocalGeometry) where {dim}
  x,y,z = aux.coord = g.coord
  if dim == 2
    aux.a = x*y + z*y
    aux.b = 2*x*y + sin(x)*y^2/2 - (z-1)^2*y^3/3
  else
    aux.a = x*z + z^2/2
    aux.b = 2*x*z + sin(x)*y*z - (1+(z-1)^3)*y^2/3
  end
end

function update_aux!(dg::DGModel, m::IntegralTestModel, Q::MPIStateArray,
                     auxstate::MPIStateArray, t::Real)
  indefinite_stack_integral!(dg, m, Q, auxstate, t)
  reverse_indefinite_stack_integral!(dg, m, auxstate, t)
end

@inline function integrate_aux!(m::IntegralTestModel, integrand::Vars,
                                state::Vars, aux::Vars)
  x,y,z = aux.coord
  integrand.a = x + z
  integrand.b = 2*x + sin(x)*y - (z-1)^2*y^2
end

using Test
function run(mpicomm, dim, ArrayType, Ne, N, FT)

  brickrange = ntuple(j->range(FT(0); length=Ne[j]+1, stop=3), dim)
  topl = StackedBrickTopology(mpicomm, brickrange,
                              periodicity=ntuple(j->true, dim))

  grid = DiscontinuousSpectralElementGrid(topl,
                                          FloatType = FT,
                                          DeviceArray = ArrayType,
                                          polynomialorder = N,
                                         )
  dg = DGModel(IntegralTestModel{dim}(),
               grid,
               Rusanov(),
               CentralNumericalFluxDiffusive(),
               CentralGradPenalty())

  Q = init_ode_state(dg, FT(0))
  dQdt = similar(Q)

  dg(dQdt, Q, nothing, 0.0)

  # Wrapping in Array ensure both GPU and CPU code use same approx
  @test Array(dg.auxstate.data[:, 1, :]) ≈ Array(dg.auxstate.data[:, 8, :])
  @test Array(dg.auxstate.data[:, 2, :]) ≈ Array(dg.auxstate.data[:, 9, :])
end

let
  MPI.Initialized() || MPI.Init()

  mpicomm = MPI.COMM_WORLD
  ll = uppercase(get(ENV, "JULIA_LOG_LEVEL", "INFO"))
  loglevel = ll == "DEBUG" ? Logging.Debug :
  ll == "WARN"  ? Logging.Warn  :
  ll == "ERROR" ? Logging.Error : Logging.Info
  logger_stream = MPI.Comm_rank(mpicomm) == 0 ? stderr : devnull
  global_logger(ConsoleLogger(logger_stream, loglevel))
  @static if haspkg("CUDAnative")
    device!(MPI.Comm_rank(mpicomm) % length(devices()))
  end

  numelem = (5, 5, 5)
  lvls = 1

  polynomialorder = 4

  @testset "$(@__FILE__)" for ArrayType in ArrayTypes
    for FT in (Float64,) #Float32)
      for dim = 2:3
        err = zeros(FT, lvls)
        for l = 1:lvls
          @info (ArrayType, FT, dim)
          run(mpicomm, dim, ArrayType, ntuple(j->2^(l-1) * numelem[j], dim),
              polynomialorder, FT)
        end
      end
    end
  end
end

nothing

max_stars_repo_name: akshaysridhar/CLIMA

max_issues_repo_name: akshaysridhar/CLIMA

max_forks_repo_name: akshaysridhar/CLIMA


text:
function add(a, b)
   try
       a + b
   catch
       println("caught exception")
       a * b
   end
end

println(add(2, 6))
println(add(1//2, 1//2))
println(add("Hello ", "world"))

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
println("Building LiterateOrg.jl")

# Bootstrap package by extracting all package code to a Julia file and
# include it.

src_file = joinpath(dirname(@__FILE__), "..", "src", "LiterateOrg.org")

deps_dir = joinpath(dirname(@__FILE__), "..", "deps")
bootstrap_file = joinpath(deps_dir, "bootstrap.jl")

start_code_pat = r"[ ]*#\+begin_src[ ]+julia(.*)"
end_code_pat = r"[ ]*#\+end_src"

println("Bootstrapping LiterateOrg.jl to $(bootstrap_file)")

code_mode = false

open(src_file) do infile
    open(bootstrap_file, "w") do outfile
        for line in readlines(infile)
            global code_mode
            if occursin(start_code_pat, lowercase(line))
                code_mode = true
                continue
            elseif occursin(end_code_pat, lowercase(line))
                code_mode = false
                continue
            end
            code_mode && write(outfile, "$(line)\n")
        end
    end
end

println("Running bootstrap file")

# This is necessary for included test expressions to work.
if VERSION.major > 0 || VERSION.minor > 6
    using Test
else
    using Base.Test
end
include(bootstrap_file)
rm(bootstrap_file)

println("Tangling LiterateOrg.jl")

tangle_package(src_file, "LiterateOrg")


text:
module ImmersedMirrors

export Grid, coarse_grid, function_to_grid, operator_to_grid
export Mirror, distance_to_mirror, mirror_image, smoothstep
export Extrapolator, make_dirichlet

using Distributed, RecipesBase
@everywhere using LinearAlgebra, SparseArrays, ProgressMeter, ForwardDiff

include("./mirror.jl")
include("./grid.jl")
include("./extrapolator.jl")

end


text:
using Interpolations
using Trapz
using Unitful

abstract type AbstractFilter{T} <: AbstractVector{T} end

abstract type DetectorType end
struct Photon <: DetectorType end
struct Energy <: DetectorType end

function Base.parse(::Type{DetectorType}, s::AbstractString)
    tok = lowercase(s)
    if tok === "photon"
        return Photon()
    elseif tok === "energy"
        return Energy()
    else
        throw(ArgumentError("failed to parse DetectorType from $s"))
    end
end

struct PhotometricFilter{T,WT,DT<:DetectorType,TT<:AbstractVector{T},ST<:Union{String,Nothing},ET} <: AbstractFilter{T}
    wave::WT
    throughput::TT
    detector::DT
    name::ST
    etp::ET
end

function PhotometricFilter(wave::AbstractVector, throughput::AbstractVector{T}; detector=Photon(), name=nothing) where T
    length(wave) == length(throughput) || error("wavelength and throughput arrays must match")
    bc = zero(T)
    etp = LinearInterpolation(ustrip(wave), throughput; extrapolation_bc=bc)
    return PhotometricFilter(wave, throughput, detector, name, etp)
end

Base.getindex(f::PhotometricFilter, idx...) = throughput(f)[idx...]

Base.show(io::IO, f::PhotometricFilter) = print(io, f.name)

function Base.show(io::IO, ::MIME"text/plain", f::PhotometricFilter{T}) where T
    N = length(f)
    min_wl = min_wave(f)
    max_wl = max_wave(f)
    piv_wl = pivot_wavelength(f)
    cen_wl = central_wavelength(f)
    eff_wl = effective_wavelength(f)
    eff_width = width(f)
    Γ = fwhm(f)
    print(io, "$N-element PhotometricFilter{$T}: ", f.name, "\n min. wave.: ", min_wl, "\n max. wave.: ", max_wl, "\n effective wave.: ", eff_wl, "\n central wave.: ", cen_wl, "\n pivot wave.: ", piv_wl, "\n eff. width: ", eff_width, "\n fwhm: ", Γ)
end

wave(f::PhotometricFilter) = f.wave
throughput(f::PhotometricFilter) = f.throughput

(f::PhotometricFilter)(wave) = f.etp(wave)

function (f::PhotometricFilter)(wave::Q) where Q <: Unitful.Length
    # strip units because Interpolations.jl doesn't work nicely with them
    wl = ustrip(unit(eltype(f.wave)), wave)
    return f.etp(wl)
end

Base.size(f::PhotometricFilter) = size(throughput(f))

"""
    effective_wavelength(::PhotometricFilter)

Return the effective wavelength using the Vega spectrum as a standard
"""
function effective_wavelength(f::PhotometricFilter)
    wvega, fvega = Vega()
    u = unit(eltype(wave(f)))
    wl = ustrip.(u, wvega)
    filt = f.(wvega) .* fvega
    norm = trapz(wl, filt)
    leff = trapz(wl, wl .* filt)
    return leff / norm * u
end

"""
    pivot_wavelength(::PhotometricFilter)

Returns the pivot wavelength of the filter, described by the equation below. Internally integration is carried out using trapezoidal integration. It can be convenient to think of this as the "center of mass" of the filter.
"""
function pivot_wavelength(f::PhotometricFilter{T,S,<:Photon}) where {T,S}
    wl = ustrip(wave(f))
    y = throughput(f) ./ wl
    norm = trapz(wl, wl .* throughput(f))
    lp2 = norm / trapz(wl, y)
    return sqrt(lp2) * unit(eltype(wave(f)))
end

function pivot_wavelength(f::PhotometricFilter{T,S,<:Energy}) where {T,S}
    wl = ustrip(wave(f))
    y = throughput(f) ./ wl.^2
    norm = trapz(wl, throughput(f))
    lp2 = norm / trapz(wl, y)
    return sqrt(lp2) * unit(eltype(wave(f)))
end

function central_wavelength(f::PhotometricFilter)
    wl = ustrip(wave(f))
    norm = trapz(wl, throughput(f))
    lt = trapz(wl, wl .* throughput(f))
    return  lt / norm * unit(eltype(wave(f)))
end

function min_wave(f::PhotometricFilter; level=0.01)
    y = throughput(f)
    idx = findfirst(q -> q / maximum(y) > level, y)
    return wave(f)[idx]
end

function max_wave(f::PhotometricFilter; level=0.01)
    y = throughput(f)
    idx = findlast(q -> q / maximum(y) > level, y)
    return wave(f)[idx]
end

"""
    apply(::PhotometricFilter, wave, flux)

Use linear interpolation to map the wavelenghts of the photometric filter to the given wavelengths and apply the filter throughput to the `flux`. The wavelengths of the filter and `wave` need to be compatible. This means if one has units, the other one needs units, too.

# See also
[`interpolator`](@ref)
"""
apply(filt::PhotometricFilter, wave, flux) = apply!(filt, wave, flux, similar(flux))

"""
    apply!(::PhotometricFilter, wave, flux, out)

In-place version of [`apply`](@ref) which modifies `out`. It should have a compatible element type with `flux`.
"""
function apply!(filt::PhotometricFilter, wave, flux, out)
    @. out = flux * filt(wave)
    return out
end

function fwhm(filt::PhotometricFilter)
    Δ = diff(sign.(filt ./ maximum(filt) .- 0.5))
    i1 = findfirst(!iszero, Δ)
    i2 = findnext(!iszero, Δ, i1 + 1)
    return wave(filt)[i2] - wave(filt)[i1]
end

function width(f::PhotometricFilter)
    norm = trapz(ustrip(wave(f)), throughput(f))
    return norm / maximum(throughput(f)) * unit(eltype(wave(f)))
end

function Base.:*(f1::PhotometricFilter, f2::PhotometricFilter)
    # find total extent and log-spacing
    w1 = wave(f1)
    w2 = wave(f2)
    e1 = extrema(w1)
    e2 = extrema(w2)
    min_wl = max(e1[1], e2[1])
    max_wl = min(e1[2], e2[2])
    max_wl ≤ min_wl && error("no overlap between the two filters")
    δwl1 = (e1[2] - e1[1]) / length(w1)
    δwl2 = (e2[2] - e2[1]) / length(w2)
    δ = min(δwl1, δwl2)

    wl = min_wl:δ:max_wl

    through = @. f1(wl) * f2(wl)
    name = "$(f1.name) * $(f2.name)"
    return PhotometricFilter(wl, through; detector=f1.detector, name=name)
end


text:
function createDensityWindow{F<:(@compat AbstractFloat)}(X::Vector{F}, Y::Vector{F};
                                               width::Int = 40,
                                               height::Int = 20,
                                               margin::Int = 3,
                                               padding::Int = 1,
                                               title::(@compat AbstractString) = "",
                                               border::Symbol = :solid,
                                               labels::Bool = true,
                                               xlim::Vector = [0.,0.],
                                               ylim::Vector = [0.,0.])
  length(xlim) == length(ylim) == 2 || throw(ArgumentError("xlim and ylim must only be vectors of length 2"))
  margin >= 0 || throw(ArgumentError("Margin must be greater than or equal to 0"))
  length(X) == length(Y) || throw(DimensionMismatch("X and Y must be the same length"))
  width = max(width, 5)
  height = max(height, 5)
  minX = float(minimum(xlim))
  maxX = float(maximum(xlim))
  if minX == 0. && maxX == 0.
    minX = minimum(X)
    maxX = maximum(X)
    minX, maxX = plottingRange(minX, maxX)
  end
  minY = float(minimum(ylim))
  maxY = float(minimum(ylim))
  if minY == 0. && maxY == 0.
    minY = minimum(Y)
    maxY = maximum(Y)
    minY, maxY = plottingRange(minY, maxY)
  end
  plotOriginX = minX; plotWidth = maxX - plotOriginX
  plotOriginY = minY; plotHeight = maxY - plotOriginY

  canvas = DensityCanvas(width, height,
                         plotOriginX = plotOriginX, plotOriginY = plotOriginY,
                         plotWidth = plotWidth, plotHeight = plotHeight)
  newPlot = Plot(canvas, title=title, margin=margin,
                 padding=padding, border=border, showLabels=labels)

  minXString=string(isinteger(minX) ? safeRound(minX) : minX)
  maxXString=string(isinteger(maxX) ? safeRound(maxX) : maxX)
  minYString=string(isinteger(minY) ? safeRound(minY) : minY)
  maxYString=string(isinteger(maxY) ? safeRound(maxY) : maxY)
  annotate!(newPlot, :l, 1, maxYString)
  annotate!(newPlot, :l, height, minYString)
  annotate!(newPlot, :bl, minXString)
  annotate!(newPlot, :br, maxXString)
  newPlot
end

function densityplot{F<:Real,R<:Real}(X::Vector{F}, Y::Vector{R};
                                      color::Symbol = :white,
                                      args...)
  X = convert(Vector{(@compat AbstractFloat)},X)
  Y = convert(Vector{(@compat AbstractFloat)},Y)
  minX = minimum(X); minY = minimum(Y)
  maxX = maximum(X); maxY = maximum(Y)
  newPlot = createDensityWindow(X, Y; args...)
  setPoint!(newPlot, X, Y, color)
end

function densityplot!{T<:Canvas,F<:Real,R<:Real}(plot::Plot{T},
                                                 X::Vector{F},
                                                 Y::Vector{R};
                                                 color::Symbol = :white,
                                                 args...)
  X = convert(Vector{(@compat AbstractFloat)},X)
  Y = convert(Vector{(@compat AbstractFloat)},Y)
  setPoint!(plot, X, Y, color)
end


text:
@testset "testing parallel solution" begin

    let royalty_rates = [0.25,],
        geology_types = 2:2,
        roy = royalty_rates[1],
        geoid = geology_types[1],
        itype = (geoid, roy,)

        rmprocs(workers())
        pids = IN_SLURM ? addprocs_slurm(parse(Int, ENV["SLURM_CPUS_PER_TASK"])) : addprocs()

        @everywhere @show pwd()
        @everywhere using ShaleDrillingModel

        sev = SharedEV(pids, prim, geology_types, royalty_rates)
        @eval @everywhere begin
            set_g_dcdp_primitives($prim)
            set_g_dcdp_tmpvars($tmpv)
            set_g_SharedEV($sev)
            set_g_ItpSharedEV(ItpSharedEV(get_g_SharedEV(), get_g_dcdp_primitives(), $σv))
        end
        s = remotecall_fetch(get_g_SharedEV, pids[2])
        fetch(s)

        zero!(sev)
        @show parallel_solve_vf_all!(sev, vcat(θt,σv), true)
        @test all(isfinite.(sev.EV))
        @test all(isfinite.(sev.dEV))
        @test all(isfinite.(sev.dEVσ))
        @test !all(sev.EV .== 0.0)
        @test !all(sev.dEV .== 0.0)
        @test !all(sev.dEVσ .== 0.0)

        let EVcopy = similar(sev.EV),
            EVcoefs = similar(sev.EV),
            sitev = get_g_ItpSharedEV()

            @test(sitev.EV.itp.coefs === sev.EV)
            # store EV
            EVcopy .= sev.EV
            #interpolate & store
            serial_prefilterByView!(sev,sitev,true)
            EVcoefs .= sev.EV
            # check interpolation happened
            @test !all(EVcopy .== EVcoefs)
            # restore coefs
            sev.EV .= EVcopy
            @test all(sitev.EV.itp.coefs .== EVcopy)
            parallel_prefilterByView!(sev,sitev,true)
            # double check that parallel does the same as serial
            @test !all(EVcopy .== sev.EV)
            @test all(EVcoefs .== sev.EV)
        end

        rmprocs(workers())
    end
end


text:
function EGraphs.preprocess(c::Chain)
    if length(c.args) > 2
        return Chain(c.args[1], EGraphs.preprocess(Chain(c.args[2:end]...)))
    end
    return c
end

TermInterface.gethead(t::Chain) = :call
TermInterface.getargs(t::Chain) = [:Chain, t.args...]
TermInterface.arity(t::Chain) = length(TermInterface.getargs(t))
TermInterface.isterm(e::Type{Chain}) = true

TermInterface.gethead(t::Gate) = :call
TermInterface.getargs(t::Gate) = [:Gate, t.operation,  t.locations]
TermInterface.arity(t::Gate) = length(TermInterface.getargs(t))
TermInterface.isterm(e::Type{Gate}) = true

TermInterface.gethead(t::Ctrl) = :call
TermInterface.getargs(t::Ctrl) = [:Ctrl, t.gate, t.ctrl]
TermInterface.arity(t::Ctrl) = length(TermInterface.getargs(t))
TermInterface.isterm(e::Type{Ctrl}) = true

function TermInterface.similarterm(x::Type{Chain}, head, args; metadata=nothing)
    @assert head == :call && args[1] == :Chain
    Chain(args[2:end]...)
end

function TermInterface.similarterm(x::Type{Gate}, head, args; metadata=nothing)
    @assert head == :call && args[1] == :Gate
    Gate(args[2], args[3])
end

function TermInterface.similarterm(x::Type{Ctrl}, head, args; metadata=nothing)
    @assert head == :call && args[1] == :Ctrl
    Ctrl(args[2], args[3])
end

postprocess(c::Chain) = c
postprocess(c::Union{Gate, Ctrl}) = Chain(c)


text:
@testset "delaydifferentialequation" begin
# Tests
# 1. x'(t) = 2x(t) - exp(1)*x(t - 1) has a double root at λ = 1
    ev = dde_stability(2, -exp(1), 1)
    @test count((abs.(real.(ev.values) .- 1) .< 1e-6) .& (abs.(imag.(ev.values)) .< 1e-6)) == 2
# 2. x'(t) = -x(t - τ) goes unstable at τ=π/2
    ev = dde_stability(0, -1, π/2 - 0.01)
    @test all(real.(ev.values) .< 0)
    ev = dde_stability(0, -1, π/2 + 0.01)
    @test any(real.(ev.values) .> 0)
# 3. x'(t) = α*A*x(t) + (1-α)*A*x(t-τ)
#    A = [
#          50  284   41   23   50   32
#        -280  -46  -19  -37  -10  -28
#          35   -1   26  143   35   17
#          5   -31 -139  -22    5  -13
#          20  -16   11   -7 -115  137
#         -10  -46  -19  -37 -145 -163
#    ]
#    with α=0.9 is stable for 3π/(270*(2*0.9 - 1)) ≤ τ ≤ 4π/270
#    that is 0.04363323129985824 ≤ τ ≤ 0.046542113386515455
    A = [
          50  284   41   23   50   32
        -280  -46  -19  -37  -10  -28
          35   -1   26  143   35   17
          5   -31 -139  -22    5  -13
          20  -16   11   -7 -115  137
         -10  -46  -19  -37 -145 -163
    ]
    α = 0.9
    ev = dde_stability(α * A, (1 - α)*A, 0.043)
    @test any(real.(ev.values) .> 0)
    ev = dde_stability(α * A, (1 - α)*A, 0.044)
    @test all(real.(ev.values) .< 0)
    ev = dde_stability(α * A, (1 - α)*A, 0.047)
    @test any(real.(ev.values) .> 0)
end


text:
##################################################################################
# This file is part of ModelBaseEcon.jl
# BSD 3-Clause License
# Copyright (c) 2020, Bank of Canada
# All rights reserved.
##################################################################################

abstract type AbstractEquation end

# equations must have these fields: expr, vinds, vsyms, eval_resid, eval_RJ
for fn in (:expr, :vinds, :vsyms, :eval_resid, :eval_RJ)
    local qnfn = QuoteNode(fn)
    quote
        @inline $fn(eqn::AbstractEquation) = getfield(eqn, $qnfn)
    end |> eval
end

const default_eqn_type = :any

# equations might have these fields. If not, we provide defaults
@inline type(eqn::AbstractEquation) = :type in fieldnames(typeof(eqn)) ? getfield(eqn, :type) : default_eqn_type
@inline doc(eqn::AbstractEquation) = :doc in fieldnames(typeof(eqn)) ? getfield(eqn, :doc) : ""

@inline eval_resid(eqn::AbstractEquation, x) = eval_resid(eqn)(x)
@inline eval_RJ(eqn::AbstractEquation, x) = eval_RJ(eqn)(x)

#
function Base.show(io::IO, eqn::AbstractEquation)
    if get(io, :compact, false)
        return print(io, expr(eqn))
    end
    docstr = isempty(doc(eqn)) ? "" : "\"$(doc(eqn))\" "
    typestr = type(eqn) === default_eqn_type ? "" : "$(type(eqn)) "
    print(io, docstr, typestr, expr(eqn))
end

Base.:(==)(e1::AbstractEquation, e2::AbstractEquation) = type(e1) == type(e1) && expr(e1) == expr(e2)
Base.hash(e::AbstractEquation, h::UInt) = hash((type(e), expr(e)), h)

abstract type AbstractModel end

# a subtype of AbstractModel is expected to have a number of fields.
# If it doesn't, the creater of the new model type must define the
# access methods that follow.

@inline variables(m::AM) where AM <: AbstractModel = getfield(m, :variables)
@inline nvariables(m::AM) where AM <: AbstractModel = length(variables(m))

@inline shocks(m::AM) where AM <: AbstractModel = getfield(m, :shocks)
@inline nshocks(m::AM) where AM <: AbstractModel = length(shocks(m))

@inline allvars(m::AM) where AM <: AbstractModel = vcat(variables(m), shocks(m))
@inline nallvars(m::AM) where AM <: AbstractModel = length(variables(m)) + length(shocks(m))

@inline sstate(m::AM) where AM <: AbstractModel = getfield(m, :sstate)

@inline parameters(m::AM) where AM <: AbstractModel = getfield(m, :parameters)

@inline equations(m::AM) where AM <: AbstractModel = getfield(m, :equations)
@inline nequations(m::AM) where AM <: AbstractModel = length(equations(m))

@inline alleqns(m::AM) where AM <: AbstractModel = getfield(m, :equations)
@inline nalleqns(m::AM) where AM <: AbstractModel = length(equations(m))

export parameters
export variables, nvariables
export shocks, nshocks
export equations, nequations
export sstate
#######

# @inline moduleof(f::Function) = parentmodule(f)
@inline moduleof(e::AbstractEquation) = parentmodule(eval_resid(e))
function moduleof(m::AbstractModel)
    eqns = equations(m)
    if isempty(eqns)
        error("Unable to determine the module containing the given model. Try adding equations to it and call `@initialize`.")
    end
    return moduleof(first(eqns))
end


text:
using jInv.Utils
using Base.Test
using OED

n = 10
N = 10
nrhs = 1

A = sdiag(1./linspace(1,N,N).^4)
x = rand(n,nrhs); x[4:7,:] = 1;
D = A*x  + 1e-8* randn(N,nrhs)/1e3

OEDparam = OEDTypeAParam(A,D,zeros(0,0))

display(OEDparam)

(Nt,nt,nrhst) = getProblemSize(OEDparam)
@test Nt==N
@test nt==n
@test nrhst ==nrhs
@test getNumberOfDesignParameters(OEDparam)==N

println("--- test getM and derivatives ---")
v = randn(n)
p0  = rand(N)
pass, = checkDerivative(p-> getM(p,OEDparam)*v, p-> getdpMx(p,v,OEDparam), p0)
@test pass

println("--- test getM^T and derivatives ---")
v = randn(N)
p0  = rand(N)
pass, = checkDerivative(p-> getM(p,OEDparam)'*v, p->getdpMTx(p,v,OEDparam), p0)
@test pass

println("--- test getD and derivatives ---")
pass, = checkDerivative(p->getD(p,OEDparam)[:,1], p->getdpD(p,1,OEDparam),p0)
@test pass


text:
using IterableTables
using NamedTuples
using JuliaDB
using DataFrames
using NullableArrays
using Base.Test

@testset "JuliaDB" begin

source_df = DataFrame(a=[1,2,3], b=[4.,5.,6.])

jdb = distribute(source_df)

@test isa(jdb, JuliaDB.DTable)

target_df = DataFrame(jdb)

@test source_df == target_df

end


text:
using SparseArrays
using LinearAlgebra
#using CUDA
#using CUDA.CUSPARSE
using Printf
using MatrixMarket

include("DiagonalSBP.jl")

const BC_DIRICHLET = 1
const BC_NEUMANN = 2
const BC_LOCKED_INTERFACE = 0
const BC_JUMP_INTERFACE   = 7

#CUDA.allowscalar(false)

⊗(A,B) = kron(A, B)

function mod_data!(δ, uf2, ge, K, vf, RS, metrics, Lw, t)

    (xf, yf) = metrics.facecoord
    coord = metrics.coord
    sJ = metrics.sJ
    μf2 = metrics.μf2

    ge .= 0

    for i in 1:4
        if i == 1
            # fault data
            vf .=  δ ./ 2
        elseif i == 2
            # face 2 data
            vf .= (uf2 .+ (RS.Vp/2 * t))
        elseif i == 3
            # face 3 data
            vf .= 0
        elseif i == 4
            # face 4 data
            vf .= 0
        end
        ge .+= K[i] * vf
    end

end

function static_data_mms!(δ, ge, K, JH, vf, MMS, B_p, RS, metrics, t)

    (xf, yf) = metrics.facecoord
    coord = metrics.coord
    sJ = metrics.sJ
    μf2 = μ(xf[2], yf[2], B_p)

    ge .= 0

    for i in 1:4
        if i == 1
            vf .=  ue(xf[1], yf[1], t, MMS)
        elseif i == 2
            # dirichlet h
            vf .= ue(xf[2], yf[2], t, MMS)
        elseif i == 3
            # neumann h
            vf .= sJ[3] .* τe(xf[3], yf[3], t, 3, B_p, MMS)
        elseif i == 4
            # neumann h
            vf .= sJ[4] .* τe(xf[4], yf[4], t, 4, B_p, MMS)
        end

        ge .+= K[i] * vf
    end

    ge .+= JH * Forcing_u(coord[1][:], coord[2][:], t, B_p, MMS)

end

function mod_data_mms!(δ, ge, K, H̃, JH, vf, MMS, B_p, RS, metrics, t)

    (xf, yf) = metrics.facecoord
    coord = metrics.coord
    sJ = metrics.sJ
    μf2 = μ(xf[2], yf[2], B_p)

    ge .= 0

    for i in 1:4
        if i == 1
            vf .=  δ ./ 2
        elseif i == 2
            # dirichlet h
            vf .= he(xf[2], yf[2], t, MMS)
        elseif i == 3
            # neumann h
            vf .= sJ[3] .* τhe(xf[3], yf[3], t, 3, B_p, MMS)
        elseif i == 4
            # neumann h
            vf .= sJ[4] .* τhe(xf[4], yf[4], t, 4, B_p, MMS)
        end

        ge .+= K[i] * vf
    end

    ge .+= JH * Forcing_h(coord[1][:], coord[2][:], t, B_p, MMS)

end

function traction(ops, metrics, f, u, û)

    HI = ops.HI[f]
    G = ops.G[f]
    Γ = ops.Γ[f]
    L = ops.L[f]
    sJ = ops.sJ[f]

    return (HI * G * u + Γ * (û - L * u)) ./ sJ

end

function operators(p, Nr, Ns, μ, ρ, R, B_p, metrics,
                     τscale = 2,
                     crr = metrics.crr,
                     css = metrics.css,
                     crs = metrics.crs)

    csr = crs
    J = metrics.J

    hr = 2/Nr
    hs = 2/Ns

    hmin = min(hr, hs)

    r = -1:hr:1
    s = -1:hs:1

    Nrp = Nr + 1
    Nsp = Ns + 1
    Np = Nrp * Nsp
    Nn = Np
    nn = Nrp
    # Derivative operators for the rest of the computation
    (Dr, HrI, Hr, r) = D1(p, Nr; xc = (-1,1))
    Qr = Hr * Dr
    QrT = sparse(transpose(Qr))

    (Ds, HsI, Hs, s) = D1(p, Ns; xc = (-1,1))
    Qs = Hs * Ds
    QsT = sparse(transpose(Qs))

    # Identity matrices for the comuptation
    Ir = sparse(I, Nrp, Nrp)
    Is = sparse(I, Nsp, Nsp)

    (_, S0e, SNe, _, _, Ae, _) = variable_D2(p, Nr, rand(Nrp))
    IArr = Array{Int64,1}(undef,Nsp * length(Ae.nzval))
    JArr = Array{Int64,1}(undef,Nsp * length(Ae.nzval))
    VArr = Array{Float64,1}(undef,Nsp * length(Ae.nzval))
    stArr = 0
    A_t = @elapsed begin
        for j = 1:Nsp
            rng = (j-1) * Nrp .+ (1:Nrp)
            (_, S0e, SNe, _, _, Ae, _) =  variable_D2(p, Nr, crr[rng])
            (Ie, Je, Ve) = findnz(Ae)
            IArr[stArr .+ (1:length(Ve))] = Ie .+ (j-1) * Nrp
            JArr[stArr .+ (1:length(Ve))] = Je .+ (j-1) * Nrp
            VArr[stArr .+ (1:length(Ve))] = Hs[j,j] * Ve
            stArr += length(Ve)
        end

        Ãrr = sparse(IArr[1:stArr], JArr[1:stArr], VArr[1:stArr], Np, Np)

        (_, S0e, SNe, _, _, Ae, _) =  variable_D2(p, Ns, rand(Nsp))
        IAss = Array{Int64,1}(undef,Nrp * length(Ae.nzval))
        JAss = Array{Int64,1}(undef,Nrp * length(Ae.nzval))
        VAss = Array{Float64,1}(undef,Nrp * length(Ae.nzval))
        stAss = 0

        for i = 1:Nrp
            rng = i .+ Nrp * (0:Ns)
            (_, S0e, SNe, _, _, Ae, _) =  variable_D2(p, Ns, css[rng])

            (Ie, Je, Ve) = findnz(Ae)
            IAss[stAss .+ (1:length(Ve))] = i .+ Nrp * (Ie .- 1)
            JAss[stAss .+ (1:length(Ve))] = i .+ Nrp * (Je .- 1)
            VAss[stAss .+ (1:length(Ve))] = Hr[i,i] * Ve
            stAss += length(Ve)

        end

        Ãss = sparse(IAss[1:stAss], JAss[1:stAss], VAss[1:stAss], Np, Np)

        Ãsr = (QsT ⊗ Ir) * sparse(1:length(crs), 1:length(crs), view(crs, :)) * (Is ⊗ Qr)
        Ãrs = (Is ⊗ QrT) * sparse(1:length(csr), 1:length(csr), view(csr, :)) * (Qs ⊗ Ir)

        Ã = Ãrr + Ãss + Ãrs + Ãsr
    end

    #@printf "Got Ã in %f seconds\n" A_t

    # volume quadrature
    H̃ = kron(Hr, Hs)
    H̃inv = spdiagm(0 => 1 ./ diag(H̃))

    # diagonal rho
    rho = ρ(metrics.coord[1], metrics.coord[2], B_p)
    rho = reshape(rho, Nrp*Nsp)
    P̃ = spdiagm(0 => rho)
    P̃inv = spdiagm(0 => (1 ./ rho))

    JI = spdiagm(0 => reshape(metrics.JI, Nrp*Nsp))

    er0 = sparse([1  ], [1], [1], Nrp, 1)
    erN = sparse([Nrp], [1], [1], Nrp, 1)
    es0 = sparse([1  ], [1], [1], Nsp, 1)
    esN = sparse([Nsp], [1], [1], Nsp, 1)

    er0T = sparse([1], [1  ], [1], 1, Nrp)
    erNT = sparse([1], [Nrp], [1], 1, Nrp)
    es0T = sparse([1], [1  ], [1], 1, Nsp)
    esNT = sparse([1], [Nsp], [1], 1, Nsp)

    cmax = maximum([maximum(crr), maximum(crs), maximum(css)])

    # Surface quadtrature matrices
    H1 = H2 = Hs
    H3 = H4 = Hr

    H = (Hs, Hs, Hr, Hr)
    HI = (HsI, HsI, HrI, HrI)
    # Volume to Face Operators (transpose of these is face to volume)
    L = (convert(SparseMatrixCSC{Float64, Int64}, kron(Ir, es0)'),
         convert(SparseMatrixCSC{Float64, Int64}, kron(Ir, esN)'),
         convert(SparseMatrixCSC{Float64, Int64}, kron(er0, Is)'),
         convert(SparseMatrixCSC{Float64, Int64}, kron(erN, Is)'))
    # coefficent matrices
    Crr1 = spdiagm(0 => crr[1, :])
    Crs1 = spdiagm(0 => crs[1, :])
    Csr1 = spdiagm(0 => crs[1, :])
    Css1 = spdiagm(0 => css[1, :])

    Crr2 = spdiagm(0 => crr[Nrp, :])
    Crs2 = spdiagm(0 => crs[Nrp, :])
    Csr2 = spdiagm(0 => crs[Nrp, :])
    Css2 = spdiagm(0 => css[Nrp, :])

    Css3 = spdiagm(0 => css[:, 1])
    Crs3 = spdiagm(0 => crs[:, 1])
    Csr3 = spdiagm(0 => crs[:, 1])
    Crr3 = spdiagm(0 => crr[:, 1])

    Css4 = spdiagm(0 => css[:, Nsp])
    Crs4 = spdiagm(0 => crs[:, Nsp])
    Csr4 = spdiagm(0 => crs[:, Nsp])
    Crr4 = spdiagm(0 => crr[:, Nsp])

    (_, S0, SN, _, _) = D2(p, Nr, xc=(-1,1))[1:5]
    S0 = sparse(Array(S0[1,:])')
    SN = sparse(Array(SN[end, :])')

    # Boundars Derivatives
    B1r =  Crr1 * kron(Is, S0)
    B1s = Crs1 * L[1] * kron(Ds, Ir)
    B2r = Crr2 * kron(Is, SN)
    B2s = Crs2 * L[2] * kron(Ds, Ir)
    B3s = Css3 * kron(S0, Ir)
    B3r = Csr3 * L[3] * kron(Is, Dr)
    B4s = Css4 * kron(SN, Ir)
    B4r = Csr4 * L[4] * kron(Is, Dr)

    (xf1, xf2, xf3, xf4) = metrics.facecoord[1]
    (yf1, yf2, yf3, yf4) = metrics.facecoord[2]

    Z̃f = (metrics.sJ[1] .* sqrt.(ρ(xf1, yf1, B_p) .* μ(xf1, yf1, B_p)),
          metrics.sJ[2] .* sqrt.(ρ(xf2, yf2, B_p) .* μ(xf2, yf2, B_p)),
          metrics.sJ[3] .* sqrt.(ρ(xf3, yf3, B_p) .* μ(xf3, yf3, B_p)),
          metrics.sJ[4] .* sqrt.(ρ(xf4, yf4, B_p) .* μ(xf4, yf4, B_p)))

    # Penalty terms

    if p == 2
        l = 2
        α = 1.0
        θ_R = 1 / 2
    elseif p == 4
        l = 4
        β = 0.2505765857
        α = 0.5776
        θ_R = 17 / 48
    elseif p == 6
        l = 7
        β = 0.1878687080
        α = 0.3697
        θ_R = 13649 / 43200
    else
        error("unknown order")
    end

    ψmin_r = reshape(crr, Nrp, Nsp)
    ψmin_s = reshape(css, Nrp, Nsp)
    @assert minimum(ψmin_r) > 0
    @assert minimum(ψmin_s) > 0

    hr = 2 / Nr
    hs = 2 / Ns

    ψ1 = ψmin_r[  1, :]
    ψ2 = ψmin_r[Nrp, :]
    ψ3 = ψmin_s[:,   1]
    ψ4 = ψmin_s[:, Nsp]

    for k = 2:l
        ψ1 = min.(ψ1, ψmin_r[k, :])
        ψ2 = min.(ψ2, ψmin_r[Nrp+1-k, :])
        ψ3 = min.(ψ3, ψmin_s[:, k])
        ψ4 = min.(ψ4, ψmin_s[:, Nsp+1-k])
    end

    τR1 = (1/(α*hr))*Is
    τR2 = (1/(α*hr))*Is
    τR3 = (1/(α*hs))*Ir
    τR4 = (1/(α*hs))*Ir

    p1 = ((crr[  1, :]) ./ ψ1)
    p2 = ((crr[Nrp, :]) ./ ψ2)
    p3 = ((css[:,   1]) ./ ψ3)
    p4 = ((css[:, Nsp]) ./ ψ4)

    P1 = sparse(1:Nsp, 1:Nsp, p1)
    P2 = sparse(1:Nsp, 1:Nsp, p2)
    P3 = sparse(1:Nrp, 1:Nrp, p3)
    P4 = sparse(1:Nrp, 1:Nrp, p4)

    # dynamic penalty matrices
    Γ = ((2/(θ_R*hr))*Is + τR1 * P1,
         (2/(θ_R*hr))*Is + τR2 * P2,
         (2/(θ_R*hs))*Ir + τR3 * P3,
         (2/(θ_R*hs))*Ir + τR4 * P4)

    JH = sparse(1:Np, 1:Np, view(J, :)) * (Hs ⊗ Hr)

    JIHP = JI * H̃inv * P̃inv

    Cf = ((Crr1, Crs1), (Crr2, Crs2), (Css3, Csr3), (Css4, Csr4))
    B = ((B1r, B1s), (B2r, B2s), (B3s, B3r), (B4s, B4r))
    nl = (-1, 1, -1, 1)
    G = (-H[1] * (B[1][1] + B[1][2]),
         H[2] * (B[2][1] + B[2][2]),
         -H[3] * (B[3][1] + B[3][2]),
         H[4] * (B[4][1] + B[4][2]))

    static_t = @elapsed begin

        # boundary data operators for quasi-static displacement conditions
        K1 = L[1]' * H[1] * Cf[1][1] * Γ[1] - G[1]'
        K2 = L[2]' * H[2] * Cf[2][1] * Γ[2] - G[2]'
        #K3 = L[3]' * H[3] * Γ[3] - G[3]'
        #K4 = L[4]' * H[4] * Γ[4] - G[4]'

        # boundary data operator for quasi-static traction-free conditions
        #K1 = L[1]' * H2]
        #K2 = L[2]' * H[2]
        K3 = L[3]' * H[3]
        K4 = L[4]' * H[4]

        # modification of second derivative operator for displacement conditions
        M̃ = copy(Ã)
        for f in 1:2
            M̃ -= L[f]' * G[f]
            M̃ += L[f]' * H[f] * Cf[f][1] * Γ[f] * L[f]
            M̃ -= G[f]' * L[f]
        end

    end

    Λ_t = @elapsed begin
        faces = [1,2,4]
        dv_u = -Ã

        for i in faces
            dv_u .+= (L[i]' * H[i] * ((1 - R[i])/2 .* (nl[i] * (B[i][1] + B[i][2]) - Cf[i][1] * Γ[i] * L[i]))) +
                nl[i] * (B[i][1]' + B[i][2]') * H[i] * L[i]

        end

        dv_v = spzeros(Nn, Nn)
        for i in faces
            dv_v .+=  L[i]' * H[i] * (-(1 - R[i])/2 .* Z̃f[i] .* L[i])
        end

        dv_û = spzeros(Nn, 4nn)
        dû_u = spzeros(4nn, Nn)
        dû_v = spzeros(4nn, Nn)

        for i in faces
            dv_û[ : , (i-1) * nn + 1 : i * nn] .=
                (L[i]' * H[i] * ((1 - R[i])/2 .* Cf[i][1] * Γ[i])) -
                nl[i] * (B[i][1]' + B[i][2]') * H[i]

            dû_u[(i-1) * nn + 1 : i * nn , : ] .=
                -(1 + R[i])/2 .* (nl[i] * (B[i][1] + B[i][2]) -
                Cf[i][1] * Γ[i] * L[i])./Z̃f[i]

            dû_v[(i-1) * nn + 1 : i * nn , : ] .= (1 + R[i])/2 .* L[i]
        end

        dû_û = spzeros(4nn, 4nn)
        for i in faces
            dû_û[(i-1) * nn + 1 : i * nn, (i-1) * nn + 1 : i * nn] .= -(1 + R[i])/2 .* (Cf[i][1] * Γ[i])./Z̃f[i]
        end

        dû_ψ = spzeros(4nn, nn)

        Λ = [ spzeros(Nn, Nn) sparse(I, Nn, Nn) spzeros(Nn, 5nn)
              dv_u dv_v dv_û spzeros(Nn, nn)
              dû_u dû_v dû_û dû_ψ
              spzeros(nn, 2Nn + 5nn) ]

        nCnΓ1 = Crr1 * Γ[1]
        HIGΓL1 = nl[1] * (B[1][1] + B[1][2]) - nCnΓ1 * L[1]
    end

    #@printf "Got Λ and friends in %f seconds\n" Λ_t

    (Λ = Λ,
     M̃ = cholesky(Symmetric(M̃)),
     K = (K1, K2, K3, K4),
     G = G,
     Crr = Crr1,
     Γ = Γ,
     HI = HI,
     P̃I = P̃inv,
     H̃ = H̃,
     H̃I = H̃inv,
     JI = JI,
     JIHP = JIHP,
     nCnΓ1 = nCnΓ1,
     HIGΓL1 = HIGΓL1,
     hmin = hmin,
     cmax = cmax,
     JH = JH,
     sJ = metrics.sJ,
     nx = metrics.nx,
     ny = metrics.ny,
     L = L,
     H = H,
     Z̃f = Z̃f)

end


text:
# Test tail sampling of TruncatedNormal

using Distributions
using Base.Test

# tests tail sampling, #343, #493
@test isfinite(rand(TruncatedNormal(0,1,35,Inf)))
@test isfinite(rand(TruncatedNormal(0,1,-Inf,-35)))


text:
using Test
using Jecco
using Jecco.AdS5_3_1

sigmavars   = tuple(ones(58)...)

axx, ayy, axy, bx, by, cc= AdS5_3_1.AH_eq_coeff(sigmavars, AdS5_3_1.Outer())

S = AdS5_3_1.AH_eq_res(sigmavars, AdS5_3_1.Outer())

print("(axx, ayy, axy, bx, by, cc, S) =",(axx, ayy, axy, bx, by, cc, S),"\n")


text:
mutable struct EventHandler
    handlers::Dict{String,Function}
end
EventHandler() = EventHandler(Dict{String,Function}())
function set!(eh::EventHandler, event_name::String, handler::Function=DefaultHandler())
    eh.handlers[event_name] = handler
end

DefaultHandler(args...) = nothing
get_handler(eh::EventHandler, event_name::String) = get(eh.handlers, event_name, DefaultHandler)


text:
# TODO: specialize to different call sites.
mutable struct UpdateContext{C <: CallSite,
                             T <: Trace,
                             K <: ConstrainedSelection,
                             P <: Parameters,
                             D <: Diff} <: ExecutionContext
    prev::C
    tr::T
    select::K
    weight::Float64
    score::Float64
    discard::HierarchicalTrace
    visited::Visitor
    params::P
    argdiffs::D

    # Re-write with dispatch for specialized vs. black box.
    UpdateContext(cl::C, select::K, argdiffs::D) where {C <: CallSite, K <: ConstrainedSelection, D <: Diff} = new{C, typeof(cl.trace), K, EmptyParameters, D}(cl, typeof(cl.trace)(), select, 0.0, 0.0, Trace(), Visitor(), Parameters(), argdiffs)
end
Update(cl, select, argdiffs) = UpdateContext(cl, select, argdiffs)

# Update has a special dynamo.
@dynamo function (mx::UpdateContext)(a...)
    ir = IR(a...)
    ir == nothing && return
    recur!(ir)
    return ir
end

# ------------ Convenience ------------ #

function update(ctx::UpdateContext, bbcs::HierarchicalCallSite, args...) where D <: Diff
    ret = ctx(bbcs.fn, args...)
    return ret, HierarchicalCallSite(ctx.tr, ctx.score, bbcs.fn, args, ret), ctx.weight, UndefinedChange(), ctx.discard
end

function update(sel::L, bbcs::HierarchicalCallSite) where L <: ConstrainedSelection
    argdiffs = NoChange()
    ctx = UpdateContext(bbcs, sel, argdiffs)
    return update(ctx, bbcs, bbcs.args...)
end

function update(sel::L, bbcs::HierarchicalCallSite, argdiffs::D, new_args...) where {L <: ConstrainedSelection, D <: Diff}
    ctx = UpdateContext(bbcs, sel, argdiffs)
    return update(ctx, bbcs, new_args...)
end

# TODO: disallowed for now.
#function update(sel::L, vcs::VectorizedCallSite{typeof(plate)}, argdiffs::D, new_args...) where {L <: ConstrainedSelection, D <: Diff}
#    addr = gensym()
#    v_sel = selection(addr => sel)
#    ctx = UpdateContext(vcs, v_sel, argdiffs)
#    ret = ctx(plate, addr, vcs.fn, new_args...)
#    return ret, VectorizedCallSite{typeof(plate)}(ctx.tr, ctx.score, vcs.fn, vcs.args, ret), ctx.weight, UndefinedChange(), ctx.discard
#end

function update(sel::L, vcs::VectorizedCallSite{typeof(plate)}) where {L <: ConstrainedSelection, D <: Diff}
    argdiffs = NoChange()
    addr = gensym()
    v_sel = selection(addr => sel)
    ctx = UpdateContext(vcs, v_sel, argdiffs)
    ret = ctx(plate, addr, vcs.fn, vcs.args)
    return ret, VectorizedCallSite{typeof(plate)}(ctx.tr, ctx.score, vcs.fn, vcs.args, ret), ctx.weight, UndefinedChange(), ctx.discard
end

function update(sel::L, vcs::VectorizedCallSite{typeof(markov)}) where {L <: ConstrainedSelection, D <: Diff}
    argdiffs = NoChange()
    addr = gensym()
    v_sel = selection(addr => sel)
    ctx = UpdateContext(vcs, v_sel, argdiffs)
    ret = ctx(markov, addr, vcs.fn, vcs.args[1], vcs.args[2]...)
    return ret, VectorizedCallSite{typeof(markov)}(ctx.tr, ctx.score, vcs.fn, vcs.args, ret), ctx.weight, UndefinedChange(), ctx.discard
end

function update(sel::L, vcs::VectorizedCallSite{typeof(markov)}, d::NoChange, len::Int) where {L <: ConstrainedSelection, D <: Diff}
    addr = gensym()
    v_sel = selection(addr => sel)
    ctx = UpdateContext(vcs, v_sel, d)
    ret = ctx(markov, addr, vcs.fn, len, vcs.args[2]...)
    return ret, VectorizedCallSite{typeof(markov)}(ctx.tr, ctx.score, vcs.fn, vcs.args, ret), ctx.weight, UndefinedChange(), ctx.discard
end

function update(sel::L, vcs::VectorizedCallSite{typeof(markov)}, len::Int) where {L <: ConstrainedSelection, D <: Diff}
    return update(sel, vcs, NoChange(), len)
end

# ------------ includes ------------ #

include("hierarchical/update.jl")
include("plate/update.jl")
include("markov/update.jl")
include("factor/update.jl")

# ------------ Documentation ------------ #

@doc(
"""
```julia
mutable struct UpdateContext{C <: CallSite,
                             T <: Trace,
                             K <: ConstrainedSelection,
                             P <: Parameters,
                             D <: Diff} <: ExecutionContext
    prev::C
    tr::T
    select::K
    weight::Float64
    score::Float64
    discard::HierarchicalTrace
    visited::Visitor
    params::P
    argdiffs::D
end
```

Inner constructor:

```julia
UpdateContext(cl::C, select::K, argdiffs::D) where {C <: CallSite, K <: ConstrainedSelection, D <: Diff} = new{C, typeof(cl.trace), K, EmptyParameters, D}(cl, typeof(cl.trace)(), select, 0.0, 0.0, Trace(), Visitor(), Parameters(), argdiffs)
```

`UpdateContext` is an execution context used for updating the value of random choices in an existing recorded call site. This context will perform corrective updates to the likehood weights and scores so that this operation produces the correct weights and scores for the original model program constrained with the `select` selection in the `UpdateContext`.
""", UpdateContext)

@doc(
"""
```julia
ret, cl, w, retdiff, d = update(ctx::UpdateContext, bbcs::HierarchicalCallSite, args...) where D <: Diff
ret, cl, w, retdiff, d = update(sel::L, bbcs::HierarchicalCallSite) where L <: ConstrainedSelection
ret, cl, w, retdiff, d = update(sel::L, bbcs::HierarchicalCallSite, argdiffs::D, new_args...) where {L <: ConstrainedSelection, D <: Diff}
ret, v_cl, w, retdiff, d = update(sel::L, vcs::VectorizedCallSite{typeof(plate)}) where {L <: ConstrainedSelection, D <: Diff}
ret, v_cl, w, retdiff, d = update(sel::L, vcs::VectorizedCallSite{typeof(markov)}) where {L <: ConstrainedSelection, D <: Diff}
ret, v_cl, w, retdiff, d = update(sel::L, vcs::VectorizedCallSite{typeof(markov)}, d::NoChange, len::Int) where {L <: ConstrainedSelection, D <: Diff}
ret, v_cl, w, retdiff, d = update(sel::L, vcs::VectorizedCallSite{typeof(markov)}, len::Int) where {L <: ConstrainedSelection, D <: Diff}
```

`update` provides an API to the `UpdateContext` execution context. You can use this function on any of the matching signatures above - it will return the return value `ret`, the updated `RecordSite` instance `cl` or `v_cl`, the updated weight `w`, a `Diff` instance for the return value `retdiff`, and a structure which contains any changed (i.e. discarded) record sites `d`.
""", update)


text:
using Parameters, JuMP
using JuMP.Containers: DenseAxisArray

# --- Utility ---

data(a::DenseAxisArray) = a.data
data(a::Any) = a

round_int(x::AbstractFloat) = Integer(round(x))
convert_int(x, ::Type{<:Integer}) = round_int(x)
convert_int(x, ::Type{Array{T, N}}) where T <: Integer where N = round_int.(x)
convert_int(x, ::Any) = x

"""The function queries values from the model to data type based on its field names. It extracts values from DenseAxisArray from its `data` field. Then, it converts the values to the corresponding field type. The function rounds integers before conversion because JuMP outputs integer variables as floats."""
function model_to_dtype(dtype::DataType, model::Model)
    fields = (
        value.(model[n]) |> data |> x -> convert_int(x, t)
        for (n, t) in zip(fieldnames(dtype), fieldtypes(dtype)))
    dtype(fields...)
end

# --- Model ---

"""Defines the ExampleModel type."""
const ExampleModel = Model

"""Specification for different model scenarios. For example, we can specify toggling on and off certain constraints and objectives.

# Arguments
- `relax_integer::Bool=false`: If true, relax integer constraints from variables.
"""
@with_kw struct Specs
    relax_integer::Bool=false
end

"""Contains indices of the model."""
@with_kw struct Indices
    m::Integer
    n::Integer
    k::Integer
end

"""Contains the parameters of the model."""
@with_kw struct Params
    a::Array{AbstractFloat, 1}
    b::Array{AbstractFloat, 1}
    c::Array{AbstractFloat, 1}
    A::Array{AbstractFloat, 2}
    B::Array{AbstractFloat, 2}
end

"""Contains the variable values of the model."""
@with_kw struct Variables
    x::Array{AbstractFloat, 1}
    y::Array{Integer, 1}
end

"""Contains the objective values of the model."""
@with_kw struct Objectives
    f::AbstractFloat
end

"""Queries the variable values from the model into Variables type."""
Variables(model::ExampleModel) = model_to_dtype(Variables, model)

"""Queries the objective values from the model into Objectives type."""
Objectives(model::ExampleModel) = model_to_dtype(Objectives, model)

"""Initializes the ExampleModel."""
function ExampleModel(specs::Specs, indices::Indices, params::Params)
    @unpack m, n, k = indices
    @unpack a, b, c, A, B = params

    # Initialize the model
    model = ExampleModel()

    # Variables
    @variable(model, x[1:n]≥0)
    @variable(model, y[1:k]≥0, Int)
    if specs.relax_integer
        unset_integer.(y)
    end

    # Objectives
    @expression(model, f,
        sum(a[i]*x[i] for i in 1:n) +
        sum(b[i]*y[i] for i in 1:k)
    )
    @objective(model, Max, f)

    # Constraints
    @constraint(model, [i = 1:m],
        sum(A[i, j]*x[j] for j in 1:n) +
        sum(B[i, j]*y[j] for j in 1:k) ≤ c[i])

    return model
end

max_stars_repo_name: nnhjy/JuMP-project-template

max_issues_repo_name: nnhjy/JuMP-project-template

max_forks_repo_name: nnhjy/JuMP-project-template