text:
export
    MultiagentNGSIMEnvVideoMaker,
    reset,
    step,
    observation_space_spec,
    action_space_spec,
    obs_names,
    render

#=
Description:
    Multiagent NGSIM env that plays NGSIM trajectories, allowing a variable
    number of agents to simultaneously control vehicles in the scene

    Raunak: This is basically a copy of multiagent_ngsim_env.jl with just
    a few additions to enable color coded video making
=#
mutable struct MultiagentNGSIMEnvVideoMaker <: Env
    trajdatas::Vector{ListRecord}
    trajinfos::Vector{Dict}
    roadways::Vector{Roadway}
    roadway::Union{Nothing, Roadway} # current roadway
    scene::Scene
    rec::SceneRecord
    ext::MultiFeatureExtractor
    egoids::Vector{Int} # current ids of relevant ego vehicles
    ego_vehs::Vector{Union{Nothing, Vehicle}} # the ego vehicles
    traj_idx::Int # current index into trajdatas
    t::Int # current timestep in the trajdata
    h::Int # current maximum horizon for egoid
    H::Int # maximum horizon
    primesteps::Int # timesteps to prime the scene
    Δt::Float64

    # multiagent type members
    n_veh::Int # number of simultaneous agents
    remove_ngsim_veh::Bool # whether to remove ngsim veh from all scenes
    features::Array{Float64}

    # metadata
    epid::Int # episode id
    render_params::Dict # rendering options
    infos_cache::Dict # cache for infos intermediate results
    function MultiagentNGSIMEnvVideoMaker(
            params::Dict;
            trajdatas::Union{Nothing, Vector{ListRecord}} = nothing,
            trajinfos::Union{Nothing, Vector{Dict}} = nothing,
            roadways::Union{Nothing, Vector{Roadway}} = nothing,
            reclength::Int = 5,
            Δt::Float64 = .1,
            primesteps::Int = 50,
            H::Int = 50,
            n_veh::Int = 20,
            remove_ngsim_veh::Bool = false,
            render_params::Dict = Dict("zoom"=>5., "viz_dir"=>"/tmp"))
        param_keys = keys(params)
        @assert in("trajectory_filepaths", param_keys)

        # optionally overwrite defaults
        reclength = get(params, "reclength", reclength)
        primesteps = get(params, "primesteps", primesteps)
        H = get(params, "H", H)
        n_veh = get(params, "n_veh", n_veh)
        remove_ngsim_veh = get(params, "remove_ngsim_veh", remove_ngsim_veh)
        for (k,v) in get(params, "render_params", render_params)
            render_params[k] = v
        end

        # load trajdatas if not provided
        if trajdatas == nothing || trajinfos == nothing || roadways == nothing
            trajdatas, trajinfos, roadways = load_ngsim_trajdatas(
                params["trajectory_filepaths"],
                minlength=primesteps + H
            )
        end

        # build components
        scene_length = max_n_objects(trajdatas)
        scene = Scene(scene_length)
        rec = SceneRecord(reclength, Δt, scene_length)
        ext = build_feature_extractor(params)
        infos_cache = fill_infos_cache(ext)
        # features are stored in row-major order because they will be transferred
        # to python; this is inefficient in julia, but don't change or it will
        # break the python side of the interaction
        features = zeros(n_veh, length(ext))
        egoids = zeros(n_veh)
        ego_vehs = [nothing for _ in 1:n_veh]
        return new(
            trajdatas,
            trajinfos,
            roadways,
            nothing,
            scene,
            rec,
            ext,
            egoids, ego_vehs, 0, 0, 0, H, primesteps, Δt,
            n_veh, remove_ngsim_veh, features,
            0, render_params, infos_cache
        )
    end
end

#=
Description:
    Reset the environment. Note that this environment maintains the following
    invariant attribute: at any given point, all vehicles currently being controlled
    will end their episode at the same time. This simplifies the rest of the code
    by enforcing synchronized restarts, but it does somewhat limit the sets of
    possible vehicles that can simultaneously interact. With a small enough minimum
    horizon (H <= 250 = 25 seconds) and number of vehicle (n_veh <= 100)
    this should not be a problem. If you need to run with larger numbers then those
    implement an environment with asynchronous resets.

Args:
    - env: env to reset
    - dones: bool vector indicating which indices have reached a terminal state
        these must be either all true or all false
=#
function reset(
        env::MultiagentNGSIMEnvVideoMaker,
        dones::Vector{Bool} = fill!(Vector{Bool}(undef, env.n_veh), true);
        offset::Int=env.H + env.primesteps,
        random_seed::Union{Nothing, Int} = nothing)
    # enforce environment invariant reset property
    # (i.e., always either all true or all false)
    @assert (all(dones) || all(.!dones))
    # first == all at this point, so if first is false, skip the reset
    if !dones[1]
        return
    end

    # sample multiple ego vehicles
    # as stated above, these will all end at the same timestep
    env.traj_idx, env.egoids, env.t, env.h = sample_multiple_trajdata_vehicle(
        env.n_veh,
        env.trajinfos,
        offset,
        rseed=random_seed
    )

    # update / reset containers
    env.epid += 1
    empty!(env.rec)
    empty!(env.scene)

    # prime
    for t in env.t:(env.t + env.primesteps)
        get!(env.scene, env.trajdatas[env.traj_idx], t)
        if env.remove_ngsim_veh
            keep_vehicle_subset!(env.scene, env.egoids)
        end
        update!(env.rec, env.scene)
    end

    # set the ego vehicle
    for (i, egoid) in enumerate(env.egoids)
        vehidx = findfirst(egoid,env.scene)
        #TODO: Check for nothing

        env.ego_vehs[i] = env.scene[vehidx]
    end
    # set the roadway
    env.roadway = env.roadways[env.traj_idx]
    # env.t is the next timestep to load
    env.t += env.primesteps + 1
    # enforce a maximum horizon
    env.h = min(env.h, env.t + env.H)
    return get_features(env)
end

#=
Description:
    Propagate a single vehicle through an otherwise predeterined trajdata

Args:
    - env: environment to be stepped forward
    - action: array of floats that can be converted into an AccelTurnrate
=#
function _step!(env::MultiagentNGSIMEnvVideoMaker, action::Array{Float64})
    # make sure number of actions passed in equals number of vehicles
    @assert size(action, 1) == env.n_veh
    ego_states = Vector{VehicleState}(undef, env.n_veh)
    # propagate all the vehicles and get their new states
    for (i, ego_veh) in enumerate(env.ego_vehs)
        # convert action into form
	# ego_action = AccelTurnrate(action[i,:]...)
	ego_action = LatLonAccel(action[i,:]...) # RpB: To work with IDM+MOBIL
        # propagate the ego vehicle
        ego_states[i] = propagate(
            ego_veh,
            ego_action,
            env.roadway,
            env.Δt
        )
        # update the ego_veh
        env.ego_vehs[i] = Entity(ego_veh, ego_states[i])
    end

    # load the actual scene, and insert the vehicles into it
    get!(env.scene, env.trajdatas[env.traj_idx], env.t)
    if env.remove_ngsim_veh
        keep_vehicle_subset!(env.scene, env.egoids)
    end
    orig_vehs = Vector{Vehicle}(undef, env.n_veh)

    for (i, egoid) in enumerate(env.egoids)

    vehidx = findfirst(egoid, env.scene)
    # TODO: Check for nothing

    # track the original vehicle for validation / infos purposes

    orig_vehs[i] = env.scene[vehidx]

	# replace the original with the controlled vehicle

        env.scene[vehidx] = env.ego_vehs[i]

	# Raunak testing how to access lane id
	# Commented out for now as not relevant to reward augmentation
	# lane4print = env.ego_vehs[i].state.posF.roadind.tag.lane
	# println("Lane number = $lane4print")

	# Raunak's failed attempt at Ghost vehicle. The below treats both original
	# and policy driven as real cars and ends up showing the both to be colliding
	#push!(env.scene,env.ego_vehs[i])
    end

    # update rec with current scene
    update!(env.rec, env.scene)

    # compute info about the step
#    step_infos = Dict{String, Vector{Float64}}(
#        "rmse_pos"=>Float64[],
#        "rmse_vel"=>Float64[],
#        "rmse_t"=>Float64[],
#        "x"=>Float64[],
#        "y"=>Float64[],
#        "s"=>Float64[],
#        "phi"=>Float64[],
#    )
#    for i in 1:env.n_veh
#        push!(step_infos["rmse_pos"], sqrt(abs2((orig_vehs[i].state.posG - env.ego_vehs[i].state.posG))))
#        push!(step_infos["rmse_vel"], sqrt(abs2((orig_vehs[i].state.v - env.ego_vehs[i].state.v))))
#        push!(step_infos["rmse_t"], sqrt(abs2((orig_vehs[i].state.posF.t - env.ego_vehs[i].state.posF.t))))
#        push!(step_infos["x"], env.ego_vehs[i].state.posG.x)
#        push!(step_infos["y"], env.ego_vehs[i].state.posG.y)
#        push!(step_infos["s"], env.ego_vehs[i].state.posF.s)
#        push!(step_infos["phi"], env.ego_vehs[i].state.posF.ϕ)
#    end

    # Raunak adds in original vehicle properties
    step_infos = Dict{String, Vector{Float64}}(
        #"rmse_pos"=>Float64[],
        #"rmse_vel"=>Float64[],
        #"rmse_t"=>Float64[],
        "x"=>Float64[],
        "y"=>Float64[],
        "s"=>Float64[],
        "phi"=>Float64[],
	"orig_x"=>Float64[],
	"orig_y"=> Float64[],
	"orig_theta"=>Float64[],
	"orig_length"=>Float64[],
	"orig_width"=>Float64[]
    )
    for i in 1:env.n_veh
        #push!(step_infos["rmse_pos"], sqrt(abs2((orig_vehs[i].state.posG - env.ego_vehs[i].state.posG))))
        #push!(step_infos["rmse_vel"], sqrt(abs2((orig_vehs[i].state.v - env.ego_vehs[i].state.v))))
        #push!(step_infos["rmse_t"], sqrt(abs2((orig_vehs[i].state.posF.t - env.ego_vehs[i].state.posF.t))))
        push!(step_infos["x"], env.ego_vehs[i].state.posG.x)
        push!(step_infos["y"], env.ego_vehs[i].state.posG.y)
        push!(step_infos["s"], env.ego_vehs[i].state.posF.s)
        push!(step_infos["phi"], env.ego_vehs[i].state.posF.ϕ)
	push!(step_infos["orig_x"], orig_vehs[i].state.posG.x)
	push!(step_infos["orig_y"], orig_vehs[i].state.posG.y)
	push!(step_infos["orig_theta"], orig_vehs[i].state.posG.θ)
	push!(step_infos["orig_length"], orig_vehs[i].def.length)
	push!(step_infos["orig_width"], orig_vehs[i].def.width)
    end

    return step_infos
end

function _extract_rewards(env::MultiagentNGSIMEnvVideoMaker, infos::Dict{String, Array{Float64}})
    rewards = zeros(env.n_veh)
    R = 0

    for i in 1:env.n_veh
        if infos["is_colliding"][i] == 1
            rewards[i] -= R
        elseif infos["is_offroad"][i] == 1
            rewards[i] -= R
        elseif infos["hard_brake"][i] == 1
            rewards[i] -= (R*0.5) # braking hard is not as bad as a collision
        end
    end
    return rewards
end

function Base.step(env::MultiagentNGSIMEnvVideoMaker, action::Array{Float64})
    step_infos = _step!(env, action)

    # compute features and feature_infos
    features = get_features(env)
    feature_infos = _compute_feature_infos(env, features)

    # combine infos
    infos = merge(step_infos, feature_infos)

    # update env timestep to be the next scene to load
    env.t += 1

    # compute terminal
    terminal = env.t > env.h ? true : false
    terminal = [terminal for _ in 1:env.n_veh]
    # vectorized sampler does not call reset on the environment
    # but expects the environment to handle resetting, so do that here
    # note: this mutates env.features in order to return the correct obs when resetting
    reset(env, terminal)
	rewards = _extract_rewards(env, infos)
    return deepcopy(env.features), rewards, terminal, infos
end

function _compute_feature_infos(env::MultiagentNGSIMEnvVideoMaker, features::Array{Float64};
                                                         accel_thresh::Float64=-3.0)
    feature_infos = Dict{String, Array{Float64}}(
                "is_colliding"=>Float64[],
                "is_offroad"=>Float64[],
                "hard_brake"=>Float64[],
		"colliding_veh_ids"=>Float64[],
		"offroad_veh_ids"=>Float64[],
		"hardbrake_veh_ids"=>Float64[]
		)

    # Raunak explains: env.n_veh will be number of policy driven cars
    # Caution: This need not be the same as number of cars in the scene
    # Because the scene contains both policy driven cars and ngsim replay cars

    # Raunak debugging
    #sizeFeatures=size(features)
    #println("Compute repoting featurs are \n $sizeFeatures")
    #something = env.infos_cache
    #println("env.infos_cache =\n $something")
    #somethingelse = env.infos_cache["is_colliding_idx"]
    #println("env.infos_cache[is_colliding_idx]=$somethingelse")

    for i in 1:env.n_veh
	# Raunak debugging
	#println("i=$i\n")
	#egoVehid = env.ego_vehs[i].id
	#println("egoVeh[i] = $egoVehid")
	#infos_cache = env.infos_cache
	#println("env.infos_cache = $infos_cache")

	is_colliding = features[i, env.infos_cache["is_colliding_idx"]]
	#println("is_colliding=$is_colliding\n")
	is_offroad = features[i, env.infos_cache["out_of_lane_idx"]]
        accel = features[i, env.infos_cache["accel_idx"]]
        push!(feature_infos["hard_brake"], accel <= accel_thresh)
        push!(feature_infos["is_colliding"], is_colliding)
        push!(feature_infos["is_offroad"], is_offroad)

	# Raunak adding list of colliding ego ids into the feature list that gets passed to render
	if is_colliding==1
		push!(feature_infos["colliding_veh_ids"],env.ego_vehs[i].id)
		#println("Collision has happened see red")
	end
	if is_offroad==1
		push!(feature_infos["offroad_veh_ids"],env.ego_vehs[i].id)
		#println("Offroad has happened see yellow")
	end
	if accel <= accel_thresh
		push!(feature_infos["hardbrake_veh_ids"],env.ego_vehs[i].id)
		#println("Hard brake has happened see some color")
	end
    end
    return feature_infos
end

function AutoRisk.get_features(env::MultiagentNGSIMEnvVideoMaker)
    for (i, egoid) in enumerate(env.egoids)
	#println("i=$i\n")
	#println("egoid = $egoid\n")
    veh_idx = findfirst(egoid, env.scene)

    # TODO: Check for nothing

	#println("veh_idx=$veh_idx")
	pull_features!(env.ext, env.rec, env.roadway, veh_idx)
        env.features[i, :] = deepcopy(env.ext.features)

	# Raunak checking whether lane_id got added as a feature
	# Does changing code in the julia folder automatically reflect
	#feature4print = env.ext.features
	#println("features are $feature4print")
    end
    return deepcopy(env.features)
end

function observation_space_spec(env::MultiagentNGSIMEnvVideoMaker)
    low = zeros(length(env.ext))
    high = zeros(length(env.ext))
    feature_infos = feature_info(env.ext)
    for (i, fn) in enumerate(feature_names(env.ext))
        low[i] = feature_infos[fn]["low"]
        high[i] = feature_infos[fn]["high"]
    end
    infos = Dict("high"=>high, "low"=>low)
    return (length(env.ext),), "Box", infos
end
action_space_spec(env::MultiagentNGSIMEnvVideoMaker) = (2,), "Box", Dict("high"=>[4.,.15], "low"=>[-4.,-.15])
obs_names(env::MultiagentNGSIMEnvVideoMaker) = feature_names(env.ext)
vectorized(env::MultiagentNGSIMEnvVideoMaker) = true
num_envs(env::MultiagentNGSIMEnvVideoMaker) = env.n_veh

# Raunak defined this render function to enable color changing based on collisions and offroads
function render(
        env::MultiagentNGSIMEnvVideoMaker;
	infos=Dict(),
        egocolor::Vector{Float64}=[0.,0.,1.],
        camtype::String="follow",
        static_camera_pos::Vector{Float64}=[0.,0.],
        camera_rotation::Float64=0.,
        canvas_height::Int=800,
        canvas_width::Int=800)
    # define colors for all the vehicles
    carcolors = Dict{Int,Colorant}()
    egocolor = ColorTypes.RGB(egocolor...)

    # Loop over all the vehicles in the scene. Note these may be both policy driven and ngsim replay
    for veh in env.scene
	# If current vehicle is a policy driven vehicle then color it blue otherwise color it green
	carcolors[veh.id] = in(veh.id, env.egoids) ? egocolor : colorant"green"

	# If the current vehicle is in the list of colliding vehicles color it red
"""
	if in(veh.id,infos["colliding_veh_ids"])
		carcolors[veh.id] = ColorTypes.RGB([1.,0.,0.]...)
	end

	# If current vehicle is in the list of offroad vehicles color it yellow
	if in(veh.id,infos["offroad_veh_ids"])
		carcolors[veh.id]=ColorTypes.RGB([1.,1.,0.]...)
	end

	# If current vehicle is in the list of hard brakers then color it light blue
	if in(veh.id,infos["hardbrake_veh_ids"])
		carcolors[veh.id]=ColorTypes.RGB([0.,1.,1.]...)
	end
"""
    end

    # define a camera following the ego vehicle
    if camtype == "follow"
        # follow the first vehicle in the scene
        cam = AutoViz.CarFollowCamera{Int}(env.egoids[1], env.render_params["zoom"])
    elseif camtype == "static"
        cam = AutoViz.StaticCamera(VecE2(static_camera_pos...), env.render_params["zoom"])
    else
        error("invalid camera type $(camtype)")
    end

    # Raunak commented this out because it was creating rays that were being used for
    # some research that Tim had been doing
    overlays = [
        CarFollowingStatsOverlay(env.egoids[1], 2),
    #    NeighborsOverlay(env.egoids[1], textparams = TextParams(x = 600, y_start=300))
    ]

    # Raunak video plotting the ghost vehicle
    # See 'OrigVehicleOverlay' defined in AutoViz/src/2d/overlays.jl
    # to understand how the ghost vehicle is being plotted
#    overlays = [
#       CarFollowingStatsOverlay(env.egoids[1], 2)
#	,OrigVehicleOverlay(infos["orig_x"][1],infos["orig_y"][1],infos["orig_theta"][1],infos["orig_length"][1],infos["orig_width"][1])
#    ]

    # rendermodel for optional rotation
    # note that for this to work, you have to comment out a line in AutoViz
    # src/overlays.jl:27 `clear_setup!(rendermodel)` in render
    rendermodel = RenderModel()
    camera_rotate!(rendermodel, deg2rad(camera_rotation))

    # render the frame
    frame = render(
        env.scene,
        env.roadway,
        overlays,
        rendermodel = rendermodel,
        cam = cam,
        car_colors = carcolors,
        canvas_height=canvas_height,
        canvas_width=canvas_width
    )
"""
    # save the frame
    if !isdir(env.render_params["viz_dir"])
        mkdir(env.render_params["viz_dir"])
    end
    ep_dir = joinpath(env.render_params["viz_dir"], "episode_$(env.epid)")
    if !isdir(ep_dir)
        mkdir(ep_dir)
    end
    filepath = joinpath(ep_dir, "step_$(env.t).png")
    AutoViz.write_to_png(frame, filepath)

    # load and return the frame as an rgb array
    img = PyPlot.imread(filepath)
    return img
"""
    return frame
end

# Raunak trying to add the original vehicle as a ghost car in the video
# Erroring out when defined here. Worked fine when defined within AutoViz/src/2d/Overlays.jl
# which is the same location where CarFollowingStatsOverlay is defined
#mutable struct OrigVehicleOverlay <: SceneOverlay
	# structure attributes here
#	x::Float64
#	y::Float64
#	theta::Float64
#	length::Float64
#	width::Float64
#	color::Colorant

	# constructor function

#	function OrigVehicleOverlay(x::Float64,y::Float64,theta::Float64,length::Float64,
#				    width::Float64,color::Colorant=colorant"white")
#	new(x,y,theta,length,width,color)
#	end
#end

# Better have this function signature exactly same as CarFollowingStatsOverlay
#function render!(rendermodel::RenderModel,overlay::OrigVehicleOverlay,scene::Scene, roadway::Roadway)
#	verbose = true
#	if verbose
#		println("render! within OrigVehicleOverlay reporting here")
#		x=overlay.x
#		y=overlay.y
#		println("x=$x,y = $y")
#	end
#   add_instruction!(rendermodel, render_vehicle, (overlay.x, overlay.y, overlay.theta,
#						   overlay.length, overlay.width, overlay.color))
#    return rendermodel
#end


text:
# ---
# title: 1481. Least Number of Unique Integers after K Removals
# id: problem1481
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array, Sort
# link: <https://leetcode.com/problems/least-number-of-unique-integers-after-k-removals/description/>
# hidden: true
# ---
#
# Given an array of integers `arr` and an integer `k`. Find the _least number of
# unique integers_  after removing **exactly** `k` elements **.**
#
#
#
# **Example 1:**
#
#
#
#     Input: arr = [5,5,4], k = 1
#     Output: 1
#     **Explanation** : Remove the single 4, only 5 is left.
#
#
# **Example 2:**
#
#
#
#     Input: arr = [4,3,1,1,3,3,2], k = 3
#     Output: 2
#     **Explanation** : Remove 4, 2 and either one of the two 1s or three 3s. 1 and 3 will be left.
#
#
#
# **Constraints:**
#
#   * `1 <= arr.length <= 10^5`
#   * `1 <= arr[i] <= 10^9`
#   * `0 <= k <= arr.length`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end

ext: jl


text:
function extract_composition(elementstring::String)
    # (C,H,N,O,S,P)
    regexes = [r"[C]\d*", r"[H]\d*", r"[N]\d*", r"[O]\d*", r"[S]\d*", r"[P]\d*"]
    chemical_indices = Vector{Int64}(undef, 6)
    for i in eachindex(regexes)
        m = collect(m.match for m in eachmatch(regexes[i], elementstring))
        if size(m,1) == 0
            chemical_indices[i] = 0
        else
            try
                chemical_indices[i] = parse(Int64, split(m[1], r"[A-Z]")[2])
            catch e
                chemical_indices[i] = 1
            end
        end
    end
    return chemical_indices
end

function get_stoich_electron_donor(elementstring::String)
    chemical_indices = extract_composition(elementstring)
    a = chemical_indices[1] #C
    b = chemical_indices[2] #H
    c = chemical_indices[3] #N
    d = chemical_indices[4] #O
    e = chemical_indices[5] #S
    f = chemical_indices[6] #P
    z = 0
    return [-1,  -(3*a+4*e-d), a, c, e, f, 5*a+b-4*c-2*d+7*e-f, -z+4*a+b-3*c-2*d+5*e-2*f, 0, 0]
end

function get_stoich_electron_acceptor()
    stoich_electron_acceptor = zeros(10)
    stoich_electron_acceptor[9] = -1  # oxygen
    stoich_electron_acceptor[7] = -4  # h+
    stoich_electron_acceptor[8] = -4  # e-
    stoich_electron_acceptor[2] =  2  # h2o
    return stoich_electron_acceptor
end

function get_stoich_electron_acceptor_N()
    stoich_electron_acceptor = zeros(10)
    stoich_electron_acceptor[3] = 1  # nh4
    stoich_electron_acceptor[7] = +6  # h+
    stoich_electron_acceptor[8] = +5  # e-
    stoich_electron_acceptor[2] =  -3  # h2o
    return stoich_electron_acceptor
end

function get_stoich_catabolic_reaction(elementstring::String)
    stoich_electron_donor = get_stoich_electron_donor(elementstring)
    stoich_electron_acceptor = get_stoich_electron_acceptor()
    yEd = stoich_electron_donor[8]
    yEa = stoich_electron_acceptor[8]
    stoich_cat_rxns = stoich_electron_donor-(yEd/yEa)*stoich_electron_acceptor
    return stoich_cat_rxns
end

function get_stoich_catabolic_reaction_N(elementstring::String)
    stoich_electron_donor = get_stoich_electron_donor(elementstring)
    stoich_electron_acceptor = get_stoich_electron_acceptor_N()
    yEd = stoich_electron_donor[8]
    yEa = stoich_electron_acceptor[8]
    stoich_cat_rxns = stoich_electron_donor-(yEd/yEa)*stoich_electron_acceptor
    return stoich_cat_rxns
end

function get_stoich_anabolic_reaction(elementstring::String, chemFormBiom)
    chemical_indices = extract_composition(elementstring)
    stoich_electron_donor = get_stoich_electron_donor(elementstring)
    stoich_electron_acceptor = get_stoich_electron_acceptor()
    a = chemical_indices[1] #C

    #chemFormBiom = [1, 1.8, 0.2, 0.5, 0, 0, 0]  # C H_1.8 N_0.2 O_0.5 made variable
    aB = chemFormBiom[1]
    bB = chemFormBiom[2]
    cB = chemFormBiom[3]
    dB = chemFormBiom[4]
    eB = chemFormBiom[5]
    fB = chemFormBiom[6]
    zB = chemFormBiom[7]

    ySource = -1
    yH2o = -(3*aB+4*eB-dB)
    yHco3 = aB
    yNh4 = cB
    yHpo4 = eB
    yHs = fB
    yH = 5*aB+bB-4*cB-2*dB+7*eB-fB
    yE = -zB+4*aB+bB-3*cB-2*dB+5*eB-2*fB

    stoichAnStarB = [ySource,yH2o,yHco3,yNh4,yHpo4,yHs,yH,yE,0,0]
    stoichAnStarB = -stoichAnStarB
    stoichAnStarB[end] = stoichAnStarB[1]
    stoichAnStarB[1] = 0

    # O2
    stoichAnStar_O2 = stoichAnStarB+(1/a)*stoich_electron_donor
    yEana = stoichAnStar_O2[8]
    if yEana > 0.0
        yEa = stoich_electron_acceptor[8]
        stoichAn_O2 = stoichAnStar_O2-yEana/yEa*stoich_electron_acceptor
    elseif yEana < 0.0
        yEd = stoich_electron_donor[8]
        stoichAn_O2 = stoichAnStar_O2-yEana/yEd*stoich_electron_donor
    else
        stoichAn_O2 = stoichAnStar_O2
    end

    # HCO3-
    yEd = stoich_electron_donor[8]
    yEa = stoichAnStarB[8]
    stoichAn_HCO3 = stoich_electron_donor-(yEd/yEa)*stoichAnStarB
    stoichAn_HCO3 = stoichAn_HCO3/stoichAn_HCO3[10]

    return stoichAn_O2, stoichAn_HCO3
end

function get_lambda(elementstring::String, chemFormBiom)
    chemical_indices = extract_composition(elementstring)
    stoich_electron_donor = get_stoich_electron_donor(elementstring)
    stoich_cat_rxns = get_stoich_catabolic_reaction(elementstring)
    stoich_anabolic_O2, stoich_anabolic_HCO3 = get_stoich_anabolic_reaction(elementstring, chemFormBiom)

    a = chemical_indices[1] #C
    b = chemical_indices[2] #H
    c = chemical_indices[3] #N
    d = chemical_indices[4] #O
    e = chemical_indices[5] #S
    f = chemical_indices[6] #P
    z = 0

    ne = -z+4*a+b-3*c-2*d+5*e-2*f  # number of electrons transferred in D
    nosc = -(ne/a) + 4
    delGcox0PerC = 60.3-28.5*nosc  # kJ/C-mol
    delGcox0 = delGcox0PerC*a*abs(stoich_electron_donor[1])
    # - estimate delGf0 for electron donor
    delGf0_D_zero = 0
    delGf0_zero = [delGf0_D_zero, -237.2, -586.8, -79.3, -1096.1, 12.1, 0, 0, 16.4, -67]
    delGcox0_zero = dot(delGf0_zero, stoich_electron_donor)
    delGf0_D_est = (delGcox0-delGcox0_zero)/stoich_electron_donor[1]
    # - finally, delGf0
    delGf0 = delGf0_zero
    delGf0[1] = delGf0_D_est

    # - standard delG at pH=0
    delGcat0 = dot(delGf0, stoich_cat_rxns)
    delGan0_O2 = dot(delGf0, stoich_anabolic_O2)
    delGan0_HCO3 = dot(delGf0, stoich_anabolic_HCO3)

    # - stadard delG at pH=7
    R = 0.008314  # kJ/(K.mol)
    T = 298  # K
    iProton = 7  # [eD,h2o,hco3-,nh4+,hpo4**2-,hs-,h+,e-,eA,biom]
    delGcox = delGcox0+R*T*stoich_electron_donor[iProton]*log(1e-7)
    delGcat = delGcat0+R*T*stoich_cat_rxns[iProton]*log(1e-7)
    delGan_O2 = delGan0_O2+R*T*stoich_anabolic_O2[iProton]*log(1e-7)
    delGan_HCO3 = delGan0_HCO3+R*T*stoich_anabolic_HCO3[iProton]*log(1e-7)

    # The Thermodynamic Electron Equivalents Model (TEEM)
    # --------
    eta = 0.43
    delGsyn = 200 # kJ/(mol.X)

    if delGan_O2 < 0
        m_O2 = 1
    else
        m_O2 = -1
    end

    if delGan_HCO3 < 0
        m_HCO3 = 1
    else
        m_HCO3 = -1
    end

    lambda_O2 = (delGan_O2*eta^m_O2+delGsyn)/(-delGcat*eta)
    lambda_HCO3 = (delGan_HCO3*eta^m_HCO3+delGsyn)/(-delGcat*eta)

    if lambda_O2 > 0
        stoichMet_O2 = lambda_O2*stoich_cat_rxns+stoich_anabolic_O2
    else
        stoichMet_O2 = stoich_anabolic_O2
    end

    if lambda_HCO3 > 0
        stoichMet_HCO3 = lambda_HCO3*stoich_cat_rxns+stoich_anabolic_HCO3
    else
        stoichMet_HCO3 = stoich_anabolic_HCO3
    end

    delGdis_O2 = dot(delGf0,stoichMet_O2) + R*T*stoichMet_O2[iProton]*log(1e-7)
    delGdis_HCO3 = dot(delGf0,stoichMet_HCO3) + R*T*stoichMet_HCO3[iProton]*log(1e-7)

    return [lambda_O2, lambda_HCO3], [delGcox0PerC,delGcox0,delGcox,delGcat0,delGcat,delGan0_O2,delGan0_HCO3, delGan_O2,delGan_HCO3,delGdis_O2,delGdis_HCO3], [m_O2, m_HCO3], stoichMet_O2, stoichMet_HCO3
end

function get_lambda(elementstring::String, chemFormBiom, V_c, Min_gen_time, Gram_stain, ρ_p)
    chemical_indices = extract_composition(elementstring)
    stoich_electron_donor = get_stoich_electron_donor(elementstring)
    stoich_cat_rxns = get_stoich_catabolic_reaction(elementstring)
    stoich_anabolic_O2, stoich_anabolic_HCO3 = get_stoich_anabolic_reaction(elementstring, chemFormBiom)

    a = chemical_indices[1] #C
    b = chemical_indices[2] #H
    c = chemical_indices[3] #N
    d = chemical_indices[4] #O
    e = chemical_indices[5] #S
    f = chemical_indices[6] #P
    z = 0

    ne = -z+4*a+b-3*c-2*d+5*e-2*f  # number of electrons transferred in D
    nosc = -(ne/a) + 4
    delGcox0PerC = 60.3-28.5*nosc  # kJ/C-mol
    delGcox0 = delGcox0PerC*a*abs(stoich_electron_donor[1])
    # - estimate delGf0 for electron donor
    delGf0_D_zero = 0
    delGf0_zero = [delGf0_D_zero, -237.2, -586.8, -79.3, -1096.1, 12.1, 0, 0, 16.4, -67]
    delGcox0_zero = dot(delGf0_zero, stoich_electron_donor)
    delGf0_D_est = (delGcox0-delGcox0_zero)/stoich_electron_donor[1]
    # - finally, delGf0
    delGf0 = delGf0_zero
    delGf0[1] = delGf0_D_est

    # - standard delG at pH=0
    delGcat0 = dot(delGf0, stoich_cat_rxns)
    delGan0_O2 = dot(delGf0, stoich_anabolic_O2)
    delGan0_HCO3 = dot(delGf0, stoich_anabolic_HCO3)

    # - stadard delG at pH=7
    R = 0.008314  # kJ/(K.mol)
    T = 298  # K
    iProton = 7  # [eD,h2o,hco3-,nh4+,hpo4**2-,hs-,h+,e-,eA,biom]
    delGcox = delGcox0+R*T*stoich_electron_donor[iProton]*log(1e-7)
    delGcat = delGcat0+R*T*stoich_cat_rxns[iProton]*log(1e-7)
    delGan_O2 = delGan0_O2+R*T*stoich_anabolic_O2[iProton]*log(1e-7)
    delGan_HCO3 = delGan0_HCO3+R*T*stoich_anabolic_HCO3[iProton]*log(1e-7)

    # The Thermodynamic Electron Equivalents Model (TEEM)
    # --------
    eta = 0.43
    #
    r_p = 1e-9
    N_A = 6.022e23
    r_c = (3*V_c/(4*pi)).^(1/3)   # cell radius
    gmax = log(2)./Min_gen_time
    N_p = @. 4*ρ_p*r_c^2/r_p^2
    N_SB = transporter_density_to_monomer_uptake_sites(V_c, ρ_p*ones(1,1), Min_gen_time, Gram_stain)[1]
    delG_porter = @. 1.09e-19*60^2*gmax*N_SB*N_A/1e3*6 # kJ/mol.X

    delGsyn = 4.18 + delG_porter[1]  # kJ/(mol.X)

    if delGan_O2 < 0
        m_O2 = 1
    else
        m_O2 = -1
    end

    if delGan_HCO3 < 0
        m_HCO3 = 1
    else
        m_HCO3 = -1
    end

    lambda_O2 = (delGan_O2*eta^m_O2+delGsyn)/(-delGcat*eta)
    lambda_HCO3 = (delGan_HCO3*eta^m_HCO3+delGsyn)/(-delGcat*eta)

    if lambda_O2 > 0
        stoichMet_O2 = lambda_O2*stoich_cat_rxns+stoich_anabolic_O2
    else
        stoichMet_O2 = stoich_anabolic_O2
    end

    if lambda_HCO3 > 0
        stoichMet_HCO3 = lambda_HCO3*stoich_cat_rxns+stoich_anabolic_HCO3
    else
        stoichMet_HCO3 = stoich_anabolic_HCO3
    end

    delGdis_O2 = dot(delGf0,stoichMet_O2) + R*T*stoichMet_O2[iProton]*log(1e-7)
    delGdis_HCO3 = dot(delGf0,stoichMet_HCO3) + R*T*stoichMet_HCO3[iProton]*log(1e-7)

    return [lambda_O2, lambda_HCO3], [delGcox0PerC,delGcox0,delGcox,delGcat0,delGcat,delGan0_O2,delGan0_HCO3, delGan_O2,delGan_HCO3,delGdis_O2,delGdis_HCO3], [m_O2, m_HCO3], stoichMet_O2, stoichMet_HCO3
end

max_stars_repo_name: giannamars/DEBMicroTrait

max_issues_repo_name: giannamars/DEBMicroTrait

max_forks_repo_name: giannamars/DEBMicroTrait


text:
module ScrewDimensions
using DataFrames
using Unitful
using UnitfulRecipes
using RecipesBase

export diagonal_length

"""
```SCREWS```

A table of screw dimensions, from

http://www.emuge-franken-bg.com/attachments/article/97/15%20Gewindetabellen.pdf

and

https://en.wikipedia.org/wiki/Unified_Thread_Standard
"""
const SCREWS = DataFrame(;
                         name=Symbol[],
                         nominal_diameter=empty([1.0]u"mm"), # d = D
                         pitch=empty([1.0]u"mm"), # P
                         average_diameter=empty([1.0]u"mm"), # d_2 = D_2
                         outer_minor=empty([1.0]u"mm"), # d_3
                         inner_minor=empty([1.0]u"mm"), # D_1
                         outer_depth=empty([1.0]u"mm"), # h_3
                         inner_depth=empty([1.0]u"mm"), # H
                         rounding=empty([1.0]u"mm"), # R
                        )

push!(SCREWS, (:M1,            1u"mm",   0.25u"mm", 0.838u"mm",  0.693u"mm",  0.729u"mm",  0.153u"mm", 0.135u"mm", 0.036u"mm"))
push!(SCREWS, (Symbol("M1.1"), 1.1u"mm", 0.25u"mm", 0.938u"mm",  0.793u"mm",  0.829u"mm",  0.153u"mm", 0.135u"mm", 0.036u"mm"))
push!(SCREWS, (Symbol("M1.2"), 1.2u"mm", 0.25u"mm", 1.038u"mm",  0.893u"mm",  0.929u"mm",  0.153u"mm", 0.135u"mm", 0.036u"mm"))
push!(SCREWS, (Symbol("M1.4"), 1.4u"mm", 0.3u"mm",  1.205u"mm",  1.032u"mm",  1.075u"mm",  0.184u"mm", 0.162u"mm", 0.043u"mm"))
push!(SCREWS, (Symbol("M1.6"), 1.6u"mm", 0.35u"mm", 1.373u"mm",  1.171u"mm",  1.221u"mm",  0.215u"mm", 0.189u"mm", 0.051u"mm"))
push!(SCREWS, (Symbol("M1.8"), 1.8u"mm", 0.35u"mm", 1.573u"mm",  1.371u"mm",  1.421u"mm",  0.215u"mm", 0.189u"mm", 0.051u"mm"))
push!(SCREWS, (:M2,            2u"mm",   0.4u"mm",  1.740u"mm",  1.509u"mm",  1.567u"mm",  0.245u"mm", 0.217u"mm", 0.058u"mm"))
push!(SCREWS, (Symbol("M2.2"), 2.2u"mm", 0.45u"mm", 1.908u"mm",  1.648u"mm",  1.713u"mm",  0.276u"mm", 0.244u"mm", 0.065u"mm"))
push!(SCREWS, (Symbol("M2.5"), 2.5u"mm", 0.45u"mm", 2.208u"mm",  1.948u"mm",  2.013u"mm",  0.276u"mm", 0.244u"mm", 0.065u"mm"))
push!(SCREWS, (:M3,            3u"mm",   0.5u"mm",  2.675u"mm",  2.387u"mm",  2.459u"mm",  0.307u"mm", 0.271u"mm", 0.072u"mm"))
push!(SCREWS, (Symbol("M3.5"), 3.5u"mm", 0.6u"mm",  3.110u"mm",  2.764u"mm",  2.850u"mm",  0.368u"mm", 0.325u"mm", 0.087u"mm"))
push!(SCREWS, (:M4,            4u"mm",   0.7u"mm",  3.545u"mm",  3.141u"mm",  3.242u"mm",  0.429u"mm", 0.379u"mm", 0.101u"mm"))
push!(SCREWS, (Symbol("M4.5"), 4.5u"mm", 0.75u"mm", 4.013u"mm",  3.580u"mm",  3.688u"mm",  0.460u"mm", 0.406u"mm", 0.108u"mm"))
push!(SCREWS, (:M5,            5u"mm",   0.8u"mm",  4.480u"mm",  4.019u"mm",  4.134u"mm",  0.491u"mm", 0.433u"mm", 0.115u"mm"))
push!(SCREWS, (:M6,            6u"mm",   1u"mm",    5.350u"mm",  4.773u"mm",  4.917u"mm",  0.613u"mm", 0.541u"mm", 0.144u"mm"))
push!(SCREWS, (:M7,            7u"mm",   1u"mm",    6.350u"mm",  5.773u"mm",  5.917u"mm",  0.613u"mm", 0.541u"mm", 0.144u"mm"))
push!(SCREWS, (:M8,            8u"mm",   1.25u"mm", 7.188u"mm",  6.466u"mm",  6.647u"mm",  0.767u"mm", 0.677u"mm", 0.180u"mm"))
push!(SCREWS, (:M9,            9u"mm",   1.25u"mm", 8.188u"mm",  7.466u"mm",  7.647u"mm",  0.767u"mm", 0.677u"mm", 0.180u"mm"))
push!(SCREWS, (:M10,           10u"mm",  1.5u"mm",  9.026u"mm",  8.160u"mm",  8.376u"mm",  0.920u"mm", 0.812u"mm", 0.217u"mm"))
push!(SCREWS, (:M11,           11u"mm",  1.5u"mm",  10.026u"mm", 9.160u"mm",  9.376u"mm",  0.920u"mm", 0.812u"mm", 0.217u"mm"))
push!(SCREWS, (:M12,           12u"mm",  1.75u"mm", 10.863u"mm", 9.853u"mm",  10.106u"mm", 1.074u"mm", 0.947u"mm", 0.253u"mm"))
push!(SCREWS, (:M14,           14u"mm",  2u"mm",    12.701u"mm", 11.546u"mm", 11.835u"mm", 1.227u"mm", 1.083u"mm", 0.289u"mm"))
push!(SCREWS, (:M16,           16u"mm",  2u"mm",    14.701u"mm", 13.546u"mm", 13.835u"mm", 1.227u"mm", 1.083u"mm", 0.289u"mm"))
push!(SCREWS, (:M18,           18u"mm",  2.5u"mm",  16.376u"mm", 14.933u"mm", 15.294u"mm", 1.534u"mm", 1.353u"mm", 0.361u"mm"))
push!(SCREWS, (:M20,           20u"mm",  2.5u"mm",  18.376u"mm", 16.933u"mm", 17.294u"mm", 1.534u"mm", 1.353u"mm", 0.361u"mm"))
push!(SCREWS, (:M22,           22u"mm",  2.5u"mm",  20.376u"mm", 18.933u"mm", 19.294u"mm", 1.534u"mm", 1.353u"mm", 0.361u"mm"))
push!(SCREWS, (:M24,           24u"mm",  3u"mm",    22.051u"mm", 20.319u"mm", 20.752u"mm", 1.840u"mm", 1.624u"mm", 0.433u"mm"))
push!(SCREWS, (:M27,           27u"mm",  3u"mm",    25.051u"mm", 23.319u"mm", 23.752u"mm", 1.840u"mm", 1.624u"mm", 0.433u"mm"))
push!(SCREWS, (:M30,           30u"mm",  3.5u"mm",  27.727u"mm", 25.706u"mm", 26.211u"mm", 2.147u"mm", 1.894u"mm", 0.505u"mm"))
push!(SCREWS, (:M33,           33u"mm",  3.5u"mm",  30.727u"mm", 28.706u"mm", 29.211u"mm", 2.147u"mm", 1.894u"mm", 0.505u"mm"))
push!(SCREWS, (:M36,           36u"mm",  4u"mm",    33.402u"mm", 31.093u"mm", 31.670u"mm", 2.454u"mm", 2.165u"mm", 0.577u"mm"))
push!(SCREWS, (:M39,           39u"mm",  4u"mm",    36.402u"mm", 34.093u"mm", 34.670u"mm", 2.454u"mm", 2.165u"mm", 0.577u"mm"))
push!(SCREWS, (:M42,           42u"mm",  4.5u"mm",  39.077u"mm", 36.479u"mm", 37.129u"mm", 2.760u"mm", 2.436u"mm", 0.650u"mm"))
push!(SCREWS, (:M45,           45u"mm",  4.5u"mm",  42.077u"mm", 39.479u"mm", 40.129u"mm", 2.760u"mm", 2.436u"mm", 0.650u"mm"))
push!(SCREWS, (:M48,           48u"mm",  5u"mm",    44.752u"mm", 41.866u"mm", 42.587u"mm", 3.067u"mm", 2.706u"mm", 0.722u"mm"))
push!(SCREWS, (:M52,           52u"mm",  5u"mm",    48.752u"mm", 45.866u"mm", 46.587u"mm", 3.067u"mm", 2.706u"mm", 0.722u"mm"))
push!(SCREWS, (:M56,           56u"mm",  5.5u"mm",  52.428u"mm", 49.252u"mm", 50.046u"mm", 3.374u"mm", 2.977u"mm", 0.794u"mm"))
push!(SCREWS, (:M60,           60u"mm",  5.5u"mm",  56.428u"mm", 53.252u"mm", 54.046u"mm", 3.374u"mm", 2.977u"mm", 0.794u"mm"))
push!(SCREWS, (:M64,           64u"mm",  6u"mm",    60.103u"mm", 56.639u"mm", 57.505u"mm", 3.681u"mm", 3.248u"mm", 0.866u"mm"))
push!(SCREWS, (:M68,           68u"mm",  6u"mm",    64.103u"mm", 60.639u"mm", 61.505u"mm", 3.681u"mm", 3.248u"mm", 0.866u"mm"))

macro uncScrew(name, D, p)
    quote
        uncScrewHelper($name, $D, $p)
    end
end
function uncScrewHelper(name, D, p)
    H = sqrt(3)/2*p
    D_min = D - 3*sqrt(3)/8*p
    D_avg = (D_min + D)/2
    push!(SCREWS,
          (name, D*u"inch", p*u"inch", D_avg*u"inch", D_min*u"inch",
           D_min*u"inch", H*u"inch", H*u"inch", 0u"mm")
         )
end
@uncScrew :UNC1                 0.0730  0.015625
@uncScrew :UNC2                 0.0860  0.017857
@uncScrew :UNC3                 0.0990  0.020833
@uncScrew :UNC4                 0.1120  0.025000
@uncScrew :UNC5                 0.1250  0.025000
@uncScrew :UNC6                 0.1380  0.031250
@uncScrew :UNC8                 0.1640  0.031250
@uncScrew :UNC10                0.1900  0.041667
@uncScrew :UNC12                0.2160  0.041667
@uncScrew Symbol("UNC1/4")      0.2500  0.050000
@uncScrew Symbol("UNC5/16")     0.3125  0.055556
@uncScrew Symbol("UNC3/8")      0.3750  0.062500
@uncScrew Symbol("UNC7/16")     0.4375  0.071428
@uncScrew Symbol("UNC1/2")      0.5000  0.076923
@uncScrew Symbol("UNC9/16")     0.5625  0.083333
@uncScrew Symbol("UNC5/8")      0.6250  0.090909
@uncScrew Symbol("UNC3/4")      0.7500  0.100000
@uncScrew Symbol("UNC7/8")      0.8750  0.111111
@uncScrew Symbol("UNC1/1")      1.0000  0.125000

# Plot recipe
@userplot ScrewPlot

@recipe function f(h::ScrewPlot)
    measured_x = length(h.args) >= 1 ? h.args[1] : nothing
    measured_y = length(h.args) >= 2 ? h.args[2] : nothing
    xguide := "D"
    yguide := "d"
    @series begin
        seriestype := :scatter
        primary := false
        seriescolor := startswith.(String.(SCREWS[!,:name]), "M")*1
        markersize := 3
        label := nothing
        SCREWS[!, :nominal_diameter], diagonal_length.(SCREWS[!, :name])
    end

    if !isnothing(measured_x)
        rows = sort(SCREWS, order(:nominal_diameter, by=t->abs2(t-measured_x)))
        xs = ustrip.(get(plotattributes, :xunit, u"mm"), rows[1:4, :nominal_diameter])
        push!(xs, ustrip(get(plotattributes, :xunit, u"mm"), measured_x))
        ys = ustrip.(get(plotattributes, :yunit, u"mm"), rows[1:4, :outer_minor])
        xlims := extrema(xs)
        xwiden := true
        if !isnothing(measured_y)
            push!(ys, ustrip(get(plotattributes, :yunit, u"mm"), measured_y))
            ylims := extrema(ys)
            ywiden := true
            @series begin
                seriestype := :scatter
                seriescolor := :red
                markersize := 5
                label := "Measured"
                [measured_x], [measured_y]
            end
        else
            ylims := extrema(ys)
            ywiden := true
            @series begin
                seriestype := :vline
                seriescolor := :red
                label := "Measured"
                [measured_x]
            end
        end
        @series begin
            seriestype := :scatter
            primary := true
            series_annotations := [(s, -45.0, "Computer Modern") for s in SCREWS[!, :name]]
            seriescolor := startswith.(String.(SCREWS[!,:name]), "M")*1
            label := nothing
            SCREWS[!, :nominal_diameter], SCREWS[!, :outer_minor]
        end
    end
end

"""
```julia
screwplot()
```

Plot all the screws in the register, minor diameter vs major. Above each point will also be
the diagonal diameter, i.e. the diameter as measured by a caliper between a point along the
thread and another half a turn later, with a smaller marker.

The major diameter of a measured screw can be supplied to zoom the ``x`` axis to the nearest
couple of suggestions, as well as mark this value in the plot.

The minor diameter can be supplied as a second argument and will be used to further pinpoint
your screw in the plot.

## Example

A screw is somewhat haphazardly measured to be ``9.8\\;\\mathrm{mm}`` in diameter. The user
runs

```julia
screwplot(9.8u"mm")
```

and sees that the closest alternatives are M10 and UNC⅜. To differentiate them, the user
measures the diagonal minor diameter to be ``8.0\\;\\mathrm{mm}``. Now

```julia
using ScrewDimensions, Plots, Unitful
screwplot(9.8u"mm", 8.0u"mm")
```

shows that the closest screw is indeed M10.
"""
function screwplot end

"""
```julia
diagonal_length(name)
```
Find the diameter of a screw, as measured between two consecutive threads, as one would with
calipers.
"""
function diagonal_length(name)
    screw = SCREWS[SCREWS.name .=== name,:][1,:]
    return sqrt((screw.pitch/2)^2 + screw.outer_minor^2)
end

end


text:
const LEGEND_KINDS = (:line, :line3d, :bar, :errorbar)

const LEGEND_LOCATIONS = Dict(
    :left => [2, 3, 6],
    :center_h => [8, 9, 10],
    :right_out => [11, 12, 13],
    :bottom => [3, 4, 8],
    :center_v => [5, 6, 7, 10, 12]
)

"""
    Legend(size::Tuple{Float64, Float64}, cursors::Vector{Tuple{Float64, Float64}})

Return a `Legend` object.

This type defines the frame where a legend is plotted.
The fields contained in a `Legend` object are a 2-tuple with the size
of the legend box in NDC (width and height, respectively), and a vector of
2-tuples with the positions of the legend items.
"""
struct Legend
    size::Tuple{Float64, Float64}
    cursors::Vector{Tuple{Float64, Float64}}
end

"""
    Legend(geoms, frame [, maxrows])

Return a `Legend` object defined by the collection of geometries that
are meant to be referred to in the legend (`geoms`), and the dimensions
(width, height) of the `frame` in which the legend should be drawn.

The geometries are used to set the number of items to be drawn in the
legend, and their labels. The frame is used to estimate the font size
of the labels.

Optionally, this constructor can take the maximum number of items that
are represented in each column of the legend. Only the items in that collection
of geometries where `label` is not empty will be included.
"""
function Legend(geoms::Array{<:Geometry}, frame, maxrows=length(geoms))
    cursors = Tuple{Float64, Float64}[]
    charsize = charheight(frame)
    row = 0
    x = 0.08         # width reserved for the guide
    y = -0.015       # vertical top margin
    labelwidth = 0.0
    w = h = 0.0      # width and height of the full box
    scale = Int(GR.inqscale())
    # GR.selntran(0)
    GR.setscale(0)
    for g in geoms
        if !isempty(g.label) && g.kind ∈ LEGEND_KINDS
            row += 1
            # New column if the limit is exceeded
            if row > maxrows
                (-y > h) && (h = -y) # increase height if needed
                y = -0.015
                x += labelwidth + 0.08
                labelwidth = 0.0
                row = 1
            end
            sz = stringsize(g.label, charsize)
            (sz[1] > labelwidth) && (labelwidth = sz[1]) # increase label width
            dy = max(sz[2] - 0.03, 0.0) # height of the item
            push!(cursors, (x, y - dy))
            y -= dy + 0.03 # add item height and vertical margin
        end
    end
    GR.setscale(scale)
    # GR.selntran(1)
    if !isempty(cursors)
        # Define width and height
        (-y > h) && (h = -y)
        w = x + labelwidth + 0.015
        return Legend((w, h), cursors)
    else
        return Legend()
    end
end

const EMPTYLEGEND = Legend(NULLPAIR, Tuple{Float64, Float64}[])
Legend() = EMPTYLEGEND

"""
    legend_box(frame, (w, h), location)

Define the rectangle for a legend box, given the main `frame` with respect to
which it is located, the width and height of the box — in the tuple `(w, h)`,
and the location of the legend — an integer that defines the location code of
[Matplotlib legends](https://matplotlib.org/3.1.1/api/_as_gen/matplotlib.pyplot.legend.html).

All the values are given in NDC.
"""
function legend_box(frame, (w, h), location)
    if location ∈ LEGEND_LOCATIONS[:right_out]
        px = frame[2] + 0.01
    elseif location ∈ LEGEND_LOCATIONS[:center_h]
        px = 0.5 * (frame[1] + frame[2] - w)
    elseif location ∈ LEGEND_LOCATIONS[:left]
        px = frame[1] + 0.03
    else
        px = frame[2] - 0.03 - w
    end
    if location ∈ LEGEND_LOCATIONS[:center_v]
        py = 0.5 * (frame[3] + frame[4] + h)
    elseif location == 13
        py = frame[3] + h
    elseif location ∈ LEGEND_LOCATIONS[:bottom]
        py = frame[3] + h + 0.03
    elseif location == 11
        py = frame[4]
    else
        py = frame[4] - 0.03
    end
    (px, px + w, py - h, py)
end

####################
## `draw` methods ##
####################

# `location` is an integer code that defines the location of the legend with
# respect to the main plot area — as defined in
# [Matplotlib legends](https://matplotlib.org/3.1.1/api/_as_gen/matplotlib.pyplot.legend.html).

function draw(lg::Legend, geoms, location=1)
    # Do not draw if the legend is empty or is not meant to be plotted
    (lg == EMPTYLEGEND || location == 0) && return nothing
    # First draw the frame
    GR.savestate()
    # GR.selntran(0)
    viewport = legend_box(GR.inqviewport(), lg.size, location)
    GR.setviewport(viewport...)
    w, h = lg.size
    window = (0, w, -h, 0)
    GR.setwindow(window...)
    GR.setscale(0)
    # Fill rectangle with background color
    GR.setfillintstyle(GR.INTSTYLE_SOLID)
    GR.setfillcolorind(0)
    GR.fillrect(window...)
    # Draw border
    GR.setlinetype(GR.LINETYPE_SOLID)
    GR.setlinecolorind(1)
    GR.setlinewidth(1)
    GR.drawrect(window...)
    # Draw the geometries
    GR.uselinespec(" ")
    GR.settextalign(GR.TEXT_HALIGN_LEFT, GR.TEXT_VALIGN_HALF)
    c = 1
    for g in geoms
        (c > length(lg.cursors)) && break
        cursor = lg.cursors[c]
        if !isempty(g.label) && g.kind ∈ LEGEND_KINDS
            guide(g, cursor[1] - 0.04, cursor[2])
            text(cursor[1], cursor[2], g.label, true)
            c += 1
        else
            GR.uselinespec("")
        end
    end
    resetcolors()
    # GR.selntran(1)
    GR.restorestate()
    return nothing
end

"""
    guide(g, x, y)

Draw a guide to the geometry `g` in the coordinates `x`, `y` of the legend box.

The kind-specific method `guide(::Val{kind}, x, y)` defines the
low-level plotting instructions for the given kind of `g`.
"""
function guide(g::Geometry, x, y)
    GR.savestate()
    GR.settransparency(get(g.attributes, :alpha, 1.0))
    guide(Val(g.kind), g, x, y)
    GR.restorestate()
end

function guide(::Val{:line}, g, x, y)
    mask = _uselinespec(g.spec, g.attributes)
    if hasline(mask)
        GR.setlinewidth(float(get(g.attributes, :linewidth, 1.0)))
        GR.polyline(x .+ [-0.03, 0.03], y .+ [0.0, 0.0])
    end
    if hasmarker(mask)
        GR.setmarkersize(2float(get(g.attributes, :markersize, 1.0)))
        GR.polymarker([x], [y])
    end
    return nothing
end
guide(::Val{:line3d}, args...) = guide(Val(:line), g, x, y)

function guide(::Val{:errorbar}, g, x, y)
    mask = _uselinespec(g.spec, g.attributes)
    if hasline(mask)
        GR.setlinewidth(float(get(g.attributes, :linewidth, 1.0)))
        GR.polyline([x-0.03, x+0.03], [y, y]) # main bar
        GR.polyline([x-0.03, x-0.03], [y-0.01, y+0.01]) # left bar
        GR.polyline([x+0.03, x+0.03], [y-0.01, y+0.01]) # right bar
    end
    if hasmarker(mask)
        GR.setmarkersize(2float(get(g.attributes, :markersize, 1.0)))
        GR.polymarker([x], [y])
    end
    return nothing
end

function guide(::Val{:bar}, g, x, y)
    if haskey(g.attributes, :color)
        colorind = colorindex(Int(g.attributes[:color]))
    else
        ind = get(COLOR_INDICES, :barfill, 0)
        ind = COLOR_INDICES[:barfill] = ind + 1
        colorind = SERIES_COLORS[ind]
    end
    GR.setfillcolorind(colorind)
    GR.setfillintstyle(GR.INTSTYLE_SOLID)
    GR.fillrect(x - 0.03, x + 0.03, y - 0.012, y + 0.012)
    GR.setfillcolorind(1)
    GR.setfillintstyle(GR.INTSTYLE_HOLLOW)
    GR.fillrect(x - 0.03, x + 0.03, y - 0.012, y + 0.012)
end

guide(kind, g, x, y) = nothing


text:
"""
	version()

Returns the current Garamond version using the `Project.toml` and `git`.
If the `Project.toml`, `git` are not available, the version defaults to
an empty string.
"""
function version()
    commit = DEFAULT_VERSION_COMMIT
    date = DEFAULT_VERSION_DATE
    ver = DEFAULT_VERSION
    try
        root_path = "/" * joinpath(split(pathof(Garamond), "/")[1:end-2]...)

        # Check that the current git repository is the garamond one
        _check_commit = split(Garamond.DEFAULT_VERSION_COMMIT, "*")[1]
        cd(root_path)
        @assert occursin("master", read(pipeline(`git branch --contains $(_check_commit)`,
                                                 stderr=devnull), String))
        # Try reading the latest commit and date
        gitstr = read(pipeline(`git show --oneline -s --format="%h%ci"`, stderr=devnull), String)
        commit = gitstr[1:7]
        date = gitstr[8:17]
        # Try reading the version from the Project.toml file
        project_file = "Project.toml"
        open(project_file) do fid
            for line in eachline(fid)
                if occursin("=", line)
                    id, ver = strip.(split(line, "="))
                    id == "version" && break
                end
            end
        end
    catch
        # do nothing
    end
	return ver, commit, date
end

"""
    printable_version()

Returns a pretty version string that includes the git commit and date.
"""
function printable_version()
    ver, commit, date = version()
    vstr = ""
    if !isempty(ver)
        vstr = vstr * "v"*ver
    end
    if !isempty(commit)
        vstr = vstr * " commit: $commit"
    end
    if !isempty(date)
        vstr = vstr * " ($date)"
    end
    return vstr
end


text:
using LinearAlgebra
using BenchmarkTools
using InteractiveUtils

# test
T = 26
nx = 2
nu = 1
nw = 0
ft = (x, u) -> 1.0 * dot(x, x) + 0.1 * dot(u, u)
fT = (x, u) -> 1.0 * dot(x, x)
ct = Cost(ft, nx, nu, [t for t = 1:T-1])
cT = Cost(fT, nx, 0, [T])
obj = [ct, cT]

grad = zeros((T - 1) * (nx + nu) + nx)

function dynamics(z, u, w)
    mass = 1.0
    lc = 1.0
    gravity = 9.81
    damping = 0.1

    [z[2], (u[1] / ((mass * lc * lc)) - gravity * sin(z[1]) / lc - damping * z[2] / (mass * lc * lc))]
end

h = 0.1
function discrete_dynamics(y, x, u, w)
    y - (x + h * dynamics(y, u, w))
end

dt = Dynamics(discrete_dynamics, nx, nx, nu, nw);
dyn = [dt for t = 1:T-1]
model = DynamicsModel(dyn)

nz = sum([t < T ? dyn[t].nx : dyn[t-1].ny for t = 1:T]) + sum([dyn[t].nu for t = 1:T-1])
nc = num_con(dyn)
nj = num_jac(dyn)
z = rand(nz)
c = zeros(nc)
j = zeros(nj)
x = [zero(z[model.idx.x[t]]) for t = 1:T]
u = [[zero(z[model.idx.u[t]]) for t = 1:T-1]..., zeros(0)]
w = [zeros(nw) for t = 1:T-1]

con = StageConstraint()
cons = [con for t = 1:T]
cons isa StageConstraints
model isa DynamicsModel
obj isa Objective

trajopt = TrajectoryOptimizationProblem(obj, model, cons, x, u, w)

p = Problem(trajopt, nz, nc,
    [-Inf * ones(nz), ones(nz)], [zeros(nc), zeros(nc)],
    sparsity(dyn, model.dim.x, model.dim.u), [], false)

@code_warntype MOI.eval_objective(p, z)
@benchmark MOI.eval_objective($p, $z)
@benchmark MOI.eval_objective_gradient($p, $grad, $z)

@code_warntype MOI.eval_constraint(p, c, z)
@benchmark MOI.eval_constraint($p, $c, $z)

@code_warntype MOI.eval_constraint_jacobian(p, j, z)
@benchmark MOI.eval_constraint_jacobian($p, $j, $z)

sparsity_jacobian(p)
sparsity_hessian_lagrangian(p)

z0 = 0.1 * randn(p.num_var)
MOI.eval_constraint(p, c, z0)

function obje(z)
    sum([ft(z[p.trajopt.model.idx.x[t]], z[p.trajopt.model.idx.u[t]]) for t = 1:T-1]) + fT(z[p.trajopt.model.idx.x[T]], zeros(0))
end
norm(obje(z0) - MOI.eval_objective(p, z0)) < 1.0e-8
grad .= 0.0
MOI.eval_objective_gradient(p, grad, z0)
norm(ForwardDiff.gradient(obje, z0) - grad)

function dyn_con(z)
    vcat([discrete_dynamics(z[p.trajopt.model.idx.x[t+1]], z[p.trajopt.model.idx.x[t]], z[p.trajopt.model.idx.u[t]], p.trajopt.w) for t = 1:T-1]...)
end
norm(c - dyn_con(z0)) < 1.0e-8

dyn_con_jac = zeros(num_con(p.trajopt.model.dyn), num_xuy(p.trajopt.model.dyn))
MOI.eval_constraint_jacobian(p, j, z0)
sp = sparsity(p.trajopt.model.dyn, p.trajopt.model.dim.x, p.trajopt.model.dim.u)
for (i, v) in enumerate(sp)
    dyn_con_jac[v...] = j[i]
end
norm(ForwardDiff.jacobian(dyn_con, z0) - dyn_con_jac) < 1.0e-8

zl = -Inf * ones(p.num_var)
zu = Inf * ones(p.num_var)
cl = zeros(p.num_con)
cu = zeros(p.num_con)

x1 = [0.0; 0.0]
xT = [π; 0.0]
x_interp = linear_interpolation(x1, xT, T)
for (t, idx) in enumerate(p.trajopt.model.idx.x)
    z0[idx] = x_interp[t]
end

zl[p.trajopt.model.idx.x[1]] = x1
zu[p.trajopt.model.idx.x[1]] = x1
zl[p.trajopt.model.idx.x[T]] = xT
zu[p.trajopt.model.idx.x[T]] = xT

nlp_bounds = MOI.NLPBoundsPair.(cl, cu)
block_data = MOI.NLPBlockData(nlp_bounds, p, true)

solver = Ipopt.Optimizer()
solver.options["max_iter"] = 100
solver.options["tol"] = 1.0e-3
solver.options["constr_viol_tol"] = 1.0e-3
# solver.options["print_level"] = mapl
# solver.options["linear_solver"] = "ma57"

z = MOI.add_variables(solver, p.num_var)

for i = 1:p.num_var
    MOI.add_constraint(solver, z[i], MOI.LessThan(zu[i]))
    MOI.add_constraint(solver, z[i], MOI.GreaterThan(zl[i]))
    MOI.set(solver, MOI.VariablePrimalStart(), z[i], z0[i])
end

# Solve the problem
MOI.set(solver, MOI.NLPBlock(), block_data)
MOI.set(solver, MOI.ObjectiveSense(), MOI.MIN_SENSE)

@time MOI.optimize!(solver)


text:
function w_partial_t(w_x, w_y, phi_x, phi_y, psi_x, psi_y, am)

    _a = -(phi_x.*w_y - phi_y.*w_x)
    _b = -(psi_x.*w_y - psi_y.*w_x)
    return _a + am*_b
end

function psi_partial_t(psi_x, psi_y, phi_x, phi_y)

    return -(phi_x.*psi_y - phi_y.*psi_x)
end

function forward_euler(x, x_t, K1, K2, dt, gd::GhostData)

    return (x + dt*(K1*x_t)) ./ (1 .+ dt*diag(K2))
end

function solve_pde(A, x0, b, gd::GhostData)

    An, _ = require_boundary_conditions(A, gd)
    bn = gd.Proj*b

    bn = transpose(An)*bn
    An = transpose(An)*An
    # x = jacobi_preconditioned_gmres(An, gd.Proj*x0, bn, 1, 1)
    x = gd.Proj*x0

    gauss_seidel!(x, An, bn; maxiter=10)
    gmres!(x, An, bn)
    return gd.R*transpose(gd.Proj)*x
end

function solve_diffusion(x, k, Dxx, Dyy, gd::GhostData)

    M = I - k*(Dxx + Dyy)
    return solve_pde(M, x, x, gd)
end

function solve_vorticity_eqn(phi, w, Dxx, Dyy, gd::GhostData)

    return solve_pde(Dxx + Dyy, phi, w, gd)
end

function leapfrog!(w, psi, phi, Dx, Dxx, Dy, Dyy, K1, K2, gd, dt, am, k_psi, k_w)

    # take derivatives
    psi_x = Dx*psi[:,2]
    psi_y = Dy*psi[:,2]
    phi_x = Dx*phi
    phi_y = Dy*phi
    w_x = Dx*w[:,2]
    w_y = Dy*w[:,2]

    w_t = w_partial_t(w_x, w_y, phi_x, phi_y, psi_x, psi_y, am)
    psi_t = psi_partial_t(psi_x, psi_y, phi_x, phi_y)

    # compute predictor step
    psi[:,3] = forward_euler(0.5*(psi[:,1] + psi[:,2]), psi_t, K1, K2, dt, gd)
    w[:,3] = forward_euler(0.5*(w[:,1] + w[:,2]), w_t, K1, K2, dt, gd)

    # compute diffusion
    psi[:,3] = solve_diffusion(psi[:,3], k_psi, Dxx, Dyy, gd)
    w[:,3] = solve_diffusion(w[:,3], k_w, Dxx, Dyy, gd)

    # solve vorticity equation
    phi[:] = solve_vorticity_eqn(phi, w[:,3], Dxx, Dyy, gd)

    # take derivatives again
    psi_x = Dx*psi[:,3]
    psi_y = Dy*psi[:,3]
    phi_x = Dx*phi
    phi_y = Dy*phi
    w_x = Dx*w[:,3]
    w_y = Dy*w[:,3]

    w_t = w_partial_t(w_x, w_y, phi_x, phi_y, psi_x, psi_y, am)
    psi_t = psi_partial_t(psi_x, psi_y, phi_x, phi_y)

    # compute corrector step
    psi[:,3] = forward_euler(psi[:,2], psi_t, K1, K2, dt, gd)
    w[:,3] = forward_euler(w[:,2], w_t, K1, K2, dt, gd)

    # compute diffusion
    psi[:,3] = solve_diffusion(psi[:,3], k_psi, Dxx, Dyy, gd)
    w[:,3] = solve_diffusion(w[:,3], k_w, Dxx, Dyy, gd)

    # solve vorticity equation
    phi[:] = solve_vorticity_eqn(phi, w[:,3], Dxx, Dyy, gd)

    # shift indices
    psi[:,1:2] = psi[:,2:3]
    w[:,1:2] = w[:,2:3]
end


text:
function C1!(champu :: Array{ComplexF64,2},
             champv :: Array{ComplexF64,2},
             m      :: MicMac,
             t      :: Float64,
             fft_u  :: Vector{ComplexF64},
             fft_v  :: Vector{ComplexF64} )

    ftau!(champu, champv, m, t, fft_u, fft_v)

    fft!(champu, 1)
    fft!(champv, 1)

    champu[1, :] .= 0.0
    champv[1, :] .= 0.0

    champu ./= (1im * m.ktau)
    champv ./= (1im * m.ktau)

    ifft!(champu, 1)
    ifft!(champv, 1)

    champu .*= m.epsilon
    champv .*= m.epsilon

    champu .+= transpose(fft_u)
    champv .+= transpose(fft_v)

end


text:
seats = reduce(vcat, permutedims.(collect.(readlines("day11/input.txt"))))

function step_seats!(seats)
    copy_seats = copy(seats)
    function execute_seat_rules(i, j)
        hi_i, hi_j = size(seats)
        function first_visible_seat(i_mod, j_mod, i_inner, j_inner)
            if 1 <= i_inner <= hi_i && 1 <= j_inner <= hi_j
                if copy_seats[i_inner, j_inner] == '.'
                    return first_visible_seat(i_mod, j_mod, i_inner+i_mod, j_inner+j_mod)
                else
                    return copy_seats[i_inner, j_inner]
                end
            end
            return ['.']
        end
        adj_seats = []
        for im in -1:1, jm in -1:1
            if im == 0 && jm == 0
                continue
            end
            push!(adj_seats, first_visible_seat(im, jm, i+im, j+jm))
        end
        if copy_seats[i, j] == 'L' && count(s -> s != '#', adj_seats) == length(adj_seats)
            seats[i,j] = '#'
        elseif copy_seats[i, j] == '#' && count(s -> s == '#', adj_seats) >= 5
            seats[i,j] = 'L'
        end
    end

    for i in 1:size(seats, 1), j in 1:size(seats, 2)
        execute_seat_rules(i, j)
    end
end

old_seats = nothing
while old_seats != seats
    global old_seats = copy(seats)
    step_seats!(seats)
end

println(count(s -> s == '#', seats))

# that was a super mess! never a bad day to code though :)


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    rescale_magnitude(localaniso; r1, r2, r3, clip=[0.05,0.95])
    rescale_magnitude!(localaniso; r1, r2, r3, clip=[0.05,0.95])

Rescale magnitude values of `localaniso` to desired limits using min-max scaling.
`r1`, `r2` and `r3` refer to the ranges along each main axis. `clip` defines
the quantiles to which min and max limits are defined. The default is set to
[0.05,0.95] to avoid outliers influence.

## Example

In the example below, a 3-D local anisotropy object is rescaled. `r1` has all the
values set to 1.0. `r2` and `r3` are rescaled between (0.2, 1.0) and (0.1, 0.5)
respectively. That means that if we associate it to a SphericalVariogram of range=10,
range1 will be equal to 10 at every place. range2 will vary between 2 and 10.
And range3 will vary between 1 and 5.

```julia
example = rescale_magnitude(localaniso3d, r1=1, r2=(0.2,1.0), r3=(0.1,0.5))
```
"""
function rescale_magnitude!(lp::LocalAnisotropy; r1=nothing, r2=nothing,
	                        r3=nothing, clip=[0.05,0.95])
    N = ndims(lp)
    m = lp.magnitude
	rvals = (r1,r2,r3)

    for i in 1:N
		rvals[i] == nothing && continue

        mi = m[i,:]
        qx = quantile(mi, clip)

		if rvals[i] isa Number
			m[i,:] .= rvals[i]
		elseif qx[1] == qx[2]
			throw(ArgumentError("clipped values are unique; inform a single value to r$i"))
		else
	        mi[mi .< qx[1]] .= qx[1]
	        mi[mi .> qx[2]] .= qx[2]
	        mi = (mi .- qx[1]) ./ (qx[2]-qx[1])
	        b1, b2 = rvals[i] isa Number ? (rvals[i],rvals[i]) : rvals[i]
	        m[i,:] .= (b1 .+ (b2-b1) .* mi)
		end
    end
	lp
end

function rescale_magnitude(lpo::LocalAnisotropy; kwargs...)
    lp = deepcopy(lpo)
	rescale_magnitude!(lp; kwargs...)
end

"""
    reference_magnitude(localaniso, axis)
    reference_magnitude!(localaniso, axis)

Rescale the magnitudes of `localaniso` by setting `axis` as a unit reference.
`axis` can be :X, :Y or :Z.

## Example

For example, consider a point with magnitude equal to [1.0, 0.2]. After scaling
it with a ExponentialVariogram of range=10m, the range in X would be 10m and in
Y would be 2m at this point. If we first call `reference_magnitude!(localaniso, :Y)`,
the magnitude will turn to [5.0, 1.0], and after scaling to the same previous
variogram, the range in X would be 50m and in Y would be 10m.
"""
function reference_magnitude!(lpar::LocalAnisotropy, ax::Symbol)
  # rescale magnitude according to reference axis
  m = lpar.magnitude
  ref = m[iaxis(ax),:]
  for i in 1:size(m, 1)
    m[i,:] ./= ref
  end
  lpar
end

function reference_magnitude(lpar::LocalAnisotropy, ax::Symbol)
  lp = deepcopy(lpar)
  reference_magnitude!(lp, ax)
end

"""
    to_vtk(vtkfile, coords, localaniso; dir=:ellips, magnitude=:ranges)

Export local anisotropies `localaniso` at `coords` to VTK format. It export local
ellipsoids by default. If `dir` is set to `:X`, `:Y`, `:Z` or `:XYZ`, local
vectors/arrows are exported on these axes. The magnitude can be expressed as
`:ranges` (default) or `:ratios`. In the case of `:ratios`, `ratio1=range2/range1`
and `ratio2=range3/range1`. Instead of an array of coordinates, `coords` can
also be just the georeferenced spatial object. The output `.vtu` file can be
loaded in Paraview or similars and processed as Glyphs to plot the directions.
If `dir = :ellips`, TensorGlyphs must be used for proper visualization (need to
disable extract eigenvalues).
"""
function to_vtk(vtkfile, coords::AbstractArray, lpars::LocalAnisotropy;
	            dir=:ellips, magnitude=:ranges)
	@assert dir in [:ellips, :XYZ, :X, :Y, :Z] "`dir` must be :ellips, :XYZ :X, :Y or :Z"
	@assert magnitude in [:ratios, :ranges] "`magnitude` must be :ratios or :ranges"

	elp = (dir == :ellips)
	dim = size(coords,1)
	n   = size(coords,2)
	axs = dir in [:XYZ, :ellips] ? (1:dim) : [iaxis(dir)]
	xyz = dim == 2 ? vcat(coords, zeros(Float64,1,n)) : coords
	vtx = [MeshCell( VTKCellTypes.VTK_VERTEX, [i]) for i in 1:n]
	ijk = elp ? zeros(Float64, 9, n, 1, 1) : [zeros(Float64, 3, n, 1, 1) for i in axs]
	zx  = (dim == 2 && elp) ? minimum(lpars.magnitude)/10 : 0.0

	# put data into vtk structure format
	for x in 1:n
		if elp
			mag = lpars.magnitude[:,x]
			dim == 2 && push!(mag, zx)
			Λ = Diagonal(SVector{3}(mag))
			P = rotmat(lpars, x)' * Λ
			for v in 1:length(P)
				ijk[v, x, 1, 1] = P[v]
			end
		else
			dcm = rotmat(lpars, x)
	        f = dcm[3,3] < 0 ? -1 : 1
			for d in axs
				ijk[d][1, x, 1, 1] = f*dcm[d,1]
				ijk[d][2, x, 1, 1] = f*dcm[d,2]
				ijk[d][3, x, 1, 1] = f*dcm[d,3]
			end
		end
	end

	# write vtk file
	outfiles = vtk_grid(vtkfile,xyz,(vtx)) do vtk
		if elp
			vtk["ellips"] = ijk
		else
			for d in axs
				vtk["dir$d"] = ijk[d]
			end
		end
		for d in 1:dim
			if magnitude == :ranges
				vtk["range$d"] = lpars.magnitude[d,:]
			elseif d > 1
				vtk["ratio$(d-1)"] = lpars.magnitude[d,:] ./ lpars.magnitude[1,:]
			end
		end
	end
end

function to_vtk(vtkfile, D::SpatialData, lpars::LocalAnisotropy; kwargs...)
	coords = reduce(hcat, [coordinates(centroid(D,x)) for x in 1:nelements(D)])
	to_vtk(vtkfile, coords, lpars; kwargs...)
end

toqmat(lp) = [qmat(rotation(lp,i),magnitude(lp,i)) for i in 1:nvals(lp)]


text:
module CollidingSphereWithBunnies

using  ModiaLang
import Modia3D
using  Modia3D.ModiaInterface

filenameBunny = joinpath(Modia3D.path, "objects", "bunny", "bunny.obj")

# Properties
mat1 = "Red"

massAndGeoSphere = MassPropertiesFromShape()

# Objects3D
ConvexPartitions = Model(
    gravField = UniformGravityField(g=9.81, n=[0, -1, 0]),
    world = Object3D(feature=Scene(gravityField=:gravField)),

    sphere = Object3D(feature=Solid(shape=Sphere(diameter=0.2),
                      visualMaterial=mat1, solidMaterial="Steel",
                      massProperties=massAndGeoSphere, collision=true)),
    free = FreeMotion(obj1=:world, obj2=:sphere, r=Var(init=[0.0, 2.0, 0.0])),

    bunnySolid = Object3D(parent=:world, feature=Solid(shape =
        FileMesh(filename=filenameBunny, scale=[0.1, 0.1, 0.1]), solidMaterial="Steel", collision = true )),

    bunnyPartitionSolid = Object3D(parent=:world, translation=[0.5, -2.0, 2.0], feature=Solid(shape =
        FileMesh(filename = filenameBunny, scale=[0.1, 0.1, 0.1], convexPartition=true), solidMaterial="Steel", collision=true)),
)

convexPartitions = @instantiateModel(buildModia3D(ConvexPartitions), unitless=true, log=false, logStateSelection=false, logCode=false)

stopTime = 1.3
requiredFinalStates = [1.6852536030611518, -1.9601748772369065, 1.553773946072609, 4.576339600748371, -3.635822903848681, -1.8278191918133224, 16.048785711292467, 0.11321559138403894, 8.101411301098336, 28.75260252729643, 8.402335107079937, 37.18798381568381]
@time simulate!(convexPartitions, stopTime=stopTime, log=true, logStates=false, logEvents=true, requiredFinalStates=requiredFinalStates)

end


text:
"""
    mutable struct OuterControl{
        A <: ActivePowerControl,
        R <: ReactivePowerControl
    } <: DynamicInverterComponent
        active_power::A
        reactive_power::R
        ext::Dict{String, Any}
        states::Vector{Symbol}
        n_states::Int
    end
Parameters of a Outer-Loop controller using a active power controller and a reactive power droop controller.

# Arguments
- `A <: ActivePowerControl`: Active power controller (typically droop or virtual inertia).
- `R <: ReactivePowerControl`: Reactive power controller (typically droop).
- `ext::Dict{String, Any}`
- `states::Vector{Symbol}`: Vector of states (will depend on the components).
- `n_states::Int`: Number of states (will depend on the components).
"""
mutable struct OuterControl{A <: ActivePowerControl, R <: ReactivePowerControl} <:
               DynamicInverterComponent
    active_power::A
    reactive_power::R
    ext::Dict{String, Any}
    states::Vector{Symbol}
    n_states::Int
end

function OuterControl(
    active_power::A,
    reactive_power::R,
    ext = Dict{String, Any}(),
) where {A <: ActivePowerControl, R <: ReactivePowerControl}
    OuterControl(
        active_power,
        reactive_power,
        ext,
        vcat(active_power.states, reactive_power.states),
        active_power.n_states + reactive_power.n_states,
    )
end

function OuterControl(;
    active_power,
    reactive_power,
    ext = Dict{String, Any}(),
    states = nothing,
    n_states = nothing,
)
    if states === nothing
        @assert n_states === nothing
        return OuterControl(active_power, reactive_power, ext)
    end
    @assert n_states !== nothing
    return OuterControl(active_power, reactive_power, ext, states, n_states)
end

get_active_power(value::OuterControl) = value.active_power
get_reactive_power(value::OuterControl) = value.reactive_power
get_ext(value::OuterControl) = value.ext
get_states(value::OuterControl) = value.states
get_n_states(value::OuterControl) = value.n_states


text:
dark_palette = [
    colorant"#E06C75", # red
    colorant"#98C379", # green
    colorant"#61AFEF",  # blue
    colorant"#C678DD",  # purple
    colorant"#E5C07B", # yellow
    colorant"#56B6C2", # teal
    colorant"#ABB2BF" # grey
]
dark_bg = colorant"#373E4B"

_themes[:onedark] = PlotTheme(
    bg = dark_bg,
    bginside = colorant"#282C34",
    fg = colorant"#ADB2B7",
    fgtext = colorant"#B1BED6",
    fgguide = colorant"#B1BED6",
    fglegend = colorant"#B1BED6",
    legendfontcolor = colorant"#B1BED6",
    legendtitlefontcolor = colorant"#B1BED6",
    titlefontcolor = colorant"#B1BED6",
    palette = expand_palette(dark_bg, dark_palette; lchoices = [57], cchoices = [100]),
    colorgradient = :fire
)


text:
include("../foamLia.jl")

using foamLia

# okay, test implementation
# case = OpenFoam("/users/a/r/areagan/scratch/run/juliatest")
# caseFolder = "/users/a/r/areagan/OpenFOAM/areagan-2.2.1/run/juliabase"
caseFolder = "/users/a/r/areagan/scratch/run/BCTest-260-370"
# caseFolder = "/users/a/r/areagan/scratch/run/2014-10-15-small-mesh-long-runs/400"
case = OpenFoam(caseFolder)

# manipulate parameters
case.controlDict["endTime"] = 1
case.controlDict["deltaT"] = .1
case.controlDict["writeInterval"] = .1
# case.T["..."] = replace(case.T["..."],Text=340,Text=Ttop)
# # ["boundaryField"]["bottominside"]["variables"] = "\"Text=340;hc=225;gradT=(Text-T)*hc\""

baseCase = "/users/a/r/areagan/OpenFOAM/areagan-2.2.1/run/juliabase"
# initCase(case,baseCase)

# (o::OpenFoam,d::OrderedDict,name::String,location::String)
# writeVolScalarField(test,defaultT,"T","0")

# cd("/users/a/r/areagan/scratch/run/juliatest")

# run(`./Allrun`)

# readMesh(case)

timeSaves = findTimes(case)
println(timeSaves)

# T = readVar(case,"0.24","T")

# println(size(T))

# Tbottom = readVarSpec(case,"0.24","T",[1,10,100,1000,1600])

# println(Tbottom)

for t in timeSaves
    if t>0
	continue
    end
end


text:
getindex(a::Number, ::Block{0}) = a

axes(A::AbstractTriangular{<:Any,<:AbstractBlockMatrix}) = axes(parent(A))
axes(A::HermOrSym{<:Any,<:AbstractBlockMatrix}) = axes(parent(A))

Base.print_matrix_row(io::IO,
        X::Union{AbstractTriangular{<:Any,<:AbstractBlockMatrix},
                 Symmetric{<:Any,<:AbstractBlockMatrix},
                 Hermitian{<:Any,<:AbstractBlockMatrix},
                 Adjoint{<:Any,<:AbstractBlockMatrix},
                 Transpose{<:Any,<:AbstractBlockMatrix}}, A::Vector,
        i::Integer, cols::AbstractVector, sep::AbstractString) =
        _blockarray_print_matrix_row(io, X, A, i, cols, sep)


text:
export FDTD, fdtd, fdtd!

immutable FDTD{T,N}
	p0::Array{T,N}
	p1::Array{T,N}
	p2::Array{T,N}
	d::Tuple{T,T,T,T}
	q::Vector{T}
	G::Array{UInt8,3}
	env::FDTDEnv

	function FDTD(p0,p1,p2,geo)
		a,b,c,λ = geo.env.scheme.a, geo.env.scheme.b, geo.env.scheme.c, geo.env.scheme.λ
		d = (λ^2*(1-4*a+4*b),
      		     λ^2*(a-2*b),
                     λ^2*b,
                     2*(1+λ^2*(-3+6*a-4*b)) )
		q = λ./(geo.ξ*2.0)
		new(p0,p1,p2,d,q,geo.G,geo.env)
	end

	FDTD(p0,p1,p2,d,q,G,env) = new(p0,p1,p2,d,q,G,env)

end

"""
`FDTD(S::Type,geo::AbstractGeometry)`

Returns a `FDTD` object that can be used in `fdtd` and `fdtd!`.
`S` must be an `AbstractFloat` and determinates the numerical precision of the FDTD simulation.
"""
function FDTD(S::Type,geo::AbstractGeometry)

	p0 = zeros(S,geo.Nx,geo.Ny,geo.Nz)
	p1 = zeros(S,geo.Nx,geo.Ny,geo.Nz)
	p2 = zeros(S,geo.Nx,geo.Ny,geo.Nz)
	FDTD{S,3}(p0,p1,p2,geo)
end

FDTD(geo::AbstractGeometry) = FDTD(Float64, geo)

function setIC!{T,N}(f::FDTD{T,N},p0::Array{T,N},p1::Array{T,N})
	f.p0 .= p0
	f.p1 .= p1
end

function resetIC!{T,N}(f::FDTD{T,N})
	f.p0 .= 0.0
	f.p1 .= 0.0
end

fdtd(xr::Tuple{Int,Int,Int}, xs::Tuple{Int,Int,Int}, args...) =
fdtd([xr], [xs], args...)

fdtd(xr::Vector{Tuple{Int,Int,Int}}, xs::Tuple{Int,Int,Int}, args...) =
fdtd(xr, [xs], args...)

fdtd(xs::Tuple{Int,Int,Int}, args...) =
fdtd([xs], args...)

"""
* `p = fdtd([xr, xs], Nt, geo, s, [p0ic, p1ic])`
  * `xs`  : microphone positions
  * `xr`  : source positions (must be a `Tuple{Int,Int,Int}` or and `Vector{Tuple}` like `xs`)
  * `Nt`  : number of time samples `Int`
  * `geo` : geometry specification, must be a `AbstractGeometry` or a `FDTD` object (the latter saves some memory allocation)
  * `s`   :
    * if `xs` is given must be an `Array` of `size(s) = (Nt, length(xs))` where `s[:,k]` is the `k`th source signal positioned at `xs[k]`
    * otherwise must be a an `Array` of size `(size(geo)...,Nt)`
  * `p0ic` and `p1ic`: `Array` of `size(geo)` consisting of the initial conditions

  * reurns `p`:
    * if `xr` is given `p` is an `Array` of size `(Nt, length(xr))` where the `p[:,k]` is the `k`th mic signal positioned at `xr[k]`
    * otherwise `p` will be an `Array` os size `(size(geo)...,Nt)`

"""
function fdtd{T}(xr::Vector{Tuple{Int,Int,Int}}, xs::Vector{Tuple{Int,Int,Int}}, Nt::Int,
		 geo::AbstractGeometry,
		 s::Array{T})
	f = FDTD(eltype(s),geo)
	return fdtd(xr,xs,Nt,f,s)
end

function fdtd{T}(xr::Vector{Tuple{Int,Int,Int}}, xs::Vector{Tuple{Int,Int,Int}}, Nt::Int,
		 f::FDTD{T,3}, s::Array{T})
	size(s,2) != length(xs) && throw(ArgumentError("size(s,2) must be equal to length(xs)"))

	#initialize output
	p_out = Array{T}(Nt,length(xr))
	fdtd!(p_out, xr, xs, Nt, f, s)
	return p_out
end

"""
`p = fdtd!(p::Array, xr, xs, Nt, f::FDTD, s, [p0ic, p1ic])`

In-place version of `fdtd` overwrites the `Array` `p`. See `fdtd` for more details.
"""
function fdtd!{T}(p_out::Array{T},
		  xr::Vector{Tuple{Int,Int,Int}}, xs::Vector{Tuple{Int,Int,Int}},Nt::Int,
		  f::FDTD{T,3}, s::Array{T})
	resetIC!(f)
	p2 = f.p2
	p1 = f.p1
	p0 = f.p0
	for n = 1:Nt
		fdtdKernel!(n,p_out,p0,p1,p2,xr,xs,s,f.d[1],f.d[2],f.d[3],f.d[4],f.q,f.G)
		p2, p1, p0 = p0, p2, p1
	end
	return p_out
end

## full with xr, null IC
function fdtd{T}(xs::Vector{Tuple{Int,Int,Int}}, Nt::Int,
		 geo::AbstractGeometry,
		 s::Array{T})
	f = FDTD(eltype(s),geo)
	return fdtd(xs,Nt,f,s)
end

function fdtd{T}(xs::Vector{Tuple{Int,Int,Int}}, Nt::Int,
		 f::FDTD{T,3}, s::Array{T})
	size(s,2) != length(xs) && throw(ArgumentError("size(s,2) must be equal to length(xs)"))

	#initialize output
	p_out = Array{T}(size(f.p0)...,Nt)
	fdtd!(p_out, xs, Nt, f, s)
	return p_out
end

function fdtd!{T}(p_out::Array{T,4},
		  xs::Vector{Tuple{Int,Int,Int}},Nt::Int,
		  f::FDTD{T,3}, s::Array{T})
	resetIC!(f)
	p2 = f.p2
	p1 = f.p1
	p0 = f.p0
	for n = 1:Nt
		fdtdKernel!(n,p_out,p0,p1,p2,xs,s,f.d[1],f.d[2],f.d[3],f.d[4],f.q,f.G)
		p2, p1, p0 = p0, p2, p1
	end
	return p_out
end

## xr with source everywhere, null IC
function fdtd{T}(xr::Vector{Tuple{Int,Int,Int}},Nt::Int,
		 geo::AbstractGeometry,
		 s::Array{T,4})
	f = FDTD(eltype(s),geo)
	return fdtd(xr,Nt,f,s)
end

function fdtd{T}(xr::Vector{Tuple{Int,Int,Int}},Nt::Int,
		 f::FDTD{T,3}, s::Array{T,4})

	#initialize output
	p_out = Array{T}(Nt,length(xr))
	fdtd!(p_out, xr, Nt, f, s)
	return p_out
end

function fdtd!{T}(p_out::Array{T},
		  xr::Vector{Tuple{Int,Int,Int}}, Nt::Int,
		  f::FDTD{T,3}, s::Array{T,4})
	resetIC!(f)
	p2 = f.p2
	p1 = f.p1
	p0 = f.p0
	for n = 1:Nt
		fdtdKernel!(n,p_out,p0,p1,p2,xr,s,f.d[1],f.d[2],f.d[3],f.d[4],f.q,f.G)
		p2, p1, p0 = p0, p2, p1
	end
	return p_out
end

## full with source everywhere, null IC
function fdtd{T}(Nt::Int,
		 geo::AbstractGeometry,
		 s::Array{T,4})
	f = FDTD(eltype(s),geo)
	return fdtd(Nt,f,s)
end

function fdtd{T}(Nt::Int,
		 f::FDTD{T,3}, s::Array{T,4})

	#initialize output
	p_out = Array{T}(size(f.p0)...,Nt)
	fdtd!(p_out, Nt, f, s)
	return p_out
end

function fdtd!{T}(p_out::Array{T,4},
		  Nt::Int,
		  f::FDTD{T,3}, s::Array{T,4})
	resetIC!(f)
	p2 = f.p2
	p1 = f.p1
	p0 = f.p0
	for n = 1:Nt
		fdtdKernel!(n,p_out,p0,p1,p2,s,f.d[1],f.d[2],f.d[3],f.d[4],f.q,f.G)
		p2, p1, p0 = p0, p2, p1
	end
	return p_out
end

## xr and xs with IC

function fdtd{T}(xr::Vector{Tuple{Int,Int,Int}}, xs::Vector{Tuple{Int,Int,Int}}, Nt::Int,
		 geo::AbstractGeometry,
		 s::Array{T},
		 p0ic::Array{T,3}, p1ic::Array{T,3})
	f = FDTD(eltype(s),geo)
	return fdtd(xr,xs,Nt,f,s,p0ic,p1ic)
end

function fdtd{T}(xr::Vector{Tuple{Int,Int,Int}}, xs::Vector{Tuple{Int,Int,Int}}, Nt::Int,
		 f::FDTD{T,3}, s::Array{T},
		 p0ic::Array{T,3}, p1ic::Array{T,3})
	size(s,2) != length(xs) && throw(ArgumentError("size(s,2) must be equal to length(xs)"))

	#initialize output
	p_out = Array{T}(Nt,length(xr))
	fdtd!(p_out, xr, xs, Nt, f, s, p0ic, p1ic)
	return p_out
end

function fdtd!{T}(p_out::Array{T},
		  xr::Vector{Tuple{Int,Int,Int}}, xs::Vector{Tuple{Int,Int,Int}},Nt::Int,
		  f::FDTD{T,3}, s::Array{T},
		  p0ic::Array{T,3}, p1ic::Array{T,3})
	setIC!(f,p0ic,p1ic)
	p2 = f.p2
	p1 = f.p1
	p0 = f.p0
	for n = 1:Nt
		fdtdKernel!(n,p_out,p0,p1,p2,xr,xs,s,f.d[1],f.d[2],f.d[3],f.d[4],f.q,f.G)
		p2, p1, p0 = p0, p2, p1
	end
	return p_out
end

#full soundfield & IC
function fdtd{T}(xs::Vector{Tuple{Int,Int,Int}}, Nt::Int,
		 geo::AbstractGeometry,
		 s::Array{T},
		 p0ic::Array{T,3}, p1ic::Array{T,3})
	f = FDTD(eltype(s),geo)
	return fdtd(xs,Nt,f,s,p0ic,p1ic)
end

function fdtd{T}(xs::Vector{Tuple{Int,Int,Int}}, Nt::Int,
		 f::FDTD{T,3}, s::Array{T},
		 p0ic::Array{T,3}, p1ic::Array{T,3})
	size(s,2) != length(xs) && throw(ArgumentError("size(s,2) must be equal to length(xs)"))

	#initialize output
	p_out = Array{T}(size(f.p0)...,Nt)
	fdtd!(p_out, xs, Nt, f, s, p0ic, p1ic)
	return p_out
end

function fdtd!{T}(p_out::Array{T,4},
		  xs::Vector{Tuple{Int,Int,Int}},Nt::Int,
		  f::FDTD{T,3}, s::Array{T},
		  p0ic::Array{T,3}, p1ic::Array{T,3})
	setIC!(f,p0ic,p1ic)
	p2 = f.p2
	p1 = f.p1
	p0 = f.p0
	for n = 1:Nt
		fdtdKernel!(n,p_out,p0,p1,p2,xs,s,f.d[1],f.d[2],f.d[3],f.d[4],f.q,f.G)
		p2, p1, p0 = p0, p2, p1
	end
	return p_out
end

## xr with source everywhere, & IC
function fdtd{T}(xr::Vector{Tuple{Int,Int,Int}},Nt::Int,
		 geo::AbstractGeometry,
		 s::Array{T,4},
		 p0ic::Array{T,3}, p1ic::Array{T,3})
	f = FDTD(eltype(s),geo)
	return fdtd(xr,Nt,f,s,p0ic,p1ic)
end

function fdtd{T}(xr::Vector{Tuple{Int,Int,Int}},Nt::Int,
		 f::FDTD{T,3}, s::Array{T,4},
		 p0ic::Array{T,3}, p1ic::Array{T,3})

	#initialize output
	p_out = Array{T}(Nt,length(xr))
	fdtd!(p_out, xr, Nt, f, s,p0ic,p1ic)
	return p_out
end

function fdtd!{T}(p_out::Array{T},
		  xr::Vector{Tuple{Int,Int,Int}}, Nt::Int,
		  f::FDTD{T,3}, s::Array{T,4},
		  p0ic::Array{T,3}, p1ic::Array{T,3})
	setIC!(f,p0ic,p1ic)
	p2 = f.p2
	p1 = f.p1
	p0 = f.p0
	for n = 1:Nt
		fdtdKernel!(n,p_out,p0,p1,p2,xr,s,f.d[1],f.d[2],f.d[3],f.d[4],f.q,f.G)
		p2, p1, p0 = p0, p2, p1
	end
	return p_out
end

## full with source everywhere, & IC
function fdtd{T}(Nt::Int,
		 geo::AbstractGeometry,
		 s::Array{T,4},
		 p0ic::Array{T,3}, p1ic::Array{T,3})
	f = FDTD(eltype(s),geo)
	return fdtd(Nt,f,s,p0ic,p1ic)
end

function fdtd{T}(Nt::Int,
		 f::FDTD{T,3}, s::Array{T,4},
		 p0ic::Array{T,3}, p1ic::Array{T,3})

	#initialize output
	p_out = Array{T}(size(f.p0)...,Nt)
	fdtd!(p_out, Nt, f, s,p0ic,p1ic)
	return p_out
end

function fdtd!{T}(p_out::Array{T,4},
		  Nt::Int,
		  f::FDTD{T,3}, s::Array{T,4},
		  p0ic::Array{T,3}, p1ic::Array{T,3})
	setIC!(f,p0ic,p1ic)
	p2 = f.p2
	p1 = f.p1
	p0 = f.p0
	for n = 1:Nt
		fdtdKernel!(n,p_out,p0,p1,p2,s,f.d[1],f.d[2],f.d[3],f.d[4],f.q,f.G)
		p2, p1, p0 = p0, p2, p1
	end
	return p_out
end

##walls
include("kernels.jl")

function air!{T,N}(l::Int,m::Int,i::Int,
		   d1::T,d2::T,d3::T,d4::T,
		   p0::Array{T,N},p1::Array{T,N},p2::Array{T,N})
	sum1,sum2,sum3,sum4 = 0.,0.,0.,0.
	if d1 != 0.
		sum1 += p1[l+1,m,i]
		sum1 += p1[l-1,m,i]
		sum1 += p1[l,m+1,i]
		sum1 += p1[l,m-1,i]
		sum1 += p1[l,m,i+1]
		sum1 += p1[l,m,i-1]
		sum1 *= d1
	end
	if d2 != 0.
		sum2 += p1[l+1,m+1,i]
		sum2 += p1[l+1,m-1,i]
		sum2 += p1[l+1,m,i+1]
		sum2 += p1[l+1,m,i-1]
		sum2 += p1[l,m+1,i+1]
		sum2 += p1[l,m+1,i-1]
		sum2 += p1[l,m-1,i+1]
		sum2 += p1[l,m-1,i-1]
		sum2 += p1[l-1,m+1,i]
		sum2 += p1[l-1,m-1,i]
		sum2 += p1[l-1,m,i+1]
		sum2 += p1[l-1,m,i-1]
		sum2 *= d2
	end
	if d3 != 0.
		sum3 += p1[l+1,m+1,i+1]
		sum3 += p1[l+1,m-1,i+1]
		sum3 += p1[l+1,m+1,i-1]
		sum3 += p1[l+1,m-1,i-1]
		sum3 += p1[l-1,m+1,i+1]
		sum3 += p1[l-1,m-1,i+1]
		sum3 += p1[l-1,m+1,i-1]
		sum3 += p1[l-1,m-1,i-1]
		sum3 *= d3
	end
	if d4 != 0.
		sum4 += d4*p1[l,m,i]
	end
	p2[l,m,i]  = sum1+sum2+sum3+sum4-p0[l,m,i]
end

##adding source

function air!{T,N}(n::Int,l::Int,m::Int,i::Int,
		   d1::T,d2::T,d3::T,d4::T,
		   p0::Array{T,N},p1::Array{T,N},p2::Array{T,N},s::Array{T,4})
	sum1,sum2,sum3,sum4 = 0.,0.,0.,0.
	if d1 != 0.
		sum1 += p1[l+1,m,i]
		sum1 += p1[l-1,m,i]
		sum1 += p1[l,m+1,i]
		sum1 += p1[l,m-1,i]
		sum1 += p1[l,m,i+1]
		sum1 += p1[l,m,i-1]
		sum1 *= d1
	end
	if d2 != 0.
		sum2 += p1[l+1,m+1,i]
		sum2 += p1[l+1,m-1,i]
		sum2 += p1[l+1,m,i+1]
		sum2 += p1[l+1,m,i-1]
		sum2 += p1[l,m+1,i+1]
		sum2 += p1[l,m+1,i-1]
		sum2 += p1[l,m-1,i+1]
		sum2 += p1[l,m-1,i-1]
		sum2 += p1[l-1,m+1,i]
		sum2 += p1[l-1,m-1,i]
		sum2 += p1[l-1,m,i+1]
		sum2 += p1[l-1,m,i-1]
		sum2 *= d2
	end
	if d3 != 0.
		sum3 += p1[l+1,m+1,i+1]
		sum3 += p1[l+1,m-1,i+1]
		sum3 += p1[l+1,m+1,i-1]
		sum3 += p1[l+1,m-1,i-1]
		sum3 += p1[l-1,m+1,i+1]
		sum3 += p1[l-1,m-1,i+1]
		sum3 += p1[l-1,m+1,i-1]
		sum3 += p1[l-1,m-1,i-1]
		sum3 *= d3
	end
	if d4 != 0.
		sum4 += d4*p1[l,m,i]
	end
	p2[l,m,i]  = sum1+sum2+sum3+sum4-p0[l,m,i]+s[l,m,i,n]
end

##walls
include("walls.jl")
##edges
include("edges.jl")
##corners
include("corners.jl")

max_stars_repo_name: nantonel/JuFdtd

max_issues_repo_name: nantonel/JuFdtd

max_forks_repo_name: nantonel/JuFdtd


text:
using Documenter, FlowFarm

makedocs(;
    modules=[FlowFarm],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
        "Tutorial" => "Tutorial.md",
        "How-To" => "How_to.md",
        "Explanation" => "Explanation.md",
        "Reference" => "Reference.md"
    ],
    repo="https://github.com/byuflowlab/FlowFarm.jl/blob/{commit}{path}#L{line}",
    sitename="FlowFarm.jl",
    authors="PJ Stanley <stanley_andrewpj@yahoo.com>",
    assets=String[],
)

deploydocs(
    repo = "github.com/byuflowlab/FlowFarm.jl.git",
)


text:
using Distributions
using Test

# to make sure that subtypes provides the required behavoir without having to add
# a dependency to InteractiveUtils
function _subtypes(m::Module, x::Type, sts=Base.IdSet{Any}(), visited=Base.IdSet{Module}())
    push!(visited, m)
    xt = Base.unwrap_unionall(x)
    if !isa(xt, DataType)
        return sts
    end
    xt = xt::DataType
    for s in names(m, all = true)
        if isdefined(m, s) && !Base.isdeprecated(m, s)
            t = getfield(m, s)
            if isa(t, DataType)
                t = t::DataType
                if t.name.name === s && supertype(t).name == xt.name
                    ti = typeintersect(t, x)
                    ti != Base.Bottom && push!(sts, ti)
                end
            elseif isa(t, UnionAll)
                t = t::UnionAll
                tt = Base.unwrap_unionall(t)
                isa(tt, DataType) || continue
                tt = tt::DataType
                if tt.name.name === s && supertype(tt).name == xt.name
                    ti = typeintersect(t, x)
                    ti != Base.Bottom && push!(sts, ti)
                end
            elseif isa(t, Module)
                t = t::Module
                in(t, visited) || _subtypes(t, x, sts, visited)
            end
        end
    end
    return sts
end

function _subtypes_in(mods::Array, x::Type)
    if !isabstracttype(x)
        # Fast path
        return Type[]
    end
    sts = Base.IdSet{Any}()
    visited = Base.IdSet{Module}()
    for m in mods
        _subtypes(m, x, sts, visited)
    end
    return sort!(collect(sts), by=string)
end

get_subtypes(m::Module, x::Type) = _subtypes_in([m], x)

get_subtypes(x::Type) = _subtypes_in(Base.loaded_modules_array(), x)

dists = get_subtypes(UnivariateDistribution)
filter!(x -> hasmethod(x, ()), dists)
filter!(x -> isbounded(x()), dists)

@testset "bound checking $dist" for dist in dists
    d = dist()
    lb,ub = float.(extrema(support(d)))
    lb = prevfloat(lb)
    ub = nextfloat(ub)
    @test iszero(cdf(d, lb))
    @test isone(cdf(d, ub))

    lb_lcdf = logcdf(d,lb)
    @test isinf(lb_lcdf) & (lb_lcdf < 0)
    @test iszero(logcdf(d, ub))

    @test isone(ccdf(d, lb))
    @test iszero(ccdf(d, ub))

    ub_lccdf = logccdf(d,ub)
    @test isinf(ub_lccdf) & (ub_lccdf < 0)
    @test iszero(logccdf(d, lb))

    @test iszero(pdf(d, lb))
    @test iszero(pdf(d, ub))

    lb_lpdf = logpdf(d, lb)
    @test isinf(lb_lpdf) & (lb_lpdf < 0)
    ub_lpdf = logpdf(d, ub)
    @test isinf(ub_lpdf) & (ub_lpdf < 0)
end


text:
module Day15

using ..ReTest, ..DocStringExtensions

import ..solve, ..PuzzleInfo, ..TEST_DATA_DIRECTORY

abstract type Day15PuzzleInfo <: PuzzleInfo end

"""
Day 15 part 1 puzzle information.

$(FIELDS)
"""
struct Part1PuzzleInfo <: Day15PuzzleInfo
    """Path to input data file."""
    datapath::String
end

"""
Day 15 part 2 puzzle information.

$(FIELDS)
"""
struct Part2PuzzleInfo <: Day15PuzzleInfo
    """Path to input data file."""
    datapath::String
end

"""
    $(FUNCTIONNAME)(input)

Read two-dimensional integer (1:9) risk level data from input stream `input`.
"""
function readrisklevels(input::IO)
    reduce(vcat, [parse(Int, c) for c in line]' for line in eachline(input))
end

"""
    $(FUNCTIONNAME)(index, shape)

Get iterator over valid neighbour indices to `index` for 2D array with size `shape`.
"""
function neighboursof(index::CartesianIndex{2}, shape::Tuple{Integer, Integer})
    i_offsets(i) = i == 1 ? (1,) : (i == shape[1] ? (-1,) : (-1, 1))
    j_offsets(j) = j == 1 ? (1,) : (j == shape[2] ? (-1,) : (-1, 1))
    Iterators.flatten((
        (index + CartesianIndex(i, 0) for i in i_offsets(index[1])),
        (index + CartesianIndex(0, j) for j in j_offsets(index[2]))
    ))
end

"""
    $(FUNCTIONNAME)(risklevels)

Expand two-dimensional risk levels array by tiling right and downward 5 times, on
each horizontal or vertical shift incrementing the risk levels by 1 and wrapping to 1-9.
"""
function expandrisklevels(risklevels::Array{T, 2}) where T
    shape = size(risklevels)
    expandedrisklevels = Array{T, 2}(undef, size(risklevels) .* 5)
    wrap(x) = mod(x - 1, 9) + 1
    for i in 0:4
        for j in 0:4
            expandedrisklevels[
                (shape[1] * i + 1):(shape[1] * (i + 1)),
                (shape[2] * j + 1):(shape[2] * (j + 1)),
            ] = wrap.(risklevels .+ (i + j))
        end
    end
    expandedrisklevels
end

"""
    $(FUNCTIONNAME)(values, condition)
Get index of minimum value of `values[condition]` within `values`.
"""
argminwhere(values, condition) = findall(condition)[argmin(values[condition])]

function solve(info::Day15PuzzleInfo)
    risklevels = open(readrisklevels, info.datapath)
    isa(info, Part2PuzzleInfo) && (risklevels = expandrisklevels(risklevels))
    visited = falses(size(risklevels))
    totalrisk = fill(typemax(Int), size(risklevels))
    current = CartesianIndex(1, 1)
    target = CartesianIndex(size(risklevels))
    totalrisk[current] = 0
    while !(current == target)
        visited[current] = true
        for neighbour in neighboursof(current, size(risklevels))
            if !visited[neighbour]
                altrisk = totalrisk[current] + risklevels[neighbour]
                altrisk < totalrisk[neighbour] && (totalrisk[neighbour] = altrisk)
            end
        end
        current = argminwhere(totalrisk, .!visited)
    end
    totalrisk[target]
end

@testset begin
    @test solve(Part1PuzzleInfo(joinpath(TEST_DATA_DIRECTORY, "day15.txt"))) == 40
    @test solve(Part2PuzzleInfo(joinpath(TEST_DATA_DIRECTORY, "day15.txt"))) == 315
end

end # module


text:
using Basic
using Test

@testset "Basic.jl" begin
    # Write your tests here.
end


text:
# Property: a thing that changes how other things are drawn.

import Base.fill

export stroke, fill, linewidth, font, fontsize, visible, clip, opacity, svgid,
       svgclass, svglink, onactive, onclick, onfocusin, onfocusout, onload,
       onmousedown, onmousemove, onmouseout, onmouseover, onmouseup, svgmask,
       svgdefmask, svgembed, svgattribute, d3embed, d3hook

# A property primitive is something can be directly applied.
abstract PropertyPrimitive

# A property is a (possibly empty) sequence of property primitives. The set of
# Properties with the combine operator forms a monoid (specifically, it can be
# thought of as a free monoid). The empty_property constant is the identity
# element.
abstract Property

# An empty (i.e., nop) property which form the identity element of the Property
# monoid.
type EmptyProperty <: Property end
const empty_property = EmptyProperty()

copy(p::EmptyProperty) = p

# A non-empty sequence of property primitives.
type PropertySeq <: Property
    primitive::PropertyPrimitive
    next::Property

    function PropertySeq(primitive::PropertyPrimitive,
                         next::Property)
        new(primitive, next)
    end

    function PropertySeq(primitive::PropertyPrimitive)
        new(primitive, empty_property)
    end

    # shallow copy constructor
    function PropertySeq(ps::PropertySeq)
        new(ps.primitive, ps.next)
    end
end

copy(ps::PropertySeq) = PropertySeq(ps)

function length(p::Property)
    n = 0
    while !is(p, empty_property)
        n += 1
        p = p.next
    end
    n
end

# Combination of properties, which is simply list concatenation.
function combine(head::Property, rest::Property...)
    a = b = head === empty_property ? head : copy(head)
    for p in rest
        if p === empty_property
            continue
        elseif a === empty_property
            a = b = copy(p)
        else
            while !is(b.next, empty_property)
                b.next = copy(b.next)
                b = b.next
            end
            b.next = p
        end
    end
    a
end

combine(head::Property, ::Nothing) = head

# Unit conversion functions.

function native_measure(property::EmptyProperty,
                        t::NativeTransform,
                        unit_box::BoundingBox,
                        box::NativeBoundingBox,
                        backend::Backend)
    property
end

# Convert a sequence of properties to a sequence of properties in native
# coordinates.
function native_measure(property::PropertySeq,
                        t::NativeTransform,
                        unit_box::BoundingBox,
                        box::NativeBoundingBox,
                        backend::Backend)
    p = property = copy(property)
    while !is(p, empty_property)
        p.primitive = native_measure(p.primitive, t, unit_box, box, backend)
        p.next = copy(p.next)
        p = p.next
    end
    property
end

# Catchall for properties that don't require unit conversion.
function native_measure(property::PropertyPrimitive,
                        t::NativeTransform,
                        unit_box::BoundingBox,
                        box::NativeBoundingBox,
                        backend::Backend)
    property
end

# A property primitive controlling the fill color (or lack of color) of a form.
type Fill <: PropertyPrimitive
    value::ColorOrNothing

    Fill(value::ColorOrNothing) = new(value)
    Fill(value::String) = new(color(value))
end

# Singleton sequence contructor.
fill(value) = PropertySeq(Fill(value))

# A property primitive controlling the strok color (or lack of color) of a form.
type Stroke <: PropertyPrimitive
    value::ColorOrNothing

    Stroke(value::ColorOrNothing) = new(value)
    Stroke(value::String) = new(color(value))
end

# Singleton sequence contructor.
stroke(value) = PropertySeq(Stroke(value))

# A property primitive controlling the widths of lines drawn in stroke
# operations.
type LineWidth <: PropertyPrimitive
    value::Measure

    LineWidth(value::MeasureOrNumber) = new(size_measure(value))
end

# Singleton sequence contructor.
linewidth(value::MeasureOrNumber) = PropertySeq(LineWidth(value))

# Unit conversion
function native_measure(property::LineWidth,
                        t::NativeTransform,
                        unit_box::BoundingBox,
                        box::NativeBoundingBox,
                        backend::Backend)
    LineWidth(native_measure(property.value, t, unit_box, box, backend))
end

# Set the visible attribute in SVG output, or just skip rendering in raster
# graphics.
type Visible <: PropertyPrimitive
    value::Bool
end

visible(value::Bool) = PropertySeq(Visible(value))

type Opacity <: PropertyPrimitive
    value::Float64
end

opacity(value::Number) = PropertySeq(Opacity(convert(Float64, value)))

# Clipping path
type Clip <: PropertyPrimitive
    points::Vector{Point}
end

function clip(points::XYTupleOrPoint...)
    PropertySeq(Clip([convert(Point, point) for point in points]))
end

function native_measure(property::Clip,
                        t::NativeTransform,
                        unit_box::BoundingBox,
                        box::NativeBoundingBox,
                        backend::Backend)
    Clip([native_measure(point, t, unit_box, box, backend)
          for point in property.points])
end

# A property primitive assigning an ID, in particular in SVG, to enable
# manipulation of portions of the graphic.
type SVGID <: PropertyPrimitive
    value::String
end

# Singleton sequence contructor.
svgid(value::String) = PropertySeq(SVGID(html_escape_string(value)))

# A property setting the class field in a form's group element.
type SVGClass <: PropertyPrimitive
    value::String
end

# Singleton sequence contructor.
svgclass(value::String) = PropertySeq(SVGClass(html_escape_string(value)))

# A SVG <a> tag.
type SVGLink <: PropertyPrimitive
    target::String
end

# Singleton sequence contructor.
svglink(target::String) = PropertySeq(SVGLink(target))

# The font property primitive.
type Font <: PropertyPrimitive
    family::String
end

# Singletone sequence constructor.
font(family::String) = PropertySeq(Font(family))

# The font size property.
type FontSize <: PropertyPrimitive
    value::Measure

    FontSize(value::MeasureOrNumber) = new(size_measure(value))
end

# Singletone sequence constructor.
fontsize(value::MeasureOrNumber) = PropertySeq(FontSize(value))

# Native unit conversion.
function native_measure(property::FontSize,
                        t::NativeTransform,
                        unit_box::BoundingBox,
                        box::NativeBoundingBox,
                        backend::Backend)
    FontSize(native_measure(property.value, t, unit_box, box, backend))
end

# Events:
# These correspond to DOM events, and associate some javascript blurb that gets
# executed when the event is triggered.

const events = (:OnActivate, :OnClick, :OnFocusIn, :OnFocusOut,
                :OnLoad, :OnMouseDown, :OnMouseMove, :OnMouseOut,
                :OnMouseOver, :OnMouseUp)

for event in events
    event_lc  = symbol(lowercase(string(event)))
    @eval begin
        type ($event) <: PropertyPrimitive
            value::String
        end

        function ($event_lc)(value::String)
            PropertySeq(($event)(value))
        end
    end
end

# SVG Mask objects

type SVGMask <: PropertyPrimitive
    id::String
end

svgmask(id::String) = PropertySeq(SVGMask(id))

type SVGDefineMask <: PropertyPrimitive
    id::String
end

svgdefmask(id::String) = PropertySeq(SVGDefineMask(id))

# A general purpose SVG attribute.

type SVGAttribute <: PropertyPrimitive
    attribute::String
    value::String
end

function svgattribute(attribute::String, value::String)
    PropertySeq(SVGAttribute(attribute, value))
end

# Embedding raw markup in SVG files.

# This is a property that when processed will be saved and then appended to the
# end of the SVG file, but ignored by default on other backends. This allows us
# to use some of SVGs more obsure features without exposing all of SVG.

type SVGEmbed <: PropertyPrimitive
    markup::String
end

svgembed(markup::String) = PropertySeq(SVGEmbed(markup))

# Embedding raw js code in D3 output.

type D3Embed <: PropertyPrimitive
    code::String
end

type D3Hook <: PropertyPrimitive
    code::String
end

d3hook(code::String) = PropertySeq(D3Hook(code))


text:
"""
# script bench_solver

- Julia version: 1.4
- Author: niallcullinane
- Date: 2020-08-20

Script to benchmark GPU and CPU implementions of model.
"""

using DrWatson
@quickactivate "MEngProject"
using MEngProject,
    # CUDA,
    DifferentialEquations,
    PyPlot,
    NNlib,
    ImageFiltering,
    Images,
    ImageIO,
    MEngProject.LaminartKernels,
    MEngProject.LaminartInitFunc,
    MEngProject.Utils,
    BenchmarkTools,
    Test

using OrdinaryDiffEq,
    ParameterizedFunctions, LSODA, Sundials, DiffEqDevTools, Noise

batch = 1000

global benchm_s = []
global benchm_sname = []
global alg_ = []
global stats_ = []
global benchm_sc = []
global benchm_snamec = []
global prob_s

tspan = (0.0f0, 800f0)

batch_ = string(batch, "_", rand(1000:9999))
mkdir(plotsdir(string("bench_solver", batch_)))
file = "kan_sq_cont_l.png"

solvers = [AutoTsit5(Rosenbrock23()), Tsit5(), BS3()lsoda(), Vern7()]
# alg=lsoda()

# GPU
p = LaminartInitFunc.parameterInit_conv_gpu(
    datadir("img", file),
    Parameters.parameters_f32,
);

u0 = cu(reshape(
    zeros(Float32, p.dim_i, p.dim_j * (5 * p.K + 2)),
    p.dim_i,
    p.dim_j,
    5 * p.K + 2,
    1,
))

arr1 = similar(@view u0[:, :, 1:2, :])
arr2 = similar(@view u0[:, :, 1:1, :])
f = LaminartFunc.LamFunction(
    arr1, #x
    similar(arr1), #m
    similar(arr1), #s
    arr2, #x_lgn,
    similar(arr1), #C,
    similar(arr1), #H_z,
    similar(arr1), # dy_temp,
    similar(arr1), # dm_temp,
    similar(arr1), # dz_temp,
    similar(arr1), # ds_temp,
    similar(arr2), # dv_temp,
    similar(arr1), # H_z_temp,
    similar(arr2), #  V_temp_1,
    similar(arr2), #  V_temp_2,
    similar(arr1), #  Q_temp,
    similar(arr1), #   P_temp
)
prob_i = ODEProblem(f, u0, tspan, p)

for alg in solvers
    @. u0 = 0.0f0
    try
        sol = solve(prob_s, alg)
        push!(alg_, sol.alg)
        push!(stats_, sol.destats)
        push!(benchm_s, @benchmark solve(prob_s, alg))
    catch err
        print(err)
    end
end

# CPU conv

# p = LaminartInitFunc.parameterInit_imfil_cpu(
#     datadir("img", file),
#     Parameters.parameters_f32,
# );
#
# u0 = reshape(
#     zeros(Float32, p.dim_i, p.dim_j * (5 * p.K + 2)),
#     p.dim_i,
#     p.dim_j,
#     5 * p.K + 2,
#     1,
# )
#
#
# f = LaminartFunc.LamFunction_imfil_cpu(
#     similar(arr1), # H_z_temp,
#     similar(arr2[:,:,1]), # v_C_temp1,
#     similar(arr2[:,:,1]), # v_C_temp2,
#     similar(arr1), # v_C_tempA,
#     similar(arr1[:,:,1]), #W_temp
# )
# prob_s = ODEProblem(f, u0, tspan, p)
# push!(benchm_i, @benchmark solve(prob_s))

# benchmark plot

# time
fig, ax = plt.subplots()
for ben in enumerate(test_name_plt)
    ax.scatter(
        ben[2],
        median(benchm_i[ben[1]].times) * 1e-9,
        color = Utils.colours[ben[1]],
        edgecolors = "none",
    )
end

ax.set_ylabel("Time (\$s\$)")
ax.set_ylim(ymin = 0)
ax.grid(true)
fig.tight_layout()
plt.savefig(plotsdir(
    string("bench_solver", batch_),
    string("bench_solver_time.png"),
))
close("all")

# memory

fig, ax = plt.subplots()
for ben in enumerate(test_name_plt)
    ax.scatter(
        ben[2],
        benchm_i[ben[1]].memory * 1e-6,
        color = Utils.colours[ben[1]],
        edgecolors = "none",
    )
end

ax.set_ylabel("Memory (\$MB\$)")
ax.set_ylim(ymin = 0)
ax.grid(true)
fig.tight_layout()
plt.savefig(plotsdir(
    string("bench_solver", batch_),
    string("bench_solver_mem.png"),
))
close("all")

# alloc

fig, ax = plt.subplots()
for ben in enumerate(test_name_plt)
    ax.scatter(
        ben[2],
        benchm_i[ben[1]].allocs,
        color = Utils.colours[ben[1]],
        edgecolors = "none",
    )
end

ax.set_ylabel("Allocations")
ax.set_ylim(ymin = 0)
ax.grid(true)
fig.tight_layout()
plt.savefig(plotsdir(
    string("bench_solver", batch_),
    string("bench_solver_alloc.png"),
))
close("all")

# benchm_i = nothing
# prob_i = nothing
# GC.gc
# CUDA.reclaim()