text:
using JuMP, Ipopt

# Create JuMP model, using Ipopt as the solver
rocket = Model(optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0))

# Constants
# Note that all parameters in the model have been normalized
# to be dimensionless. See the COPS3 paper for more info.
h_0 = 1    # Initial height
v_0 = 0    # Initial velocity
m_0 = 1    # Initial mass
g_0 = 1    # Gravity at the surface

T_c = 3.5  # Used for thrust
h_c = 500  # Used for drag
v_c = 620  # Used for drag
m_c = 0.6  # Fraction of initial mass left at end

c     = 0.5 * sqrt(g_0 * h_0)  # Thrust-to-fuel mass
m_f   = m_c * m_0            # Final mass
D_c   = 0.5 * v_c * m_0 / g_0    # Drag scaling
T_max = T_c * g_0 * m_0        # Maximum thrust

n = 800   # Time steps

@variables(rocket, begin
    Δt ≥ 0, (start = 1/n) # Time step
    # State variables
    v[1:n] ≥ 0            # Velocity
    h[1:n] ≥ h_0          # Height
    m_f ≤ m[1:n] ≤ m_0    # Mass
    # Control
    0 ≤ T[1:n] ≤ T_max    # Thrust
end)

# Objective: maximize altitude at end of time of flight
@objective(rocket, Max, h[n])

# Initial conditions
@constraints(rocket, begin
    v[1] == v_0
    h[1] == h_0
    m[1] == m_0
    m[n] == m_f
end)

# Forces
# Drag(h,v) = Dc v^2 exp( -hc * (h - h0) / h0 )
@NLexpression(rocket, drag[j = 1:n], D_c * (v[j]^2) * exp(-h_c * (h[j] - h_0) / h_0))
# Grav(h)   = go * (h0 / h)^2
@NLexpression(rocket, grav[j = 1:n], g_0 * (h_0 / h[j])^2)
# Time of flight
@NLexpression(rocket, t_f, Δt * n)

# Dynamics
for j in 2:n
    # h' = v

    # Rectangular integration
    # @NLconstraint(rocket, h[j] == h[j - 1] + Δt * v[j - 1])

    # Trapezoidal integration
    @NLconstraint(rocket,
        h[j] == h[j - 1] + 0.5 * Δt * (v[j] + v[j - 1]))

    # v' = (T-D(h,v))/m - g(h)

    # Rectangular integration
    # @NLconstraint(rocket, v[j] == v[j - 1] + Δt *(
    #                 (T[j - 1] - drag[j - 1]) / m[j - 1] - grav[j - 1]))

    # Trapezoidal integration
    @NLconstraint(rocket,
        v[j] == v[j-1] + 0.5 * Δt * (
            (T[j] - drag[j] - m[j] * grav[j]) / m[j] +
            (T[j - 1] - drag[j - 1] - m[j - 1] * grav[j - 1]) / m[j - 1]))

    # m' = -T/c

    # Rectangular integration
    # @NLconstraint(rocket, m[j] == m[j - 1] - Δt * T[j - 1] / c)

    # Trapezoidal integration
    @NLconstraint(rocket,
        m[j] == m[j - 1] - 0.5 * Δt * (T[j] + T[j-1]) / c)
end

# Solve for the control and state
println("Solving...")
status = optimize!(rocket)

# Display results
# println("Solver status: ", status)
println("Max height: ", objective_value(rocket))

# Can visualize the state and control variables
using Gadfly

h_plot = plot(x = (1:n) * value.(Δt), y = value.(h)[:], Geom.line,
                Guide.xlabel("Time (s)"), Guide.ylabel("Altitude"))
m_plot = plot(x = (1:n) * value.(Δt), y = value.(m)[:], Geom.line,
                Guide.xlabel("Time (s)"), Guide.ylabel("Mass"))
v_plot = plot(x = (1:n) * value.(Δt), y = value.(v)[:], Geom.line,
                Guide.xlabel("Time (s)"), Guide.ylabel("Velocity"))
T_plot = plot(x = (1:n) * value.(Δt), y = value.(T)[:], Geom.line,
                Guide.xlabel("Time (s)"), Guide.ylabel("Thrust"))
draw(SVG(6inch, 6inch), vstack(hstack(h_plot, m_plot), hstack(v_plot, T_plot)))


text:
function helper_fn(z, attempt)
	single_data_name, test_fold = z
	proc_id = myid();

	helper_start_time = time()

	transaction_counts = convert(Array{Int64}, readcsv( string("../tcm_data/IRI_data_assortments/", single_data_name, "_transaction_counts.csv") ) )
	assortments = convert(Array{Int64}, readcsv( string("../tcm_data/IRI_data_assortments/", single_data_name, "_assortments.csv") ) )
	folds  = convert(Array{Int64}, readcsv( string("../tcm_data/IRI_data_assortments/", single_data_name, "_folds.csv") ) )

	train_assortments = assortments[ find(folds .!= test_fold), :]
	test_assortments = assortments[ find(folds .== test_fold), :]

	train_transaction_counts = transaction_counts[ find(folds .!= test_fold), :]
	test_transaction_counts = transaction_counts[ find(folds .== test_fold), :]

	num_train_transactions = sum(train_transaction_counts)
	num_test_transactions = sum(test_transaction_counts)

	num_train_assortments = size(train_assortments, 1)
	num_test_assortments = size(test_assortments, 1)

	N = 10

	time_limit = 1200

	perTol = 1e-6

	## RANKING-BASED MODEL ##
	CS_size_vec = [2,3,4,5,6,7,8,9]

	for CS_size in CS_size_vec

		ranking_lambda, orderings, ranking_loglik, ranking_elapsed_time = tcm_estimateRanking_CS(N, train_assortments, train_transaction_counts, time_limit, perTol, CS_size)
		ranking_K = length(ranking_lambda)

		test_ranking_predict, test_ranking_predict_by_class = tcm_predictRanking(N, ranking_K, ranking_lambda, orderings, test_assortments)
		ranking_KL, ranking_noinf_count, ranking_KL_avg = tcm_evaluateKL(N, test_ranking_predict, test_transaction_counts)

		avg_CS_size, weighted_avg_CS_size, max_CS_size = tcm_rankingStats(N, ranking_lambda, orderings)

		row_string = string(single_data_name,  ",",
								test_fold, ",",
								num_train_assortments, ",",
								num_test_assortments, ",",
								num_train_transactions, ",",
								num_test_transactions, ",",
								proc_id, ",",
								attempt, ",",
								"ranking", ",",
								CS_size, ",",
								ranking_loglik, ",",
								ranking_K, ",",
								perTol, ",",
								ranking_KL, ",",
								ranking_noinf_count, ",",
								ranking_KL_avg, ",",
								ranking_elapsed_time, ",",
								avg_CS_size, ",",
								weighted_avg_CS_size, ",",
								max_CS_size)
		row_string = string(row_string, "\n");

		@spawnat 1 begin
			print(ranking_outcsvhandle, row_string);
			flush(ranking_outcsvhandle);
		end
	end

	## RANKING-CV MODEL ##

	numInitializations = 5
	perTol_CV = 1e-6;
	CS_list = [2,3,4,5,6,7,8,9]
	nfolds = 4
	train_folds = folds[ find(folds .!= test_fold)]
	best_CS_size, best_KL_cv, cv_elapsed_time = tcm_estimateRanking_kcv_weighted_ass(N, CS_list, nfolds, train_folds, train_assortments, train_transaction_counts, numInitializations, time_limit, perTol_CV)

	row_string = string(single_data_name,  ",",
							test_fold, ",",
							num_train_assortments, ",",
							num_test_assortments, ",",
							num_train_transactions, ",",
							num_test_transactions, ",",
							proc_id, ",",
							attempt, ",",
							"ranking-CV", ",",
							best_CS_size, ",",
							perTol_CV, ",",
							numInitializations, ",",
							best_KL_cv, ",",
							cv_elapsed_time)
	row_string = string(row_string, "\n");

	@spawnat 1 begin
		print(ranking_cv_outcsvhandle, row_string);
		flush(ranking_cv_outcsvhandle);
	end

	@show best_CS_size

	helper_end_time = time()

	return helper_end_time - helper_start_time;

end

ext: jl

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_licenses: MIT Unlicense


text:
import BLASBenchmarksCPU
import StatsPlots
@testset "Interface" begin
    benchmark_result = BLASBenchmarksCPU.runbench(Float64; sizes = [1, 2, 5, 10, 20, 50, 100, 200], threaded=false) #test that threads=false at least doesn't throw somewhere.
    dfmin = BLASBenchmarksCPU.benchmark_result_df(benchmark_result) # minimum
    dfmedian = BLASBenchmarksCPU.benchmark_result_df(benchmark_result, :median)
    dfmean = BLASBenchmarksCPU.benchmark_result_df(benchmark_result, :mean)
    dfmax = BLASBenchmarksCPU.benchmark_result_df(benchmark_result, :maximum)
    @test_throws ArgumentError  BLASBenchmarksCPU.benchmark_result_df(benchmark_result, :foobar)
    @test dfmin isa BLASBenchmarksCPU.DataFrame
    @test dfmedian isa BLASBenchmarksCPU.DataFrame
    @test dfmean isa BLASBenchmarksCPU.DataFrame
    @test dfmax isa BLASBenchmarksCPU.DataFrame
    for df ∈ (dfmin,dfmedian,dfmean,dfmax)
        df[!, :Size] = Float64.(df[!, :Size]);
        df[!, :GFLOPS] = Float64.(df[!, :GFLOPS]);
        df[!, :Seconds] = Float64.(df[!, :Seconds]);
        p = StatsPlots.@df df StatsPlots.plot(:Size, :GFLOPS; group = :Library, legend = :bottomright)
        @test p isa StatsPlots.Plots.Plot
    end
    @test all(dfmin[!, :GFLOPS] .≥ dfmedian[!, :GFLOPS])
    @test all(dfmin[!, :GFLOPS] .≥ dfmean[!, :GFLOPS])
    @test all(dfmin[!, :GFLOPS] .≥ dfmax[!, :GFLOPS])
    @test any(dfmin[!, :GFLOPS] .≠ dfmedian[!, :GFLOPS])
    @test any(dfmin[!, :GFLOPS] .≠ dfmean[!, :GFLOPS])
    @test any(dfmin[!, :GFLOPS] .≠ dfmax[!, :GFLOPS])
    @test any(dfmedian[!, :GFLOPS] .≥ dfmax[!, :GFLOPS])
    @test any(dfmean[!, :GFLOPS] .≥ dfmax[!, :GFLOPS])
    @test any(dfmedian[!, :GFLOPS] .≠ dfmax[!, :GFLOPS])
    @test any(dfmean[!, :GFLOPS] .≠ dfmax[!, :GFLOPS])
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/HeatTransfer.jl/blob/master/LICENSE

using Pkg, Documenter, Literate, Dates, JuliaFEMDocs

PAGES = Any[]

root = abspath(pathof(JuliaFEMDocs), "..", "..", "docs", "src")

# First we get main page from JuliaFEM.jl/README.md

using JuliaFEM
JuliaFEMDocs.copydocs(JuliaFEM)
cp(joinpath(root, "packages", "JuliaFEM", "README.md"), joinpath(root, "index.md"); force=true)
push!(PAGES, "Home" => "index.md")

# Add packages

PACKAGES = Any[]

# HeatTransfer.jl
using HeatTransfer
JuliaFEMDocs.copydocs(HeatTransfer)
HeatTransferDocs = Any[]
push!(HeatTransferDocs, "Home" => "packages/HeatTransfer/README.md")
push!(HeatTransferDocs, "Docs" => "packages/HeatTransfer/index.md")
push!(PACKAGES, "HeatTransfer.jl" => HeatTransferDocs)

push!(PAGES, "Packages" => PACKAGES)

# Done. Make docs.

makedocs(modules=[JuliaFEMDocs],
         format = Documenter.HTML(),
         checkdocs = :all,
         sitename = "JuliaFEM.jl",
         authors = "Jukka Aho",
         pages = PAGES)


text:
import ContinuumArrays: apply, MulQuasiMatrix

@testset "$(rpad("Bernstein Basis Tests",80))" begin

    l = Bernstein(2)
    d = Derivative(axes(l,1))

    @test apply(*,d,l) isa BernsteinDerivative
    @test d*l isa BernsteinDerivative
    @test l' isa BernsteinDerivative

    @test basis(l) == l.b
    @test nbasis(l) == 2
    @test eachindex(l) == 0:1
    @test order(l) == 2
    @test degree(l) == 1

    l1 = Bernstein(Float64, 2)
    l2 = Bernstein(Integer, 2)
    l3 = Bernstein(Float64, 3)

    @test hash(l) == hash(l1)
    @test hash(l) == hash(l2)
    @test hash(l) != hash(l3)

    @test l == l1
    @test l == l2
    @test l != l3

    @test  isequal(l, l1)
    @test !isequal(l, l2)
    @test !isequal(l, l3)

    y = rand(5)

    z1 = parent([ l[y[i], j] for i in eachindex(y), j in eachindex(l)])
    z2 =   hcat([ l[y[i], :] for i in eachindex(y)]...)'
    z3 =   hcat([ l[y,    j] for j in eachindex(l)]...)
    z4 = parent(  l[y,    :] )

    @test z1 == z2 == z3 == z4

    z1 = parent([ (d*l)[y[i], j] for i in eachindex(y), j in eachindex(l)])
    z2 =   hcat([ (d*l)[y[i], :] for i in eachindex(y)]...)'
    z3 =   hcat([ (d*l)[y,    j] for j in eachindex(l)]...)
    z4 = parent(  (d*l)[y,    :] )

    @test z1 == z2 == z3 == z4

    @test l(0.0, 0) == 1.0
    @test l(0.5, 0) == 0.5
    @test l(1.0, 0) == 0.0

    @test l(0.0, 1) == 0.0
    @test l(0.5, 1) == 0.5
    @test l(1.0, 1) == 1.0

    @test l[0.0, 0] == 1.0
    @test l[0.5, 0] == 0.5
    @test l[1.0, 0] == 0.0

    @test l[0.0, 1] == 0.0
    @test l[0.5, 1] == 0.5
    @test l[1.0, 1] == 1.0

    @test l[0, 0] == +1.0
    @test l[1, 0] ==  0.0
    @test l[2, 0] == -1.0

    @test l[0, 1] == 0.0
    @test l[1, 1] == 1.0
    @test l[2, 1] == 2.0

    @test (d*l)[0.0, 0] == -1.0
    @test (d*l)[0.5, 0] == -1.0
    @test (d*l)[1.0, 0] == -1.0

    @test (d*l)[0.0, 1] == +1.0
    @test (d*l)[0.5, 1] == +1.0
    @test (d*l)[1.0, 1] == +1.0

    @test (d*l)[0, 0] == -1.0
    @test (d*l)[1, 0] == -1.0
    @test (d*l)[2, 0] == -1.0

    @test (d*l)[0, 1] == +1.0
    @test (d*l)[1, 1] == +1.0
    @test (d*l)[2, 1] == +1.0

end


text:
function is_joinpathable(x)
    if x isa CSTParser.EXPR{CSTParser.Call} && x.args[1] isa CSTParser.IDENTIFIER && CSTParser.str_value(x.args[1]) == "joinpath"
        joinpathable = true
        for i = 2:length(x.args)
            arg = x.args[i]
            if !(arg isa CSTParser.PUNCTUATION)
                joinpathable &= CSTParser.is_lit_string(arg)
            end
        end
        return joinpathable
    end
    return false
end

function get_path(x::CSTParser.EXPR{CSTParser.Call})
    if length(x.args) == 4
        parg = x.args[3]
        if CSTParser.is_lit_string(parg)
            path = CSTParser.str_value(parg)
            return path
        elseif parg isa CSTParser.EXPR{CSTParser.Call} && parg.args[1] isa CSTParser.IDENTIFIER && CSTParser.str_value(parg.args[1]) == "joinpath"
            path = ""
            for i = 2:length(parg.args)
                arg = parg.args[i]
                if arg isa CSTParser.PUNCTUATION
                elseif CSTParser.is_lit_string(arg)
                    path = string(path, "/", CSTParser.str_value(arg))
                else
                    return ""
                end
            end
            return path
        end
    end
    return ""
end

function follow_include(x, s, S::State{FileSystem})
    path = get_path(x)
    path = isabspath(path) ? path : joinpath(dirname(S.loc.path), path)
    !isfile(path) && return

    parent = S.includes[S.loc.path]
    f = File(path, (parent, S.loc.offset + x.span), [])
    push!(parent.children, f)
    S.includes[path] = f

    x1 = CSTParser.parse(read(path, String), true)
    old_Sloc = S.loc
    S.loc = Location(path, 0)
    trav(x1, s, S)
    S.loc = old_Sloc
end


text:
###############################################################################
#
# Steps of non-steady state stomatal controls
#
###############################################################################
#=
0. Initialize plant and environment with minimal stomatal opening
1. Update transpiration rate
2. Update canopy temperature
3. Update plant hydraulic profile every 10 minutes
4. Calculate assimilation rate
5. Calculate dH₂O/dt and dCO₂/dt for the environment
6. Calculate dGs/dt for the given environment
7. Update stomatal conductance
8. Goto step 1

This function will be splited into many functions later!
=#

function test_diurnal(
            node::SPACMono{FT},
            Δt::FT=FT(10)
) where {FT<:AbstractFloat}
    # 0.1 create variables required
    @unpack n_canopy, wl_set = node;
    canopy_rt = node.canopy_rt;
    fraction_sl::Array{FT,1} = repeat(canopy_rt.lidf, outer=[ length(canopy_rt.lazitab) ]) / length(canopy_rt.lazitab);
    n_sl = length(canopy_rt.lidf) * length(canopy_rt.lazitab);

    # change node.in_rad to alter direct and diffuse fraction
    # sum of direct  energy in W:
    # e_dir = sum( (wl_set.swl[2:end] .- wl_set.swl[1:end-1]) .* node.in_rad.E_direct  ) / 1000
    # sum of diffuse energy in W:
    # e_dif = sum( (wl_set.swl[2:end] .- wl_set.swl[1:end-1]) .* node.in_rad.E_diffuse ) / 1000
    # target e_dir is tar_dir, do it this way:
    # node.in_rad.E_direct  .*= tar_dir / e_dir
    # target e_dif is tar_dif, do it this way:
    # node.in_rad.E_diffuse .*= tar_dif / e_dif

    # 0.2 A diurnal cycle for radiation and Tair
    Tmean  = FT(295.15)
    DeltaT = 3
    omega  = 2 * π / (24*3600)
    t      = range(0, stop=24*3600*2, step=Δt)
    phi_t  = omega .* t .- π .* ones(FT, size(t)) ./ 2
    PARr_t = zeros(FT, size(t))
    Tair_t = zeros(FT, size(t))
    Tlef_t = zeros(FT, size(t))
    N      = length(PARr_t)

    # 0.3 customize dirnal cycles of temperatures and PAR
    for i_tim in 1:N
        PARr_t[i_tim] = max( sin(phi_t[i_tim]), FT(0.0) )
        Tair_t[i_tim] = Tmean + DeltaT * sin(phi_t[i_tim] - π/3)
        Tlef_t[i_tim] = Tmean + DeltaT * sin(phi_t[i_tim] - π/3)
    end

    # 0.4 initialize the canopy temperature
    for i_can in 1:n_canopy
        # Assume RH = 50%
        photosystem_temperature_dependence!(node.plant_ps[i_can].ps.PSM, node.envirs[i_can], Tlef_t[1]);
        node.plant_ps[i_can].ps.p_H₂O_sat = saturation_vapor_pressure(Tlef_t[1]);
        node.envirs[i_can].p_H₂O = node.envirs[i_can].p_H₂O_sat / 2;
        # latent heat flux as well
        node.plant_ps[i_can].LV = latent_heat_vapor(Tlef_t[i_can]) * 1000 / 18;
        # vc dependecy as well
        temperature_effects!(node.plant_hs.leaves[i_can], Tlef_t[i_can]);
    end

    # run dirnal cycle
    # NPP per ground area
    Agpp_t = zeros(FT, size(t));
    # NPP per ground area
    Asum_t = zeros(FT, size(t));
    # trasnpiration rate per ground area
    Esum_t = zeros(FT, size(t));
    for i_tim in 1:N
        # 0. recalculate the canopy RT stuff
        zenith = zenith_angle(node.latitude, 180 + t[i_tim]/86400 );
        zenith = min(88, zenith);
        node.angles.tts = zenith;
        canopy_geometry!(canopy_rt, node.angles, node.can_opt, node.rt_con);
        canopy_matrices!(node.leaves_rt, node.can_opt);
        short_wave!(canopy_rt, node.can_opt, node.can_rad, node.in_rad, node.soil_opt, node.rt_con);
        canopy_fluxes!(canopy_rt, node.can_opt, node.can_rad, node.in_rad, node.soil_opt, node.leaves_rt, wl_set, node.rt_con);
        thermal_fluxes!(node.leaves_rt, node.can_opt, node.can_rad, canopy_rt, node.soil_opt, [FT(400.0)], wl_set);
        SIF_fluxes!(node.leaves_rt, node.can_opt, node.can_rad, canopy_rt, node.soil_opt, wl_set, node.rt_con, node.rt_dim);

        # 2. Update transpiration rate
        a_gpp = FT(0);
        a_sum = FT(0);
        e_sum = FT(0);
        for i_can in 1:n_canopy
            rt_layer = node.n_canopy + 1 - i_can;

            # set the diffuse par to all the leaves
            #node.plant_ps[i_can].APAR           .= PARr_t[i_tim]  * node.can_rad.absPAR_shadeCab[rt_layer] * FT(1e6);
            #node.plant_ps[i_can].APAR[1:end-1] .+= PARr_t[i_tim] .* reshape(node.can_rad.absPAR_sunCab[:,:,rt_layer],(:,1))[:,1] .* FT(1e6);

            # calculate the fraction of sunlit and shaded leaves
            f_view = (node.can_opt.Ps[rt_layer] + node.can_opt.Ps[rt_layer+1]) / 2;
            #la_new = [f_view .* fraction_sl; 1-f_view];
            #node.plant_ps[i_can].LAIx .= la_new;
            ##=
            for i_leaf in 1:n_sl
                node.plant_ps[i_can].APAR[i_leaf] = PARr_t[i_tim] * node.can_rad.absPAR_sunCab[(rt_layer-1)*n_sl + i_leaf] * FT(1e6);
                node.plant_ps[i_can].LAIx[i_leaf] = f_view * fraction_sl[i_leaf];
            end
            node.plant_ps[i_can].APAR[end] = PARr_t[i_tim] * node.can_rad.absPAR_shadeCab[rt_layer] * FT(1e6);
            node.plant_ps[i_can].LAIx[end] = 1 - f_view;
            # =#

            # update leaf temperature
            # update the parameters from stomatal model
            node.plant_ps[i_can].T = Tlef_t[i_tim];
            update_leaf_TP!(node.plant_ps[i_can], node.plant_hs.leaves[i_can], node.envirs[i_can]);
            # uncomment this for Sperry, Eller, and Gentine models
            #update_leaf_AK!(node.plant_ps[i_can], node.plant_hs.leaves[i_can], node.envirs[i_can]);

            #canopyi.t_list[1:end-1] .= reshape(node.can_rad.T_sun3D[:,:,rt_layer],(:,1))[:,1]
            #canopyi.t_list[end]      = node.can_rad.T_shade[rt_layer]
            # Assume RH = 50%
            node.envirs[i_can].t         = Tair_t[i_tim];
            node.envirs[i_can].p_H₂O_sat = saturation_vapor_pressure(Tair_t[i_tim]);
            node.envirs[i_can].p_H₂O     = node.envirs[i_can].p_H₂O_sat / 2;

            # update the flow rates
            g_i_can = FT(0);
            a_i_can = FT(0);
            e_i_can = FT(0);
            for i_leaf in 1:(n_sl+1)
                g_i_can += node.plant_ps[i_can].Ag[i_leaf] * node.plant_ps[i_can].LAIx[i_leaf] * node.plant_ps[i_can].LA;
                a_i_can += node.plant_ps[i_can].An[i_leaf] * node.plant_ps[i_can].LAIx[i_leaf] * node.plant_ps[i_can].LA;
                e_i_can += node.plant_ps[i_can].g_lw[i_leaf] * (node.plant_ps[i_can].p_sat - node.envirs[i_can].p_H₂O) / node.envirs[i_can].P_AIR * node.plant_ps[i_can].LAIx[i_leaf] *
                           node.plant_ps[i_can].LA;
            end
            #g_i_can = sum( node.plant_ps[i_can].Ag .* node.plant_ps[i_can].LAIx ) * node.plant_ps[i_can].LA;
            #a_i_can = sum( node.plant_ps[i_can].An .* node.plant_ps[i_can].LAIx ) * node.plant_ps[i_can].LA;
            #e_i_can = sum( node.plant_ps[i_can].g_lw .* (node.plant_ps[i_can].p_sat - node.envirs[i_can].p_H₂O) ./ node.envirs[i_can].P_AIR .* node.plant_ps[i_can].LAIx ) * node.plant_ps[i_can].LA;

            # calculate the photosynthetic rates
            gas_exchange!(node.plant_ps[i_can], node.envirs[i_can], GswDrive());
            gsw_control!(node.plant_ps[i_can], node.envirs[i_can]);

            # use the ball-berry model here for now as the ∂A/∂E and ∂A/∂Θ functions are not yet ready
            gsw_ss = stomatal_conductance(node.stomata_model, node.plant_ps[i_can], node.envirs[i_can], FT(1));

            # assume τ = 10 minutes
            for i_leaf in 1:(n_sl+1)
                node.plant_ps[i_can].g_sw[i_leaf] += (gsw_ss[i_leaf] - node.plant_ps[i_can].g_sw[i_leaf]) / 600 * Δt;
            end
            #node.plant_ps[i_can].g_sw .+= (gsw_ss .- node.plant_ps[i_can].g_sw) / 600 * Δt;

            a_gpp += g_i_can;
            a_sum += a_i_can;
            e_sum += e_i_can;
        end

        Agpp_t[i_tim] = a_gpp / node.ga;
        Asum_t[i_tim] = a_sum / node.ga;
        Esum_t[i_tim] = e_sum / node.ga;

        # update pressure profile every 10 minutes
        # TO BE DONE within PlantHydraulics

        _par = PARr_t[i_tim]  * node.can_rad.absPAR_shadeCab[1] * FT(1e6);
    end

    #=
    # plot the results
    figure(1)
    clf()
    plot(t ./ 3600, Agpp_t, "g-");
    plot(t ./ 3600, Asum_t, "k-");

    figure(2)
    clf()
    plot(t ./ 3600, Esum_t, "k-");
    =#

    return nothing
end


text:
export Line, Wire, NWN, WireProp, NWN_JDA, len

using StaticArrays, LinearAlgebra, LightGraphs, SimpleWeightedGraphs

"""
A structure representing a line with end points
p₁ and p₂. Parameterized by the datatype of the
coordinates and the length (i.e. dimension) N.
"""
struct Line{T<:Real,N}
    p₁::SVector{N,T}
    p₂::SVector{N,T}
end

struct WireProp{T<:Number}
    ρ::T # Resistivity
    D::T # Wire diameter
end

mutable struct NWN{S<:Real,T<:Number,N,M,P}
    lines::Vector{Line{S,N}}
    props::Vector{WireProp{T}}
    elecs::Vector{Line{S,N}}
    volts::Vector{T}
    grnds::Vector{Line{S,N}}
    dims::SVector{N,S}
    Rⱼ::T
    Rₑ::T
    graph::SimpleWeightedGraph{Int, T}
    function NWN{S,T,N,M,P}(
        lines::Vector{Line{S,N}},props::Vector{WireProp{T}},
        elecs::Vector{Line{S,N}},volts::Vector{T},grnds::Vector{Line{S,N}},
        dims::SVector{N,S},Rⱼ::T,Rₑ::T
        ) where {S,T,N,M,P}

        nwn = new{S,T,N,M,P}(lines, props, elecs, volts, grnds, dims, Rⱼ, Rₑ)
        nwn.graph = graph(nwn)
        nwn
    end
end

function graph(nwn::NWN{S,T,N,:JDA,P}) where {S,T,N,P}
    ww_src, ww_dst = connections(nwn.lines, nwn.dims)
    we_src, we_dst = connections(nwn.lines, nwn.elecs, nwn.dims)
    wg_src, wg_dst = connections(nwn.lines, nwn.grnds, nwn.dims)
    nₗ, nₑ = length(nwn.lines), length(nwn.elecs)
    mₗ, mₑ, m₀ = length(ww_src), length(we_src), length(wg_src)
    src = vcat(ww_src, we_src, wg_src)
    dst = vcat(ww_dst, we_dst.+nₗ, wg_dst.+(nₗ+nₑ))
    wgt = vcat(ones(mₗ)./nwn.Rⱼ,ones(mₑ+m₀)./nwn.Rₑ)
    g = SimpleWeightedGraph(src, dst, wgt)
    g
end


text:
##TODO: the question: do I follow the tree strucutre of matplotlib, or try to emulate
# the more logical structure of SVG? (meaning I would need to do mappings like
# facecolor -> fill, edgecolor -> stroke, linewidth -> stroke_width etc)

# axes style
type AxesStyle
    stroke_width::Real # linewidth
    stroke::RGB # edgecolor
    #grid
    #facecolor::Color
    #labelcolor::Color
    #axisbelow
    function AxesStyle(;stroke_width, stroke)
        AxesStyle(stroke_width, stroke)
    end
end

# grid style

# legend style

# tick style

# Container for all styles
type Style
    axes::AxesStyle
end

function style(;axes)
end

function style!(s::Style)

end

max_stars_repo_name: gabrielgellner/OdinSon

max_issues_repo_name: gabrielgellner/OdinSon

max_forks_repo_name: gabrielgellner/OdinSon


text:
# these do not require decoding to work correctly

typemin(::Type{Float65}) = typemin(Float64)
typemax(::Type{Float65}) = typemin(Float64)
typemin(x::Float65) = typemin(Float65)
typemax(x::Float65) = typemax(Float65)

realmin(::Type{Float65}) = nextfloat(Tiny(Float64))
realmax(::Type{Float65}) = prevfloat(Huge(Float64))
realmin(x::Float65) = realmin(Float65)
realmax(x::Float65) = realmax(Float65)

zero(::Float65)  = convert(Float65, 0.0)
one(::Float65)   = convert(Float65, 1.0)
zero(x::Float65) = zero(Float65)
one(x::Float65)  = one(Float65)

sizeof(::Type{Float65}) = 8
sizeof(x::Float65) = sizeof(Float65)

hash(x::Float65, h::UInt) = hash(x.fp, h)


text:
# runtests.jl

using Test: @test, @testset, detect_ambiguities
using Sinograms

@testset "Sinograms" begin
    # todo
    include("fbp-par.jl")
    include("bdd_2d.jl")

    @test length(detect_ambiguities(Sinograms)) == 0
end


text:
using BNPRegressionGGA2021
using CSV
using DataStructures
using DataFrames
using Random
const BNP = BNPRegressionGGA2021

# Get the simulated data
data = CSV.read("data/simulated-data-erlang-nph.csv", DataFrame)

# Get the best γ for 100 simulated samples (using the MAP estimator)
begin
    Niter = 100
    Random.seed!(1);
    df = zeros(Bool, Niter, 6)
    gdata = groupby(data, :iter)
    for iter in 1:Niter
        c0 = gdata[iter][!, r"c"] |> Matrix |> vec
        y0 = gdata[iter][!, r"y"] |> Matrix |> vec
        X0 = gdata[iter][!, r"x"] |> Matrix
        y1 = zeros(0)
        X1 = zeros(0, size(X0, 2))
        m = BNP.DGPMErlang(; c0, y0, X0, y1, X1);
        _, _, _, chaing = BNP.sample!(m; mcmcsize = 10000, burnin = 5000);
        df[iter, :] .= begin
            chaing |>
            DataStructures.counter |>
            collect |>
            x -> sort(x, by = x -> x[2], rev = true) |>
            x -> getindex(x, 1) |>
            x -> getindex(x, 1)
        end
        println(iter)
    end
    filename = "data/simulated-best-gamma-erlang-nph-bnp.csv"
    CSV.write(filename, DataFrame(df, :auto))
end


text:
# example for README

using DiffWrappers
using DiffResults
f(x) = sum(x .* linspace(0, 1, length(x)));
g = ForwardGradientWrapper(f, 4); # will use length 4 vectors
length(g)
x = ones(4);
gx = g(x);
DiffResults.value(gx)
DiffResults.gradient(gx)


text:
# This is a `mutable struct` ATM so that it can be replaced with a 64 bit CAS.
# TODO: make it an immutable `struct` and compute `log2length` on-the-fly?
mutable struct CircularVector{T} <: AbstractVector{T}
    log2length::Int
    data::Vector{T}
end

function CircularVector{T}(log2length::Int) where {T}
    @assert log2length >= 0
    data = Vector{T}(undef, 1 << log2length)
    Nothing <: T && fill!(data, nothing)  # [^speculative_load]
    T <: Signed && fill!(data, typemin(T))
    return CircularVector{T}(log2length, data)
end

Base.size(A::CircularVector{T}) where {T} = size(A.data)

@noinline function _check_data_length(A::CircularVector)
    @assert 1 << A.log2length == length(A.data)
end

Base.@propagate_inbounds function indexof(A::CircularVector, i::Int)
    @boundscheck _check_data_length(A)
    return (i - 1) & ((1 << A.log2length) - 1) + 1
end

Base.@propagate_inbounds Base.getindex(A::CircularVector, i::Int) = A.data[indexof(A, i)]

Base.@propagate_inbounds function Base.setindex!(A::CircularVector, v, i::Int)
    v = convert(eltype(A), v)
    A.data[indexof(A, i)] = v
end

Base.pointer(A::CircularVector, i::Integer) = pointer(A.data, indexof(A, i))

function tryresize(A::CircularVector, log2inc::Integer, indices)
    log2length = A.log2length + log2inc
    n = 1 << log2length
    n < length(indices) && return nothing
    B = CircularVector{eltype(A)}(log2length)
    for i in indices
        @inbounds B[i] = A[i]
    end
    return B
end

mutable struct WorkStealingDeque{T,S}
    @atomic buffer::CircularVector{S}
    @atomic top::Int
    @atomic bottom::Int
    # TODO: pad
end

function WorkStealingDeque{T}() where {T}
    if is_pointerfree_type(T)
        S = T
    else
        # Any element type that may contain a boxed object uses
        # `CircularVector{Any}` which initialize all locations with `nothing`.
        # This is for supporting "blind" `buffer[top]`.  [^speculative_load]
        S = Any
    end
    return WorkStealingDeque{T,S}(CircularVector{S}(4), 1, 1)
end

Base.eltype(::Type{WorkStealingDeque{T}}) where {T} = T

function Base.length(deque::WorkStealingDeque)
    bottom = @atomic deque.bottom
    top = @atomic deque.top
    return bottom - top
end

function try_resize_buffer!(
    deque::WorkStealingDeque,
    buffer::CircularVector,
    top::Int,
    bottom::Int,
    log2inc::Integer,
)
    target_size = 1 << (buffer.log2length + log2inc)
    current_size = bottom - top
    if target_size <= max(current_size, 4)
        # too large to shrink
        return nothing
    end

    buffer = tryresize(buffer, log2inc, top:bottom-1)
    if buffer === nothing
        return nothing
    else
        @atomic deque.buffer = buffer
        return buffer
    end
end

grow!(deque::WorkStealingDeque, buffer, top, bottom) =
    something(try_resize_buffer!(deque, buffer, top, bottom, 1))
tryshrink!(deque::WorkStealingDeque, buffer, top, bottom) =
    try_resize_buffer!(deque, buffer, top, bottom, -1)

function Base.push!(deque::WorkStealingDeque, v)
    v = convert(eltype(deque), v)
    bottom = @atomic deque.bottom
    top = @atomic deque.top
    buffer = @atomic deque.buffer
    current_size = bottom - top
    if current_size >= length(buffer) - 1
        buffer = grow!(deque, buffer, top, bottom)
        @atomic deque.buffer = buffer
    end
    # TODO: Technically, this should use atomic store. However, there is no way
    # to do this in a GC-compatible way at the moment.
    buffer[bottom] = v
    bottom += 1
    @atomic deque.bottom = bottom
    return deque
end

function ConcurrentCollections.maybepop!(deque::WorkStealingDeque)
    bottom = @atomic deque.bottom
    buffer = @atomic deque.buffer
    bottom -= 1
    @atomic deque.bottom = bottom
    top = @atomic deque.top
    next_size = bottom - top
    if next_size < 0
        @atomic deque.bottom = top
        return nothing
    end
    # TODO: Technically, this should also use atomic load
    r = Some(buffer[bottom])
    if next_size > 0
        tryshrink!(deque, buffer, top, bottom)
        return r
    end
    bottom = top + 1
    if !@atomicreplace(deque.top, top => top + 1)[2]
        r = nothing
    end
    @atomic deque.bottom = bottom
    return r
end

function ConcurrentCollections.maybepopfirst!(deque::WorkStealingDeque{T}) where {T}
    top = @atomic deque.top
    bottom = @atomic deque.bottom
    buffer = @atomic deque.buffer
    current_size = bottom - top
    if current_size <= 0
        return nothing
    end
    # TODO: Technically, this should be an atomic load.  See below for some
    # discussions.  [^speculative_load]
    y = buffer[top]
    if @atomicreplace(deque.top, top => top + 1)[2]
        y = y::T
        return Some{T}(y)
    else
        return nothing
    end
end
# [^speculative_load]: The difficulty here is that we don't know the validity of
# `y = buffer[top]` until the CAS below but we still need to make sure `y` is a
# ("somewhat") valid Julia object since it has to be rooted before the CAS.
# This is because successful CAS means that `buffer` may not root `y` anymore.
# This can be an issue if the compiler inserts a safepoint between `buffer[top]`
# and the CAS.  The current implementation relies on the stop-the-world behavior
# of the GC (the "somewhat" part).  That is to say, the object `y` is not
# *always* properly accessible at the (potential) safepoint before the CAS.
# However, since the GC (currently) starts traversing the object graph after it
# knows all worker threads reach the safepoint, the implementation above may
# actually be OK for now.
#
# Another issue is that `buffer[top]` may not be assigned and accessing it can
# throw if the element type is a boxed value.  The current workaround is to fill
# `buffer` with `nothing`. Alternatively, checking the buffer with `isassigned`
# should also work. However, since dequeue can happen many times while the
# buffer is initialized less frequently (upon construction and resizes),
# the `isassigned`-based strategy is not used.
#
# TODO: The current implementation also expects that it is "OK" to load a
# pointerfree value non-atomically (it's a data race in C++20 memory model hence
# an UB). A proper solution may be to store/load these values as a sequence of
# `UInt`s with relaxed atomics and then type-pun `NTuple{_,UInt}` into the Julia
# values.  See also the discussion in the C++ standard: [p0690r1: Tearable
# Atomics](http://www.open-std.org/jtc1/sc22/wg21/docs/papers/2018/p0690r1.html)

function Base.pop!(deque::WorkStealingDeque)
    r = maybepop!(deque)
    if r === nothing
        error("deque is empty")
    else
        return something(r)
    end
end

function Base.popfirst!(deque::WorkStealingDeque)
    r = maybepopfirst!(deque)
    if r === nothing
        error("deque is empty")
    else
        return something(r)
    end
end

function Base.sizehint!(deque::WorkStealingDeque, n::Integer)
    bottom = @atomic deque.bottom
    top = @atomic deque.top
    buffer = @atomic deque.buffer
    log2inc = ceillog2(n) - buffer.log2length
    if log2inc != 0
        try_resize_buffer!(deque, buffer, top, bottom, log2inc)
    end
    return deque
end

#=
struct WorkStealingDequeSet{T}
    deques::Vector{WorkStealingDeque{T}}
end

WorkStealingDequeSet{T}() where {T} =
    WorkStealingDequeSet{T}([WorkStealingDeque{T}() for _ in 1:Threads.nthreads()])
=#


text:
# OpenCL.Context

const _ctx_reference_count = Dict{CL_context, Int}()

function create_jl_reference!(ctx_id::CL_context)
    if haskey(_ctx_reference_count, ctx_id) # for the first jl reference, we already have a refcount of 1
        @check api.clRetainContext(ctx_id) # increase internal refcount, if creating an additional reference
    end
    refcount = get!(_ctx_reference_count, ctx_id, 0)
    _ctx_reference_count[ctx_id] = refcount + 1
    return
end
function free_jl_reference!(ctx_id::CL_context)
    if !haskey(_ctx_reference_count, ctx_id)
        error("Freeing unknown context")
    end
    refcount = _ctx_reference_count[ctx_id]
    if refcount == 0
        error("Double free of context id: ", ctx_id)
    elseif refcount == 1
        delete!(_ctx_reference_count, ctx_id)
        return
    end
    _ctx_reference_count[ctx_id] = refcount - 1
    return
end

type Context <: CLObject
    id :: CL_context
    # If created from ctx_id already, we need to increase the reference count
    # because then we give out multiple context references with multiple finalizers to the world
    # TODO should we make it in a way, that you can't overwrite it?
    function Context(ctx_id::CL_context; retain = false)
        retain && @check api.clRetainContext(ctx_id)
        if !is_ctx_id_alive(ctx_id)
            error("ctx_id not alive: ", ctx_id)
        end
        ctx = new(ctx_id)
        create_jl_reference!(ctx_id)
        finalizer(ctx, c -> begin
            retain || _deletecached!(c);
            if c.id != C_NULL
                release_ctx_id(c.id)
                free_jl_reference!(c.id)
                c.id = C_NULL
            end
        end )
        return ctx
    end
end

number_of_references(ctx::Context) = number_of_references(ctx.id)
function number_of_references(ctx_id::CL_context)
    refcounts = Ref{CL_uint}()
    @check api.clGetContextInfo(
        ctx_id, CL_CONTEXT_REFERENCE_COUNT,
        sizeof(CL_uint), refcounts, C_NULL
    )
    return refcounts[]
end

function is_ctx_id_alive(ctx_id::CL_context)
    number_of_references(ctx_id) > 0
end
function release_ctx_id(ctx_id::CL_context)
    if is_ctx_id_alive(ctx_id)
        @check api.clReleaseContext(ctx_id)
    else
        error("Double free for context: ", ctx_id)
    end
    return
end

Base.pointer(ctx::Context) = ctx.id

function Base.show(io::IO, ctx::Context)
    dev_strs = [replace(d[:name], r"\s+", " ") for d in devices(ctx)]
    devs_str = join(dev_strs, ",")
    ptr_val = convert(UInt, Base.pointer(ctx))
    ptr_address = "0x$(hex(ptr_val, Sys.WORD_SIZE>>2))"
    print(io, "OpenCL.Context(@$ptr_address on $devs_str)")
end

immutable _CtxErr
    handle :: Ptr{Void}
    err_info :: Ptr{Cchar}
    priv_info :: Ptr{Void}
    cb :: Csize_t
end

const io_lock = ReentrantLock()
function log_error(message...)
    @async begin
        lock(STDERR)
        lock(io_lock)
        print(STDERR, string(message..., "\n"))
        unlock(io_lock)
        unlock(STDERR)
    end
end

function ctx_notify_err(
        err_info::Ptr{Cchar}, priv_info::Ptr{Void},
        cb::Csize_t, func::Ptr{Void}
    )
    ccall(func, Void, (Ptr{Cchar}, Ptr{Void}, Csize_t), err_info, priv_info, cb)
    return
end

ctx_callback_ptr() = cfunction(ctx_notify_err, Void,
                               Tuple{Ptr{Cchar}, Ptr{Void}, Csize_t, Ptr{Void}})

function raise_context_error(err_info, private_info, cb)
    log_error("OpenCL Error: | ", unsafe_string(err_info), " |")
    return
end

function Context(devs::Vector{Device};
                         properties=nothing,
                         callback::Union{Function, Void} = nothing)
    if isempty(devs)
        ArgumentError("No devices specified for context")
    end
    if properties !== nothing
        ctx_properties = _parse_properties(properties)
    else
        ctx_properties = C_NULL
    end

    n_devices = length(devs)
    device_ids = Vector{CL_device_id}(n_devices)
    for (i, d) in enumerate(devs)
        device_ids[i] = d.id
    end

    err_code = Ref{CL_int}()
    payload = callback == nothing ? raise_context_error : callback
    f_ptr = cfunction(payload, Void, Tuple{Ptr{Cchar}, Ptr{Void}, Csize_t})
    ctx_id = api.clCreateContext(
        ctx_properties, n_devices, device_ids,
        ctx_callback_ptr(), f_ptr, err_code
    )
    if err_code[] != CL_SUCCESS
        throw(CLError(err_code[]))
    end
    return Context(ctx_id)
end

Context(d::Device; properties=nothing, callback=nothing) =
        Context([d], properties=properties, callback=callback)

function Context(dev_type::CL_device_type;
                 properties=nothing, callback=nothing)
    if properties !== nothing
        ctx_properties = _parse_properties(properties)
    else
        ctx_properties = C_NULL
    end
    if callback !== nothing
        ctx_user_data = callback
    else
        ctx_user_data = raise_context_error
    end
    err_code = Ref{CL_int}()
    ctx_id = api.clCreateContextFromType(ctx_properties, dev_type,
                                         ctx_callback_ptr(), ctx_user_data, err_code)
    if err_code[] != CL_SUCCESS
        throw(CLError(err_code[]))
    end
    return Context(ctx_id)
end

function Context(dev_type::Symbol;
                 properties=nothing, callback=nothing)
    Context(cl_device_type(dev_type),
            properties=properties, callback=callback)
end

function properties(ctx_id::CL_context)
    nbytes = Ref{Csize_t}(0)
    @check api.clGetContextInfo(ctx_id, CL_CONTEXT_PROPERTIES, 0, C_NULL, nbytes)

    # Calculate length of storage array
    # At nbytes[] the size of the properties array in bytes is stored
    # The length of the property array is then nbytes[] / sizeof(CL_context_properties)
    # Note: nprops should be odd since it requires a C_NULL terminated array
    nprops = div(nbytes[], sizeof(CL_context_properties))

    props = Vector{CL_context_properties}(nprops)
    @check api.clGetContextInfo(ctx_id, CL_CONTEXT_PROPERTIES,
                                nbytes[], props, C_NULL)
    #properties array of [key,value..., C_NULL]
    result = Any[]
    for i in 1:2:nprops
        key = props[i]
        value = i < nprops ? props[i+1] : nothing

        if key == CL_CONTEXT_PLATFORM
            push!(result, (key, Platform(cl_platform_id(value))))
        elseif key == CL_GL_CONTEXT_KHR ||
           key == CL_EGL_DISPLAY_KHR ||
           key == CL_GLX_DISPLAY_KHR ||
           key == CL_WGL_HDC_KHR ||
           key == CL_CGL_SHAREGROUP_KHR
            push!(result, (key, value))
        elseif is_apple() ? (key == CL_CONTEXT_PROPERTY_USE_CGL_SHAREGROUP_APPLE) : false
            push!(result, (key, value))
        elseif key == 0
            if i != nprops
                warn("Encountered OpenCL.Context property key == 0 at position $i")
            end
            break
        else
            warn("Unknown OpenCL.Context property key encountered $key")
        end
    end
    return result
end

function properties(ctx::Context)
    properties(ctx.id)
end

#Note: properties list needs to be terminated with a NULL value!
function _parse_properties(props)
    if isempty(props)
        return C_NULL
    end
    cl_props = CL_context_properties[]
    for prop_tuple in props
        if length(prop_tuple) != 2
            throw(ArgumentError("Context property tuples must be of type (key, value)"))
        end
        prop, val = prop_tuple
        push!(cl_props, cl_context_properties(prop))
        if prop == CL_CONTEXT_PLATFORM
            isa(val, Platform) && (val = val.id)
            push!(cl_props, cl_context_properties(val))
        elseif prop == CL_WGL_HDC_KHR
            push!(cl_props, cl_context_properties(val))
        elseif is_apple() ? (prop == CL_CONTEXT_PROPERTY_USE_CGL_SHAREGROUP_APPLE) : false
            push!(cl_props, cl_context_properties(val))
        elseif prop == CL_GL_CONTEXT_KHR ||
            prop == CL_EGL_DISPLAY_KHR ||
            prop == CL_GLX_DISPLAY_KHR ||
            prop == CL_CGL_SHAREGROUP_KHR
            push!(cl_props, cl_context_properties(val))
        else
            throw(OpenCLException("Invalid OpenCL Context property"))
        end
    end
    push!(cl_props, cl_context_properties(C_NULL))
    return cl_props
end

function num_devices(ctx::Context)
    ndevices = Ref{CL_uint}()
    @check api.clGetContextInfo(ctx.id, CL_CONTEXT_NUM_DEVICES,
                                sizeof(CL_uint), ndevices, C_NULL)
    return ndevices[]
end

function devices(ctx::Context)
    n = num_devices(ctx)
    if n == 0
        return []
    end
    dev_ids = Vector{CL_device_id}(n)
    @check api.clGetContextInfo(ctx.id, CL_CONTEXT_DEVICES,
                                n * sizeof(CL_device_id), dev_ids, C_NULL)
    return [Device(id) for id in dev_ids]
end

function create_some_context()
    if isempty(platforms())
        throw(OpenCLException("No OpenCL.Platform available"))
    end
    gpu_devices = devices(:gpu)
    if !isempty(gpu_devices)
        for dev in gpu_devices
            local ctx::Context
            try
                ctx = Context(dev)
            catch
                continue
            end
            return ctx
        end
    end
    cpu_devices = devices(:cpu)
    if !isempty(cpu_devices)
        for dev in cpu_devices
            local ctx::Context
            try
                ctx = Context(dev)
            catch
                continue
            end
            return ctx
        end
    end
    if isempty(gpu_devices) && isempty(cpu_devices)
        throw(OpenCLException("Unable to create any OpenCL.Context, no available devices"))
    else
        throw(OpenCLException("Unable to create any OpenCL.Context, no devices worked"))
    end
end


text:
# work around for julia issue #28948, Gurobi.jl issue #152
if VERSION ≤ v"0.7-"
    sparse_transpose(A) = sparse(transpose(A))
else
    sparse_transpose(A) = SparseMatrixCSC(transpose(A))
end

## Add Linear constraints

function constrainttype(rel)
    if rel == Cchar('>') || rel == XPRS_GEQ
        return XPRS_GEQ
    elseif rel == Cchar('<') || rel == XPRS_LEQ
        return XPRS_LEQ
    end
    # elseif rel == Cchar('=')
    return XPRS_EQ
end

function add_constr!(model::Model, inds::Vector{Cint}, coeffs::Vector{Float64}, rel::Cchar, rhs::Float64)

    rel = constrainttype(rel)

    length(inds) == length(coeffs) || error("Inconsistent argument dimensions.")

    ret = @xprs_ccall(addrows, Cint,(
        Ptr{Nothing}, #prob
        Cint,    #number new rowws
        Cint, #number nonzeros
        Ptr{Cchar}, #cstr type
        Ptr{Float64}, # rhs
        Ptr{Float64}, # range (size new row)
        Ptr{Cint}, # start (size newrow)
        Ptr{Cint}, # ind (size new nz)
        Ptr{Float64} # val (size newnz)
        ),
        model.ptr_model, 1, length(inds), Cchar[rel], Float64[rhs], C_NULL, Cint[0], inds .- Cint(1), coeffs
    )

    if ret != 0
        throw(XpressError(model))
    end

    return nothing
end

"""
    add_constr!{I<:Integer,R<:Real}(model::Model, inds::Vector{I}, coeffs::Vector{R}, rel::GChars, rhs::Real)

Adds a single constraint  to the model.
"""
function add_constr!(model::Model, inds::Vector{I}, coeffs::Vector{R}, rel::GChars, rhs::Real) where {I<:Integer,R<:Real}
    add_constr!(model, ivec(inds), fvec(coeffs), cchar(rel), Float64(rhs))
end

"""
    add_constr!(model::Model, coeffs::Vector, rel::GChars, rhs::Real)

Adds a constraint based on a dense vector os coefficients `coeffs`
"""
function add_constr!(model::Model, coeffs::Vector, rel::GChars, rhs::Real)
    inds = Compat.findall(!iszero, coeffs)
    vals = coeffs[inds]
    add_constr!(model, inds, vals, rel, rhs)
end

"""
    add_constrs!(model::Model, cbeg::Vector, inds::Vector, coeffs::Vector, rel::GCharOrVec, rhs::Vector)

Adds multiple rows in sparse(not julia standard) form (see Xpress manual)

    add_constrs!(model::Model, A::CoeffMat, rel::GCharOrVec, b::Vector{Float64})

Adds multiple rows in dense or SparseMatrixCSC form (see Xpress manual)
"""
function add_constrs!(model::Model, cbegins::Vector{Cint}, inds::Vector{Cint}, coeffs::Vector{Float64}, rel::Vector{Cchar}, rhs::Vector{Float64})

    for i in 1:length(rel)
        rel[i] = constrainttype(rel[i])
    end

    m = length(cbegins)
    nnz = length(inds)
    (m == length(rel) == length(rhs) && nnz == length(coeffs)) || error("Inconsistent argument dimensions.")

    if m > 0
        ret = @xprs_ccall(addrows, Cint,(
            Ptr{Nothing}, #prob
            Cint,    #number new rowws
            Cint, #number nonzeros
            Ptr{Cchar}, #cstr type
            Ptr{Float64}, # rhs
            Ptr{Float64}, # range (size new row) - no range constraints
            Ptr{Cint}, # start (size newrow) - offsets of coefficientes by each new row
            Ptr{Cint}, # ind (size new nz)
            Ptr{Float64} # val (size newnz)
            ),
            model.ptr_model, m, nnz, rel, rhs, C_NULL, cbegins .- Cint(1), inds .- Cint(1), coeffs
        )

        if ret != 0
            throw(XpressError(model))
        end
    end
    nothing
end
function add_constrs!(model::Model, cbeg::Vector, inds::Vector, coeffs::Vector, rel::GCharOrVec, rhs::Vector)
    add_constrs!(model, ivec(cbeg), ivec(inds), fvec(coeffs), cvecx(rel, length(cbeg)), fvec(rhs))
end
function add_constrs!(model::Model, A::CoeffMat, rel::GCharOrVec, b::Vector{Float64})
    m, n = size(A)
    (m == length(b) && n == num_vars(model)) || error("Incompatible argument dimensions.")
    add_constrs_t!(model, sparse_transpose(A), rel, b)
end
@static if VERSION > v"0.7-"
function add_constrs!(model::Model, A::Adjoint{Float64, MAT}, rel::GCharOrVec, b::Vector{Float64}) where MAT <: CoeffMat
    m, n = size(A)
    (m == length(b) && n == num_vars(model)) || error("Incompatible argument dimensions.")
    add_constrs_t!(model, sparse_transpose(A), rel, b)
end
end

"""
    add_constrs_t!(model::Model, At::SparseMatrixCSC{Float64}, rel::GCharOrVec, b::Vector)

add multiple constrints given in transpose sparse form.

    add_constrs_t!(model::Model, At::Matrix{Float64}, rel::GCharOrVec, b::Vector)

add multiple constrints given in transpose dense form.
"""
function add_constrs_t!(model::Model, At::SparseMatrixCSC{Float64}, rel::GCharOrVec, b::Vector)
    n, m = size(At)
    (m == length(b) && n == num_vars(model)) || error("Incompatible argument dimensions.")
    add_constrs!(model, At.colptr[1:At.n], rowvals(At), nonzeros(At), rel, b)
end
function add_constrs_t!(model::Model, At::AbstractMatrix{Float64}, rel::GCharOrVec, b::Vector)
    n, m = size(At)
    (m == length(b) && n == num_vars(model)) || error("Incompatible argument dimensions.")
    add_constrs_t!(model, sparse(At), rel, b)
end

"""
    add_rangeconstr!(model::Model, inds::Vector, coeffs::Vector, lb::Real, ub::Real)

Adds single range constraint in sparse format

    add_rangeconstr!(model::Model, coeffs::Vector, lb::Real, ub::Real)

Adds single range constraint in dense format
"""
function add_rangeconstr!(model::Model, inds::IVec, coeffs::FVec, lb::Float64, ub::Float64)
    # b = ub
    r = ub - lb
    ret = @xprs_ccall(addrows, Cint,(
        Ptr{Nothing}, #prob
        Cint,    #number new rowws
        Cint, #number nonzeros
        Ptr{Cchar}, #cstr type
        Ptr{Float64}, # rhs
        Ptr{Float64}, # range (size new row)
        Ptr{Cint}, # start (size newrow)
        Ptr{Cint}, # ind (size new nz)
        Ptr{Float64} # val (size newnz)
        ),
        model.ptr_model, 1, length(inds), Cchar['R'], Float64[ub], Float64[r], Cint[0], inds .- Cint(1), coeffs
    )
    if ret != 0
        throw(XpressError(model))
    end
    nothing
end
function add_rangeconstr!(model::Model, inds::Vector, coeffs::Vector, lb::Real, ub::Real)
    add_rangeconstr!(model, ivec(inds), fvec(coeffs), Float64(lb), Float64(ub))
end
function add_rangeconstr!(model::Model, coeffs::Vector, lb::Real, ub::Real)
    inds = find(coeffs)
    vals = coeffs[inds]
    add_rangeconstr!(model, inds, vals, lb, ub)
end

"""
    add_rangeconstrs!(model::Model, cbeg::Vector, inds::Vector, coeffs::Vector, lb::Vector, ub::Vector)

Adds multiple rows in sparse(not julia standard) form (see Xpress manual)

    add_rangeconstrs!(model::Model, A::CoeffMat, lb::Vector, ub::Vector)

Adds multiple rows in dense or SparseMatrixCSC form (see Xpress manual)
"""
function add_rangeconstrs!(model::Model, cbegins::IVec, inds::IVec, coeffs::FVec, lb::FVec, ub::FVec)
    m = length(cbegins)
    nnz = length(inds)
    (m == length(lb) == length(ub) && nnz == length(coeffs)) || error("Incompatible argument dimensions.")

    if m > 0
        r = ub .- lb
        ret = @xprs_ccall(addrows, Cint,(
            Ptr{Nothing}, #prob
            Cint,    #number new rowws
            Cint, #number nonzeros
            Ptr{Cchar}, #cstr type
            Ptr{Float64}, # rhs
            Ptr{Float64}, # range (size new row)
            Ptr{Cint}, # start (size newrow)
            Ptr{Cint}, # ind (size new nz)
            Ptr{Float64} # val (size newnz)
            ),
            model.ptr_model, m, nnz, cvecx( convert(Cchar,'R') , m), ub, r, cbegins .- Cint(1), inds .- Cint(1), coeffs
        )

        if ret != 0
            throw(XpressError(model))
        end
    end
    nothing
end
function add_rangeconstrs!(model::Model, cbeg::Vector, inds::Vector, coeffs::Vector, lb::Vector, ub::Vector)
    add_rangeconstrs!(model, ivec(cbeg), ivec(inds), fvec(coeffs), fvec(lb), fvec(ub))
end
function add_rangeconstrs!(model::Model, A::CoeffMat, lb::Vector, ub::Vector)
    m, n = size(A)
    (m == length(lb) == length(ub) && n == num_vars(model)) || error("Incompatible argument dimensions.")
    add_rangeconstrs_t!(model, sparse_transpose(A), lb, ub)
end

"""
    add_rangeconstrs_t!(model::Model, At::SparseMatrixCSC{Float64}, lb::Vector, ub::Vector)

add multiple constrints given in transpose sparse form.

    add_rangeconstrs_t!(model::Model, At::Matrix{Float64}, lb::Vector, ub::Vector)

add multiple constrints given in transpose dense form.
"""
function add_rangeconstrs_t!(model::Model, At::SparseMatrixCSC{Float64}, lb::Vector, ub::Vector)
    add_rangeconstrs!(model, At.colptr[1:At.n], rowvals(rowval), At.nzval, lb, ub)
end

function add_rangeconstrs_t!(model::Model, At::Matrix{Float64}, lb::Vector, ub::Vector)
    add_rangeconstrs_t!(model, sparse(At), lb, ub)
end

"""
    get_constrmatrix(model::Model)

Return constraint matrix (A) in SparseMatrixCSC form.
"""
function get_constrmatrix(model::Model)
    nnz = num_cnzs(model)
    m = Cint(num_constrs(model))
    n = Cint(num_vars(model))
    numnzP = Array{Cint}(undef,  1)
    cbeg = zeros(Cint, m+1)
    cind = Array{Cint}(undef,  nnz)
    cval = Array{Float64}(undef,  nnz)
    ret = @xprs_ccall(getrows, Cint, (
                     Ptr{Nothing},
                     Ptr{Cint},
                     Ptr{Cint},
                     Ptr{Float64},
                     Cint,
                     Ptr{Cint},
                     Cint,
                     Cint
                     ),
                     model.ptr_model, cbeg, cind, cval, Cint(nnz), numnzP, Cint(0), m.-Cint(1))
    if ret != 0
        throw(XpressError(model))
    end
    cbeg[end] = nnz
    I = Array{Int}(undef,  nnz)
    J = Array{Int}(undef,  nnz)
    V = Array{Float64}(undef,  nnz)
    for i in 1:length(cbeg)-1
        for j in (cbeg[i]+1):cbeg[i+1]
            I[j] = i
            J[j] = cind[j]+1
            V[j] = cval[j]
        end
    end
    return sparse(I, J, V, m, n)
end
function get_rows_nnz(model::Model, first::Integer, last::Integer)
    numnzP = Array{Cint}(undef,  1)
    ret = @xprs_ccall(getrows, Cint, (
                     Ptr{Nothing},
                     Ptr{Cint},
                     Ptr{Cint},
                     Ptr{Float64},
                     Cint,
                     Ptr{Cint},
                     Cint,
                     Cint
                     ),
                     model.ptr_model, C_NULL, C_NULL, C_NULL, Cint(0), numnzP, Cint(first-1), Cint(last-1))

    if ret != 0
        throw(XpressError(model))
    end
    return Int(numnzP[1])
end
function get_rows(model::Model, first::Integer, last::Integer)
    nnz = get_rows_nnz(model, first, last)
    m = Cint(last-first+1) #nrows to get coefs
    n = Cint(num_vars(model))
    numnzP = Array{Cint}(undef,  1)
    cbeg = zeros(Cint, m+1)
    cind = Array{Cint}(undef,  nnz)
    cval = Array{Float64}(undef,  nnz)
    ret = @xprs_ccall(getrows, Cint, (
                     Ptr{Nothing},
                     Ptr{Cint},
                     Ptr{Cint},
                     Ptr{Float64},
                     Cint,
                     Ptr{Cint},
                     Cint,
                     Cint
                     ),
                     model.ptr_model, cbeg, cind, cval, Cint(nnz), numnzP, Cint(first-1), Cint(last-1))
    if ret != 0
        throw(XpressError(model))
    end
    cbeg[end] = nnz
    I = Array{Int}(undef,  nnz)
    J = Array{Int}(undef,  nnz)
    V = Array{Float64}(undef,  nnz)
    for i in 1:length(cbeg)-1
        for j in (cbeg[i]+1):cbeg[i+1]
            I[j] = i
            J[j] = cind[j]+1
            V[j] = cval[j]
        end
    end
    # @show (I, J, V, m, n)
    return sparse(I, J, V, m, n)
end

"""
    add_sos!(model::Model, sostype::Symbol, idx::Vector{Cint}, weight::Vector{Float64})

Add SOS constraint of type `sostype`
Options are `:SOS1` and `:SOS2`
"""
add_sos!(model::Model, sostype::Symbol, idx::Vector{I}, weight::Vector{R}) where {I<:Integer, R<:Real} = add_sos!(model, sostype, ivec(idx), fvec(weight))
function add_sos!(model::Model, sostype::Symbol, idx::Vector{Cint}, weight::Vector{Float64})
    ((nelem = length(idx)) == length(weight)) || error("Index and weight vectors of unequal length")
    (sostype == :SOS1) ? (typ = XPRS_SOS_TYPE1) : ( (sostype == :SOS2) ? (typ = XPRS_SOS_TYPE2) : error("Invalid SOS constraint type") )
    ret = @xprs_ccall(addsets, Cint, (
                     Ptr{Nothing},
                     Cint,
                     Cint,
                     Ptr{Cchar},
                     Ptr{Cint},
                     Ptr{Cint},
                     Ptr{Float64}
                     ),
                     model.ptr_model,
                     convert(Cint, 1),
                     convert(Cint, nelem),
                     Cchar[typ],
                     Cint[0,0],
                     convert(Vector{Cint}, idx .- Cint(1)),
                     weight)
    if ret != 0
        throw(XpressError(model))
    end
end
function del_sos!(model::Model, idx::Vector{Cint})
    # int XPRS_CC XPRSdelsets(XPRSprob prob, int ndelsets, const int mindex[])
    numdel = length(idx)
    ret = @xprs_ccall(delsets, Cint, (
                     Ptr{Nothing},
                     Cint,
                     Ptr{Cint}),
                     model.ptr_model, convert(Cint,numdel), idx .- Cint(1))
    if ret != 0
        throw(XpressError(model))
    end
end

function get_sos_matrix(m::Model)
    nnz = num_setmembers(m)
    sets = num_sos(m)

    n = Cint(num_vars(m))

    settypes = Array{Cchar}(undef, sets)
    setstart = Array{Cint}(undef, sets+1)
    setcols = Array{Cint}(undef, nnz)
    setvals = Array{Float64}(undef, nnz)

    intents = Array{Cint}(undef, 1)
    nsets = Array{Cint}(undef, 1)

    # int XPRS_CC XPRSgetglobal(XPRSprob prob, int*nglents, int*sets,
    #char qgtype[], int mgcols[], double dlim[], char qstype[],
    #int msstart[],int mscols[], double dref[]);
    ret = @xprs_ccall(getglobal, Cint, (
        Ptr{Nothing},
        Ptr{Cint}, # nglents
        Ptr{Cint}, # sets
        Ptr{Cchar}, # qgtype
        Ptr{Cint}, # mgcols
        Ptr{Float64}, # dlim
        Ptr{Cchar}, # qstype
        Ptr{Cint}, # msstart
        Ptr{Cint}, # mscols
        Ptr{Float64}, # dref
        ),
        m.ptr_model, intents, nsets, C_NULL, C_NULL, C_NULL,
        settypes, setstart, setcols, setvals)

    if ret != 0
        throw(XpressError(m))
    end
    setstart[end] = nnz
    I = Array{Int}(undef,  nnz)
    J = Array{Int}(undef,  nnz)
    V = Array{Float64}(undef,  nnz)
    for i in 1:length(setstart)-1
        for j in (setstart[i]+1):setstart[i+1]
            I[j] = i
            J[j] = setcols[j]+1
            V[j] = setvals[j]
        end
    end
    # @show (I, J, V, sets, n)
    return sparse(I, J, V, sets, n), settypes
end

"""
    del_constrs!{T<:Real}(model::Model, idx::Vector{T})

Delete constraintes indexed in `idx`
"""
del_constrs!(model::Model, idx::T) where {T<:Real} = del_constrs!(model, Cint[idx])
del_constrs!(model::Model, idx::Vector{T}) where {T<:Real} = del_constrs!(model, convert(Vector{Cint},idx))
function del_constrs!(model::Model, idx::Vector{Cint})
    numdel = length(idx)
    ret = @xprs_ccall(delrows, Cint, (
                     Ptr{Nothing},
                     Cint,
                     Ptr{Cint}),
                     model.ptr_model, convert(Cint,numdel), idx .- Cint(1))
    if ret != 0
        throw(XpressError(model))
    end
end

"""
    chg_coeffs!{T<:Real, S<:Real}(model::Model, cidx::Vector{T}, vidx::Vector{T}, val::Vector{S})

Change multiple coefficients of the `A` matrix given constraints `cidx`, variables `vidx` and values `val`
"""
chg_coeffs!(model::Model, cidx::T, vidx::T, val::S) where {T<:Real, S<:Real} = chg_coeffs!(model, Cint[cidx], Cint[vidx], Float64[val])
chg_coeffs!(model::Model, cidx::Vector{T}, vidx::Vector{T}, val::Vector{S}) where {T<:Real, S<:Real} = chg_coeffs!(model, convert(Vector{Cint},cidx), convert(Vector{Cint},vidx), fvec(val))
function chg_coeffs!(model::Model, cidx::Vector{Cint}, vidx::Vector{Cint}, val::FVec)

    (length(cidx) == length(vidx) == length(val)) || error("Inconsistent argument dimensions.")

    numchgs = length(cidx)
    ret = @xprs_ccall(chgmcoef, Cint, (
                     Ptr{Nothing},
                     Cint,
                     Ptr{Cint},
                     Ptr{Cint},
                     Ptr{Float64}),
                     model, convert(Cint,numchgs), cidx .- Cint(1), vidx .- Cint(1), val)
    if ret != 0
        throw(XpressError(model))
    end
end

function chg_rhsrange!(model::Model, cidx::Vector{Cint}, val::FVec)
    # XPRSchgrhsrange(XPRSprob prob, int nels, const int mindex[], const double rng[])

    (length(cidx) == length(val)) || error("Inconsistent argument dimensions.")

    numchgs = length(cidx)
    ret = @xprs_ccall(chgrhsrange, Cint, (
                        Ptr{Nothing},
                        Cint,
                        Ptr{Cint},
                        Ptr{Float64}),
                        model, convert(Cint,numchgs), cidx .- Cint(1), val)
    if ret != 0
        throw(XpressError(model))
    end
end

function get_rhsrange!(model::Model, out::Vector{Float64}, rowb::Integer, rowe::Integer)

    _chklen(out, rowe-rowb+1)

    ret = @xprs_ccall(getrhsrange, Cint, (
                    Ptr{Nothing},
                    Ptr{Float64},
                    Cint,
                    Cint),
                    model, out, Cint(rowb-1), Cint(rowe-1))

    if ret != 0
        throw(XpressError(model))
    end

    return nothing
end

function get_rhsrange(model::Model, rowb::Integer, rowe::Integer)

    out = Array{Float64}(undef, rowe-rowb+1)

    get_rhsrange!(model, out, rowb, rowe)

    return out
end


text:
include("gym_tests.jl")
include("atari_tests.jl")
include("gridworld_tests.jl")
include("lavaworld_tests.jl")
include("pendulum_tests.jl")
include("continuumworld_tests.jl")
include("adversarial_mdp_tests.jl")
include("cartpole_tests.jl")

max_forks_repo_name: ancorso/POMDPGym


text:
#function to generate a Bernoulli(p) random variable
rand_bit(p::Real = 0.5) = rand() ≤ p ? 1 : 0   #\le + [TAB]

function queue_join_with_empty_check(time::Float64, current_station::Int, state::State, new_timed_events::Array{TimedEvent}, scenario::NetworkParameters, arrival_time::Float64)::Bool
    state.queues[current_station] += 1
    state.queues[current_station] == 1 && push!(new_timed_events, TimedEvent(EndOfServiceEvent(), time + rand(Gamma(1/scenario.gamma_shape, scenario.gamma_shape/scenario.μ_vector[current_station])), current_station, arrival_time))
    return true
end

"""
    new_timed_events = process_event(time, state, event)

Generate an array of 0 or more new `TimedEvent`s based on the current `event` and `state`.
"""
function process_event end # This defines a function with zero methods (to be added later)

"""
Terminates the simulation.
"""
function process_event(time::Float64, state::State, es_event::EndSimEvent)
    println("Ending simulation at time $time.")
    return []
end

"""
Prints out the system state at the current time after accepting a LogStateEvent as an arg.
"""
function process_event(time::Float64, state::State, ls_event::LogStateEvent)
    println("Logging state at time $time.")
    println(state)
    return []
end

"""
Handles jobs entering into the system
"""

# Process an arrival event
function process_event(time::Float64, state::State, location_ID, ::ArrivalEvent, scenario::NetworkParameters, arrival_time::Float64 = 0.0)
    # Increase number in system
    state.number_in_system += 1
    new_timed_events = TimedEvent[]

    # Set time ID
    arrival_time = time

    # Prepare next arrival
    push!(new_timed_events,TimedEvent(ArrivalEvent(),time + rand(Gamma(1/scenario.gamma_shape, scenario.gamma_shape/scenario.λ)),0, arrival_time))

    # Using probability vector, find the first station that the new arrival heads to
    current_station = sample(scenario.p_e)

    # Job enters initial station, if it is not full, the job either queues or is being actively serviced.
    if scenario.K[current_station] != -1
        state.queues[current_station] < scenario.K[current_station] && queue_join_with_empty_check(time, current_station, state, new_timed_events, scenario, arrival_time) && return new_timed_events
    else
        queue_join_with_empty_check(time, current_station, state, new_timed_events, scenario, arrival_time) && return new_timed_events
    end

    # If the job doesnt leave the system (Determined by the overflow matrix Q), then a new OverflowEvent is created
    overflow_state_control(time, state, scenario.Q, scenario, arrival_time, current_station, new_timed_events)

    return new_timed_events
end

# Process an EndOfServiceEvent event
function process_event(time::Float64, state::State, location_ID, ::EndOfServiceEvent, scenario::NetworkParameters, arrival_time::Float64)
    new_timed_events = TimedEvent[]

    current_station = location_ID

    state.queues[current_station] -= 1

    overflow_state_control(time, state, scenario.P, scenario, arrival_time, current_station, new_timed_events)

    @assert state.queues[location_ID] > -1
    state.queues[location_ID] > 0 && push!(new_timed_events, TimedEvent(EndOfServiceEvent(), time + rand(Gamma(1/scenario.gamma_shape, scenario.gamma_shape/scenario.μ_vector[current_station])), location_ID, arrival_time))

    return new_timed_events
end

# Process an Orbit and OverflowEvent event, here, the job has finished orbiting and must either now leave the system, join the queue, or join another queue.
function process_event(time::Float64, state::State, location_ID, ::OverflowEvent, scenario::NetworkParameters, arrival_time::Float64)
    new_timed_events = TimedEvent[]
    current_station = location_ID

    # Job enters initial station, if it is not full, the job either queues or is being actively serviced.
    state.orbiting_jobs -= 1

    if scenario.K[current_station] != -1
        state.queues[current_station] < scenario.K[current_station] && queue_join_with_empty_check(time, current_station, state, new_timed_events, scenario, arrival_time) && return new_timed_events
    else
        queue_join_with_empty_check(time, current_station, state, new_timed_events, scenario, arrival_time) && return new_timed_events
    end

    overflow_state_control(time, state, scenario.Q, scenario, arrival_time, current_station, new_timed_events)

    return new_timed_events
end

# Process an EndSimEvent event
function process_event(time::Float64, state::State, location_ID, es_event::EndSimEvent, scenario::NetworkParameters, arrival_time::Float64)
    return []
end

function overflow_state_control(time::Float64, state::State, prob_matrix::Matrix{Float64}, scenario::NetworkParameters, arrival_time::Float64, current_station, new_timed_events::Array{TimedEvent})::Nothing
    if rand_bit(1-sum(prob_matrix[current_station,:])) == 1
        state.number_in_system_decreased = true
        state.number_in_system -= 1
        @assert state.number_in_system >= 0
    else
        push!(new_timed_events, TimedEvent(OverflowEvent(), time + rand(Gamma(1/scenario.gamma_shape, scenario.gamma_shape/scenario.η)), sample(AnalyticWeights(prob_matrix[current_station,:])), arrival_time))
        state.orbiting_jobs += 1
    end
    return nothing
end


text:
using NiLang, NiLang.AD
using Test, LinearAlgebra

@testset "inorm2, dot" begin
    out = 0.0
    vec = [1.0, 2.0, 3.0]
    vec2 = [1.0, 2.0, 5.0]
    @instr inorm2(out, vec)
    @test out ≈ norm(vec)^2
    @test check_inv(inorm2, (out, vec))

    out = 0.0
    vec = [1.0, 2.0, 3.0]
    vec2 = [1.0, 2.0, 5.0]
    @test check_grad(inorm2, (Loss(out), vec); verbose=true)

    out = 0.0
    @instr idot(out, vec, vec2)
    @test out ≈ dot(vec, vec2)
    @test check_inv(idot, (out, vec, vec2))

    @test check_grad(idot, (Loss(0.0), vec, vec2); verbose=true)

    m = randn(4,4)
    n = randn(4,4)
    out = 0.0
    @instr idot(out, m[:,2], n[:,4])
    @test out ≈ dot(m[:,2], n[:,4])
    @test check_inv(idot, (out, m[:,2], n[:,4]))

    @test check_grad(idot, (Loss(0.0), vec, vec2); verbose=true)
end


text:
# Ост - "Е6" - восток, Норд - "N8" - север, Запад - "W4" - запад, Юг - "S2" - юг.
function main_function12!(r::Robot, cell_size::Int)
    path = go_to_west_south_corner_and_return_path!(r)
    x=0; y=0
    horisontalDirection = Ost

    while !(isborder(r, Nord) && isborder(r, Ost))
        marker_special!(r, x, y, cell_size)
        if move_up_condition(r)
            move!(r, Nord)
            y += 1
            marker_special!(r, x, y, cell_size)
            horisontalDirection = inverse_side(horisontalDirection)
        end

        move!(r,horisontalDirection)
        (horisontalDirection == Ost) ? x += 1 : x -= 1
    end

    marker_special!(r, x, y, cell_size)

    go_to_west_south_corner_and_return_path!(r)
    go_by_path!(r, path)
end

function marker_special!(r, x, y, cell_size)
    if (mod(x, 2 * cell_size)) < cell_size && (mod(y, 2 * cell_size)) < cell_size ||
        (mod(x + cell_size, 2 * cell_size)) < cell_size && (mod(y, 2 * cell_size)) >= cell_size
        putmarker!(r)
    end
end

function move_up_condition(r)
    return isborder(r, Ost) || isborder(r, West) && !(isborder(r, Sud) && isborder(r, West))
end
function clockwise_side(side::HorizonSide)::HorizonSide
    return HorizonSide(mod(Int(side)-1,4))
end
function inverse_side(side::HorizonSide)::HorizonSide
    return HorizonSide(mod(Int(side)+2, 4))
end
function go_to_west_south_corner_and_return_path!(r::Robot; go_around_barriers::Bool = false, markers = false)::Array{Tuple{HorizonSide,Int64},1}
    my_ans = []
    a = go_to_border_and_return_path!(r, West; go_around_barriers, markers)
    b = go_to_border_and_return_path!(r, Sud; go_around_barriers, markers)

    for i in a
        push!(my_ans, i)
    end
    for i in b
        push!(my_ans, i)
    end
    return my_ans
end
function go_by_path!(r::Robot, path::Array{Tuple{HorizonSide,Int64},1})
    new_path = reverse(path)
    for i in new_path
        go!(r, i[1]; steps = i[2])
    end
end


text:
function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ArgumentConstructStage,
    model::DeviceModel{R, D},
    ::Type{S},
) where {
    R <: PSY.RenewableGen,
    D <: AbstractRenewableDispatchFormulation,
    S <: PM.AbstractPowerModel,
}
    devices = get_available_components(R, sys)

    add_variables!(container, ActivePowerVariable, devices, D())
    add_variables!(container, ReactivePowerVariable, devices, D())

    add_parameters!(container, ActivePowerTimeSeriesParameter, devices, model)

    add_expressions!(container, ProductionCostExpression, devices, model)

    # Expression
    add_to_expression!(
        container,
        ActivePowerBalance,
        ActivePowerVariable,
        devices,
        model,
        S,
    )
    add_to_expression!(
        container,
        ReactivePowerBalance,
        ReactivePowerVariable,
        devices,
        model,
        S,
    )
    if has_service_model(model)
        add_to_expression!(
            container,
            ActivePowerRangeExpressionLB,
            ActivePowerVariable,
            devices,
            model,
            S,
        )
        add_to_expression!(
            container,
            ActivePowerRangeExpressionUB,
            ActivePowerVariable,
            devices,
            model,
            S,
        )
    end
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ModelConstructStage,
    model::DeviceModel{R, D},
    ::Type{S},
) where {
    R <: PSY.RenewableGen,
    D <: AbstractRenewableDispatchFormulation,
    S <: PM.AbstractPowerModel,
}
    devices = get_available_components(R, sys)

    if has_service_model(model)
        add_constraints!(
            container,
            ActivePowerVariableLimitsConstraint,
            ActivePowerRangeExpressionUB,
            devices,
            model,
            S,
        )
    else
        add_constraints!(
            container,
            ActivePowerVariableLimitsConstraint,
            ActivePowerVariable,
            devices,
            model,
            S,
        )
    end

    add_constraints!(
        container,
        ReactivePowerVariableLimitsConstraint,
        ReactivePowerVariable,
        devices,
        model,
        S,
    )
    add_feedforward_constraints!(container, model, devices)

    objective_function!(container, devices, model, S)

    add_constraint_dual!(container, sys, model)
    return
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ArgumentConstructStage,
    model::DeviceModel{R, D},
    ::Type{S},
) where {
    R <: PSY.RenewableGen,
    D <: AbstractRenewableDispatchFormulation,
    S <: PM.AbstractActivePowerModel,
}
    devices = get_available_components(R, sys)

    add_variables!(container, ActivePowerVariable, devices, D())

    add_parameters!(container, ActivePowerTimeSeriesParameter, devices, model)

    add_expressions!(container, ProductionCostExpression, devices, model)

    # Expression
    add_to_expression!(
        container,
        ActivePowerBalance,
        ActivePowerVariable,
        devices,
        model,
        S,
    )
    if has_service_model(model)
        add_to_expression!(
            container,
            ActivePowerRangeExpressionLB,
            ActivePowerVariable,
            devices,
            model,
            S,
        )
        add_to_expression!(
            container,
            ActivePowerRangeExpressionUB,
            ActivePowerVariable,
            devices,
            model,
            S,
        )
    end
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ModelConstructStage,
    model::DeviceModel{R, D},
    ::Type{S},
) where {
    R <: PSY.RenewableGen,
    D <: AbstractRenewableDispatchFormulation,
    S <: PM.AbstractActivePowerModel,
}
    devices = get_available_components(R, sys)

    if has_service_model(model)
        add_constraints!(
            container,
            ActivePowerVariableLimitsConstraint,
            ActivePowerRangeExpressionUB,
            devices,
            model,
            S,
        )
    else
        add_constraints!(
            container,
            ActivePowerVariableLimitsConstraint,
            ActivePowerVariable,
            devices,
            model,
            S,
        )
    end
    add_feedforward_constraints!(container, model, devices)

    objective_function!(container, devices, model, S)

    add_constraint_dual!(container, sys, model)

    return
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ArgumentConstructStage,
    model::DeviceModel{R, FixedOutput},
    ::Type{S},
) where {R <: PSY.RenewableGen, S <: PM.AbstractPowerModel}
    devices = get_available_components(R, sys)

    add_parameters!(container, ActivePowerTimeSeriesParameter, devices, model)
    add_parameters!(container, ReactivePowerTimeSeriesParameter, devices, model)

    add_to_expression!(
        container,
        ActivePowerBalance,
        ActivePowerTimeSeriesParameter,
        devices,
        model,
        S,
    )
    add_to_expression!(
        container,
        ReactivePowerBalance,
        ReactivePowerTimeSeriesParameter,
        devices,
        model,
        S,
    )
    return
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ArgumentConstructStage,
    model::DeviceModel{R, FixedOutput},
    ::Type{S},
) where {R <: PSY.RenewableGen, S <: PM.AbstractActivePowerModel}
    devices = get_available_components(R, sys)

    add_parameters!(container, ActivePowerTimeSeriesParameter, devices, model)
    add_to_expression!(
        container,
        ActivePowerBalance,
        ActivePowerTimeSeriesParameter,
        devices,
        model,
        S,
    )
    return
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ModelConstructStage,
    model::DeviceModel{R, FixedOutput},
    ::Type{S},
) where {R <: PSY.RenewableGen, S <: PM.AbstractPowerModel}
    # FixedOutput doesn't add any constraints to the model. This function covers
    # AbstractPowerModel and AbtractActivePowerModel
    return
end


text:
using Glob, Images
let
    files = glob("./scripts/output/*.png")
    for i in 1:length(files)
        img_source = load(files[i])
        if size(img_source)[1] <700 || size(img_source)[2] <700
            percentage_scale = 0.6
        else
            percentage_scale = 0.33
        end
        new_size = trunc.(Int, size(img_source) .* percentage_scale)
        img_rescaled = imresize(img_source, new_size)
        save("./thumbnails/"*files[i][18:end-4]* "_thumb"*".png", img_rescaled)
    end
end

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
##################################################
#   Override default handling (default = inline) #
##################################################

get_latex_function(ode::DiffEqBase.AbstractParameterizedFunction) = latexalign
function get_latex_function(x::AbstractArray{T}) where T <: DiffEqBase.AbstractParameterizedFunction
    return latexalign
end

get_md_function(args::DiffEqBase.AbstractParameterizedFunction) = mdalign
function get_md_function(x::AbstractArray{T}) where T <: DiffEqBase.AbstractParameterizedFunction
    return mdalign
end

###############################################
#         Overload environment functions      #
###############################################

function latexraw(ode::DiffEqBase.AbstractParameterizedFunction; kwargs...)
    lhs = ["\\frac{d$x}{dt} = " for x in ode.syms]
    rhs = latexraw(ode.funcs; kwargs...)
    return lhs .* rhs
end

function latexinline(ode::DiffEqBase.AbstractParameterizedFunction; kwargs...)
    lhs = ["\\frac{d$x}{dt} = " for x in ode.syms]
    rhs = latexraw(ode.funcs; kwargs...)
    return latexstring.( lhs .* rhs )
end

function latexalign(ode::DiffEqBase.AbstractParameterizedFunction; field::Symbol=:funcs, bracket=false, kwargs...)
    lhs = [Meta.parse("d$x/dt") for x in ode.syms]
    rhs = getfield(ode, field)
    if bracket
        rhs = add_brackets(rhs, ode.syms)
        lhs = [:(d[$x]/dt) for x in ode.syms]
    end
    return latexalign(lhs, rhs; kwargs...)
end

"""
latexalign(odearray; field=:funcs)

Display ODEs side-by-side.
"""
function latexalign(odearray::AbstractVector{T}; field::Symbol=:funcs, kwargs...) where T<:DiffEqBase.AbstractParameterizedFunction
    a = []
    maxrows = maximum(length.(getfield.(odearray, :syms)))

    blank = LaTeXString("")
    for ode in odearray
        nr_eq = length(ode.syms)

        lhs = [Meta.parse("d$x/dt") for x in ode.syms]
        rhs = getfield(ode, field)
        first_separator = fill(LaTeXString(" &= "), nr_eq)
        second_separator = fill(LaTeXString(" & "), nr_eq)
        if nr_eq < maxrows
            ### This breaks type-safety, but I don't think that it will be a bottle neck for anyone.
            lhs = [lhs; fill(blank, maxrows - nr_eq)]
            rhs = [rhs; fill(blank, maxrows - nr_eq)]
            first_separator = [first_separator; fill(LaTeXString(" & "), maxrows-nr_eq)]
            second_separator = [second_separator; fill(LaTeXString(" & "), maxrows-nr_eq)]
        end

        append!(a, [lhs, first_separator, rhs, second_separator])
    end
    a = reduce(hcat, a)
    return latexalign(a; separator=" ", kwargs...)
end


text: const SECRET_TOKEN = "671dd48fd32aab34ca3be3997aef33a5731a5fa5f6d09bec7a4f7be58967ada9"


text:
using StatsBase

include("functions_align.jl")
include("functions_encoding.jl")

function get_difference_matrix(target_A, query_A)

    difference_array = Array{Float64,2}(undef, size(target_A, 2), size(query_A, 2))

    Threads.@threads for query_index in 1:size(query_A, 2)
        for target_index in 1:size(target_A, 2)
            differences = 0
            denominator = 0
            for r in 1:size(target_A, 1)

                @inbounds x = target_A[r,target_index]
                @inbounds y = query_A[r,query_index]

                different = (x & y) < 16
                same = (x & 8 == 8) && x == y

                if different
                    differences += 1
                end

                if different || same
                    denominator += 1
                end
            end

            difference_array[target_index, query_index] = convert(AbstractFloat, differences) / convert(AbstractFloat, denominator)

        end
    end

    return difference_array
end

function get_SNP_distance_matrix(target_A, query_A, ref_A)
    difference_array = Array{Int64,2}(undef, size(target_A, 2), size(query_A, 2))

    query_difference_list = get_list_of_differences(query_A, ref_A)
    target_difference_list = get_list_of_differences(target_A, ref_A)

    Threads.@threads for query_index in 1:size(query_A, 2)
        for target_index in 1:size(target_A, 2)
            differences = 0
            for r in query_difference_list[query_index]

                @inbounds x = target_A[r,target_index]
                @inbounds y = query_A[r,query_index]

                different = (x & y) < 16

                if different
                    differences += 1
                end

            end

            for r in target_difference_list[target_index]
                if r in query_difference_list[query_index]
                    continue
                end

                @inbounds x = target_A[r,target_index]
                @inbounds y = query_A[r,query_index]

                different = (x & y) < 16

                if different
                    differences += 1
                end

            end

            difference_array[target_index, query_index] = differences
        end
    end

    return difference_array
end

function get_min_float_indices(array_of_scores)

    min::Float64 = array_of_scores[1]
    indices::Array{Int64,1} = [1]

    for (i, score) in enumerate(array_of_scores)
        if i == 1
            min = score
            indices = [i]
        elseif score == min
            indices = vcat(indices, i)
        elseif score < min
            min = score
            indices = [i]
        end
    end

    return min, indices
end

function get_max_int_indices(array_of_scores)

    max::Int64 = array_of_scores[1]
    indices::Array{Int64,1} = [1]

    for (i, score) in enumerate(array_of_scores)
        if i == 1
            max = score
            indices = [i]
        elseif score == max
            indices = vcat(indices, i)
        elseif score > max
            max = score
            indices = [i]
        end
    end

    return max, indices
end

function get_best_target_index(differences_V, completeness_V)
    distance_min, distance_min_indices = get_min_float_indices(differences_V)

    if length(distance_min_indices) > 1
        completenesses = completeness_V[distance_min_indices]
        completeness_max, completeness_max_indices = get_max_int_indices(completenesses)
        indx = distance_min_indices[completeness_max_indices[1]]
    else
        indx = distance_min_indices[1]
    end

    return indx
end

function get_SNPs(target_V, query_V)
    nuc_dict = make_nuc_dict()

    SNPs = Array{String, 1}()
    for r in 1:length(target_V)
        @inbounds x = target_V[r]
        @inbounds y = query_V[r]
        different = (x & y) < 16
        if different
            nuc_T = nuc_dict[x]
            nuc_Q = nuc_dict[y]
            snp = string(r) * nuc_T * nuc_Q
            push!(SNPs, snp)
        end
    end

    return join(SNPs, "|")
end

function closest(target_file, query_file, outfile)
    @time T_A, T_names = populate_byte_array_get_names(target_file)
    @time Q_A, Q_names = populate_byte_array_get_names(query_file)

    @time differences = get_difference_matrix(T_A, Q_A)

    @time target_completeness = score_alignment2(T_A)

    @time open(outfile, "w") do io

        println(io, "query,closest,distance,SNPs")

        for query_index in 1:size(differences, 2)
            best_indx = get_best_target_index(differences[:,query_index], target_completeness)
            distance = differences[best_indx, query_index]
            SNPs = get_SNPs(Q_A[:,query_index], T_A[:,best_indx])
            Q_name = Q_names[query_index]
            T_name = T_names[best_indx]

            println(io, Q_name * "," * T_name * "," * string(round(distance, digits = 9)) * "," * SNPs)
        end
    end
end

function fast_closest(target_file, query_file, ref_file, outfile)
    R_A, R_names = populate_byte_array_get_names(ref_file)
    @time T_A, T_names = populate_byte_array_get_names(target_file)
    @time Q_A, Q_names = populate_byte_array_get_names(query_file)

    @time differences = get_SNP_distance_matrix(T_A, Q_A, R_A)

    @time target_completeness = score_alignment2(T_A)

    @time open(outfile, "w") do io

        println(io, "query,closest,SNP_distance,SNPs")

        for query_index in 1:size(differences, 2)
            best_indx = get_best_target_index(differences[:,query_index], target_completeness)
            distance = differences[best_indx, query_index]
            SNPs = get_SNPs(Q_A[:,query_index], T_A[:,best_indx])
            Q_name = Q_names[query_index]
            T_name = T_names[best_indx]

            println(io, Q_name * "," * T_name * "," * string(distance) * "," * SNPs)
        end
    end
end

max_stars_repo_name: cov-ert/julialign

max_issues_repo_name: cov-ert/julialign

max_forks_repo_name: cov-ert/julialign


text:
# Discord delegate functions so that the bot can be tested more easily with the Pretend framework.

@mockable function discord_channel(c::Client, channel_id::UInt64)
    return @discord retrieve(c, DiscordChannel, channel_id)
end

@mockable function discord_upload_file(
    c::Client, channel::DiscordChannel, filename::AbstractString; kwargs...
)
    return @discord upload_file(c, channel, filename; kwargs...)
end

@mockable function discord_reply(c::Client, m::Message, content::AbstractString)
    return @discord reply(c, m, content)
end


text:
function acqr_em_starting(y)
    """
    computing the starting point for the em algorithm for acQR model

    y: time series

    javier.cara@upm.es, 2020-05
    """

    #
    nt = length(y)

    # least squares for c
    # ya = c*yb + v
    ya = y[2:nt]
    yb = y[1:nt-1]
    c = sum(ya.*yb)/sum(yb.^2)

    # noise
    v = ya - c*yb
    r = var(v)

    #
    a = 0.5 # must be < 1

    # noise
    w = ya - a*yb
    q = var(w)

    # first state
    m1 = y[1]/c
    P1 = 1

    return acqr(a,c,q,r,m1,P1)
end


text:
"""
    MenuAdventures.Talking

A sub-module that enables talking to other players.

```jldoctest
julia> using MenuAdventures

julia> using MenuAdventures.Testing

julia> using MenuAdventures.Talking

julia> import MenuAdventures: ever_possible

julia> @universe struct Universe <: AbstractUniverse
        end;

julia> @noun struct Room <: AbstractRoom
            already_lit::Bool = true
        end;

julia> @noun struct Person <: Noun
        end;

julia> @noun struct NPC <: Noun
            dialog::Dict{String, Answer} = Dict{String, Answer}()
        end;

julia> ever_possible(::Say, ::Visible, ::NPC) = true;

julia> cd(joinpath(pkgdir(MenuAdventures), "test", "Talking")) do
            check_choices() do interface
                you = Person(
                    "Brandon",
                    grammatical_person = second_person,
                    indefinite_article = "",
                )
                friend = NPC(
                    "your friend",
                    indefinite_article = "",
                    dialog = Dict("Hello" => Answer("Hello", (universe, person) -> nothing))
                )
                room = Room("room")
                universe = Universe(you, interface = interface)
                universe[room, you] = Containing()
                universe[room, friend] = Containing()
                universe
            end
        end
true
```
"""
module Talking

using MenuAdventures: Action, Answer, Domain, ever_possible, println_wrapped, Sentence, subject_to_verb, Verb, Visible
import MenuAdventures: argument_domains, find_in_domain, possible_now, print_sentence

"""
    Dialog()

Things a player might [`Say`](@ref).

A [`Domain`](@ref).
"""
struct Dialog <: Domain end

export Dialog

# TODO: why won't it let me use an _ for universe?
function find_in_domain(universe, sentence, ::Dialog; lit = true)
    [Answer(repr(text), text) for text in keys(sentence.argument_answers[1].object.dialog)]
end

"""
    Say()

Say, to someone [`Visible`](@ref), some [`Dialog`](@ref).

An [`Action`](@ref).
Note the first argument is the addressee, and the second is the dialog.
If `possible(::Say, ::Dialog, addressee)`, then `addressee` must have a dialog field containing a dialog dictionary.
The keys of this dictionary should be text the player says, and the values should be [`Answer`](@ref)s containing a response from the `addressee` and a trigger.
After a dialog option is chosen, the response trigger will be called with two arguments: `universe` and `addressee`.
Then the dialog option will be removed.
"""
struct Say <: Action end

export Say

function possible_now(_, sentence::Sentence{Say}, domain::Visible, person)
    ever_possible(sentence.action, domain, person) && !isempty(person.dialog)
end

function argument_domains(::Say)
    Visible(), Dialog()
end

const REPLY = Verb("reply"; third_person_singular_present = "replies")

function (::Say)(universe, person, text)
    interface = universe.interface
    dialog = person.dialog
    answer = dialog[text]
    delete!(dialog, text)
    buffer = IOBuffer()
    show(IOContext(buffer, :capitalize => true, :is_subject => true), person)
    print(buffer, ' ')
    print(buffer, subject_to_verb(person, REPLY))
    print(buffer, " \"")
    print(buffer, answer.text)
    print(buffer, '"')
    println_wrapped(interface, String(take!(buffer)))
    answer.object(universe, person)
    return false
end

function print_sentence(io, ::Say, person_answer, text_answer)
    print(io, "Say ")
    print(io, text_answer.text)
    print(io, " to ")
    print(io, person_answer.text)
end

end


text:
# Contents: Functions for hashing passwords.

"""
Type for storing user-specific salt and hashed password.

The password is hashed by the constructor using the following algorithm:
1) Generate a 16 byte (128 bit) salt using a cryptographically secure RNG.
2) Hash the salted password using PBKDF2 with:
    - SHA-512 as the pseudorandom function
    - 5000 iterations
    - A 512-bit derived key length
"""
immutable StoredPassword
    salt::Array{UInt8, 1}
    hashed_password::Array{UInt8, 1}

    function StoredPassword(password::AbstractString)
	salt            = csrng(16)
	hashed_password = compute_hashed_password(salt, password)
	new(salt, hashed_password)
    end
end

"""
Returns true if hash(salt * password) == stored hashed password.

The algorithm:
1) Compute hash(salt, password), where the salt is supplied in a StoredPassword.
2) Return true if hash(password) == stored hashed password
"""
function password_is_valid(password::AbstractString, salt::Vector{UInt8}, hashed_pwd::Vector{UInt8})
    result = false
    hp     = compute_hashed_password(salt, password)
    if hp == hashed_pwd
	result = true
    end
    result
end
password_is_valid(password::AbstractString, sp::StoredPassword) = password_is_valid(password, sp.salt, sp.hashed_password)

"Selects and applies the password hashing algorithm."
function compute_hashed_password(salt::Array{UInt8, 1}, password::AbstractString)
    pbkdf2(salt, password, 5000, 64)    # 5000 iterations; dklen = 64 bytes = 512 bits
end

### EOF


text:
#! /usr/bin/env julia

println("Start test 'tests/lesson.jl'")

include("../src/DCellC.jl")
using DCellC

lmgs = [ synthesize(512, 512, (300, 500)) for i in 1:3 ]

region() = (rand(1:50), rand(1:50), rand(260:400), rand(260:400))
selections = Selection[]

for lmg in lmgs
  r = region()
  lmg = crop(lmg, r...)
  push!(selections, (lmg..., r))
end

lesson1 = Lesson(FCRNA,
                 imgtype = GreyscaleImage,
                 selections = selections,
                 patchmode = 2,
                 epochs = 2,
                 patchsize = 128)

lesson2 = Lesson(Multiscale3(GreyscaleImage, bn = true),
                 patchmode = 1,
                 epochs = 2,
                 patchsize = 64)

println("Test lesson.jl: Creating instance of Lesson works")

lname = tempname()
lessonsave(lname, lesson1)
println("Test lesson.jl: Saving Lessons works")

lesson3 = lessonload(lname)
println("Test lesson.jl: Loading Lessons works")

model  = train(lesson1, log=false)
model2 = train(lesson2, log=false)

println("Test lesson.jl: Learned the lesson")

println("Test lesson.jl completed")

max_stars_repo_name: tscode/DCellC

max_issues_repo_name: tscode/DCellC


text:
@testset "Ambiguities" begin
    if VERSION.prerelease == () && !Sys.iswindows() && VERSION < v"1.7.0"
        mbs = Test.detect_ambiguities(ManifoldsBase)
        # Interims solution until we follow what was proposed in
        # https://discourse.julialang.org/t/avoid-ambiguities-with-individual-number-element-identity/62465/2
        fmbs = filter(x -> !any(has_type_in_signature.(x, Identity)), mbs)
        FMBS_LIMIT = 22
        @test length(fmbs) <= FMBS_LIMIT
        if length(fmbs) > FMBS_LIMIT
            for amb in fmbs
                println(amb)
                println()
            end
        end
        ms = Test.detect_ambiguities(Manifolds)
        # Interims solution until we follow what was proposed in
        # https://discourse.julialang.org/t/avoid-ambiguities-with-individual-number-element-identity/62465/2
        fms = filter(x -> !any(has_type_in_signature.(x, Identity)), ms)
        FMS_LIMIT = 23
        if length(fms) > FMS_LIMIT
            for amb in fms
                println(amb)
                println()
            end
        end
        @test length(fms) <= FMS_LIMIT
        # this test takes way too long to perform regularly
        # @test length(our_base_ambiguities()) <= 4
    else
        @info "Skipping Ambiguity tests for pre-release versions"
    end
end


text:
using DataFrames, MixedModels,StatLib
dfx = readtable("/mapr/mapr04p/analytics0001/analytic_users/jmdl/v3.0/inputs/mctargetgrus00115.csv")
dfs=readtable("/mapr/mapr04p/analytics0001/analytic_users/jmdl/v3.0/sas_output/fixed1_grus00115.csv")
#mctargetgrus00147.csv  mctargetgrus00247.csv  mctargetgrus00256.csv
#Sample 1:for ppgid :115 region: US outlet: GR
f = logvolume ~ 1 + logregprice + logpriceindex + displayonly + featuredisplay + featureonly + logseason + absprice19 + lift19 + absprice220 + lift220 + absprice178 + lift178 + absprice132 + lift132 + absprice131 + lift131 + absprice191 + lift191 + absprice95 + lift95 + absprice140 + lift140 + absprice139 + lift139 + absprice31 + lift31 + absprice233 + lift233 + absprice289 + lift289 + absprice275 + lift275 + absprice236 + lift236 + absprice230 + lift230 + absprice235 + lift235 + absprice246 + lift246 + absprice265 + lift265 + absprice290 + lift290 + absprice239 + lift239 + absprice7 + lift7 + absprice305 + lift305 + absprice322 + lift322 + absprice302 + lift302 + absprice299 + lift299 + absprice229 + lift229 + absprice309 + lift309 + absprice231 + lift231 + absprice298 + lift298 + absprice311 + lift311 +
( ( 0 + logregprice ) | store ) +
( ( 0 + logpriceindex ) | store ) +
( ( 0 + displayonly ) | store ) +
( ( 0 + featuredisplay ) | store ) +
( ( 0 + featureonly ) | store ) +
( ( 0 + logseason ) | store )

lowercase!(dfx)
for col in[:weekending, :outlet, :region] dfx[col] = pool(dfx[col]) end

m=fit!(lmm(f, dfx), true)

mx = lmm(f, dfx)
reml!(mx,true)
m = fit(mx)

@time m1 = fit!(lmm(Y ~ 1 + (1|G) + (1|H), ml1m))

fit!(lmm(f, dfx), true)

gmm1 = fit!(glmm(f, dfd[(dfd[:buyer_pos_p1].==1),cols] , m[:dist], m[:lnk]), verbose=false, fast=true)
glm(f, dfd[(dfd[:buyer_pos_p1].==1),cols] , m[:dist], m[:lnk])

modelREML = lmm(f, dfx)
reml!(modelREML,true)
lmeModel = fit(modelREML)

Sample 2:for ppgid :247 region: US outlet: GR
logvolume ~ 1 + logregprice + logpriceindex + displayonly + featuredisplay + featureonly + logseason + absprice263 + lift263 + absprice258 + lift258 + absprice259 + lift259 + absprice254 + lift254 + absprice261 + lift261 + absprice265 + lift265 + absprice248 + lift248 + absprice260 + lift260 + absprice266 + lift266 + absprice255 + lift255 + absprice276 + lift276 + absprice271 + lift271 + absprice289 + lift289 + absprice228 + lift228 + absprice236 + lift236 + absprice242 + lift242 + absprice279 + lift279 + absprice282 + lift282 + absprice239 + lift239 + absprice230 + lift230 + absprice311 + lift311 + absprice10 + lift10 + lift189 + absprice115 + lift115 + absprice20 + lift20 + absprice33 + lift33 + absprice293 + lift293 + absprice12 + lift12 + absprice298 + lift298 + lift45 + ( ( 0 + logregprice ) | store ) + ( ( 0 + logpriceindex ) | store ) + ( ( 0 + displayonly ) | store ) + ( ( 0 + featuredisplay ) | store ) + ( ( 0 + featureonly ) | store ) + ( ( 0 + logseason ) | store )

Sample 3: for ppgid :147 region: US outlet: GR
logvolume ~ 1 + logregprice + logpriceindex + displayonly + logseason + absprice88 + lift88 + absprice48 + lift48 + absprice95 + lift95 + absprice194 + lift194 + absprice89 + lift89 + absprice227 + lift227 + absprice49 + lift49 + absprice65 + lift65 + absprice225 + lift225 + absprice33 + lift33 + absprice235 + lift235 + absprice272 + lift272 + absprice276 + lift276 + absprice288 + lift288 + absprice278 + lift278 + absprice265 + lift265 + absprice274 + lift274 + absprice228 + lift228 + absprice322 + lift322 + absprice13 + lift13 + absprice327 + lift327 + absprice236 + lift236 + absprice305 + lift305 + absprice299 + lift299 + absprice248 + lift248 + absprice324 + lift324 + absprice275 + lift275 + absprice289 + lift289 + absprice309 + lift309 + absprice266 + lift266 + ( ( 0 + logregprice ) | store ) + ( ( 0 + logpriceindex ) | store ) + ( ( 0 + displayonly ) | store ) + ( ( 0 + logseason ) | store )

Sample 4: for ppgid :256 region: US outlet: GR
logvolume ~ 1 + logregprice + logpriceindex + featuredisplay + featureonly + logseason + lift265 + absprice248 + lift248 + absprice278 + lift278 + lift229 + lift238 + absprice230 + lift230 + lift281 + lift239 + absprice246 + lift246 + lift236 + absprice275 + lift275 + absprice274 + lift274 + lift33 + lift178 + lift308 + absprice95 + lift95 + absprice20 + lift20 + absprice300 + lift300 + absprice12 + lift12 + absprice305 + lift305 + absprice107 + lift107 + absprice199 + lift199 + absprice106 + lift106 + absprice7 + absprice3 + lift3 + lift189 + lift298 + absprice297 + lift297 + absprice303 + lift303 + absprice5 + lift5 + ( ( 0 + logregprice ) | store ) + ( ( 0 + logpriceindex ) | store ) + ( ( 0 + featuredisplay ) | store ) + ( ( 0 + featureonly ) | store ) + ( ( 0 + logseason ) | store )

max_stars_repo_name: werbac/jlib

max_issues_repo_name: werbac/jlib

max_forks_repo_name: werbac/jlib


text:
using ThreeBodyDecay
using Test

@testset "LS amplitude integer spin" begin
    tbs = ThreeBodySystem(2.0, 1.0, 1.5; m0=6.0)
    dpp = randomPoint(tbs)
    #
    dc = DecayChainLS(1, σ->1/(4.1^2-σ-0.1im);
        two_s=2,
        parity='-',
        Ps = ['+','+','+','+'],
        tbs=tbs)
    @test sum(reim(amplitude(dpp, dc)) .≈ 0.0) == 0
end

@testset "LS amplitude half-integer spin" begin
    tbs = ThreeBodySystem( 2.0, 1.0, 1.5; m0=6.0,
        two_js = ThreeBodySpins(1, 0, 0; two_h0=1))  # 1/2+ 0- 0- 1/2+
    σs = randomPoint(tbs.ms)
    dpp = DalitzPlotPoint(;σs=σs,two_λs=[1,0,0,1])
    #
    dc = DecayChainLS(3, σ->1/(4.1^2-σ-0.1im); two_s=3, parity='-', Ps=['+','-','-','+'], tbs=tbs)
    # @show amplitude(dpp, dc)
    @test sum(reim(amplitude(dpp, dc)) .≈ 0.0) == 0
    # testing something else?
end

max_stars_repo_name: mmikhasenko/ThreeBodyDecay

max_issues_repo_name: mmikhasenko/ThreeBodyDecay

max_forks_repo_name: mmikhasenko/ThreeBodyDecay


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    precompile(FastTransforms.pochhammer, (Float64, Int,))
end