text:
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the compressible Euler equations

equations = CompressibleEulerEquations1D(1.4)

"""
    initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations1D)

The Sedov blast wave setup based on Flash
- http://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel/node184.html#SECTION010114000000000000000
"""
function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations1D)
  # Set up polar coordinates
  inicenter = SVector(0.0)
  x_norm = x[1] - inicenter[1]
  r = abs(x_norm)

  # Setup based on http://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel/node184.html#SECTION010114000000000000000
  r0 = 0.21875 # = 3.5 * smallest dx (for domain length=4 and max-ref=6)
  # r0 = 0.5 # = more reasonable setup
  E = 1.0
  p0_inner = 6 * (equations.gamma - 1) * E / (3 * pi * r0)
  p0_outer = 1.0e-5 # = true Sedov setup

  # Calculate primitive variables
  rho = 1.0
  v1  = 0.0
  p   = r > r0 ? p0_outer : p0_inner

  return prim2cons(SVector(rho, v1, p), equations)
end

initial_condition = initial_condition_sedov_blast_wave

surface_flux = flux_lax_friedrichs
volume_flux  = flux_ranocha
basis = LobattoLegendreBasis(3)
shock_indicator_variable = density_pressure
indicator_sc = IndicatorHennemannGassner(equations, basis,
                                         alpha_max=1.0,
                                         alpha_min=0.001,
                                         alpha_smooth=true,
                                         variable=shock_indicator_variable)
volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
                                                 volume_flux_dg=volume_flux,
                                                 volume_flux_fv=surface_flux)
solver = DGSEM(basis, surface_flux, volume_integral)

cells_per_dimension = (64,)
coordinates_min = (-2.0,)
coordinates_max = ( 2.0,)
mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max, periodicity=true)

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)

###############################################################################
# ODE solvers, callbacks etc.

tspan = (0.0, 12.5)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

analysis_interval = 1000
analysis_callback = AnalysisCallback(semi, interval=analysis_interval)

alive_callback = AliveCallback(analysis_interval=analysis_interval)

save_solution = SaveSolutionCallback(interval=100,
                                     save_initial_solution=true,
                                     save_final_solution=true,
                                     solution_variables=cons2prim)

stepsize_callback = StepsizeCallback(cfl=0.5)

callbacks = CallbackSet(summary_callback,
                        analysis_callback,
                        alive_callback,
                        save_solution,
                        stepsize_callback)

###############################################################################
# run the simulation

sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt=stepsize_callback(ode), # solve needs some value here but it will be overwritten by the stepsize_callback
            save_everystep=false, callback=callbacks);
summary_callback() # print the timer summary

ext: jl


text:
function state2array(s::State)
    st = s
    res = State[]
    while !isnull(st.prev)
        unshift!(res, st)
        st = st.prev
    end
    unshift!(res, st)
    res
end

function stack2array(s::State)
    st = s
    res = Int[]
    while !isnull(st.left)
        unshift!(res, st.top)
        st = st.left
    end
    unshift!(res, st.top)
    res
end

# prints [ a/NN b/VB ][ c/NN d/PP ]
function Base.print(io::IO, s::State)
    stack = map(stack2array(s)) do i
        i == 0 ? "ROOT/ROOT" :
        getkey(s.parser.words, s.tokens[i].word) * "/" *
        getkey(s.parser.tags, s.tokens[i].tag)
    end
    stack = join(stack, " ")
    buffer = map(s.right:length(s.tokens)) do i
        getkey(s.parser.words, s.tokens[i].word) * "/" *
        getkey(s.parser.tags, s.tokens[i].tag)
    end
    buffer = join(buffer, " ")
    print(io, "[", stack, "][", buffer, "]")
end

function actstr(s::State{Labeled})
    act = s.prevact
    labels(id) = getkey(s.parser.labels, id)
    return acttype(act) == SHIFT ? "SHIFT" :
           acttype(act) == REDUCEL ? "REDUCE-lEFT($(labels((tolabel(act)))))" :
           acttype(act) == REDUCER ? "REDUCE-RIGHT($(labels((tolabel(act)))))" :
           throw("Invalid action: $(act).")
end

function actstr(s::State{Unlabeled})
    act = s.prevact
    return acttype(act) == SHIFT ? "SHIFT" :
           acttype(act) == REDUCEL ? "REDUCE-LEFT" :
           acttype(act) == REDUCER ? "REDUCE-RIGHT" :
           throw("Invalid action: $(act).")
end

function stacktrace{T}(io::IO, s::State{T})
    ss = state2array(s)
    for i in 2:length(ss)
        i > 1 && println(io, actstr(ss[i]))
        print(io, ss[i])
        println(io)
    end
end
stacktrace(s::State) = stacktrace(STDOUT, s)

function toconll(io::IO, s::State)
    pred = heads(s)
    for i in 1:length(s.tokens)
        t = s.tokens[i]
        items = [i, getkey(s.parser.words, t.word), "-",
                    getkey(s.parser.tags, t.tag), pred[i],
                    t.head, getkey(s.parser.labels, t.label)]
        println(io, join(items, "\t"))
    end
    println(io, "")
end
toconll(s::State) = toconll(STDOUT, s)

max_stars_repo_name: hshindo/JukaiNLP

max_issues_repo_name: hshindo/JukaiNLP

max_forks_repo_name: hshindo/JukaiNLP


text:
module EbitSolver
using Main.EbitODEMessages
using Main.EbitCloudSpatialExtend

export solve_ode, pack_ode_result

using DifferentialEquations
using ProtoBuf
using ParameterizedFunctions
using Memoize
using SpecialFunctions

_ret_codes_ = Dict(
    :Default => 0,
    :Success => 1,
    :MaxIters => 2,
    :DtLessThanMin => 3,
    :Unstable => 4,
    :InitialFailure => 5,
    :ConvergenceFailure => 6,
    :Failure => 7
)

function create_matrix(dim, sparse_values, f = (v, matrix) -> matrix[v.row,v.column] += v.value)
    matrix = zeros(dim, dim)
    map(v -> f(v,matrix), sparse_values)
    return matrix
end

struct EbitParameters
    qVₑ::Array{Float64,1}
    qVₜ::Array{Float64,1}
    A::Array{Float64,1}
    ϕ::Array{Float64,1}
    q::Array{Float64,1}

    source_n::Array{Float64,1}
    source_kT::Array{Float64,1}

    χ::Array{Float64,2}
    dN::Array{Float64,2}
    CX::Array{Float64,2}

    τ::Array{Float64,1}

    no_dimensions::UInt32
    min_N::Float64
    report_function
    ion_ion_overlap_function
    electron_beam_ion_overlap_function
    one_over_ion_cloud_vol_function
    heat_capacity_function

    function EbitParameters(dep::EbitODEMessages.DiffEqParameters, report_function,
                            ion_ion_overlap_function,
                            electron_beam_ion_overlap_function,
                            one_over_ion_cloud_vol_function,
                            heat_capacity_function)
        dim = dep.no_dimensions
        qVₑ = dep.qVe
        qVₜ = dep.qVt
        q = dep.q
        A = dep.mass_number
        ϕ = dep.spitzer_divided_by_overlap
        χ = create_matrix(dim, dep.inverted_collision_constant)
        dN = create_matrix(dim, dep.rate_of_change_divided_by_N)
        CX = create_matrix(dim, dep.dCharge_ex_divided_by_N_times_tau)

        return new(qVₑ, qVₜ, A, ϕ, q,
                   dep.source_terms_n, dep.source_terms_kt, χ, dN, CX,
                   Array{Float64}(undef, dep.no_dimensions),
                   dep.no_dimensions, dep.minimum_N,
                   report_function,
                   ion_ion_overlap_function,
                   electron_beam_ion_overlap_function,
                   one_over_ion_cloud_vol_function,
                   heat_capacity_function)
    end
end

function du(du::Array{Float64, 1}, u::Array{Float64,1}, p::EbitParameters, t::Float64)
    @inbounds begin
        N = view(u, 1:p.no_dimensions)
        Nτ = view(u, p.no_dimensions+1:p.no_dimensions*2)

        dN = view(du, 1:p.no_dimensions)
        dNτ = view(du, p.no_dimensions+1:p.no_dimensions*2)

        for i in 1:p.no_dimensions
            if N[i] > p.min_N
                p.τ[i] = Nτ[i] / (1.5 * N[i])
            else
                p.τ[i] = 0.0
            end
        end

        @simd for i in 1:p.no_dimensions
            ν = 0.0
            R_exchange_sum_j = 0.0
            dN[i] = p.source_n[i]
            dNτ[i] = p.source_kT[i]
            beta_i = p.q[i]/p.τ[i]
            @simd for j in 1:p.no_dimensions
                # calculate everything that is based on interaction of two species
                if (N[i] > p.min_N && N[j] > p.min_N && p.τ[j] > 0.0 && p.τ[i] > 0.0)
                    beta_j = p.q[j]/p.τ[j]
                    fᵢⱼ = p.ion_ion_overlap_function(beta_i, beta_j)
                    nⱼ = N[j] * p.one_over_ion_cloud_vol_function(beta_j)
                    arg = (p.τ[i]/p.A[i] + p.τ[j]/p.A[j])
                    Σ = p.χ[i,j] * nⱼ * arg^(-1.5)
                    ν += Σ
                    R_exchange_sum_j += fᵢⱼ * Σ * (p.τ[j] - p.τ[i])

                    dN[i] += p.CX[i,j]*N[j]*sqrt(p.τ[j])
                    dNτ[i] += p.CX[i,j]*Nτ[j]*sqrt(p.τ[j])
                end

                dN[i] += p.dN[i,j]*N[j]
                dNτ[i] += p.dN[i,j]*Nτ[j]
            end

            dNτ[i] += R_exchange_sum_j*N[i]

            if p.τ[i] > 0.0
                ɷ = p.qVₜ[i] / p.τ[i]
                dNi_dt = - N[i] * ν * ( exp(-ɷ) / ɷ - sqrt(ɷ)*(erf(ɷ) - 1))
                dNτ[i] += N[i] * p.electron_beam_ion_overlap_function(beta_i) * p.ϕ[i]
                dN[i] +=  dNi_dt
                dNτ[i] += - ( 2/3*N[i]*ν*exp(-ɷ) - dNi_dt)*p.τ[i]
            end

        end

        return du
    end
end

function create_initial_values(initial_values, dimensions)
    ret = zeros(2*dimensions)

    map((iv) ->
        begin
        ret[iv.index] = iv.number_of_particles
        ret[dimensions + iv.index] = iv.temperature_in_ev
        end,
        initial_values)

    return ret
end

@memoize function create_heat_capacity(r_stop, r_e, V_0)
    @info "Creating heat capacity approximation"
    EbitCloudSpatialExtend.create_interpolation(
        x -> EbitCloudSpatialExtend.heat_capacitance(x, r_stop, r_e, V_0))
end

@memoize function create_ion_ion_overlap_function(r_stop, r_e, V_0)
    @info "Creating ion/ion overlap approximation"
    EbitCloudSpatialExtend.create_2d_interpolation(
        (x,y) -> EbitCloudSpatialExtend.ion_ion_overlap(x, y, r_stop, r_e, V_0))
end

@memoize function create_electron_ion_overlap_function(r_stop, r_e, V_0)
    @info "Creating electron/ion overlap approximation"
    EbitCloudSpatialExtend.create_interpolation(
        x -> EbitCloudSpatialExtend.electron_ion_overlap(x, r_stop, r_e, V_0))
end

@memoize function create_one_over_ion_cloud_volume_function(r_stop, r_e, V_0, l_dt)
    @info "Creating ion density approximation"
    EbitCloudSpatialExtend.create_interpolation(
        x -> 1 / (pi * l_dt * (EbitCloudSpatialExtend.effective_radius(x, r_stop, r_e, V_0))^2))

end

@noinline function create_diffeq_prob(problem, report_function)
    if problem.problem_parameters.problem_type == EbitODEMessages.ProblemType.ODEProblem
        @debug "Creating ODEProblem"
        dp = problem.diff_eq_parameters
        p = EbitParameters(dp, report_function,
                           create_ion_ion_overlap_function(dp.r_dt, dp.r_e, dp.V_0),
                           create_electron_ion_overlap_function(dp.r_dt, dp.r_e, dp.V_0),
                           create_one_over_ion_cloud_volume_function(dp.r_dt, dp.r_e, dp.V_0, dp.l_dt),
                           create_heat_capacity(dp.r_dt, dp.r_e, dp.V_0))

        @debug "Created differential equation parameters" p.dN

        tspan = (problem.problem_parameters.time_span.start, problem.problem_parameters.time_span.stop)

        initial_values = create_initial_values(
            problem.diff_eq_parameters.initial_values,
            problem.diff_eq_parameters.no_dimensions
        )

        @debug "Created initial values from list" initial_values

        global last_p = p
        global last_initial_values = initial_values

        return ODEProblem(du, initial_values, tspan, p)
    else
        throw(ErrorException("Unknown Problem type"))
    end
end

function pack_values(solution, nuclides)
    s = solution
    len_t = length(s.t)
    len_i = length(nuclides)

    ret_n = [EbitODEMessages.ValuesPerNuclide(nuclide=nuclide, values=Array{Float64}(undef, len_t))
             for nuclide in nuclides]
    ret_kT = [EbitODEMessages.ValuesPerNuclide(nuclide=nuclide, values=Array{Float64}(undef, len_t))
              for nuclide in nuclides]

    for t in 1:len_t
        for i in 1:len_i
            ret_n[i].values[t] = s.u[t][i]
            ret_kT[i].values[t] = s.u[t][len_i + i]
        end
    end

    return ret_n, ret_kT
end

function pack_ode_result(solution, problem, start, stop)
    ns, kts = pack_values(solution, problem.nuclides)
    return EbitODEMessages.Result(problem=problem, start_time=start, stop_time=stop,
                                  return_code=get(_ret_codes_, solution.retcode, 0),
                                  times=solution.t, n=ns, kT=kts)
end

function pack_ode_msg(res)
    return EbitODEMessages.Message(msg_type=EbitODEMessages.MessageType.ODEResult, ode_result = res)
end

function callback(u,t,integrator)
    if integrator.p.report_function(t)
        terminate!(integrator)
        @info "Terminating integration" t
    end
end

function solve_ode(problem, report_progress)
    start_ = time()
    @info "Starting solver"
    sol = solve(create_diffeq_prob(problem, report_progress),
                # alghints=[:stiff],
                saveat=problem.solver_parameters.saveat,
                # abstol=1e-3,
                callback=FunctionCallingCallback(callback, funcat=problem.solver_parameters.saveat),
                )
    stop = time()
    @info "Finished integration" start_ stop
    return pack_ode_msg(pack_ode_result(sol, problem, start_, stop))
end

end


text:
str = """Hello,
        world.
         """

print(str)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module TestKmerAnalysis

using KmerAnalysis, BioSequences, Test

@testset "Count Mode" begin
    iterres = BioSequences.MerIterResult(1, mer"TCGGAAGACTAAGGAGCCTATCAGATATGTC", mer"GACATATCTGATAGGCTCCTTAGTCTTCCGA")
    @test CANONICAL(iterres) === mer"GACATATCTGATAGGCTCCTTAGTCTTCCGA"
    @test NONCANONICAL(iterres) === mer"TCGGAAGACTAAGGAGCCTATCAGATATGTC"
end

include("MerCount.jl")
include("Counters.jl")

end # module


text:
@testset "send_window.jl" begin
    @testset "Initial receiver advertised window is respected" begin
        isn = rand(UInt32(0):typemax(UInt32))

        sender = TCPSender(; fixed_isn=WrappingInt32(isn))
        fill_window!(sender)

        expect_seg(sender; no_flag=true, syn=true, payload_size=0, seqno=WrappingInt32(isn))
        ack_received_test(sender, WrappingInt32(isn + 1), 4)
        expect_no_seg(sender)
        write_bytes!(sender, "abcdefg")
        expect_seg(sender; no_flag=true, data = "abcd")
        expect_no_seg(sender)
    end

    @testset "Immediate window is respected" begin
        isn = rand(UInt32(0):typemax(UInt32))

        sender = TCPSender(; fixed_isn=WrappingInt32(isn))
        fill_window!(sender)

        expect_seg(sender; no_flag=true, syn=true, payload_size=0, seqno=WrappingInt32(isn))
        ack_received_test(sender, WrappingInt32(isn + 1), 6)
        expect_no_seg(sender)
        write_bytes!(sender, "abcdefg")
        expect_seg(sender; no_flag=true, data = "abcdef")
        expect_no_seg(sender)
    end
    @testset "Window" begin
        N_REPS = 1000
        MIN_WIN = 5
        MAX_WIN = 100
        for i in 1:N_REPS

            isn = rand(UInt32(0):typemax(UInt32))
            len = rand(MIN_WIN:MAX_WIN)

            sender = TCPSender(; fixed_isn=WrappingInt32(isn))
            fill_window!(sender)
            expect_seg(sender; no_flag=true, syn=true, payload_size=0, seqno=WrappingInt32(isn))
            ack_received_test(sender, WrappingInt32(isn + 1), len)
            expect_no_seg(sender)
            write_bytes!(sender, "a"^(2N_REPS))
            expect_seg(sender; payload_size=len, no_flag=true)
            expect_no_seg(sender)
        end
    end
    @testset "Window growth is exploited" begin
        isn = rand(UInt32(0):typemax(UInt32))

        sender = TCPSender(; fixed_isn=WrappingInt32(isn))
        fill_window!(sender)
        expect_seg(sender; no_flag=true, syn=true, payload_size=0, seqno=WrappingInt32(isn))
        ack_received_test(sender, WrappingInt32(isn + 1), 4)
        expect_no_seg(sender)
        write_bytes!(sender, "0123456789")
        expect_seg(sender; no_flag=true, data="0123")
        ack_received_test(sender, WrappingInt32(isn+5), 5)
        expect_seg(sender; no_flag=true, data="45678")
        expect_no_seg(sender)
    end
    @testset "FIN flag occupies space in window" begin
        isn = rand(UInt32(0):typemax(UInt32))

        sender = TCPSender(; fixed_isn=WrappingInt32(isn))
        fill_window!(sender)
        expect_seg(sender; no_flag=true, syn=true, payload_size=0, seqno=WrappingInt32(isn))
        ack_received_test(sender, WrappingInt32(isn + 1), 7)
        expect_no_seg(sender)
        write_bytes!(sender, "1234567")
        sender_close(sender)
        expect_seg(sender; no_flag=true, data="1234567")
        expect_no_seg(sender)

        ack_received_test(sender, WrappingInt32(isn+8), 1)
        expect_seg(sender; fin=true, data="")
        expect_no_seg(sender)
    end
    @testset "FIN flag occupies space in window (part II)" begin
        isn = rand(UInt32(0):typemax(UInt32))

        sender = TCPSender(; fixed_isn=WrappingInt32(isn))
        fill_window!(sender)
        expect_seg(sender; no_flag=true, syn=true, payload_size=0, seqno=WrappingInt32(isn))
        ack_received_test(sender, WrappingInt32(isn + 1), 7)
        expect_no_seg(sender)
        write_bytes!(sender, "1234567")
        sender_close(sender)
        expect_seg(sender; no_flag=true, data="1234567")
        expect_no_seg(sender)

        ack_received_test(sender, WrappingInt32(isn+1), 8)
        expect_seg(sender; fin=true, data="")
        expect_no_seg(sender)
    end
    @testset "Piggyback FIN in segment when space is available" begin
        isn = rand(UInt32(0):typemax(UInt32))

        sender = TCPSender(; fixed_isn=WrappingInt32(isn))
        fill_window!(sender)
        expect_seg(sender; no_flag=true, syn=true, payload_size=0, seqno=WrappingInt32(isn))
        ack_received_test(sender, WrappingInt32(isn + 1), 3)
        expect_no_seg(sender)
        write_bytes!(sender, "1234567")
        sender_close(sender)
        expect_seg(sender; no_flag=true, data="123")
        expect_no_seg(sender)

        ack_received_test(sender, WrappingInt32(isn+1), 8)
        expect_seg(sender; fin=true, data="4567")
        expect_no_seg(sender)
    end

end

max_stars_repo_name: AquaIndigo/JLSponge

max_issues_repo_name: AquaIndigo/JLSponge

max_forks_repo_name: AquaIndigo/JLSponge


text:
module DataTests

using Photon, Test, JLD2

resetContext()
ctx = getContext()

getimages(s=224) = KorA(randn(ctx.dtype,s,s,3,4))

function test_dataset()
    ds = TestDataset((28,28,1), (10,), 50)
	@test length(ds) == 50
	sample = ds[1]
	@test  size(sample[1]) == (28,28,1)
	@test  size(sample[2]) == (10,)
end

function test_arraydataset()
	X = [randn(100,100,3) for _ in 1:10]
	Y = [rand(0:9) for _ in 1:10]

    ds = ArrayDataset(X, Y)
	@test length(ds) == 10
	sample = ds[1]
	@test  size(sample[1]) == (100,100,3)
end

function test_jld()
	rm("example.jld2", force=true)

	jldopen("example.jld2", "w") do file
	    file["image4"] = (randn(Float32,28,28,1), rand(0:9,1))
	    file["image5"] = (randn(Float32,28,28,1), rand(0:9,1))
	end

	f = jldopen("example.jld2", "r")
	ds = JLDDataset(f)
	(X,Y) = ds[1]
	close(f)
	@assert size(X) == (28,28,1)
	@assert size(Y) == (1,)
	rm("example.jld2", force=true)
end

function test_output()

	samples = 9
	X1 = [zeros(Float32,28,28,1) for _ = 1:samples]
	Y1 = [ones(2) for _ = 1:samples]
	data = ArrayDataset(X1,Y1) |> MiniBatch(4)

	@assert data.batchsize == 4
	@assert length(data) == samples ÷ 4

	for (X,Y) in data
		@assert size(X,4) == data.batchsize
		@assert size(Y,2) == data.batchsize
		@assert sum(X) == 0.0
		@assert sum(Y) == 1.0 * 2 * data.batchsize
	end

end

function test_minibatch()
    data = TestDataset((28,28,1), (10,), 50) |> MiniBatch(8)
	sample = first(data)
	@test  size(sample[1]) == (28,28,1,8)
	@test  size(sample[2]) == (10,8)

	data = TestDataset((28,28,1), (10,), 50) |> MiniBatch(8, shuffle=false)
	sample = first(data)
end

function test_training()
    data = TestDataset((28,28,1), (10,), 100) |> MiniBatch()

	model = Sequential(
        Conv2D(16, 3, relu),
        Conv2D(32, 3, relu),
        MaxPool2D(),
        Dense(32, relu),
        Dense(10)
    )

	workout = Workout(model, MSELoss())
	fit!(workout, data, epochs=2)
	@assert workout.steps == (100÷8)*2

end

function test_threading(sleep)
    data = TestDataset((28,28,1), (10,), 100; sleep=sleep) |> MiniBatch(32)

	model = Sequential(
        Conv2D(16, 3, relu),
        Conv2D(32, 3, relu),
        MaxPool2D(),
        Dense(32, relu),
        Dense(10)
    )

	workout = Workout(model, MSELoss())
	@time fit!(workout, data, epochs=2)
	@assert workout.steps == (100÷32)*2
end

function test_transformers()
	ds = TestDataset((10,),(1,),100)
	ds = ds |> NoisingTransfomer()
	sample = ds[1]
	@assert size(sample[1]) == (10,)

	ds = TestDataset((100,100,3),(20,20,1),100)
	ds = ds |> ImageCrop((50,50),(10,10))
	sample = ds[10]
	@assert size(sample[1]) == (50,50,3)
	@assert size(sample[2]) == (10,10,1)

	X = [randn(100,100,3) for _ in 1:10]
	Y = [rand(0:9) for _ in 1:10]

	ds = ArrayDataset(X,Y)
	ds = ds |> OneHotEncoder(0:9)
	sample = ds[5]
	@assert size(sample[2]) == (10,)
end

function test_split()
	ds = TestDataset((30,30,3),(10,),100)

	t,v = ds |> Split()
	@assert length(t) == 80
	@assert length(v) == 20

	X,Y = t[4]
	@assert size(X) == (30,30,3)
	@assert size(Y) == (10,)

	X,Y = v[10]
	@assert size(X) == (30,30,3)
	@assert size(Y) == (10,)
end

@testset "Data" begin
    test_dataset()
	test_arraydataset()
	test_jld()
	test_output()
	test_minibatch()
	test_training()
	test_split()
	test_threading(0.01)
	test_threading(0.02)
	test_transformers()
end

end


text:
using Test
using Random
using CellLists

function brute_force(p::Array{T, 2}, r::T) where T <: AbstractFloat
    ps = Vector{Tuple{Int, Int}}()
    n, d = size(p)
    for i in 1:(n-1)
        for j in (i+1):n
            if distance_condition(p[i, :], p[j, :], r)
                push!(ps, (i, j))
            end
        end
    end
    return ps
end

function cell_lists(p::Array{T, 2}, r::T) where T <: AbstractFloat
    c = CellList(p, r)
    return near_neighbors(c, p, r)
end

function test_cell_lists_serial(rng::AbstractRNG, ns::Vector{Int}, ds::Vector{Int}, rs::Vector{<:AbstractFloat}, iterations::Int)
    for (n, d, r) in Iterators.product(ns, ds, rs)
        @info "Testing cell lists serial: n: $n | d: $d | r: $r"
        for i in 1:iterations
            p = 2 .* rand(rng, n, d) .- 1.0
            a = Set(Set.(brute_force(p, r)))
            b = Set(Set.(cell_lists(p, r)))
            @test b == a
        end
    end
end

function test_merge(rng::AbstractRNG, ns, ds, rs)
    for (n, d, r) in Iterators.product(ns, ds, rs)
        @info "Testing merge: n: $n | d: $d | r: $r"
        p = 2 .* rand(rng, n, d) .- 1.0
        l = fld(n, 2)
        c1 = CellList(p, r)
        c2 = CellList(p[1:l, :], r)
        c3 = CellList(p[(l+1):end, :], r; offset=l)
        c4 = merge(c1, c2)
        @test Set(keys(c1.data)) == Set(keys(c4.data))
        for cell in keys(c1.data)
            @test Set(c1.data[cell]) == Set(c4.data[cell])
        end
    end
end

function test_cell_list_constructor_threads(rng::AbstractRNG, ns, ds, rs)
    for (n, d, r) in Iterators.product(ns, ds, rs)
        @info "Testing cell list constructor threads: n: $n | d: $d | r: $r"
        p = 2 .* rand(rng, n, d) .- 1.0
        c1 = CellList(p, r)
        c2 = CellList(p, r, Val(:threads))
        @test Set(keys(c1.data)) == Set(keys(c2.data))
        for cell in keys(c1.data)
            @test Set(c1.data[cell]) == Set(c2.data[cell])
        end
    end
end

function test_near_neighbors_threads(rng::AbstractRNG, ns::Vector{Int}, ds::Vector{Int}, rs::Vector{<:AbstractFloat}, iterations::Int)
    for (n, d, r) in Iterators.product(ns, ds, rs)
        @info "Testing near neighbors threads: n: $n | d: $d | r: $r"
        for i in 1:iterations
            p = 2 .* rand(rng, n, d) .- 1.0
            c = CellList(p, r)
            a = near_neighbors(c, p, r)
            b = near_neighbors(c, p, r, Val(:threads))
            @test Set(Set.(b)) == Set(Set.(a))
        end
    end
end

function test_near_neighbors_threads_large(rng::AbstractRNG, n::Int, d::Int, r::Float64)
    @info "Testing near neighbors threads: n: $n | d: $d | r: $r"
    p = 2 .* rand(rng, n, d) .- 1.0
    c = CellList(p, r)
    a = near_neighbors(c, p, r)
    b = near_neighbors(c, p, r, Val(:threads))
    @test Set(Set.(b)) == Set(Set.(a))
end

const rng = MersenneTwister(894)
const ns = [0, 1, 2, 10, 100]
const ds = [1, 2, 3]
const rs = [0.1, 0.2, 0.3, 0.5, 1.0, 2.0]
test_cell_lists_serial(rng, ns, ds, rs, 20)
test_merge(rng, ns, ds, rs)
test_cell_list_constructor_threads(rng, ns, ds, rs)
test_near_neighbors_threads(rng, ns, ds, rs, 20)
test_near_neighbors_threads_large(rng, 20000, 2, 0.01)
test_near_neighbors_threads_large(rng, 30000, 3, 0.01)


text:
using Test
include("./computer.jl")
input_data = open("input09.txt") do file
    data = readline(file)
    entries = split(data, ",")
    numbers = map(x->parse(Int64, x), entries)
    numbers
end

@test run_program([104,1125899906842624,99]).output == [1125899906842624]

max_stars_repo_name: buotex/aoc

max_issues_repo_name: buotex/aoc

max_forks_repo_name: buotex/aoc


text:
using BlockArrays, LinearAlgebra, Test

include("test_blockindices.jl")
include("test_blockarrays.jl")
include("test_blockviews.jl")
include("test_blockrange.jl")
include("test_blockarrayinterface.jl")

include("test_blockbroadcast.jl")

include("../docs/make.jl")


text:
module Advection

export
    div_𝐯u, div_𝐯v, div_𝐯w, div_Uc,

    momentum_flux_uu,
    momentum_flux_uv,
    momentum_flux_uw,
    momentum_flux_vu,
    momentum_flux_vv,
    momentum_flux_vw,
    momentum_flux_wu,
    momentum_flux_wv,
    momentum_flux_ww,
    advective_tracer_flux_x,
    advective_tracer_flux_y,
    advective_tracer_flux_z,

    CenteredSecondOrder,
    UpwindBiasedFirstOrder,
    UpwindBiasedThirdOrder,
    UpwindBiasedFifthOrder,
    CenteredFourthOrder,
    WENO5,
    VectorInvariant,
    EnergyConservingScheme,
    EnstrophyConservingScheme

using DocStringExtensions

using Oceananigans.Grids
using Oceananigans.Operators

import Oceananigans.Grids: required_halo_size

abstract type AbstractAdvectionScheme{Buffer} end
abstract type AbstractCenteredAdvectionScheme{Buffer} <: AbstractAdvectionScheme{Buffer} end
abstract type AbstractUpwindBiasedAdvectionScheme{Buffer} <: AbstractAdvectionScheme{Buffer} end

required_halo_size(scheme::AbstractAdvectionScheme{Buffer}) where Buffer = Buffer + 1

include("centered_advective_fluxes.jl")
include("upwind_biased_advective_fluxes.jl")
include("flat_advective_fluxes.jl")

include("upwind_biased_first_order.jl")
include("centered_second_order.jl")
include("upwind_biased_third_order.jl")
include("centered_fourth_order.jl")
include("upwind_biased_fifth_order.jl")
include("weno_fifth_order.jl")
include("vector_invariant_advection.jl")

include("topologically_conditional_interpolation.jl")

include("momentum_advection_operators.jl")
include("tracer_advection_operators.jl")

end # module


text:
immutable JustStrategy{T} <: AbstractStrategy
    x::T
end

just(x) = JustStrategy(x)

Base.minimum(x::JustStrategy) = x.x
arbitrary(strategy::JustStrategy) = strategy.x
shrink{T}(x::T, strategy::JustStrategy{T}) = Task(() -> nothing)

max_stars_repo_name: cpritcha/PropertyCheck

max_issues_repo_name: cpritcha/PropertyCheck

max_forks_repo_name: cpritcha/PropertyCheck


text:
connectivity(network::WattsStrogatzNetwork) = network.k

"""
	numshortcuts(network::WattsStrogatzNetwork)

Return the number of connections that were rewired.
"""
numshortcuts(network::WattsStrogatzNetwork) = network._numShortcuts

#---------- Calculators ----------
calcNumConnections!(network::WattsStrogatzNetwork) =
	network._props.numConnections = Int(network.k * network.N / (isdirected(network) ? 1 : 2))

function generateWSAdjMat(
	N::Integer,
	k::Integer,
	β::Real;
	directed::Bool,
	seed::Integer
)
	mat = adjMat(RegularNetwork(N, k), sparse=true)
	(β <= 0) && (return mat, 0)

	numShortcuts = 0
	rng = Random.MersenneTwister(seed)
	edges = allEdges(mat, directed)

	for edge in edges
		if rand(rng) < β
			# Find a valid reconnection
		    i = rand(rng, 1:N); j = rand(rng, 1:N)
			while mat[i,j] || i == j
				i = rand(rng, 1:N)
		    	j = rand(rng, 1:N)
			end

			mat[i,j] = 1
			mat[edge[1], edge[2]] = 0

			if !directed
				mat[j,i] = 1
				mat[edge[2], edge[1]] = 0
			end

			numShortcuts += 1
		end
	end

	return SparseArrays.dropzeros(mat), numShortcuts
end

function generateWSAdjMat(
	N::Integer,
	k::Integer,
	numShortcuts::Integer;
	directed::Bool,
	seed::Integer
)
	mat = adjMat(RegularNetwork(N, k), sparse=true)
	(numShortcuts <= 0) && (return mat)

	rng = Random.MersenneTwister(seed)
	edges = allEdges(mat, directed)

	# Connections that will be rewired
	rewiredEdges = Array{Tuple{Int,Int}}(undef, numShortcuts)
	for i in 1:numShortcuts
		rewiredEdges[i] = rand(rng, edges)
		edges = filter(x -> x != rewiredEdges[i], edges)
	end

	for edge in rewiredEdges
		# Find a valid reconnection
	    i = rand(rng, 1:N); j = rand(rng, 1:N)
		while mat[i,j] || i == j
			i = rand(rng, 1:N)
	    	j = rand(rng, 1:N)
		end

		mat[i,j] = 1
		mat[edge[1], edge[2]] = 0

		if !directed
			mat[j,i] = 1
			mat[edge[2], edge[1]] = 0
		end
	end

	return SparseArrays.dropzeros(mat)
end


text:
using JuliaGSGP
using CSV, DataFrames
using ScikitLearn: CrossValidation.train_test_split

# df = CSV.read("path/to/dataset", DataFrame)
df = df[2:end]
data = convert(Matrix, df)
X = data[:, 1:end-1]
y = data[:, end]
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=1)
X_train = X_train' # columns corresponds to one instances
X_test = X_test'

aq(a, b) = a / sqrt(1+b^2)

pset = PrimitiveSet()
add_function!(pset, +, 2) #(pset, function, arity)
add_function!(pset, -, 2)
add_function!(pset, *, 2)
add_function!(pset, aq, 2)

for i in 1:size(X)[2]
    add_variable!(pset, i)
end

cfg = Config(population_size=100, num_generations=100, p_crossover=0.5, p_mutation=0.5, mutation_step=0.01, min_depth=1, max_depth=6, tournament_size=4, pset=pset, create_rand_tree=generate_grow, create_population=ramped_half_and_half, fitness=rmse)

best_idx, rec = gsgp(cfg, X_train, y_train, X_test, y_test);

# `reconstruct_tree()` will take long time to execute.
# tree = reconstruct_tree(best_idx, cfg.mutation_step)

max_stars_repo_name: kogad/JuliaGSGP

max_issues_repo_name: kogad/JuliaGSGP

max_forks_repo_name: kogad/JuliaGSGP


text:
"Explicit Runge-Kutta integrator."
immutable IntegratorERK{T,F} <: Integrator{T}
    equation::ODE{T,F}
    tableau::TableauERK{T}
    Δt::T

    x::Array{T,1}
    y::Array{T,1}
    X::Array{T,2}
    Y::Array{T,2}
    F::Array{T,2}
    tX::Array{T,1}
    tF::Array{T,1}

    function IntegratorERK(equation, tableau, Δt)
        D = equation.d
        S = tableau.s
        new(equation, tableau, Δt, zeros(T,D), zeros(T,D), zeros(T,D,S), zeros(T,D,S), zeros(T,D,S), zeros(T,D), zeros(T,D))
    end
end

function IntegratorERK{T,F}(equation::ODE{T,F}, tableau::TableauERK{T}, Δt::T)
    IntegratorERK{T,F}(equation, tableau, Δt)
end

"Integrate ODE with explicit Runge-Kutta integrator."
function integrate!{T,F}(int::IntegratorERK{T,F}, sol::SolutionODE{T})
    local tᵢ::T

    # loop over initial conditions
    for m in 1:sol.n0
        # copy initial conditions from solution
        simd_copy_xy_first!(int.x, sol, 0, m)

        # loop over time steps
        for n in 1:sol.ntime
            # compute internal stages
            fill!(int.Y, zero(T))
            for i in 1:int.tableau.s
                for k in 1:sol.nd
                    for j = 1:i-1
                        int.Y[k,i] += int.tableau.a[i,j] * int.F[k,j]
                    end
                    int.X[k,i] = int.x[k] + int.Δt * int.Y[k,i]
                end
                tᵢ = sol.t[n] + int.Δt * int.tableau.c[i]
                simd_copy_xy_first!(int.tX, int.X, i)
                int.equation.f(tᵢ, int.tX, int.tF)
                simd_copy_yx_first!(int.tF, int.F, i)
            end

            # compute final update
            simd_mult!(int.y, int.F, int.tableau.b)
            simd_axpy!(int.Δt, int.y, int.x)

            # copy to solution
            if mod(n, sol.nsave) == 0
                simd_copy_yx_first!(int.x, sol, div(n, sol.nsave), m)
            end
        end
    end
    nothing
end


text:
# Recombinations
# ==============
function average{T <: Vector}(population::Vector{T})
    obj = zeros(eltype(T), length(population[1]))
    l = length(population)
    for i in 1:l
        obj += population[i]
    end
    return obj./l
end

function marriage{T <: Vector}(population::Vector{T})
    s = length(population)
    l = length(population[1])
    obj = zeros(eltype(T), l)
    for i in 1:l
        obj[i] = population[rand(1:s)][i]
    end
    return obj
end

# Strategy recombinations
# =======================
function averageSigma1{S <: Strategy}(ss::Vector{S})
    s = copy(ss[1])
    σ = 0.0
    l = length(ss)
    for i in 1:l
        σ += ss[i][:σ]
    end
    s[:σ] = σ/l
    return s
end

function averageSigmaN{S <: Strategy}(ss::Vector{S})
    s = copy(ss[1])
    σ = zeros(length(ss[1][:σ]))
    l = length(ss)
    for i in 1:l
        σ += ss[i][:σ]
    end
    s[:σ] = σ./l
    return s
end

# ==================
# Genetic algorithms
# ==================

# Binary crossovers
# -----------------

# Single point crossover
function singlepoint{T <: Vector}(v1::T, v2::T)
    l = length(v1)
    c1 = copy(v1)
    c2 = copy(v2)
    pos = rand(1:l)
    for i in pos:l
        vswap!(c1, c2, i)
    end
    return c1, c2
end

# Two point crossover
function twopoint{T <: Vector}(v1::T, v2::T)
    l = length(v1)
    c1 = copy(v1)
    c2 = copy(v2)
    from, to = rand(1:l, 2)
    from, to = from > to ? (to, from)  : (from, to)
    for i in from:to
        vswap!(c1, c2, i)
    end
    return c1, c2
end

# Uniform crossover
function uniform{T <: Vector}(v1::T, v2::T)
    l = length(v1)
    c1 = copy(v1)
    c2 = copy(v2)
    xch = randbool(l)
    for i in 1:l
        if xch[i]
            vswap!(c1, c2, i)
        end
    end
    return c1, c2
end

# Real valued crossovers
# ----------------------

function discrete{T <: Vector}(v1::T, v2::T)
    l = length(v1)
    c1 = similar(v1)
    c2 = similar(v2)
    sltc = randbool(2,l)
    for i in 1:l
        c1[i] = sltc[1,i] ? v1[i] : v2[i]
        c2[i] = sltc[2,i] ? v2[i] : v1[i]
    end
    return c1, c2
end

# Weighted arithmetic mean
function waverage(w::Vector{Float64})
    function wavexvr{T <: Vector}(v1::T, v2::T)
        c1 = (v1+v2)./w
        return c1, copy(c1)
    end
    return wavexvr
end

# Intermediate recombination
function intermediate(d::Float64 = 0.0)
    function intermxvr{T <: Vector}(v1::T, v2::T)
        l = length(v1)
        α = (1.0+2d)*rand(l)-d
        c1 = v2 .+ α .* (v1 - v2)
        α = (1.0+2d)*rand(l)-d
        c2 = v1 .+ α .* (v2 - v1)
        return c1, c2
    end
    return intermxvr
end

# Line recombination
function line(d::Float64 = 0.0)
    function linexvr{T <: Vector}(v1::T, v2::T)
        α1, α2 = (1.0+2d)*rand(2)-d
        c1 = v2 .+ α2 * (v1 - v2)
        c2 = v1 .+ α1 * (v2 - v1)
        return c1, c2
    end
    return linexvr
end

# Permutation crossovers
# ----------------------

# Partially mapped crossover
function pmx{T <: Vector}(v1::T, v2::T)
    s = length(v1)
    from, to = rand(1:s, 2)
    from, to = from > to ? (to, from)  : (from, to)
    c1 = similar(v1)
    c2 = similar(v2)

    # Swap
    c1[from:to] = v2[from:to]
    c2[from:to] = v1[from:to]

    # Fill in from parents
    for i in vcat(1:from-1, to+1:s)
        # Check conflicting offspring
        in1 = inmap(v1[i], c1, from, to)
        if in1 == 0
            c1[i] = v1[i]
        else
            tmpin = in1
            while tmpin > 0
                tmpin = inmap(c2[in1], c1, from, to)
                in1 = tmpin > 0 ? tmpin : in1
            end
            c1[i] = v1[in1]
        end

        in2 = inmap(v2[i], c2, from, to)
        if in2 == 0
            c2[i] = v2[i]
        else
            tmpin = in2
            while tmpin > 0
                tmpin = inmap(c1[in2], c2, from, to)
                in2 = tmpin > 0 ? tmpin : in2
            end
            c2[i] = v2[in2]
        end
    end
    return c1, c2
end

# Order crossover
function ox1{T <: Vector}(v1::T, v2::T)
    s = length(v1)
    from, to = rand(1:s, 2)
    from, to = from > to ? (to, from)  : (from, to)
    c1 = similar(v1)
    c2 = similar(v2)
    # Swap
    c1[from:to] = v2[from:to]
    c2[from:to] = v1[from:to]
    # Fill in from parents
    return c1, c2
end

# Cycle crossover
function cx{T <: Vector}(v1::T, v2::T)
    # TODO
end

# Order-based crossover
function ox2{T <: Vector}(v1::T, v2::T)
    # TODO
end

# Position-based crossover
function pos{T <: Vector}(v1::T, v2::T)
    # TODO
end

# Utils
# =====
function vswap!{T <: Vector}(v1::T, v2::T, idx::Int)
    val = v1[idx]
    v1[idx] = v2[idx]
    v2[idx] = val
end

function inmap{T}(v::T, c::Vector{T}, from::Int, to::Int)
    exists = 0
    for j in from:to
        if exists == 0 && v == c[j]
            exists = j
        end
    end
    return exists
end


text:
import BenchmarkTools
import ConcurrentCollectionsBenchmarks
import JSON

function git_info(dir = @__DIR__)
    git(cmd) = strip(read(setenv(`git $cmd`; dir), String))
    return (;
        revision = git(`rev-parse HEAD`),
        status = git(`status --short --untracked-files=no --porcelain`),
    )
end

function julia_info()
    return (
        version = string(VERSION),
        git = (
            commit = Base.GIT_VERSION_INFO.commit,
            branch = Base.GIT_VERSION_INFO.branch,
        ),
        is_debugbuild = ccall(:jl_is_debugbuild, Cint, ()) != 0,
        libllvm_version = string(Base.libllvm_version),
        Sys = (
            WORD_SIZE = Sys.WORD_SIZE,
            JIT = Sys.JIT,
            # CPU_NAME = Sys.CPU_NAME,
            # CPU_THREADS = Sys.CPU_THREADS,
        ),
        env = Dict(k => v for (k, v) in ENV if startswith(k, "JULIA_")),
    )
end

function main(args = ARGS)
    output = get(args, 1, joinpath(@__DIR__, "build", "results.json"))
    mkpath(dirname(output))

    info = (; git = git_info(), julia = julia_info())
    open(joinpath(dirname(output), "info.json"), write = true) do io
        JSON.print(io, info)
    end

    suite = ConcurrentCollectionsBenchmarks.BenchDictHistogram.setup(
        ntasks_list = 1:Threads.nthreads(),
        nkeys_list = [2^13, 2^16, 2^19, 2^25],
    )
    results = run(suite; verbose = true)
    BenchmarkTools.save(output, results)
    return results
end

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


text:
using BioSequences
using Random

function addition(x,y)
	sum = x+y
	return sum
end

println(addition(2,7))
println(addition(-1,5))

f(x) = 3*x^2

println(f(4))

println("""the sum of 2 and 7 is $(addition(2,7))""")
#prints the sum of 2 and 7 is 9

seq = randdnaseq(50)
println(seq)
println(typeof(seq))

max_stars_repo_name: KorfLab/julie

max_issues_repo_path: AlanZ/playground backup 1.jl

max_issues_repo_name: KorfLab/julie

max_forks_repo_path: AlanZ/playground backup 1.jl

max_forks_repo_name: KorfLab/julie


text:
# This file is part of the IntervalArithmetic.jl package; MIT licensed

"""
    in(x, a)
    ∈(x, a)

Checks if the number `x` is a member of the interval `a`, treated as a set.
Corresponds to `isMember` in the ITF-1788 Standard.
"""
function in(x::T, a::Interval) where T<:Real
    isinf(x) && return false
    a.lo <= x <= a.hi
end

in(x::Interval, y::Interval) = throw(ArgumentError("$x ∈ $y is not defined"))

in(x::T, a::Complex{<:Interval}) where T<:Real = x ∈ real(a) && 0 ∈ imag(a)
in(x::Complex{T}, a::Complex{<:Interval}) where T<:Real = real(x) ∈ real(a) && imag(x) ∈ imag(a)

"""
    issubset(a,b)
    ⊆(a,b)

Checks if all the points of the interval `a` are within the interval `b`.
"""
function ⊆(a::Interval, b::Interval)
    isempty(a) && return true
    b.lo ≤ a.lo && a.hi ≤ b.hi
end

function ⊂(a::Interval, b::Interval)
    a == b && return false
    a ⊆ b
end

function ⊇(a::Interval, b::Interval)
    b ⊆ a
end

function ⊃(a::Interval, b::Interval)
    b ⊂ a
end

# isinterior
function isinterior(a::Interval, b::Interval)
    isempty(a) && return true
    islessprime(b.lo, a.lo) && islessprime(a.hi, b.hi)
end
const ⪽ = isinterior  # \subsetdot

# Disjoint:
function isdisjoint(a::Interval, b::Interval)
    (isempty(a) || isempty(b)) && return true
    islessprime(b.hi, a.lo) || islessprime(a.hi, b.lo)
end

function isdisjoint(a::Complex{<:Interval}, b::Complex{<:Interval})
    isdisjoint(real(a),real(b)) || isdisjoint(imag(a),imag(b))
end

# Intersection
"""
    intersect(a, b)
    ∩(a,b)

Returns the intersection of the intervals `a` and `b`, considered as
(extended) sets of real numbers. That is, the set that contains
the points common in `a` and `b`.
"""
function intersect(a::Interval{T}, b::Interval{T}) where T

    x = Interval(max(a.lo, b.lo), min(a.hi, b.hi))

    if x.lo > x.hi
        return emptyinterval(T)
    end

    return x
end
# Specific promotion rule for intersect:
intersect(a::Interval{T}, b::Interval{S}) where {T,S} =
    intersect(promote(a, b)...)

function intersect(a::Complex{Interval{T}},b::Complex{Interval{T}}) where T
    isdisjoint(a,b) && return complex(emptyinterval(T),emptyinterval(T)) # for type stability

    complex(intersect(real(a),real(b)),intersect(imag(a),imag(b)))
end
intersect(a::Interval{Complex{T}}, b::Interval{Complex{S}}) where {T,S} =
    intersect(promote(a, b)...)

"""
    intersect(a::Interval{T}...) where T

Returns the n-ary intersection of its arguments.

This function is applicable to any number of input intervals, as in
`intersect(a1, a2, a3, a4)` where `ai` is an interval.
If your use case needs to splat the input, as in `intersect(a...)`, consider
`reduce(intersect, a)` instead, because you save the cost of splatting.
"""
function intersect(a::Interval{T}...) where T
    low = maximum(broadcast(ai -> ai.lo, a))
    high = minimum(broadcast(ai -> ai.hi, a))

    !is_valid_interval(low, high) && return emptyinterval(T)
    return Interval(low, high)
end

## Hull
"""
    hull(a, b)

Returns the "interval hull" of the intervals `a` and `b`, considered as
(extended) sets of real numbers, i.e. the smallest interval that contains
all of `a` and `b`.
"""
hull(a::Interval, b::Interval) = Interval(min(a.lo, b.lo), max(a.hi, b.hi))
#
# hull{T,S}(a::Interval{T}, b::Interval{S}) = hull(promote(a, b)...)
hull(a::Complex{<:Interval},b::Complex{<:Interval}) =
    complex(hull(real(a),real(b)),hull(imag(a),imag(b)))
hull(a...) = reduce(hull, a)
hull(a::Vector{Interval{T}}) where {T} = reduce(hull, a)

"""
    union(a, b)
    ∪(a,b)

Returns the union (convex hull) of the intervals `a` and `b`; it is equivalent
to `hull(a,b)`.
"""
union(a::Interval, b::Interval) = hull(a, b)
#
# union(a::Interval, b::Interval) = union(promote(a, b)...)
union(a::Complex{<:Interval},b::Complex{<:Interval}) = hull(a, b)
union(a...) = reduce(union, a)

"""
    setdiff(x::Interval, y::Interval)

Calculate the set difference `x ∖ y`, i.e. the set of values
that are inside the interval `x` but not inside `y`.

Returns an array of intervals.
The array may:

- be empty if `x ⊆ y`;
- contain a single interval, if `y` overlaps `x`
- contain two intervals, if `y` is strictly contained within `x`.
"""
function setdiff(x::Interval, y::Interval)
    intersection = x ∩ y

    isempty(intersection) && return [x]
    intersection == x && return typeof(x)[]  # x is subset of y; setdiff is empty

    x.lo == intersection.lo && return [Interval(intersection.hi, x.hi)]
    x.hi == intersection.hi && return [Interval(x.lo, intersection.lo)]

    return [Interval(x.lo, y.lo),
            Interval(y.hi, x.hi)]

end


text:
#=
the IDX file format is a simple format for vectors and multidimensional matrices
of various numerical types.
The basic format is

magic number
size in dimension 0
size in dimension 1
size in dimension 2
.....
size in dimension N
data

The magic number is an integer (MSB first). The first 2 bytes are always 0.

The third byte codes the type of the data:
0x08: unsigned byte
0x09: signed byte
0x0B: short (2 bytes)
0x0C: int (4 bytes)
0x0D: float (4 bytes)
0x0E: double (8 bytes)

The 4-th byte codes the number of dimensions of the vector/matrix: 1 for
vectors, 2 for matrices....

The sizes in each dimension are 4-byte integers (MSB first, high endian, like in
most non-Intel processors).

The data is stored like in a C array, i.e. the index in the last dimension
changes the fastest.
=#

module MNIST

export read_file_from_index

using CodecZlib

TRAINING_SAMPLES_PATH = "datasets/mnist/train-images-idx3-ubyte.gz"
TRAINING_LABELS_PATH = "datasets/mnist/train-labels-idx1-ubyte.gz"

TEST_SAMPLES_PATH = "datasets/mnist/t10k-images-idx3-ubyte.gz"
TEST_LABELS_PATH = "datasets/mnist/t10k-labels-idx1-ubyte.gz"

"""
    read_file_from_index(filename; sample_index, num_samples)

Reads samples from `filename`.

- `sample_index`: index of the first sample to read, or the beginning of the
  file if not provided.
- `num_samples`: the number of samples to read, or the entire file if not
  provided.

Returns a vector of tensors, each of which represents a single sample.
"""
function read_file_from_index(filename::String;
                              sample_index=missing, num_samples=missing)
    open(GzipDecompressorStream, filename, "r") do stream
        tensor_dims, current_pos = read_header(stream)

        # tensor_dims[1] is assumed to represent the number of samples, with the
        # remaining fields (of which there may be none) representing the number
        # of dimensions and the size along each dimension.
        total_num_samples = tensor_dims[1]
        sample_dims = tensor_dims[2:end]

        if ismissing(sample_index) && ismissing(num_samples)
            # Read all samples from the beginning of the data.
            samples, current_pos = read_samples(stream, sample_dims, total_num_samples)

            # Make sure that we've actually read the whole file.
            seekend(stream)
            end_pos = position(stream)
            @assert current_pos == end_pos
        elseif !ismissing(sample_index) && !ismissing(num_samples)
            # Read num_samples from the offset, if valid.
            @assert 1 <=
                sample_index <=
                sample_index + num_samples - 1 <=
                total_num_samples

            target_offset = (sample_index - 1) * prod(sample_dims)
            skip(stream, target_offset)

            samples, current_pos = read_samples(stream, sample_dims, num_samples)
        else
            error("sample_index and num_samples must both be provided, or neither at all.")
        end

        return samples
    end
end

"""
    read_header(stream)

Reads the IDX header from a file stream.

Returns the tensor dimensions of the data, and the stream position of the
beginning of the data; also seeks the stream to that position.
"""
function read_header(stream::IO)
    seekstart(stream)

    magic_number = maybe_bswap(read(stream, UInt32))
    @assert (magic_number & 0xFFFF0000) == 0

    datatype = (magic_number >>> 8) & 0x000000FF
    tensor_order = magic_number & 0x000000FF

    tensor_dims = Vector{UInt32}(undef, tensor_order)
    for i in 1:tensor_order
        tensor_dims[i] = maybe_bswap(read(stream, UInt32))
    end

    return Tuple(tensor_dims), position(stream)
end

"""
    read_samples(stream, sample_dims, num_samples)

Reads `num_samples` samples from `stream`; does not perform bounds checking.

- `stream`: Stream from which to read, assumed to already be at the correct
  offset.
- `sample_dims`: Tuple containing the size of each sample along each axis;
  length(sample_dims) represents the number of axes.
- `num_samples`: Number of samples to read.

Returns a vector of samples, and the updated stream position after reading the
samples. Also seeks the stream to that position.
"""
function read_samples(stream::IO, sample_dims::Tuple, num_samples::Integer)
    # Each sample has tensor order length(sample_dims).
    samples = Vector{Array{UInt8, length(sample_dims)}}(undef, num_samples)

    for i in 1:num_samples
        samples[i] = read_as_sample(stream, sample_dims)
    end

    return samples, position(stream)
end

"""
    read_as_sample(stream, sample_dims)

Reads a single sample from `stream`; does not perform bounds checking.

- `stream`: Stream from which to read, assumed to already be at the correct
  offset; seeking on each individual sample would be expensive, since
  decompressed streams may not support arbitrary seeking.
- `sample_dims`: Tuple containing the size along each axis; length(sample_dims)
  represents the number of axes.
"""
function read_as_sample(stream::IO, sample_dims::Tuple)
    # Base.prod(...) nicely handles empty tuples by returning 1 (multiplicative
    # identity).
    return reshape(read(stream, prod(sample_dims)), # Vector
                   sample_dims) # Reshaped to a length(sample_dims)-order tensor
end

"""
    maybe_bswap(bytes)

Return `bytes` in reversed order if on a little-endian machine, otherwise return
`bytes`

Guarantees that the output is big-endian.
"""
function maybe_bswap(bytes)
    return ENDIAN_BOM == 0x04030201 ? bswap(bytes) : bytes
end

"""
    get_datatype(bytes)

Returns a string representing the datatype of the data.

Throws an error if the mapping is invalid.
"""
function get_datatype(bytes)
    datatypes = Dict(
        0x08=>"unsigned byte",
        0x09=>"signed byte",
        0x0B=>"short",
        0x0C=>"int",
        0x0D=>"float",
        0x0E=>"double"
    )
    return datatypes[bytes]
end

end # module

max_stars_repo_name: brianxie/automaton-edification

max_issues_repo_name: brianxie/automaton-edification

max_forks_repo_name: brianxie/automaton-edification


text:
abstract type pseudoGPUArray <: SolidStateDetectors.GPUArrays.AbstractGPUArray{T where T, N where N} end

@testset "Simulate example detector: Inverted Coax" begin
    sim = Simulation{T}(SSD_examples[:InvertedCoax])

    calculate_electric_potential!(
        sim,
        device_array_type = pseudoGPUArray,
        convergence_limit = 1e-5,
        refinement_limits = [0.2, 0.1],
        depletion_handling = true
    )

    W = SolidStateDetectors.calculate_stored_energy(sim)
    @test isapprox(ustrip(W), 2.41e-4, atol = 1e-5)
end

@testset "Simulate example detector: Inverted Coax" begin
    sim = Simulation{T}(SSD_examples[:CGD])

    calculate_electric_potential!(
        sim,
        device_array_type = pseudoGPUArray,
        convergence_limit = 1e-5,
        refinement_limits = [0.2],
        depletion_handling = true
    )

    W = SolidStateDetectors.calculate_stored_energy(sim)
    @test isapprox(ustrip(W), 4.39e-6, atol = 1e-5)
end


text:
function get_tblis_data(N, C)

    open(joinpath(@__DIR__, "N$N", "tblis_c$(C)_t$(N).out"), "r") do lines

        for line = eachline(lines)
            m = match(r"Average time per iteration\s+?(\d+?\.\d+)", line)
            if m !== nothing
                return parse(Float64, m.captures[1])
            end
        end
    end
end


text:
"""
The BoundingBox type holds two Points, `corner1` and `corner2`.

    BoundingBox(;centered=true)     # the bounding box of the Drawing
    BoundingBox(s::AbstractString)  # the bounding box of a text string
    BoundingBox(pt::Array)          # the bounding box of a polygon

`BoundingBox(;centered=true)` returns a BoundingBox the same size and position as the current drawing, assuming
the origin (0, 0) is at the center.

The `centered` option defaults to `true`, and assumes the drawing is currently
centered. If `false`, the function assumes that the origin is at the top left
of the drawing. So this function doesn't really work if the current matrix has
been modified (by `translate()`, `scale()`, `rotate()` etc.)
"""
mutable struct BoundingBox
   corner1::Point
   corner2::Point
end

function Base.show(io::IO, bbox::BoundingBox)
  print(io, " ⤡ ", bbox.corner1, " : ", bbox.corner2)
end

function BoundingBox(;centered=true)
    if centered
        # ignore current matrix
        r = BoundingBox(
        Point(-current_width()/2, -current_height()/2),
        Point(current_width()/2, current_height()/2))
    else
        b = getmatrix()
        setmatrix([1.0, 0.0, 0.0, 1.0, 0.0, 0.0])
        r = BoundingBox(Point(0, 0), Point(current_width(), current_height()))
        setmatrix(b)
    end
    return r
end

"""
    BoundingBox(pointlist::Array)

Return the BoundingBox of a polygon (array of points).
"""
function BoundingBox(pointlist::Array{Point, 1})
    lowx, lowy = pointlist[1].x, pointlist[1].y
    highx, highy = pointlist[end].x, pointlist[end].y
    for p in pointlist
        p.x < lowx  && (lowx  = p.x)
        p.y < lowy  && (lowy  = p.y)
        p.x > highx && (highx = p.x)
        p.y > highy && (highy = p.y)
    end
    return BoundingBox(Point(lowx, lowy), Point(highx, highy))
end

"""
    BoundingBox(str::AbstractString)

Return a BoundingBox that just encloses a text string, given the current font
selection.
"""
function BoundingBox(str::AbstractString)
    xbearing, ybearing, width, height, xadvance, yadvance = textextents(str)
    lcorner = Point(xbearing, ybearing)
    ocorner = Point(lcorner.x + width, lcorner.y + height)
    return BoundingBox(lcorner, ocorner)
end

#state is either 1 or 2!
function Base.iterate(bbox::BoundingBox)
    return (bbox.corner1, 2)
end

function Base.iterate(bbox::BoundingBox, state)
    if state > 2
        return
    end
    state == 1 ? (bbox.corner1, state + 1) : state == 2 ? (bbox.corner2, state + 1) : error("this should never happen")
end

Base.last(bbox::BoundingBox) = bbox.corner2

Base.size(bbox::BoundingBox) = (2)

function Base.length(bbox::BoundingBox)
    return 2
end

unsafe_length(bbox::BoundingBox) = length(bbox)

Base.eltype(::Type{BoundingBox}) = Point

Base.lastindex(bbox::BoundingBox) = 2

function Base.getindex(bbox::BoundingBox, i::Int)
    if i == 1
        return bbox.corner1
    else
        return bbox.corner2
    end
end

Base.getindex(bbox::BoundingBox, i::Number) = bbox[convert(Int, i)]
Base.getindex(bbox::BoundingBox, I) = [bbox[i] for i in I]

function Base.setindex(bbox::BoundingBox, v::Point, i::Number)
    if i == 1
        return BoundingBox(v, bbox.corner2)
    else
        return BoundingBox(bbox.corner1, v)
    end
end

# basics
# adding a number to a box makes the box larger all round
+(z::Number, bb1::BoundingBox)              = BoundingBox(bb1.corner1 - z,          bb1.corner2 + z)
+(bb1::BoundingBox, z::Number)              = BoundingBox(bb1.corner1 - z,          bb1.corner2 + z)

# subtracting a number makes the box smaller all round
-(z::Number, bb1::BoundingBox)              = BoundingBox(bb1.corner1 + z,          bb1.corner2 - z)
-(bb1::BoundingBox, z::Number)              = BoundingBox(bb1.corner1 + z,          bb1.corner2 - z)

function (*)(bb::BoundingBox, s::Real)
    c = midpoint(bb...)
    newcorner1 = between(c, bb.corner1, s)
    newcorner2 = between(c, bb.corner2, s)
    return BoundingBox(newcorner1, newcorner2)
end

(*)(s::Real, bb::BoundingBox) = bb * s

function (/)(bb::BoundingBox, s::Real)
    c = midpoint(bb...)
    newcorner1 = between(c, bb.corner1, 1/s)
    newcorner2 = between(c, bb.corner2, 1/s)
    return BoundingBox(newcorner1, newcorner2)
end

(/)(s::Real, bb::BoundingBox) = bb / s

# modifying with tuples shifts the box in x and y
+(bb1::BoundingBox, shift::NTuple{2, Real}) = BoundingBox(
    Point(bb1.corner1.x + shift[1], bb1.corner1.y + shift[2]),
    Point(bb1.corner2.x + shift[1], bb1.corner2.y + shift[2]))

-(bb1::BoundingBox, shift::NTuple{2, Real}) = BoundingBox(
    Point(bb1.corner1.x - shift[1], bb1.corner1.y - shift[2]),
    Point(bb1.corner2.x - shift[1], bb1.corner2.y - shift[2]))

# if you add two bounding boxes, you get the biggest enclosing box
function +(bb1::BoundingBox, bb2::BoundingBox)
    pts = [bb1.corner1, bb1.corner2, bb2.corner1, bb2.corner2]
    lox, hix = sort!(pts, by = a -> a.x, lt  = (a, b) -> a < b)[[1, end]]
    loy, hiy = sort!(pts, by = a -> a.y, lt  = (a, b) -> a < b)[[1, end]]
    return BoundingBox(Point(lox.x, loy.y), Point(hix.x, hiy.y))
end

"""
    boxwidth(bb::BoundingBox)

Return the width of bounding box `bb`.
"""
boxwidth(bb::BoundingBox)       = abs(bb.corner1.x - bb.corner2.x)

"""
    boxheight(bb::BoundingBox)

Return the height of bounding box `bb`.
"""
boxheight(bb::BoundingBox)      = abs(bb.corner1.y - bb.corner2.y)

"""
    boxdiagonal(bb::BoundingBox)

Return the length of the diagonal of bounding box `bb`.
"""
boxdiagonal(bb::BoundingBox)    = hypot(boxwidth(bb), boxheight(bb))

"""
    boxaspectratio(bb::BoundingBox)

Return the aspect ratio (the height divided by the width) of bounding box `bb`.
"""
boxaspectratio(bb::BoundingBox) = boxheight(bb)/boxwidth(bb)

"""
    convert(Point, bbox::BoundingBox)

Convert a BoundingBox to a four-point clockwise polygon.

    convert(Vector{Point}, BoundingBox())
"""
function Base.convert(::Type{Vector{Point}}, bbox::BoundingBox)
    return [bbox.corner1,
            Point(bbox.corner2.x, bbox.corner1.y),
            bbox.corner2,
            Point(bbox.corner1.x, bbox.corner2.y)]
end

"""
    box(bbox::BoundingBox, :action;
            vertices=false)

Make a box using the bounds in `bbox`.

Use `vertices=true` to return an array of the four corner points: bottom left,
top left, top right, bottom right.
"""
function box(bbox::BoundingBox, action::Symbol=:nothing;
        vertices=false)
    if vertices
        botleft  = Point(bbox.corner1.x, bbox.corner2.y)
        topleft  = bbox.corner1
        topright = Point(bbox.corner2.x, bbox.corner1.y)
        botright = bbox.corner2
        return [botleft, topleft, topright, botright]
    else
        box(bbox.corner1, bbox.corner2, action)
    end
end

"""
    poly(bbox::BoundingBox, :action; kwargs...)

Make a polygon around the BoundingBox in `bbox`.
"""
poly(bbox::BoundingBox, action::Symbol=:nothing; kwargs...) =
    poly(convert(Vector{Point}, bbox), action; kwargs...)

"""
    prettypoly(bbox::BoundingBox, :action; kwargs...)

Make a decorated polygon around the BoundingBox in `bbox`. The vertices are in
the order: bottom left, top left, top right, and bottom right.
"""
prettypoly(bbox::BoundingBox, action::Symbol=:nothing; kwargs...) =
    prettypoly(box(bbox, vertices=true), action; kwargs...)

"""
    boundingboxesintersect(bbox1::BoundingBox, bbox2::BoundingBox)
    boundingboxesintersect(acorner1::Point, acorner2::Point, bcorner1::Point, bcorner2::Point)

Return true if the two bounding boxes intersect.
"""
function boundingboxesintersect(acorner1::Point, acorner2::Point, bcorner1::Point, bcorner2::Point)
    minax, maxax = minmax(acorner1.x, acorner2.x)
    minay, maxay = minmax(acorner1.y, acorner2.y)

    minbx, maxbx = minmax(bcorner1.x, bcorner2.x)
    minby, maxby = minmax(bcorner1.y, bcorner2.y)

    maxax < minbx && return false
    maxbx < minax && return false

    maxay < minby && return false
    maxby < minay && return false
    return true # boxes overlap
end

boundingboxesintersect(bbox1::BoundingBox, bbox2::BoundingBox) =
    boundingboxesintersect(bbox1.corner1, bbox1.corner2, bbox2.corner1, bbox2.corner2)

"""
    intersectboundingboxes(bb1::BoundingBox, bb2::BoundingBox)

Return a BoundingBox that's an intersection of the two bounding boxes.
"""
function intersectboundingboxes(bb1::BoundingBox, bb2::BoundingBox)
    !boundingboxesintersect(bb1, bb2) && return BoundingBox(O, O)
    minax, maxax = minmax(bb1.corner1.x, bb1.corner2.x)
    minay, maxay = minmax(bb1.corner1.y, bb1.corner2.y)

    minbx, maxbx = minmax(bb2.corner1.x, bb2.corner2.x)
    minby, maxby = minmax(bb2.corner1.y, bb2.corner2.y)
    return BoundingBox(Point(max(minax, minbx),
                             max(minay, minby)),
                       Point(min(maxax, maxbx),
                             min(maxay, maxby)))
end

"""
    isinside(p::Point, bb:BoundingBox)

Returns `true` if `pt` is inside bounding box `bb`.
"""
isinside(p::Point, bb::BoundingBox) = (bb.corner1.x <= p.x <= bb.corner2.x) &&
    (bb.corner1.y <= p.y <= bb.corner2.y)

"""
    midpoint(bb::BoundingBox)

Returns the point midway between the two points of the BoundingBox. This should
also be the center, unless I've been very stupid...
"""
midpoint(bb::BoundingBox) = midpoint(bb...)

"""
    between(bb::BoundingBox, x)

Find a point between the two corners of a BoundingBox corresponding to `x`,
where `x` is typically between 0 and 1.
"""
between(bb::BoundingBox, k=0.5) = between(bb[1], bb[2], k)

"""
    boxtopleft(bb::BoundingBox)

Return the point at the top left of the BoundingBox `bb`.

```
• ⋅ ⋅
⋅ ⋅ ⋅
⋅ ⋅ ⋅
```

"""
boxtopleft(bb::BoundingBox)        = bb[1]

"""
    boxtopcenter(bb::BoundingBox)

Return the point at the top center of the BoundingBox `bb`.
```
⋅ • ⋅
⋅ ⋅ ⋅
⋅ ⋅ ⋅
```

"""
boxtopcenter(bb::BoundingBox) = midpoint(bb.corner1, bb.corner2) - (0, boxheight(bb)/2)

"""
    boxtopright(bb::BoundingBox)

Return the point at the top right of the BoundingBox `bb`.
```
⋅ ⋅ •
⋅ ⋅ ⋅
⋅ ⋅ ⋅
```

"""
boxtopright(bb::BoundingBox)       = Point(bb[2].x, bb[1].y)

"""
    boxmiddleleft(bb::BoundingBox)

Return the point at the middle left of the BoundingBox `bb`.
```
⋅ ⋅ ⋅
• ⋅ ⋅
⋅ ⋅ ⋅
```

"""
boxmiddleleft(bb::BoundingBox)     = Point(bb[1].x, midpoint(bb[1], bb[2]).y)

"""
    boxmiddlecenter(bb::BoundingBox)

Return the point at the center of the BoundingBox `bb`.
```
⋅ ⋅ ⋅
⋅ • ⋅
⋅ ⋅ ⋅
```
"""
boxmiddlecenter(bb::BoundingBox)   = midpoint(bb[1], bb[2])

"""
    boxmiddleright(bb::BoundingBox)

Return the point at the midde right of the BoundingBox `bb`.
```
⋅ ⋅ ⋅
⋅ ⋅ •
⋅ ⋅ ⋅
```
"""
boxmiddleright(bb::BoundingBox)    = Point(bb[2].x, midpoint(bb[1], bb[2]).y)

"""
    boxbottomleft(bb::BoundingBox)

Return the point at the bottom left of the BoundingBox `bb`.
```
⋅ ⋅ ⋅
⋅ ⋅ ⋅
• ⋅ ⋅
```
"""
boxbottomleft(bb::BoundingBox)     = Point(bb[1].x, bb[2].y)

"""
    boxbottomcenter(bb::BoundingBox)

Return the point at the bottom center of the BoundingBox `bb`.
```
⋅ ⋅ ⋅
⋅ ⋅ ⋅
⋅ • ⋅
```
"""
boxbottomcenter(bb::BoundingBox)   = midpoint(bb.corner1, bb.corner2) + (0, boxheight(bb)/2)

"""
    boxbottomright(bb::BoundingBox)

Return the point at the bottom right of the BoundingBox `bb`.
```
⋅ ⋅ ⋅
⋅ ⋅ ⋅
⋅ ⋅ •
```
"""
boxbottomright(bb::BoundingBox)    = bb[2]

# legacy defs will be deprecated
boxtop(bb::BoundingBox) = boxtopcenter(bb)
boxbottom(bb::BoundingBox) = boxbottomcenter(bb)

"""
    pointcrossesboundingbox(pt, bbox::BoundingBox)

Find and return the point where a line from the center of bounding box
`bbox` to point `pt` would, if continued, cross the edges of the box.
"""
function pointcrossesboundingbox(pt, bbox::BoundingBox)
    minpt, maxpt = extrema(bbox)
    mp = midpoint(minpt, maxpt)
    midX, midY = mp.x, mp.y

    m = (midY - pt.y) / (midX - pt.x)
    if pt.x <= midX # left
        min_Xy = m * (minpt.x - pt.x) + pt.y
        if minpt.y <= min_Xy && min_Xy <= maxpt.y
            return Point(minpt.x, min_Xy)
        end
    end
    if pt.x >= midX # right
        max_Xy = m * (maxpt.x - pt.x) + pt.y
        if minpt.y <= max_Xy && max_Xy <= maxpt.y
            return Point(maxpt.x, max_Xy)
        end
    end
    if pt.y <= midY # top
        min_Yx = (minpt.y - pt.y) / m + pt.x
        if minpt.x <= min_Yx && min_Yx <= maxpt.x
            return Point(min_Yx, minpt.y)
        end
    end
    if pt.y >= midY # bottom
        max_Yx = (maxpt.y - pt.y) / m + pt.x
        if minpt.x <= max_Yx && max_Yx <= maxpt.x
            return Point(max_Yx, maxpt.y)
        end
    end
    if pt.x == midX && pt.y == midY
        return Point(pt.x, pt.y)
    end
    error("oh dear, something has gone wrong but I don't know what")
end


text:
import Base: quantile, convert
import Distributions.pnormalize!

abstract type InferenceAlgorithm end
abstract type SamplingAlgorithm <: InferenceAlgorithm end
abstract type MCMCAlgorithm <: SamplingAlgorithm end

## Types of algorithm
# - Preimage Partition: pre_bfs
# - Exact Preimage Sample: ___
# - Approximate Preimage Sample: pre_tlmh
# - Exact Conditional Abstract Sample: __
# - Approximate Conditional Abstract Sample: ___
# - Exact Conditional Point Sample
# - Approximate Conditional Point Sample

for finame in ["prob.jl",
               "rand.jl",
               "model.jl"]
    include(joinpath("query", finame))
end


text:
push!(LOAD_PATH, "../src/")

using Documenter
using BioEnergeticFoodWebs

makedocs(
         modules = [BioEnergeticFoodWebs],
         sitename = "BioEnergeticFoodWebs.jl"
         )

deploydocs(
        deps   = Deps.pip("pygments", "python-markdown-math"),
        repo   = "github.com/PoisotLab/BioEnergeticFoodWebs.jl.git",
        devbranch = "next",
        push_preview = true
        )


text:
import Base: copy, +

#*********************************#
#       COST FUNCTION CLASS       #
#*********************************#

"""
Abstract type that represents a scalar-valued function that accepts a state and control
at a single knot point.
"""
abstract type CostFunction end

import RobotDynamics: diffmethod, ForwardAD, FiniteDifference
diffmethod(::CostFunction) = ForwardAD()

# Automatic differentiation methods for generic costs
@inline function gradient!(E, cost::CostFunction, z::AbstractKnotPoint, cache=ExpansionCache(cost))
    _gradient!(diffmethod(cost), E, cost, z, cache)
end

@inline function hessian!(E, cost::CostFunction,  z::AbstractKnotPoint, cache=ExpansionCache(cost))
    _hessian!(diffmethod(cost), E, cost, z, cache)
end

# ForwardDiff methods
function _gradient!(::ForwardAD, E, cost::CostFunction, z::AbstractKnotPoint, cache=nothing)
    if is_terminal(z)
        costfun_term(x) = stage_cost(cost, x)
        ForwardDiff.gradient!(E.x, costfun_term, state(z))
    else
        ix,iu = z._x, z._u
        costfun(z) = stage_cost(cost, z[ix], z[iu])
        ForwardDiff.gradient!(E.grad, costfun, z.z)
    end
    return false
end

function _hessian!(::ForwardAD, E, cost::CostFunction, z::AbstractKnotPoint, cache=nothing)
    if is_terminal(z)
        costfun_term(x) = stage_cost(cost, x)
        ForwardDiff.hessian!(E.xx, costfun_term, state(z))
    else
        ix,iu = z._x, z._u
        costfun(z) = stage_cost(cost, z[ix], z[iu])
        ForwardDiff.hessian!(E.hess, costfun, z.z)
    end
    return false
end

# FiniteDiff methods
function _gradient!(::FiniteDifference, E, cost::CostFunction, z::AbstractKnotPoint,
        cache=ExpansionCache(cost))
    if is_terminal(z)
        cache = cache[3]
        costfun_term(x) = stage_cost(cost, x)
        cache.c3 .= state(z)
        FiniteDiff.finite_difference_gradient!(E.x.data, costfun_term, cache.c3, cache)
    else
        cache = cache[1]
        ix,iu = z._x, z._u
        costfun(z) = stage_cost(cost, z[ix], z[iu])
        cache.c3 .= z.z
        FiniteDiff.finite_difference_gradient!(E.grad, costfun, cache.c3, cache)
    end
    return false
end

function _hessian!(::FiniteDifference, E, cost::CostFunction, z::AbstractKnotPoint,
        cache=ExpansionCache(cost))
    if is_terminal(z)
        cache = cache[4]
        costfun_term(x) = stage_cost(cost, x)
        x = state(z)
        cache.xmm .= x
        cache.xmp .= x
        cache.xpm .= x
        cache.xpp .= x
        FiniteDiff.finite_difference_hessian!(E.xx.data, costfun_term, cache.xmm, cache)
    else
        cache = cache[2]
        ix,iu = z._x, z._u
        costfun(z) = stage_cost(cost, z[ix], z[iu])
        cache.xmm .= z.z
        cache.xmp .= z.z
        cache.xpm .= z.z
        cache.xpp .= z.z
        FiniteDiff.finite_difference_hessian!(E.hess, costfun, cache.xmm, cache)
    end
    return false
end

"""
    ExpansionCache(cost)

Allocate the cache needed to evaluate the gradient and Hessian of the cost function.
Returns a 4-element vector with the following elements: `[grad, hess, grad_term, hess_term]`,
which are the caches need for the gradient and Hessian for the state and terminal cost
functions.
"""
@inline ExpansionCache(cost::CostFunction) =
    ExpansionCache(diffmethod(cost), state_dim(cost), control_dim(cost))
@inline ExpansionCache(::ForwardAD, n::Int, m::Int) = (nothing,nothing,nothing,nothing)
function ExpansionCache(::FiniteDifference, n::Int, m::Int)
    grad_cache0 = FiniteDiff.GradientCache(zeros(n), zeros(n), Val(:forward))
    hess_cache0 = FiniteDiff.HessianCache(zeros(n))
    grad_cache  = FiniteDiff.GradientCache(zeros(n+m), zeros(n+m), Val(:forward))
    hess_cache  = FiniteDiff.HessianCache(zeros(n+m))
    return (grad_cache, hess_cache, grad_cache0, hess_cache0)
end


text:
for (func, typ) in [(:clblasCherk, cl.CL_float),
                    (:clblasZherk, cl.CL_double)]

    @eval @blasfun $func(order::clblasOrder, uplo::clblasUplo,
                              trans::clblasTranspose,
                              N::Csize_t, K::Csize_t,
                              alpha::$typ,
                              A::cl.CL_mem, offa::Csize_t, lda::Csize_t,
                              beta::$typ,
                              C::cl.CL_mem, offc::Csize_t, ldc::Csize_t,
                              n_queues::cl.CL_uint,
                              queues::Ptr{cl.CL_command_queue},
                              n_events_in_wait_list::cl.CL_uint,
                              event_wait_list::Ptr{cl.CL_event},
                              events::Ptr{cl.CL_event})

    @eval @blasfun2 $func(order::clblasOrder, uplo::clblasUplo,
                               trans::clblasTranspose,
                               N::Csize_t, K::Csize_t,
                               alpha::$typ,
                               A::cl.CL_mem, offa::Csize_t, lda::Csize_t,
                               beta::$typ,
                               C::cl.CL_mem, offc::Csize_t, ldc::Csize_t)

end


text:
function varargout = correlation(varargin)
#  CORRELATION MATLAB code for correlation.fig
#       CORRELATION, by itself, creates a new CORRELATION or raises the existing
#       singleton*.
#
#       H = CORRELATION returns the handle to a new CORRELATION or the handle to
#       the existing singleton*.
#
#       CORRELATION('CALLBACK',hObject,eventData,handles,...) calls the local
#       function named CALLBACK in CORRELATION.M with the given input arguments.
#
#       CORRELATION('Property','Value',...) creates a new CORRELATION or raises the
#       existing singleton*.  Starting from the left, property value pairs are
#       applied to the GUI before correlation_OpeningFcn gets called.  An
#       unrecognized property name or invalid value makes property application
#       stop.  All inputs are passed to correlation_OpeningFcn via varargin.
#
#       *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
#       instance to run (singleton)".
#
#  See also: GUIDE, GUIDATA, GUIHANDLES

#  Edit the above text to modify the response to help correlation

#  Last Modified by GUIDE v2.5 11-Nov-2013 16:38:24

#  Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @correlation_OpeningFcn, ...
                   'gui_OutputFcn',  @correlation_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
#  End initialization code - DO NOT EDIT

#  --- Executes just before correlation is made visible.
function correlation_OpeningFcn(hObject, ~, handles, varargin)
#  This function has no output args, see OutputFcn.
#  hObject    handle to figure
#  eventdata  reserved - to be defined in a future version of MATLAB
#  handles    structure with handles and user data (see GUIDATA)
#  varargin   command line arguments to correlation (see VARARGIN)

#  Choose default command line output for correlation
handles.output = hObject;

#  Update handles structure
guidata(hObject, handles);

#  UIWAIT makes correlation wait for user response (see UIRESUME)
#  uiwait(handles.figure1);

# File: correlation.jl
# Author:
# Date: 29-June-2020

# import XL data
[names,~,returns] = getXLData('hfma_matlab_data.xlsx','10 Hedge Funds');
# calculate correlation matrix
correl = corr(returns);
# format data
data = num2cell(correl);
    for i = 1:numel(data)
        data{i} = sprintf('# 8.3f',data{i});
    end
# set size column widths
cwidth = num2cell(repmat(50,1,10));
# set table output
set(handles.table_correlation_matrix,'ColumnName',names,'ColumnWidth',cwidth,'RowName',names,'Data',data);

#  --- Outputs from this function are returned to the command line.
function varargout = correlation_OutputFcn(~, ~, handles)
#  varargout  cell array for returning output args (see VARARGOUT);
#  hObject    handle to figure
#  eventdata  reserved - to be defined in a future version of MATLAB
#  handles    structure with handles and user data (see GUIDATA)

#  Get default command line output from handles structure
varargout{1} = handles.output;


text:
#fix iterations
function lowest_rank_TAME_test(A::ThirdOrderSymTensor, B::ThirdOrderSymTensor,
                          U_0::Array{F,2},V_0::Array{F,2},
                          β::F, max_iter::Int,tol::F,α::F;profile=false) where {F <:AbstractFloat}

	@assert size(U_0,2) == size(V_0,2)

    dimension = minimum((A.n,B.n))

    best_triangle_count = -Inf
    best_x = zeros(size(U_0,1),size(V_0,1))
    best_index = -1

	r = 1
	X_k = U_0 * V_0'
	normalization_factor = norm(X_k)
	#TODO: fix factorization
	X_k ./= normalization_factor

	X_0 = copy(X_k)
	U_0 ./= normalization_factor
#	V_0 ./= sqrt(normalization_factor)

	U_k = deepcopy(U_0)
	V_k = deepcopy(V_0)

	All_vectors = zeros(A.n,B.n,max_iter)

	i = 1
    lambda = Inf

    if profile
        experiment_profile = Dict(
            "ranks"=>[],
            "contraction_timings"=>[],
            "svd_timings"=>[],
            "matching_timings"=>[],
            "scoring_timings"=>[]
        )
    end

    for _ in 1:max_iter

#		println("U_k",U_k[1:2,:])
		println("size of U_0",size(U_0))
		println("size of U_k:",size(U_k))

		#println("|U_k - U_0|_F = ",norm(U_k - U_0))
        if profile
            #X_k_1,t = @timed kron_contract(A,B,U_k,V_k)
		    (A_comps, B_comps),t = @timed get_kron_contract_comps(A,B,U_k,V_k)
            push!(experiment_profile["contraction_timings"],t)
        else
            A_comps, B_comps = get_kron_contract_comps(A,B,U_k,V_k)
			# = size(U_k,2)^2 + size(U_k,2) + size(U_0,2)
			X_k_1 = kron_contract(A,B,sparse(X_k),r)
			println("|X_{k_1} - T_A*T_B'|_F=",norm(X_k_1 - A_comps*B_comps')/norm(X_k_1))
        end

		if α != 1.0 && β != 0.0
			U_temp = hcat(sqrt(α) * A_comps, sqrt(α) * sqrt(β) * U_k, sqrt(1-α) * U_0)
			V_temp = hcat(sqrt(α) * B_comps, sqrt(α) * sqrt(β) * V_k, sqrt(1-α) * V_0)
#			U_temp = hcat(α * A_comps, α * β * U_k, (1-α) * U_0)
#			V_temp = hcat(B_comps, V_k, V_0)
		elseif α != 1.0
			U_temp = hcat(sqrt(α)*A_comps, sqrt(1-α)*U_0)
			V_temp = hcat(sqrt(α)*B_comps, sqrt(1-α)*V_0)
		elseif β != 1.0
			U_temp = hcat(A_comps, sqrt(β) * U_k)
			V_temp = hcat(B_comps, sqrt(β) * V_k)
		else
			U_temp = A_comps
			V_temp = B_comps
		end

		println("U_temp size:",size(U_temp))
		println("U_temp rank:",rank(U_temp))
		X_k_1 = α*(X_k_1 + β*X_k) + (1-α)*X_0
	#	if i > 1
	#		return X_k_1 - U_temp*V_temp'
	#	end

		println("|(α*(X_k_1 + β*X_k) + (1-α)*X_0) - U_temp*V_temp'|_F=",norm(X_k_1 - U_temp*V_temp')/norm(X_k_1))

       # Q_A,R_A = qr(U_temp)
       # Q_B,R_B = qr(V_temp)

		#@show svdvals(U_temp*V_temp')
		#@show  svdvals(R_A*R_B')

    #    _,test_singvals,_ = svd(R_A*R_B')
     #   _,normalized_singvals,_ = svd((R_A*R_B')/norm(R_A*R_B'))

        #println("singular values of R_A':",test_singvals)
    #    println("singular values of R_A'R_B:",normalized_singvals)
#		println("trace of R_A'R_B:",sum(normalized_singvals))
       #println("singular values of R_B':",RS_B)

		 if profile
            #X_k_1,t = @timed kron_contract(A,B,U_k,V_k)
			(result_A),t_A = @timed svd(U_temp)
			#println(typeof(result_A))
			A_U,A_S,A_Vt = result_A.U,result_A.S,result_A.Vt
			(result_B),t_B = @timed svd(V_temp)
			B_U,B_S,B_Vt = result_B.U,result_B.S,result_B.Vt
            push!(experiment_profile["svd_timings"],t_A + t_B)
        else
			A_U,A_S,A_Vt = svd(U_temp)
			B_U,B_S,B_Vt = svd(V_temp)
		end

		#println("OG A singular values:",A_S)
		#println("OG B singular values:",B_S)
		singular_indexes_A = [i for i in 1:length(A_S) if A_S[i] > A_S[1]*eps(Float64)*size(U_temp,2)]
		singular_indexes_B = [i for i in 1:length(B_S) if B_S[i] > B_S[1]*eps(Float64)*size(V_temp,2)]

	#	println("length of singular_indexes_A:",length(singular_indexes_A))
		U_k_1 = A_U[:,singular_indexes_A]
#		U_k_1 = A_U*(diagm(A_S)*A_Vt')

		V_k_1 = V_temp*(A_Vt[:,singular_indexes_A]*diagm(A_S[singular_indexes_A]))

#		println("U_k_1 cond:",cond(U_k_1))
#		println("V_k_1 cond:",cond(V_k_1))
	#	V_k_1 ./= norm(V_k_1)

		#=
		core = diagm(A_S[singular_indexes_A])*(A_Vt[:,singular_indexes_A]'*B_Vt[:,singular_indexes_B])*diagm(B_S[singular_indexes_B])
		#println("core size is:",size(core))
		core ./= norm(core)  #normalizing the core is cheaper than the other components
		U_k_1 = A_U[:,singular_indexes_A]
		#V_k_1 = V_temp*(A_Vt[:,singular_indexes_A]*diagm(A_S[singular_indexes_A]))
		V_k_1 = B_U[:,singular_indexes_B]*core'
		=#
		#new_lambda = dot(V_k_1'*V_k, U_k_1'*U_k)

		#U_k_1 = A_U[:,singular_indexes_A]*(diagm(A_S[singular_indexes_A])*A_Vt[:,singular_indexes_A]')
		#V_k_1 = B_U[:,singular_indexes_B]*(diagm(B_S[singular_indexes_B])*B_Vt[:,singular_indexes_B]')

		#println(norm(X_k_1))

#		normalization_factor = norm(X_k_1)
		#X_k_1 = U_k_1*V_k_1'

		C = norm(X_k_1)
		println("|X_k_1 - U_k_1*V_k_1'|_F=",norm(X_k_1 - U_k_1*V_k_1')/C)
		X_k_1 ./= C

		C = norm(U_k_1*V_k_1')
		U_k_1 ./= sqrt(C)
		V_k_1 ./= sqrt(C)

		All_vectors[:,:,i] = U_k_1*V_k_1'

		r = rank(X_k_1)
		println("rank of X_k_1:       ",r)
		println("rank of U_k_1*V_k_1':",rank(All_vectors[:,:,i]))
#		V_k_1 ./= sqrt(normalization_factor)
#		X_k_1 ./= normalization_factor

		lam = dot(X_k_1,kron_contract(A,B,U_k_1,V_k_1))
		#println("old lambda:$new_lambda")
#

	#   NOTE:   triangle_matching is turned off for speed
# 		println(norm(X_k_1))
	#	sparse_X_k_1 = sparse(X_k_1)

        if profile
            triangles, gaped_triangles, matching_time, scoring_time = TAME_score(A,B,sparse_X_k_1;return_timings=true)
            push!(experiment_profile["matching_timings"],matching_time)
            push!(experiment_profile["scoring_timings"],scoring_time)
        else
     #       triangles, gaped_triangles = TAME_score(A,B,sparse_X_k_1)
        end

      #=  println("finished iterate $(i):tris:$(triangles) -- gaped_t:$(gaped_triangles)")

        if triangles > best_triangle_count
            best_x = deepcopy(X_k_1)
            best_triangle_count = triangles
            best_iterate = i
        end
		=#
        println("λ_i: $(lam)")

        if abs(lam - lambda) < tol || i >= max_iter
            if profile
                return best_x, best_triangle_count, experiment_profile
            else
                return best_x, best_triangle_count, All_vectors
            end
        end

		#get the low rank factorization for the next one
		U_k = copy(U_k_1)
		V_k = copy(V_k_1)
		X_k = copy(X_k_1)

		lambda = lam
		i += 1

		println("U_k:",U_k)
		println("norm of U_k*V_k' at end:",norm(U_k*V_k'))
		#println("norm of V_k at end:",norm(V_k))
		println("norm of X_k at end:",norm(X_k))

    end

end

function low_rank_TAME_contraction(A::ThirdOrderSymTensor, B::ThirdOrderSymTensor,
                                   U_k::Array{F,2},V_k::Array{F,2},α::Float64,β::Float64,
                                   U_0::Array{F,2},V_0::Array{F,2};profile=false) where {F <:AbstractFloat}

	if profile
		#X_k_1,t = @timed kron_contract(A,B,U_k,V_k)
		(A_comps, B_comps),t = @timed get_kron_contract_comps(A,B,U_k,V_k)
#		push!(experiment_profile["contraction_timings"],t)
	else
		A_comps, B_comps = get_kron_contract_comps(A,B,U_k,V_k)
	end

	#add in shift and mixing terms
	if α != 1.0 && β != 0.0
		U_temp = hcat(sqrt(α)*A_comps,sqrt(α*β)*U_k, sqrt(1-α)*U_0)
		#V_temp = hcat(B_comps,V_k,V_0)
		V_temp = hcat(sqrt(α)*B_comps,sqrt(α*β)*V_k, sqrt(1-α)*V_0)
	elseif α != 1.0
		U_temp = hcat(sqrt(α)*A_comps, sqrt(1-α)*U_0)
		V_temp = hcat(sqrt(α)*B_comps, sqrt(1-α)*V_0)
	elseif β != 1.0
		U_temp = hcat(A_comps, sqrt(β)*U_k)
		V_temp = hcat(B_comps, sqrt(β)*V_k)
	else
		U_temp = A_comps
		V_temp = B_comps
	end

	return U_temp , V_temp
	if profile
		#X_k_1,t = @timed kron_contract(A,B,U_k,V_k)
		(result_A),t_A = @timed svd(U_temp)
		A_U,A_S,A_Vt = result_A.U,result_A.S,result_A.Vt
		(result_B),t_B = @timed svd(V_temp)
		B_U,B_S,B_Vt = result_B.U,result_B.S,result_B.Vt
	#	push!(experiment_profile["svd_timings"],t_A + t_B)
	else
		A_U,A_S,A_Vt = svd(U_temp)
		#B_U,B_S,B_Vt = svd(V_temp)
	end
	#println(A_S)
	#println(B_S)
	singular_indexes_A = [i for i in 1:length(A_S) if A_S[i] > A_S[1]*eps(Float64)*A.n]
	#singular_indexes_B = [i for i in 1:length(B_S) if B_S[i] > B_S[1]*eps(Float64)*B.n]
	println(length(singular_indexes_A))
	#println(length(singular_indexes_B))

	#core = diagm(A_S[singular_indexes_A])*(A_Vt[:,singular_indexes_A]'*B_Vt[:,singular_indexes_B])*diagm(B_S[singular_indexes_B])
	#core ./= norm(core) #normalize
	U_k_1 = A_U[:,singular_indexes_A]
	V_k_1 = V_temp*(A_Vt[:,singular_indexes_A]*diagm(A_S[singular_indexes_A]))

	return U_k_1, V_k_1
end

#TODO: would be nice to load in processes as needed
#TODO: need a more robust way to load in the functions,
#          currently only works while running w/in directory
#TODO: check type stability
function lowest_rank_TAME_parallel_test(A::ThirdOrderSymTensor, B::ThirdOrderSymTensor,
                          U_0::Array{F,2},V_0::Array{F,2},
                          β::F, max_iter::Int,tol::F,α::F;profile=false) where {F <:AbstractFloat}

	@assert size(U_0,2) == size(V_0,2)

	matching_proc = addprocs(1)
	futures = [] #TODO: can we type this?
	#load the functions on process
	@everywhere include_string(Main,$(read("LambdaTAME.jl",String)),"LambdaTAME.jl")

    dimension = minimum((A.n,B.n))

	X_k = U_0 * V_0'
	X_0 = copy(X_k)

	U_k = copy(U_0)
	V_k = copy(V_0)

	U_k ./= norm(U_k)
	V_k ./= norm(V_k)

	i = 1
    lambda = Inf

    if profile
        experiment_profile = Dict(
            "ranks"=>[],
            "contraction_timings"=>[],
            "svd_timings"=>[],
            "matching_timings"=>[],
            "scoring_timings"=>[]
        )
    end

    while true

        if profile
            #X_k_1,t = @timed kron_contract(A,B,U_k,V_k)
		    (A_comps, B_comps),t = @timed get_kron_contract_comps(A,B,U_k,V_k)
            push!(experiment_profile["contraction_timings"],t)
        else
            A_comps, B_comps = kron_contract_test2(A,B,U_k,V_k)
        end

		U_temp = hcat(sqrt(α)*A_comps,sqrt(α*β)*U_k, sqrt(1-α)*U_0)
		V_temp = hcat(sqrt(α)*B_comps,sqrt(α*β)*V_k, sqrt(1-α)*V_0)

		 if profile
			#X_k_1,t = @timed kron_contract(A,B,U_k,V_k)
			(result_A),t_A = @timed svd(U_temp)
			#println(typeof(result_A))
			A_U,A_S,A_Vt = result_A.U,result_A.S,result_A.Vt
			(result_B),t_B = @timed svd(V_temp)
			B_U,B_S,B_Vt = result_B.U,result_B.S,result_B.Vt
			push!(experiment_profile["svd_timings"],t_A + t_B)
        else
			A_U,A_S,A_Vt = svd(U_temp)
			B_U,B_S,B_Vt = svd(V_temp)
		end
		singular_indexes_A = [i for i in 1:length(A_S) if A_S[i] > A_S[1]*eps(Float64)*dimension]
		singular_indexes_B = [i for i in 1:length(B_S) if B_S[i] > B_S[1]*eps(Float64)*dimension]

		U_k_1 = A_U[:,singular_indexes_A]*diagm(A_S[singular_indexes_A])
		V_k_1 = (B_U[:,singular_indexes_B]*(diagm(B_S[singular_indexes_B])*B_Vt[:,singular_indexes_B]'*A_Vt[:,singular_indexes_A]))

		U_k_1 ./= norm(U_k_1)
		V_k_1 ./= norm(V_k_1)

		#run Triangle matching on another process
		f = @spawnat matching_proc[1] produce_scored_matching(A,B,U_k_1,U_k_1;return_timings=true)
		#TODO: check memory usage
		push!(futures,f)

		new_lambda = dot(V_k_1'*V_k, U_k_1'*U_k)

        println("finished iterate $(i):λ_i = $new_lambda")

        println("λ: $(new_lambda)")

        if abs(new_lambda - lambda) < tol || i >= max_iter
			break
        else
			#get the low rank factorization for the next one
			U_k = copy(U_k_1)
			V_k = copy(V_k_1)

			X_k = copy(X_k_1)
			lambda = new_lambda
			i += 1
        end

    end

	#check and return the best matching
	results = [fetch(future) for future in futures]

	best_triangle_count = -Inf
	best_index = -1

	for (i,(triangle_count, gaped_triangles, bipartite_matching_time, scoring_time, U, V)) in enumerate(results)
		if profile
			push!(experiment_profile["matching_timings"],bipartite_matching_time)
			push!(experiment_profile["scoring_timings"],scoring_timings)
		end
		if triangle_count > best_triangle_count
			best_index = i
			best_triangle_count = triangle_count
		end
	end

	(_,_,_,_, best_U, best_V) = results[i]

	best_X = best_U*best_V'

	addprocs(1)
	if profile
		return best_U, best_V, best_triangle_count, experiment_profile
	else
		return best_U, best_V, best_triangle_count
	end

end

function produce_scored_matching(A::ThirdOrderSymTensor,B::ThirdOrderSymTensor,
                                 U::Array{Float64,2},V::Array{Float64,2};
								 kwargs...)

    #unsure if low rank option exists for exact solution.
	X= U*V'
	sparseX = sparse(X)
	return TAME_score(A,B,sparseX;kwargs...), U , V

end

#old code
#runs TAME, but reduces down to lowest rank form first
function lowest_rank_TAME_original_svd(A::ThirdOrderSymTensor, B::ThirdOrderSymTensor,
                          U_0::Array{F,2},V_0::Array{F,2},
                          β::F, max_iter::Int,tol::F,α::F;profile=false) where {F <:AbstractFloat}

	@assert size(U_0,2) == size(V_0,2)

    dimension = minimum((A.n,B.n))

    best_triangle_count = -Inf
    best_x = zeros(size(U_0,1),size(V_0,1))
    best_index = -1

	X_k = U_0 * V_0'
	X_0 = copy(X_k)
#    X_k = copy(W)
 #   X_k_1 = copy(W)
 	rank_0 = size(U_0,2)
	#All_vectors = zeros(A.n,B.n,max_iter)
    i = 1
    lambda = Inf

    if profile
        experiment_profile = Dict(
            "ranks"=>[],
            "contraction_timings"=>[],
            "svd_timings"=>[],
            "matching_timings"=>[],
            "scoring_timings"=>[]
        )
    end

	U_k = copy(U_0)
	V_k = copy(V_0)

    while true

        if profile
            X_k_1,t = @timed kron_contract(A,B,U_k,V_k)
            push!(experiment_profile["contraction_timings"],t)
        else
            X_k_1 = kron_contract(A,B,U_k,V_k)
        end

        new_lambda = dot(X_k_1,X_k)

        if β != 0.0
            X_k_1 .+= β * X_k
        end

        if α != 1.0
            X_k_1 = α * X_k_1 + (1 - α) * X_0
        end

        X_k_1 ./= norm(X_k_1)
		sparse_X_k_1 = sparse(X_k_1)

        if profile
            triangles, gaped_triangles, matching_time, scoring_time = TAME_score(A,B,sparse_X_k_1;return_timings=true)
            push!(experiment_profile["matching_timings"],matching_time)
            push!(experiment_profile["scoring_timings"],scoring_time)
        else
            triangles, gaped_triangles = TAME_score(A,B,sparse_X_k_1)
        end

        println("finished iterate $(i):tris:$(triangles) -- gaped_t:$(gaped_triangles)")

        if triangles > best_triangle_count
            best_x = copy(X_k_1)
            best_triangle_count = triangles
            best_iterate = i
        end

        println("λ: $(new_lambda)")

        if abs(new_lambda - lambda) < tol || i >= max_iter
            if profile
                return best_x, best_triangle_count, experiment_profile
            else
                return best_x, best_triangle_count
            end
        else

			#get the low rank factorization for the next one

			rank_k = size(U_k,2)
			rank_k_1 = rank_k^2 + rank_k + rank_0

			if profile
				(result,t) = @timed svds(sparse_X_k_1,nsv = rank_k_1)
				push!(experiment_profile["svd_timings"],t)
			else
				result = svds(sparse_X_k_1,nsv = rank_k_1)
			end

			U,S,VT = result[1]
			singular_indexes = [i for i in 1:length(S) if S[i] > S[1]*eps(Float64)*dimension]

			println("rank is",length(singular_indexes))
			if profile
				push!(experiment_profile["ranks"],length(singular_indexes))
			end

			U_k = U[:,singular_indexes]
			V_k = VT[:,singular_indexes]*diagm(S[singular_indexes])

			X_k = copy(X_k_1)
			lambda = new_lambda
			i += 1
        end

    end

end

function lowest_rank_TAME_original_svd_for_ranks(A::ThirdOrderSymTensor, B::ThirdOrderSymTensor,
                          U_0::Array{F,2},V_0::Array{F,2},
                          β::F, max_iter::Int,tol::F,α::F) where {F <:AbstractFloat}

	@assert size(U_0,2) == size(V_0,2)

    dimension = minimum((A.n,B.n))

    best_triangle_count = -Inf
    best_x = zeros(size(U_0,1),size(V_0,1))
    best_index = -1

	X_k = U_0 * V_0'
	X_0 = copy(X_k)
#    X_k = copy(W)
 #   X_k_1 = copy(W)
 	rank_0 = size(U_0,2)
	#All_vectors = zeros(A.n,B.n,max_iter)
    i = 1
    lambda = Inf

	ranks = []

	U_k = copy(U_0)
	V_k = copy(V_0)

    while true

		X_k_1 = kron_contract(A,B,U_k,V_k)

        new_lambda = dot(X_k_1,X_k)

        if β != 0.0
            X_k_1 .+= β * X_k
        end

        if α != 1.0
            X_k_1 = α * X_k_1 + (1 - α) * X_0
        end

        X_k_1 ./= norm(X_k_1)
		sparse_X_k_1 = sparse(X_k_1)

		#println("λ: $(new_lambda)")

        if abs(new_lambda - lambda) < tol || i >= max_iter
    		return best_x, ranks
        else

			#get the low rank factorization for the next one
			rank_k = size(U_k,2)
			rank_k_1 = rank_k^2 + rank_k + rank_0
			result = svds(sparse_X_k_1,nsv = rank_k_1)

			U,S,VT = result[1]
			singular_indexes = [i for i in 1:length(S) if S[i] > S[1]*eps(Float64)*dimension]

			#println("rank is",length(singular_indexes))
			push!(ranks,length(singular_indexes))

			U_k = U[:,singular_indexes]
			V_k = VT[:,singular_indexes]*diagm(S[singular_indexes])

			X_k = copy(X_k_1)
			lambda = new_lambda
			i += 1
        end

    end

end

function low_rank_TAME(A::COOTen, B::COOTen,W::Array{F,2},rank::Int,
                       β::F, max_iter::Int,tol::F,α::F) where {F <:AbstractFloat}

    #low rank factorization
    #U,S,VT = svd(W)
    #U = U[:,1:rank]
    #V = diagm(S)*VT'[:,1:rank]
    #TODO: fix this to compute top k vectors
    U,V = NMF.randinit(W,rank)
    V = Matrix(V')

    best_triangle_count,_ = TAME_score(A,B,U,V)
    best_U = copy(U)
    best_V = copy(V)
    best_index = 1

    #w = reshape(W, A.cubical_dimension,B.cubical_dimension)
    #get the low rank factors of W

    x_k = copy(W)
    i = 1
    lambda = Inf
    while true

        x_k_1 = kron_contract(A,B,U,V)
        new_lambda = dot(x_k_1,x_k)

        if β != 0.0
            println(size(x_k))
            println(size(x_k_1))
            x_k_1 .+= β * x_k
        end

        if α != 1.0
            x_k_1 = α * x_k_1 + (1 - α) * W
        end

        x_k_1 ./= norm(x_k_1)

        #low rank factorization
        #TODO: fix this to compute top k vectors

        U,V = NMF.randinit(x_k_1,rank)
        V = Matrix(V')
#        U,S,VT = svd(reshape(x_k_1,A.cubical_dimension,B.cubical_dimension))
#        U = U[:,1:rank]
#        V = diagm(S)*VT'[:,1:rank]

        triangles, gaped_triangles = TAME_score(A,B,U,V)

      #  println("finished iterate $(i):tris:$(triangles) -- gaped_t:$(gaped_triangles)")

        if triangles > best_triangle_count
            best_x = copy(x_k_1)
            best_triangle_count = triangles
            best_iterate = i
            best_U = copy(U)
            best_V = copy(V)
        end

        println("λ: $(new_lambda)")

        if abs(new_lambda - lambda) < tol || i >= max_iter
            return x_k_1, best_U, best_V, best_triangle_count
        else

            lambda = new_lambda
            i += 1
        end

    end

end

#TODO: needs to be fixed
function align_tensors(A::COOTen,B::COOTen,rank::Int,method="ΛTAME")
    iter =15
    tol=1e-6
    alphas = [.15,.5,.85]
    betas =[1000.0,100.0,10.0,1.0,0.0,0.1,0.01,0.001]

    max_triangle_match = min(size(A.indices,1),size(B.indices,1))

    TAME_timings = Array{Float64,1}(undef,length(alphas)*length(betas))
    Krylov_Search_timings = Array{Float64,1}(undef,length(alphas)*length(betas))

    U = Array{Float64,2}(undef,A.cubical_dimension,iter)
    V = Array{Float64,2}(undef,B.cubical_dimension,iter)

    best_TAME_PP_tris = -1
    best_i  = -1
    best_j = -1
    exp_index = 1

    for α in alphas
        for beta in betas

            if method == "ΛTAME"
                U,V = TAME(A,B,beta,iter,tol,α)

            else
                ((_, best_U, best_V, best_triangle_count),runtime) = @timed low_rank_TAME(A, B, ones(A.cubical_dimension,B.cubical_dimension),rank,beta, iter,tol,α)
            end
            TAME_timings[exp_index] = runtime
            search_tris, _  = TAME_score(A,B,best_U,Matrix(best_V'))

            if search_tris > best_TAME_PP_tris
                best_TAME_PP_tris = search_tris
                best_i = i
                best_j = j
            end
            exp_index += 1
        end

    end

    avg_TAME_timings = sum(TAME_timings)/length(TAME_timings)
    avg_Krylov_timings = sum(Krylov_Search_timings)/length(Krylov_Search_timings)

    return best_TAME_PP_tris,max_triangle_match, avg_TAME_timings, avg_Krylov_timings
end

#  lowest_rank_TAME old code
function lowest_rank_TAME(A::ThirdOrderSymTensor, B::ThirdOrderSymTensor,W::Array{F,2},
                       β::F, max_iter::Int,tol::F,α::F;profile=false) where {F <:AbstractFloat}

    dimension = minimum((A.n,B.n))
    ranks = []

    if profile
        experiment_profile = Dict(
            "ranks"=>[],
            "contraction_timings"=>[],
            "svd_timings"=>[],
            "matching_timings"=>[],
            "hungarian_timings"=>[]
        )
    end

    best_triangle_count = -Inf
    best_x = copy(W)
    best_index = -1
    X_k = copy(W)
    X_k_1 = copy(W)
    i = 1
    lambda = Inf

    while true
        if profile
            (U,S,VT),t = @timed svd(X_k)
            singular_indexes = [i for i in 1:length(S) if S[i] > S[1]*eps(Float64)*dimension]
            push!(experiment_profile["ranks"],length(singular_indexes))
            push!(experiment_profile["svd_timings"],t)
        else
            (U,S,VT) = svd(X_k)
            singular_indexes = [i for i in 1:length(S) if S[i] > S[1]*eps(Float64)*dimension]
        end

        U = U[:,singular_indexes]
        V = VT[:,singular_indexes]*diagm(S[singular_indexes])

        if profile
            X_k_1,t = @timed kron_contract(A,B,U,V)
            push!(experiment_profile["contraction_timings"],t)
        else
            X_k_1 = kron_contract(A,B,U,V)
        end

        new_lambda = dot(X_k_1,X_k)

        if β != 0.0
            X_k_1 .+= β * X_k
        end

        if α != 1.0
            X_k_1 = α * X_k_1 + (1 - α) * W
        end

        X_k_1 ./= norm(X_k_1)

        if profile
            triangles, gaped_triangles, hungarian_time, matching_time = TAME_score(A,B,sparse(X_k_1);return_timings=true)
            push!(experiment_profile["hungarian_timings"],hungarian_time)
            push!(experiment_profile["matching_timings"],matching_time)
        else
            triangles, gaped_triangles =  TAME_score(A,B,X_k_1)
        end

        println("finished iterate $(i):tris:$(triangles) -- gaped_t:$(gaped_triangles)")

        if triangles > best_triangle_count
            best_x = copy(X_k_1)
            best_triangle_count = triangles
            best_iterate = i
            best_U = copy(U)
            best_V = copy(V)
        end
        println("λ: $(new_lambda)")

        if abs(new_lambda - lambda) < tol || i >= max_iter
            if profile
                return best_x, best_triangle_count, ranks, experiment_profile
            else
                return best_x, best_triangle_count, ranks
            end
        else
            X_k = copy(X_k_1)
            lambda = new_lambda
            i += 1
        end

    end

    #compute the number of triangles matched in the last iterate

end


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

mutable struct PrivateEndpointServiceConnectionStatus <: SwaggerModel

    function PrivateEndpointServiceConnectionStatus(;)
        o = new()
        o
    end
end # type PrivateEndpointServiceConnectionStatus

const _property_map_PrivateEndpointServiceConnectionStatus = Dict{Symbol,Symbol}()
const _property_types_PrivateEndpointServiceConnectionStatus = Dict{Symbol,String}()
Base.propertynames(::Type{ PrivateEndpointServiceConnectionStatus }) = collect(keys(_property_map_PrivateEndpointServiceConnectionStatus))
Swagger.property_type(::Type{ PrivateEndpointServiceConnectionStatus }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_PrivateEndpointServiceConnectionStatus[name]))}
Swagger.field_name(::Type{ PrivateEndpointServiceConnectionStatus }, property_name::Symbol) =  _property_map_PrivateEndpointServiceConnectionStatus[property_name]

function check_required(o::PrivateEndpointServiceConnectionStatus)
    true
end

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