text:
"""
$(TYPEDEF)

Defines constant system bath coupling operators.

# Fields

$(FIELDS)
"""
struct ConstantCouplings <: AbstractCouplings
    "1-D array for independent coupling operators"
    mats::Vector{AbstractMatrix}
    "String representation for the coupling (for display purpose)"
    str_rep::Union{Vector{String},Nothing}
end

(c::ConstantCouplings)(t) = c.mats
Base.iterate(c::ConstantCouplings, state = 1) =
    state > length(c.mats) ? nothing : ((x) -> c.mats[state], state + 1)
Base.getindex(c::ConstantCouplings, inds...) = (x) -> getindex(c.mats, inds...)
Base.length(c::ConstantCouplings) = length(c.mats)
Base.eltype(c::ConstantCouplings) = typeof(c.mats[1])
Base.size(c::ConstantCouplings) = size(c.mats[1])
Base.size(c::ConstantCouplings, d) = size(c.mats[1], d)

"""
    function ConstantCouplings(mats::Union{Vector{Matrix{T}},Vector{SparseMatrixCSC{T,Int}}}; unit=:h) where {T<:Number}

Constructor of `ConstantCouplings` object. `mats` is 1-D array of matrices. `str_rep` is the optional string representation of the coupling terms. `unit` is the unit one -- `:h` or `:ħ`. The `mats` will be scaled by ``2π`` is unit is `:h`.
"""
function ConstantCouplings(
    mats::Union{Vector{Matrix{T}},Vector{SparseMatrixCSC{T,Int}}};
    unit = :h,
) where {T<:Number}
    msize = size(mats[1])
    if msize[1] <= 10
        if issparse(mats[1])
            @warn "For matrices smaller than 10×10, use StaticArrays by default."
            mats = Array.(mats)
        end
        mats = [SMatrix{msize[1],msize[2]}(unit_scale(unit) * m) for m in mats]
    else
        mats = unit_scale(unit) .* mats
    end
    ConstantCouplings(mats, nothing)
end

"""
    function ConstantCouplings(c::Vector{T}; sp = false, unit=:h) where T <: AbstractString

If the first argument is a 1-D array of strings. The constructor will automatically construct the matrics represented by the string representations.
"""
function ConstantCouplings(
    c::Vector{T};
    sp = false,
    unit = :h,
) where {T<:AbstractString}
    mats = q_translate.(c, sp = sp)
    msize = size(mats[1])
    if msize[1] <= 10
        if sp
            @warn "For matrices smaller than 10×10, use StaticArrays by default."
            mats = Array.(mats)
        end
        mats = [SMatrix{msize[1],msize[2]}(unit_scale(unit) * m) for m in mats]
    else
        mats = unit_scale(unit) .* mats
    end
    ConstantCouplings(mats, c)
end

"""
    function collective_coupling(op, num_qubit; sp = false, unit = :h)

Create `ConstantCouplings` object with operator `op` on each qubits. `op` is the string representation of one of the Pauli matrices. `num_qubit` is the total number of qubits. `sp` set whether to use sparse matrices. `unit` set the unit one -- `:h` or `:ħ`.
"""
function collective_coupling(op, num_qubit; sp = false, unit = :h)
    res = Vector{String}()
    for i = 1:num_qubit
        temp = "I"^(i - 1) * uppercase(op) * "I"^(num_qubit - i)
        push!(res, temp)
    end
    ConstantCouplings(res; sp = sp, unit = unit)
end

Base.summary(C::ConstantCouplings) = string(
    TYPE_COLOR,
    nameof(typeof(C)),
    NO_COLOR,
    " with ",
    TYPE_COLOR,
    eltype(C.mats[1]),
    NO_COLOR,
)

function Base.show(io::IO, C::ConstantCouplings)
    println(io, summary(C))
    print(io, "and string representation: ")
    show(io, C.str_rep)
end

"""
$(TYPEDEF)

Defines a single time dependent system bath coupling operator. It is defined as ``S(s)=∑f(s)×M``.  Keyword argument `unit` set the unit one -- `:h` or `:ħ`.

# Fields

$(FIELDS)

# Examples
```julia-repl
julia> TimeDependentCoupling([(s)->s], [σz], unit=:ħ)
```
"""
struct TimeDependentCoupling
    "1-D array of time dependent functions"
    funcs::Any
    "1-D array of constant matrics"
    mats::Any

    function TimeDependentCoupling(funcs, mats; unit = :h)
        new(funcs, unit_scale(unit) * mats)
    end
end

(c::TimeDependentCoupling)(t) = sum((x) -> x[1](t) * x[2], zip(c.funcs, c.mats))
Base.size(c::TimeDependentCoupling) = size(c.mats[1])
Base.size(c::TimeDependentCoupling, d) = size(c.mats[1], d)

abstract type AbstractTimeDependentCouplings <: AbstractCouplings end

"""
$(TYPEDEF)

Defines an 1-D array of time dependent system bath coupling operators.

# Fields

$(FIELDS)
"""
struct TimeDependentCouplings <: AbstractTimeDependentCouplings
    "A tuple of single `TimeDependentCoupling` operators"
    coupling::Tuple

    function TimeDependentCouplings(args...)
        new(args)
    end
end

(c::TimeDependentCouplings)(t) = [x(t) for x in c.coupling]
Base.size(c::TimeDependentCouplings) = size(c.coupling[1])
Base.size(c::TimeDependentCouplings, d) = size(c.coupling[1], d)

"""
$(TYPEDEF)

`CustomCouplings` is a container for any user defined coupling operators.

# Fields

$(FIELDS)
"""
struct CustomCouplings <: AbstractTimeDependentCouplings
    "A 1-D array of callable objects that returns coupling matrices"
    coupling::Any
    "Size of the coupling operator"
    size::Any
end

"""
$(SIGNATURES)

Create a `CustomCouplings` object from a list of functions `funcs`.
"""
function CustomCouplings(funcs; unit = :h)
    mat = funcs[1](0.0)
    if unit == :h
        funcs = [(s) -> 2π * f(s) for f in funcs]
    elseif unit != :ħ
        throw(ArgumentError("The unit can only be :h or :ħ."))
    end
    CustomCouplings(funcs, size(mat))
end

(c::CustomCouplings)(s) = [x(s) for x in c.coupling]
Base.size(c::CustomCouplings) = c.size
Base.size(c::CustomCouplings, d) = c.size[d]
Base.getindex(c::AbstractTimeDependentCouplings, inds...) = getindex(c.coupling, inds...)
Base.iterate(c::AbstractTimeDependentCouplings, state = 1) =
    Base.iterate(c.coupling, state)
Base.length(c::AbstractTimeDependentCouplings) = length(c.coupling)
Base.eltype(c::AbstractTimeDependentCouplings) = typeof(c.coupling[1])


text:
type Const <: Op; init; initialized; out; Const(;init=nothing,o...)=new(init,false); end
_KENV[:const]=Const
_KENV[:arr]=Const               # old name, to be deprecated
infersize(a::Const,ysize)=tuple(size(a.init))
ninputs(::Const)=0
canoverwrite(::Const)=false
back_reads_x(::Const)=false
back_reads_y(::Const)=false
back(::Const,dy;o...)=nothing

function forw(a::Const,y;o...)
    if !a.initialized
        a.out=copysync!(y, a.init)
        a.initialized=true
    end
    a.out===y || error("Constant modified")
    return y
end

ext: jl


text:
using Core: LineInfoNode
const NullLineInfo = LineInfoNode(@__MODULE__, Symbol(""), Symbol(""), 0, 0)

if false
    import Base: Base, @show
else
    macro show(s)
        return :(println(stdout, $(QuoteNode(s)), " = ", $(esc(s))))
    end
end

include("compiler/ssair/ir.jl")
include("compiler/ssair/domtree.jl")
include("compiler/ssair/slot2ssa.jl")
include("compiler/ssair/queries.jl")
include("compiler/ssair/passes.jl")
include("compiler/ssair/inlining2.jl")
include("compiler/ssair/verify.jl")
include("compiler/ssair/legacy.jl")
#@isdefined(Base) && include("compiler/ssair/show.jl")

function normalize_expr(stmt::Expr)
    if stmt.head === :gotoifnot
        return GotoIfNot(stmt.args[1], stmt.args[2]::Int)
    elseif stmt.head === :return
        return (length(stmt.args) == 0) ? ReturnNode(nothing) : ReturnNode(stmt.args[1])
    elseif stmt.head === :unreachable
        return ReturnNode()
    else
        return stmt
    end
end

function normalize(@nospecialize(stmt), meta::Vector{Any}, table::Vector{LineInfoNode}, loc::RefValue{Int})
    if isa(stmt, Expr)
        if stmt.head == :meta
            args = stmt.args
            if length(args) > 0
                a1 = args[1]
                if a1 === :push_loc
                    let
                        current = loc[]
                        filename = args[2]::Symbol
                        methodname = NullLineInfo.method
                        mod = table[current].mod
                        line = 0
                        for i = 3:length(args)
                            ai = args[i]
                            if ai isa Symbol
                                methodname = ai
                            elseif ai isa Int32
                                line = Int(ai)
                            elseif ai isa Int64
                                line = Int(ai)
                            elseif ai isa Module
                                mod = ai
                            end
                        end
                        push!(table, LineInfoNode(mod, methodname, filename, line, current))
                        loc[] = length(table)
                    end
                elseif a1 === :pop_loc
                    n = (length(args) > 1) ? args[2]::Int : 1
                    for i in 1:n
                        current = loc[]
                        current = table[current].inlined_at
                        current == 0 && break
                        loc[] = current
                    end
                else
                    push!(meta, stmt)
                end
            end
            return nothing
        elseif stmt.head === :line
            return nothing # deprecated - we shouldn't encounter this
        else
            return normalize_expr(stmt)
        end
    elseif isa(stmt, LabelNode)
        return nothing
    elseif isa(stmt, LineNumberNode)
        let # need to expand this node so that it is source-location independent
            current = loc[]
            info = table[current]
            methodname = info.method
            mod = info.mod
            file = stmt.file
            file isa Symbol || (file = info.file)
            line = stmt.line
            push!(table, LineInfoNode(mod, methodname, file, line, info.inlined_at))
            loc[] = length(table)
        end
        return nothing
    end
    return stmt
end

function just_construct_ssa(ci::CodeInfo, code::Vector{Any}, nargs::Int, linetable::Vector{LineInfoNode})
    mod = linetable[1].mod
    # Go through and add an unreachable node after every
    # Union{} call. Then reindex labels.
    idx = 1
    while idx <= length(code)
        stmt = code[idx]
        if isexpr(stmt, :(=))
            stmt = stmt.args[2]
        end
        if isa(stmt, Expr) && stmt.typ === Union{}
            if !(idx < length(code) && isexpr(code[idx+1], :unreachable))
                insert!(code, idx + 1, ReturnNode())
                idx += 1
            end
        end
        idx += 1
    end
    reindex_labels!(code) # update labels changed above

    inbounds_depth = 0 # Number of stacked inbounds
    meta = Any[]
    lines = fill(0, length(code))
    flags = fill(0x00, length(code))
    let loc = RefValue(1)
        for i = 1:length(code)
            stmt = code[i]
            if isexpr(stmt, :inbounds)
                arg1 = stmt.args[1]
                if arg1 === true # push
                    inbounds_depth += 1
                elseif arg1 === false # clear
                    inbounds_depth = 0
                elseif inbounds_depth > 0 # pop
                    inbounds_depth -= 1
                end
                stmt = nothing
            else
                stmt = normalize(stmt, meta, linetable, loc)
            end
            code[i] = stmt
            if !(stmt === nothing)
                lines[i] = loc[]
                if inbounds_depth > 0
                     flags[i] |= IR_FLAG_INBOUNDS
                end
            end
        end
    end
    code = strip_trailing_junk!(code, lines, flags)
    cfg = compute_basic_blocks(code)
    defuse_insts = scan_slot_def_use(nargs, ci, code)
    @timeit "domtree 1" domtree = construct_domtree(cfg)
    ir = let code = Any[nothing for _ = 1:length(code)]
             argtypes = ci.slottypes[1:(nargs+1)]
            IRCode(code, Any[], lines, flags, cfg, linetable, argtypes, mod, meta)
        end
    @timeit "construct_ssa" ir = construct_ssa!(ci, code, ir, domtree, defuse_insts, nargs)
    return ir
end

function run_passes(ci::CodeInfo, nargs::Int, linetable::Vector{LineInfoNode}, sv::OptimizationState)
    ir = just_construct_ssa(ci, copy(ci.code), nargs, linetable)
    #@Base.show ("after_construct", ir)
    # TODO: Domsorting can produce an updated domtree - no need to recompute here
    @timeit "compact 1" ir = compact!(ir)
    #@timeit "verify 1" verify_ir(ir)
    @timeit "Inlining" ir = ssa_inlining_pass!(ir, linetable, sv)
    #@timeit "verify 2" verify_ir(ir)
    @timeit "domtree 2" domtree = construct_domtree(ir.cfg)
    ir = compact!(ir)
    #@Base.show ("before_sroa", ir)
    @timeit "SROA" ir = getfield_elim_pass!(ir, domtree)
    #@Base.show ir.new_nodes
    #@Base.show ("after_sroa", ir)
    ir = adce_pass!(ir)
    #@Base.show ("after_adce", ir)
    @timeit "type lift" ir = type_lift_pass!(ir)
    @timeit "compact 3" ir = compact!(ir)
    #@Base.show ir
    @timeit "verify 3" (verify_ir(ir); verify_linetable(linetable))
    return ir
end

max_stars_repo_name: ninjin/julia

max_issues_repo_name: ninjin/julia

max_forks_repo_name: ninjin/julia


text:
using JuMP, Cbc, DelimitedFiles

function distances_from_coords(node_coords::Array{Float32,2})
    """
    Recieves a 2D array of the X,Y coordinates for each of the n nodes.
    E.g.
    node_coords[node,1] : Latitude of the node

    Returns a n,n matrix of the euclidean distances between the nodes
    """
    node_number = size(node_coords,1)
    node_distances = zeros(node_number,node_number)
    for row in range(1,length=node_number)
        for column in range(1,length=node_number)
            #Applying Euclidean distance
            node_distances[row,column] = sqrt(
                Float32(node_coords[column,1]-node_coords[row,1])^2 +
                Float32(node_coords[column,2]-node_coords[row,2])^2 )
        end
    end

    return node_distances
end

function tsp_main(args)
    """
    Use Cbc to calculate Travelling Salesperson solution for a given graph.
    Args: Filename of the parsed graph file.
        The file should follow this standard
        http://www.math.uwaterloo.ca/tsp/world/countries.html#DJ

    """

    node_coords = DelimitedFiles.readdlm(args[1],' ',Float32,'\n')

    """
    nxn Matrix of floats representing the distances between the nodes
    """
    node_distances = distances_from_coords(node_coords)

    node_number = size(node_distances)[1]

    tsp_model = Model(with_optimizer(Cbc.Optimizer , threads=4))

    #Limits execution time
    set_time_limit_sec(tsp_model , 7200)

    #Decision variable - binary nxn Matrix representing the edges taken between nodes
    @variable(tsp_model, edges_taken[1:node_number, 1:node_number], Bin)
	@variable(tsp_model,  0 <= u[1:node_number]  <= (node_number-1),Int)

    #Min Sum ( node_distances[i][j] * edges_taken[i][j] for all nodes)
    @objective(tsp_model,Min, sum( map(*,node_distances,edges_taken) ) )

    fix(u[1], 0,force=true) # change the 1 for the position of initial node

	for i in 1:node_number
		@constraint(tsp_model,sum(edges_taken[i,j] for j in 1:node_number if i != j ) == 1)
	end
	# Restrições para garantir a chegada em cada vértice
	for j in 1:node_number
		@constraint(tsp_model,sum(edges_taken[i,j] for i in 1:node_number if i != j) == 1)
	end
	for i in 1:node_number
		for j in 2:node_number
			if i != j
				@constraint(tsp_model,u[i] - u[j] + node_number * edges_taken[i,j] <= node_number - 1)
			end
		end
	end
	for i in 2:node_number
		@constraint(tsp_model, 0 <= u[i] <= node_number - 1)
	end

    optimize!(tsp_model)

	#Código para facilitar a impressão das cidades percorridas
	auxU = value.(u)
	ordem = zeros(UInt64 , (node_number,1))
	for i in 1:node_number
		ordem[Int64(round(auxU[i]+1))] = i
	end

	print("O valor do caminho otimo é: ")
	println(objective_value(tsp_model))
	print("A ordem dos pontos percorridos: ")
	for i in 1:node_number
		print(ordem[i])
		print(" ")
	end
end

tsp_main(ARGS)

max_stars_repo_name: jvlsg/julia-tsp

max_issues_repo_name: jvlsg/julia-tsp

max_forks_repo_name: jvlsg/julia-tsp


text:
using Test

const EXAMPLES = ["AJPR14e54.jl"]

@testset "run_examples.jl" begin
    @testset "$(example)" for example in EXAMPLES
        include(example)
    end
end


text:
function hasfieldcount(@nospecialize(dt))
    try
        fieldcount(dt)
    catch
        return false
    end
    return true
end


text:
@testset "Arp.jl" begin
    filename = "data\\test.arp"
    df = DDR2import.Arp.read(filename)
    @test df[1, :AIRPORT] == "EGSX"
    @test df[2, :LAT_DEG] ≈ 53.3941667 atol = 0.001
    @test df[3, :LON_DEG] ≈ -1.7 atol = 0.001
    @test df[4, :FIR] == "EGTT_FIR"
    #LAST TEST
    @test df[5, :AIRPORT] == "EGTB"
end


text:
using .JSort
import CSV
import DelimitedFiles: writedlm
using Polynomials

const E_MIN = 5  # E Threshold

function sortevent!(matrices, event::Event, parameters::Parameters)
    #=
    Algorithm:
    For each ΔE, look for an E in the associated E-detector. If no E is found,
    or multiple correlations are made, ignore this count.
    =#
    minievent = correlate(event)

    # increment!(matrices[:emul], length(minievents))
    # increment!(matrices[:Δemul], length(minievents))

    # If there were more than 1 correlated event
    # in the time window, it is impossible to tell which
    # LaBr corresponds to which SiRi detection
    isnothing(minievent) && return

    # Calibrate the events
    e = calibrate(minievent, :e, parameters)
    Δe = calibrate(minievent, :Δe, parameters)

    # Reject event if under threshold
    e ≤ E_MIN || Δe ≤ E_MIN && return

    # Look for events within the banana defined by Qkinz
    thick::Float64 = parameters.range(e+Δe) - parameters.range(e)
    thickdev::Float64 = parameters.thickrange[2] + e*parameters.thickrange[3]

    # Increment the matrices
    b, f = minievent.back, minievent.front
    increment!(matrices[:front], Δe, minievent.frontchannel)
    increment!(matrices[:back], e, minievent.backchannel)
    increment!(matrices[:eΔe], e, Δe)
    # increment!(matrices[Symbol("eΔeb",b, "f", f)], e, Δe)
    increment!(matrices[:thick], thick)

    # Check if the event comes within the banana
    within_banana = abs(thick-parameters.thickrange[1]) < thickdev
    !within_banana && return

    # Increment even more matrices
    increment!(matrices[:eΔethick], e, Δe)
    increment!(matrices[:heΔe], e+Δe)
    # increment!(matrices[Symbol("heΔerf", f)], e+Δe)
    # increment!(matrices[Symbol("heΔeb", b, "f", f)], e+Δe)

    ex = calibrate(e+Δe, minievent, :ex, parameters)
    increment!(matrices[:hex], ex)
    # increment!(matrices[Symbol("hexrf", f)], ex)

    # TODO Find a way to support time evolution
    # Time evolution

    # LaBr
    increment!(matrices[:labrmul], event.labrnum)
    for (channel, labr) in eachlabr(event)
        # Check that both adc and tdc have been set
        !islabrgood(labr) && return
        labre = calibrate(labr.adc, :labre, parameters)
        # TODO why this magic /8?
        labrtime = calibrate(labr.tdc/8, :labrtime, parameters)
        labrtime_c = tLaBr(labrtime, labre, e)

        increment!(matrices[:labre], labre, channel)
        increment!(matrcies[:labrtimec], channel)

    end
    # for (channel, labr) in event.labr
        # labr.adc ≤ 0 && continue
        # e = 0
    # end
    return
    # e = ecalibrator(channel_b, channel_f)(si_e[channel_b])
    # Do thick stuff here
end

function oldmakeeΔe(events::Array{TrueEvent, 1}, parameters::Parameters)
    meΔe = OMatrix{Array{Int64, 2}}((NUMBINS_E, NUMBINS_E), max=(MAX_E, MAX_ΔE))
    meΔebf = [OMatrix{Array{Int64, 2}}((NUMBINS_E, NUMBINS_E), max=(MAX_E, MAX_ΔE))
              for f in 1:8, b in 1:8]
    meΔebf_raw = [[] for f in 1:8, b in 1:8]

    for (i, event) in enumerate(events)
        e = calibrate(event.e, 0.0, 5.0)
        Δe = calibrate(event.Δe, 0.0, 2.5)
        push!(meΔebf_raw[event.b, event.f], (e, Δe))
        increment!(meΔe, e, Δe)
        increment!(meΔebf[event.b, event.f], e, Δe)
    end
    #save(meΔe, "mede", parameters.savepath)
    for b in 1:8, f in 1:8
        #save(meΔebf[b, f], "medeb$(b)f$f", parameters.savepath)
        path = joinpath(getpath(parameters, :savepath), "/medeb$(b)f$f.csv")
        CSV.write(path, meΔebf_raw[b, f])
    end
    meΔe, meΔebf
end

function makeeΔe(events::AbstractVector{TrueEvent}, parameters::Parameters; docalibrate=true)
    eΔe = [Tuple{Float32, Float32}[] for f in 1:8, b in 1:8]
    ecalibrators  = parameters.calibrator[:e]
    Δecalibrators = parameters.calibrator[:Δe]
    e = Δe = 0.0
    for (i, event) in enumerate(events)
        e  = ecalibrators[event.f, event.b](event.e)
        Δe = Δecalibrators[event.f, event.b](event.Δe)
        push!(eΔe[event.f, event.b], (e, Δe))
    end
    for f in 1:8
        total = Tuple{Float32, Float32}[]
        for b in 1:8
            path = joinpath(getpath(parameters, :savepath), "edeb$(b)f$(f).csv")
            # Stupid writedlm syntax
            data = []
            push!(data, ("e", "de"))
            append!(data, eΔe[f, b])
            open(path, "w") do io
                writedlm(io, data, ',')
            end
            append!(total, eΔe[f, b])
        end
        path = joinpath(getpath(parameters, :savepath), "edef$(f).csv")
        data = []
        push!(data, ("e", "de"))
        append!(data, total)
        open(path, "w") do io
            writedlm(io, data, ',')
        end
    end
end

function makeeΔebin(events::AbstractVector{TrueEvent}, parameters::Parameters; T=Float32)
    eΔe = [Tuple{T, T}[] for f in 1:8, b in 1:8]
    ecalibrators  = parameters.calibrator[:e]
    Δecalibrators = parameters.calibrator[:Δe]
    e = Δe = 0.0
    for (i, event) in enumerate(events)
        e  = ecalibrators[event.f, event.b](event.e)
        Δe = Δecalibrators[event.f, event.b](event.Δe)
        push!(eΔe[event.f, event.b], (e, Δe))
    end
    for f in 1:8
        path = joinpath(getpath(parameters, :savepath), "edef$(f).bin")
        iototal = open(path, "w")
        # Write the total length
        totallength = length.(eΔe[f, :]) |> sum |> T
        write(iototal, totallength)
        for b in 1:8
            path = joinpath(getpath(parameters, :savepath), "edeb$(b)f$(f).bin")
            data = eΔe[f, b]
            # Stupid writedlm syntax
            open(path, "w") do io
                write(io, length(data) |> T)
                write(io, data)
            end
            # Append to the summed data
            write(iototal, data)
        end
        close(iototal)
    end
end

function makelabr(events::AbstractVector{TrueEvent}, parameters::Parameters; docalibrate=true)
    # TODO: OMatrices are *almost* useless
    # Use long arrays instead
    #gamma_e = OMatrix{Array{Int64, 2}}((2000, 32), max=(14000, 32))
    #gamma_t = OMatrix{Array{Int64, 2}}((500, 32))
    channels = 1:33
    Γe = Dict(i => Float32[] for i in channels)
    Γt = Dict(i => Float32[] for i in channels)

    ecalibrators = parameters.calibrator[:γ]
    tcalibrators = parameters.calibrator[:t]
    @debug "Making LaBr gamma events"
    e = 0.0
    t = 0.0
    for event in events
        for γ in event.labr
            e = ecalibrators[γ.channel+1](γ.adc)
            t = tcalibrators[γ.channel+1](γ.tdc/8)
            #increment!(gamma_e, e, gamma.channel)
            push!(Γe[γ.channel+1], e)
            push!(Γt[γ.channel+1], t)
            #increment!(gamma_t, t, gamma.channel)
        end
    end
    #save(gamma_e, "gammae", getpath(parameters, :savepath))
    #save(gamma_t, "gammat", getpath(parameters, :savepath))
    return Γe, Γt
end

max_stars_repo_name: Caronthir/JSort

max_issues_repo_name: Caronthir/JSort

max_forks_repo_name: Caronthir/JSort


text:
# @copyright (c) 2019-2020 RWTH Aachen. All rights reserved.
#
# @file test/runtests.jl
# @version 1.0.0
# @author Aleksandr Mikhalev
# @date 2020-02-06

using Maxvol, LinearAlgebra, Random, Test

# Size and number of test matrices
nrows = 100
ncols = 20
ntests = 50

# Parameters for maxvol
tol = 1.05
maxiters = 100
# Set random seed
rng = Random.MersenneTwister(100)

# Lapack-based Maxvol
println("Testing LAPACK-based Maxvol")
types = [Float32, Float64, ComplexF32, ComplexF64]
for T in types
    println("    ", T)
    # Generate orthogonal/unitary input
    A = Matrix(qr(randn(rng, T, nrows, ncols)).Q)
    B = copy(A)
    piv, niters = Maxvol.maxvol!(B, tol, maxiters)
    @test A ≈ B*A[piv,:]
    # Lapack-based maxvol in complex case uses izamax(), which returns maximum
    # of |re|+|im|, not |re|^2+|im|^2. So we do not check Inf norm, but
    # absolute value of matrix entry with maximum 1-norm.
    if T <: AbstractFloat
        @test norm(B, Inf) <= tol
    else
        C = abs.(real(B)) + abs.(imag(B))
        maxv, maxi = findmax(C)
        @test abs(B[maxi]) <= tol
    end
end

# Lapack-based Maxvol does not exist for certain types
types = [Float16, BigFloat, ComplexF16, Complex{BigFloat}, Int]
for T in types
    println("    ", T)
    A = zeros(T, nrows, ncols)
    @test_throws MethodError Maxvol.maxvol!(A, tol, maxiters)
end

# Lapack-based Maxvol throws ArgumentError if input is singular matrix
types = [Float32, Float64, ComplexF32, ComplexF64]
for T in types
    println("    ", T)
    A = ones(T, nrows, ncols)
    @test_throws ArgumentError Maxvol.maxvol!(A, tol, maxiters)
end

# Generic Maxvol
println("Testing Generic Maxvol")
types = [Float16, Float32, Float64, BigFloat, ComplexF16, ComplexF32,
         ComplexF64, Complex{BigFloat}]
for T in types
    println("    ", T)
    A = Matrix(qr(rand(rng, T, nrows, ncols)).Q)
    B = copy(A)
    piv, niters = Maxvol.maxvol_generic!(B, tol, maxiters)
    @test A ≈ B*A[piv,:]
    # Check that matrix B is upperbounded
    @test norm(B, Inf) <= tol
end

# Generic Maxvol throws ArgumentError if input is singular matrix
types = [Float16, Float32, Float64, BigFloat, ComplexF16, ComplexF32,
         ComplexF64, Complex{BigFloat}]
for T in types
    println("    ", T)
    A = ones(T, nrows, ncols)
    @test_throws ArgumentError Maxvol.maxvol_generic!(A, tol, maxiters)
end

# Parameters for rect_maxvol
tol = 1.0
# Set random seed
rng = Random.MersenneTwister(100)

# Lapack-based Rect_maxvol
println("Testing LAPACK-based Rect_maxvol")
types = [Float32, Float64, ComplexF32, ComplexF64]
for T in types
    println("    ", T)
    # Generate orthogonal/unitary input
    A = Matrix(qr(randn(rng, T, nrows, ncols)).Q)
    piv, C = Maxvol.rect_maxvol(A, tol)
    @test A ≈ C*A[piv,:]
    # Check that norm of each row of C is less than tol
    @test all([norm(C[i,:]) <= tol for i=1:nrows])
end

# Lapack-based Rect_maxvol does not exist for certain types
types = [Float16, BigFloat, ComplexF16, Complex{BigFloat}, Int]
for T in types
    println("    ", T)
    A = zeros(T, nrows, ncols)
    @test_throws MethodError Maxvol.rect_maxvol(A, tol)
end

# Lapack-based Rect_maxvol throws ArgumentError if input is singular matrix
types = [Float32, Float64, ComplexF32, ComplexF64]
for T in types
    println("    ", T)
    A = ones(T, nrows, ncols)
    @test_throws ArgumentError Maxvol.rect_maxvol(A, tol)
end

# Generic Rect_maxvol
println("Testing Generic Rect_maxvol")
types = [Float16, Float32, Float64, BigFloat, ComplexF16, ComplexF32,
         ComplexF64, Complex{BigFloat}]
for T in types
    println("    ", T)
    A = Matrix(qr(rand(rng, T, nrows, ncols)).Q)
    piv, C = Maxvol.rect_maxvol_generic(A, tol)
    @test A ≈ C*A[piv,:]
    # Check that norm of each row of C is less than tol
    @test all([norm(C[i,:]) <= tol for i=1:nrows])
end

# Generic Rect_maxvol throws ArgumentError if input is singular matrix
types = [Float16, Float32, Float64, BigFloat, ComplexF16, ComplexF32,
         ComplexF64, Complex{BigFloat}]
for T in types
    println("    ", T)
    A = ones(T, nrows, ncols)
    @test_throws ArgumentError Maxvol.rect_maxvol_generic(A, tol)
end


text:
Failed to fetch UNSO time:
raw html:
 <!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 3.2 Final"//EN>

<TITLE>What time is it?</TITLE>

 US Naval Observatory Master Clock Time

Apr. 20, 17:55:31 UTC               Universal Time

Apr. 20, 01:55:31 PM EDT            Eastern Time

Apr. 20, 12:55:31 PM CDT            Central Time

Apr. 20, 11:55:31 AM MDT            Mountain Time

Apr. 20, 10:55:31 AM PDT            Pacific Time

Apr. 20, 09:55:31 AM AKDT   Alaska Time

Apr. 20, 07:55:31 AM HAST   Hawaii-Aleutian Time

 US Naval Observatory

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
struct ConditionalModel{A0,A,B,M} <: Distribution{MixedVariate, MixedSupport}
    model::Model{A,B,M}
    args::A0
end

(jd::ConditionalModel)(nt::NamedTuple) = jd.model(merge(jd.args, nt))

# function (m::Model)(nt::NamedTuple)
#     badargs = setdiff(keys(nt), variables(m))
#     isempty(badargs) || @error "Unused arguments $badargs"

#     m = predictive(m, keys(nt)...)
#     return ConditionalModel(m, nt)
# end

(m::Model)(;args...)= m((;args...))

(m::Model{A,B,M})(nt::NamedTuple) where {A,B,M} = ConditionalModel(m,nt)

function Base.show(io::IO, d :: ConditionalModel)
    m = d.model
    println(io, "Joint Distribution")
    print(io, "    Bound arguments: [")
    join(io, fieldnames(arguments(d)), ", ")
    println(io, "]")
    print(io, "    Variables: [")
    join(io, setdiff(toposort(m),arguments(m)), ", ")
    println(io, "]\n")
    println(io, convert(Expr, m))
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

# Base.require is the implementation for the `import` statement

# Cross-platform case-sensitive path canonicalization

if is_unix() && !is_apple()
    # assume case-sensitive filesystems, don't have to do anything
    isfile_casesensitive(path) = isfile(path)
elseif is_windows()
    # GetLongPathName Win32 function returns the case-preserved filename on NTFS.
    function isfile_casesensitive(path)
        isfile(path) || return false  # Fail fast
        basename(Filesystem.longpath(path)) == basename(path)
    end
elseif is_apple()
    # HFS+ filesystem is case-preserving. The getattrlist API returns
    # a case-preserved filename. In the rare event that HFS+ is operating
    # in case-sensitive mode, this will still work but will be redundant.

    # Constants from <sys/attr.h>
    const ATRATTR_BIT_MAP_COUNT = 5
    const ATTR_CMN_NAME = 1
    const BITMAPCOUNT = 1
    const COMMONATTR = 5
    const FSOPT_NOFOLLOW = 1  # Don't follow symbolic links

    const attr_list = zeros(UInt8, 24)
    attr_list[BITMAPCOUNT] = ATRATTR_BIT_MAP_COUNT
    attr_list[COMMONATTR] = ATTR_CMN_NAME

    # This essentially corresponds to the following C code:
    # attrlist attr_list;
    # memset(&attr_list, 0, sizeof(attr_list));
    # attr_list.bitmapcount = ATTR_BIT_MAP_COUNT;
    # attr_list.commonattr = ATTR_CMN_NAME;
    # struct Buffer {
    #    u_int32_t total_length;
    #    u_int32_t filename_offset;
    #    u_int32_t filename_length;
    #    char filename[max_filename_length];
    # };
    # Buffer buf;
    # getattrpath(path, &attr_list, &buf, sizeof(buf), FSOPT_NOFOLLOW);
    function isfile_casesensitive(path)
        isfile(path) || return false
        path_basename = String(basename(path))
        local casepreserved_basename
        const header_size = 12
        buf = Vector{UInt8}(length(path_basename) + header_size + 1)
        while true
            ret = ccall(:getattrlist, Cint,
                        (Cstring, Ptr{Void}, Ptr{Void}, Csize_t, Culong),
                        path, attr_list, buf, sizeof(buf), FSOPT_NOFOLLOW)
            systemerror(:getattrlist, ret ≠ 0)
            filename_length = unsafe_load(
              convert(Ptr{UInt32}, pointer(buf) + 8))
            if (filename_length + header_size) > length(buf)
                resize!(buf, filename_length + header_size)
                continue
            end
            casepreserved_basename =
              view(buf, (header_size+1):(header_size+filename_length-1))
            break
        end
        # Hack to compensate for inability to create a string from a subarray with no allocations.
        Vector{UInt8}(path_basename) == casepreserved_basename && return true

        # If there is no match, it's possible that the file does exist but HFS+
        # performed unicode normalization. See  https://developer.apple.com/library/mac/qa/qa1235/_index.html.
        isascii(path_basename) && return false
        Vector{UInt8}(normalize_string(path_basename, :NFD)) == casepreserved_basename
    end
else
    # Generic fallback that performs a slow directory listing.
    function isfile_casesensitive(path)
        isfile(path) || return false
        dir, filename = splitdir(path)
        any(readdir(dir) .== filename)
    end
end

function load_hook(prefix::String, name::String, ::Void)
    name_jl = "$name.jl"
    path = joinpath(prefix, name_jl)
    isfile_casesensitive(path) && return abspath(path)
    path = joinpath(prefix, name_jl, "src", name_jl)
    isfile_casesensitive(path) && return abspath(path)
    path = joinpath(prefix, name, "src", name_jl)
    isfile_casesensitive(path) && return abspath(path)
    return nothing
end
load_hook(prefix::String, name::String, path::String) = path
load_hook(prefix, name::String, ::Any) =
    throw(ArgumentError("unrecognized custom loader in LOAD_PATH: $prefix"))

_str(x::AbstractString) = String(x)
_str(x) = x

# `wd` is a working directory to search. defaults to current working directory.
# if `wd === nothing`, no extra path is searched.
function find_in_path(name::String, wd::Union{Void,String})
    isabspath(name) && return name
    base = name
    if endswith(name,".jl")
        base = name[1:end-3]
    else
        name = string(base,".jl")
    end
    if wd !== nothing
        isfile_casesensitive(joinpath(wd,name)) && return joinpath(wd,name)
    end
    path = nothing
    path = _str(load_hook(_str(Pkg.dir()), base, path))
    for dir in LOAD_PATH
        path = _str(load_hook(_str(dir), base, path))
    end
    return path
end
find_in_path(name::AbstractString, wd::AbstractString = pwd()) =
    find_in_path(String(name), String(wd))

function find_in_node_path(name::String, srcpath, node::Int=1)
    if myid() == node
        return find_in_path(name, srcpath)
    else
        return remotecall_fetch(find_in_path, node, name, srcpath)
    end
end

function find_source_file(file::String)
    (isabspath(file) || isfile(file)) && return file
    file2 = find_in_path(file)
    file2 !== nothing && return file2
    file2 = joinpath(JULIA_HOME, DATAROOTDIR, "julia", "base", file)
    return isfile(file2) ? file2 : nothing
end

function find_all_in_cache_path(mod::Symbol)
    name = string(mod)
    paths = String[]
    for prefix in LOAD_CACHE_PATH
        path = joinpath(prefix, name*".ji")
        if isfile_casesensitive(path)
            push!(paths, path)
        end
    end
    return paths
end

# these return either the array of modules loaded from the path / content given
# or an Exception that describes why it couldn't be loaded
function _include_from_serialized(content::Vector{UInt8})
    return ccall(:jl_restore_incremental_from_buf, Any, (Ptr{UInt8}, Int), content, sizeof(content))
end
function _include_from_serialized(path::String)
    return ccall(:jl_restore_incremental, Any, (Cstring,), path)
end

# returns an array of modules loaded, or an Exception that describes why it failed
# and also attempts to load the same file across all nodes (if toplevel_node and myid() == master)
# and it reconnects the Base.Docs.META
function _require_from_serialized(node::Int, mod::Symbol, path_to_try::String, toplevel_load::Bool)
    local restored = nothing
    local content::Vector{UInt8}
    if toplevel_load && myid() == 1 && nprocs() > 1
        # broadcast top-level import/using from node 1 (only)
        if node == myid()
            content = open(read, path_to_try)
        else
            content = remotecall_fetch(open, node, read, path_to_try)
        end
        restored = _include_from_serialized(content)
        isa(restored, Exception) && return restored

        results = sizehint!(Vector{Tuple{Int,Any}}(), nprocs())
        @sync for p in procs()
            if p != myid()
                @async begin
                    result = remotecall_fetch(p) do
                        let m = try
                                    _include_from_serialized(content)
                                catch ex
                                    isa(ex, Exception) ? ex : ErrorException(string(ex))
                                end
                            isa(m, Exception) ? m : nothing
                        end
                    end
                    push!(results, (p, result))
                end
            end
        end
        for (id, m) in results
            if m !== nothing
                warn("Node state is inconsistent: node $id failed to load cache from $path_to_try. Got:")
                warn(m, prefix="WARNING: ")
            end
        end
    elseif node == myid()
        restored = _include_from_serialized(path_to_try)
    else
        content = remotecall_fetch(open, node, read, path_to_try)
        restored = _include_from_serialized(content)
    end

    if !isa(restored, Exception)
        for M in restored::Vector{Any}
            if isdefined(M, Base.Docs.META)
                push!(Base.Docs.modules, M)
            end
        end
    end
    return restored
end

# returns `true` if require found a precompile cache for this mod, but couldn't load it
# returns `false` if the module isn't known to be precompilable
# returns the set of modules restored if the cache load succeeded
function _require_search_from_serialized(node::Int, mod::Symbol, sourcepath::String, toplevel_load::Bool)
    if node == myid()
        paths = find_all_in_cache_path(mod)
    else
        paths = @fetchfrom node find_all_in_cache_path(mod)
    end

    for path_to_try in paths::Vector{String}
        if node == myid()
            if stale_cachefile(sourcepath, path_to_try)
                continue
            end
        else
            if @fetchfrom node stale_cachefile(sourcepath, path_to_try)
                continue
            end
        end
        restored = _require_from_serialized(node, mod, path_to_try, toplevel_load)
        if isa(restored, Exception)
            if isa(restored, ErrorException) && endswith(restored.msg, " uuid did not match cache file.")
                # can't use this cache due to a module uuid mismatch,
                # defer reporting error until after trying all of the possible matches
                DEBUG_LOADING[] && info("JL_DEBUG_LOADING: Failed to load $path_to_try because $(restored.msg)")
                continue
            end
            warn("Deserialization checks failed while attempting to load cache from $path_to_try.")
            throw(restored)
        else
            return restored
        end
    end
    return !isempty(paths)
end

# this value is set by `require` based on whether JULIA_DEBUG_LOADING
# is presently defined as an environment variable
# and makes the logic in this file noisier about what it is doing and why
const DEBUG_LOADING = Ref(false)

# to synchronize multiple tasks trying to import/using something
const package_locks = Dict{Symbol,Condition}()

# to notify downstream consumers that a module was successfully loaded
# Callbacks take the form (mod::Symbol) -> nothing.
# WARNING: This is an experimental feature and might change later, without deprecation.
const package_callbacks = Any[]

# used to optionally track dependencies when requiring a module:
const _concrete_dependencies = Any[] # these dependency versions are "set in stone", and the process should try to avoid invalidating them
const _require_dependencies = Any[] # a list of (path, mtime) tuples that are the file dependencies of the module currently being precompiled
const _track_dependencies = Ref(false) # set this to true to track the list of file dependencies
function _include_dependency(_path::AbstractString)
    prev = source_path(nothing)
    path = (prev === nothing) ? abspath(_path) : joinpath(dirname(prev), _path)
    if myid() == 1 && _track_dependencies[]
        apath = abspath(path)
        push!(_require_dependencies, (apath, mtime(apath)))
    end
    return path, prev
end

"""
    include_dependency(path::AbstractString)

In a module, declare that the file specified by `path` (relative or absolute) is a
dependency for precompilation; that is, the module will need to be recompiled if this file
changes.

This is only needed if your module depends on a file that is not used via `include`. It has
no effect outside of compilation.
"""
function include_dependency(path::AbstractString)
    _include_dependency(path)
    return nothing
end

# We throw PrecompilableError(true) when a module wants to be precompiled but isn't,
# and PrecompilableError(false) when a module doesn't want to be precompiled but is
struct PrecompilableError <: Exception
    isprecompilable::Bool
end
function show(io::IO, ex::PrecompilableError)
    if ex.isprecompilable
        print(io, "Declaring __precompile__(true) is only allowed in module files being imported.")
    else
        print(io, "Declaring __precompile__(false) is not allowed in files that are being precompiled.")
    end
end
precompilableerror(ex::PrecompilableError, c) = ex.isprecompilable == c
precompilableerror(ex::WrappedException, c) = precompilableerror(ex.error, c)
precompilableerror(ex, c) = false

# Call __precompile__ at the top of a file to force it to be precompiled (true), or
# to be prevent it from being precompiled (false).  __precompile__(true) is
# ignored except within "require" call.
"""
    __precompile__(isprecompilable::Bool=true)

Specify whether the file calling this function is precompilable. If `isprecompilable` is
`true`, then `__precompile__` throws an exception when the file is loaded by
`using`/`import`/`require` *unless* the file is being precompiled, and in a module file it
causes the module to be automatically precompiled when it is imported. Typically,
`__precompile__()` should occur before the `module` declaration in the file.

If a module or file is *not* safely precompilable, it should call `__precompile__(false)` in
order to throw an error if Julia attempts to precompile it.

`__precompile__()` should *not* be used in a module unless all of its dependencies are also
using `__precompile__()`. Failure to do so can result in a runtime error when loading the module.
"""
function __precompile__(isprecompilable::Bool=true)
    if (myid() == 1 &&
        JLOptions().use_compilecache != 0 &&
        isprecompilable != (0 != ccall(:jl_generating_output, Cint, ())) &&
        !(isprecompilable && toplevel_load::Bool))
        throw(PrecompilableError(isprecompilable))
    end
end

function require_modname(name::AbstractString)
    # This function can be deleted when the deprecation for `require`
    # is deleted.
    # While we could also strip off the absolute path, the user may be
    # deliberately directing to a different file than what got
    # cached. So this takes a conservative approach.
    if Bool(JLOptions().use_compilecache)
        if endswith(name, ".jl")
            tmp = name[1:end-3]
            for prefix in LOAD_CACHE_PATH
                path = joinpath(prefix, tmp*".ji")
                if isfile(path)
                    return tmp
                end
            end
        end
    end
    return name
end

"""
    reload(name::AbstractString)

Force reloading of a package, even if it has been loaded before. This is intended for use
during package development as code is modified.
"""
function reload(name::AbstractString)
    if isfile(name) || contains(name,Filesystem.path_separator)
        # for reload("path/file.jl") just ask for include instead
        error("use `include` instead of `reload` to load source files")
    else
        # reload("Package") is ok
        require(Symbol(require_modname(name)))
    end
end

# require always works in Main scope and loads files from node 1
toplevel_load = true

"""
    require(module::Symbol)

This function is part of the implementation of `using` / `import`, if a module is not
already defined in `Main`. It can also be called directly to force reloading a module,
regardless of whether it has been loaded before (for example, when interactively developing
libraries).

Loads a source file, in the context of the `Main` module, on every active node, searching
standard locations for files. `require` is considered a top-level operation, so it sets the
current `include` path but does not use it to search for files (see help for `include`).
This function is typically used to load library code, and is implicitly called by `using` to
load packages.

When searching for files, `require` first looks for package code under `Pkg.dir()`,
then tries paths in the global array `LOAD_PATH`. `require` is case-sensitive on
all platforms, including those with case-insensitive filesystems like macOS and
Windows.
"""
function require(mod::Symbol)
    _require(mod::Symbol)
    # After successfully loading notify downstream consumers
    for callback in package_callbacks
        invokelatest(callback, mod)
    end
end

function _require(mod::Symbol)
    # dependency-tracking is only used for one top-level include(path),
    # and is not applied recursively to imported modules:
    old_track_dependencies = _track_dependencies[]
    _track_dependencies[] = false
    DEBUG_LOADING[] = haskey(ENV, "JULIA_DEBUG_LOADING")

    global toplevel_load
    loading = get(package_locks, mod, false)
    if loading !== false
        # load already in progress for this module
        wait(loading)
        return
    end
    package_locks[mod] = Condition()

    last = toplevel_load::Bool
    try
        toplevel_load = false
        # perform the search operation to select the module file require intends to load
        name = string(mod)
        path = find_in_node_path(name, nothing, 1)
        if path === nothing
            throw(ArgumentError("Module $name not found in current path.\nRun `Pkg.add(\"$name\")` to install the $name package."))
        end

        # attempt to load the module file via the precompile cache locations
        doneprecompile = false
        if JLOptions().use_compilecache != 0
            doneprecompile = _require_search_from_serialized(1, mod, path, last)
            if !isa(doneprecompile, Bool)
                return # success
            end
        end

        # if the module being required was supposed to have a particular version
        # but it was not handled by the precompile loader, complain
        for (concrete_mod, concrete_uuid) in _concrete_dependencies
            if mod === concrete_mod
                warn("""Module $mod with uuid $concrete_uuid is missing from the cache.
                     This may mean module $mod does not support precompilation but is imported by a module that does.""")
                if JLOptions().incremental != 0
                    # during incremental precompilation, this should be fail-fast
                    throw(PrecompilableError(false))
                end
            end
        end

        if doneprecompile === true || JLOptions().incremental != 0
            # spawn off a new incremental pre-compile task from node 1 for recursive `require` calls
            # or if the require search declared it was pre-compiled before (and therefore is expected to still be pre-compilable)
            cachefile = compilecache(mod)
            m = _require_from_serialized(1, mod, cachefile, last)
            if isa(m, Exception)
                warn("The call to compilecache failed to create a usable precompiled cache file for module $name. Got:")
                warn(m, prefix="WARNING: ")
                # fall-through, TODO: disable __precompile__(true) error so that the normal include will succeed
            else
                return # success
            end
        end

        # just load the file normally via include
        # for unknown dependencies
        try
            # include on node 1 first to check for PrecompilableErrors
            Base.include_from_node1(Main, path)

            if last && myid() == 1 && nprocs() > 1
                # broadcast top-level import/using from node 1 (only)
                @sync begin
                    for p in filter(x -> x != 1, procs())
                        @async remotecall_fetch(p) do
                            Base.include_from_node1(Main, path)
                            nothing
                        end
                    end
                end
            end
        catch ex
            if doneprecompile === true || JLOptions().use_compilecache == 0 || !precompilableerror(ex, true)
                rethrow() # rethrow non-precompilable=true errors
            end
            # the file requested `__precompile__`, so try to build a cache file and use that
            cachefile = compilecache(mod)
            m = _require_from_serialized(1, mod, cachefile, last)
            if isa(m, Exception)
                warn(m, prefix="WARNING: ")
                # TODO: disable __precompile__(true) error and do normal include instead of error
                error("Module $mod declares __precompile__(true) but require failed to create a usable precompiled cache file.")
            end
        end
    finally
        toplevel_load = last
        loading = pop!(package_locks, mod)
        notify(loading, all=true)
        _track_dependencies[] = old_track_dependencies
    end
    nothing
end

# remote/parallel load

"""
    include_string(m::Module, code::AbstractString, filename::AbstractString="string")

Like `include`, except reads code from the given string rather than from a file. Since there
is no file path involved, no path processing or fetching from node 1 is done.
"""
include_string(m::Module, txt::String, fname::String) =
    ccall(:jl_load_file_string, Any, (Ptr{UInt8}, Csize_t, Cstring, Any),
          txt, sizeof(txt), fname, m)

include_string(m::Module, txt::AbstractString, fname::AbstractString="string") =
    include_string(m, String(txt), String(fname))

function source_path(default::Union{AbstractString,Void}="")
    t = current_task()
    while true
        s = t.storage
        if s !== nothing && haskey(s, :SOURCE_PATH)
            return s[:SOURCE_PATH]
        end
        if t === t.parent
            return default
        end
        t = t.parent
    end
end

function source_dir()
    p = source_path(nothing)
    p === nothing ? pwd() : dirname(p)
end

include_from_node1(mod::Module, path::AbstractString) = include_from_node1(mod, String(path))
function include_from_node1(mod::Module, _path::String)
    path, prev = _include_dependency(_path)
    tls = task_local_storage()
    tls[:SOURCE_PATH] = path
    local result
    try
        if myid()==1
            # sleep a bit to process file requests from other nodes
            nprocs()>1 && sleep(0.005)
            result = Core.include(mod, path)
            nprocs()>1 && sleep(0.005)
        else
            result = include_string(mod, remotecall_fetch(readstring, 1, path), path)
        end
    finally
        if prev === nothing
            delete!(tls, :SOURCE_PATH)
        else
            tls[:SOURCE_PATH] = prev
        end
    end
    result
end

"""
    include(m::Module, path::AbstractString)

Evaluate the contents of the input source file into module `m`. Returns the result
of the last evaluated expression of the input file. During including, a task-local include
path is set to the directory containing the file. Nested calls to `include` will search
relative to that path. All paths refer to files on node 1 when running in parallel, and
files will be fetched from node 1. This function is typically used to load source
interactively, or to combine files in packages that are broken into multiple source files.
"""
include # defined in sysimg.jl

"""
    evalfile(path::AbstractString, args::Vector{String}=String[])

Load the file using [`include`](@ref), evaluate all expressions,
and return the value of the last one.
"""
function evalfile(path::AbstractString, args::Vector{String}=String[])
    return eval(Module(:__anon__),
                Expr(:toplevel,
                     :(const ARGS = $args),
                     :(eval(x) = $(Expr(:core, :eval))(__anon__, x)),
                     :(eval(m, x) = $(Expr(:core, :eval))(m, x)),
                     :(include(x) = $(Expr(:top, :include))(__anon__, x)),
                     :(include($path))))
end
evalfile(path::AbstractString, args::Vector) = evalfile(path, String[args...])

function create_expr_cache(input::String, output::String, concrete_deps::Vector{Any})
    rm(output, force=true)   # Remove file if it exists
    code_object = """
        while !eof(STDIN)
            eval(Main, deserialize(STDIN))
        end
        """
    io = open(pipeline(detach(`$(julia_cmd()) -O0
                              --output-ji $output --output-incremental=yes
                              --startup-file=no --history-file=no
                              --color=$(have_color ? "yes" : "no")
                              --eval $code_object`), stderr=STDERR),
              "w", STDOUT)
    in = io.in
    try
        serialize(in, quote
                  empty!(Base.LOAD_PATH)
                  append!(Base.LOAD_PATH, $LOAD_PATH)
                  empty!(Base.LOAD_CACHE_PATH)
                  append!(Base.LOAD_CACHE_PATH, $LOAD_CACHE_PATH)
                  empty!(Base.DL_LOAD_PATH)
                  append!(Base.DL_LOAD_PATH, $DL_LOAD_PATH)
                  empty!(Base._concrete_dependencies)
                  append!(Base._concrete_dependencies, $concrete_deps)
                  Base._track_dependencies[] = true
                  end)
        source = source_path(nothing)
        if source !== nothing
            serialize(in, quote
                      task_local_storage()[:SOURCE_PATH] = $(source)
                      end)
        end
        serialize(in, :(Base.include(Main, $(abspath(input)))))
        if source !== nothing
            serialize(in, :(delete!(task_local_storage(), :SOURCE_PATH)))
        end
        close(in)
    catch ex
        close(in)
        process_running(io) && Timer(t -> kill(io), 5.0) # wait a short time before killing the process to give it a chance to clean up on its own first
        rethrow(ex)
    end
    return io
end

compilecache(mod::Symbol) = compilecache(string(mod))

"""
    Base.compilecache(module::String)

Creates a precompiled cache file for
a module and all of its dependencies.
This can be used to reduce package load times. Cache files are stored in
`LOAD_CACHE_PATH[1]`, which defaults to `~/.julia/lib/VERSION`. See
[Module initialization and precompilation](@ref)
for important notes.
"""
function compilecache(name::String)
    myid() == 1 || error("can only precompile from node 1")
    # decide where to get the source file from
    path = find_in_path(name, nothing)
    path === nothing && throw(ArgumentError("$name not found in path"))
    path = String(path)
    # decide where to put the resulting cache file
    cachepath = LOAD_CACHE_PATH[1]
    if !isdir(cachepath)
        mkpath(cachepath)
    end
    cachefile::String = abspath(cachepath, name*".ji")
    # build up the list of modules that we want the precompile process to preserve
    concrete_deps = copy(_concrete_dependencies)
    for existing in names(Main)
        if isdefined(Main, existing)
            mod = getfield(Main, existing)
            if isa(mod, Module) && !(mod === Main || mod === Core || mod === Base)
                mod = mod::Module
                if module_parent(mod) === Main && module_name(mod) === existing
                    push!(concrete_deps, (existing, module_uuid(mod)))
                end
            end
        end
    end
    # run the expression and cache the result
    if isinteractive() || DEBUG_LOADING[]
        if isfile(cachefile)
            info("Recompiling stale cache file $cachefile for module $name.")
        else
            info("Precompiling module $name.")
        end
    end
    if success(create_expr_cache(path, cachefile, concrete_deps))
        # append checksum to the end of the .ji file:
        open(cachefile, "a+") do f
            data = Mmap.mmap(f, Vector{UInt8}, filesize(f), 0)
            checksum = crc32c(data)
            finalize(data)
            write(f, hton(checksum))
        end
    else
        error("Failed to precompile $name to $cachefile.")
    end
    return cachefile
end

module_uuid(m::Module) = ccall(:jl_module_uuid, UInt64, (Any,), m)

isvalid_cache_header(f::IOStream) = 0 != ccall(:jl_read_verify_header, Cint, (Ptr{Void},), f.ios)

function parse_cache_header(f::IO)
    modules = Dict{Symbol,UInt64}()
    while true
        n = ntoh(read(f, Int32))
        n == 0 && break
        sym = Symbol(read(f, n)) # module symbol
        uuid = ntoh(read(f, UInt64)) # module UUID (mostly just a timestamp)
        modules[sym] = uuid
    end
    totbytes = ntoh(read(f, Int64)) # total bytes for file dependencies
    # read the list of files
    files = Tuple{String,Float64}[]
    while true
        n = ntoh(read(f, Int32))
        n == 0 && break
        totbytes -= 4 + n + 8
        @assert n >= 0 "EOF while reading cache header" # probably means this wasn't a valid file to be read by Base.parse_cache_header
        push!(files, (String(read(f, n)), ntoh(read(f, Float64))))
    end
    @assert totbytes == 4 "header of cache file appears to be corrupt"
    # read the list of modules that are required to be present during loading
    required_modules = Dict{Symbol,UInt64}()
    while true
        n = ntoh(read(f, Int32))
        n == 0 && break
        sym = Symbol(read(f, n)) # module symbol
        uuid = ntoh(read(f, UInt64)) # module UUID
        required_modules[sym] = uuid
    end
    return modules, files, required_modules
end

function parse_cache_header(cachefile::String)
    io = open(cachefile, "r")
    try
        !isvalid_cache_header(io) && throw(ArgumentError("Invalid header in cache file $cachefile."))
        return parse_cache_header(io)
    finally
        close(io)
    end
end

function cache_dependencies(f::IO)
    defs, files, modules = parse_cache_header(f)
    return modules, files
end

function cache_dependencies(cachefile::String)
    io = open(cachefile, "r")
    try
        !isvalid_cache_header(io) && throw(ArgumentError("Invalid header in cache file $cachefile."))
        return cache_dependencies(io)
    finally
        close(io)
    end
end

function stale_cachefile(modpath::String, cachefile::String)
    io = open(cachefile, "r")
    try
        if !isvalid_cache_header(io)
            DEBUG_LOADING[] && info("JL_DEBUG_LOADING: Rejecting cache file $cachefile due to it containing an invalid cache header.")
            return true # invalid cache file
        end
        modules, files, required_modules = parse_cache_header(io)

        # Check if transitive dependencies can be fullfilled
        for mod in keys(required_modules)
            if mod == :Main || mod == :Core || mod == :Base
                continue
            # Module is already loaded
            elseif isbindingresolved(Main, mod)
                continue
            end
            name = string(mod)
            path = find_in_node_path(name, nothing, 1)
            if path === nothing
                return true # Won't be able to fullfill dependency
            end
        end

        # check if this file is going to provide one of our concrete dependencies
        # or if it provides a version that conflicts with our concrete dependencies
        # or neither
        for (mod, uuid_req) in _concrete_dependencies
            uuid = get(modules, mod, UInt64(0))
            if uuid !== UInt64(0)
                if uuid === uuid_req
                    return false # this is the file we want
                end
                DEBUG_LOADING[] && info("JL_DEBUG_LOADING: Rejecting cache file $cachefile because it provides the wrong uuid (got $uuid) for $mod (want $uuid_req).")
                return true # cachefile doesn't provide the required version of the dependency
            end
        end

        # now check if this file is fresh relative to its source files
        if !samefile(files[1][1], modpath)
            DEBUG_LOADING[] && info("JL_DEBUG_LOADING: Rejecting cache file $cachefile because it is for file $(files[1][1])) not file $modpath.")
            return true # cache file was compiled from a different path
        end
        for (f, ftime_req) in files
            # Issue #13606: compensate for Docker images rounding mtimes
            # Issue #20837: compensate for GlusterFS truncating mtimes to microseconds
            ftime = mtime(f)
            if ftime != ftime_req && ftime != floor(ftime_req) && ftime != trunc(ftime_req, 6)
                DEBUG_LOADING[] && info("JL_DEBUG_LOADING: Rejecting stale cache file $cachefile (mtime $ftime_req) because file $f (mtime $ftime) has changed.")
                return true
            end
        end

        # finally, verify that the cache file has a valid checksum
        data = Mmap.mmap(io, Vector{UInt8}, filesize(io), 0)
        # checksum = UInt32 read in bigendian format from the last 4 bytes:
        checksum = UInt32(data[end]) + UInt32(data[end-1])<<8 + UInt32(data[end-2])<<16 + UInt32(data[end-3])<<24
        crc = crc32c(@view(data[1:end-4]))
        finalize(data)
        if checksum != crc
            DEBUG_LOADING[] && info("JL_DEBUG_LOADING: Rejecting cache file $cachefile because it has an invalid checksum.")
            return true
        end

        return false # fresh cachefile
    finally
        close(io)
    end
end

"""
    @__LINE__ -> Int

`@__LINE__` expands to the line number of the location of the macrocall.
Returns `0` if the line number could not be determined.
"""
macro __LINE__()
    return __source__.line
end

"""
    @__FILE__ -> AbstractString

`@__FILE__` expands to a string with the path to the file containing the
macrocall, or an empty string if evaluated by `julia -e <expr>`.
Returns `nothing` if the macro was missing parser source information.
Alternatively see [`PROGRAM_FILE`](@ref).
"""
macro __FILE__()
    __source__.file === nothing && return nothing
    return String(__source__.file)
end

"""
    @__DIR__ -> AbstractString

`@__DIR__` expands to a string with the absolute path to the directory of the file
containing the macrocall.
Returns the current working directory if run from a REPL or if evaluated by `julia -e <expr>`.
"""
macro __DIR__()
    __source__.file === nothing && return nothing
    return abspath(dirname(String(__source__.file)))
end

max_stars_repo_name: caodi/julia

max_issues_repo_name: caodi/julia

max_forks_repo_name: caodi/julia


text:
#  Copyright 2022 California Institute of Technology
#
#  Licensed under the Apache License, Version 2.0 (the "License");
#  you may not use this file except in compliance with the License.
#  You may obtain a copy of the License at
#
#      http://www.apache.org/licenses/LICENSE-2.0
#
#  Unless required by applicable law or agreed to in writing, software
#  distributed under the License is distributed on an "AS IS" BASIS,
#  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#  See the License for the specific language governing permissions and
#  limitations under the License.
#
# Author: Philip G. Brodrick, philip.g.brodrick@jpl.nasa.gov

using JuMP
using NLopt
using LinearAlgebra

function opt_solve(A::Matrix{Float64}, b::Vector{Float64}, x0::Vector{Float64}, lb::Vector{Float64}, ub::Vector{Float64})

    #mle = Model(optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0))
    mle = Model(NLopt.Optimizer)
    set_optimizer_attribute(mle, "algorithm", :LD_SLSQP) #LD_SLSQP

    x0[x0 .< 0] .= 0
    x0[x0 .> 1] .= 1

    @variable(mle, lb[i] <= x[i=1:length(x0)] <= ub[i])
    #@constraint(mle, sum(x) == 1)
    for n in 1:length(x0) set_start_value(x[n], x0[n]) end

    #@NLexpression(mle, mse ,  sum( (b .- sum(x[i] .* A[:,i] for i in 1:length(x) )).^2    ) )
    @expression(mle, mse ,  sum( (b .- sum(x[i] .* A[:,i] for i in 1:length(x) )).^2    ) )

    @NLobjective(mle, Min, mse)

    JuMP.optimize!(mle)
    return value.(mle[:x]) , exp(objective_value(mle))
end

function dolsq(A, b; method::String="default")
    if method == "default"
        x = A \ b
    elseif method == "pinv"
        x = pinv(A)*b
    elseif method == "qr"
        #Q,R = qr(A)
        #x = inv(R)*Q'*b
        qrA = qr(A)
        x = qrA\b
    end
    return x
end

function bvls(A, b, x_lsq, lb, ub, tol::Float64, max_iter::Int64, verbose::Int64, inverse_method::String)
    n_iter = 0
    m, n = size(A)

    x = x_lsq
    on_bound = zeros(n)

    mask = x .<= lb
    x[mask] = lb[mask]
    on_bound[mask] .= -1

    mask = x .>= ub
    x[mask] = ub[mask]
    on_bound[mask] .= 1

    free_set = on_bound .== 0
    active_set = .!free_set
    free_set = (1:size(free_set)[1])[free_set .!= 0]

    r = A * x - b
    cost = 0.5 * dot(r , r)
    initial_cost = cost
    g = A' * r

    cost_change = nothing
    step_norm = nothing
    iteration = 0

    while size(free_set)[1] > 0
        if verbose == 2
            optimality = compute_kkt_optimality(g, on_bound)
        end

        iteration += 1
        x_free_old = x[free_set]

        A_free = A[:, free_set]
        b_free = b - A * (x .* active_set)
        z = dolsq(A_free, b_free, method=inverse_method)

        lbv = z .< lb[free_set]
        ubv = z .> ub[free_set]

        v = lbv .| ubv

        if any(lbv)
            ind = free_set[lbv]
            x[ind] = lb[ind]
            active_set[ind] .= true
            on_bound[ind] .= -1
        end

        if any(ubv)
            ind = free_set[ubv]
            x[ind] = ub[ind]
            active_set[ind] .= true
            on_bound[ind] .= 1
        end

        ind = free_set[.!v]
        x[ind] = z[.!v]

        r = A * x -b
        cost_new = 0.5 * dot(r , r)
        cost_change = cost - cost_new
        cost = cost_new
        g = A' * r
        step_norm = sum((x[free_set] .- x_free_old).^2)

        if any(v)
            free_set = free_set[.!v]
        else
            break
        end
    end

    if isnothing(max_iter)
        max_iter = n
    end
    max_iter += iteration

    termination_status = nothing

    optimality = compute_kkt_optimality(g, on_bound)
    for iteration in iteration:max_iter
        if optimality < tol
            termination_status = 1
        end

        if !isnothing(termination_status)
            break
        end

        move_to_free = argmax(g .* on_bound)
        on_bound[move_to_free] = 0

        x_free = copy(x)
        x_free_old = copy(x)
        while true

            free_set = on_bound .== 0
            sum(free_set)
            active_set = .!free_set
            free_set = (1:size(free_set)[1])[free_set .!= 0]

            x_free = x[free_set]
            x_free_old = copy(x_free)
            lb_free = lb[free_set]
            ub_free = ub[free_set]

            A_free = A[:, free_set]
            b_free = b - A * (x .* active_set)
            z = dolsq(A_free, b_free, method=inverse_method)

            lbv = (1:size(free_set)[1])[ z .< lb_free]
            ubv = (1:size(free_set)[1])[ z .> ub_free]
            v = cat(lbv, ubv, dims=1)

            if size(v)[1] > 0
                alphas = cat(lb_free[lbv] - x_free[lbv], ub_free[ubv] - x_free[ubv],dims=1) ./ (z[v] - x_free[v])

                i = argmin(alphas)
                i_free = v[i]
                alpha = alphas[i]

                x_free .*= (1 .- alpha)
                x_free .+= (alpha .* z)
                x[free_set] = x_free

                vsize = size(lbv)
                if i <= size(lbv)[1]
                    on_bound[free_set[i_free]] = -1
                else
                    on_bound[free_set[i_free]] = 1
                end
            else
                x_free = z
                x[free_set] = x_free
                @goto start
            end
        end #while
        @label start
        step_norm = sum((x_free .- x_free_old).^2)

        r = A * x - b
        cost_new = 0.5 * dot(r , r)
        cost_change = cost - cost_new

        combo = tol * cost
        if cost_change < tol * cost
            termination_status = 2
        end
        cost = cost_new

        g = A' * r
        optimality = compute_kkt_optimality(g, on_bound)
    end #iteration

    if isnothing(termination_status)
        termination_status = 0
    end

    x[x .< 1e-5] .= 0
    return x, cost
end

function compute_kkt_optimality(g::Vector{Float64}, on_bound::Vector)
  g_kkt = g .* on_bound
  free_set = on_bound .== 0
  g_kkt[free_set] = broadcast(abs, g[free_set])

  return maximum(g_kkt)

end

max_stars_repo_name: fotxoa-geo/SpectralUnmixing

max_issues_repo_name: fotxoa-geo/SpectralUnmixing

max_forks_repo_name: fotxoa-geo/SpectralUnmixing


text:
module Core

import ..MarketData
import ..Currencies
import ..BusinessDays
import ..InterestRates
import ..CAL_BRL
import ..ForwardDiff
import ..OptimizingIR
import ..BS

using Dates
using Printf

include("types.jl")
include("observables_algebra.jl")
include("exch.jl")
include("contract_attributes.jl")
include("fixed_income.jl")
include("events.jl")
include("builtin_contracts.jl")
include("scenario.jl")
include("pricer.jl")
include("print.jl")
include("project_cashflows.jl")
include("compiler/compiler.jl")
include("docs.jl")

end # module Core


text:
using GeometricInterpretation
using Documenter

DocMeta.setdocmeta!(GeometricInterpretation, :DocTestSetup, :(using GeometricInterpretation); recursive=true)

makedocs(;
    modules=[GeometricInterpretation],
    authors="Shalin Patel <shalin.kp.patel@gmail.com> and contributors",
    repo="https://github.com/shalinkpatel/GeometricInterpretation.jl/blob/{commit}{path}#{line}",
    sitename="GeometricInterpretation.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://shalinkpatel.github.io/GeometricInterpretation.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/shalinkpatel/GeometricInterpretation.jl",
    devbranch = "main"
)


text:
include("Junction.jl")
include("Arc.jl")
include("Pattern.jl")
include("Curve.jl")
include("EPANETPumpParams.jl")
include("NormPumpParams.jl")
include("PumpParams.jl")
include("Pump.jl")
include("OpenPipe.jl")
include("GatePipe.jl")
include("CVPipe.jl")
include("StaticDemand.jl")
include("CylindricalTank.jl")
include("Reservoir.jl")

export get__forecasts
export get_arc
export get_available
export get_base_demand
export get_base_price
export get_coordinates
export get_diameter
export get_effncy_bep
export get_elevation
export get_epnt_data
export get_epnt_efficiency
export get_epnt_head
export get_epnt_power
export get_epnt_type
export get_flow
export get_flow_bep
export get_flowlimits
export get_from
export get_head
export get_head_bep
export get_headgain
export get_headloss
export get_internal
export get_junction
export get_length
export get_level
export get_level_limits
export get_multipliers
export get_name
export get_norm_coefs
export get_open_status
export get_operating
export get_pattern_name
export get_points
export get_power_bep
export get_pumpparams
export get_roughness
export get_to
export get_type


text:
using Xpress

const XPR_OPTIMIZER = Xpress.Optimizer()
if !Sys.iswindows()
    MOI.set(XPR_OPTIMIZER, MOI.Silent(), true)
end
const XPR_CACHE = MOIU.UniversalFallback(MOIU.Model{Float64}())
const XPR_CACHED = MOIU.CachingOptimizer(XPR_CACHE, XPR_OPTIMIZER)
const XPR_BRIDGED = MOI.Bridges.full_bridge_optimizer(XPR_CACHED, Float64)

push!(solvers, (opt = XPR_BRIDGED, mode = BilevelJuMP.PositiveSOS1Mode()))
push!(solvers_sos, (opt = XPR_BRIDGED, mode = BilevelJuMP.PositiveSOS1Mode()))
push!(solvers_quad, (opt = XPR_BRIDGED, mode = BilevelJuMP.PositiveSOS1Mode()))
push!(solvers_sos_quad, (opt = XPR_BRIDGED, mode = BilevelJuMP.PositiveSOS1Mode()))


text:
module Missings

export allowmissing, disallowmissing, ismissing, missing, missings,
       Missing, MissingException, levels, coalesce, passmissing, nonmissingtype

using Base: ismissing, missing, Missing, MissingException
using Base: @deprecate

@static if VERSION < v"1.3.0-alpha.121"
    nonmissingtype(::Type{S}) where {S} = Core.Compiler.typesubtract(S, Missing)
else
    using Base: nonmissingtype
end

import DataAPI

@deprecate T nonmissingtype false

# vector constructors
missings(dims::Dims) = fill(missing, dims)
missings(::Type{T}, dims::Dims) where {T >: Missing} = fill!(Array{T}(undef, dims), missing)
missings(::Type{T}, dims::Dims) where {T} = fill!(Array{Union{T, Missing}}(undef, dims), missing)
missings(dims::Integer...) = missings(dims)
missings(::Type{T}, dims::Integer...) where {T} = missings(T, dims)

"""
    allowmissing(x::AbstractArray)

Return an array equal to `x` allowing for [`missing`](@ref) values,
i.e. with an element type equal to `Union{eltype(x), Missing}`.

When possible, the result will share memory with `x` (as with [`convert`](@ref)).

See also: [`disallowmissing`](@ref)
"""
allowmissing(x::AbstractArray{T}) where {T} = convert(AbstractArray{Union{T, Missing}}, x)

"""
    disallowmissing(x::AbstractArray)

Return an array equal to `x` not allowing for [`missing`](@ref) values,
i.e. with an element type equal to `nonmissingtype(eltype(x))`.

When possible, the result will share memory with `x` (as with [`convert`](@ref)).
If `x` contains missing values, a `MethodError` is thrown.

See also: [`allowmissing`](@ref)
"""
disallowmissing(x::AbstractArray{T}) where {T} = convert(AbstractArray{nonmissingtype(T)}, x)

# Iterators
"""
    Missings.replace(itr, replacement)

Return an iterator wrapping iterable `itr` which replaces [`missing`](@ref) values with
`replacement`. When applicable, the size of `itr` is preserved.
If the type of `replacement` differs from the element type of `itr`,
it will be converted.

See also: [`skipmissing`](@ref), [`Missings.fail`](@ref)

# Examples
```jldoctest
julia> collect(Missings.replace([1, missing, 2], 0))
3-element Array{Int64,1}:
 1
 0
 2

julia> collect(Missings.replace([1 missing; 2 missing], 0))
2×2 Array{Int64,2}:
 1  0
 2  0

```
"""
replace(itr, replacement) = EachReplaceMissing(itr, convert(eltype(itr), replacement))
struct EachReplaceMissing{T, U}
    x::T
    replacement::U
end
Base.IteratorSize(::Type{<:EachReplaceMissing{T}}) where {T} =
    Base.IteratorSize(T)
Base.IteratorEltype(::Type{<:EachReplaceMissing{T}}) where {T} =
    Base.IteratorEltype(T)
Base.length(itr::EachReplaceMissing) = length(itr.x)
Base.size(itr::EachReplaceMissing) = size(itr.x)
Base.eltype(itr::EachReplaceMissing) = nonmissingtype(eltype(itr.x))

@inline function Base.iterate(itr::EachReplaceMissing)
    st = iterate(itr.x)
    st === nothing && return nothing
    v, s = st
    return (v isa Missing ? itr.replacement : v, s)
end

@inline function Base.iterate(itr::EachReplaceMissing, state)
    st = iterate(itr.x, state)
    st === nothing && return nothing
    v, s = st
    return (v isa Missing ? itr.replacement : v, s)
end

"""
    Missings.fail(itr)

Return an iterator wrapping iterable `itr` which will throw a [`MissingException`](@ref)
if a [`missing`](@ref) value is found.

Use [`collect`](@ref) to obtain an `Array` containing the resulting values.
If `itr` is an array, the resulting array will have the same dimensions.

See also: [`skipmissing`](@ref), [`Missings.replace`](@ref)

# Examples
```jldoctest
julia> collect(Missings.fail([1 2; 3 4]))
2×2 Array{Int64,2}:
 1  2
 3  4

julia> collect(Missings.fail([1, missing, 2]))
ERROR: MissingException: missing value encountered by Missings.fail
[...]
```
"""
fail(itr) = EachFailMissing(itr)
struct EachFailMissing{T}
    x::T
end
Base.IteratorSize(::Type{EachFailMissing{T}}) where {T} =
    Base.IteratorSize(T)
Base.IteratorEltype(::Type{EachFailMissing{T}}) where {T} =
    Base.IteratorEltype(T)
Base.length(itr::EachFailMissing) = length(itr.x)
Base.size(itr::EachFailMissing) = size(itr.x)
Base.eltype(itr::EachFailMissing) = nonmissingtype(eltype(itr.x))

@inline function Base.iterate(itr::EachFailMissing)
    st = iterate(itr.x)
    st === nothing && return nothing
    v, s = st
    ismissing(v) && throw(MissingException("missing value encountered by Missings.fail"))
    return (v::eltype(itr), s)
end

@inline function Base.iterate(itr::EachFailMissing, state)
    st = iterate(itr.x, state)
    st === nothing && return nothing
    v, s = st
    ismissing(v) && throw(MissingException("missing value encountered by Missings.fail"))
    return (v::eltype(itr), s)
end

const levels = DataAPI.levels

struct PassMissing{F} <: Function
    f::F
end

function (f::PassMissing{F})(x) where {F}
    if @generated
        return x === Missing ? missing : :(f.f(x))
    else
        return x === missing ? missing : f.f(x)
    end
end

function (f::PassMissing{F})(xs...) where {F}
    if @generated
        for T in xs
            T === Missing && return missing
        end
        return :(f.f(xs...))
    else
        return any(ismissing, xs) ? missing : f.f(xs...)
    end
end

"""
    passmissing(f)

Return a function that returns `missing` if any of its positional arguments
are `missing` (even if their number or type is not consistent with any of the
methods defined for `f`) and otherwise applies `f` to these arguments.

`passmissing` does not support passing keyword arguments to the `f` function.

# Examples
```jldoctest
julia> passmissing(sqrt)(4)
2.0

julia> passmissing(sqrt)(missing)
missing

julia> passmissing(sqrt).([missing, 4])
2-element Array{Union{Missing, Float64},1}:
  missing
   2.0

julia> passmissing((x,y)->"\$x \$y")(1, 2)
"1 2"

julia> passmissing((x,y)->"\$x \$y")(missing)
missing
"""
passmissing(f) = PassMissing{Core.Typeof(f)}(f)

end # module


text:
struct ValidateError <: Exception end

const rgroups = (:POINT, :ANALOG)

const descriptives = (:LABELS, :DESCRIPTIONS, :UNITS)

const rpoint = (:USED, :DATA_START, :FRAMES)
const ratescale = (:SCALE, :RATE)

const ranalog = (:USED, :GEN_SCALE, :OFFSET) ∪ ratescale
const bitsformat = (:BITS, :FORMAT)

const rforceplatf = (:TYPE, :ZERO, :CORNERS, :ORIGIN, :CHANNEL, :CAL_MATRIX)

const pointsigncheck = ((:POINT, :USED),
                   (:POINT, :DATA_START),
                   (:POINT, :FRAMES))

# const analogsigncheck = (:ANALOG, :USED)
# const fpsigncheck = (:FORCE_PLATFORM, :ZERO)

function validatec3d(header::Header, groups::Dict{Symbol,Group}; complete=false)
    # The following if-else ensures the minimum set of information needed to succesfully read a C3DFile
    if !(rgroups ⊆ keys(groups))
        if !haskey(groups, :ANALOG)
            groups[:ANALOG] = Group(0, Int8(6), false, Int8(0), "ANALOG", :ANALOG, Int16(0), UInt8(22), "Analog data parameters", Dict{Symbol,AbstractParameter}())
            groups[:ANALOG].params[:USED] = ScalarParameter(0, Int8(5), false, Int8(0), "USED", :USED, Int16(0), zero(Int16), UInt8(30), "Number of analog channels used")
        else
            d = setdiff(rgroups, keys(groups))
            msg = "Required group(s)"
            for p in d
                msg *= " :"*string(p)
            end
            msg *= " are missing"
            throw(ErrorException(msg))
        end
    end

    # Validate the :POINT group
    pointkeys = keys(groups[:POINT].params)
    if !(rpoint ⊆ pointkeys)
        # The minimum set of parameters in :POINT is rpoint
        d = setdiff(rpoint, pointkeys)
        msg = ":POINT is missing required parameter(s)"
        for p in d
            msg *= " :"*string(p)
        end
        throw(ErrorException(msg))
    end

    # Fix the sign for any point parameters that are likely to need it
    for (group, param) in pointsigncheck
        if any(signbit.(groups[group].params[param].data))
            groups[group].params[param] = unsigned(groups[group].params[param])
        end
    end

    if groups[:POINT].USED != 0 # There are markers
        if !(ratescale ⊆ pointkeys) # If there are markers, the additional set of required parameters is ratescale
            if !(:RATE ∈ pointkeys) && groups[:ANALOG].USED == 0
                # If there is no analog data, POINT:RATE isn't technically required
            else
                d = setdiff(rpoint, pointkeys)
                msg = ":POINT is missing required parameter(s)"
                for p in d
                    msg *= " :"*string(p)
                end
                throw(ErrorException(msg))
            end
        end

        if !(descriptives ⊆ pointkeys) # Check that the descriptive parameters exist
            if !haskey(groups[:POINT].params, :LABELS)
                # While the C3D file can technically be read in the absence of a LABELS parameter,
                # this implementation requires LABELS (for indexing)
                @debug ":POINT is missing parameter :LABELS"
                labels = [ "M"*string(i, pad=3) for i in 1:groups[:POINT].USED ]
                groups[:POINT].params[:LABELS] =
                      StringParameter(0, Int8(0), false, abs(groups[:POINT].gid), "LABELS", :LABELS, Int16(0), labels, UInt8(13), "Marker labels")
            elseif !haskey(groups[:POINT].params, :DESCRIPTIONS)
                @debug ":POINT is missing parameter :DESCRIPTIONS"
            elseif !haskey(groups[:POINT].params, :UNITS)
                @debug ":POINT is missing parameter :UNITS"
            end
        elseif groups[:POINT].params[:LABELS] isa ScalarParameter # ie There is only one used marker (or the others are unlabeled)
            groups[:POINT].params[:LABELS] = StringParameter(groups[:POINT].params[:LABELS])
        end

        # Valid labels are required for each marker by the C3DFile constructor
        if any(isempty.(groups[:POINT].LABELS)) ||
           length(groups[:POINT].LABELS) < groups[:POINT].USED # Some markers don't have labels
            i = 2
            while length(groups[:POINT].LABELS) < groups[:POINT].USED
                if haskey(groups[:POINT].params, Symbol("LABEL",i)) # Check for the existence of a runoff labels group
                    append!(groups[:POINT].LABELS, groups[:POINT].params[Symbol("LABEL",i)].data)
                    i += 1
                else
                    push!(groups[:POINT].LABELS, "")
                end
            end

            idx = findall(isempty, groups[:POINT].LABELS)
            labels = [ "M"*string(i, pad=3) for i in 1:length(idx) ]
            groups[:POINT].LABELS[idx] .= labels
        end

        if length(unique(groups[:POINT].LABELS)) !== length(groups[:POINT].LABELS)
            dups = String[]
            for i in 1:groups[:POINT].USED
                if !in(groups[:POINT].LABELS[i], dups)
                    push!(dups, groups[:POINT].LABELS[i])
                else
                    m = match(r"_(?<num>\d+)$", groups[:POINT].LABELS[i])

                    if m == nothing
                        groups[:POINT].LABELS[i] *= "_2"
                        push!(dups, groups[:POINT].LABELS[i])
                    else
                        newlabel = groups[:POINT].LABELS[i][1:(m.offset - 1)]*string('_',tryparse(Int,m[:num])+1)
                        groups[:POINT].LABELS[i] = newlabel
                        push!(dups, groups[:POINT].LABELS[i])
                    end
                end
            end
        end
    end # End validate :POINT

    # Validate the :ANALOG group
    analogkeys = keys(groups[:ANALOG].params)
    if !haskey(groups[:ANALOG].params, :USED)
        msg = ":ANALOG is missing required parameter :USED"
        throw(ErrorException(msg))
    end

    if signbit(groups[:ANALOG].USED)
        groups[:ANALOG].params[:USED] = unsigned(groups[:ANALOG].params[:USED])
    end

    if groups[:ANALOG].USED != 0 # There are analog channels

        @label analogkeychanged
        if !(ranalog ⊆ analogkeys) # If there are analog channels, the required set of parameters is ranalog
            if :OFFSETS ∈ analogkeys
                groups[:ANALOG].params[:OFFSET] = groups[:ANALOG].params[:OFFSETS]
                delete!(groups[:ANALOG].params, :OFFSETS)
                @goto analogkeychanged # OFFSETS might not be the only missing parameter
            else
                d = setdiff(ranalog, analogkeys)
                msg = ":ANALOG is missing required parameter(s)"
                for p in d
                    msg *= " :"*string(p)
                end
                throw(ErrorException(msg))
            end
        elseif !(descriptives ⊆ analogkeys) # Check that the descriptive parameters exist
            if !haskey(groups[:ANALOG].params, :LABELS)
                @debug ":ANALOG is missing parameter :LABELS"
                labels = [ "A"*string(i, pad=3) for i in 1:groups[:ANALOG].USED ]
                groups[:ANALOG].params[:LABELS] =
                      StringParameter(0, Int8(0), false, abs(groups[:ANALOG].gid), "LABELS", :LABELS, Int16(0), labels, UInt8(14), "Channel labels")
            elseif !haskey(groups[:ANALOG].params, :DESCRIPTIONS)
                @debug ":ANALOG is missing parameter :DESCRIPTIONS"
            elseif !haskey(groups[:ANALOG].params, :UNITS)
                @debug ":ANALOG is missing parameter :UNITS"
            end
        elseif groups[:ANALOG].params[:LABELS] isa ScalarParameter
            groups[:ANALOG].params[:LABELS] = StringParameter(groups[:ANALOG].params[:LABELS])
        end

        # Pad scale and offset if shorter than :USED
        l = length(groups[:ANALOG].SCALE)
        if l < groups[:ANALOG].USED
            append!(groups[:ANALOG].params[:SCALE].data, fill(Float32(1.0), groups[:ANALOG].USED - l))
        end

        l = length(groups[:ANALOG].OFFSET)
        if l < groups[:ANALOG].USED
            append!(groups[:ANALOG].params[:OFFSET].data, fill(Float32(1.0), groups[:ANALOG].USED - l))
        end

        if any(isempty.(groups[:ANALOG].LABELS)) ||
           length(groups[:ANALOG].LABELS) < groups[:ANALOG].USED # Some markers don't have labels
            i = 2
            while length(groups[:ANALOG].LABELS) < groups[:ANALOG].USED
                if haskey(groups[:ANALOG].params, Symbol("LABEL",i)) # Check for the existence of a runoff labels group
                    append!(groups[:ANALOG].LABELS, groups[:ANALOG].params[Symbol("LABEL",i)].data)
                    i += 1
                else
                    push!(groups[:ANALOG].LABELS, "")
                end
            end

            idx = findall(isempty, groups[:ANALOG].LABELS)
            labels = [ "A"*string(i, pad=3) for i in 1:length(idx) ]
            groups[:ANALOG].LABELS[idx] .= labels
        end

        if length(unique(groups[:ANALOG].LABELS)) !== length(groups[:ANALOG].LABELS)
            dups = String[]
            for i in 1:groups[:ANALOG].USED
                if !in(groups[:ANALOG].LABELS[i], dups)
                    push!(dups, groups[:ANALOG].LABELS[i])
                else
                    m = match(r"_(?<num>\d+)$", groups[:ANALOG].LABELS[i])

                    if m == nothing
                        groups[:ANALOG].LABELS[i] *= "_2"
                        push!(dups, groups[:ANALOG].LABELS[i])
                    else
                        newlabel = groups[:ANALOG].LABELS[i][1:(m.offset - 1)]*string('_',tryparse(Int,m[:num])+1)
                        groups[:ANALOG].LABELS[i] = newlabel
                        push!(dups, groups[:ANALOG].LABELS[i])
                    end
                end
            end
        end
    end # End if analog channels exist

    nothing
end


text:
# printing
using Printf, Statistics, JLD2
# Oceananigans
using Oceananigans
using Oceananigans.BoundaryConditions
using Oceananigans.Fields
using Oceananigans.OutputWriters
using Oceananigans.Diagnostics
using Oceananigans.Utils
using Oceananigans.AbstractOperations
using Oceananigans.Advection
# CUDA
using CUDA

include(pwd() * "/scripts/ic.jl")

max_stars_repo_name: sandreza/Mesoscale

max_issues_repo_name: sandreza/Mesoscale

max_forks_repo_name: sandreza/Mesoscale


text:
using Distributions:Random
using AdaptiveABC
using QuantEcon
using Plots
using Plots.PlotMeasures
using Random
using Distributions
using Optim
using NLopt

# Set Simulator Parameters and Seed
Random.seed!(37162)
ϕ = 0.7
θ = 0.3
σ = 0.01
parameter_true = [ϕ, θ, σ]
n = 1000
p0 = 100
parameter_names = ["phi", "theta", "sigma"]

# Create model, simulate to obtain observation vector and then calculate summary statistics
true_arma = ARMA(ϕ, θ, σ)
obs = simulation(true_arma, ts_length=n)
obs_summary_stats = get_summary_stats(obs)

prices = zeros(n + 1)
prices[1] = p0
for i in 1:n
    prices[i + 1] = exp(obs[i]) * prices[i]
end

# l = @layout [a b]

# p1 = plot(1:n, obs, legend=false,
#     xlab="Time", ylab="ARMA Log Returns"
# );

# p2 = plot(1:n + 1, prices, legend=false,
# xlab="Time", ylab="ARMA Mid-Price"
# );
# p3 = plot(p1, p2, layout=l, tickfontsize=6, guidefontsize=8,
#     titlefontsize=10, right_margin=5mm, size=(1000, 400)) ;

# savefig(p3, "examples/arma/arma_lr_mp_plot.pdf")

#### Calibration Assessment

#####
# Utility Functions
#####

# Define Summary function used in AdaptiveABC. Takes in a vector of parameters and outputs a success boolean and the
# simulated summart statistics
function summary_fn(parameters, n_summary_stats, n_replications)
    success = true
    model = ARMA(parameters...)
    sim_obs = Array{Float64,2}(undef, n, n_replications)
    summary_stats = Array{Float64,2}(undef, n_summary_stats, n_replications)
    for i in 1:n_replications
        sim_obs[:, i] = simulation(model, ts_length=n)
        summary_stats[:, i] = get_summary_stats(obs, sim_obs[:, i])
    end
    return success, summary_stats
end

prior = Prior([
    Uniform(0.05, 0.95),
    Uniform(0.05, 0.95),
    Uniform(0.005, 0.1),
])

#####
# Calibration Technique: ABC Rejection BBWM
#####

# Distance Function to measure distance between observed summary stats and simulated summary stats
boot_weight_matrix = BlockBootstrapWeightMatrix(1312, obs, get_summary_stats, 100, 10000)
weighted_bootstrap = WeightedBootstrap(obs, get_summary_stats, boot_weight_matrix)

abc_input_wb = ABCInput(
    prior,
    parameter_names,
    obs_summary_stats,
    summary_fn,
    17,
    weighted_bootstrap
)

abc_reject_wb = ABCRejection(
    abc_input_wb,
    100_000,
    10,
    250.0,
    parallel=true,
    seed=817221
)
posterior_means = AdaptiveABC.parameter_means(abc_reject_wb)
println()
println("ABC rejection BBWM: ϕ, θ, σ")
print(round.(posterior_means[:,end], digits=3))
println()

#####
# Calibration Technique: ABC Rejection MADWE
#####
abc_input_we = ABCInput(
    prior,
    parameter_names,
    obs_summary_stats,
    summary_fn,
    17,
    WeightedEuclidean(obs, get_summary_stats, ones(17), "MAD")
)

@time abc_reject_we = ABCRejection(
    abc_input_we,
    100_000,
    10,
    12.0,
    parallel=true,
    seed=2125179095
)
posterior_means = AdaptiveABC.parameter_means(abc_reject_we)
println()
println("ABC rejection MADWE: ϕ, θ, σ")
print(round.(posterior_means[:,end], digits=3))
println()
#####
# Calibration Technique: ABC-PMC BBWM
#####

# Distance Function to measure distance between observed summary stats and simulated summary stats
boot_weight_matrix = BlockBootstrapWeightMatrix(1312, obs, get_summary_stats, 100, 10000)
weighted_bootstrap = WeightedBootstrap(obs, get_summary_stats, boot_weight_matrix)

abc_input_wb = ABCInput(
    prior,
    parameter_names,
    obs_summary_stats,
    summary_fn,
    17,
    weighted_bootstrap
)

abc_pmc_out_wb = ABC_PMC(
    abc_input_wb,
    500,
    10,
    0.15,
    10000,
    9;
    parallel=true,
    batch_size=250,
    seed=2077971236
)
posterior_means = AdaptiveABC.parameter_means(abc_pmc_out_wb)
print(round.(posterior_means[:,end], digits=3))

#####
# Calibration Technique: ABC-PMC MADWE
#####
abc_input_we = ABCInput(
    prior,
    parameter_names,
    obs_summary_stats,
    summary_fn,
    17,
    WeightedEuclidean(obs, get_summary_stats, ones(17), "MAD")
)

abc_pmc_out_we = ABC_PMC(
    abc_input_we,
    100,
    10,
    0.15,
    100_000,
    9;
    parallel=true,
    batch_size=250,
    seed=589181349
)
posterior_means = AdaptiveABC.parameter_means(abc_pmc_out_we)
println("")
println("ABC-PMC MADWE: ϕ, θ, σ")
print(round.(posterior_means[:,end], digits=3))
println("")

#####
# Calibration Technique: Nelder-Mead BBWM
#####

# Distance Function to measure distance between observed summary stats and simulated summary stats
boot_weight_matrix = BlockBootstrapWeightMatrix(1312, obs, get_summary_stats, 100, 10000)
weighted_bootstrap = WeightedBootstrap(obs, get_summary_stats, boot_weight_matrix)

function f(x, weight, reps)
    success, sum_stats = summary_fn(x, 17, reps)
    return weight(sum_stats)
end

f_wb(x, grad) = f(x, weighted_bootstrap, 50)
lower = [0.05, 0.05, 0.005]
upper = [0.95, 0.95, 0.1]

nm_reps = 3
results = zeros(3, nm_reps)
Random.seed!(831838178)
for i in 1:nm_reps
    println("Iteration $i")
    opt = Opt(:LN_NELDERMEAD, 3)
    opt.lower_bounds = lower
    opt.upper_bounds = upper
    opt.xtol_rel = 1e-6
    opt.min_objective = f_wb
    init_params = rand(Int.(rand(MersenneTwister(), UInt32)), prior)

    (minf, minx, ret) = NLopt.optimize(opt, init_params)

    results[:,i] = minx
end
println()
println("Nelder-Mead BBWM: ")
round.(mean(results, dims=2), digits=3)
results
println()

#####
# Calibration Technique: Nelder-Mead MADWE
#####

weighted_mad = WeightedEuclidean(obs, get_summary_stats, ones(17), "MAD")
sim_sum_stats = zeros(17, 1000, 10)
Random.seed!(156271221)
for i in 1:1000
    params = rand(Int.(rand(MersenneTwister(), UInt32)), prior)
    success, sum_stats = summary_fn(params, 17, 10)
    sim_sum_stats[:,i,:] = sum_stats
end

inv_weights = [AdaptiveABC.MAD(sim_sum_stats[i, :, :]) for i in 1:17]  # For each statistic, calculate the MAD across all n_sim, n_rep values
weights = 1.0 ./ inv_weights

weighted_mad.weights = weights

f_mad(x, grad) = f(x, weighted_mad, 50)
lower = [0.05, 0.05, 0.005]
upper = [0.95, 0.95, 0.1]

nm_reps = 3
results = zeros(3, nm_reps)
Random.seed!(11019388531)
for i in 1:nm_reps
    println("Iteration $i")
    opt = Opt(:LN_NELDERMEAD, 3)
    opt.lower_bounds = lower
    opt.upper_bounds = upper
    opt.xtol_rel = 1e-6
    opt.min_objective = f_mad
    init_params = rand(Int.(rand(MersenneTwister(), UInt32)), prior)

    (minf, minx, ret) = NLopt.optimize(opt, init_params)

    results[:,i] = minx
end
println()
println("Nelder-Mead BBWM: ")
round.(mean(results, dims=2), digits=3)
results
println()

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

max_stars_repo_name: GantZA/CALOBMTMD

max_issues_repo_name: GantZA/CALOBMTMD

max_forks_repo_name: GantZA/CALOBMTMD


text:
using Cassette
using Documenter

makedocs(;
    modules = [Cassette],
    authors = "Jarrett Revels",
    sitename = "Cassette",
    format = Documenter.HTML(;
        prettyurls = get(ENV, "CI", "false") == "true",
        ansicolor = true,
    ),
    pages = Any[
        "Introduction" => "index.md",
        "Disclaimers" => "disclaimers.md",
        "Why Cassette?" => "whycassette.md",
        "The Overdubbing Mechanism" => "overdub.md",
        "Contextual Dispatch" => "contextualdispatch.md",
        "Contextual Compiler Pass Injection" => "contextualpass.md",
        "Contextual Tagging of Values" => "contextualtagging.md",
        "Cassette API Documentation" => "api.md",
        "Related Work" => "relatedwork.md",
    ],
    strict = false,
)

deploydocs(; repo = "github.com/jrevels/Cassette.jl",
             push_preview = true,
             )


text:
using ModelingToolkit
using Base.Test

@IVar t
@DVar x(t) y(t) z(t)

null_op = 0*t
@test simplify_constants(null_op) == Constant(0)

one_op = 1*t
@test simplify_constants(one_op) == t

identity_op = Operation(identity,[x])
@test simplify_constants(identity_op) == x

minus_op = -x
@test simplify_constants(minus_op) == -1*x
simplify_constants(minus_op)