text:
using Distributions
using Statistics
using DelimitedFiles
using Random
#using Plots
using PyPlot

function load_data(path, file, shufffle = false)
    arti_rep1 = readdlm(path * file, ',')
    arti_seq = arti_rep1[2:end,1]
    arti_seq = convert(Array{String,1}, arti_seq)

    arti_exp = arti_rep1[2:end, end]
    for i in 1:length(arti_exp)
        if arti_exp[i] == "NA"
            arti_exp[i] = missing
        end
    end
    arti_exp = convert(Array{Union{Float64,Missing},1}, arti_exp)
    rng = MersenneTwister(1234);
    if shufffle == true
        shuffle(rng, arti_exp)
    end
    return arti_seq, arti_exp
end

function mutational_neightbors(seq1, base = ['A','U','C','G'])::Array{String,1}
    mut_neighbors::Array{String,1} = ["" for i in 1:((length(base)-1)*length(seq1))]
    index = 1
    for i in 1:length(seq1)
        if seq1[i] != 'N'
            for ch in base
                if ch != seq1[i]
                    mut_neighbors[index] = seq1[1:i-1]*ch*seq1[i+1:end]
                    index += 1
                    #push!(mut_neighbors, seq1[1:i-1]*ch*seq1[i+1:end])
                end
            end
        end
    end
    return mut_neighbors
end

function pearson_correlation(x,y)
    ind_x = [!ismissing(i) for i in x]
    ind_y = [!ismissing(i) for i in y]
    ind = ind_x .& ind_y
    x2::Array{Float64,1} = x[ind]
    y2::Array{Float64,1} = y[ind]
    return cor(x2,y2), length(x2),x2,y2
end

arti_seq, arti_array = load_data("/Users/imac/Desktop/Julia_repo/", "arti_SDR_union_count25.csv")
#arti_dict = Dict(zip(arti_seq,arti_array))
fepb_seq, fepb_array = load_data("/Users/imac/Desktop/Julia_repo/", "fepB_SDR_union_count25.csv")
#fepb_dict = Dict(zip(fepb_seq,fepb_array))
dmsc_seq, dmsc_array = load_data("/Users/imac/Desktop/Julia_repo/", "dmsC_SDR_union_count25.csv")

yfp_seq, yfp_array = load_data("/Users/imac/Desktop/Julia_repo/", "yfp_SDR_rep1_count25.csv")

function neighbor_center_correlation(seq_list::Array{String,1},exp_list::Array{Union{Missing, Float64},1},show::Bool = true)::Float64
    lib_dict = Dict{String,Union{Missing, Float64}}(zip(seq_list,exp_list))
    center::Array{Float64,1} = []
    neighbor_mean_value_array::Array{Float64,1} = []
    for seq in seq_list
        neighbors = mutational_neightbors(seq)
        if (!ismissing(lib_dict[seq])) & all([!isequal(lib_dict[ne],missing) for ne in neighbors])
            push!(center, lib_dict[seq])
            push!(neighbor_mean_value_array, mean([lib_dict[ne] for ne in neighbors]))
        end
    end
    # Calculate Pearson Correlation
    r = cor(center, neighbor_mean_value_array)
    if show == true
        println("Pearson correlation is ", r)
    end
    return r
end

#@time r_arti = neighbor_center_correlation(arti_seq,arti_array)

function MonteCarlo_permutation_test(lib_seq,lib_exp, sample_size = 1000, tail = 2 ,filename="")
    r_from_data::Float64 = neighbor_center_correlation(lib_seq,lib_exp)
    println("r from sample is $r_from_data")
    greater_count::Int64 = 0
    r_distribution::Array{Float64,1} = []
    rng = MersenneTwister(1234);
    for i in 1:sample_size
        lib_exp_new = shuffle(rng, lib_exp)
        r = neighbor_center_correlation(lib_seq,lib_exp_new, false)
        push!(r_distribution,r)
        #println(r)
        if tail == 1
            if r_from_data > 0
                if r > r_from_data
                    greater_count += 1
                end
            elseif r_from_data < 0
                if r < r_from_data
                    greater_count += 1
                end
            end

        elseif tail == 2
            if abs(r) > abs(r_from_data)
                greater_count += 1
            end
        end
    end
    println("Write!")
    writedlm("/Users/imac/Desktop/Julia_repo/" * filename * "_" * "r_distribution" * "_$sample_size" * ".csv", r_distribution, ",")
    println("P-value under $sample_size sampling = ", greater_count/sample_size)
    #histogram(r_distribution)
    #savefig("Users/imac/Desktop/Julia_repo/" * filename * "_" * "histogramOFr.svg")
end

function plot_r_distribution(filename = "")
    r_distribution = readdlm("/Users/imac/Desktop/Julia_repo/" * filename * "_" * "r_distribution_old.csv")
    #fitted_dist = fit_mle(Normal, r_distribution)
    #histogram(r_distribution, legend=false, xlims=(-0.05,0.05), fillcolor = :green)
    #plot!(fitted_dist, legend = false, xlims=(-1,1))

    h = plt.hist(r_distribution, 50, color = "green")
    xlim(-0.03,0.03)
    ylim(0,750)
    xticks(-0.03:0.01:0.03)
    yticks(0:250:750)
    #xlabel("r")
    #ylabel("Counts")
    savefig(filename * "_" * "histogramOFr.svg")
    close()
end

function cal_p_value_from_rlist(list::Array{Float64,1}, target::Float64, tail = 2)
    len = length(list)
    if tail == 1
        if target > 0
            count = sum([i > target for i in list])
        elseif target < 0
            count = sum([i < target for i in list])
        end
    elseif tail == 2
        count = sum([abs(i) > abs(target) for i in list])
    end
    p_value = count/len
    println(p_value)
    return p_value
end

# sum(skipmissing([1, missing]))

#plot_neighbor_vs_center(X,Y)

#println(">>>>>> arti >>>>>>")
#@time MonteCarlo_permutation_test(arti_seq,arti_array, 10000, "arti")
#println(">>>>>> dmsC >>>>>>")
#@time MonteCarlo_permutation_test(dmsc_seq,dmsc_array, 10000, "dmsC")
#println(">>>>>> fepB >>>>>>")
#@time MonteCarlo_permutation_test(fepb_seq,fepb_array, 10000, "fepB")
#println(">>>>>> yfp >>>>>>")
#@time MonteCarlo_permutation_test(yfp_seq,yfp_array, 10000, "yfp")

function load_r(filename = "")
    r_distribution = readdlm("/Users/imac/Desktop/Julia_repo/" * filename * "_" * "r_distribution_old.csv")
    return r_distribution[:,1]
end

arti_r = load_r("arti")
fepb_r = load_r("fepb")
dmsc_r = load_r("dmsc")
yfp_r = load_r("yfp")

println(">>> Original data")
cal_p_value_from_rlist(fepb_r, 0.929)
cal_p_value_from_rlist(arti_r, 0.921)
cal_p_value_from_rlist(dmsc_r, 0.883)
cal_p_value_from_rlist(yfp_r, 0.922)

println(">>> Shuffled data")
cal_p_value_from_rlist(fepb_r, -0.00751, 2)
cal_p_value_from_rlist(arti_r, -0.003329, 2)
cal_p_value_from_rlist(dmsc_r, -0.0032, 2)
#cal_p_value_from_rlist(yfp_r, -0.000445)

plot_r_distribution("arti")
plot_r_distribution("dmsC")
plot_r_distribution("fepB")
plot_r_distribution("yfp")


text:
if !isdir("../data") mkdir("../data") end
if !isdir("../data/horizons") mkdir("../data/horizons") end
println("build.jl running from ", pwd())

#include("build_julia_specific_conda.jl")
include("pyimports.jl")

ext: jl


text:
#--------------------------------------------------------------------------------------------------
using Test
using Mahler

#--------------------------------------------------------------------------------------------------

@test Mahler.sgn(2) == 1
@test Mahler.sgn(7) == -1

@test Mahler.fac(5)  == 120
@test Mahler.fac(20) == factorial(20)
@test Mahler.fac(20) == factorial(big(20))
@test Mahler.fac(21) == factorial(big(21))

@test Mahler.bin(4,3)    == 4
@test Mahler.bin(3,4)    == 0
@test Mahler.bin(33,4)   == binomial(33,4)
@test Mahler.bin(133,44) == binomial(big(133),big(44))
@test Mahler.bin(4.0,3)  == 4.0
@test Mahler.bin(3.0,4)  == 0.0
@test Mahler.bin(33.0,4)  ≈ Mahler.bin(33,4)
@test Mahler.bin(133.0,4) ≈ Mahler.bin(133,4)

#--------------------------------------------------------------------------------------------------

max_stars_repo_name: ayarodionov/Mahler

max_issues_repo_name: ayarodionov/Mahler

max_forks_repo_name: ayarodionov/Mahler


text:
using Documenter
using MultipleTesting

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

makedocs(
    modules = [MultipleTesting],
    sitename = "MultipleTesting.jl",
    format = Documenter.HTML(disable_git = true,  # disable source and edit links to github
        canonical = "https://juliangehring.github.io/MultipleTesting.jl/stable/"),
    doctest = true,
    checkdocs = :exports,
    linkcheck = true,
    pages = [
        "Home" => "index.md",
        "Library" => Any[
            "adjustment.md",
            "combination.md",
            "pi0.md",
            "higher-criticism.md",
            "models.md",
        ]
    ]
)

deploydocs(
    repo = "github.com/juliangehring/MultipleTesting.jl.git",
    target = "build",
    devbranch = "develop",
    push_preview = true
)


text:
using Plots

include("src/pwl_approx.jl")

## Define a submodular cardinality based function via Combined gadget
k = 2
w = random_scb_function(k::Int64)

w = 10*w/maximum(w)
pl = scatter(0:k,w,legend = false, xticks = 0:k)
epsi = 10.0
z0, zk, a, b, cgf = SubCardFun_to_CGF_weights(w,epsi,false)
Jtrue = length(a) + 1
J = Jtrue
xs = 0:0.01:k
ys = cgf(xs)
plot!(pl,xs,ys,color= :blue)

# Refine to get better approximation
z0l, zkl, al, bl, cgfl, best_eps = Refined_SCB_to_CGF(w,epsi)
xs = 0:0.01:k
ys = cgfl(xs)
plot!(pl,xs,ys,color= :red)

max_stars_repo_name: nveldt/SparseCardDSFM

max_issues_repo_name: nveldt/SparseCardDSFM

max_forks_repo_name: nveldt/SparseCardDSFM


text:
"""
    http_get(url::String)

Perform GET request on url and return the response.
"""
function http_get(url::String)
    return HTTP.get(url)
end

"""
    parse_body(response::HTTP.Messages.Response)

Parse the body of a HTTP response.
"""
function parse_body(response::HTTP.Messages.Response)
    return JSON.parse(String(response.body))
end

"""
    querry_api(url::String, action::String, parameters::Dict)

Querry an API with the given parameters, for a given CKAN action.
"""
function querry_api(url::String, action::String, parameters::Pair...)

    url *= "action/" * action
    url *= generate_request_query(parameters...)

    @info "Querrying API with url $(shorten_url(url))"

    response = parse_body(http_get(url))

    if !(response["success"])
        error("Querry failed")
    end

    return response
end

"""
    generate_request_query(parameters::Pair...)

Copied from Mangal.jl, formats a query string for use in a url, from a set of
Pairs.
"""
function generate_request_query(parameters::Pair...)
    query = ""
    for (i, pair) in enumerate(promote(parameters...))
        i == 1 || (query *= "&")
        i == 1 && (query *= "?")
        query *= pair.first * "=" * string(pair.second)
    end
    return replace(query, " " => "+")
end


text:
# Empirical estimation of CDF and PDF

## Empirical CDF

function ecdf{T<:Real}(X::RealVector{T})
    Xs = sort(X)
    n = length(X)

    ef(x::Real) = searchsortedlast(Xs, x) / n

    function ef(v::RealVector)
        ord = sortperm(v)
        m = length(v)
        r = Array(T, m)
        r0 = 0

        i = 1
        for x in Xs
            while i <= m && x > v[ord[i]]
                r[ord[i]] = r0
                i += 1
            end
            r0 += 1
            if i > m
            	break
            end
        end
        while i <= m
            r[ord[i]] = n
            i += 1
        end
        return r / n
    end

    return ef
end


text:
# Include the known function
using Random
using Distributions
using IntervalSets
using SparseArrays

using CautiousEngine

# Get the absolute path of the experiment dir
EXPERIMENT_DIR = @__DIR__
experiment_type = "verification-example-refinement"
global_exp_dir = "$EXPERIMENT_DIR/$experiment_type/global"
specification_file = "foo"

#======================================================================================================
0. Define the Unknown Function, Process Noise, and the Full Dynamics Function
======================================================================================================#
system_mat = [0.9 0.5; -0.5 0.4]
f = (x) -> system_mat*x

system = "mode1"
known_part = nothing

# Define the zero-mean process noise
σ_proc = 0.01
measurement_noise_dist = Normal(0., σ_proc)
process_noise_dist = nothing                                # Process noise is only used in synthesis, not verification

#======================================================================================================
1.  Perform GP regression and construct IMDP for each unknown mode.
======================================================================================================#
exp_dir = global_exp_dir
res_dirs = Array{String,1}()                        # Array that holds the single-mode result directories
X = Dict("x1" => [-2., 2.], "x2" => [-2., 2.])      # Defines the compact safe set
grid_sizes = Dict("x1" => 1.0, "x2" => 1.0)       # Specifies the discretization coarseness in each dimension
safety_dims = ["x1", "x2"]                          # Dimensions that have safety bounds
dependency_dims = Dict("x1" => [1, 1],              # Defines dependencies on other states
                       "x2" => [1, 1])
number_of_datapoints = 400                          # Number of data points per mode to generate
run_exps_flag = true                                # Set to true to rerun the single-mode constructions
m_opt = -1                                          # Unused
bound_type = "rkhs-tight"                           # GP bound type
epsilon = -1.                                       # Distance bound parameter, set to -1 for automatic determination

Lf_bound = 1.5                                 # Lipschitz constant bound for the mode over X
random_seed = 11                             	# Random seed for the random twister
eta = -1.                                       # η parameter set in a naiive manner

# Setup experiment parameter structures (defined in CautiousEngine.jl)
system_params = SystemParameters(system, f, known_part, measurement_noise_dist, process_noise_dist, Lf_bound, dependency_dims, 2, 2)
data_params = DataParameters(number_of_datapoints, m_opt, bound_type, σ_proc, epsilon, eta, safety_dims, -1., -1.)
experiment_params = ExperimentParameters(exp_dir, experiment_type, specification_file, X, grid_sizes, "foo", random_seed, system_params, data_params)

_, trans_mats = CautiousEngine.end_to_end_transition_bounds(experiment_params, reuse_gps_flag=true, reuse_regions_flag=true)
results_path = experiment_params.experiment_directory

#======================================================================================================
2.  Construct the IMDP Labels and call the appropriate verification solver wrapper.
======================================================================================================#
# Wrapper function
imdp = CautiousEngine.create_simple_imdp(trans_mats["minPr"], trans_mats["maxPr"])
refinement_steps = 3

k_vals = [1,]
for k in k_vals
    global imdp
    verification_result_mat = CautiousEngine.Globally(imdp, "safe", k, "$results_path/imdp.txt")
    CautiousEngine.save_legacy_mats(verification_result_mat, results_path, k)
    CautiousEngine.plot_2d_verification_results(results_path, prob_plots=true, state_outlines_flag=true)
    working_dir = results_path

    # TODO: Reusing regions only works for a single instance of k. Handle this in a smarter way.
    refine_result_dirs = CautiousEngine.safety_based_refinement(experiment_params, results_path, verification_result_mat, refinement_steps, horizon=k, reuse_regions_flag=false, reuse_transition_mats_flag=false)

    for res in refine_result_dirs
        CautiousEngine.plot_2d_verification_results(res, prob_plots=true, state_outlines_flag=true)
    end
end


text:
module Postgres

  using Logging, GnuTLS, Compat

  include("message.jl")
  include("io.jl")
  include("error.jl")
  include("utils.jl")
  include("types.jl")

export pg_connect
end # module


text:
module DoubleUAVs_module

push!(LOAD_PATH, "./")
disc = "mid3"
include("./discretizations/" * disc * "_disc.jl")
include("./discretizations/common.jl")

using Distributions

export states, advisories, weights, advisory_to_action_dict, disc, state_dim
export sigma_str, sigma_v, sigma_alert_deg, sigma_coc_deg
export terminal_state, terminal_state_var
export num_states, num_actions, num_sigmas
export get_reward, get_next_state
export is_terminal
export norm_angle
export rho_max, sensing_range, v_min, v_max, collision_range
export rho_ind, theta_ind, phi_ind, v_own_ind, v_int_ind
export grid_states
export actual_weight

"""
	tells if the state is terminal
"""
function is_terminal(state::Vector{Float64})
	rho = state[rho_ind]
	if rho > rho_max
		return true
	else
		return false
	end
end

"""
	keeps the angles in [0, 2*pi]
"""
function norm_angle(angle::Float64)
	return ((angle % (2 * pi)) + 2 * pi) % (2 * pi)
end

"""
	state update
	assuming the intruder keeps the straight path
"""
function get_next_state(state::Vector{Float64}, advisory::Symbol; isigma::Int64=1, if_sample::Bool=false)
	if is_terminal(state)
		return terminal_state
	else
		if if_sample == false
			rho = state[rho_ind]
			theta = state[theta_ind] + sigmas[theta_sigma_ind, isigma]
			phi = state[phi_ind] + sigmas[phi_sigma_ind, isigma]
			v_own = state[v_own_ind]
			v_int = state[v_int_ind]

			if advisory == :COC
				turn_rate = sign(sigmas[turn_rate_sigma_ind, isigma]) * sigma_coc
			else
				turn_rate = deg2rad(advisory_to_action_dict[advisory]) + sigmas[turn_rate_sigma_ind, isigma]
			end
			turn_angle = turn_rate * dt

			v_own += sigmas[v_own_sigma_ind, isigma]
			v_int += sigmas[v_int_sigma_ind, isigma]

		else # if_sample == true
			# rho = state[rho_ind]
			# theta = state[theta_ind]
			# phi = state[phi_ind] + clamp(rand(Normal(0, sigma_phi)), -sigma_phi, sigma_phi)
			# v_own = state[v_own_ind]
			# v_int = state[v_int_ind]

			# if advisory == :COC
			# 	turn_rate = clamp(rand(Normal(0, sigma_coc)), -sigma_coc, sigma_coc)
			# else
			# 	turn_rate = deg2rad(advisory_to_action_dict[advisory]) +
			# 		clamp(rand(Normal(0, sigma_alert)), -sigma_alert, sigma_alert)
			# end
			# turn_angle = turn_rate * dt

			# v_own += clamp(rand(Normal(0, sigma_v)), -sigma_v, sigma_v)
			# v_int += clamp(rand(Normal(0, sigma_v)), -sigma_v, sigma_v)
			nothing
		end
		###################################################################
		# Using the body fixed frame origined at the ownship.
		# The positive x-direction is coincide with the heading of the ownship.
		# desperate

		x_own = 0.
		y_own = 0.
		x_int = rho * cos(theta)
		y_int = rho * sin(theta)

		next_x_own = x_own + v_own * cos(turn_angle) * dt
		next_y_own = y_own + v_own * sin(turn_angle) * dt
		next_x_int = x_int + v_int * cos(phi) * dt
		next_y_int = y_int + v_int * sin(phi) * dt

		next_rho = sqrt((next_x_own - next_x_int)^2 + (next_y_own - next_y_int)^2)

		if next_rho > rho_max
			return terminal_state
		end

		temp_var = norm_angle(atan2(next_y_int - next_y_own, next_x_int - next_x_own))
		next_theta = norm_angle(temp_var - turn_angle)
		next_phi = norm_angle(phi - turn_angle)

		next_state = Vector{Float64}(state_dim)
		next_state[rho_ind] = next_rho
		next_state[theta_ind] = next_theta
		next_state[phi_ind] = next_phi
		next_state[v_own_ind] = v_own
		next_state[v_int_ind] = v_int

		return next_state
	end
end

"""
	reward function
"""
function get_reward(state::Vector{Float64}, advisory::Symbol;
	penalty_action=0.02,
	penalty_closeness=10.0,
	penalty_nmac=1000.0,
	penalty_conflict=1.0)

	possible_max_rew = 0.0
	possible_min_rew = - penalty_conflict - penalty_action * 10.^2 - penalty_nmac - penalty_closeness * e

	action = advisory_to_action_dict[advisory] # [deg]
	reward = 0.0
	if advisory == :COC
		action = 0.0
	else
		reward -= penalty_conflict
	end

	if is_terminal(state)
		reward -= penalty_action * action^2
		reward /= (possible_max_rew - possible_min_rew)
		if abs(reward) > 1
			println(reward)
		end
		return reward
	else
		reward -= penalty_action * action^2
		rho = state[rho_ind]
		if rho < collision_range
			reward -= penalty_nmac
		end
		reward -= penalty_closeness * exp(-(rho - collision_range) / collision_range)
	end

	# scale the magnitude of every step reward between 0 and 1
	reward = reward / (possible_max_rew - possible_min_rew)
	if abs(reward) > 1
		println(reward)
	end
	return reward
end

end # module

max_stars_repo_name: parachutel/DeepCorrectionCAS


text:
# Model
include_model("cybertruck")

# ramp
jump_slope = 0.3
jump_length = 1.0

function ϕ_func(::CYBERTRUCK, q)
    if q[1] < jump_length && q[1] > 0.0
        return @SVector [q[3] - jump_slope * q[1]]
    else
        return @SVector [q[3]]
    end
end

# Horizon
T = 21

# Time step
tf = 2.0
h = tf / (T - 1)

# Bounds
_uu = Inf * ones(model.m)
_uu[model.idx_u] = [Inf; 1.0]
_ul = zeros(model.m)
_ul[model.idx_u] .= [-1.0; -1.0]

ul, uu = control_bounds(model, T, _ul, _uu)

# Initial and final states
q1 = [-1.0, 0.0, 0.0, 0.0]
qT = [3.0, 0.0, 0.0, 0.0]

x1 = [q1; q1]
xT = [qT; qT]

xl, xu = state_bounds(model, T,
    x1 = x1, xT = xT)

# Objective
Qq = Diagonal(ones(model.nq))
Q = cat(0.5 * Qq, 0.5 * Qq, dims = (1, 2))
QT = cat(10.0 * Qq, 1000.0 * Diagonal(ones(model.nq)), dims = (1, 2))
R = Diagonal([1.0e-1 * ones(model.nu)..., zeros(model.m - model.nu)...])

obj_tracking = quadratic_tracking_objective(
    [t < T ? Q : QT for t = 1:T],
    [R for t = 1:T-1],
    [xT for t = 1:T],
    [zeros(model.m) for t = 1:T])

obj_penalty = PenaltyObjective(1.0e5, model.m)
obj = MultiObjective([obj_tracking, obj_penalty])

# Constraints
include_constraints(["contact", "control_complementarity"])

n_cc_stage = 2
n_cc_con = n_cc_stage * (T - 1)
con_ctrl_comp = ControlComplementarity(n_cc_con, (1:n_cc_con), n_cc_stage)

con_contact = contact_constraints(model, T)

con = multiple_constraints([con_contact, con_ctrl_comp])

# Problem
prob = trajectory_optimization_problem(model,
               obj,
               T,
               h = h,
               xl = xl,
               xu = xu,
               ul = ul,
               uu = uu,
               con = con)

# Trajectory initialization
q_ref = linear_interpolation(q1, qT, T)
x_ref = configuration_to_state(q_ref)
x0 = deepcopy(x_ref) # linear interpolation on state
u0 = [1.0e-5 * rand(model.m) for t = 1:T-1] # random controls

# Pack trajectories into vector
z0 = pack(x0, u0, prob)

#NOTE: may need to run examples multiple times to get good trajectories
# Solve nominal problem
include_snopt()
@time z̄ , info = solve(prob, copy(z0),
    nlp = :SNOPT7,
    tol = 1.0e-3, c_tol = 1.0e-3)

check_slack(z̄, prob)
x̄, ū = unpack(z̄, prob)

include(joinpath(pwd(), "models/visualize.jl"))
vis = Visualizer()
open(vis)
default_background!(vis)

obj_path = joinpath(pwd(), "models/cybertruck/cybertruck.obj")
mtl_path = joinpath(pwd(), "models/cybertruck/cybertruck.mtl")

ctm = ModifiedMeshFileObject(obj_path, mtl_path, scale=0.1)
setobject!(vis["cybertruck"], ctm)
settransform!(vis["cybertruck"], LinearMap(RotZ(pi) * RotX(pi / 2.0)))

Ns = 100
height = range(0, stop = jump_slope * jump_length, length = Ns)
width = range(0, stop = jump_length, length = Ns)
wid = jump_length / Ns
for i = 1:Ns
    setobject!(vis["stair$i"],
        Rect(Vec(0., 0.0, 0.0), Vec(0.01, 1.0, height[i])),
            MeshPhongMaterial(color = RGBA(0.5, 0.5, 0.5, 1.0)))
    settransform!(vis["stair$i"], Translation(width[i], -0.5, 0))
end

q = state_to_configuration([[x̄[1] for i = 1:10]...,x̄..., [x̄[end] for i = 1:10]...])
anim = MeshCat.Animation(convert(Int,floor(1.0 / h)))
for t = 1:length(q)
    MeshCat.atframe(anim, t) do
        settransform!(vis["cybertruck"],
            compose(Translation((q[t][1:3] + ((q[t][1] < jump_length + 1.0 && q[t][1] > -0.1) ? [0.0; 0.0; 0.01] : zeros(3))) ...),
            LinearMap(RotZ(q[t][4] + pi) * RotY((q[t][1] < jump_length + 0.0 && q[t][1] > -0.25) ? tan(jump_slope) : 0.0) * RotX(pi / 2.0))))
    end
end
MeshCat.setanimation!(vis, anim)
settransform!(vis["/Cameras/default"],
    compose(LinearMap(RotZ(-pi / 2.0)), Translation(0.0, 3.0, 1.5)))


text:
include("../src/Transit.jl")
# To run this, Pkg.clone the Benchmarks.jl package at:
#
# https://github.com/johnmyleswhite/Benchmarks.jl/
#
# with:
#
# > Pkg.clone("https://github.com/johnmyleswhite/Benchmarks.jl")

import Benchmarks
import Transit

const SEATTLEDIR = "../transit-format/examples/0.8/"
const JSONFILE = string(SEATTLEDIR, "example.json")
const JSONVERBOSEFILE = string(SEATTLEDIR, "example.verbose.json")

function run_tests(data::AbstractString)
    io = IOBuffer(data)
    Transit.parse(io)
end

function read_json()
    readall(open(JSONFILE))
end

function read_json_verbose()
    readall(open(JSONVERBOSEFILE))
end

println("Transit JSON parse")
println(@Benchmarks.benchmark run_tests(read_json()))
println("\n")
println("Transit JSON verbose parse")
println(@Benchmarks.benchmark run_tests(read_json_verbose()))


text:
function generate_problem_dictionary(path_to_parameters_file::String)::Dict{String,Any}

    # initialize -
    problem_dictionary = Dict{String,Any}()

    try

        # load the TOML parameters file -
        toml_dictionary = TOML.parsefile(path_to_parameters_file)["biophysical_constants"]

        # setup the initial condition array -
        initial_condition_array = [
            0.0 ;   # 1 mRNA
            5.0 ;   # gene concentration in nM (5 nM = 0.005 micromolar)
            0.0 ;   # 3 I = we'll fill this in the execute script #QN: WHAT DOES THIS MEAN?
        ]
        #trying to fix error message: "mRNA_half_life_in_min undefined"
        mRNA_half_life_in_min = 13.5                           # units: min
        #mRNA_half_life_in_sec = mRNA_half_life_in_min*60       #units: sec
        mRNA_half_life_in_sec = 810                            #units: sec
        transcription_elongation_rate=12                       # units: nt/s
        gene_length_in_nt=690                                   # units: nt
        RNAPII_concentration=60                                 # units: nM
        #RNAPII_concentration=0.6                                 # units: nM #test value

        # TODO: calculate the mRNA_degradation_constant in SECONDS
        mRNA_degradation_constant = log(2)/(mRNA_half_life_in_sec)

        # calculate the elongation constant (1/s)
        elongation_constant = transcription_elongation_rate/gene_length_in_nt

        # TODO: VMAX for transcription (nM/s)
        VMAX = elongation_constant*RNAPII_concentration

        # --- PUT STUFF INTO problem_dictionary ----
        problem_dictionary["transcription_time_constant"] = toml_dictionary["transcription_time_constant"]
        problem_dictionary["transcription_saturation_constant"] = toml_dictionary["transcription_saturation_constant"]
        problem_dictionary["E1"] = toml_dictionary["energy_promoter_state_1"]
        problem_dictionary["E2"] = toml_dictionary["energy_promoter_state_2"]
        problem_dictionary["inducer_dissociation_constant"] = toml_dictionary["inducer_dissociation_constant"]
        problem_dictionary["ideal_gas_constant_R"] = 0.008314 # kJ/mol-K
        problem_dictionary["temperature_K"] = (273.15+37)
        problem_dictionary["initial_condition_array"] = initial_condition_array
#trying to fix an error message
        problem_dictionary["inducer_cooperativity_parameter"]=toml_dictionary["inducer_cooperativity_parameter"]

        # TODO: If we want to use a value later e.g., VMAX or mRNA_degradation_constant you need to put them in the Dictionary
        problem_dictionary["mRNA_degradation_constant"]=mRNA_degradation_constant
        problem_dictionary["elongation_constant"]=elongation_constant
        problem_dictionary["VMAX"]=VMAX

        # return -
        return problem_dictionary
    catch error
        throw(error)
    end
end


text:
module Bar1
    include("foo.jl")
    using .Foo: @foo
    m() = r
    q() = z
    @foo y
end

include("foo.jl")
module Bar2
    using ..Foo: @foo, @foo_y, @multifoo
    @foo u
    @foo v
    @foo_y w
    @multifoo q
end


text:
module RocketMergedObservableTest

using Test
using Rocket

include("../test_helpers.jl")

@testset "MergeObservable" begin

    println("Testing: merged")

    @testset begin
        @test_throws Exception merged([ of(1), of(2.0) ])
    end

    @testset begin
        source = merged((of(1), of(2.0)))

        io = IOBuffer()

        show(io, source)

        printed = String(take!(io))

        @test occursin("MergeObservable", printed)
        @test occursin(string(eltype(source)), printed)

        subscription = subscribe!(source, void())

        show(io, subscription)

        printed = String(take!(io))

        @test occursin("MergeSubscription", printed)

        unsubscribe!(subscription)
    end

    run_testset([
        (
            source = merged((of(1), of(2.0))),
            values = @ts([ 1, 2.0, c ]),
            source_type = Union{Int, Float64}
        ),
        (
            source = merged((of(2.0), from("Hello"), from("World") |> async(0))) |> take(10),
            values = @ts([ 2.0, 'H', 'e', 'l', 'l', 'o' ] ~ [ 'W' ] ~ [ 'o' ] ~ [ 'r' ] ~ [ 'l', c ]),
            source_type = Union{Float64, Char}
        )
    ])

end

end


text:
"""
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
-                                QUANTUM STATES                                -
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
"""

qubits(N::Int) = siteinds("Qubit", N)

qudits(N::Int; dim::Int = 3) = siteinds("Qudit",N; dim = dim)

#
# State-like gates, used to define product input states
#

# TODO: add an arbitrary state specified by angles

state(sn::String) = state(StateName(sn))
state(sn::String, dim::Int) = state(StateName(sn), dim)

function state(sn::String, i::Index)
  st = hastags(i, "Qubit") ? state(sn) : state(sn, dim(i))
  return ITensors.itensor(st, i)
end

# Pauli eingenstates
state(::StateName"X+") = [
  1 / sqrt(2)
  1 / sqrt(2)
]

state(::StateName"X-") = [
  1 / sqrt(2)
  -1 / sqrt(2)
]

state(::StateName"Y+") = [
  1 / sqrt(2)
  im / sqrt(2)
]

state(::StateName"Y-") = [
  1 / sqrt(2)
  -im / sqrt(2)
]

state(::StateName"Z+") = [
  1
  0
]

state(::StateName"Z-") = [
  0
  1
]

state(::StateName"0") = state("Z+")
state(::StateName"1") = state("Z-")

# SIC-POVMs

state(::StateName"T1") = state("Z+")
state(::StateName"T2") = [
  1/√3
  √2/√3
]
state(::StateName"T3") = [
  1/√3
  √2/√3 * exp(im*2π/3)
]
state(::StateName"T4") = [
  1/√3
  √2/√3 * exp(im*4π/3)
]

function state(::StateName{N}, dim::Int) where {N}
  n = parse(Int, String(N))
  st = zeros(Int64,dim)
  st[n + 1] = 1
  return st
end

"""
    productstate(N::Int)

    productstate(sites::Vector{<:Index})

Initialize qubits to an MPS wavefunction in the 0 state (`|ψ⟩ = |0⟩ ⊗ |0⟩ ⊗ …`).
"""
# TODO: add dimension to use qudit instead of qubit
function productstate(N::Int; dim::Int = 2, sitetype::String = "Qubit")
  dim > 2 && return productstate(siteinds("Qudit", N; dim = dim))
  return productstate(siteinds(sitetype,N))
end

productstate(sites::Vector{<:Index}) = productMPS(sites, "0")

"""
    productstate(M::Union{MPS,MPO,LPDO})

Initialize qubits on the Hilbert space of a reference state,
given as `MPS`, `MPO` or `LPDO`.
"""
productstate(M::Union{MPS,MPO,LPDO}) = productstate(hilbertspace(M))

"""
    productstate(N::Int, states::Vector{T})

    productstate(sites::Vector{<:Index}, states::Vector{T})

    productstate(M::Union{MPS,MPO,LPDO}, states::Vector{T})

Initialize the qubits to a given product state, where the state `T` can be specified either
with a Vector of states specified as Strings or bit values (0 and 1).
"""
function productstate(N::Int, states::Vector; dim::Int = 2, sitetype::String = "Qubit")
  dim > 2 && return productstate(siteinds("Qudit", N; dim = dim), states)
  return productstate(siteinds(sitetype, N), states)
end

function productstate(sites::Vector{<:Index}, states::Vector{<:Integer})
  return MPS(state.(string.(Int.(states)), sites))
end

function productstate(sites::Vector{<:Index}, states::Vector)
  return MPS(state.(states, sites))
end

function productstate(M::Union{MPS,MPO,LPDO}, states::Vector)
  return productstate(hilbertspace(M), states)
end

function productstate(sites::Vector{<:Index}, state::Union{String,Integer})
  return productstate(sites, fill(state, length(sites)))
end

function productstate(sites::Vector{<:Index}, states::Function)
  return productstate(sites, map(states, 1:length(sites)))
end

"""
    productoperator(N::Int)

    productoperator(sites::Vector{<:Index})

Initialize an MPO that is a product of identity operators.
"""
function productoperator(N::Int; dim::Int = 2, sitetype::String = "Qubit")
  dim > 2 && return productoperator(siteinds("Qudit", N; dim = dim))
  return productoperator(siteinds(sitetype, N))
end

productoperator(M::Union{MPS,MPO,LPDO}) = productoperator(hilbertspace(M))

productoperator(sites::Vector{<:Index}) = MPO([op("Id", s) for s in sites])


text:
#####
##### AtmosModel
#####

module DycomsModel
##### PlanetParameters
# Physical constants
gas_constant =     8.3144598                      #          "Universal gas constant (J/mol/K)"
light_speed =      2.99792458e8                   #          "Speed of light in vacuum (m/s)"
h_Planck =         6.626e-34                      #          "Planck constant (m^2 kg/s)"
k_Boltzmann =      1.381e-23                      #          "Boltzmann constant (m^2 kg/s^2/K)"
Stefan =           5.670e-8                       #          "Stefan-Boltzmann constant (W/m^2/K^4)"
astro_unit =       1.4959787e11                   #          "Astronomical unit (m)"
k_Karman =         0.4                            #          "Von Karman constant (1)"

# Properties of dry air
molmass_dryair =   28.97e-3                       #          "Molecular weight dry air (kg/mol)"
R_d =              gas_constant/molmass_dryair    #          "Gas constant dry air (J/kg/K)"
kappa_d =          2//7                           #          "Adiabatic exponent dry air"
cp_d =             R_d/kappa_d                    #          "Isobaric specific heat dry air"
cv_d =             cp_d - R_d                     #          "Isochoric specific heat dry air"

# Properties of water
ρ_cloud_liq =      1e3                            #          "Density of liquid water (kg/m^3)"
ρ_cloud_ice =      916.7                          #          "Density of ice water (kg/m^3)"
molmass_water =    18.01528e-3                    #          "Molecular weight (kg/mol)"
molmass_ratio =    molmass_dryair/molmass_water   #          "Molar mass ratio dry air/water"
R_v =              gas_constant/molmass_water     #          "Gas constant water vapor (J/kg/K)"
cp_v =             1859                           #          "Isobaric specific heat vapor (J/kg/K)"
cp_l =             4181                           #          "Isobaric specific heat liquid (J/kg/K)"
cp_i =             2100                           #          "Isobaric specific heat ice (J/kg/K)"
cv_v =             cp_v - R_v                     #          "Isochoric specific heat vapor (J/kg/K)"
cv_l =             cp_l                           #          "Isochoric specific heat liquid (J/kg/K)"
cv_i =             cp_i                           #          "Isochoric specific heat ice (J/kg/K)"
T_freeze =         273.15                         #          "Freezing point temperature (K)"
T_min =            150.0                          #          "Minimum temperature guess in saturation adjustment (K)"
T_max =            1000.0                         #          "Maximum temperature guess in saturation adjustment (K)"
T_icenuc =         233.00                         #          "Homogeneous nucleation temperature (K)"
T_triple =         273.16                         #          "Triple point temperature (K)"
T_0 =              T_triple                       #          "Reference temperature (K)"
LH_v0 =            2.5008e6                       #          "Latent heat vaporization at T_0 (J/kg)"
LH_s0 =            2.8344e6                       #          "Latent heat sublimation at T_0 (J/kg)"
LH_f0 =            LH_s0 - LH_v0                  #          "Latent heat of fusion at T_0 (J/kg)"
e_int_v0 =         LH_v0 - R_v*T_0                #          "Specific internal energy of vapor at T_0 (J/kg)"
e_int_i0 =         LH_f0                          #          "Specific internal energy of ice at T_0 (J/kg)"
press_triple =     611.657                        #          "Triple point vapor pressure (Pa)"

# Properties of sea water
ρ_ocean =          1.035e3                        #          "Reference density sea water (kg/m^3)"
cp_ocean =         3989.25                        #          "Specific heat sea water (J/kg/K)"

# Planetary parameters
planet_radius =    6.371e6                        #          "Mean planetary radius (m)"
day =              86400                          #          "Length of day (s)"
Omega =            7.2921159e-5                   #          "Ang. velocity planetary rotation (1/s)"
grav =             9.81                           #          "Gravitational acceleration (m/s^2)"
year_anom =        365.26*day                     #          "Length of anomalistic year (s)"
orbit_semimaj =    1*astro_unit                   #          "Length of semimajor orbital axis (m)"
TSI =              1362                           #          "Total solar irradiance (W/m^2)"
MSLP =             1.01325e5                      #          "Mean sea level pressure (Pa)"

#####

using StaticArrays
abstract type AtmosConfiguration end
struct AtmosLESConfiguration <: AtmosConfiguration end

abstract type TurbulenceClosure end

abstract type ReferenceState end
struct HydrostaticState{P,F} <: ReferenceState
  temperatureprofile::P
  relativehumidity::F
end

abstract type PrecipitationModel end
struct NoPrecipitation <: PrecipitationModel end

abstract type TemperatureProfile end

struct LinearTemperatureProfile{FT} <: TemperatureProfile
  "minimum temperature (K)"
  T_min::FT
  "surface temperature (K)"
  T_surface::FT
  "lapse rate (K/m)"
  Γ::FT
end

struct SmagorinskyLilly{FT} <: TurbulenceClosure
  "Smagorinsky Coefficient [dimensionless]"
  C_smag::FT
end

abstract type MoistureModel end
Base.@kwdef struct EquilMoist <: MoistureModel
  maxiter::Int = 3
end

abstract type RadiationModel end
struct DYCOMSRadiation{FT} <: RadiationModel
  "mass absorption coefficient `[m^2/kg]`"
  κ::FT
  "Troposphere cooling parameter `[m^(-4/3)]`"
  α_z::FT
  "Inversion height `[m]`"
  z_i::FT
  "Density"
  ρ_i::FT
  "Large scale divergence `[s^(-1)]`"
  D_subsidence::FT
  "Radiative flux parameter `[W/m^2]`"
  F_0::FT
  "Radiative flux parameter `[W/m^2]`"
  F_1::FT
end

abstract type Orientation end
struct SphericalOrientation <: Orientation end
abstract type Source end
struct Gravity <: Source end

struct GeostrophicForcing{FT} <: Source
  f_coriolis::FT
  u_geostrophic::FT
  v_geostrophic::FT
end
struct RayleighSponge{FT} <: Source
  "Maximum domain altitude (m)"
  z_max::FT
  "Altitude at with sponge starts (m)"
  z_sponge::FT
  "Sponge Strength 0 ⩽ α_max ⩽ 1"
  α_max::FT
  "Relaxation velocity components"
  u_relaxation::SVector{3,FT}
  "Sponge exponent"
  γ::FT
end
struct Subsidence{FT} <: Source
  D::FT
end
struct GeostrophicForcing{FT} <: Source
  f_coriolis::FT
  u_geostrophic::FT
  v_geostrophic::FT
end

abstract type BoundaryCondition
end
struct DYCOMS_BC{FT} <: BoundaryCondition
  "Drag coefficient"
  C_drag::FT
  "Latent Heat Flux"
  LHF::FT
  "Sensible Heat Flux"
  SHF::FT
end

function init_dycoms!(bl, state, aux, (x,y,z), t)
    FT = eltype(state)
    return nothing
end

abstract type BalanceLaw end # PDE part
struct AtmosModel{FT,O,RS,T,M,P,R,S,BC,IS} <: BalanceLaw
  orientation::O
  ref_state::RS
  turbulence::T
  moisture::M
  precipitation::P
  radiation::R
  source::S
  # TODO: Probably want to have different bc for state and diffusion...
  boundarycondition::BC
  init_state::IS
end

struct FlatOrientation <: Orientation
  # for Coriolis we could add latitude?
end

function AtmosModel{FT}(::Type{AtmosLESConfiguration};
                         orientation::O=FlatOrientation(),
                         ref_state::RS=HydrostaticState(LinearTemperatureProfile(FT(200),
                                                                                FT(280),
                                                                                FT(grav) / FT(cp_d)),
                                                                                FT(0)),
                         turbulence::T=SmagorinskyLilly{FT}(0.21),
                         moisture::M=EquilMoist(),
                         precipitation::P=NoPrecipitation(),
                         radiation::R=NoRadiation(),
                         source::S=( Gravity(),
                                     Coriolis(),
                                     GeostrophicForcing{FT}(7.62e-5, 0, 0)),
                         # TODO: Probably want to have different bc for state and diffusion...
                         boundarycondition::BC=NoFluxBC(),
                         init_state::IS=nothing) where {FT<:AbstractFloat,O,RS,T,M,P,R,S,BC,IS}
  @assert init_state ≠ nothing

  atmos = (
        orientation,
        ref_state,
        turbulence,
        moisture,
        precipitation,
        radiation,
        source,
        boundarycondition,
        init_state,
       )

  return AtmosModel{FT,typeof.(atmos)...}(atmos...)
end

function get_model()
    FT = Float32
    xmax = 1000
    ymax = 1000
    zmax = 2500
    # Reference state
    T_min   = FT(289)
    T_s     = FT(290.4)
    Γ_lapse = FT(grav/cp_d)
    T       = LinearTemperatureProfile(T_min, T_s, Γ_lapse)
    rel_hum = FT(0)
    ref_state = HydrostaticState(T, rel_hum)

    # Radiation model
    κ             = FT(85)
    α_z           = FT(1)
    z_i           = FT(840)
    ρ_i           = FT(1.13)
    D_subsidence  = FT(0) # 0 for stable testing, 3.75e-6 in practice
    F_0           = FT(70)
    F_1           = FT(22)
    radiation = DYCOMSRadiation{FT}(κ, α_z, z_i, ρ_i, D_subsidence, F_0, F_1)

    # Sources
    f_coriolis    = FT(1.03e-4)
    u_geostrophic = FT(7.0)
    v_geostrophic = FT(-5.5)
    w_ref         = FT(0)
    u_relaxation  = SVector(u_geostrophic, v_geostrophic, w_ref)
    # Sponge
    c_sponge = 1
    # Rayleigh damping
    zsponge = FT(1500.0)
    rayleigh_sponge = RayleighSponge{FT}(zmax, zsponge, c_sponge, u_relaxation, 2)
    # Geostrophic forcing
    geostrophic_forcing = GeostrophicForcing{FT}(f_coriolis, u_geostrophic, v_geostrophic)

    # Boundary conditions
    # SGS Filter constants
    C_smag = FT(0.21) # 0.21 for stable testing, 0.18 in practice
    C_drag = FT(0.0011)
    LHF    = FT(115)
    SHF    = FT(15)
    bc = DYCOMS_BC{FT}(C_drag, LHF, SHF)
    ics = init_dycoms!
    source = (Gravity(),
              rayleigh_sponge,
              Subsidence{FT}(D_subsidence),
              geostrophic_forcing)

    model = AtmosModel{FT}(AtmosLESConfiguration;
                           ref_state=ref_state,
                          turbulence=SmagorinskyLilly{FT}(C_smag),
                            moisture=EquilMoist(5),
                           radiation=radiation,
                              source=source,
                   boundarycondition=bc,
                          init_state=ics)
    return model
end

model = get_model()

end


text:
# Portions translated from SLICOT-Reference distribution
# Copyright (c) 2002-2020 NICONET e.V.
function run_mb03pd(datfile, io=stdout)
    ZERO = 0.0e0
    ONE = 1.0e0
    NIN = 5
    NOUT = 6
    NMAX = 10
    MMAX = 10
    LDA = NMAX
    LDTAU = min(MMAX,NMAX)
    LDWORK = 3*MMAX
    A = Array{Float64,2}(undef, LDA,NMAX)
    JPVT = Array{BlasInt,1}(undef, MMAX)
    TAU = Array{Float64,1}(undef, LDTAU)
    SVAL = Array{Float64,1}(undef, 3)
    DWORK = Array{Float64,1}(undef, LDWORK)
    f = open(datfile,"r")
    readline(f)
    vs = split(readline(f))
    M = parse(BlasInt, vs[1])
    N = parse(BlasInt, vs[2])
    JOBRQ =  vs[3][1]
    RCOND = parse(Float64, replace(vs[4],'D'=>'E'))
    SVLMAX = parse(Float64, replace(vs[5],'D'=>'E'))
    if ( N<0 || N>NMAX )
        @error "Illegal N=$N"
    end
    if ( M<0 || M>MMAX )
        @error "Illegal M=$M"
    end

    vs = String[]
    _isz,_jsz = (M,N)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       A[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end
    close(f)

    JPVT[1:M] .= 0
# interp call 1

    RANK, INFO = SLICOT.mb03pd!(JOBRQ, M, N, A, JPVT, RCOND, SVLMAX, TAU, SVAL, LDWORK)
    @test INFO == 0
    println(io, "RANK = $RANK")

# interp output 1
    println(io, "JPVT:")
    _nr = M
    show(io, "text/plain", JPVT[1:_nr])
    println(io)

# interp output 2
    println(io, "SVAL:")
    _nr = 3
    show(io, "text/plain", SVAL[1:_nr])
    println(io)

end # run_mb03pd()


text:
##########
# This is more of a test case than a proper working example
##########

using ConstrainedDynamics

# Parameters
ex = [1.;0.;0.]

l1 = 1.0 * 0.3
d = .05
b1 = Box(l1, d, d, l1, color = RGBA(1., 1., 0.))
b2 = Box(.1, 1., 1., 2., color = RGBA(1., 0., 0.))

vert11 = [l1 / 2;0.;0.]
vert12 = -vert11

# Initial orientation
phi1 = 0.
q1 = Quaternion(RotX(phi1))

# Links
origin = Origin{Float64}()
link1 = Body(b1)
link2 = Body(b2)

# Constraints
socket0to1 = EqualityConstraint(Spherical(origin, link1, zeros(3), vert11))
joint1to5 = EqualityConstraint(Revolute(link1, link2, vert12, zeros(3), ex))

links = [link1;link2]
constraints = [socket0to1;joint1to5]
shapes = [b1;b2]

mech = Mechanism(origin, links, constraints;tend = 10.0,Δt = 0.001, shapes = shapes)
setPosition!(mech,origin,link1,p2 = vert11,Δq = q1)
setPosition!(mech,link1,link2,p1 = vert12,Δq = q1)
setVelocity!(mech,link1,ω = [50.;0;0])
setVelocity!(mech,link1,link2,p1 = vert12)

simulate!(mech,save = true)
visualize!(mech)


text:
using Libui
using ApplicationBuilderAppUtils

#const progressbar::Ptr{uiProgressBar} = C_NULL
#const spinbox::Ptr{uiProgressBar} = C_NULL
#const slider::Ptr{uiProgressBar} = C_NULL
#

function controlgallery()
	o = Libui.uiInitOptions(0)
	err = uiInit(Ref(o))

	menu = uiNewMenu("File")
	item = uiMenuAppendItem(menu, "Open")
	uiMenuItemOnClicked(item, cfunction(openClicked, Nothing, (Ptr{uiWindow},)), C_NULL)
	item = uiMenuAppendItem(menu, "Save")
	uiMenuItemOnClicked(item, cfunction(saveClicked, Nothing, (Ptr{uiWindow},)), C_NULL)
	item = uiMenuAppendQuitItem(menu)
	uiOnShouldQuit(cfunction(shouldQuit, Ptr{Nothing}, (Ptr{uiWindow},)), C_NULL)

	menu = uiNewMenu("Edit")
	item = uiMenuAppendCheckItem(menu, "Checkable Item")
	uiMenuAppendSeparator(menu)
	item = uiMenuAppendItem(menu, "Disabled Item")
	uiMenuItemDisable(item)
	item = uiMenuAppendPreferencesItem(menu)

	menu = uiNewMenu("Help")
	item = uiMenuAppendItem(menu, "Help")
	item = uiMenuAppendAboutItem(menu)

	mainwin = uiNewWindow("libui Control Gallery", 640, 480, 1)
	uiWindowSetMargined(mainwin, 1)
	uiWindowOnClosing(mainwin, cfunction(onClosing, Ptr{Nothing}, (Ptr{uiWindow},)), C_NULL)

	box = uiNewVerticalBox()
	uiBoxSetPadded(box, 1)
	uiWindowSetChild(mainwin, convert(Ptr{uiControl}, box))

	hbox = uiNewHorizontalBox()
	uiBoxSetPadded(hbox, 1)
	uiBoxAppend(box, convert(Ptr{uiControl}, hbox), 1)

	group = uiNewGroup("Basic Controls")
	uiGroupSetMargined(group, 1)
	uiBoxAppend(hbox, uiControl_(group), 0)

	inner = uiNewVerticalBox()
	uiBoxSetPadded(inner, 1)
	uiGroupSetChild(group, uiControl_(inner))

	uiBoxAppend(inner, uiControl_(uiNewButton("Button")), 0)
	uiBoxAppend(inner, uiControl_(uiNewCheckbox("Checkbox")), 0)
	entry = uiNewEntry()
	uiEntrySetText(entry, "Entry")
	uiBoxAppend(inner, uiControl_(entry), 0)
	uiBoxAppend(inner, uiControl_(uiNewLabel("Label")), 0)

	uiBoxAppend(inner, uiControl_(uiNewHorizontalSeparator()), 0)
	uiBoxAppend(inner, uiControl_(uiNewDatePicker()), 0)
	uiBoxAppend(inner, uiControl_(uiNewTimePicker()), 0)
	uiBoxAppend(inner, uiControl_(uiNewDateTimePicker()), 0)
	uiBoxAppend(inner, uiControl_(uiNewFontButton()), 0)
	uiBoxAppend(inner, uiControl_(uiNewColorButton()), 0)

	inner2 = uiNewVerticalBox()
	uiBoxSetPadded(inner2, 1)
	uiBoxAppend(hbox, uiControl_(inner2), 1)

	group = uiNewGroup("Numbers")
	uiGroupSetMargined(group, 1)
	uiBoxAppend(inner2, uiControl_(group), 0)

	inner = uiNewVerticalBox()
	uiBoxSetPadded(inner, 1)
	uiGroupSetChild(group, uiControl_(inner))

	uiSpinboxOnChanged(spinbox, cfunction(onSpinboxChanged, Nothing, (Ptr{uiSpinbox},)), C_NULL)
	uiBoxAppend(inner, uiControl_(spinbox), 0)

	uiSliderOnChanged(slider, cfunction(onSliderChanged, Nothing, (Ptr{uiSlider},)), C_NULL)
	uiBoxAppend(inner, uiControl_(slider), 0)

	uiBoxAppend(inner, uiControl_(progressbar), 0)

	group = uiNewGroup("Lists")
	uiGroupSetMargined(group, 1)
	uiBoxAppend(inner2, uiControl_(group), 0)

	inner = uiNewVerticalBox()
	uiBoxSetPadded(inner, 1)
	uiGroupSetChild(group, uiControl_(inner))

	cbox = uiNewCombobox()
	uiComboboxAppend(cbox, "Combobox Item 1")
	uiComboboxAppend(cbox, "Combobox Item 2")
	uiComboboxAppend(cbox, "Combobox Item 3")
	uiBoxAppend(inner, uiControl_(cbox), 0)

	ecbox = uiNewEditableCombobox()
	uiEditableComboboxAppend(ecbox, "Editable Item 1")
	uiEditableComboboxAppend(ecbox, "Editable Item 2")
	uiEditableComboboxAppend(ecbox, "Editable Item 3")
	uiBoxAppend(inner, uiControl_(ecbox), 0)

	rb = uiNewRadioButtons()
	uiRadioButtonsAppend(rb, "Radio Button 1")
	uiRadioButtonsAppend(rb, "Radio Button 2")
	uiRadioButtonsAppend(rb, "Radio Button 3")
	uiBoxAppend(inner, uiControl_(rb), 1)

	tab = uiNewTab()
	uiTabAppend(tab, "Page 1", uiControl_(uiNewHorizontalBox()))
	uiTabAppend(tab, "Page 2", uiControl_(uiNewHorizontalBox()))
	uiTabAppend(tab, "Page 3", uiControl_(uiNewHorizontalBox()))
	uiBoxAppend(inner2, uiControl_(tab), 1)

	uiControlShow(convert(Ptr{uiControl}, mainwin))
	uiMain()
	uiUninit()
end

# ---------------------------------------------------------------------------------------------
function onSpinboxChanged(spinbox::Ptr{uiSpinbox})
	update(uiSpinboxValue(spinbox))
	return
end

# ---------------------------------------------------------------------------------------------
function onSliderChanged(slider::Ptr{uiSlider})
	update(uiSliderValue(slider))
	return nothing
end

# ---------------------------------------------------------------------------------------------
function openClicked(mainwin::Ptr{uiWindow})
	filename = uiOpenFile(mainwin)
	if (filename == C_NULL)
		uiMsgBoxError(mainwin, "No file selected", "Don't be alarmed!")
		return
	end
	uiMsgBox(mainwin, "File selected", filename)
	uiFreeText(filename)
	return nothing
end

# ---------------------------------------------------------------------------------------------
function saveClicked(mainwin::Ptr{uiWindow})
	filename = uiSaveFile(mainwin)
	if (filename == C_NULL)
		uiMsgBoxError(mainwin, "No file selected", "Don't be alarmed!")
		return
	end
	uiMsgBox(mainwin, "File selected (don't worry, it's still there)", filename)
	uiFreeText(filename)
end

# ---------------------------------------------------------------------------------------------
function update(value::Integer)
	uiSpinboxSetValue(spinbox, value)
	uiSliderSetValue(slider, value)
	uiProgressBarSetValue(progressbar, value)
end

Base.@ccallable function julia_main(ARGS::Vector{String})::Cint
    ApplicationBuilderAppUtils.cd_to_bundle_resources()
    global progressbar = uiNewProgressBar()
    global spinbox = uiNewSpinbox(0, 100)
    global slider = uiNewSlider(0, 100)

    controlgallery()

    return 0
end


text:
"TruthTables exports the macro @truth_table"
module TruthTables

## For pretty printing
using NamedArrays: NamedArray

export @truth_table

¬(p::Bool) = !p
∧(p::Bool, q::Bool) = p && q
∨(p::Bool, q::Bool) = p || q
⟹(p::Bool, q::Bool) = if p; q else true end
⟺(p::Bool, q::Bool) = p === q

"""
    @truth_table predicate

Return the truth table of `predicate`.

The `predicate` can contain built-in boolean operators, as well as user-defined boolean-valued functions.
The truth table is a NamedArray whose columns are variables and sub-expressions of the input predicate.
"""
macro truth_table(expr)
    sub_exprs, vars = get_sub_exprs(expr)
    col_names = string.([vars..., sub_exprs...])
    c = gensym(:c)
    vals = gensym(:vals)

    ## escape user-defined functions
    sub_exprs = escape_calls.(sub_exprs)

    ## Form expression
    ex = :($vals[$c+=1, :] = [$(vars...), $(sub_exprs...)])
    for i in reverse(vars)
        ex = :(for $i in [true, false]; $ex end)
    end

    quote
        $c = 0
        $vals = Matrix{Bool}(undef, 2^length($vars), length($vars) + length($sub_exprs))
        $ex
        NamedArray($vals, (1:2^length($vars), $col_names), ("", ""))
    end
end

function get_sub_exprs(expr)
    sub_exprs = Expr[]
    vars = Symbol[]

    ## Traverse expression using multiple dispatch
    function walk!(ex::Expr)
        n = ex.head == :call ? 2 : 1    # Ignore function names
        walk!.(ex.args[n:end])          # Mind the dot
        ex ∉ sub_exprs && push!(sub_exprs, ex)
    end
    walk!(var::Symbol) = var ∉ vars && push!(vars, var)
    walk!(::Bool) = nothing

    walk!(expr)
    return sub_exprs, vars
end

escape_calls(s) = s
function escape_calls(expr::Expr)
    ex = expr
    if ex.head == :call && ex.args[1] isa Symbol && !isdefined(@__MODULE__, ex.args[1])
        ex.args[1] = esc(ex.args[1])
    end
    ex.args[2:end] = escape_calls.(ex.args[2:end])
    ex
end

end #module

max_stars_repo_name: savq/TruthTables

max_issues_repo_name: savq/TruthTables

max_forks_repo_name: savq/TruthTables


text:
using CompScienceMeshes
using Test

p1 = point(1,0,0)
p2 = point(0,1,0)
p3 = point(0,0,1)
p4 = point(0,0,0)

tet = simplex(p1,p2,p3,p4)
c = neighborhood(tet, point(1,1,1)/4)

@test cartesian(c) ≈ point(0.25, 0.25, 0.25)
@test parametric(c) ≈ point(0.25, 0.25, 0.25)

@test length(c) == 3
@test c[1] ≈ 0.25
@test c[2] ≈ 0.25
@test c[3] ≈ 0.25

@test jacobian(1.234) ≈ 1

@test tangents(c,1) ≈ point(1,0,0)
@test tangents(c,2) ≈ point(0,1,0)
@test tangents(c,3) ≈ point(0,0,1)

tri = simplex(p1,p2,p4)
c_tri = neighborhood(tri, point(1,1)/3)

@test normal(c_tri) ≈ point(0,0,1)


text:
using Distributed, Plots, StatsPlots

N = 4
if nworkers() < N
    addprocs(N, exeflags="--project=$(Base.active_project())")
end

@everywhere begin
    using DataFrames, CSV
    using Streamfall, BlackBoxOptim

    HERE = @__DIR__
    DATA_PATH = joinpath(HERE, "../test/data/cotter/")

    # Load observations
    date_format = "YYYY-mm-dd"
    obs_data = CSV.File(joinpath(DATA_PATH, "climate/CAMELS-AUS_410730.csv"),
                        comment="#",
                        dateformat=date_format) |> DataFrame

    hist_streamflow = obs_data[:, ["Date", "410730_Q"]]
    climate_data = obs_data[:, ["Date", "410730_P", "410730_PET"]]
    climate = Climate(climate_data, "_P", "_PET")

    burn_in = 366  # 1 year burn-in period

    # Create objective function to minimize (here we use Normalized KGE')
    func = (obs, sim) -> 1.0 - Streamfall.NmKGE(obs[burn_in:end], sim[burn_in:end])
    opt_func = (node) -> calibrate!(node, climate, hist_streamflow[:, "410730_Q"]; metric=func, MaxTime=180)
end

# Create individual nodes
hymod_node = create_node(SimpleHyModNode, "410730", 129.2)
gr4j_node = create_node(GR4JNode, "410730", 129.2)
symhyd_node = create_node(SYMHYDNode, "410730", 129.2)
ihacres_node = create_node(BilinearNode, "410730", 129.2)

# Calibrate each node separately using multiprocessing
node_names = ["HyMod", "GR4J", "SYMHYD", "IHACRES"]
node_list = [hymod_node, gr4j_node, symhyd_node, ihacres_node]
result = pmap(opt_func, node_list)

# Create comparison plot
Qo = hist_streamflow[:, "410730_Q"]
Qo_burn = Qo[burn_in:end]
res_plots = []

for ((res, opt), node, n_name) in zip(result, node_list, node_names)
    update_params!(node, best_candidate(res)...)
    reset!(node)
    run_node!(node, climate)

    node_burn = node.outflow[burn_in:end]

    res_plot = quickplot(Qo, node, climate, n_name; burn_in=366)
    push!(res_plots, res_plot)
end

combined_plot = plot(
    [rplt for rplt in res_plots]...,
    layout=(length(node_list),1),
    size=(950,450*length(node_list)),
)

display(combined_plot)

# savefig("multi_model_comparison.png")


text:
function onlineplot!(fig, axis_dict, stats::AbstractVector, iter, data, variable, i)
    for (j, stat) in enumerate(stats)
        axis_dict[(variable, stat)] = fig[i, j] = Axis(fig, title="$(name(stat))")
        limits!(axis_dict[(variable, stat)], 0.0, 10.0, -1.0, 1.0)
        onlineplot!(axis_dict[(variable, stat)], stat, iter, data[variable], data[:iter], i, j)
        tight_ticklabel_spacing!(axis_dict[(variable, stat)])
    end
end

function onlineplot!(axis, stat::Symbol, args...)
    onlineplot!(axis, Val(stat), args...)
end

onlineplot!(axis, ::Val{:mean}, args...) = onlineplot!(axis, Mean(), args...)

onlineplot!(axis, ::Val{:var}, args...) = onlineplot!(axis, Variance(), args...)

onlineplot!(axis, ::Val{:autocov}, args...) = onlineplot!(axis, AutoCov(20), args...)

onlineplot!(axis, ::Val{:hist}, args...) = onlineplot!(axis, KHist(50, Float32), args...)

# Generic fallback for OnlineStat objects
function onlineplot!(axis, stat::T, iter, data, iterations, i, j) where {T<:OnlineStat}
    window = data.b
    @eval TStat = $(nameof(T))
    stat = Observable(TStat(Float32))
    on(iter) do _
        stat[] = fit!(stat[], last(value(data)))
    end
    statvals = Observable(MovingWindow(window, Float32))
    on(stat) do s
        statvals[] = fit!(statvals[], Float32(value(s)))
    end
    statpoints = map!(Observable(Point2f0.([0], [0])), statvals) do v
        Point2f0.(value(iterations), value(v))
    end
    lines!(axis, statpoints, color = std_colors[i], linewidth = 3.0)
end

function onlineplot!(axis, ::Val{:trace}, iter, data, iterations, i, j)
    trace = map!(Observable([Point2f0(0, 0)]), iter) do _
        Point2f0.(value(iterations), value(data))
    end
    lines!(axis, trace, color = std_colors[i]; linewidth = 3.0)
end

function onlineplot!(axis, stat::KHist, iter, data, iterations, i, j)
    nbins = stat.k
    stat = Observable(KHist(nbins, Float32))
    on(iter) do _
        stat[] = fit!(stat[], last(value(data)))
    end
    hist_vals = Node(Point2f0.(collect(range(0f0, 1f0, length=nbins)), zeros(Float32, nbins)))
    on(stat) do h
        edges, weights = OnlineStats.xy(h)
        weights = nobs(h) > 1 ? weights / OnlineStats.area(h) : weights
        hist_vals[] = Point2f0.(edges, weights)
    end
    barplot!(axis, hist_vals; color=std_colors[i])
    # barplot!(axis, rand(4), rand(4))
end

function expand_extrema(xs)
    xmin, xmax = xs
    diffx = xmax - xmin
    xmin = xmin - 0.1 * abs(diffx)
    xmax = xmax + 0.1 * abs(diffx)
    return (xmin, xmax)
end

function onlineplot!(axis, ::Val{:kde}, iter, data, iterations, i, j)
    interpkde = Observable(InterpKDE(kde([1f0])))
    on(iter) do _
        interpkde[] = InterpKDE(kde(value(data)))
    end
    xs = Observable(range(0, 2, length=10))
    on(iter) do _
        xs[] = range(expand_extrema(extrema(value(data)))..., length = 200)
    end
    kde_pdf = lift(xs) do xs
        pdf.(Ref(interpkde[]), xs)
    end
    lines!(axis, xs, kde_pdf, color = std_colors[i], linewidth = 3.0)
end

name(s::Val{:histkde}) = "Hist. + KDE"

function onlineplot!(axis, ::Val{:histkde}, iter, data, iterations, i, j)
    onlineplot!(axis, KHist(50), iter, data, iterations, i, j)
    onlineplot!(axis, Val(:kde), iter, data, iterations, i, j)
end

function onlineplot!(axis, stat::AutoCov, iter, data, iterations, i, j)
    b = length(stat.cross)
    stat = Observable(AutoCov(b, Float32))
    on(iter) do _
        stat[] = fit!(stat[], last(value(data)))
    end
    statvals = map!(Observable(zeros(Float32, b + 1)), stat) do s
        value(s)
    end
    scatter!(axis, Point2f0.([0.0, b], [-0.1, 1.0]), markersize = 0.0, color = RGBA(0.0, 0.0, 0.0, 0.0)) # Invisible points to keep limits fixed
    lines!(axis, 0:b, statvals, color = std_colors[i], linewidth = 3.0)
end


text:
function tugbot_visualizer!(vis::Visualizer, mech::Mechanism;)
    material = MeshPhongMaterial(color=Colors.RGBA(0.0, 0.0, 0.0, 1.0))
    drone = MeshFileObject(joinpath(dirname(pathof(DojoEnvironments)), "tugbot/deps/drone.obj"))
    setobject!(vis["robot"]["bodies"]["drone__id_1"][:mesh], drone)
    settransform!(vis["robot"]["bodies"]["drone__id_1"][:mesh], MeshCat.LinearMap(RotX(π/2)))
end


text:
## Mapping, until Mimi can handle string parameters:
## String -> Int64
## "zero" -> 0, "rcp26" -> 26, "rcp45" -> 45, "rcp85" -> 85, "rcpw" -> 10, "rcp26extra" -> 260
## "ssp1" -> 1, "ssp2" -> 2, "ssp5" -> 5, "sspw" -> 10, "ssp234" -> 234
## $(rcp) -> rcp$(rcp), $(ssp) -> ssp$(ssp)

@defcomp RCPSSPScenario begin
    region = Index()

    rcp = Parameter{Int64}() # like rcp26
    ssp = Parameter{Int64}() # like ssp1

    y_year = Parameter(index=[time], unit="year")
    weight_scenarios = Parameter(unit="%") # from -100% to 100%, only used for sspw, rcpw

    extra_abate_rate = Parameter(unit="%/year") # only used for rcp26extra
    extra_abate_start = Parameter(unit="year")
    extra_abate_end = Parameter(unit="year")

    # RCP scenario values
    er_CO2emissionsgrowth = Variable(index=[time,region], unit="%")
    er_CH4emissionsgrowth = Variable(index=[time,region], unit="%")
    er_N2Oemissionsgrowth = Variable(index=[time,region], unit="%")
    er_LGemissionsgrowth = Variable(index=[time,region], unit="%")
    pse_sulphatevsbase = Variable(index=[time, region], unit="%")
    exf_excessforcing = Variable(index=[time], unit="W/m2")

    extra_abate_compound = Variable(index=[time])

    # SSP scenario values
    popgrw_populationgrowth = Variable(index=[time, region], unit="%/year") # From p.32 of Hope 2009
    grw_gdpgrowthrate = Variable(index=[time, region], unit="%/year") # From p.32 of Hope 2009

    function init(p, v, d)
        # Set the RCP values
        if p.rcp == 10
            v.er_CO2emissionsgrowth[:, :] =
                weighted_scenario(readpagedata(nothing, "data/rcps/rcp26_co2.csv"), readpagedata(nothing, "data/rcps/rcp45_co2.csv"),
                                  readpagedata(nothing, "data/rcps/rcp85_co2.csv"), p.weight_scenarios)
            v.er_CH4emissionsgrowth[:, :] =
                weighted_scenario(readpagedata(nothing, "data/rcps/rcp26_ch4.csv"), readpagedata(nothing, "data/rcps/rcp45_ch4.csv"),
                                  readpagedata(nothing, "data/rcps/rcp85_ch4.csv"), p.weight_scenarios)
            v.er_N2Oemissionsgrowth[:, :] =
                weighted_scenario(readpagedata(nothing, "data/rcps/rcp26_n2o.csv"), readpagedata(nothing, "data/rcps/rcp45_n2o.csv"),
                                  readpagedata(nothing, "data/rcps/rcp85_n2o.csv"), p.weight_scenarios)
            v.er_LGemissionsgrowth[:, :] =
                weighted_scenario(readpagedata(nothing, "data/rcps/rcp26_lin.csv"), readpagedata(nothing, "data/rcps/rcp45_lin.csv"),
                                  readpagedata(nothing, "data/rcps/rcp85_lin.csv"), p.weight_scenarios)
            v.pse_sulphatevsbase[:, :] =
                weighted_scenario(readpagedata(nothing, "data/rcps/rcp26_sulph.csv"), readpagedata(nothing, "data/rcps/rcp45_sulph.csv"),
                                  readpagedata(nothing, "data/rcps/rcp85_sulph.csv"), p.weight_scenarios)
            v.exf_excessforcing[:] =
                weighted_scenario(readpagedata(nothing, "data/rcps/rcp26_excess.csv"), readpagedata(nothing, "data/rcps/rcp45_excess.csv"),
                                  readpagedata(nothing, "data/rcps/rcp85_excess.csv"), p.weight_scenarios)
        elseif p.rcp == 260
            # Fill in within run_timestep
        elseif p.rcp == 0
            v.er_CO2emissionsgrowth[:, :] .= 0.
            v.er_CH4emissionsgrowth[:, :] .= 0.
            v.er_N2Oemissionsgrowth[:, :] .= 0.
            v.er_LGemissionsgrowth[:, :] .= 0.
            v.pse_sulphatevsbase[:, :] .= 0.
            v.exf_excessforcing[:] .= 0.
        else
            v.er_CO2emissionsgrowth[:, :] = readpagedata(nothing, "data/rcps/rcp$(p.rcp)_co2.csv")
            v.er_CH4emissionsgrowth[:, :] = readpagedata(nothing, "data/rcps/rcp$(p.rcp)_ch4.csv")
            v.er_N2Oemissionsgrowth[:, :] = readpagedata(nothing, "data/rcps/rcp$(p.rcp)_n2o.csv")
            v.er_LGemissionsgrowth[:, :] = readpagedata(nothing, "data/rcps/rcp$(p.rcp)_lin.csv")
            v.pse_sulphatevsbase[:, :] = readpagedata(nothing, "data/rcps/rcp$(p.rcp)_sulph.csv")
            v.exf_excessforcing[:] = readpagedata(nothing, "data/rcps/rcp$(p.rcp)_excess.csv")
        end

        # Set the SSP values
        if p.ssp == 234 || p.ssp == 10
            v.popgrw_populationgrowth[:, :] = (readpagedata(nothing, "data/ssps/ssp2_pop_rate.csv") +
                                   readpagedata(nothing, "data/ssps/ssp3_pop_rate.csv") +
                                   readpagedata(nothing, "data/ssps/ssp4_pop_rate.csv")) / 3
            v.grw_gdpgrowthrate[:, :] = (readpagedata(nothing, "data/ssps/ssp2_gdp_rate.csv") +
                                         readpagedata(nothing, "data/ssps/ssp3_gdp_rate.csv") +
                                         readpagedata(nothing, "data/ssps/ssp4_gdp_rate.csv")) / 3
            if p.ssp == 10
                v.popgrw_populationgrowth[:, :] =
                    weighted_scenario(readpagedata(nothing, "data/ssps/ssp1_pop_rate.csv"), v.popgrw_populationgrowth[:, :],
                                      readpagedata(nothing, "data/ssps/ssp5_pop_rate.csv"), p.weight_scenarios)
                v.grw_gdpgrowthrate[:, :] =
                    weighted_scenario(readpagedata(nothing, "data/ssps/ssp1_gdp_rate.csv"), v.grw_gdpgrowthrate[:, :],
                                      readpagedata(nothing, "data/ssps/ssp5_gdp_rate.csv"), p.weight_scenarios)
            end
        else
            v.popgrw_populationgrowth[:, :] = readpagedata(nothing, "data/ssps/ssp$(p.ssp)_pop_rate.csv")
            v.grw_gdpgrowthrate[:, :] = readpagedata(nothing, "data/ssps/ssp$(p.ssp)_gdp_rate.csv")
        end
    end

    function run_timestep(p, v, d, t)
        # Only used for rcp26extra
        if p.rcp == 260
            if is_first(t)
                duration = 5
            else
                duration = p.y_year[t] - p.y_year[t - 1]
            end

            extra_abate_period = ifelse(p.y_year[t] <= p.extra_abate_start || p.y_year[t] > p.extra_abate_end, 1.,
                                        (1 - p.extra_abate_rate / 100.)^duration)
            if is_first(t)
                v.extra_abate_compound[t] = extra_abate_period
            else
                v.extra_abate_compound[t] = extra_abate_period * v.extra_abate_compound[t - 1]
            end

            er_rcp26_CO2emissionsgrowth = readpagedata(nothing, "data/rcps/rcp26_co2.csv")
            er_rcp26_CH4emissionsgrowth = readpagedata(nothing, "data/rcps/rcp26_ch4.csv")
            er_rcp26_N2Oemissionsgrowth = readpagedata(nothing, "data/rcps/rcp26_n2o.csv")
            er_rcp26_LGemissionsgrowth = readpagedata(nothing, "data/rcps/rcp26_lin.csv")
            pse_rcp26_sulphatevsbase = readpagedata(nothing, "data/rcps/rcp26_sulph.csv")
            exf_rcp26_excessforcing = readpagedata(nothing, "data/rcps/rcp26_excess.csv")

            v.er_CO2emissionsgrowth[t, :] = (er_rcp26_CO2emissionsgrowth[t.t, :] - er_rcp26_CO2emissionsgrowth[p.y_year[:] .== 2100, :][:]) * v.extra_abate_compound[t] .+ er_rcp26_CO2emissionsgrowth[p.y_year[:] .== 2100, :][:]
            v.er_CH4emissionsgrowth[t, :] = er_rcp26_CH4emissionsgrowth[t.t, :] * v.extra_abate_compound[t]
            v.er_N2Oemissionsgrowth[t, :] = er_rcp26_N2Oemissionsgrowth[t.t, :] * v.extra_abate_compound[t]
            v.er_LGemissionsgrowth[t, :] = er_rcp26_LGemissionsgrowth[t.t, :] * v.extra_abate_compound[t]
            v.pse_sulphatevsbase[t, :] = pse_rcp26_sulphatevsbase[t.t, :] * v.extra_abate_compound[t]
            v.exf_excessforcing[t, :] = exf_rcp26_excessforcing[t.t, :] * v.extra_abate_compound[t]
        end
    end
end

function weighted_scenario(lowscen, medscen, highscen, weight)
    lowscen * .25 * (1 - weight / 100)^2 +
        medscen * .5 * (1 - (weight / 100)^2) +
        highscen * .25 * (1 + weight / 100)^2
end

function addrcpsspscenario(model::Model, scenario::String)
    rcpsspscenario = add_comp!(model, RCPSSPScenario)

    # Default parameters
    rcpsspscenario[:y_year] = Mimi.dim_keys(model.md, :time)
    rcpsspscenario[:extra_abate_rate] = 0
    rcpsspscenario[:extra_abate_start] = 2015
    rcpsspscenario[:extra_abate_end] = 2100
    rcpsspscenario[:weight_scenarios] = 0.

    if scenario == "Zero Emissions & SSP1"
        rcpsspscenario[:rcp] = 0
        rcpsspscenario[:ssp] = 1
    elseif scenario == "1.5 degC Target"
        rcpsspscenario[:rcp] = 260
        rcpsspscenario[:ssp] = 1
        rcpsspscenario[:extra_abate_rate] = 4.053014079712271
        rcpsspscenario[:extra_abate_start] = 2020
        rcpsspscenario[:extra_abate_end] = 2100
    elseif scenario == "2 degC Target"
        rcpsspscenario[:rcp] = 260
        rcpsspscenario[:ssp] = 1
        rcpsspscenario[:extra_abate_rate] = 0.2418203462401034
        rcpsspscenario[:extra_abate_start] = 2020
        rcpsspscenario[:extra_abate_end] = 2100
    elseif scenario == "2.5 degC Target"
        rcpsspscenario[:rcp] = 10
        rcpsspscenario[:ssp] = 10
        rcpsspscenario[:weight_scenarios] = -69.69860118334117
    elseif scenario == "NDCs"
        rcpsspscenario[:rcp] = 10
        rcpsspscenario[:ssp] = 10
        rcpsspscenario[:weight_scenarios] = -14.432092365610856
    elseif scenario == "NDCs Partial"
        rcpsspscenario[:rcp] = 10
        rcpsspscenario[:ssp] = 10
        rcpsspscenario[:weight_scenarios] = 10.318413035360622
    elseif scenario == "BAU"
        rcpsspscenario[:rcp] = 10
        rcpsspscenario[:ssp] = 10
        rcpsspscenario[:weight_scenarios] = 51.579276825415874
    elseif scenario == "RCP2.6 & SSP1"
        rcpsspscenario[:rcp] = 26
        rcpsspscenario[:ssp] = 1
    elseif scenario == "RCP4.5 & SSP2"
        rcpsspscenario[:rcp] = 45
        rcpsspscenario[:ssp] = 2
    elseif scenario == "RCP8.5 & SSP5"
        rcpsspscenario[:rcp] = 85
        rcpsspscenario[:ssp] = 5
    elseif scenario == "RCP8.5 & SSP2"
        rcpsspscenario[:rcp] = 85
        rcpsspscenario[:ssp] = 2
    else
        error("Unknown scenario")
    end

    rcpsspscenario
end


text:
# A versatile graph type
#
# It implements edge_list, adjacency_list and incidence_list
#

type GenericGraph{V,E,VList,EList,IncList} <: AbstractGraph{V,E}
    is_directed::Bool
    vertices::VList     # an indexable container of vertices
    edges::EList        # an indexable container of edges
    finclist::IncList   # forward incidence list
    binclist::IncList   # backward incidence list
end

@graph_implements GenericGraph vertex_list edge_list vertex_map edge_map
@graph_implements GenericGraph bidirectional_adjacency_list bidirectional_incidence_list

# SimpleGraph:
#   V:          Int
#   E:          IEdge
#   VList:      Range1{Int}
#   EList:      Vector{IEdge}
#   AdjList:    Vector{Vector{Int}}
#   IncList:    Vector{Vector{IEdge}}
#
typealias SimpleGraph GenericGraph{Int,IEdge,Range1{Int},Vector{IEdge},Vector{Vector{IEdge}}}

typealias Graph{V,E} GenericGraph{V,E,Vector{V},Vector{E},Vector{Vector{E}}}

# construction

simple_graph(n::Integer; is_directed::Bool=true) =
    SimpleGraph(is_directed,
                1:int(n),  # vertices
                IEdge[],   # edges
                multivecs(IEdge, n), # finclist
                multivecs(IEdge, n)) # binclist

function graph{V,E}(vs::Vector{V}, es::Vector{E}; is_directed::Bool=true)
    n = length(vs)
    g = Graph{V,E}(is_directed, vs, E[], multivecs(E, n), multivecs(E, n))
    for e in es
        add_edge!(g, e)
    end
    return g
end

# required interfaces

is_directed(g::GenericGraph) = g.is_directed

num_vertices(g::GenericGraph) = length(g.vertices)
vertices(g::GenericGraph) = g.vertices

num_edges(g::GenericGraph) = length(g.edges)
edges(g::GenericGraph) = g.edges

vertex_index{V<:ProvidedVertexType}(v::V, g::GenericGraph{V}) = vertex_index(v)
edge_index{V,E}(e::E, g::GenericGraph{V,E}) = edge_index(e)

out_edges{V}(v::V, g::GenericGraph{V}) = g.finclist[vertex_index(v, g)]
out_degree{V}(v::V, g::GenericGraph{V}) = length(out_edges(v, g))
out_neighbors{V}(v::V, g::GenericGraph{V}) = TargetIterator(g, out_edges(v, g))

in_edges{V}(v::V, g::GenericGraph{V}) = g.binclist[vertex_index(v, g)]
in_degree{V}(v::V, g::GenericGraph{V}) = length(in_edges(v, g))
in_neighbors{V}(v::V, g::GenericGraph{V}) = SourceIterator(g, in_edges(v, g))

# mutation

function add_vertex!{V}(g::GenericGraph{V}, v::V)
    push!(g.vertices, v)
    push!(g.finclist, Int[])
    push!(g.binclist, Int[])
    v
end
add_vertex!{V}(g::GenericGraph{V}, x) = add_vertex!(g, make_vertex(g, x))

function add_edge!{V,E}(g::GenericGraph{V,E}, u::V, v::V, e::E)
    # add an edge e between u and v
    ui = vertex_index(u, g)::Int
    vi = vertex_index(v, g)::Int

    push!(g.edges, e)
    push!(g.finclist[ui], e)
    push!(g.binclist[vi], e)

    if !g.is_directed
        rev_e = revedge(e)
        push!(g.finclist[vi], rev_e)
        push!(g.binclist[ui], rev_e)
    end
    e
end

add_edge!{V,E}(g::GenericGraph{V,E}, e::E) = add_edge!(g, source(e, g), target(e, g), e)
add_edge!{V,E}(g::GenericGraph{V,E}, u::V, v::V) = add_edge!(g, u, v, make_edge(g, u, v))


text:
# Decompressor Codecs
# ====================

abstract type DecompressorCodec <: TranscodingStreams.Codec end

function Base.show(io::IO, codec::DecompressorCodec)
    print(io, summary(codec), "(windowbits=$(codec.windowbits))")
end

# Gzip
# ----

struct GzipDecompressor <: DecompressorCodec
    zstream::ZStream
    windowbits::Int
end

"""
    GzipDecompressor(;windowbits=$(Z_DEFAULT_WINDOWBITS), gziponly=false)

Create a gzip decompressor codec.

If `gziponly` is `false`, this codec can decompress the zlib format as well.

Arguments
---------
- `windowbits`: size of history buffer (8..15)
- `gziponly`: flag to inactivate data format detection
"""
function GzipDecompressor(;windowbits::Integer=Z_DEFAULT_WINDOWBITS, gziponly::Bool=false)
    if !(8 ≤ windowbits ≤ 15)
        throw(ArgumentError("windowbits must be within 8..15"))
    end
    return GzipDecompressor(ZStream(), windowbits+(gziponly ? 16 : 32))
end

const GzipDecompressorStream{S} = TranscodingStream{GzipDecompressor,S} where S<:IO

"""
    GzipDecompressorStream(stream::IO; kwargs...)

Create a gzip decompression stream (see `GzipDecompressor` for `kwargs`).
"""
function GzipDecompressorStream(stream::IO; kwargs...)
    x, y = splitkwargs(kwargs, (:windowbits, :gziponly))
    return TranscodingStream(GzipDecompressor(;x...), stream; y...)
end

# Zlib
# ----

struct ZlibDecompressor <: DecompressorCodec
    zstream::ZStream
    windowbits::Int
end

"""
    ZlibDecompressor(;windowbits=$(Z_DEFAULT_WINDOWBITS))

Create a zlib decompression codec.

Arguments
---------
- `windowbits`: size of history buffer (8..15) or (-15..-8) or 0
"""
function ZlibDecompressor(;windowbits::Integer=Z_DEFAULT_WINDOWBITS)
    if !(8 ≤ abs(windowbits) ≤ 15) || windowbits == 0
        throw(ArgumentError("|windowbits| must be within 8..15 or 0"))
    end
    return ZlibDecompressor(ZStream(), windowbits)
end

const ZlibDecompressorStream{S} = TranscodingStream{ZlibDecompressor,S} where S<:IO

"""
    ZlibDecompressorStream(stream::IO; kwargs...)

Create a deflate decompression stream (see `ZlibDecompressor` for `kwargs`).
"""
function ZlibDecompressorStream(stream::IO; kwargs...)
    x, y = splitkwargs(kwargs, (:windowbits,))
    return TranscodingStream(ZlibDecompressor(;x...), stream; y...)
end

# Deflate
# -------

struct DeflateDecompressor <: DecompressorCodec
    zstream::ZStream
    windowbits::Int
end

"""
    DeflateDecompressor(;windowbits=$(Z_DEFAULT_WINDOWBITS))

Create a deflate decompression codec.

Arguments
---------
- `windowbits`: size of history buffer (8..15)
"""
function DeflateDecompressor(;windowbits::Integer=Z_DEFAULT_WINDOWBITS)
    if !(8 ≤ windowbits ≤ 15)
        throw(ArgumentError("windowbits must be within 8..15"))
    end
    return DeflateDecompressor(ZStream(), -Int(windowbits))
end

const DeflateDecompressorStream{S} = TranscodingStream{DeflateDecompressor,S} where S<:IO

"""
    DeflateDecompressorStream(stream::IO; kwargs...)

Create a deflate decompression stream (see `DeflateDecompressor` for `kwargs`).
"""
function DeflateDecompressorStream(stream::IO; kwargs...)
    x, y = splitkwargs(kwargs, (:windowbits,))
    return TranscodingStream(DeflateDecompressor(;x...), stream; y...)
end

# Methods
# -------

function TranscodingStreams.initialize(codec::DecompressorCodec)
    code = inflate_init!(codec.zstream, codec.windowbits)
    if code != Z_OK
        zerror(codec.zstream, code)
    end
    return
end

function TranscodingStreams.finalize(codec::DecompressorCodec)
    zstream = codec.zstream
    if zstream.state != C_NULL
        code = inflate_end!(zstream)
        if code != Z_OK
            zerror(zstream, code)
        end
    end
    return
end

function TranscodingStreams.startproc(codec::DecompressorCodec, ::Symbol, error::Error)
    code = inflate_reset!(codec.zstream)
    if code == Z_OK
        return :ok
    else
        error[] = ErrorException(zlib_error_message(codec.zstream, code))
        return :error
    end
end

function TranscodingStreams.process(codec::DecompressorCodec, input::Memory, output::Memory, error::Error)
    zstream = codec.zstream
    zstream.next_in = input.ptr
    zstream.avail_in = input.size
    zstream.next_out = output.ptr
    zstream.avail_out = output.size
    code = inflate!(zstream, Z_NO_FLUSH)
    Δin = Int(input.size - zstream.avail_in)
    Δout = Int(output.size - zstream.avail_out)
    if code == Z_OK
        return Δin, Δout, :ok
    elseif code == Z_STREAM_END
        return Δin, Δout, :end
    else
        error[] = ErrorException(zlib_error_message(zstream, code))
        return Δin, Δout, :error
    end
end


text:
# Stopgaps, should have proper treatment later
Random.randn(ω::Ω, ::Type{Float64} = Float64) = normal(ω, 0, 1)
Random.randexp(ω::Ω, ::Type{Float64} = Float64) = exponential(ω, 1)


text:
using Distributions
using PyPlot

plt[:hist](randn(1000))


text:
module BoosterWhite

using ..Ahorn, Maple

@mapdef Entity "brokemiahelper/whitebooster" WhiteBooster(x::Integer, y::Integer)

const placements = Ahorn.PlacementDict(
    "Booster (White) (Brokemia Helper)" => Ahorn.EntityPlacement(
        WhiteBooster
    )
)

function boosterSprite(entity::WhiteBooster)
    return "objects/booster/booster00"
end

function Ahorn.selection(entity::WhiteBooster)
    x, y = Ahorn.position(entity)
    sprite = boosterSprite(entity)

    return Ahorn.getSpriteRectangle(sprite, x, y)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::WhiteBooster, room::Maple.Room)
    sprite = boosterSprite(entity)

    Ahorn.drawSprite(ctx, sprite, 0, 0)
end

end

max_stars_repo_name: Brokemia/BrokemiaHelper

max_issues_repo_name: Brokemia/BrokemiaHelper

max_forks_repo_name: Brokemia/BrokemiaHelper


text:
import Test
import DataFrames
import PredictMD

my_vector = Vector{Any}(undef, 5)
my_vector[1] = Float64(1.1)
my_vector[2] = Float64(2.2)
my_vector[3] = DataFrames.missing
my_vector[4] = Float64(4.4)
my_vector[5] = DataFrames.missing
Test.@test(eltype(my_vector) == Any)
my_vector_fixed = PredictMD.fix_type(my_vector)
Test.@test(eltype(my_vector_fixed) == Union{Float64, DataFrames.Missing})

Test.@test( PredictMD.fix_type(nothing) == nothing )

Test.@test( PredictMD.fix_type(3.14) == 3.14 )

dict_1 = Dict()
dict_1[Symbol(:x)] = Float64(1.1)
dict_1[Symbol(:y)] = Float64(2.2)
dict_1[Symbol(:z)] = Float64(3.3)
Test.@test(typeof(dict_1) <: Dict{Any, Any})
Test.@test(length(dict_1) == 3)

dict_2 = PredictMD.fix_type(dict_1)
Test.@test(typeof(dict_2) <: Dict{Symbol, Float64})
Test.@test(length(dict_2) == 3)
Test.@test(dict_1[:x] == 1.1)
Test.@test(dict_1[:y] == 2.2)
Test.@test(dict_1[:z] == 3.3)

dict_3 = PredictMD.fix_type(dict_2)
Test.@test(typeof(dict_3) <: Dict{Symbol, Float64})
Test.@test(length(dict_3) == 3)
Test.@test(dict_1[:x] == 1.1)
Test.@test(dict_1[:y] == 2.2)
Test.@test(dict_1[:z] == 3.3)

dict_4 = PredictMD.fix_type(dict_3)
Test.@test(typeof(dict_4) <: Dict{Symbol, Float64})
Test.@test(length(dict_4) == 3)
Test.@test(dict_1[:x] == 1.1)
Test.@test(dict_1[:y] == 2.2)
Test.@test(dict_1[:z] == 3.3)

Test.@test( PredictMD.fix_type(nothing) == nothing )
Test.@test( PredictMD.fix_type(nothing) == nothing )

vector_1 = []
push!(vector_1, Float64(1.0))
push!(vector_1, Float64(2.0))
push!(vector_1, Float64(3.0))
Test.@test(typeof(vector_1) <: Vector{Any})
Test.@test(length(vector_1) == 3)
Test.@test(size(vector_1) == (3,))

vector_2 = PredictMD.fix_type(vector_1)
Test.@test(typeof(vector_2) <: Vector{Float64})
Test.@test(length(vector_2) == 3)
Test.@test(size(vector_2) == (3,))
Test.@test(vector_1[1] == 1.0)
Test.@test(vector_1[2] == 2.0)
Test.@test(vector_1[3] == 3.0)

vector_3 = PredictMD.fix_type(vector_2)
Test.@test(typeof(vector_3) <: Vector{Float64})
Test.@test(length(vector_3) == 3)
Test.@test(size(vector_3) == (3,))
Test.@test(vector_1[1] == 1.0)
Test.@test(vector_1[2] == 2.0)
Test.@test(vector_1[3] == 3.0)

vector_4 = PredictMD.fix_type(vector_3)
Test.@test(length(vector_4) == 3)
Test.@test(size(vector_4) == (3,))
Test.@test(typeof(vector_4) <: Vector{Float64})
Test.@test(vector_1[1] == 1.0)
Test.@test(vector_1[2] == 2.0)
Test.@test(vector_1[3] == 3.0)

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

array_1 = Array{Any}(undef, 2,3,4)
array_1[1,1,1] = Float64(10)
array_1[1,1,2] = Float64(20)
array_1[1,1,3] = Float64(30)
array_1[1,1,4] = Float64(40)
array_1[1,2,1] = Float64(50)
array_1[1,2,2] = Float64(60)
array_1[1,2,3] = Float64(70)
array_1[1,2,4] = Float64(80)
array_1[1,3,1] = Float64(90)
array_1[1,3,2] = Float64(100)
array_1[1,3,3] = Float64(110)
array_1[1,3,4] = Float64(120)
array_1[2,1,1] = Float64(130)
array_1[2,1,2] = Float64(140)
array_1[2,1,3] = Float64(150)
array_1[2,1,4] = Float64(160)
array_1[2,2,1] = Float64(170)
array_1[2,2,2] = Float64(180)
array_1[2,2,3] = Float64(190)
array_1[2,2,4] = Float64(200)
array_1[2,3,1] = Float64(210)
array_1[2,3,2] = Float64(220)
array_1[2,3,3] = Float64(230)
array_1[2,3,4] = Float64(240)
Test.@test(typeof(array_1) <: Array{Any, 3})
Test.@test(length(array_1) == 24)
Test.@test(size(array_1) == (2,3,4,))

array_2 = PredictMD.fix_type(array_1)
Test.@test(typeof(array_2) <: Array{Float64, 3})
Test.@test(length(array_2) == 24)
Test.@test(size(array_2) == (2,3,4,))
Test.@test(array_2[1,1,1] == 10)
Test.@test(array_2[1,1,2] == 20)
Test.@test(array_2[1,1,3] == 30)
Test.@test(array_2[1,1,4] == 40)
Test.@test(array_2[1,2,1] == 50)
Test.@test(array_2[1,2,2] == 60)
Test.@test(array_2[1,2,3] == 70)
Test.@test(array_2[1,2,4] == 80)
Test.@test(array_2[1,3,1] == 90)
Test.@test(array_2[1,3,2] == 100)
Test.@test(array_2[1,3,3] == 110)
Test.@test(array_2[1,3,4] == 120)
Test.@test(array_2[2,1,1] == 130)
Test.@test(array_2[2,1,2] == 140)
Test.@test(array_2[2,1,3] == 150)
Test.@test(array_2[2,1,4] == 160)
Test.@test(array_2[2,2,1] == 170)
Test.@test(array_2[2,2,2] == 180)
Test.@test(array_2[2,2,3] == 190)
Test.@test(array_2[2,2,4] == 200)
Test.@test(array_2[2,3,1] == 210)
Test.@test(array_2[2,3,2] == 220)
Test.@test(array_2[2,3,3] == 230)
Test.@test(array_2[2,3,4] == 240)

array_3 = PredictMD.fix_type(array_2)
Test.@test(typeof(array_3) <: Array{Float64, 3})
Test.@test(length(array_3) == 24)
Test.@test(size(array_3) == (2,3,4,))
Test.@test(array_3[1,1,1] == 10)
Test.@test(array_3[1,1,2] == 20)
Test.@test(array_3[1,1,3] == 30)
Test.@test(array_3[1,1,4] == 40)
Test.@test(array_3[1,2,1] == 50)
Test.@test(array_3[1,2,2] == 60)
Test.@test(array_3[1,2,3] == 70)
Test.@test(array_3[1,2,4] == 80)
Test.@test(array_3[1,3,1] == 90)
Test.@test(array_3[1,3,2] == 100)
Test.@test(array_3[1,3,3] == 110)
Test.@test(array_3[1,3,4] == 120)
Test.@test(array_3[2,1,1] == 130)
Test.@test(array_3[2,1,2] == 140)
Test.@test(array_3[2,1,3] == 150)
Test.@test(array_3[2,1,4] == 160)
Test.@test(array_3[2,2,1] == 170)
Test.@test(array_3[2,2,2] == 180)
Test.@test(array_3[2,2,3] == 190)
Test.@test(array_3[2,2,4] == 200)
Test.@test(array_3[2,3,1] == 210)
Test.@test(array_3[2,3,2] == 220)
Test.@test(array_3[2,3,3] == 230)
Test.@test(array_3[2,3,4] == 240)

array_4 = PredictMD.fix_type(array_3)
Test.@test(typeof(array_4) <: Array{Float64, 3})
Test.@test(length(array_4) == 24)
Test.@test(size(array_4) == (2,3,4,))
Test.@test(array_4[1,1,1] == 10)
Test.@test(array_4[1,1,2] == 20)
Test.@test(array_4[1,1,3] == 30)
Test.@test(array_4[1,1,4] == 40)
Test.@test(array_4[1,2,1] == 50)
Test.@test(array_4[1,2,2] == 60)
Test.@test(array_4[1,2,3] == 70)
Test.@test(array_4[1,2,4] == 80)
Test.@test(array_4[1,3,1] == 90)
Test.@test(array_4[1,3,2] == 100)
Test.@test(array_4[1,3,3] == 110)
Test.@test(array_4[1,3,4] == 120)
Test.@test(array_4[2,1,1] == 130)
Test.@test(array_4[2,1,2] == 140)
Test.@test(array_4[2,1,3] == 150)
Test.@test(array_4[2,1,4] == 160)
Test.@test(array_4[2,2,1] == 170)
Test.@test(array_4[2,2,2] == 180)
Test.@test(array_4[2,2,3] == 190)
Test.@test(array_4[2,2,4] == 200)
Test.@test(array_4[2,3,1] == 210)
Test.@test(array_4[2,3,2] == 220)
Test.@test(array_4[2,3,3] == 230)
Test.@test(array_4[2,3,4] == 240)


text:
using BenchmarkTools, SegmentIntersections

random_coordinates = 100 .* rand(500, 4);
segments = Segment{Float64}[]
for datap in eachrow(random_coordinates)
    push!(segments, Segment(datap...))
end

println("Running brute force...")
@btime find_intersections_brute(segments);

println("Running Bentley-Ottman...")
@btime find_intersections(segments);


text:
abstract Initializer

type NormalInit <: Initializer
    mu :: AbstractFloat
    sd :: AbstractFloat
    function NormalInit(mu=0.0, sd=1.0)
        new(mu, sd)
    end
end

type UniformInit <: Initializer
    min :: AbstractFloat
    max :: AbstractFloat
    function UniformInit(min=0.0, max=1.0)
        new(min, max)
    end
end

function initialize(T::DataType, mapdim::Tuple, init::UniformInit)
    data = rand(T, mapdim)
    data -= T(init.min)
    data /= T(init.max-init.min)
end

function initialize(T::DataType, mapdim::Tuple, init::NormalInit)
    data = convert(Array{T},randn(mapdim))
    data *= T(init.sd)
    data += T(init.mu)
end


text:
# This can probably be cleaned up and simplifed with a little thought
@inline wavelength{T<:Unitful.Frequency}(s::T) = c0/s
@inline wavelength{T<:Unitful.Length}(s::T) = s
@inline wavelength(attr::OpticalAttr) = wavelength.(attr.source)
@inline wavelength(laser::AbstractLaser) = wavelength(laser.frequency)
@inline frequency{T<:Unitful.Frequency}(s::T) = s
@inline frequency{T<:Unitful.Length}(s::T) = c0/s
@inline frequency(attr::OpticalAttr) = frequency.(attr.source)
@inline frequency(laser::AbstractLaser) = frequency(laser.frequency)
@inline frequency(model::Model) = frequency(model.laser)
@inline frequency(model::ToyModel) = frequency(model.ω0)

function get_attr(attr::OpticalAttr) attr.property end

function get_label(attr::OpticalAttr) attr.label end
"""
OpticalAttr is callable. Giving a source source will return the
value of the OpticalAttr interpolated at that point
"""
(attr::OpticalAttr)(source) = attr.fit_func(frequency(source))
"""
Alias for der
"""
function der(attr::OpticalAttr, query; order=1)
    der(attr.fit_func, query; order=order)
end

function fit(attr::OpticalAttr, poly_order=12)
    ω = frequency(attr)
    μ = mean(ω)
    σ = std(ω)
    attr_unit = oneunit(first(get_attr(attr)))
    p = attr_unit*polyfit((ω-μ)/σ, get_attr(attr)./attr_unit, poly_order)
    attr.fit_func = ScaledFit(μ, σ, p)
    return attr
end