text:
### Run Simulation ###
function initialize_sim_and_run(;m::MapData
				, routes_path::OrderedDict{Symbol, OrderedDict{Int64,Int64}}
				, routes_distances::OrderedDict{Symbol, OrderedDict{Int64,Float64}}
				, routes_types_both_dir::Dict{Symbol,Symbol}
				, orig_map_nodes_num::Int
				, when_to_run_people::Float64 = sim_params.when_to_run_people
				, when_to_run_wagons::Float64 = sim_params.when_to_run_wagons
				, TTC_freq::Int = sim_params.TTC_car_freq
				, max_load_subway::Int = sim_params.max_load_subway
				, max_load_streetcar::Int = sim_params.max_load_streetcar
				, max_load_gov_restr::Float64 =  sim_params.max_load_gov_restr
				, N_agents::Int =  sim_params.N_agents
				, sim_run_time_in_secs::Int = sim_params.sim_run_time_in_secs
				, serialize_finished_sim_object::Bool = false
				, p0::Float64 =  sim_params.p0
				, parset_id::Int64 = 1)
	"""
	Description:
		Initialization of the simulation with a given parameters. Next the simulation will be running for the specified number of seconds.
	Returns:
		Simulation objects with the results.
	"""

	println("By default parameters from `sim_params` object will be used")
	println("--------------------")
	simul = Simulation(m
						, orig_map_nodes_num
						, routes_path
						, routes_distances
						, routes_types_both_dir
						, N_agents
						, sim_run_time_in_secs
						, TTC_freq
						, max_load_subway
						, max_load_streetcar
						, max_load_gov_restr
						, p0)

	println("when_to_run_people: ",when_to_run_people)
	println("when_to_run_wagons: ",when_to_run_wagons)
	enqueue!.(Ref(simul.events), simul.agents, tuple.(when_to_run_people,rand(length(simul.agents))))
	enqueue!.(Ref(simul.events), simul.wagons, tuple.(when_to_run_wagons,rand(length(simul.wagons))))
	enqueue!(simul.events, simul.artificial_agent_stats, tuple(when_to_run_people, rand()))

	println("--------------------")
	println("Final simul object before run: \n", simul)
	println("--------------------")
	current_time = 0

	while sim_run_time_in_secs > current_time && length(simul.events) > 0
		event, t = dequeue_pair!(simul.events)
		step!(simul,  t,  event, orig_map_nodes_num)
		current_time = deepcopy(t[1])
		if simul.infected_agents_count[end] >= 0.95*N_agents
			println("Simulation has ended ealier, because 95% of agents are already infected. Current time: ",current_time)
			break
		end
	end

	if serialize_finished_sim_object == true
		path_simul = "./data_output/simulation_finished_parset_id_$parset_id.bin"
		serialize(path_simul, simul)
		println("Simulation object with results was saved in: ", path_simul)
	end

	println("Simulation has already ended. You can check the results.")
    return simul
end
### Run Simulation N times###
function initialize_sim_and_run_N_times_and_gather_results(;m::MapData
				, routes_path::OrderedDict{Symbol, OrderedDict{Int64,Int64}}
				, routes_distances::OrderedDict{Symbol, OrderedDict{Int64,Float64}}
				, routes_types_both_dir::Dict{Symbol,Symbol}
				, orig_map_nodes_num::Int
				, when_to_run_people::Float64 = sim_params.when_to_run_people
				, when_to_run_wagons::Float64 = sim_params.when_to_run_wagons
				, TTC_freq::Int = sim_params.TTC_car_freq
				, max_load_subway::Int = sim_params.max_load_subway
				, max_load_streetcar::Int = sim_params.max_load_streetcar
				, max_load_gov_restr::Float64 =  sim_params.max_load_gov_restr
				, N_agents::Int =  sim_params.N_agents
				, sim_run_time_in_secs::Int = sim_params.sim_run_time_in_secs
				, serialize_finished_sim_object::Bool = false
				, parset_id::Int64 = 1
				, p0::Float64 =  sim_params.p0
				, N_times_to_run::Int64 = 30)
	"""
	Description:
		Run the simulation with specified parameters N times and then gather related results into dict.
	Returns:
		Dict object with mean results of N runs.
	"""

	res_gathering = Dict{String,Any}()
	total_cnt = Vector{Int}[]
	ttc_cnt = Vector{Int}[]
	street_cnt = Vector{Int}[]

	prct_ag_used_TTC = Float64[]
	mean_time_spent_in_TTC = Float64[]

	max_load_achieved_streetcars = Float64[]
	max_load_achieved_subway = Float64[]

	all_TTC_trips_count_streetcars = Float64[]
	all_TTC_trips_count_subway = Float64[]

	max_pass_per_wagon_subway = Vector{Float64}[]
	max_pass_per_wagon_streetcars = Vector{Float64}[]

	println("Your simulation will be run ", N_times_to_run, " times")
	for i in 1:N_times_to_run

		println("Is is a run number $i for parameters set number $parset_id.")
		@time s = initialize_sim_and_run(m = m
						, routes_path = routes_path
						, routes_distances = routes_distances
						, routes_types_both_dir = routes_types_both_dir
						, orig_map_nodes_num = orig_map_nodes_num
						, when_to_run_people = when_to_run_people
						, when_to_run_wagons = when_to_run_wagons
						, TTC_freq = TTC_freq
						, max_load_subway = max_load_subway
						, max_load_streetcar = max_load_streetcar
						, max_load_gov_restr = max_load_gov_restr
						, N_agents = N_agents
						, sim_run_time_in_secs = sim_run_time_in_secs
						, serialize_finished_sim_object =serialize_finished_sim_object
						, p0 = p0
						, parset_id = parset_id);
		push!(total_cnt, s.infected_agents_count)
		push!(ttc_cnt, s.infected_agents_wagon)
		push!(street_cnt, s.infected_agents_street)

		time_spend_in_TTC, TTC_types_used, time_spend_in_TTC_comp_to_total_path_time = TTC_usage_by_agents(s,orig_map_nodes_num);
		push!(prct_ag_used_TTC, round((sum([1 for i in time_spend_in_TTC if i > 0]) / N_agents)*100))
		push!(mean_time_spent_in_TTC, mean([i for i in time_spend_in_TTC if i > 0]))

		loop_max_load_achieved_streetcar = [1 for (k,v) in s.max_passengers_per_TTC_car if (routes_types_both_dir[k[2]] == :streetcar) & (v >= max_load_streetcar)]
		loop_max_load_achieved_subway = [1 for (k,v) in s.max_passengers_per_TTC_car if (routes_types_both_dir[k[2]] == :subway) & (v >= max_load_subway)]
		push!(max_load_achieved_streetcars,sum(loop_max_load_achieved_streetcar))
		push!(max_load_achieved_subway,sum(loop_max_load_achieved_subway))

		loop_all_TTC_trips_count_streetcars = [1 for (k,v) in s.max_passengers_per_TTC_car if (routes_types_both_dir[k[2]] == :streetcar)]
		loop_all_TTC_trips_count_subway = [1 for (k,v) in s.max_passengers_per_TTC_car if (routes_types_both_dir[k[2]] == :subway)]
		push!(all_TTC_trips_count_streetcars,sum(loop_all_TTC_trips_count_streetcars))
		push!(all_TTC_trips_count_subway,sum(loop_all_TTC_trips_count_subway))

		loop_max_pass_per_wagon_streetcars = [v for (k,v) in s.max_passengers_per_TTC_car if (routes_types_both_dir[k[2]] == :streetcar)]
		loop_max_pass_per_wagon_subway = [v for (k,v) in s.max_passengers_per_TTC_car if (routes_types_both_dir[k[2]] == :subway)]
		push!(max_pass_per_wagon_streetcars, [mean(loop_max_pass_per_wagon_streetcars),std(loop_max_pass_per_wagon_streetcars),maximum(loop_max_pass_per_wagon_streetcars)])
		push!(max_pass_per_wagon_subway, [mean(loop_max_pass_per_wagon_subway),std(loop_max_pass_per_wagon_subway),maximum(loop_max_pass_per_wagon_subway)])

	end

	res_gathering["total_infected"] = (floor.(mean(total_cnt)), floor.(std(total_cnt)))
	res_gathering["TTC_infected"] = (floor.(mean(ttc_cnt)), floor.(std(ttc_cnt)))
	res_gathering["street_infected"] = (floor.(mean(street_cnt)), floor.(std(street_cnt)))

	res_gathering["prct_of_agents_used_TTC"] = (floor.(mean(prct_ag_used_TTC)), floor.(std(prct_ag_used_TTC)))
	res_gathering["mean_sec_spent_in_TTC_by_those_who_used"] = (floor.(mean(mean_time_spent_in_TTC)), floor.(std(mean_time_spent_in_TTC)))

	res_gathering["times_max_load_achieved_streetcars"] = (floor.(mean(max_load_achieved_streetcars)), floor.(std(max_load_achieved_streetcars)))
	res_gathering["times_max_load_achieved_subway"] = (floor.(mean(max_load_achieved_subway)), floor.(std(max_load_achieved_subway)))

	res_gathering["all_TTC_trips_count_streetcars"] = (floor.(mean(all_TTC_trips_count_streetcars)),floor.(std(all_TTC_trips_count_streetcars)))
	res_gathering["all_TTC_trips_count_subway"] = (floor.(mean(all_TTC_trips_count_subway)),floor.(std(all_TTC_trips_count_subway)))

	res_gathering["max_pass_per_wagon_subway"] = (floor.(mean(max_pass_per_wagon_subway)),floor.(maximum(max_pass_per_wagon_subway)))
	res_gathering["max_pass_per_wagon_streetcars"] = (floor.(mean(max_pass_per_wagon_streetcars)),floor.(maximum(max_pass_per_wagon_subway)))

	println("All runs completed!")
	return res_gathering
end

### Create parameters grid ###
function create_df_with_sim_params(;TTC_car_freq::Union{Vector{Int64},Int64}
                                    , max_load_gov_restr::Union{Vector{Float64},Float64}
                                    , when_to_run_people::Union{Vector{Float64},Float64}
                                    , when_to_run_wagons::Union{Vector{Float64},Float64}
                                    , sim_run_time_in_secs::Union{Vector{Int64},Int64}
                                    , N_agents::Union{Vector{Int64},Int64}
									, agents_speed_in_sec_per_m::Union{Vector{Float64},Float64}
									, max_load_subway::Union{Vector{Int64},Int64}
									, max_load_streetcar::Union{Vector{Int64},Int64}
									, p0::Union{Vector{Float64},Float64})
	"""
	Description:
		Make all possible combinations of given values.
	Returns:
		DataFrame with all possible combinations of given parameters where one row means one set of parameters for Simulation.
	"""
    df = DataFrame(Iterators.product(TTC_car_freq
                                        , max_load_gov_restr
                                        , when_to_run_people
                                        , when_to_run_wagons
                                        , sim_run_time_in_secs
                                        , N_agents
										, agents_speed_in_sec_per_m
										, max_load_subway
										, max_load_streetcar
										, p0))
    colnames = ["TTC_car_freq", "max_load_gov_restr", "when_to_run_people", "when_to_run_wagons", "sim_run_time_in_secs", "N_agents", "agents_speed_in_sec_per_m", "max_load_subway", "max_load_streetcar", "p0"]
    rename!(df,Symbol.(colnames))
	df[!,"parset_id"] = [i for i in 1:nrow(df)]
	df = df[:,["parset_id","TTC_car_freq", "max_load_gov_restr", "when_to_run_people", "when_to_run_wagons", "sim_run_time_in_secs", "N_agents", "agents_speed_in_sec_per_m", "max_load_subway", "max_load_streetcar", "p0"]]
    return df
end

println("Run sim functions were read")

max_stars_repo_name: NykPol/EpidemicInUrbanNetworkToronto

max_issues_repo_name: NykPol/EpidemicInUrbanNetworkToronto

max_forks_repo_name: NykPol/EpidemicInUrbanNetworkToronto


text:
export Sheet, ElectricResponseStyle, MagneticResponseStyle, Impedance, Admittance

"""
    Sheet

An abstract type to dispatch methods for 2D structures embedded in a homogenous
medium modelled with Generalized Sheet Transition Conditions (GSTCs).
See:
Kuester, Edward et al.
"Averaged transition conditions for electromagnetic fields at a metafilm"
https://doi.org/10.1109/TAP.2003.817560
"""
abstract type Sheet end

"""
    ResponseStyle

Based on `IndexStyle` in Base, provides an interface for user-defined types to
choose between a formulation of the same problem in terms of their preferred
material response parameters.
"""
abstract type ResponseStyle end

"""
    Impedance()

`ResponseStyle` to describe materials whose response is given by an impedance.
"""
struct Impedance <: ResponseStyle end

"""
    Admittance()

`ResponseStyle` to describe materials whose response is given by an admittance.
"""
struct Admittance <: ResponseStyle end

"""
    ElectricResponseStyle(sheet)
    ElectricResponseStyle(typeof(sheet))
"""
ElectricResponseStyle(sheet::Sheet) = ElectricResponseStyle(typeof(sheet))
ElectricResponseStyle(::Type{<:Sheet}) = Impedance()

"""
    MagneticResponseStyle(sheet)
    MagneticResponseStyle(typeof(sheet))
"""
MagneticResponseStyle(sheet::Sheet) = MagneticResponseStyle(typeof(sheet))
MagneticResponseStyle(::Type{<:Sheet}) = Admittance()

# Define default behavior of impedance/admittance parameters
for R in (:Z, :Y), s in (:ₑ, :ₘ), i in (:ˣ, :ʸ), j in (:ˣ, :ʸ)
    @eval begin
        export $(Symbol(R, s, i, j))
        $(Symbol(R, s, i, j))(::Sheet, x⃗) = zero(eltype(x⃗))
    end
end

"""
    _sMatrix(::FieldStyle, ::PlaneWaves, ::Sheet, ::UniformMedium, ::UniformMedium)
"""
_sMatrix(fs, pw::PlaneWaves{1}, sheet::Sheet, um₁, um₂) = _sMatrix(fs, PolarizationStyle(pw), pw, sheet, um₁, um₂)
_sMatrix(fs::EField, pol::TE, pw, sheet, um₁, um₂) = __sMatrix(fs, pol, pw, ElectricResponseStyle(sheet), MagneticResponseStyle(sheet), sheet, um₁, um₂)
_sMatrix(fs::HField, pol::TM, pw, sheet, um₁, um₂) = __sMatrix(fs, pol, pw, ElectricResponseStyle(sheet), MagneticResponseStyle(sheet), sheet, um₁, um₂)
function __sMatrix(fs, pol, pw, rsₑ, rsₘ, sheet, um₁, um₂)
    R̃ₑ = diagFT(_get_electric_response(pol, rsₑ, pw, sheet))
    R̃ₘ = diagFT(_get_magnetic_response(pol, rsₘ, pw, sheet))
    kz₁ = Diagonal(_get_kz(pw, um₁))
    kz₂ = Diagonal(_get_kz(pw, um₂))
    ωη₁ = pw.ω * _get_η(pol, um₁)
    ωη₂ = pw.ω * _get_η(pol, um₂)
    out, inc = _gstcMatrix(fs,
        _set_electric_terms(pol, rsₑ, R̃ₑ, kz₁, kz₂, ωη₁, ωη₂)...,
        _set_magnetic_terms(pol, rsₘ, R̃ₘ, kz₁, kz₂, ωη₁, ωη₂)...,
    )
    out\inc
end
function __sMatrix(fs, pol, pw, rsₑ, rsₘ, sheet, um::T, ::T) where T<:UniformMedium
    R̃ₑ = diagFT(_get_electric_response(pol, rsₑ, pw, sheet))
    R̃ₘ = diagFT(_get_magnetic_response(pol, rsₘ, pw, sheet))
    kz = Diagonal(_get_kz(pw, um))
    ωη = pw.ω * _get_η(pol, um)
    out, inc = _gstcMatrix(fs,
        _set_electric_terms(pol, rsₑ, R̃ₑ, kz, ωη)...,
        _set_magnetic_terms(pol, rsₘ, R̃ₘ, kz, ωη)...,
    )
    out\inc
end

_get_electric_response(::TM, ::Admittance, pw, sheet) = Yₑˣˣ.(Ref(sheet), Iterators.product(pw.x⃗...))
_get_electric_response(::TM, ::Impedance,  pw, sheet) = Zₑˣˣ.(Ref(sheet), Iterators.product(pw.x⃗...))
_get_magnetic_response(::TM, ::Admittance, pw, sheet) = Yₘʸʸ.(Ref(sheet), Iterators.product(pw.x⃗...))
_get_magnetic_response(::TM, ::Impedance,  pw, sheet) = Zₘʸʸ.(Ref(sheet), Iterators.product(pw.x⃗...))

_get_electric_response(::TE, ::Admittance, pw, sheet) = Yₑʸʸ.(Ref(sheet), Iterators.product(pw.x⃗...))
_get_electric_response(::TE, ::Impedance,  pw, sheet) = Zₑʸʸ.(Ref(sheet), Iterators.product(pw.x⃗...))
_get_magnetic_response(::TE, ::Admittance, pw, sheet) = Yₘˣˣ.(Ref(sheet), Iterators.product(pw.x⃗...))
_get_magnetic_response(::TE, ::Impedance,  pw, sheet) = Zₘˣˣ.(Ref(sheet), Iterators.product(pw.x⃗...))

"""
    diagFT(A::AbstractArray{<:Number, ndim}) where ndim

Performs a Fourier transform of a diagonal operator/scalar field over an Nd grid
"""
function diagFT(A::AbstractArray{<:Number, ndim}) where ndim
    n⃗ = size(A)
    n = prod(n⃗)
    Ã = Array(reshape(Diagonal(reshape(A, n)), (n⃗..., n⃗...)))
    Ã = fft(ifft(Ã, 1:ndim), (ndim+1):2ndim)
    reshape(Ã, (n, n))
end

_get_η(::TE, ::UniformMedium{ϵ, μ} where ϵ) where μ = μ
_get_η(::TM, ::UniformMedium{ϵ, μ} where μ) where ϵ = ϵ

_set_electric_terms(::TM, ::Impedance,  R̃, kz₁, kz₂, ωη₁, ωη₂) = (-R̃, kz₁ / 2ωη₁, kz₂ / 2ωη₂)
_set_electric_terms(::TM, ::Admittance, R̃, kz₁, kz₂, ωη₁, ωη₂) = (-I, R̃*(kz₁ / 2ωη₁), R̃*(kz₂ / 2ωη₂))
_set_magnetic_terms(::TM, ::Impedance,  R̃, kz₁, kz₂, ωη₁, ωη₂) = (R̃*(kz₁ / 0.5ωη₁), R̃*(kz₂ / 0.5ωη₂), -I)
_set_magnetic_terms(::TM, ::Admittance, R̃, kz₁, kz₂, ωη₁, ωη₂) = (kz₁ / 0.5ωη₁, kz₂ / 0.5ωη₂, -R̃)

_set_electric_terms(::TM, ::Impedance,  R̃, kz, ωη) = (-R̃, kz / 2ωη)
_set_electric_terms(::TM, ::Admittance, R̃, kz, ωη) = (-I, R̃*(kz / 2ωη))
_set_magnetic_terms(::TM, ::Impedance,  R̃, kz, ωη) = (R̃*(kz / 0.5ωη), -I)
_set_magnetic_terms(::TM, ::Admittance, R̃, kz, ωη) = (kz / 0.5ωη, -R̃)

_set_electric_terms(::TE, args...) = _set_magnetic_terms(TM(), args...)
_set_magnetic_terms(::TE, args...) = _set_electric_terms(TM(), args...)

# compute E field for TM mode by using TM relations to convert H to E
function _sMatrix(::EField, pol::TM, pw, sheet, um₁, um₂)
    S = _sMatrix(HField(), pol, pw, sheet, um₁, um₂)
    η₁ = _get_η(pol, um₁)
    η₂ = _get_η(pol, um₂)
    kz₁ = _get_kz(pw, um₁)
    kz₂ = _get_kz(pw, um₂)
    HtoE = Diagonal(vcat(-kz₁ / pw.ω*η₁, kz₂ / pw.ω*η₂))
    HtoE * S * inv(-HtoE)
end
# compute H field for TE mode by using TE relations to convert E to H
function _sMatrix(::HField, pol::TE, pw, sheet, um₁, um₂)
    S = _sMatrix(EField(), pol, pw, sheet, um₁, um₂)
    η₁ = _get_η(pol, um₁)
    η₂ = _get_η(pol, um₂)
    kz₁ = _get_kz(pw, um₁)
    kz₂ = _get_kz(pw, um₂)
    EtoH = Diagonal(vcat(kz₁ / pw.ω*η₁, -kz₂ / pw.ω*η₂))
    EtoH * S * inv(-EtoH)
end

function _sBlockMatrix(fs, pw::PlaneWaves{N}, sheet::Sheet, um₁, um₂) where N
    n = length(pw)
    BlockMatrix(_sMatrix(fs, pw, sheet, um₁, um₂), [N*n, N*n], [N*n, N*n])
end

"""
    _sLinearMap(pw, sheet, pol)

Return a function that will compute the matrix-vector product with the
scattering matrix of a Sheet in a matrix-free fashion.
"""
function _sLinearMap(fs, pw::PlaneWaves{1}, sheet::Sheet, um₁, um₂)
    rsₑ = ElectricResponseStyle(sheet)
    rsₘ = MagneticResponseStyle(sheet)
    pol = PolarizationStyle(pw)
    n = length(pw)
    R̃ₘ = LinearMap(ifft ∘ (x -> Diagonal(_get_magnetic_response(pol, rsₘ, pw, sheet)) * x) ∘ fft, n)
    R̃ₑ = LinearMap(ifft ∘ (x -> Diagonal(_get_electric_response(pol, rsₑ, pw, sheet)) * x) ∘ fft, n)
    kz₁ = LinearMap(Diagonal(_get_kz(pw, um₁)))
    kz₂ = LinearMap(Diagonal(_get_kz(pw, um₂)))
    ωη₁ = pw.ω * _get_η(pol, um₁)
    ωη₂ = pw.ω * _get_η(pol, um₂)
    out, inc = _gstcLinearMap(fs,
        _set_electric_terms(pol, rsₑ, R̃ₑ, kz₁, kz₂, ωη₁, ωη₂)...,
        _set_magnetic_terms(pol, rsₘ, R̃ₘ, kz₁, kz₂, ωη₁, ωη₂)...,
    )
    LinearMap(x -> gmres(out, inc(x)), 2n)
end

"""
    smatrixBlockMap(pw, sheet, pol)

Return a function that will compute the matrix-vector product with the
scattering matrix of a Sheet in a matrix-free fashion.
"""
function _sBlockMap(fs, pw::PlaneWaves{1}, sheet::Sheet, um₁, um₂)
    N = length(pw)
    rsₑ = ElectricResponseStyle(sheet)
    rsₘ = MagneticResponseStyle(sheet)
    pol = PolarizationStyle(pw)
    _fft = LinearMap{ComplexF64}(fft, ifft, N)
    R̃ₘ = _fft' * LinearMap(Diagonal(_get_magnetic_response(pol, rsₘ, pw, sheet))) * _fft
    R̃ₑ = _fft' * LinearMap(Diagonal(_get_electric_response(pol, rsₑ, pw, sheet))) * _fft
    kz₁ = LinearMap(Diagonal(_get_kz(pw, um₁)))
    kz₂ = LinearMap(Diagonal(_get_kz(pw, um₂)))
    ωη₁ = pw.ω * _get_η(pol, um₁)
    ωη₂ = pw.ω * _get_η(pol, um₂)
    out, inc = _gstcMatrix(fs,
        _set_electric_terms(pol, rsₑ, R̃ₑ, kz₁, kz₂, ωη₁, ωη₂)...,
        _set_magnetic_terms(pol, rsₘ, R̃ₘ, kz₁, kz₂, ωη₁, ωη₂)...,
    )
    out\inc
end

function Base.:\(A, B::BlockMap)
    B₁₁, B₁₂, B₂₁, B₂₂ = B.maps
    N = size(B₁₁, 1)
    # Build the matrix inversion block-by-block
    invA_B₁₁ = LinearMap(x -> gmres(A, vcat(B₁₁ * x, B₂₁ * x))[1:N], N)
    invA_B₁₂ = LinearMap(x -> gmres(A, vcat(B₁₂ * x, B₂₂ * x))[1:N], N)
    invA_B₂₁ = LinearMap(x -> gmres(A, vcat(B₁₁ * x, B₂₁ * x))[(N+1):2N], N)
    invA_B₂₂ = LinearMap(x -> gmres(A, vcat(B₁₂ * x, B₂₂ * x))[(N+1):2N], N)
    [
        invA_B₁₁    invA_B₁₂;
        invA_B₂₁    invA_B₂₂;
    ]
end

_sMatrix(fs, pw::PlaneWaves{2}, sheet::Sheet, um₁, um₂) = __sMatrix(fs, pw, ElectricResponseStyle(sheet), MagneticResponseStyle(sheet), sheet, um₁, um₂)
function __sMatrix(fs, pw, rsₑ, rsₘ, sheet, um₁, um₂)
    Rₑ = _get_electric_response_components(rsₑ, pw, sheet)
    Rₘ = _get_magnetic_response_components(rsₘ, pw, sheet)
    R̃ₑ = [
        diagFT(vec(Rₑ.xx))   diagFT(vec(Rₑ.xy))
        diagFT(vec(Rₑ.yx))   diagFT(vec(Rₑ.yy))
    ]
    R̃ₘ = [
        diagFT(vec(Rₘ.xx))   diagFT(vec(Rₘ.xy))
        diagFT(vec(Rₘ.yx))   diagFT(vec(Rₘ.yy))
    ]
    ωϵ₁, ωμ₁ = pw.ω .* _get_ϵμ(um₁)
    ωϵ₂, ωμ₂ = pw.ω .* _get_ϵμ(um₂)
    kz₁ = Diagonal(repeat(vec(_get_kz(pw, um₁)), 2))
    kz₂ = Diagonal(repeat(vec(_get_kz(pw, um₂)), 2))
    _K = _get_K_components(pw)
    K = [
        Diagonal(vec(_K.xx)) Diagonal(vec(_K.xy))
        Diagonal(vec(_K.yx)) Diagonal(vec(_K.yy))
    ]
    Ω₁ = K - (ωϵ₁ * ωμ₁)*I
    Ω₂ = K - (ωϵ₂ * ωμ₂)*I
    n = length(pw)
    C = [
        (0I)(n) (-I)(n)
        (1I)(n) (0I)(n)
    ]
    out, inc = _gstcMatrix(fs,
        _set_electric_terms(fs, rsₑ, R̃ₑ, kz₁, kz₂, ωϵ₁, ωμ₁, ωϵ₂, ωμ₂, Ω₁, Ω₂, C)...,
        _set_magnetic_terms(fs, rsₘ, R̃ₘ, kz₁, kz₂, ωϵ₁, ωμ₁, ωϵ₂, ωμ₂, Ω₁, Ω₂, C)...,
    )
    out\inc
end
function __sMatrix(fs, pw, rsₑ, rsₘ, sheet, um::T, ::T) where T<:UniformMedium
    Rₑ = _get_electric_response_components(rsₑ, pw, sheet)
    Rₘ = _get_magnetic_response_components(rsₘ, pw, sheet)
    R̃ₑ = [
        diagFT(vec(Rₑ.xx))   diagFT(vec(Rₑ.xy))
        diagFT(vec(Rₑ.yx))   diagFT(vec(Rₑ.yy))
    ]
    R̃ₘ = [
        diagFT(vec(Rₘ.xx))   diagFT(vec(Rₘ.xy))
        diagFT(vec(Rₘ.yx))   diagFT(vec(Rₘ.yy))
    ]
    ωϵ, ωμ = pw.ω .* _get_ϵμ(um)
    kz = Diagonal(repeat(vec(_get_kz(pw, um)), 2))
    _K = _get_K_components(pw)
    K = [
        Diagonal(vec(_K.xx)) Diagonal(vec(_K.xy))
        Diagonal(vec(_K.yx)) Diagonal(vec(_K.yy))
    ]
    Ω = K - (ωϵ * ωμ)*I
    n = length(pw)
    C = [
        (0I)(n) (-I)(n)
        (1I)(n) (0I)(n)
    ]
    out, inc = _gstcMatrix(fs,
        _set_electric_terms(fs, rsₑ, R̃ₑ, kz, ωϵ, ωμ, Ω, C)...,
        _set_magnetic_terms(fs, rsₘ, R̃ₘ, kz, ωϵ, ωμ, Ω, C)...,
    )
    out\inc
end

for (kind, ksym) in ((:electric, :ₑ), (:magnetic, :ₘ)), (rs, rsym) in ((:Impedance, :Z), (:Admittance, :Y))
    @eval function $(Symbol(:_get_, kind, :_response_components))(::$rs, pw, sheet)
        (
            xx = $(Symbol(rsym, ksym, :ˣˣ)).(Ref(sheet), Iterators.product(pw.x⃗...)),
            xy = $(Symbol(rsym, ksym, :ˣʸ)).(Ref(sheet), Iterators.product(pw.x⃗...)),
            yx = $(Symbol(rsym, ksym, :ʸˣ)).(Ref(sheet), Iterators.product(pw.x⃗...)),
            yy = $(Symbol(rsym, ksym, :ʸʸ)).(Ref(sheet), Iterators.product(pw.x⃗...)),
        )
    end
end

function _get_K_components(pw::PlaneWaves{2})
    (
        xx = [k⃗[1] .* k⃗[1] for k⃗ in Iterators.product(pw.k⃗...)],
        xy = [k⃗[1] .* k⃗[2] for k⃗ in Iterators.product(pw.k⃗...)],
        yx = [k⃗[2] .* k⃗[1] for k⃗ in Iterators.product(pw.k⃗...)],
        yy = [k⃗[2] .* k⃗[2] for k⃗ in Iterators.product(pw.k⃗...)],
    )
end

_set_electric_terms(::HField, ::Impedance,  R̃, kz₁, kz₂, ωϵ₁, ωμ₁, ωϵ₂, ωμ₂, Ω₁, Ω₂, C) = (Ω₁ * Ω₂ * -C * R̃ * C, Ω₂ * (0.5ωμ₁ * kz₁), Ω₁ * (0.5ωμ₂ * kz₂))
_set_electric_terms(::HField, ::Admittance, R̃, kz₁, kz₂, ωϵ₁, ωμ₁, ωϵ₂, ωμ₂, Ω₁, Ω₂, C) = (I, -C * R̃ * C * Ω₁ * inv(2ωϵ₁ * kz₁), -C * R̃ * C * Ω₁ * inv(2ωϵ₂ * kz₂))
_set_magnetic_terms(::HField, ::Impedance,  R̃, kz₁, kz₂, ωϵ₁, ωμ₁, ωϵ₂, ωμ₂, Ω₁, Ω₂, C) = (R̃ * -C * Ω₁ * C * inv(2ωϵ₁ * kz₁), R̃ * -C * Ω₂ * C * inv(0.5ωϵ₂ * kz₂), I)
_set_magnetic_terms(::HField, ::Admittance, R̃, kz₁, kz₂, ωϵ₁, ωμ₁, ωϵ₂, ωμ₂, Ω₁, Ω₂, C) = (Ω₂ * (2ωμ₁ * kz₁), Ω₁ * (2ωμ₂ * kz₂), Ω₁ * Ω₂ * R̃)

_set_electric_terms(::HField, ::Impedance,  R̃, kz, ωϵ, ωμ, Ω, C) = (Ω * -C * R̃ * C, 0.5ωμ * kz)
_set_electric_terms(::HField, ::Admittance, R̃, kz, ωϵ, ωμ, Ω, C) = (I, -C * R̃ * C * Ω * inv(2ωϵ * kz))
_set_magnetic_terms(::HField, ::Impedance,  R̃, kz, ωϵ, ωμ, Ω, C) = (R̃ * -C * Ω * C * inv(0.5ωϵ * kz), I)
_set_magnetic_terms(::HField, ::Admittance, R̃, kz, ωϵ, ωμ, Ω, C) = (2ωμ * kz, Ω * R̃)

_set_electric_terms(::EField, rs, R̃, kz, ωϵ₁, ωμ₁, ωϵ₂, ωμ₂, Ω₁, Ω₂, C) = _set_magnetic_terms(HField(), rs, R̃, kz, ωμ₁, ωϵ₁, ωμ₂, ωϵ₂, Ω₁, Ω₂, C)
_set_magnetic_terms(::EField, rs, R̃, kz, ωϵ₁, ωμ₁, ωϵ₂, ωμ₂, Ω₁, Ω₂, C) = _set_electric_terms(HField(), rs, R̃, kz, ωμ₁, ωϵ₁, ωμ₂, ωϵ₂, Ω₁, Ω₂, C)

_set_electric_terms(::EField, rs, R̃, kz, ωϵ, ωμ, Ω, C) = _set_magnetic_terms(HField(), rs, R̃, kz, ωμ, ωϵ, Ω, C)
_set_magnetic_terms(::EField, rs, R̃, kz, ωϵ, ωμ, Ω, C) = _set_electric_terms(HField(), rs, R̃, kz, ωμ, ωϵ, Ω, C)

ext: jl


text:
@inline near(I::CartesianIndex,a=0) = (2I-2oneunit(I)):(2I-oneunit(I)-δ(a,I))

function restrictML(b::Array{Float64,m}) where m
    N = ntuple(i-> i==m ? m-1 : 1+size(b,i)÷2, m)
    a = zeros(N)
    restrictL!(a,b)
    Poisson(a)
end
@fastmath function restrictL!(a::Array{Float64,m},b::Array{Float64,m}) where m
    @inbounds for i ∈ 1:m-1, I ∈ inside(size(a)[1:m-1])
        a[I,i] = 0.5sum(@inbounds(b[J,i]) for J ∈ near(I,i))
    end
end

@fastmath restrict!(a::Array{Float64},b::Array{Float64}) = @inbounds @simd for I ∈ inside(a)
    a[I] = sum(@inbounds(b[J]) for J ∈ near(I))
end

prolongate!(a::Array{Float64},b::Array{Float64}) = @inbounds for I ∈ inside(b)
    @simd for J ∈ near(I)
        a[J] = b[I]
end;end

@inline divisible(N) = mod(N,2)==0 && N>4

struct MultiLevelPoisson{N,M} <: AbstractPoisson{N,M}
    levels :: Vector{Poisson{N,M}}
    n :: Vector{Int16}
    function MultiLevelPoisson(L::Array{Float64,n}) where n
        levels = [Poisson(L)]
        while all(size(levels[end].x) .|> divisible)
            push!(levels,restrictML(levels[end].L))
        end
        text = "MultiLevelPoisson requires size=a2ⁿ, where a<31, n>2"
        @assert (length(levels)>2 && all(size(levels[end].x).<31)) text
        new{n-1,n}(levels,[])
    end
end
function update!(ml::MultiLevelPoisson,L::Array)
    update!(ml.levels[1],L)
    for l ∈ 2:length(ml.levels)
        restrictL!(ml.levels[l].L,ml.levels[l-1].L)
        set_diag!(ml.levels[l])
    end
end

function Vcycle!(ml::MultiLevelPoisson;l=1)
    fine,coarse = ml.levels[l],ml.levels[l+1]
    # set up coarse level
    GS!(fine,it=0)
    restrict!(coarse.r,fine.r)
    fill!(coarse.x,0.)
    # solve coarse (with recursion if possible)
    l+1<length(ml.levels) && Vcycle!(ml,l=l+1)
    GS!(coarse,it=2)
    # correct fine
    prolongate!(fine.ϵ,coarse.x)
    increment!(fine)
end

mult(ml::MultiLevelPoisson,x) = mult(ml.levels[1],x)

function solver!(x::Array{Float64,n},ml::MultiLevelPoisson{n},b::Array{Float64,n};log=false,tol=1e-3,itmx=32) where n
    p = ml.levels[1]
    p.x .= x
    residual!(p,b); r₂ = L₂(p.r)
    log && (res = [r₂])
    nᵖ=0
    while r₂>tol && nᵖ<itmx
        Vcycle!(ml)
        GS!(p,it=2); r₂ = L₂(p.r)
        5tol>r₂>tol && (GS!(p,it=2); r₂ = L₂(p.r))
        log && push!(res,r₂)
        nᵖ+=1
    end
    x .= p.x
    push!(ml.n,nᵖ)
    log && return res
end


text:
module entity

import SynapseClient: @entityfunction, @standalonefunction, synapsecall, synapseclient

export
	split_entity_namespaces,
	is_container

@entityfunction split_entity_namespaces
@entityfunction is_container

end


text:
println("Compile libraries on main processor..."); flush(stdout)
include("imports.jl")
println("Finished loading libraries."); flush(stdout)

using Distributed
rmprocs(workers())
addprocs(32)
println("Compile libraries on workers ..."); flush(stdout)
@everywhere include("imports.jl")
println("Finished loading libraries."); flush(stdout)

Ks = collect(2:9)
snums = collect(1:4)
skewtmix = [true, false]

sims = dict_list(Dict(:K => Ks, :snum => snums, :skewtmix => skewtmix))

# Run simulations in parallel.
@time res = pmap(run, sims, on_error=identity)

# Print results status.
foreach(z -> println("$(z[1]) => $(z[2])"), zip(sims, res))

# Post process
pp_res = pmap(postprocess, sims)
# postprocess(Dict(:beta=>1, :K=>5, :snum=>1))

for snum in [1,2,3,4]
  imdir = mkpath(joinpath(Info.resultsdir_simstudy, simname, "img"))
  plot(size=plotsize)
  xlabel!("K")
  ylabel!("DIC")
  for stm in skewtmix
    dics = Float64[]
    for K in Ks
      info_path = joinpath(Info.resultsdir_simstudy, simname,
                           savename(Dict(:snum => snum, :K => K, :skewtmix => stm)),
                           "info.txt")
      model_info = open(f->read(f, String), info_path)
      _dic = parse_dic(model_info)
      append!(dics, _dic)
    end
    modelname = stm ? "Skew-t mixture" : "Normal mixture"
    plot!(Ks, dics, marker=:square, ms=8, label=modelname, legend=:topright)
  end
  savefig(joinpath(imdir, "dic_snum=$(snum).pdf"))
  closeall()
end

# Send results to S3
Util.s3sync(from="$(Info.resultsdir_simstudy)/$(simname)",
            to="$(Info.awsbucket_simstudy)/$(simname)",
            tags=`--exclude '*.nfs'`)


text:
"This module contains all inner parts."
module Internal

export Problem, Options, solve

using LinearAlgebra
using UnPack

"""
```
Problem(
    aᵢⱼ::Function
)
```

A function `aᵢⱼ(i, j)` for generating a symmetric matrix (upper triangle).
"""
struct Problem
    aᵢⱼ::Function
end

"""
`Options(n::Integer)`

Options to provide to the solver.

Includes:
- order `n` of the matrix;
- maximum number of iterations `iterₘₐₓ` for the Jacobi eigenvalue algorithm;
- error `ϵₚ` for power methods;
- error `ϵⱼ` for the Jacobi eigenvalue algorithm;
- function `get_x₀(A)` to get a vector of initial values.
"""
struct Options
    n::Integer
    iterₘₐₓ::Integer
    ϵₚ::Float64
    ϵⱼ::Float64
    get_x₀::Function
end

"""
The default set of options.

Sets:
- `n` to `10`;
- `iterₘₐₓ` to `1000000`;
- `ϵₚ` to `1e-16`;
- `ϵⱼ` to `1e-9`;
- `get_x₀(A)` to `(A) -> ones(size(A, 1))`.
"""
DefaultOptions = Options(10, 1000000, 1e-16, 1e-9, (A) -> ones(size(A, 1)))

include("Functions/solve.jl")

end


text:
# to simplify change boa(lo) to be constant
# boa(lo) = boac
function foa(γ_OA, γ′_OA, x1b, x1va, η)
    v_wr_g = (x1b - x1va) * sin(γ_OA)
    a = γ′_OA
    b = -v_wr_g
    c = boac
    fv =
        c *
        (lomax - lomin) *
        (
            a^2 * lomin^2 +
            a^2 * lomin * lomax +
            a^2 * lomax^2 +
            3 * a * b * lomin +
            3 * a * b * lomax +
            3 * b^2
        ) / 3
    0.25rhoh2o * cdoa * (1 + sign(γ′_OA)) * fv * η
end


text:
struct FractionalFactorial <: GSAMethod
end

"""
Code based on the theory presetned in:
Saltelli, A. (2008). Global sensitivity analysis: The primer. Chichester: Wiley, pp. 71-76.
"""

using LinearAlgebra

function _recursive_hadamard(k::Integer)
     """
     Generate a hadamard matrix via recursion.
     """
     # base case
     if k == 2
          return [1 1; 1 -1]
     else
          h = _recursive_hadamard(k ÷ 2)
     end

     return hcat(vcat(h, h), vcat(h, -h))
end

function _expanding_window_hadamard(k::Integer)
     """
     Generate hadamard matrix of size k using expanding window approach.
     """
     @assert ispow2(k)

     # intialize
     h = ones(Int64, k, k)
     h[2,2] = -1

     let
        bot_row = 2
        right_col = 2
        while bot_row < k

             @inbounds cur = h[1:bot_row, 1:right_col]

             new_bot = bot_row+1:bot_row*2
             new_right = right_col+1:right_col*2
             # update right
             @inbounds h[1:bot_row, new_right] = cur

             # update below
             @inbounds h[new_bot, 1:right_col] = cur

             # update diagonal
             @inbounds h[new_bot, new_right] = - cur

             # update window for each to 'copy' from.
             bot_row *= 2
             right_col *= 2
        end
     end
     return h
end

function _generate_hadamard(k::Integer)
     @assert ispow2(k)
     return _expanding_window_hadamard(k)
end

function generate_ff_design_matrix(num_parameters::Integer)
     """
     param: num_parameters
          The number of parameters to be sampled for.

     return:
          A 2-level fractional factorial design matrix of -1s and 1s with a design of resolution IV.

          The number of rows is 2*number of parameters and the number of columns is the number of parameters.

          See equation (2.31) in Saltelli, A. (2008) for details.
     """

     # If k is not a power of 2, get the next larger number that is
     # so the hadamard matrices can be computed.
     k = Integer(round(2^ceil(log2(num_parameters))))

     s = _generate_hadamard(k)
     design_matrix = vcat(s, -s) # design of resolution IV.

     return design_matrix
end

function generate_ff_sample_matrix(design_matrix::Array{Int64,2},  levels_list=nothing)
      """

      Convert the desgin matrix to a matrix whose rows can
      be used as inputs to a model by replacing low and high values for each parameter
      with the the low and high end of each parameters range, as specified in levels_list.

      param: design_matrix

            The matrix outputted by 'ff_design' i.e a 2-level fractional
            factorial design matrix of -1s and 1s with a design of resolution IV.
            The number of rows is 2*number of parameters and the
            number of columns is the number of parameters.

      param: levels_list (optional)
          A list of length equal to the number of parameters, where the elements of
          the list are 2-tuples.

          The ith tuple corresponds to the the low-value and the high-value of the ith parameter.
          i.e. levels_list[2][1] corresponds to the low-value of the 2nd parameter,
          and levels_list[2][2] corresponds to the high-value of the 2nd parameter.

          If no list is provided, the low value for each parameter
          is assumed to be 0 and the high value is assumed to be 1.

      returns:
          A matrix with the same dimensions as the design matrix, but with the '-1's replaced with
          the low value specified for each parameter in level_list, and similar the '1's are replaced
          with the high values.

          Each row corresponds to a single input vector to be used as input to a model/function.
      """
      sample_matrix = copy(design_matrix)
      if !isnothing(levels_list)
           for (col_index, (low_value, high_value)) in enumerate(levels_list)
                design_col = @view design_matrix[:,col_index]
                sample_col = @view sample_matrix[:,col_index]

                sample_col[design_col .== -1] .= low_value
                sample_col[design_col .== 1] .= high_value
            end
      end
      return sample_matrix
end

function run_model(sample_matrix::AbstractArray, model)
      """

      param: sample_matrix
            A matrix where each row corresponds to a input vector for "model"
      param: model
            The function/model to be executed on each row of the sample_matrix.
            Must return a scalar.
      return:
            A vector containing the result of running each row of the
            sample_matrix through model.
      """

      N = size(sample_matrix)[1]
      y_out = zeros(N)
      Threads.@threads for i in 1:N
          y_out[i] = model(@view sample_matrix[i,:])
      end
      return y_out
end

function ff_main_effects(design_matrix::AbstractArray, response_values::AbstractArray)
      """
           Computes the main effect for each parameter

           param: design_matrix

                The matrix outputted by 'ff_design' i.e a 2-level fractional
                factorial design matrix of -1s and 1s with a design of resolution IV.
                The number of rows is 2*number of parameters and the
                number of columns is the number of parameters.

           param: response_values

                An array of values corresponding to the result of running the
                model on each row of a sample matrix generated by the design_matrix

           return: main_effects

               A vector of the main effects of each parameter on the output.

            See Saltelli, A. equation 2.32 for details.
      """

      num_rows, num_cols = size(design_matrix)
      @assert num_rows == length(response_values) "Number of rows in design matrix must match number of responses"

      main_effects = zeros(num_cols)
      for column in 1:num_cols
            contrast_vector = @view design_matrix[:,column]
            main_effects[column] = dot(response_values, contrast_vector)
      end
      return main_effects ./ num_rows
end

function gsa(f, method::FractionalFactorial; num_params, N, p_range = nothing, kwargs...)
    design_matrix = generate_ff_design_matrix(num_params)
    sample_matrix = generate_ff_sample_matrix(design_matrix, p_range)

    response_vec = run_model(sample_matrix, f)

    main_effects = ff_main_effects(design_matrix, response_vec)
    return main_effects, main_effects.^2
end


text:
export iauAtciqz
"""
Quick ICRS to CIRS transformation, given precomputed star-
independent astrometry parameters, and assuming zero parallax and
proper motion.

Use of this function is appropriate when efficiency is important and
where many star positions are to be transformed for one date.  The
star-independent parameters can be obtained by calling one of the
functions iauApci[13], iauApcg[13], iauApco[13] or iauApcs[13].

The corresponding function for the case of non-zero parallax and
proper motion is iauAtciq.

This function is part of the International Astronomical Union's
SOFA (Standards of Fundamental Astronomy) software collection.

Status:  support function.

Given:
   rc,dc  double     ICRS astrometric RA,Dec (radians)
   astrom iauASTROM* star-independent astrometry parameters:
      pmt    double       PM time interval (SSB, Julian years)
      eb     double[3]    SSB to observer (vector, au)
      eh     double[3]    Sun to observer (unit vector)
      em     double       distance from Sun to observer (au)
      v      double[3]    barycentric observer velocity (vector, c)
      bm1    double       sqrt(1-|v|^2): reciprocal of Lorenz factor
      bpn    double[3][3] bias-precession-nutation matrix
      along  double       longitude + s' (radians)
      xpl    double       polar motion xp wrt local meridian (radians)
      ypl    double       polar motion yp wrt local meridian (radians)
      sphi   double       sine of geodetic latitude
      cphi   double       cosine of geodetic latitude
      diurab double       magnitude of diurnal aberration vector
      eral   double       "local" Earth rotation angle (radians)
      refa   double       refraction constant A (radians)
      refb   double       refraction constant B (radians)

Returned:
   ri,di  double     CIRS RA,Dec (radians)

Note:

   All the vectors are with respect to BCRS axes.

References:

   Urban, S. & Seidelmann, P. K. (eds), Explanatory Supplement to
   the Astronomical Almanac, 3rd ed., University Science Books
   (2013).

   Klioner, Sergei A., "A practical relativistic model for micro-
   arcsecond astrometry in space", Astr. J. 125, 1580-1597 (2003).

Called:
   iauS2c       spherical coordinates to unit vector
   iauLdsun     light deflection due to Sun
   iauAb        stellar aberration
   iauRxp       product of r-matrix and p-vector
   iauC2s       p-vector to spherical
   iauAnp       normalize angle into range +/- pi

This revision:   2013 October 9

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""
function iauAtciqz(rc::Real, dc::Real,  astrom::iauASTROM)

   # Allocate return value
   ref_astrom = Ref{iauASTROM}(astrom)
   ref_ri     = Ref{Float64}(0.0)
   ref_di     = Ref{Float64}(0.0)

   ccall((:iauAtciqz, libsofa_c), Cvoid,
            (Cdouble, Cdouble, Ref{iauASTROM},
            Ref{Cdouble}, Ref{Cdouble}),
            convert(Float64, rc), convert(Float64, dc),
            ref_astrom, ref_ri, ref_di)

   return ref_ri[], ref_di[]
end


text:
using Venti, Test
# implement your tests here
@test 1 == 1


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

module ApplicationGatewayClient

using Random
using Dates
using HTTP
using Swagger
import Swagger: field_name, property_type, hasproperty, validate_property, SwaggerApi, SwaggerModel
import Base: convert, propertynames

include("modelincludes.jl")

include("api_ApplicationGatewaysApi.jl")

# export models
export convert, ApplicationGatewayAuthenticationCertificatePropertiesFormat
export convert, ApplicationGatewayAvailableWafRuleSetsResult
export convert, ApplicationGatewayBackendAddress
export convert, ApplicationGatewayBackendAddressPoolPropertiesFormat
export convert, ApplicationGatewayBackendHealth
export convert, ApplicationGatewayBackendHealthHttpSettings
export convert, ApplicationGatewayBackendHealthPool
export convert, ApplicationGatewayBackendHealthServer
export convert, ApplicationGatewayBackendHttpSettingsPropertiesFormat
export convert, ApplicationGatewayConnectionDraining
export convert, ApplicationGatewayFirewallDisabledRuleGroup
export convert, ApplicationGatewayFirewallRule
export convert, ApplicationGatewayFirewallRuleGroup
export convert, ApplicationGatewayFirewallRuleSetPropertiesFormat
export convert, ApplicationGatewayFrontendIPConfigurationPropertiesFormat
export convert, ApplicationGatewayFrontendPortPropertiesFormat
export convert, ApplicationGatewayHttpListenerPropertiesFormat
export convert, ApplicationGatewayIPConfigurationPropertiesFormat
export convert, ApplicationGatewayListResult
export convert, ApplicationGatewayPathRulePropertiesFormat
export convert, ApplicationGatewayProbePropertiesFormat
export convert, ApplicationGatewayPropertiesFormat
export convert, ApplicationGatewayRequestRoutingRulePropertiesFormat
export convert, ApplicationGatewaySku
export convert, ApplicationGatewaySslCertificatePropertiesFormat
export convert, ApplicationGatewaySslPolicy
export convert, ApplicationGatewayUrlPathMapPropertiesFormat
export convert, ApplicationGatewayWebApplicationFirewallConfiguration
export convert, Resource
export convert, SubResource
export convert, ApplicationGateway
export convert, ApplicationGatewayAuthenticationCertificate
export convert, ApplicationGatewayBackendAddressPool
export convert, ApplicationGatewayBackendHttpSettings
export convert, ApplicationGatewayFirewallRuleSet
export convert, ApplicationGatewayFrontendIPConfiguration
export convert, ApplicationGatewayFrontendPort
export convert, ApplicationGatewayHttpListener
export convert, ApplicationGatewayIPConfiguration
export convert, ApplicationGatewayPathRule
export convert, ApplicationGatewayProbe
export convert, ApplicationGatewayRequestRoutingRule
export convert, ApplicationGatewaySslCertificate
export convert, ApplicationGatewayUrlPathMap
export convert, NetworkInterfaceIPConfiguration

# export operations
export convert, ApplicationGatewaysApi

export check_required, field_name, property_type, hasproperty, propertynames, validate_property, convert

end


text:
using JULES.Operators

using Oceananigans.BoundaryConditions

using Oceananigans.Fields: datatuple
using Oceananigans.TimeSteppers: tick!

import Oceananigans.TimeSteppers: time_step!
import Oceananigans.Simulations: ab2_or_rk3_time_step!

ab2_or_rk3_time_step!(model::CompressibleModel, Δt; euler) = time_step!(model, Δt)

function time_step!(model::CompressibleModel, Δt)
    arch = model.architecture
    total_density  = model.total_density
    momenta = model.momenta
    tracers = model.tracers
    diffusivities  = model.diffusivities

    slow_source_terms   = model.time_stepper.slow_source_terms
    fast_source_terms   = model.time_stepper.fast_source_terms
    intermediate_fields = model.time_stepper.intermediate_fields

    momenta_names = propertynames(momenta)
    tracers_names = propertynames(tracers)

    intermediate_momenta_fields = [getproperty(intermediate_fields, ρu)         for ρu in momenta_names]
    intermediate_tracers_fields = [getproperty(intermediate_fields.tracers, ρc) for ρc in tracers_names]

    intermediate_momenta = NamedTuple{momenta_names}(intermediate_momenta_fields)
    intermediate_tracers = NamedTuple{tracers_names}(intermediate_tracers_fields)

    first_stage_Δt  = Δt / 3
    second_stage_Δt = Δt / 2
    third_stage_Δt  = Δt

    #####
    ##### Compute slow source terms
    #####

    density_update_event =
        launch!(arch, model.grid, :xyz, update_total_density!,
                datatuple(total_density), model.grid, model.gases, datatuple(tracers),
                dependencies=Event(device(arch)))

    wait(device(arch), density_update_event)

    fill_halo_regions!(merge((Σρ=total_density,), momenta, tracers), model.architecture, model.clock, nothing)
    fill_halo_regions!(momenta.ρw, model.architecture, model.clock, nothing)
    fill_halo_regions!(intermediate_momenta.ρw, model.architecture, model.clock, nothing)

    compute_slow_source_terms!(
        slow_source_terms, arch, model.grid, model.thermodynamic_variable, model.gases, model.gravity,
        model.coriolis, model.closure, total_density, momenta, tracers, diffusivities, model.forcing, model.clock)

    fill_halo_regions!(slow_source_terms.ρw, model.architecture, model.clock, nothing)

    #####
    ##### Stage 1
    #####

    density_update_event =
        launch!(arch, model.grid, :xyz, update_total_density!,
                datatuple(total_density), model.grid, model.gases, datatuple(tracers),
                dependencies=Event(device(arch)))

    wait(device(arch), density_update_event)

    fill_halo_regions!(merge((Σρ=total_density,), momenta, tracers), model.architecture, model.clock, nothing)
    fill_halo_regions!(momenta.ρw, model.architecture, model.clock, nothing)
    fill_halo_regions!(intermediate_momenta.ρw, model.architecture, model.clock, nothing)

    compute_fast_source_terms!(
        fast_source_terms, arch, model.grid, model.thermodynamic_variable, model.gases, model.gravity,
        model.advection, total_density, momenta, tracers, slow_source_terms)

    calculate_boundary_tendency_contributions!(fast_source_terms, arch, momenta, tracers, model.clock, nothing)

    advance_state_variables!(intermediate_fields, arch, model.grid, momenta, tracers, fast_source_terms, Δt=first_stage_Δt)

    tick!(model.clock, 0, stage=true)

    #####
    ##### Stage 2
    #####

    density_update_event =
        launch!(arch, model.grid, :xyz, update_total_density!,
                datatuple(total_density), model.grid, model.gases, datatuple(intermediate_tracers),
                dependencies=Event(device(arch)))

    wait(device(arch), density_update_event)

    fill_halo_regions!(merge((Σρ=total_density,), intermediate_momenta, intermediate_tracers), model.architecture, model.clock, nothing)
    fill_halo_regions!(momenta.ρw, model.architecture, model.clock, nothing)
    fill_halo_regions!(intermediate_momenta.ρw, model.architecture, model.clock, nothing)

    compute_fast_source_terms!(
        fast_source_terms, arch, model.grid, model.thermodynamic_variable, model.gases, model.gravity,
        model.advection, total_density, intermediate_momenta, intermediate_tracers, slow_source_terms)

    calculate_boundary_tendency_contributions!(fast_source_terms, arch, intermediate_momenta,
                                               intermediate_tracers, model.clock, nothing)

    advance_state_variables!(intermediate_fields, arch, model.grid, momenta, tracers, fast_source_terms, Δt=second_stage_Δt)

    tick!(model.clock, 0, stage=true)

    #####
    ##### Stage 3
    #####

    density_update_event =
        launch!(arch, model.grid, :xyz, update_total_density!,
                datatuple(total_density), model.grid, model.gases, datatuple(intermediate_tracers),
                dependencies=Event(device(arch)))

    wait(device(arch), density_update_event)

    fill_halo_regions!(merge((Σρ=total_density,), intermediate_momenta, intermediate_tracers), model.architecture, model.clock, nothing)
    fill_halo_regions!(momenta.ρw, model.architecture, model.clock, nothing)
    fill_halo_regions!(intermediate_momenta.ρw, model.architecture, model.clock, nothing)

    compute_fast_source_terms!(
        fast_source_terms, arch, model.grid, model.thermodynamic_variable, model.gases, model.gravity,
        model.advection, total_density, intermediate_momenta, intermediate_tracers, slow_source_terms)

    calculate_boundary_tendency_contributions!(fast_source_terms, arch, intermediate_momenta,
                                               intermediate_tracers, model.clock, nothing)

    state_variables = (momenta..., tracers = tracers)
    advance_state_variables!(state_variables, arch, model.grid, momenta, tracers, fast_source_terms, Δt=third_stage_Δt)

    tick!(model.clock, Δt)

    return nothing
end


text:
# GFF3 File Format
# ================

module GFF3

import Automa
import Automa.RegExp: @re_str
import BGZFStreams
import BioCore.Exceptions: missingerror
import BioSequences
import BufferedStreams
import GenomicFeatures: GenomicFeatures, Interval, IntervalCollection
import URIParser
using BioCore

include("record.jl")
include("reader.jl")
include("writer.jl")

end


text:
#=
Copyright (c) 2016 Saurav Sachidanand

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in
all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
THE SOFTWARE.
=#

function type_from_fields(::Type{T}, fields) where T
    ccall(:jl_new_structv, Any, (Any, Ptr{Cvoid}, UInt32), T, fields, length(fields))
end

function peekbyte(io::IO)
    mark(io)
    byte = read(io, UInt8)
    reset(io)
    return byte
end

struct UndefIter{IO, F}
    f::F
    io::IO
end
Base.IteratorSize(::Type{<: UndefIter}) = Base.SizeUnknown()

function Base.iterate(x::UndefIter, state = nothing)
    peekbyte(x.io) == BREAK_INDEF && return nothing
    return x.f(x.io), nothing
end

function decode_ntimes(f, io::IO)
    first_byte = peekbyte(io)
    if (first_byte & ADDNTL_INFO_MASK) == ADDNTL_INFO_INDEF
        skip(io, 1) # skip first byte
        return UndefIter(f, io)
    else
        return (f(io) for i in 1:decode_unsigned(io))
    end
end

function decode_unsigned(io::IO)
    addntl_info = read(io, UInt8) & ADDNTL_INFO_MASK
    if addntl_info < SINGLE_BYTE_UINT_PLUS_ONE
        return addntl_info
    elseif addntl_info == ADDNTL_INFO_UINT8
        return bswap(read(io, UInt8))
    elseif addntl_info == ADDNTL_INFO_UINT16
        return bswap(read(io, UInt16))
    elseif addntl_info == ADDNTL_INFO_UINT32
        return bswap(read(io, UInt32))
    elseif addntl_info == ADDNTL_INFO_UINT64
        return bswap(read(io, UInt64))
    else
        error("Unknown Int type")
    end
end

decode_internal(io::IO, ::Val{TYPE_0}) = decode_unsigned(io)

function decode_internal(io::IO, ::Val{TYPE_1})
    data = signed(decode_unsigned(io))
    if (i = Int128(data) + one(data)) > typemax(Int64)
        return -i
    else
        return -(data + one(data))
    end
end

"""
Decode Byte Array
"""
function decode_internal(io::IO, ::Val{TYPE_2})
    if (peekbyte(io) & ADDNTL_INFO_MASK) == ADDNTL_INFO_INDEF
        skip(io, 1)
        result = IOBuffer()
        while peekbyte(io) !== BREAK_INDEF
            write(result, decode_internal(io))
        end
        return take!(result)
    else
        return read(io, decode_unsigned(io))
    end
end

"""
Decode String
"""
function decode_internal(io::IO, ::Val{TYPE_3})
    if (peekbyte(io) & ADDNTL_INFO_MASK) == ADDNTL_INFO_INDEF
        skip(io, 1)
        result = IOBuffer()
        while peekbyte(io) !== BREAK_INDEF
            write(result, decode_internal(io))
        end
        return String(take!(result))
    else
        return String(read(io, decode_unsigned(io)))
    end
end

"""
Decode Vector of arbitrary elements
"""
function decode_internal(io::IO, ::Val{TYPE_4})
    return map(identity, decode_ntimes(decode_internal, io))
end

"""
Decode Dict
"""
function decode_internal(io::IO, ::Val{TYPE_5})
    return Dict(decode_ntimes(io) do io
        decode_internal(io) => decode_internal(io)
    end)
end

"""
Decode Tagged type
"""
function decode_internal(io::IO, ::Val{TYPE_6})
    tag = decode_unsigned(io)
    data = decode_internal(io)
    if tag in (POS_BIG_INT_TAG, NEG_BIG_INT_TAG)
        big_int = parse(
            BigInt, bytes2hex(data), base = HEX_BASE
        )
        if tag == NEG_BIG_INT_TAG
            big_int = -(big_int + 1)
        end
        return big_int
    end

    if tag == CUSTOM_LANGUAGE_TYPE # Type Tag
        name = data[1]
        object_serialized = data[2]
        if startswith(name, "Julia/") # Julia Type
            return deserialize(IOBuffer(object_serialized))
        end
    end
    # TODO implement other common tags!
    return Tag(tag, data)
end

function decode_internal(io::IO, ::Val{TYPE_7})
    first_byte = read(io, UInt8)
    addntl_info = first_byte & ADDNTL_INFO_MASK
    if addntl_info < SINGLE_BYTE_SIMPLE_PLUS_ONE + 1
        simple_val = if addntl_info < SINGLE_BYTE_SIMPLE_PLUS_ONE
            addntl_info
        else
            read(io, UInt8)
        end
        if simple_val == SIMPLE_FALSE
            return false
        elseif simple_val == SIMPLE_TRUE
            return true
        elseif simple_val == SIMPLE_NULL
            return nothing
        elseif simple_val == SIMPLE_UNDEF
            return Undefined()
        else
            return Simple(simple_val)
        end
    else
        if addntl_info == ADDNTL_INFO_FLOAT64
            return reinterpret(Float64, ntoh(read(io, UInt64)))
        elseif addntl_info == ADDNTL_INFO_FLOAT32
            return reinterpret(Float32, ntoh(read(io, UInt32)))
        elseif addntl_info == ADDNTL_INFO_FLOAT16
            return reinterpret(Float16, ntoh(read(io, UInt16)))
        else
            error("Unsupported Float Type!")
        end
    end
end

function decode_internal(io::IO)
    # leave startbyte in io
    first_byte = peekbyte(io)
    typ = first_byte & TYPE_BITS_MASK
    typ == TYPE_0 && return decode_internal(io, Val(TYPE_0))
    typ == TYPE_1 && return decode_internal(io, Val(TYPE_1))
    typ == TYPE_2 && return decode_internal(io, Val(TYPE_2))
    typ == TYPE_3 && return decode_internal(io, Val(TYPE_3))
    typ == TYPE_4 && return decode_internal(io, Val(TYPE_4))
    typ == TYPE_5 && return decode_internal(io, Val(TYPE_5))
    typ == TYPE_6 && return decode_internal(io, Val(TYPE_6))
    typ == TYPE_7 && return decode_internal(io, Val(TYPE_7))
end


text:
push!(LOAD_PATH, pwd())
import .MasslistFunctions
H = MasslistFunctions.createCompound(H=1)
CO2 = MasslistFunctions.createCompound(C=1, O=2)
NH3 = MasslistFunctions.createCompound(N=1, H=3)
PYRIDINE = MasslistFunctions.createCompound(N=1,H=5,C=5)
ACETONE = MasslistFunctions.createCompound(C=3,H=6,O=1)
APINENE = MasslistFunctions.createCompound(C=10, H=16)
BCARY = MasslistFunctions.createCompound(C=15, H=24)
ISOPRENE = MasslistFunctions.createCompound(C=5, H=8)
HEXANONE =  MasslistFunctions.createCompound(C=6, H=12, O=1)
PINONALDEHYDE = MasslistFunctions.createCompound(C=10, H=16, O=2)
PINONALDEHYDEPAN = MasslistFunctions.createCompound(C=10, H=15, O=6, N=1)
PINONICACID = MasslistFunctions.createCompound(C=10, H=16, O=3) #C10H16O3
PINICACID = MasslistFunctions.createCompound(C=9, H=14, O=4) #C9H14O4
ACETIC = MasslistFunctions.createCompound(C=2, H=4, O=2)
ACETICFRAG = MasslistFunctions.createCompound(C=2, H=2, O=1)
NH3 = MasslistFunctions.createCompound(N=1, H=3)
DMA = MasslistFunctions.createCompound(C=2, H=7, N=1)
TMA = MasslistFunctions.createCompound(C=3, H=9, N=1)
TMB = MasslistFunctions.createCompound(C=9, H=12)

ACETONITRILE = MasslistFunctions.createCompound(C=2, H=3, N=1)

OrgNitratNO = MasslistFunctions.createCompound(C=10, H=15, O=5, N=1)
NORPINONALDEHYDE = MasslistFunctions.createCompound(C=9, H=14, O=2)
NORPINONALDEHYDEPAN = MasslistFunctions.createCompound(C=9, H=13, O=6, N=1)
H3O = MasslistFunctions.createCompound(H=2, O=1)
H3OH2O = MasslistFunctions.createCompound(H=4, O=2)
H3OH2OH2O = MasslistFunctions.createCompound(H=6, O=3)
H3OH2OH2OH2O = MasslistFunctions.createCompound(H=8, O=4)
O2 = MasslistFunctions.createCompound(O=2, Hplus = 0)

########## Sulphur ###############
DMS = MasslistFunctions.createCompound(C=2, H=6, S=1)
DMSO = MasslistFunctions.createCompound(C=2, H=6, S=1, O=1)
DMSO2 = MasslistFunctions.createCompound(C=2, H=6, S=1, O=2)
MSIA = MasslistFunctions.createCompound(C=1, H=4, S=1, O=2)
MSA = MasslistFunctions.createCompound(C=1, H=4, S=1, O=3)

########## BCARY Specific ########
#Fast
C15H22O2 = MasslistFunctions.createCompound(C=15, H=22, O=2)
C15H24O2 = MasslistFunctions.createCompound(C=15, H=24, O=2)
C15H24O3 = MasslistFunctions.createCompound(C=15, H=24, O=3)
C15H26O4 = MasslistFunctions.createCompound(C=15, H=26, O=4)
#Slow / Sticky
C14H20O3 = MasslistFunctions.createCompound(C=14, H=20, O=3)
C13H20O4 = MasslistFunctions.createCompound(C=13, H=20, O=4)
C14H22O4 = MasslistFunctions.createCompound(C=14, H=22, O=4)
C15H22O4 = MasslistFunctions.createCompound(C=15, H=22, O=4)
C15H24O4 = MasslistFunctions.createCompound(C=15, H=24, O=4)
C14H25NO4 = MasslistFunctions.createCompound(C=14, H=25, O=4, N=1)
C15H29NO6 = MasslistFunctions.createCompound(C=15, H=29, O=6, N=1)
#Nitrates
C15H23NO4 = MasslistFunctions.createCompound(C=15, H=23, O=4, N=1)
C13H19NO6 = MasslistFunctions.createCompound(C=13, H=19, O=6, N=1)
C15H25NO4 = MasslistFunctions.createCompound(C=15, H=25, O=4, N=1)
C15H21NO5 = MasslistFunctions.createCompound(C=15, H=21, O=5, N=1)
C15H23NO5 = MasslistFunctions.createCompound(C=15, H=23, O=5, N=1)
C15H23NO6 = MasslistFunctions.createCompound(C=15, H=23, O=6, N=1)
C15H17NO7 = MasslistFunctions.createCompound(C=15, H=17, O=7, N=1)

###### END BCARY ###############

###### NAPHTHA #################

NAPHTHA = MasslistFunctions.createCompound(C=10,H=8)

###### END NAPHTHA #############

compound2 = MasslistFunctions.createCompound(C=11, H=12, O=1)
compound3 = MasslistFunctions.createCompound(C=4, H=6, O=1)
compound4 = MasslistFunctions.createCompound(C=5, H=8, O=1)
compound5 = MasslistFunctions.createCompound(C=4, H=6, O=2)
compound6 = MasslistFunctions.createCompound(C=5, H=10, O=1)
compound7 = MasslistFunctions.createCompound(C=6, H=10, O=1)
compound8 = MasslistFunctions.createCompound(C=7, H=12, O=1)
compound9 = MasslistFunctions.createCompound(C=6, H=10, O=2)
compound10 = MasslistFunctions.createCompound(C=5, H=8, O=3)

# Dimers

dimer1 = MasslistFunctions.createCompound(C=19, H=28, O=9)
dimer2 = MasslistFunctions.createCompound(C=18, H=28, O=10)
dimer3 = MasslistFunctions.createCompound(C=20, H=32, O=9)
dimer4 = MasslistFunctions.createCompound(C=19, H=28, O=11)
dimer5 = MasslistFunctions.createCompound(C=19, H=26, O=19)

dimer6 = MasslistFunctions.createCompound(C=20, H=31, N=1, O=7)
dimer7 = MasslistFunctions.createCompound(C=19, H=29, N=1, O=8)
dimer8 = MasslistFunctions.createCompound(C=20, H=31, N=1, O=9)
dimer9 = MasslistFunctions.createCompound(C=20, H=31, N=1, O=13)


text:
module CutShortURL

# Dependencies
using Requests

# Function for shortening the urls, using tinyurl service.
@doc "Shortens the supplied URL" -> function cutshorturl(url::String)
 try
  robject = Requests.get("http://tinyurl.com/api-create.php?url=" * url)
  convert(String, robject.data)
 catch Exception
  "Error! Could not shorten $url"
 end
end

# Export functions to be accessed!
export
 cutshorturl
end


text:
function total_makespan(LS::LoadingSequence, CTS::Constants)
    makespan = 0
    flag = Array{Int,1}()
    while Set(collect(1:CTS.Q)) != Set(flag)
        for t in reverse(LS.order)
            for q = 1:CTS.Q
                if t.qc == q && !(q in flag)
                    push!(flag, q)
                    if t.start_time + 2*t.task.t > makespan
                        makespan = t.start_time + 2*t.task.t
                        break
                    end
                end
            end
        end
    end
    return(makespan)
end

#get QC moves and times
function get_qc_moves(LS::LoadingSequence, CTS::Constants)
    makespan = total_makespan(LS, CTS)
    QC_MOVES=Dict{Int, Array{NamedTuple{(:bay, :end_time, :start_time, :time),Tuple{Int64,Int64,Int64,Int64}}, 1}}()
    for q = 1:CTS.Q
        l=Array{NamedTuple{(:bay, :end_time, :start_time, :time),Tuple{Int64,Int64,Int64,Int64}}, 1}()
        push!(l, (bay=init_bay(q, CTS), end_time=0, start_time=0, time=0))
        for t in LS.order
            if t.qc == q
                push!(l, (bay=t.task.b, end_time=t.start_time + 2*t.task.t, start_time=t.start_time, time=t.task.t))
            end
        end
        push!(l,(bay=l[end].bay, end_time=makespan, start_time=makespan, time=0))
        QC_MOVES[q]=l
    end
    return(QC_MOVES)
end

#get QC travel times
function get_qc_travel(QC_MOVES::Dict{Int, Array{NamedTuple{(:bay, :end_time, :start_time, :time),Tuple{Int64,Int64,Int64,Int64}}, 1}}, CTS::Constants)
    QC_TRAVEL=Dict{Int, Array{NamedTuple{(:start_bay, :end_bay, :start_time, :end_time),Tuple{Int64,Int64,Int64,Int64}}, 1}}()
    for q = 1:CTS.Q
        l=Array{NamedTuple{(:start_bay, :end_bay, :start_time, :end_time),Tuple{Int64,Int64,Int64,Int64}}, 1}()
        for t = 1:length(QC_MOVES[q])-1
            if QC_MOVES[q][t].bay != QC_MOVES[q][t+1].bay
                start_time = QC_MOVES[q][t].end_time
                end_time = QC_MOVES[q][t+1].start_time
                start_bay = QC_MOVES[q][t].bay
                end_bay = QC_MOVES[q][t+1].bay
                if end_time - start_time <= abs(start_bay - end_bay)*CTS.tt
                    push!(l, (start_bay = start_bay, end_bay = end_bay, start_time = start_time, end_time = end_time))
                else
                    push!(l, (start_bay = start_bay, end_bay = start_bay, start_time = start_time, end_time=end_time - abs(start_bay - end_bay)*CTS.tt))
                    push!(l, (start_bay = start_bay, end_bay = end_bay, start_time = end_time - abs(start_bay - end_bay)*CTS.tt, end_time = end_time))
                end
            end
        end
        if length(QC_MOVES[q]) > 1 && QC_MOVES[q][end].bay == QC_MOVES[q][end-1].bay
            push!(l, (start_bay = QC_MOVES[q][end].bay, end_bay = QC_MOVES[q][end-1].bay, start_time = QC_MOVES[q][end].end_time, end_time=QC_MOVES[q][end-1].end_time))
        end
        QC_TRAVEL[q]=l
    end
    return(QC_TRAVEL)
end

function check_loaded_filled(LS::LoadingSequence, task::LTask)
    if task.p in LS.filled_pos
        #println("position already filled.")
        return(false)
    elseif task.c in LS.loaded_cont
        #println("container already loaded.")
        return(false)
    else
        return(true)
    end
end

#check position precedences
function check_prec(LS::LoadingSequence, task::LTask, prec::Dict{Int, Array})
    if length(prec[task.p]) > 0
        for req in prec[task.p]
            if !(req in LS.filled_pos)
                #print("trying to add position "*string(task.p)*", but still missing positions: ")
                #print(setdiff(precedences[task.p], LS.filled_pos))
                #print("\n")
                return(false)
            end
        end
    end
    return(true)
end

#check if all positions/containers are filled/loaded without repetitions
function check_tasks(LS::LoadingSequence, CTS::Constants)
    deleteat!(LS.filled_pos, findall(x->x==0, LS.filled_pos))
    deleteat!(LS.loaded_cont, findall(x->x==0, LS.loaded_cont))
    if LS.tasks_left == 0 && length(setdiff(collect(1:CTS.P), LS.filled_pos)) == 0 && length(setdiff(collect(1:CTS.C), LS.loaded_cont)) == 0
        return(true)
    else
        return(false)
    end
end

#check clearance to know if a QC can go to a certain bay
function check_clearance(QC::Array{QuayCrane, 1}, q::Int, target_bay::Int, CTS::Constants)
    if target_bay == QC[q].current_bay
        return(true)
    end

    m=0
    if q == 1
        for n = 1:CTS.Q-1
            if (QC[q+n].next_bay - target_bay) > CTS.delta && (QC[q+n].current_bay - target_bay) > CTS.delta
                if length(QC[q+n].task_buffer) > 0
                    for t_buffer in QC[q+n].task_buffer
                        if (t_buffer.b - target_bay) <= CTS.delta
                            return(false)
                        end
                    end
                end
                m+=1
            end
        end
        if m==CTS.Q-1
            return(true)
        end
    elseif q == CTS.Q
        for n = 1:CTS.Q-1
            if (target_bay - QC[q-1].next_bay) > CTS.delta && (target_bay - QC[q-1].current_bay) > CTS.delta
                if length(QC[q-1].task_buffer) > 0
                    for t_buffer in QC[q-1].task_buffer
                        if (target_bay - t_buffer.b) <= CTS.delta
                            return(false)
                        end
                    end
                end
                m+=1
            end
        end
        if m==CTS.Q-1
            return(true)
        end
    else
        if target_bay < QC[q].current_bay
            for n = 1:q-1
                if (target_bay - QC[q-1].next_bay)> CTS.delta && (target_bay - QC[q-1].current_bay) > CTS.delta
                    if length(QC[q-1].task_buffer) > 0
                        for t_buffer in QC[q-1].task_buffer
                            if (target_bay - t_buffer.b) <= CTS.delta
                                return(false)
                            end
                        end
                    end
                    m+=1
                end
            end
            if m==CTS.Q-1
                return(true)
            end
        elseif target_bay > QC[q].current_bay
            for n = q+1:CTS.Q
                if (QC[q+1].next_bay - target_bay) > CTS.delta && (QC[q+1].current_bay - target_bay) > CTS.delta
                    if length(QC[q+1].task_buffer) > 0
                        for t_buffer in QC[q+1].task_buffer
                            if (t_buffer.b - target_bay) <= CTS.delta
                                return(false)
                            end
                        end
                    end
                    m+=1
                end
            end
            if m==CTS.Q-1
                return(true)
            end
        end
    end
    return(false)
end

function check_double_clearance(QC::Array{QuayCrane, 1}, q::Int, target_bay::Int, CTS::Constants)
    if CTS.Q > 2
        if q <= 2 && target_bay > QC[q].current_bay
            if QC[q+1].status == "idle" && (QC[q+2].next_bay - target_bay) > CTS.delta + 2
                    return(true, q+1, target_bay + CTS.delta + 1, "idle")
            elseif QC[q+1].status == "moving" && (QC[q+1].next_bay - target_bay) > CTS.delta
                    return(true, q+1, target_bay + CTS.delta + 1, "moving")
            end
        elseif q >= CTS.Q - 1 && target_bay < QC[q].current_bay
            if QC[q-1].status == "idle" && (target_bay - QC[q-2].next_bay) > CTS.delta + 2
                    return(true, q-1, target_bay - CTS.delta - 1, "idle")
            elseif QC[q-1].status == "moving" && (target_bay - QC[q-1].next_bay) > CTS.delta
                    return(true, q-1, target_bay - CTS.delta - 1, "moving")
            end
        elseif 2 < q && q < CTS.Q - 1
            if target_bay < QC[q].current_bay
                if QC[q-1].status == "idle" && (target_bay - QC[q-2].current_bay) > CTS.delta + 2 && (target_bay - QC[q-2].next_bay) > CTS.delta + 2
                        return(true, q-1, target_bay - CTS.delta - 1, "idle")
                elseif QC[q-1].status == "moving" && (target_bay - QC[q-1].next_bay) > CTS.delta
                        return(true, q-1, target_bay - CTS.delta - 1, "moving")
                end
            elseif target_bay > QC[q].current_bay
                if QC[q+1].status == "idle" && (QC[q+2].next_bay - target_bay) > CTS.delta + 2
                        return(true, q+1, target_bay + CTS.delta + 1, "idle")
                elseif QC[q+1].status == "moving" && (QC[q+1].next_bay - target_bay) > CTS.delta
                        return(true, q+1, target_bay + CTS.delta + 1, "moving")
                end
            end
        end
    else
        if q == 1 && target_bay > QC[q].current_bay
            #if crane is idle, look if there is enough space to move it
            if QC[q+1].status == "idle" && target_bay < CTS.J - CTS.delta
                return(true, q+1, target_bay + CTS.delta + 1, "idle")
            #if crane is moving, look if it will move enough
            elseif QC[q+1].status == "moving" && (QC[q+1].next_bay - target_bay) > CTS.delta
                if target_bay < CTS.J - CTS.delta
                    return(true, q+1, target_bay + CTS.delta + 1, "moving")
                end
            end
        elseif q == CTS.Q && target_bay < QC[q].current_bay
            #if crane is idle, look if there is enough space to move it
            if QC[q-1].status == "idle" && target_bay > CTS.delta + 1
                return(true, q-1, target_bay - CTS.delta - 1, "idle")
            #if crane is moving, look if it will move enough
            elseif QC[q-1].status == "moving" && (target_bay - QC[q-1].next_bay) > CTS.delta
                if target_bay > 1 + CTS.delta
                    return(true, q-1, target_bay - CTS.delta - 1, "moving")
                end
            end
        end
    end
    return(false, 0, 0, "false")
end

function check_triple_clearance(QC::Array{QuayCrane, 1}, q::Int, target_bay::Int, CTS::Constants)
    if CTS.Q > 2
        if q == 1 && target_bay > QC[q].current_bay
            if QC[q+1].status == "idle" && QC[q+2].status == "idle" && (QC[q+3].next_bay - target_bay) > CTS.delta + 4
                    return(true, q+1, target_bay + CTS.delta + 1, q+2, target_bay + CTS.delta + 3)
            end
        elseif q == CTS.Q && target_bay < QC[q].current_bay
            if QC[q-1].status == "idle" && QC[q-2].status == "idle" && (target_bay - QC[q-3].next_bay) > CTS.delta + 4
                    return(true, q-1, target_bay - CTS.delta - 1, q-2, target_bay - CTS.delta - 3)
            end
        elseif q == 2 && target_bay > QC[q].current_bay
            if QC[q+1].status == "idle" && QC[q+2].status == "idle" && (CTS.J - target_bay) > CTS.delta + 4
                    return(true, q+1, target_bay + CTS.delta + 1, q+2, target_bay + CTS.delta + 3)
            end
        elseif q == 3 && target_bay < QC[q].current_bay
            if QC[q-1].status == "idle" && QC[q-2].status == "idle" && (target_bay - 1) > CTS.delta + 4
                    return(true, q-1, target_bay - CTS.delta - 1, q-2, target_bay - CTS.delta - 3)
            end
        end
    end
    return(false, 0, 0, 0, 0)
end

function is_crane_done(q::Int, QC::Array{QuayCrane,1}, CTS::Constants)
    work_load = 0
    for j in QC[q].available_bays
        work_load += remaining_bay_time(j, tasks_by_w, bj, LS, CTS)
    end
    if work_load == 0
        return(true)
    else
        return(false)
    end
end

#check if a solution is correct
function check_solution(prec::Dict{Int, Array}, LS::LoadingSequence, CTS::Constants)
    #check positions and containers
    if check_tasks(LS, CTS) == false
        # println("missing tasks")
        return(false)
    end

    #check precedences are followed
    for task in LS.order
        if task[1].p != 0
            if check_prec(LS, task[1], prec) == false
                println("wrong precedences")
                return(false)
            end
        end
    end

    #check the QC constraints
    QC_MOVES=get_qc_moves(LS, CTS)
    QC_TRAVEL=get_qc_travel(QC_MOVES, CTS)
    for clock = 1:total_makespan(LS, CTS)
        # println()
        # println("|||||||||||   "*string(clock)*"   ||||||||||||||")
        # println()
        current_cranes=Array{Tuple{Int, Number}, 1}()
        # println("LOADING CRANE")
        for (key, value) in QC_MOVES
            if clock == 1
                for i = 1: length(value)-1
                    #check loading + travel time
                    if value[i].end_time + travel_time(value[i].bay, value[i+1].bay, CTS) > value[i+1].start_time
                        # println("wrong task/travel time")
                        # println("QUAY CRANE: "*string(key))
                        # println(value[i])
                        # println(value[i+1])
                        # println("-------------- ")
                        return(false)
                    end
                    #check that cranes respect bays
                    if !(value[i].bay in subset_bay(CTS, key))
                        println("wrong bays")
                        return(false)
                    end
                end
                if !(value[length(value)].bay in subset_bay(CTS, key))
                    println("wrong bays")
                    return(false)
                end
            end

            for tuple in value
                if tuple.start_time <= clock && clock <= tuple.start_time + 2*tuple.time
                    push!(current_cranes, (key, tuple.bay))
                    # println("time: "*string(clock)*", crane: "*string(key)*", bay: "*string(tuple.bay))
                end
            end
        end

        # println("-----------------------------")
        # println("MOVING CRANE")
        for (key_travel, value_travel) in QC_TRAVEL
            for tuple in value_travel
                if tuple.start_time <= clock && clock <= tuple.end_time && tuple.start_bay != tuple.end_bay
                    current_bay = (clock-tuple.end_time)*(tuple.start_bay-tuple.end_bay)/(tuple.start_time-tuple.end_time) + tuple.end_bay
                    push!(current_cranes, (key_travel, current_bay))
                    # println("time: "*string(clock)*", crane: "*string(key_travel)*", bay: "*string(current_bay))
                end
            end
        end
        # println("-----------------------------")

        #check clearance
        sort(current_cranes, by = first)
        unique!(current_cranes)
        #println(current_cranes)

        for b = 1:length(current_cranes)-1
            if abs(current_cranes[b][2] - current_cranes[b+1][2]) <= CTS.delta
                # println(QC)
                # println("wrong clearance")
                return(false)
            end
        end
    end
    return(true)
end

#plot the solution
using Gadfly
using DataFrames
Gadfly.set_default_plot_size(25cm, 10cm)
function plot_solution(LS::LoadingSequence, makespan::Int, CTS::Constants)
    QC_MOVES=get_qc_moves(LS, CTS)
    QC_TRAVEL=get_qc_travel(QC_MOVES, CTS)

    list_moves=[]
    list_travel=[]
    for n=1:CTS.Q
        push!(list_moves, DataFrame(Bay = [x.bay for x in QC_MOVES[n]], Start = [x.start_time for x in QC_MOVES[n]], End = [x.end_time for x in QC_MOVES[n]], QuayCrane = map(i -> "QC - " * string(i), n)))
        push!(list_travel, DataFrame(start_bay = [x.start_bay for x in QC_TRAVEL[n]], end_bay = [x.end_bay for x in QC_TRAVEL[n]], start_time = [x.start_time for x in QC_TRAVEL[n]], end_time = [x.end_time for x in QC_TRAVEL[n]], QuayCrane = map(i -> "QC - " * string(i), n)))
    end
    df_moves=vcat(list_moves...)
    df_travel=vcat(list_travel...)

    #plot scale
    xsc  = Scale.x_continuous(minvalue=0, maxvalue=makespan)
    ysc  = Scale.y_continuous(minvalue=0, maxvalue=CTS.J+1)
    xtk = collect(0:60:makespan)
    push!(xtk, makespan)
    ytk = collect(1:CTS.J)
    #plot(df, x = :Start, xend = :Time, y = :Bay, yend = :Bay, color = :QuayCrane, Geom.segment(filled=false), Theme(line_width=1cm, major_label_font="CMU Serif",minor_label_font="CMU Serif",
          #     major_label_font_size=16pt,minor_label_font_size=14pt), xsc, ysc)

    layer1 = layer(df_moves, x = :Start, xend = :End, y = :Bay, yend = :Bay, color = :QuayCrane, Geom.segment(filled=false), Theme(alphas=[0.5], line_width=0.25cm, major_label_font="CMU Serif",minor_label_font="CMU Serif",
               major_label_font_size=16pt,minor_label_font_size=14pt))
    layer2 = layer(df_travel, x = :start_time, xend = :end_time, y = :start_bay, yend = :end_bay, color = :QuayCrane, Geom.segment(filled=false), Theme(line_style=[:dot] ,line_width=0.05cm, major_label_font="CMU Serif",minor_label_font="CMU Serif",
               major_label_font_size=16pt,minor_label_font_size=14pt))

    plot(layer1, layer2, xsc, ysc, Guide.xticks(ticks=xtk), Guide.yticks(ticks=ytk))

end

max_stars_repo_licenses: Zlib MIT

max_issues_repo_licenses: Zlib MIT

max_forks_repo_licenses: Zlib MIT


text:
using OrdinaryDiffEq, DiffEqSensitivity, Flux
using ComponentArrays, LinearAlgebra, Optimization, Test

const nknots = 10
const h = 1.0/(nknots+1)
x = range(0, step=h, length=nknots)
u0 = sin.(π*x)

@inline function f(du,u,p,t)
  du .= zero(eltype(u))
  u₃ = @view u[3:end]
  u₂ = @view u[2:end-1]
  u₁ = @view u[1:end-2]
  @. du[2:end-1] = p.k*((u₃ - 2*u₂ + u₁)/(h^2.0))
  nothing
end

p_true = ComponentArray(k=0.42)
jac_proto = Tridiagonal(similar(u0,nknots-1), similar(u0), similar(u0, nknots-1))
prob = ODEProblem(ODEFunction(f,jac_prototype=jac_proto), u0, (0.0,1.0), p_true)
@time sol_true = solve(prob, Rodas4P(), saveat=0.1)

function loss(p)
  _prob = remake(prob, p=p)
  sol = solve(_prob, Rodas4P(autodiff=false), saveat=0.1, sensealg=ForwardDiffSensitivity())
  sum((sol .- sol_true).^2)
end

p0 = ComponentArray(k=1.0)

optf = Optimization.OptimizationFunction((x,p) -> loss(x), Optimization.AutoZygote())
optprob = Optimization.OptimizationProblem(optf, p0)
res = Optimization.solve(optprob, ADAM(0.01), maxiters = 100)

@test res.u.k ≈ 0.42461977305259074 rtol=1e-1


text:
"""
    disjunction(m, args...)

Add disjunction macro.
"""
macro disjunction(m, args...)
    #get disjunction (pos_args) and keyword arguments
    pos_args, kw_args, _ = Containers._extract_kw_args(args)
    @assert length(pos_args) > 1 "At least 2 disjuncts must be included. If there is an empty disjunct, use `nothing`."

    #get kw_args and set defaults if missing
    reformulation = filter(i -> i.args[1] == :reformulation, kw_args)
    if !isempty(reformulation)
        reformulation = reformulation[1].args[2]
        reformulation_kind = eval(reformulation)
        @assert reformulation_kind in [:big_m, :hull] "Invalid reformulation method passed to keyword argument `:reformulation`. Valid options are :big_m (Big-M Reformulation) and :hull (Hull Reformulation)."
        if reformulation_kind == :big_m
            M = filter(i -> i.args[1] == :M, kw_args)
            param = !isempty(M) ? M[1].args[2] : :(missing)
        elseif reformulation_kind == :hull
            ϵ = filter(i -> i.args[1] == :ϵ, kw_args)
            param = !isempty(ϵ) ? ϵ[1].args[2] : :(1e-6)
        else
        end
    else
        throw(UndefKeywordError(:reformulation))
    end
    name = filter(i -> i.args[1] == :name, kw_args)
    if !isempty(name)
        name = name[1].args[2]
        disj_name = Symbol("disj_",eval(name))
    else
        name = :(Symbol("disj",gensym()))
        disj_name = eval(name)
    end

    #create constraints for each disjunction
    disj_names = [Symbol("$(disj_name)[$i]") for i in eachindex(pos_args)]
    disjunction = []
    for (d,dname) in zip(pos_args,disj_names)
        if Meta.isexpr(d, :tuple)
            for (j,di) in enumerate(d.args)
                i = findfirst(x -> x == d, pos_args)
                dname_j = Symbol("$(disj_name)[$i,$j]")
                d.args[j] = add_disjunction_constraint(m, di, dname_j)
            end
            push!(disjunction, d)
        else
            push!(disjunction, add_disjunction_constraint(m, d, dname))
        end
    end

    # #XOR constraint name
    # xor_con = Symbol("XOR($disj_name)")

    #build disjunction
    code = quote
        @assert !in($name, keys(object_dictionary($m))) "The disjunction name $name is already registered in the model. Specify new name."
        $m[$name] = @variable($m, [eachindex($disjunction)], Bin, base_name = string($name), lower_bound = 0, upper_bound = 1)
        # @constraint($m, $xor_con, sum($m[$name]) == 1)
        reformulate_disjunction($m, $(disjunction...); bin_var = $name, reformulation = $reformulation, param = $param)
    end

    return esc(code)
end

"""
    add_disjunction_constraint(m, d, dname)

Add disjunction constraint with name dname.
"""
function add_disjunction_constraint(m, d, dname)
    if Meta.isexpr(d, :block)
        d = quote
            try
                @constraints($m,$d)
            catch e
                if e isa ErrorException
                    @NLconstraints($m,$d)
                else
                    throw(e)
                end
            end
        end
    elseif Meta.isexpr(d, (:call, :comparison))
        d = quote
            try
                @constraint($m,$dname,$d)
            catch e
                if e isa ErrorException
                    @NLconstraint($m,$d)
                else
                    throw(e)
                end
            end
        end
    end

    return d
end

"""
    add_disjunction!(m::Model,disj...;reformulation::Symbol,M=missing,ϵ=1e-6,name=missing)

Add disjunction and reformulate.
"""
function add_disjunction!(m::Model,disj...;reformulation::Symbol,M=missing,ϵ=1e-6,name=missing)
    #run checks
    @assert reformulation in [:big_m, :hull] "Invalid reformulation method passed to keyword argument `:reformulation`. Valid options are :big_m (Big-M Reformulation) and :hull (Hull Reformulation)."
    @assert length(disj) > 1 "At least 2 disjuncts must be included. If there is an empty disjunct, use `nothing`."

    #get kw_args and set defaults if missing
    param = reformulation == :big_m ? M : ϵ
    bin_var = ismissing(name) ? Symbol("disj",gensym()) : name
    disj_name = ismissing(name) ? bin_var : Symbol("disj_",name)

    # #XOR constraint name
    # xor_con = "XOR($disj_name)"

    #apply reformulation
    @assert !in(bin_var, keys(object_dictionary(m))) "The disjunction name $bin_var is already registered in the model. Specify new name."
    #create indicator variable
    m[bin_var] = @variable(m, [eachindex(disj)], Bin, base_name = string(bin_var), lower_bound = 0, upper_bound = 1)
    # #add xor constraint on binary variable
    # m[Symbol(xor_con)] = @constraint(m, sum(m[bin_var]) == 1, base_name = xor_con)
    #reformulate disjunction
    reformulate_disjunction(m, disj...; bin_var, reformulation, param)
end

"""
    proposition(m, expr)

Add logical proposition macro.
"""
macro proposition(m, expr)
    #get args
    expr = QuoteNode(expr)
    code = :(DisjunctiveProgramming.to_cnf!($m, $expr))

    return esc(code)
end


text:
import Flux: gelu

@init @require CuArrays="3a865a2d-5b23-5a0f-bc46-62713ec82fae" begin
    import .CuArrays
    using .CuArrays: CuArray
    using .CuArrays.CUDAnative

    CuArrays.@cufunc function gelu(x)
        λ = oftype(x/1, √(2/π))
        α = oftype(x/1, 0.044715)
        h = oftype(x/1, 0.5)
        h * x * (one(x) + tanh(λ * (x + α * x^3)))
    end
end


text:
module FluxTraining

using LightGraphs
using BSON: @load, @save
using Flux
using Flux: Params, onecold
using Flux.Optimise: update!
using Glob
import OnlineStats
using OnlineStats: EqualWeight, Mean, OnlineStat
using Parameters
using ProgressMeter: Progress, next!
using Statistics: mean
using UUIDs
using Zygote
using Animations
using TensorBoardLogger: TBLogger, log_value, log_image, log_text, log_histogram, tb_overwrite
using Zygote: Grads, gradient
using ValueHistories
using DataStructures: DefaultDict
using PrettyTables

# functional
include("./functional/metrics.jl")

# callback system
include("./callbacks/protect.jl")
include("./callbacks/phases.jl")
include("./callbacks/events.jl")
include("./callbacks/callback.jl")
include("./callbacks/graph.jl")
include("./callbacks/execution.jl")

# logging
include("./callbacks/logging/Loggables.jl")
include("./callbacks/logging/logger.jl")
include("./callbacks/logging/tensorboard.jl")
include("./callbacks/logging/checkpointer.jl")

# callback implementations
include("./callbacks/conditional.jl")
include("./callbacks/callbacks.jl")
include("./callbacks/custom.jl")
include("./callbacks/metrics.jl")
include("./callbacks/recorder.jl")

# hyperparameter scheduling
include("./callbacks/hyperparameters.jl")
include("./callbacks/scheduler.jl")

# learner
include("./learner.jl")

include("./train.jl")

export AbstractCallback,
    Loss,
    ConditionalCallback,
    CancelBatchException,
    CancelEpochException,
    CancelFittingException,
    Checkpointer,
    CustomCallback,
    EarlyStopping,
    ToGPU,
    GarbageCollect,
    Learner,
    Metric,
    Recorder,
    ProgressPrinter,
    Metrics,
    MetricsPrinter,
    TrainingPhase,
    ValidationPhase,
    Schedule,
    Scheduler,
    LogMetrics,
    LogHistograms,
    LogHyperParams,
    LogVisualization,
    TensorBoardBackend,
    StopOnNaNLoss,
    LearningRate,
    throttle,
    accuracy,
    fit!,
    onecycle,
    loadmodel,
    savemodel
end # module


text:
#EasyPlot: A quick/easy way to generate, save, & display plots.
#-------------------------------------------------------------------------------
module EasyPlot

import ..CMDimData: rootpath

demofilelist() =
	[joinpath(rootpath, "sample", "plots", "demo$i.jl") for i in 1:5]

using Colors
using MDDatasets
import NumericIO
using Graphics: BoundingBox

import ..NullOr

#==Implementation
===============================================================================#
include("core.jl")
include("attributes.jl")
include("colors.jl")
include("text.jl")
include("cartesian.jl")
include("grids.jl")
include("annotation.jl")
include("base.jl")
include("multistrip.jl")
include("eyediag.jl")
include("datamd.jl")
include("themes.jl")
include("show.jl")
include("builder.jl")
include("defaults.jl")
include("io.jl")
include("display.jl")
include("doc.jl")

#==`cons`truct interface (minimize namespace pollution)
===============================================================================#
function cons(::DS{T}, args...; kwargs...) where T
	mlist = string(methods(cons, (DS,)))
msg = """
Cannot construct object of type :$T.
Supported methods are:\n
"""
	throw(ArgumentError(msg * mlist))
end
cons(s::Symbol, args...; kwargs...) = cons(DS(s), args...; kwargs...)

#==Exported interface
===============================================================================#
export cons #Construct objects (minimize namespace pollution)
export set #Set PlotCollection/Plot/Waveform/... attributes

#==Unexported tools available to backend-interfacing modules:
================================================================================
	addwfrm(ax::AbstractWfrmBuilder, ...) #
	buildeye() #Builds multi-dimensional DataEye objects from DataF1 Leaf elements.
	getcolor()

	#Constants:
	COLORSCHEME[{:default/}] #List of color schemes
==#

#==Initialization
===============================================================================#
function __init__()
	global defaults
	_initialize(defaults)
	return
end

end #EasyPlot

#Last line


text:
using ZMQ
const ctx = Context()
const sock = Socket(ctx, REP)
ZMQ.bind(sock, "tcp://127.0.0.1:*")
endpt = split(ZMQ.get_last_endpoint(sock), "\0")[1]
println("ZMQ bound to endpoint: $(endpt)")

while true
   print_with_color(:cyan,"Waiting for changes...\n")
   msg = ZMQ.recv(sock)
   reload(bytestring(msg))
   ZMQ.send(sock, "SUCCESS")
end

max_stars_repo_name: jbnicolai/jlbox

max_issues_repo_name: jbnicolai/jlbox

max_forks_repo_name: jbnicolai/jlbox


text:
function force(x, y, z)
    r2 = x^2 + y^2 + z^2
    n = sqrt(r2)
    return [x/(n*r2), y/(n*r2), z/(n*r2)]
end

function force_sum(pts)
    s = zeros(3)
    @time for i in 1:size(pts, 1) # ~0.3 seconds
        s += force(pts[i,:]...)
    end
    return s
end


text:
using Base.Test
using PValueAdjust

# is there a mocking package?
pValues = [1, 0.5, 0.1, 0.001, 0.0001, 0.002]
@test_approx_eq padjust(pValues) pValues
@test_throws DomainError padjust([10])


text:
using ForwardDiff
using AMDGPU

function f(x)
    x = x.^2
    nrm2 = eltype(x)(0.0)
    for v in x
        nrm2 += v
    end
    return nrm2
end

n = 10

# Works
t1s{N} = ForwardDiff.Dual{Nothing,Float64, N} where N
x = rand(n)
adx = Array{t1s{n}}(undef, n)
ForwardDiff.seed!(adx, x, ForwardDiff.construct_seeds(ForwardDiff.Partials{n,Float64}))
rocadx = ROCVector{t1s{n}}(adx)
rocadg = f(rocadx)

# Doesn't works
x = ROCVector{Float64}(undef, n)
adx = ROCVector{t1s{n}}(undef, n)
ForwardDiff.seed!(adx, x, ForwardDiff.construct_seeds(ForwardDiff.Partials{n,Float64}))
adg = adf(x)

# Doesn't work
df = x -> ForwardDiff.jacobian(f,x)
g = df(x) # gradient of F


text:
#=
Basic example of 1D localization against a scalar field
=#

using Images, FileIO           # load elevation from PNG
using DataInterpolations       # handles the terrain
using Distributions
using ImageFiltering
using Plots
using Caesar, RoME
gr()

## Terrain
# load dem (18x18km span)
img = load(joinpath(dirname(dirname(pathof(Caesar))), "examples/dev/scalar/dem.png")) .|> Gray
h = 1e3*Float64.( @view img[512,:])
# x = range(-9000,9000,length = length(h))
x = range(-100,100,length = length(h))
terrain = LinearInterpolation(h, x)
plot(terrain)

## Measurement likelihood
function measurementLikelihood( meas_model::SamplableBelief,
                                terrain::LinearInterpolation)
    # meas_model = Normal(y, sigma_y)
    l = pdf.(meas_model, terrain.u)
    return l./sum(l)
end

##
function propagate(prior, model)
    imfilter(prior, reflect(centered(model)), Fill(0) )
end

function neval(dist, x)
    d = pdf.(dist, x)
    d./sum(d)
end

xt = -70    # true position
delta = 10; # motion increment
sigma_y = 5 # altitude measurement uncertainty
motion_kernel = neval(Normal(delta, 1), x)
p = neval(Normal(xt, 1), x )  # prior

##

# struct ScalarField1D
#     # scalar field over one variable f:R->R (i.e. h= f(x))
#     # [a special case of ScalarField, h = f(·)]
#     # ...
# end

## implements a scalar measurement factor (e.g. elevation against known scalar map)
# addFactor!(fg, [:x1, :topography], ScalarFieldFactor(y))
# FIXME: objective is to <: IIF.AbstractRelativeConstant (final name TBD)
struct ScalarFieldFactor{T<:SamplableBelief} <: IIF.AbstractPrior
    elevationMeasModel::T
end

function IIF.getSample(s::CalcFactor{<:ScalarFieldFactor}, N::Int=1)
    # 1. generate measurement model -
    #   1a. (high-level): we get a single value (e.g. 20m elevation) - Normal(y, sigma)
    #   1b. (low-level):  we get a vector-valued measurement (e.g. back-scatter strength) - Normal(argmax(y),sigma) [classic] or BSS(y) [full measurement]
    #   NOTE: 1b is not the n-dimensional vector field scenario; it is just a more complicated (un-processed) measurement
    # 2. generate likelihood (eval measurement model over field)
    # 3. draw samples from likelihood

    # assuming user addConstant!(fg, :terrain, ScalarField1D(dem))
    global terrain # should come from s.constants[:terrain] ?????

    meas_model = s.factor.elevationMeasModel
    # `l` is likelihood of location over grid terrain.u, given elevMeasModel
    l = measurementLikelihood(meas_model, terrain) # gridded data

    # AliasingSS only works for 1D at this time
    # NOTE, if you stray out of region of support, things blow up!
    bss = AliasingScalarSampler([terrain.t;], l)

    return (reshape(rand(bss,N),1,N),)
end

# # keep track of IIF #1051
# # parameters are both in map/trajectory space X
# function (s::CalcFactor{<:ScalarFieldFactor})(likelihoodSample, pose)
#     return likelihoodSample - pose
# end

##

# create empty factor graph
fg = initfg()
# getSolverParams(fg).inflation = 3.0

addVariable!(fg, :x0, ContinuousEuclid{1})

# add terrain to factor graph
# addConstant!(fg, :em, ScalarField1D(dem))

# simulate for a few poses
xt = -70
delta = 10
for i=1:10
    xt = xt + delta   # true position
    yt = terrain(xt)  # true elevation

    # p = propagate(p, motion_kernel) # prior
    # plot!(x,p)

    addVariable!(fg, Symbol("x$(i)"), ContinuousEuclid{1}) # BoundedEuclid{1}

    # add propagation
    addFactor!(fg, [Symbol("x$j") for j in (i-1):i], LinearRelative(Normal(delta, 5)),graphinit=false)

    em = ScalarFieldFactor( Normal(yt, 20) ) # sigma_y
    # addFactor!(fg, [:terrain, Symbol("y$(i-1)")], em)
    # the hacky version
    addFactor!(fg, [Symbol("x$(i)");], em, graphinit=false)
end

## talk some solving

# still experimental, treeinit
getSolverParams(fg).graphinit = false
getSolverParams(fg).treeinit = true

solveTree!(fg)

##

smplOnMap = approxConv(fg, :x10f1, :x10)

histogram(smplOnMap[:], bins=100)

plot(terrain)
plot!([-100; 100],[436; 436])

##

using RoMEPlotting
Gadfly.set_default_plot_size(25cm,20cm)

##

pl1 = plotKDE(fg, sortDFG(ls(fg)))

pl1.coord = Coord.cartesian(xmin=-100,xmax=85)

pl2 = Gadfly.plot(x=terrain.t,y=terrain.u, Geom.line)

vstack(pl1, pl2)

##