text:
using PyPlot

	# 2D approximation example
function plot2D()

	ndims = 2

	doplot = false

	# bounds
	lb = [-2.1,-2.5]
	ub = [2.2,2.6]

	# number of eval points and basis functions:
	# we require square basis matrices!
	npoints = [11,12]

	# number of basis funcs
	nbasis = npoints

	# splien degrees
	degs = [3,3]

	# implies a number of knots for each spline
	# remember the restriction that nknots == ncoefs
	nknots = Dict(i => nbasis[i] - degs[i] + 1 for i=1:ndims)

	# eval points
	points = Dict(i => collect(linspace(lb[i],ub[i],npoints[i])) for i=1:ndims)

	# set up BSplines
	bsp = Dict(i => BSpline(nknots[i],degs[i],lb[i],ub[i]) for i=1:ndims)

	# set of basis functions
	d = Dict{Integer,Array{Float64,2}}()
	for i=1:ndims
		d[i] = full(getBasis(points[i],bsp[i]))
	end

	# set of INVERSE basis functions
	id = Dict{Integer,Array{Float64,2}}()
	for k in collect(keys(d))
		id[k] = inv(d[k])
	end

	#  get a function
	function f(x,y)
		sin(sqrt(x^2+y^2))
	end

	# compute function values so that
	# k is fastest varying index
	y = Float64[f(i,j) for i in points[1], j in points[2]]

	yvec = y[:]

	# get coefs using the function
	mycoef = getTensorCoef(id,yvec)

	 # testing tolerance
	 tol = 1e-2

	 # make a plot
	 # ===========

	new_npoints = [17,16]
	# new_npoints = npoints
	new_points = Dict(i => linspace(lb[i],ub[i],new_npoints[i]) for i=1:ndims)
	# set of new basis functions
	nd = Dict{Integer,Array{Float64,2}}()
	for i=1:ndims
		nd[i] = full(getBasis(new_points[i],bsp[i]))
	end

	pred = BSplines.evalTensor2(nd[2],nd[1],mycoef)

	BSplines.figure()
	BSplines.subplot(121,projection="3d")
	BSplines.mesh(points[2],points[1],y)
	BSplines.title("true values on grid")

	BSplines.subplot(122,projection="3d")
	pp = reshape(pred,new_npoints[1],new_npoints[2])
	BSplines.mesh(new_points[2],new_points[1],pp)
	BSplines.title("Approximation off grid")

	BSplines.suptitle("2D TEST")
end

function plot3D()
	# 3D approximation example

	ndims = 3
	doplot = false

	# bounds
	lb = [-2.1,-2.5,0.9]
	ub = [2.2,2.6,5.0]

	# number of eval points and basis functions:
	# we require square basis matrices!
	npoints = [12,16,13]

	# number of basis funcs
	nbasis = npoints

	# splien degrees
	degs = [3,3,3]

	# implies a number of knots for each spline
	# remember the restriction that nknots == ncoefs
	nknots = Dict(i => nbasis[i] - degs[i] + 1 for i=1:ndims)

	# eval points
	points = Dict(i => linspace(lb[i],ub[i],npoints[i]) for i=1:ndims)

	# set up BSplines
	bsp = Dict(i => BSpline(nknots[i],degs[i],lb[i],ub[i]) for i=1:ndims)

	# set of basis functions
	d = Dict{Integer,Array{Float64,2}}()
	for i=1:ndims
		d[i] = full(getBasis(points[i],bsp[i]))
	end

	# set of INVERSE basis functions
	id = Dict{Integer,Array{Float64,2}}()
	for k in collect(keys(d))
		id[k] = inv(d[k])
	end

	#  get a function
	function f(x,y,z)
		sin(sqrt(x^2+y^2)) - z^2
	end

	# compute function values so that
	# k is fastest varying index
	y = Float64[f(i,j,k) for i in points[1], j in points[2], k in points[3]]

	yvec = y[:]

	# get coefs using the function
	mycoef = getTensorCoef(id,yvec)

		 # make a plot
		 # ===========

		new_npoints = [17,16,5]
		# new_npoints = npoints
		new_points = Dict(i => linspace(lb[i],ub[i],new_npoints[i]) for i=1:ndims)
		# set of new basis functions
		nd = Dict{Integer,Array{Float64,2}}()
		for i=1:ndims
			nd[i] = full(getBasis(new_points[i],bsp[i]))
		end

		pred = BSplines.evalTensor3(nd[3],nd[2],nd[1],mycoef)

		BSplines.figure()
		BSplines.subplot(221,projection="3d")
		BSplines.mesh(points[2],points[1],y[:,:,1])
		BSplines.title("true values at lowest 3D state")

		BSplines.subplot(222,projection="3d")
		pp = reshape(pred,new_npoints[1],new_npoints[2],new_npoints[3])
		BSplines.mesh(new_points[2],new_points[1],pp[:,:,1])
		BSplines.title("Approximation")

		BSplines.subplot(223,projection="3d")
		BSplines.mesh(points[2],points[1],y[:,:,npoints[3]])
		BSplines.title("true values at highest 3D state")

		BSplines.subplot(224,projection="3d")
		pp = reshape(pred,new_npoints[1],new_npoints[2],new_npoints[3])
		BSplines.mesh(new_points[2],new_points[1],pp[:,:,new_npoints[3]])
		BSplines.title("Approximation")
		BSplines.suptitle("3D TEST")

end

function plot4D()
	ndims = 4

	# bounds
	lb = [-1.1,-1.5,-0.9,-1.0]
	ub = [1.2,1.6,0.9,1]

	# number of eval points and basis functions:
	# we require square basis matrices!
	npoints = [13,13,13,13]

	# number of basis funcs
	nbasis = npoints

	# splien degrees
	degs = [3,3,3,3]

	# implies a number of knots for each spline
	# remember the restriction that nknots == ncoefs
	nknots = Dict(i => nbasis[i] - degs[i] + 1 for i=1:ndims)

	# eval points
	points = Dict(i => linspace(lb[i],ub[i],npoints[i]) for i=1:ndims)

	# set up BSplines
	bsp = Dict(i => BSpline(nknots[i],degs[i],lb[i],ub[i]) for i=1:ndims)

	# set of basis functions
	d = Dict{Integer,Array{Float64,2}}()
	for i=1:ndims
		d[i] = full(getBasis(points[i],bsp[i]))
	end

	# set of INVERSE basis functions
	id = Dict{Integer,Array{Float64,2}}()
	for k in collect(keys(d))
		id[k] = inv(d[k])
	end

	#  get a function
	function f(x,y,z,w)
		sin(sqrt(x^2+y^2)) + (z-w)^3
	end

	# compute function values so that
	# k is fastest varying index
	y = Float64[f(i,j,k,w) for i in points[1], j in points[2], k in points[3], w in points[4]]

	yvec = y[:]

	# get coefs using the function
	mycoef = getTensorCoef(id,yvec)

	 # testing tolerance
	 tol = 5e-3

		 # make a plot
		 # ===========

		new_npoints = [17,16,5,7]
		# new_npoints = npoints
		new_points = Dict(i => linspace(lb[i],ub[i],new_npoints[i]) for i=1:ndims)
		# set of new basis functions
		nd = Dict{Integer,Array{Float64,2}}()
		for i=1:ndims
			nd[i] = full(getBasis(new_points[i],bsp[i]))
		end

		# predict everywhere
		pred = BSplines.evalTensor4(nd[4],nd[3],nd[2],nd[1],mycoef)

		BSplines.figure()
		BSplines.subplot(221,projection="3d")
		BSplines.mesh(points[2],points[1],y[:,:,1,1])
		BSplines.title("true values at lowest state (:,:,1,1)")

		BSplines.subplot(222,projection="3d")
		pp = reshape(pred,new_npoints[1],new_npoints[2],new_npoints[3],new_npoints[4])
		BSplines.mesh(new_points[2],new_points[1],pp[:,:,1,1])
		BSplines.title("Approximation")

		BSplines.subplot(223,projection="3d")
		BSplines.mesh(points[2],points[1],y[:,:,npoints[3],npoints[4]])
		BSplines.title("true values at state (:,:,end,end)")

		BSplines.subplot(224,projection="3d")
		pp = reshape(pred,new_npoints[1],new_npoints[2],new_npoints[3],new_npoints[4])
		BSplines.mesh(new_points[2],new_points[1],pp[:,:,new_npoints[3],new_npoints[4]])
		BSplines.title("Approximation")
		BSplines.suptitle("4D TEST")
end


text:
"""
Records the test procedure I performed for data of stage II.
Notes,
- it has been demonstrated that DRP, and EBV are not good for training.
- The source codes for above can be found in the data directory.
- This package starts from using deregressed EBV (GSE, 2009, 41:55).
"""
function workflow()
    sum_stage_1()
end

ext: jl


text:
# Author: Lukasz Bratos

# Funkcja f wyliczajaca wartosc dla danego x
function f(x :: Float64)
    return sqrt(x^2 + one(Float64)) - one(Float64)
end

# Funkcja g wyliczajaca wartosc dla danego x
function g(x :: Float64)
    return x^2 / (sqrt(x^2 + one(Float64)) + one(Float64))
end

# Wypisywanie wartości w formacie pod LaTeXa
for i in 1:10
    println(i, " & ", f(8.0 ^ (-i)), " & ", g(8.0 ^ (-i)), " \\\\\n\\hline")
end

for i in 1:10
    println(20 * i, " & ", f(8.0 ^ (-20 * i)), " & ", g(8.0 ^ (-20 * i)), " \\\\\n\\hline")
end


text:
const olm_pomc_c = 1e5 # exploration parameter for MCTS
const olm_pomc_iterations = 300 # exploration parameter for MCTS
const olm_pomc_dt = 0.5 # exploration parameter for MCTS

mutable struct PomcNode
  x::AbstractArray{Float64, 1}
  u::AbstractArray{Float64, 1}
  b::AbstractArray{Float64, 1}
  adv_x::AbstractArray{Float64, 1}
  c::pomdp.Collision
  q::Float64
  N::Int
  Ns::Int

  PomcNode() = new(Float64[], Float64[], Float64[], Float64[],
                   pomdp.Collision(), 0.0, 1, 1)
end

function plan_pomc(x::AbstractArray{Float64, 1}, b::AbstractArray{Float64, 1},
                   agent::Agent, adv_agent::Agent, world::World,
                   reward::Function, ctrl_d::Discretization, depth::Int)
  #backup variables
  adv_agent_controller! = adv_agent.controller!
  adv_agent_custom = adv_agent.custom

  x = copy(x)
  x[1] = mod(x[1], world.road.path.S[end])

  value = PomcNode()
  value.x = x
  value.b = b
  value.adv_x = adv_agent.x
  adv_agent.controller! = pomdp.adv_controller!
  (c, _) = adv.predict_collision(adv_agent.x, agent.x, world)
  value.c = c
  node = GroupTree(value)

  num_iterations = round(Int, olm_pomc_iterations / 4 * depth)

  visited = Set{PomcNode}()
  for i in 1:num_iterations
    simulate_pomc(node, agent, adv_agent, world, reward,
                  ctrl_d, visited, depth)
  end

  us = Float64[]
  qs = -Inf
  ret = nothing
  (adv_c, _) = predict_collision(adv_agent.x, x, world)
  P = map(o -> reduce(+, map(i -> node.value.b[i] *
                             pomdp.P_adv(o, adv_c, pomdp.DRIVERS[i]),
                             1:length(pomdp.DRIVERS))), pomdp.ACTIONS)
  for a in 1:length(node.next)
    q = reduce(+, map(o -> node.next[a][o].value.q * P[o],
                      1:length(pomdp.ACTIONS)))
    if q > qs
      qs = q
      us = node.next[a][1].value.u
      ret = node
    end
  end

  # revert variables
  adv_agent.controller! = adv_agent_controller!
  adv_agent.custom = adv_agent_custom

  return (us, node)
end

function simulate_pomc(node::GroupTree, agent::Agent, adv_agent::Agent,
                       world::World, reward::Function, ctrl_d::Discretization,
                       visited::Set{PomcNode}, depth::Int)
  if depth <= 0
    return 0.0
  end

  node.value.x[1] = mod(node.value.x[1], world.road.path.S[end])
  (adv_c, _) = predict_collision(node.value.adv_x, node.value.x, world)

  if !(node.value in visited)
    push!(visited, node.value)

    # predict adv agent doing three possible actions
    adv_NX = Array{Array{Float64, 1}, 1}(undef, length(pomdp.ACTIONS))
    for o in pomdp.ACTIONS
      oidx = Int(o)

      adv_nx = copy(node.value.adv_x)
      adv_agent.custom[3] = o
      sim.advance!(sim.default_dynamics!, adv_nx, Pair(adv_agent, world), 0.0,
                   olm_pomc_dt, olm_h)
      adv_nx[1] = mod(adv_nx[1], world.road.path.S[end])

      adv_NX[oidx] = adv_nx
    end

    # allocate next level of the tree
    la_len = ctrl_d.thr[end] * ctrl_d.pt[end] # number of actions to survey
    node.next = Array{Array{GroupTree, 1}, 1}(undef, la_len)
    for aidx in 1:la_len
      node.next[aidx] = Array{GroupTree, 1}(undef, length(pomdp.ACTIONS))
    end

    for la in 0:(la_len - 1)
      u = dis.ls2x(ctrl_d, la)
      agent.custom = u
      nx = copy(node.value.x)
      sim.advance!(sim.default_dynamics!, nx, Pair(agent, world), 0.0,
                   olm_pomc_dt, olm_h)
      nx[1] = mod(nx[1], world.road.path.S[end])

      ra = 0.0
      for o in pomdp.ACTIONS
        oidx = Int(o)
        adv_nx = adv_NX[oidx]
        (nc, _) = adv.predict_collision(nx, adv_nx, world)

        custom = agent.custom
        agent.custom = nc
        r = reward(node.value.x, u, nx, agent, world)
        agent.custom = custom

        bp = pomdp.update_belief(node.value.b, o, adv_c)

        value = PomcNode()
        value.x = nx
        value.adv_x = adv_nx
        value.u = u
        value.b = bp
        value.c = nc
        value.q = r

        next_node = GroupTree(value)
        node.next[la + 1][oidx] = next_node
      end
    end
    return 0.0
  end

  # probability of each observation
  P = map(o -> reduce(+, map(i -> node.value.b[i] *
                             pomdp.P_adv(o, adv_c, pomdp.DRIVERS[i]),
                             1:length(pomdp.DRIVERS))), pomdp.ACTIONS)
  # find the best action for value and exploration
  as = -1
  vs = -Inf
  for a in 1:length(node.next)
    # expected reward
    ra = reduce(+, map(o -> node.next[a][o].value.q * P[o],
                       1:length(pomdp.ACTIONS)))
    v = ra + olm_pomc_c * sqrt(log(node.value.Ns) /
                               reduce(+, map(o -> node.next[a][o].value.N,
                                             1:length(pomdp.ACTIONS))))
    if v > vs
      vs = v
      as = a
    end
  end

  # sample a transition
  cP = cumsum(P)
  rnd = rand()
  oidx = findfirst(i -> cP[i] >= rnd, 1:length(cP))

  x = node.value.x
  adv_nx = node.next[as][oidx].value.adv_x
  u = node.next[as][oidx].value.u
  nx = node.next[as][oidx].value.x
  bp = node.next[as][oidx].value.b

  (nc, _) = adv.predict_collision(nx, adv_nx, world)
  custom = agent.custom
  agent.custom = nc
  r = reward(node.value.x, u, nx, agent, world)
  agent.custom = custom

  q = r + olm_gamma * simulate_pomc(node.next[as][oidx], agent, adv_agent,
                                    world, reward, ctrl_d, visited, depth - 1)

  # update Qsa for fuck's sake ---------------------------------------------- #
  dq = (q - node.next[as][oidx].value.q)
  node.next[as][oidx].value.q += dq / node.next[as][oidx].value.N

  # this works best BY FAR
  #node.next[as].value.q = max(q, node.next[as].value.q)

  #alpha = 0.5
  #Qsa = node.next[as].value.q
  #Qsa = Qsa * (1 - alpha) + q * alpha
  #node.next[as].value.q = Qsa
  # ------------------------------------------------------------------------- #

  node.next[as][oidx].value.N += 1
  node.value.Ns += 1

  return q
end

function rollout(x::AbstractArray{Float64, 1},
                 adv_x::AbstractArray{Float64, 1},
                 b::AbstractArray{Float64, 1},
                 agent::Agent, adv_agent::Agent,
                 world::World, reward::Function,
                 ctrl_d::Discretization, depth::Int)
  if depth <= 0
    return 0.0
  end

  return 0.0

  # choose middle action, generally reasonable
end

function controller_pomc!(u::AbstractArray{Float64},
                          x::AbstractArray{Float64},
                          dx::AbstractArray{Float64},
                          agent_world::Pair{Agent, World}, t::Float64)
  agent = agent_world.first
  u[1] = agent.custom[1]
  u[2] = agent.custom[2]
end


text:
using LEM1D
using TimerOutputs, LoopVectorization

function main(;
    Δx = 30.0, Δt = 100, U0=0.0001, βz= 2e-6, sea_level_file = Bintanja,
    path_sea_level = "/home/albert/Dropbox/Riverini/Pluto/",
    solver_opt = TerraceOnly,
    slope = -1/50
    )

    uplift_type = "constant" # "constant" or "armel"
    solver = ModelOpts{solver_opt}
        # Options:
        #  - TerraceOnly
        #  - TerraceKnickpoint
        #  - TerraceDiffusion
        #  - TerraceKnickpointDiffusion
    sea_level_file  = sea_level_file
        # no need to use  "" anymore. Options:
        #  - Spratt2016-450
        #  - Spratt2016-800
        #  - Bintanja
    sealevelopts = SeaLevelOpts{sea_level_file}(path_sea_level)

    # INITIAL GEOMETRY ========================================================
    slope = slope
    # slope = -tand(5)#1/50
    elementype = "linear"
    x, z0, _ ,nn, nel, e2n = mesher(Δx,slope,elementype)
    nnodel = elementype == "linear" ? 2 : 3
    z1 = copy(z0)
    River = Profile{RiverProfile}(x, z0, e2n, nn, Δx)
    Terrace = Profile{TerraceProfile}(x, z1, e2n, nn, Δx)
    # =========================================================================

    # FEM =====================================================================
    Scratch_FEM = sparsitymatrix(e2n, nel, nnodel)
    ScratchTerrace_FEM = ScratchTerraceFEM{Float64,Int64}(nn, Scratch_FEM.K)
    # =========================================================================

    # PHYSICAL PARAMETERS - Subaerial diffusion ===============================
    # --- Rivers
    T = vanAlbada
    RiverPhysics = KnickPointsParameters{Float64, T}(
        Kr = 1e-6, # whipple and tucker 1999_ dimensional coeff of erosion
        h = 1.92, # whipple and tucker 1999
        kappa_a = 4.6071, # whipple and tucker 1999_ area-length coeff
        m = 0.5,
        n = 1.0,
        φ = FluxLimiter{T},
    )

    RiverArrays = KnickPointsArrays{typeof(x)}(River, nn, RiverPhysics)
    # --- Terrace
    D = 4.4e-4 # Fernandes and Dietrich, 1997 (Hillslope evolution...)
    # =========================================================================

    # PHYSICAL PARAMETERS - below sea level (from Malatesta for wave erosion) =
    TerracePhysics = TerraceParameters{Float64}(
        βz = βz,
        h_wb = 15.0,
        P0 = 5e-5, # shallowest
        P_off = 5e-2, # offshore
    )
    # =========================================================================

    # PHYSICAL PARAMETERS - diffusion equation ================================
    HillSlopePhysics = HillSlopeParameters{Float64}(
        κa = 5e-5, # check miguels paper
        diffusion = 0.1, # subaerial hill diffusion (K, m2/yr)
        alpha_dif = 1.0, # precipitation rate (m/yr)
    )
    # =========================================================================

    # PHYSICAL PARAMETERS - submarine diffusion (from Miguelito) ==============
    SubmarinePhysics = SubmarineParameters{Float64}(
        K_s = 1e2, # submarine diffusion coefficient
        λ = 5e-4, # submarine diffusion decay coefficient
    )
    # =========================================================================

    # SEA LEVEL ===============================================================
    starting_time  = sea_level_starting_time(sealevelopts)

    # SEA LEVEL ===============================================================
    sea_lvl_curve, sea_age  = get_sea_level(sealevelopts)
    sea_lvl_curve, sea_age, fsea = sea_level_corrections(sea_lvl_curve, sea_age, starting_time)
    # =========================================================================

    # SOLVER ==================================================================
    Δt = Δt
    break_time = maximum(sea_age.-Δt)
    t = 0.0    # initialise time

    # SOME PREALLOCATIONS
    nit = Int64(floor(break_time / Δt))
    to = TimerOutput()
    t1 = time()
    buffer = similar(z0)
    id_shore = (
        terrace = (i1 = Terrace.nnod÷2, i2 = 1),
        river = (i1 = River.nnod÷2, i2 = 1)
    )

    terrace_age = zeros(nn)

    @timeit to "main loop" for it = 1:nit

        # GET SEA LVL AND UPLIFT RATE =========================================
        @timeit to "interpolate sea level" h_sea = fsea(t)

        # -- find uplift rate at current time
        # @timeit to "interpolate uplift" U = (uplift_type == "constant" ? U0 : fU(t))::Float64
        @timeit to "interpolate uplift" U = U0

        # GET NODE BELOW SEA LEVEL===========================================
        @timeit to "id_shore" id_shore = find_shore_id(h_sea, River, Terrace, TerracePhysics, id_shore)

        # SOLVE EQUATIONS =====================================================
        @timeit to "solver" Terrace, River, to = solve(solver,
                TerracePhysics, RiverPhysics, RiverArrays, HillSlopePhysics, SubmarinePhysics,
                Δt, River, Terrace, h_sea, id_shore, buffer,
                ScratchTerrace_FEM, Scratch_FEM, to)

        # BACKGROUND UPLIFT ===================================================
        @timeit to "add uplift" uplift!(Terrace.z, River.z, U*Δt)

        # ERODED VOLUME =======================================================
        # push!(erosion, erodedvolume(sol0,sol,x))
        # erosionOut[it] = erodedvolume(sol0,River.z,x)

        t +=  Δt

        # TERRACE AGE =========================================================
        @turbo for i in id_shore.terrace.i1:length(x)
            terrace_age[i] = t
        end

        if t > break_time
            break
        end

    end

    t2 = time()
    ftime = t2-t1
    println(to)
    println("\n Finished after $ftime seconds ")

    return Terrace, River, terrace_age

end

function runall(δt, δx)
    # to  = TimerOutput()
    for dt in float(50:50:500), dx in δx

        Terrace,  = main(Δx = dx, Δt = dt)
        fname = string("output_dt_old/Potato_dt_",dt,"_dx_", dx,".h5")

        h5open(fname, "w") do file
                g           = create_group(file, "Terrace")       # create a group
                p           = create_group(file, "Parameters")  # create a group
                g["x"]      = Terrace.x                     # create a scalar dataset inside the group
                g["z"]      = Terrace.z                     # create a scalar dataset inside the group
                p["dt"]     = dt
                p["dx"]     = dx
        end
    end
    # to
end


text:
@doc raw"""
# Planar Point Vortices with linear one-form

"""
module PointVorticesLinear

    using GeometricEquations

    export point_vortices_ode, point_vortices_iode, point_vortices_dg,
           point_vortices_formal_lagrangian,
           hamiltonian, angular_momentum, ϑ1, ϑ2, ϑ3, ϑ4,
           compute_energy, compute_energy_error, compute_angular_momentum_error,
           compute_momentum_error, compute_one_form

    const γ₁ = 4.0
    const γ₂ = 2.0
    const d  = 1.0
    const q₀ = [γ₂*d/(γ₁+γ₂), 0.0, -γ₁*d/(γ₁+γ₂), 0.0]

    function hamiltonian(t,q)
        γ₁ * γ₂ * log( (q[1] - q[3])^2 + (q[2] - q[4])^2 ) / (4π)
    end

    function angular_momentum(t,q)
        # γ₁ * (q[1]^2 + q[2]^2) * S(q[1],q[2]) +
        # γ₂ * (q[3]^2 + q[4]^2) * S(q[3],q[4])
        q[1] * ϑ2(t, q) - q[2] * ϑ1(t, q) +
        q[3] * ϑ4(t, q) - q[4] * ϑ3(t, q)
    end

    function ϑ1(t, q)
        - γ₁ * q[2] / 2
    end

    function ϑ2(t, q)
        + γ₁ * q[1] / 2
    end

    function ϑ3(t, q)
        - γ₂ * q[4] / 2
    end

    function ϑ4(t, q)
        + γ₂ * q[3] / 2
    end

    function dϑ1d1(t, q)
        zero(eltype(q))
    end

    function dϑ1d2(t, q)
        - γ₁ / 2
    end

    function dϑ1d3(t, q)
        zero(eltype(q))
    end

    function dϑ1d4(t, q)
        zero(eltype(q))
    end

    function dϑ2d1(t, q)
        + γ₁ / 2
    end

    function dϑ2d2(t, q)
        zero(eltype(q))
    end

    function dϑ2d3(t, q)
        zero(eltype(q))
    end

    function dϑ2d4(t, q)
        zero(eltype(q))
    end

    function dϑ3d1(t, q)
        zero(eltype(q))
    end

    function dϑ3d2(t, q)
        zero(eltype(q))
    end

    function dϑ3d3(t, q)
        zero(eltype(q))
    end

    function dϑ3d4(t, q)
        - γ₂ / 2
    end

    function dϑ4d1(t, q)
        zero(eltype(q))
    end

    function dϑ4d2(t, q)
        zero(eltype(q))
    end

    function dϑ4d3(t, q)
        + γ₂ / 2
    end

    function dϑ4d4(t, q)
        zero(eltype(q))
    end

    function ϑ(t, q, p)
        p[1] = ϑ1(t, q)
        p[2] = ϑ2(t, q)
        p[3] = ϑ3(t, q)
        p[4] = ϑ4(t, q)
    end

    function ω(t, q, Ω)
        Ω[1,1] = 0
        Ω[1,2] = dϑ1d2(t,q) - dϑ2d1(t,q)
        Ω[1,3] = dϑ1d3(t,q) - dϑ3d1(t,q)
        Ω[1,4] = dϑ1d4(t,q) - dϑ4d1(t,q)

        Ω[2,1] = dϑ2d1(t,q) - dϑ1d2(t,q)
        Ω[2,2] = 0
        Ω[2,3] = dϑ2d3(t,q) - dϑ3d2(t,q)
        Ω[2,4] = dϑ2d4(t,q) - dϑ4d2(t,q)

        Ω[3,1] = dϑ3d1(t,q) - dϑ1d3(t,q)
        Ω[3,2] = dϑ3d2(t,q) - dϑ2d3(t,q)
        Ω[3,3] = 0
        Ω[3,4] = dϑ3d4(t,q) - dϑ4d3(t,q)

        Ω[4,1] = dϑ4d1(t,q) - dϑ1d4(t,q)
        Ω[4,2] = dϑ4d2(t,q) - dϑ2d4(t,q)
        Ω[4,3] = dϑ4d3(t,q) - dϑ3d4(t,q)
        Ω[4,4] = 0

        nothing
    end

    function f1(t, q, v)
        + γ₁ * v[2] / 2
    end

    function f2(t, q, v)
        - γ₁ * v[1] / 2
    end

    function f3(t, q, v)
        + γ₂ * v[4] / 2
    end

    function f4(t, q, v)
        - γ₂ * v[3] / 2
    end

    function dHd1(t, q)
        + γ₁ * γ₂ * (q[1] - q[3]) / ( (q[1] - q[3])^2 + (q[2] - q[4])^2 ) / (2π)
    end

    function dHd2(t, q)
        + γ₁ * γ₂ * (q[2] - q[4]) / ( (q[1] - q[3])^2 + (q[2] - q[4])^2 ) / (2π)
    end

    function dHd3(t, q)
        - γ₁ * γ₂ * (q[1] - q[3]) / ( (q[1] - q[3])^2 + (q[2] - q[4])^2 ) / (2π)
    end

    function dHd4(t, q)
        - γ₁ * γ₂ * (q[2] - q[4]) / ( (q[1] - q[3])^2 + (q[2] - q[4])^2 ) / (2π)
    end

    function dH(t, q, dH)
        dH[1] = dHd1(t, q)
        dH[2] = dHd2(t, q)
        dH[3] = dHd3(t, q)
        dH[4] = dHd4(t, q)
        nothing
    end

    function point_vortices_p₀(q₀, t₀=0)
        p₀ = zero(q₀)
        tq = zeros(eltype(q₀), size(q₀,1))
        tp = zeros(eltype(p₀), size(p₀,1))

        if ndims(q₀) == 1
            ϑ(t₀, q₀, p₀)
        else
            for i in 1:size(q₀,2)
                simd_copy_xy_first!(tq, q₀, i)
                ϑ(t₀, tq, tp)
                simd_copy_yx_first!(tp, p₀, i)
            end
        end
        p₀
    end

    function point_vortices_ode_v(t, q, v)
        v[1] = - dHd2(t,q) / γ₁
        v[2] = + dHd1(t,q) / γ₁
        v[3] = - dHd4(t,q) / γ₂
        v[4] = + dHd3(t,q) / γ₂
        nothing
    end

    function point_vortices_ode(q₀=q₀)
        ODE(point_vortices_ode_v, q₀)
    end

    function point_vortices_iode_ϑ(t, q, v, p)
        p[1] = ϑ1(t,q)
        p[2] = ϑ2(t,q)
        p[3] = ϑ3(t,q)
        p[4] = ϑ4(t,q)
        nothing
    end

    function point_vortices_iode_f(t, q, v, f)
        f[1] = f1(t,q,v) - dHd1(t,q)
        f[2] = f2(t,q,v) - dHd2(t,q)
        f[3] = f3(t,q,v) - dHd3(t,q)
        f[4] = f4(t,q,v) - dHd4(t,q)
        nothing
    end

    function point_vortices_iode_g(t, q, λ, g)
        g[1] = f1(t,q,λ)
        g[2] = f2(t,q,λ)
        g[3] = f3(t,q,λ)
        g[4] = f4(t,q,λ)
        nothing
    end

    function point_vortices_iode_v(t, q, p, v)
        point_vortices_ode_v(t, q, v)
    end

    function point_vortices_iode(q₀=q₀)
        p₀ = point_vortices_p₀(q₀)
        IODE(point_vortices_iode_ϑ, point_vortices_iode_f,
             point_vortices_iode_g, q₀, p₀;
             v̄=point_vortices_iode_v)
    end

    function point_vortices_dg(q₀=q₀)
        IODE(point_vortices_iode_ϑ, point_vortices_iode_f,
             point_vortices_iode_g, q₀, q₀;
             v̄=point_vortices_iode_v)
    end

    function point_vortices_formal_lagrangian(q₀=q₀)
        p₀ = point_vortices_p₀(q₀)
        LODE(ϑ, point_vortices_iode_f, point_vortices_iode_g, q₀, p₀;
             v̄=point_vortices_iode_v, Ω=ω, ∇H=dH)
    end

    function compute_energy(t, q)
        h = zeros(q.nt+1)
        for i in 1:(q.nt+1)
            h[i] = hamiltonian(t.t[i], q.d[:,i])
        end
        return h
    end

    function compute_energy_error(t, q)
        h = zeros(q.nt+1)
        for i in 1:(q.nt+1)
            h[i] = hamiltonian(t.t[i], q.d[:,i])
        end
        h_error = (h .- h[1]) / h[1]
    end

    function compute_angular_momentum_error(t, q)
        P = zeros(q.nt+1)
        for i in 1:(q.nt+1)
            P[i] = angular_momentum(t.t[i], q.d[:,i])
        end
        P_error = (P .- P[1]) / P[1]
    end

    function compute_momentum_error(t, q, p)
        p1_error = zeros(q.nt+1)
        p2_error = zeros(q.nt+1)
        p3_error = zeros(q.nt+1)
        p4_error = zeros(q.nt+1)

        for i in 1:(q.nt+1)
            p1_error[i] = p.d[1,i] - ϑ1(t.t[i], q.d[:,i])
            p2_error[i] = p.d[2,i] - ϑ2(t.t[i], q.d[:,i])
            p3_error[i] = p.d[3,i] - ϑ3(t.t[i], q.d[:,i])
            p4_error[i] = p.d[4,i] - ϑ4(t.t[i], q.d[:,i])
        end

        (p1_error, p2_error, p3_error, p4_error)
    end

    function compute_one_form(t, q)
        p1 = zeros(q.nt+1)
        p2 = zeros(q.nt+1)
        p3 = zeros(q.nt+1)
        p4 = zeros(q.nt+1)

        for i in 1:(q.nt+1)
            p1[i] = ϑ1(t.t[i], q.d[:,i])
            p2[i] = ϑ2(t.t[i], q.d[:,i])
            p3[i] = ϑ3(t.t[i], q.d[:,i])
            p4[i] = ϑ4(t.t[i], q.d[:,i])
        end

        (p1, p2, p3, p4)
    end

end


text:
## Is In Alphabet ? ##

function in{T}(s::AbstractArray{T},y::Union(IntSet,Set,Associative))
  res = BitArray(size(s))
  for i in 1:length(s)
    res[i] = has(y,s[i])
  end
  res
end

in(s::ASCIIString,y)= in(s.data,y)

## Check Alphabet ##

function check{T}(s::AbstractArray{T},y::Union(IntSet,Set,Associative))
  for x in s
    if !has(y,x)
      error("$x not in this $(typeof(y))")
    end
  end
end

check(s::ASCIIString,y)= check(s.data,y)

## Interchange between Dic values ##

function swap!{T<:Number}(s::AbstractArray{T},y::Dict{T,T})
  for i in 1:length(s)
    s[i] = get(y,s[i],s[i])
  end
  s
end

function swap{T,Tk,Tv}(s::AbstractArray{T},y::Dict{Tk,Tv})
  res = similar(s,Tv)
  for i in 1:length(s)
    res[i] = get(y,s[i],s[i])
  end
  res
end

swap(s::ASCIIString,y)= ASCIIString(convert(Vector{Uint8},swap(s.data,y)))

## TODO: rand


text:
#=

A backend for Radiance

To visualize RAD files, check
https://www.ladybug.tools/spider-rad-viewer/rad-viewer/r7/rad-viewer.html

=#
export radiance,
       save_rad,
       @save_rad,
       radiance_material,
       radiance_plastic_material,
       radiance_metal_material,
       radiance_glass_material

write_rad_primitive(io::IO, modifier, typ, identifier, strings, ints, reals) =
  begin
    print_elems(elems) =
      begin
        print(io, length(elems))
        for e in elems print(io, " ", e) end
        println(io)
        # Some viewers are (incorrectly) restrictive regarding the RAD format
        #for i in 1:length(elems)
        #  print(io, " ", elems[i])
        #  if i%3 == 0
        #    println(io)
        #  end
        #end
      end
    println(io, modifier, " ", typ, " ", identifier)
    print_elems(strings)
    print_elems(ints)
    print_elems(reals)
    println(io)
  end

write_rad_polygon(io::IO, modifier, id, vertices) =
  begin
    with(current_backend, autocad) do
      polygon(vertices)
    end
    println(io, modifier, " ", "polygon", " ", id)
    println(io, 0) #0 strings
    println(io, 0) #0 ints
    println(io, 3*length(vertices))
    for v in vertices println(io, " ", v.x, " ", v.y, " ", v.z) end
  end

write_rad_sphere(io::IO, modifier, id, center, radius) =
  write_rad_primitive(io, modifier, "sphere", id, [], [],
    [center.x, center.y, center.z, radius])

write_rad_cone(io::IO, modifier, id, bot, bot_radius, top, top_radius) =
  write_rad_primitive(io, modifier, "cone", id, [], [],
    [bot.x, bot.y, bot.z,
     top.x, top.y, top.z,
     bot_radius, top_radius])

write_rad_cup(io::IO, modifier, id, bot, bot_radius, top, top_radius) =
  write_rad_primitive(io, modifier, "cup", id, [], [],
    [bot.x, bot.y, bot.z,
     top.x, top.y, top.z,
     bot_radius, top_radius])

write_rad_cylinder(io::IO, modifier, id, bot, radius, top) =
  write_rad_primitive(io, modifier, "cylinder", id, [], [],
    [bot.x, bot.y, bot.z,
     top.x, top.y, top.z,
     radius])

write_rad_tube(io::IO, modifier, id, bot, radius, top) =
  write_rad_primitive(io, modifier, "tube", id, [], [],
    [bot.x, bot.y, bot.z,
     top.x, top.y, top.z,
     radius])

write_rad_ring(io::IO, modifier, id, center, in_radius, out_radius, normal) =
   write_rad_primitive(io, modifier, "ring", id, [], [],
     [center.x, center.y, center.z,
      normal.x, normal.y, normal.z,
      in_radius, out_radius])

write_rad_quad(io::IO, modifier, id, sub_id, v0, v1, v2, v3) =
  begin
    println(io, modifier, " ", "polygon", " ", id, sub_id)
    println(io, 0) #0 strings
    println(io, 0) #0 ints
    println(io, 12)
    println(io, " ", v0.x, " ", v0.y, " ", v0.z)
    println(io, " ", v1.x, " ", v1.y, " ", v1.z)
    println(io, " ", v2.x, " ", v2.y, " ", v2.z)
    println(io, " ", v3.x, " ", v3.y, " ", v3.z)
  end

#=
Higher level primitives
=#

write_rad_box(io::IO, modifier, id, p0, l, w, h) =
  let p1 = p0 + vx(l),
      p2 = p0 + vxy(l, w),
      p3 = p0 + vy(w),
      p4 = p0 + vz(h),
      p5 = p4 + vx(l),
      p6 = p4 + vxy(l, w),
      p7 = p4 + vy(w)
    write_rad_quad(io, modifier, id, "face0", p0, p1, p5, p4)
    write_rad_quad(io, modifier, id, "face1", p1, p2, p6, p5)
    write_rad_quad(io, modifier, id, "face2", p2, p3, p7, p6)
    write_rad_quad(io, modifier, id, "face3", p3, p0, p4, p7)
    write_rad_quad(io, modifier, id, "face4", p3, p2, p1, p0)
    write_rad_quad(io, modifier, id, "face5", p2, p5, p6, p7)
  end

#=
Materials
=#

struct RadianceMaterial
  name::String
  type::String
  red::Real
  green::Real
  blue::Real
  specularity::Union{Real, Nothing}
  roughness::Union{Real, Nothing}
  transmissivity::Union{Real, Nothing}
  transmitted_specular::Union{Real, Nothing}
end

radiance_string(m::RadianceMaterial) =
  isnothing(m.transmissivity) && isnothing(m.transmitted_specular) ?
    isnothing(m.specularity) && isnothing(m.roughness) ?
      "void $(m.type) $(m.name)\n0\n0\n3 $(m.red) $(m.green) $(m.blue)\n" :
      "void $(m.type) $(m.name)\n0\n0\n5 $(m.red) $(m.green) $(m.blue) $(m.specularity) $(m.roughness)\n" :
    "void $(m.type) $(m.name)\n0\n0\n7 $(m.red) $(m.green) $(m.blue) $(m.specularity) $(m.roughness) $(m.transmissivity) $(m.transmitted_specular)\n"

write_rad_material(io::IO, material::RadianceMaterial) =
  write(io, radiance_string(material))

#Base.show(io::IO, mat::RadianceMaterial) =
#  write_rad_material(io, mat)

radiance_material(name::String, type::String;
                  gray::Real=0.3,
                  red::Real=gray, green::Real=gray, blue::Real=gray,
                  specularity=nothing, roughness=nothing,
                  transmissivity=nothing, transmitted_specular=nothing) =
  RadianceMaterial(name, type,
                   red, green, blue,
                   specularity, roughness,
                   transmissivity, transmitted_specular)

radiance_light_material(name::String; args...) =
  radiance_material(name, "light"; args...)

radiance_plastic_material(name::String; args...) =
  radiance_material(name, "plastic"; specularity=0, roughness=0, args...)

radiance_metal_material(name::String; args...) =
  radiance_material(name, "metal"; specularity=0, roughness=0, args...)

radiance_glass_material(name::String; args...) =
  radiance_material(name, "glass"; args...)

#Some pre-defined materials
radiance_material_white = radiance_plastic_material("white", gray=1.0) #
radiance_generic_ceiling_70 = radiance_plastic_material("GenericCeiling_70", gray=0.7)
radiance_generic_ceiling_80 = radiance_plastic_material("GenericCeiling_80", gray=0.8)
radiance_generic_ceiling_90 = radiance_plastic_material("HighReflectanceCeiling_90", gray=0.9)
radiance_generic_floor_20 = radiance_plastic_material("GenericFloor_20", gray=0.2)
radiance_generic_interior_wall_50 = radiance_plastic_material("GenericInteriorWall_50", gray=0.5)
radiance_generic_interior_wall_70 = radiance_plastic_material("GenericInteriorWall_70", gray=0.7)
radiance_generic_furniture_50 = radiance_plastic_material("GenericFurniture_50", gray=0.5)
radiance_outside_facade_30 = radiance_plastic_material("OutsideFacade_30", gray=0.3)
radiance_outside_facade_35 = radiance_plastic_material("OutsideFacade_35", gray=0.35)
radiance_generic_glass_80 = radiance_glass_material("Glass_80", gray=0.8)
radiance_generic_metal = radiance_metal_material("SheetMetal_80", gray=0.8)

export radiance_material_white,
       radiance_generic_ceiling_70,
       radiance_generic_ceiling_80,
       radiance_generic_ceiling_90,
       radiance_generic_floor_20,
       radiance_generic_interior_wall_50,
       radiance_generic_interior_wall_70,
       radiance_generic_furniture_50,
       radiance_outside_facade_30,
       radiance_outside_facade_35,
       radiance_generic_glass_80,
       radiance_generic_metal

#=
Simulations need to be done on a temporary folder, so that we can have multiple
simulations running at the same time.
=#

radiance_simulation_path() =
 let (path, io) = mktemp(mktempdir(tempdir(), prefix="Radiance_"))
   close(io)
   path
 end

radiance_cmd(cmd::AbstractString) =
  # By default Radiance in Windows is installed in C:\Program Files\Radiance\bin
  let radiance_folder = "C:/Program Files/Radiance/bin/"
    radiance_folder * cmd
  end

path_replace_suffix(path::String, suffix::String) =
  let (base, old_suffix) = splitext(path)
    base * suffix
  end

radiance_oconv(radpath) =
  let octpath = path_replace_suffix(radpath, ".oct")
    run(pipeline(`$(radiance_cmd("oconv")) $radpath`, stdout=octpath))
    octpath
  end

radiance_oconv(radpath, matpath, skypath) =
  let octpath = path_replace_suffix(radpath, ".oct")
    run(pipeline(`$(radiance_cmd("oconv")) $matpath $skypath $radpath`, stdout=octpath))
    octpath
  end

radiance_rview(path_oct, camera, target, light=(1,1,1)) =
  let p = camera,
      v = target-camera
    run(`$(radiance_cmd("rvu")) -vp $(p.x) $(p.y) $(p.z) -vd $(v.x) $(v.y) $(v.z) -av $(light[1]) $(light[2]) $(light[3]) $path_oct`)
  end

#=
# Test1
room_rad = radiance_simulation_path()

open(room_rad, "w") do out
  write_rad_material(out, radiance_light_material("bright", gray=100))
  write_rad_material(out, radiance_plastic_material("red_plastic", red=.7, green=.05, blue=.05, specularity=.5, roughness=.5))
  write_rad_sphere(out, "bright", "fixture", xyz(2, 1, 1.5), .125)
  write_rad_sphere(out, "red_plastic", "ball", xyz(.7, 1.125, .625), .125)
end
radiance_rview(radiance_oconv(room_rad), xyz(2.25, .375, 1), xyz(2.25, .375, 1) + vxyz(-.25, .125, -.125))

# Test2

open(room_rad, "w") do out
  write_rad_material(out, radiance_light_material("bright", gray=100))
  write_rad_material(out, radiance_plastic_material("red_plastic", red=.7, green=.05, blue=.05, specularity=.5, roughness=.5))
  write_rad_material(out, radiance_plastic_material("gray_paint", gray=.5))
  write_rad_sphere(out, "bright", "fixture", xyz(2, 1, 1.5), .125)
  write_rad_sphere(out, "red_plastic", "ball", xyz(.7, 1.125, .625), .125)
  write_rad_quad(out, "gray_paint", "room", "1", xyz(0,0,0), xyz(0,0,1.75), xyz(3,0,1.75), xyz(3,0,0))
  write_rad_quad(out, "gray_paint", "room", "2", xyz(0,0,0), xyz(0,2,0), xyz(0,2,1.75), xyz(0,0,1.75))
  write_rad_quad(out, "gray_paint", "room", "3", xyz(0,0,0), xyz(3,0,0), xyz(3,2,0), xyz(0,2,0))
  write_rad_quad(out, "gray_paint", "room", "4", xyz(3,2,1.75), xyz(0,2,1.75), xyz(0,2,0), xyz(3,2,0))
  write_rad_quad(out, "gray_paint", "room", "5", xyz(3,2,1.75), xyz(3,2,0), xyz(3,0,0), xyz(3,0,1.75))
  write_rad_quad(out, "gray_paint", "room", "6", xyz(3,2,1.75), xyz(3,0,1.75), xyz(0,0,1.75), xyz(0,2,1.75))
end

radiance_rview(radiance_oconv(room_rad), xyz(2.25, .375, 1), xyz(2.25, .375, 1) + vxyz(-.25, .125, -.125), (.5, .5, .5))

# Test2

open(room_rad, "w") do out
  write_rad_material(out, radiance_light_material("bright", gray=100))
  write_rad_material(out, radiance_plastic_material("red_plastic", red=.7, green=.05, blue=.05, specularity=.5, roughness=.5))
  write_rad_material(out, radiance_plastic_material("gray_paint", gray=.5))
  write_rad_sphere(out, "bright", "fixture", xyz(2, 1, 1.5), .125)
  write_rad_sphere(out, "red_plastic", "ball", xyz(.7, 1.125, .625), .125)
  write_rad_quad(out, "gray_paint", "room", "1", xyz(0,0,0), xyz(0,0,1.75), xyz(3,0,1.75), xyz(3,0,0))
  write_rad_quad(out, "gray_paint", "room", "2", xyz(0,0,0), xyz(0,2,0), xyz(0,2,1.75), xyz(0,0,1.75))
  write_rad_quad(out, "gray_paint", "room", "3", xyz(0,0,0), xyz(3,0,0), xyz(3,2,0), xyz(0,2,0))
  write_rad_quad(out, "gray_paint", "room", "4", xyz(3,2,1.75), xyz(0,2,1.75), xyz(0,2,0), xyz(3,2,0))
  write_rad_quad(out, "gray_paint", "room", "5", xyz(3,2,1.75), xyz(3,2,0), xyz(3,0,0), xyz(3,0,1.75))
  write_rad_quad(out, "gray_paint", "room", "6", xyz(3,2,1.75), xyz(3,0,1.75), xyz(0,0,1.75), xyz(0,2,1.75))
  write_rad_material(out, radiance_plastic_material("blue_plastic", red=.1, green=.1, blue=.6, specularity=.05, roughness=.1))
  write_rad_quad(out, "blue_plastic", "box", "1540", xyz(0.98,0.88,0), xyz(0.98,0.88,0.5), xyz(0.5,0.75,0.5), xyz(0.5,0.75,0))
  write_rad_quad(out, "blue_plastic", "box", "4620", xyz(0.5,0.75,0.5), xyz(0.37,1.23,0.5), xyz(0.37,1.23,0), xyz(0.5,0.75,0))
  write_rad_quad(out, "blue_plastic", "box", "2310", xyz(0.37,1.23,0), xyz(0.85,1.36,0), xyz(0.98,0.88,0), xyz(0.5,0.75,0))
  write_rad_quad(out, "blue_plastic", "box", "3267", xyz(0.85,1.36,0), xyz(0.37,1.23,0), xyz(0.37,1.23,0.5), xyz(0.85,1.36,0.5))
  write_rad_quad(out, "blue_plastic", "box", "5137", xyz(0.98,0.88,0.5), xyz(0.98,0.88,0), xyz(0.85,1.36,0), xyz(0.85,1.36,0.5))
  write_rad_quad(out, "blue_plastic", "box", "6457", xyz(0.37,1.23,0.5), xyz(0.5,0.75,0.5), xyz(0.98,0.88,0.5), xyz(0.85,1.36,0.5))
  write_rad_material(out, radiance_metal_material("chrome", gray=.8, specularity=.9, roughness=0))
  write_rad_cylinder(out, "chrome", "fixture_support", xyz(2,1,1.5), .05, xyz(2,1,1.75))
end

radiance_rview(radiance_oconv(room_rad), xyz(2.25, .375, 1), xyz(2.25, .375, 1) + vxyz(-.25, .125, -.125), (.5, .5, .5))

=#

#=

We need to discretize paths so that we can extract the vertices
We will use some sort of tolerance to deal with curved paths

=#

abstract type RadianceKey end
const RadianceId = Int
const RadianceRef = GenericRef{RadianceKey, RadianceId}
const RadianceNativeRef = NativeRef{RadianceKey, RadianceId}
const RadianceUnionRef = UnionRef{RadianceKey, RadianceId}
const RadianceSubtractionRef = SubtractionRef{RadianceKey, RadianceId}

mutable struct RadianceBackend{K,T,LOD} <: LazyBackend{K,T}
  shapes::Shapes
  buffer::LazyParameter{IOBuffer}
  count::Integer
  materials::Dict
  camera::Loc
  target::Loc
  lens::Real
end

const Radiance{LOD} = RadianceBackend{RadianceKey, RadianceId, LOD}

#=
The Radiance backend cannot realize shapes immediately, only when requested.
=#

void_ref(b::Radiance) = RadianceNativeRef(-1)

create_radiance_buffer() = IOBuffer()

const radiance = Radiance{500}(Shape[], LazyParameter(IOBuffer, create_radiance_buffer),
  0, Dict(), xyz(10,10,10), xyz(0,0,0), 35)
const radiance_lod100 = Radiance{100}(Shape[], LazyParameter(IOBuffer, create_radiance_buffer),
  0, Dict(), xyz(10,10,10), xyz(0,0,0), 35)

buffer(b::Radiance) = b.buffer()
next_id(b::Radiance) =
    begin
        b.count += 1
        b.count -1
    end
next_modifier(b::Radiance, key) =
    get!(b.materials, key, key.name)

save_rad(path::String, b::Radiance=current_backend()) =
  let buf = b.buffer()
    take!(buf)
    for s in b.shapes
      realize(b, s)
    end
    open(path, "w") do out
      write(out, String(take!(buf)))
    end
  end

macro save_rad()
  :(save_rad(splitext($(string(__source__.file)))[1]*".rad"))
end

#

set_view(camera::Loc, target::Loc, lens::Real, b::Radiance) =
  begin
    set_view(camera, target, lens, autocad)
    b.camera = camera
    b.target = target
    b.lens = lens
  end

get_view(b::Radiance) =
  b.camera, b.target, b.lens

#current_backend(radiance)

delete_all_shapes(b::Radiance) =
  begin
    delete_all_shapes(autocad)
    (empty!(b.shapes); empty!(b.materials); nothing)
  end

realize(b::Radiance, s::Sphere) =
  let id = next_id(b),
      mod = next_modifier(b, missing),
      kind = "sphere",
      buf = buffer(b)
    write_rad_sphere(buf, "mat$(kind)$(mod)", id, in_world(s.center), s.radius)
    id
  end

#=
backend(radiance)
delete_all_shapes()
sphere()
@save_rad()
=#

  #=
realize(b::ACAD, s::Torus) =
  ACADTorus(connection(b), s.center, vz(1, s.center.cs), s.re, s.ri)
realize(b::ACAD, s::Cuboid) =
  ACADIrregularPyramidFrustum(connection(b), [s.b0, s.b1, s.b2, s.b3], [s.t0, s.t1, s.t2, s.t3])
realize(b::ACAD, s::RegularPyramidFrustum) =
    ACADIrregularPyramidFrustum(connection(b),
                                regular_polygon_vertices(s.edges, s.cb, s.rb, s.angle, s.inscribed),
                                regular_polygon_vertices(s.edges, add_z(s.cb, s.h), s.rt, s.angle, s.inscribed))
realize(b::ACAD, s::RegularPyramid) =
  ACADIrregularPyramid(connection(b),
                          regular_polygon_vertices(s.edges, s.cb, s.rb, s.angle, s.inscribed),
                          add_z(s.cb, s.h))
realize(b::ACAD, s::IrregularPyramid) =
  ACADIrregularPyramid(connection(b), s.bs, s.t)
realize(b::ACAD, s::RegularPrism) =
  let ps = regular_polygon_vertices(s.edges, s.cb, s.r, s.angle, s.inscribed)
    ACADIrregularPyramidFrustum(connection(b),
                                   ps,
                                   map(p -> add_z(p, s.h), ps))
  end
realize(b::ACAD, s::IrregularPyramidFrustum) =
    ACADIrregularPyramidFrustum(connection(b), s.bs, s.ts)

realize(b::ACAD, s::IrregularPrism) =
  ACADIrregularPyramidFrustum(connection(b),
                              s.bs,
                              map(p -> (p + s.v), s.bs))
## FIXME: deal with the rotation angle
realize(b::ACAD, s::RightCuboid) =
  ACADCenteredBox(connection(b), s.cb, s.width, s.height, s.h)
=#
realize(b::Radiance, s::Box) =
  let id = next_id(b),
      mod = next_modifier(b, missing),
      kind = "box",
      buf = buffer(b)
    write_rad_box(buf, "mat$(kind)$(mod)", id, in_world(s.c), s.dx, s.dy, s.dz)
    id
  end

#=
open(room_rad, "w") do out
  write_rad_material(out, radiance_light_material("bright", gray=100))
  write_rad_material(out, radiance_plastic_material("red_plastic", red=.7, green=.05, blue=.05, specularity=.5, roughness=.5))
  write_rad_material(out, radiance_plastic_material("gray_paint", gray=.5))
  write_rad_sphere(out, "bright", "fixture", xyz(2, 1, 1.5), .125)
  write_rad_sphere(out, "red_plastic", "ball", xyz(.7, 1.125, .625), .125)
  write_rad_box(out, "gray_paint", "1", xyz(0,0,0), 3, 2, 1.75)
  write_rad_material(out, radiance_plastic_material("blue_plastic", red=.1, green=.1, blue=.6, specularity=.05, roughness=.1))
  write_rad_quad(out, "blue_plastic", "box", "1540", xyz(0.98,0.88,0), xyz(0.98,0.88,0.5), xyz(0.5,0.75,0.5), xyz(0.5,0.75,0))
  write_rad_quad(out, "blue_plastic", "box", "4620", xyz(0.5,0.75,0.5), xyz(0.37,1.23,0.5), xyz(0.37,1.23,0), xyz(0.5,0.75,0))
  write_rad_quad(out, "blue_plastic", "box", "2310", xyz(0.37,1.23,0), xyz(0.85,1.36,0), xyz(0.98,0.88,0), xyz(0.5,0.75,0))
  write_rad_quad(out, "blue_plastic", "box", "3267", xyz(0.85,1.36,0), xyz(0.37,1.23,0), xyz(0.37,1.23,0.5), xyz(0.85,1.36,0.5))
  write_rad_quad(out, "blue_plastic", "box", "5137", xyz(0.98,0.88,0.5), xyz(0.98,0.88,0), xyz(0.85,1.36,0), xyz(0.85,1.36,0.5))
  write_rad_quad(out, "blue_plastic", "box", "6457", xyz(0.37,1.23,0.5), xyz(0.5,0.75,0.5), xyz(0.98,0.88,0.5), xyz(0.85,1.36,0.5))
  write_rad_material(out, radiance_metal_material("chrome", gray=.8, specularity=.9, roughness=0))
  write_rad_cylinder(out, "chrome", "fixture_support", xyz(2,1,1.5), .05, xyz(2,1,1.75))
end

  radiance_rview(radiance_oconv(room_rad), xyz(2.25, .375, 1), xyz(2.25, .375, 1) + vxyz(-.25, .125, -.125), (.5, .5, .5))
=#

#=
realize(b::ACAD, s::Cone) =
  ACADCone(connection(b), add_z(s.cb, s.h), s.r, s.cb)
realize(b::ACAD, s::ConeFrustum) =
  ACADConeFrustum(connection(b), s.cb, s.rb, s.cb + vz(s.h, s.cb.cs), s.rt)
=#
realize(b::Radiance, s::Cylinder) =
  let bot_id = next_id(b),
      top_id = next_id(b),
      side_id = next_id(b),
      mod = next_modifier(b, missing),
      kind = "cylinder",
      buf = buffer(b),
      bot = in_world(s.cb),
      top = in_world(s.cb + vz(s.h, s.cb.cs)),
      normal = unitized(top-bot)
    write_rad_ring(buf, "mat$(kind)bot$(mod)", bot_id, bot, 0, s.r, -normal)
    write_rad_cylinder(buf, "mat$(kind)side$(mod)", side_id, bot, s.r, top)
    write_rad_ring(buf, "mat$(kind)top$(mod)", top_id, top, 0, s.r, normal)
    bot_id
  end

#=
realize(b::Radiance, s::EmptyShape) =
    EmptyRef{RadianceId}()
realize(b::Radiance, s::UniversalShape) =
    UniversalRef{RadianceId}()

realize(b::Radiance, s::Move) =
    let r = map_ref(s.shape) do r
                RadianceMove(connection(b), r, s.v)
                r
            end
        mark_deleted(s.shape)
        r
    end

realize(b::Radiance, s::Scale) =
    let r = map_ref(s.shape) do r
                RadianceScale(connection(b), r, s.p, s.s)
                r
            end
        mark_deleted(s.shape)
        r
    end

realize(b::Radiance, s::Rotate) =
    let r = map_ref(s.shape) do r
                RadianceRotate(connection(b), r, s.p, s.v, s.angle)
                r
            end
        mark_deleted(s.shape)
        r
    end

=#

# BIM

realize_pyramid_fustrum(b::Radiance, top::String, bot::String, side::String, bot_path::Path, top_path::Path, closed=true) =
  realize_pyramid_fustrum(b, top, bot, side, path_vertices(bot_path), path_vertices(top_path), closed)

# This requires three different keys, so that gets into the materials dict.
realize_pyramid_fustrum(b::Radiance, top::String, bot::String, side::String, bot_vs, top_vs, closed=true) =
  let buf = buffer(b)
    if closed
      write_rad_polygon(buf, top, next_id(b), reverse(bot_vs))
      write_rad_polygon(buf, bot, next_id(b), top_vs)
    end
    for vs in zip(bot_vs, circshift(bot_vs, 1), circshift(top_vs, 1), top_vs)
      write_rad_polygon(buf, side, next_id(b), vs)
    end
  end

realize_polygon(b::Radiance, key, kind::String, path::Path, acw=true) =
  realize_polygon(b, key, kind, path_vertices(path), acw)

realize_polygon(b::Radiance, key, kind::String, vs, acw=true) =
  let mod = next_modifier(b, key),
      buf = buffer(b),
      side = acw ? "top" : "bot"
    write_rad_polygon(buf, "mat$(kind)$(side)$(mod)", next_id(b), acw ? vs : reverse(vs))
  end

#=

Upon daysim processing, we need to compute the useful daylight illumination, which is the fraction of time that a sensor is within a given range of illumination

To do that, we start by reading the .ill file generated

ill = CSV.read(path, delim=' ', datarow=1)

The file looks like this:

8760×136 DataFrame. Omitted printing of 109 columns
│ Row  │ Column1 │ Column2 │ Column3  │ Column4 │ Column5 │ Column6 │ Column7 │ Column8 │ Column9 │ Column10 │ Column11 │ Column12 │ Column13 │ Column14 │ Column15 │ Column16 │ Column17 │ Column18 │ Column19 │ Column20 │ Column21 │ Column22 │ Column23 │ Column24 │ Column25 │ Column26 │ Column27 │
│      │ Int64⍰  │ Int64⍰  │ Float64⍰ │ Missing │ Int64⍰  │ Int64⍰  │ Int64⍰  │ Int64⍰  │ Int64⍰  │ Int64⍰   │ Int64⍰   │ Int64⍰   │ Int64⍰   │ Int64⍰   │ Int64⍰   │ Int64⍰   │ Int64⍰   │ Int64⍰   │ Int64⍰   │ Int64⍰   │ Int64⍰   │ Int64⍰   │ Int64⍰   │ Int64⍰   │ Int64⍰   │ Int64⍰   │ Int64⍰   │
├──────┼─────────┼─────────┼──────────┼─────────┼─────────┼─────────┼─────────┼─────────┼─────────┼──────────┼──────────┼──────────┼──────────┼──────────┼──────────┼──────────┼──────────┼──────────┼──────────┼──────────┼──────────┼──────────┼──────────┼──────────┼──────────┼──────────┼──────────┤
│ 1    │ 1       │ 1       │ 0.5      │ missing │ 0       │ 0       │ 0       │ 0       │ 0       │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │
│ 2    │ 1       │ 1       │ 1.5      │ missing │ 0       │ 0       │ 0       │ 0       │ 0       │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │
│ 3    │ 1       │ 1       │ 2.5      │ missing │ 0       │ 0       │ 0       │ 0       │ 0       │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │
│ 4    │ 1       │ 1       │ 3.5      │ missing │ 0       │ 0       │ 0       │ 0       │ 0       │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │
│ 5    │ 1       │ 1       │ 4.5      │ missing │ 0       │ 0       │ 0       │ 0       │ 0       │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │
│ 6    │ 1       │ 1       │ 5.5      │ missing │ 0       │ 0       │ 0       │ 0       │ 0       │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │ 0        │
⋮

Then, for each sensor (i.e., starting at column 5) we compute the number of times that its illuminance was within predefined limits and we normalize for all possible times

=#

#=
#path = "C:\\Users\\aml\\Dropbox\\AML\\Projects\\rosetta\\tests\\xyz.ill"

path = "C:\\Users\\aml\\Downloads\\Geometry_0_0_0629.ill"

ill = CSV.read(path, delim=' ', datarow=1)

udi_in_range(df, min, max) =
  let occupied_hours_per_year = nrow(df)
    colwise(col->
      round(Int, count(i -> min < i < max, col)/occupied_hours_per_year*100),
      df[2:end,5:end])
  end

res = udi_in_range(ill, 299.9, 3000.1)
maximum(res)
minimum(res)

test = DataFrame([
0 0 0 0 0 0 10;
0 0 0 0 0 0 10;
0 0 0 0 0 0 10;
0 0 0 0 0 0 10;
0 0 0 0 0 0 10;
0 0 0 0 0 0 10;
0 0 0 0 0 0 10;
0 0 0 0 0 0 10;
0 0 0 0 0 20 10;
0 0 0 0 0 20 10;
0 0 0 0 0 20 10;
0 0 0 0 0 20 10;
0 0 0 0 0 20 10;
0 0 0 0 0 20 10;
0 0 0 0 0 20 10;
0 0 0 0 0 20 10;
])

udi_in_range(test, -1, 5)
=#

const extra_sky_rad_contents = """
skyfunc glow sky_mat
0
0
4 1 1 1 0

sky_mat source sky
0
0
4 0 0 1 180

skyfunc glow ground_glow
0
0
4 1 .8 .5 0

ground_glow source ground
0
0
4 0 0 -1 180
"""

#
const CIE_Overcast_Sky_rad_contents = """
!gensky 12 21 12.00 -c -a 42.300 -o 71.100 -m 75.000 -B 100

skyfunc glow sky_mat
0
0
4 1 1 1 0

sky_mat source sky
0
0
4 0 0 1 180

skyfunc glow ground_glow
0
0
4 1 1 1  0

ground_glow source ground
0
0
4 0 0 -1 180
"""

export_CIE_Overcast_Sky(path::String) =
  let skypath = path_replace_suffix(path, "_sky.rad")
    open(skypath, "w") do out
      write(out, CIE_Overcast_Sky_rad_contents)
    end
    skypath
  end

#=
Radiance families need to know the different kinds of materials
that go on each surface.
In some cases it might be the same material, but in others, such
as slabs, outside walls, etc, we will have different materials.
=#

abstract type RadianceFamily <: Family end

struct RadianceMaterialFamily <: RadianceFamily
  material::RadianceMaterial
end

radiance_material_family(mat::RadianceMaterial) =
  RadianceMaterialFamily(mat)

struct RadianceSlabFamily <: RadianceFamily
  top_material::RadianceMaterial
  bottom_material::RadianceMaterial
  side_material::RadianceMaterial
end

radiance_slab_family(top::RadianceMaterial, bot::RadianceMaterial=top, side::RadianceMaterial=bot) =
  RadianceSlabFamily(top, bot, side)

struct RadianceWallFamily <: RadianceFamily
  right_material::RadianceMaterial
  left_material::RadianceMaterial
end

radiance_wall_family(right::RadianceMaterial, left::RadianceMaterial=right) =
  RadianceWallFamily(right, left)

backend_get_family_ref(b::Radiance, f::Family, rf::RadianceFamily) = rf

export radiance_material_family,
       radiance_slab_family,
       radiance_wall_family

set_backend_family(default_wall_family(), radiance,
  radiance_wall_family(radiance_generic_interior_wall_70))
set_backend_family(default_slab_family(), radiance,
  radiance_slab_family(radiance_generic_floor_20, radiance_generic_ceiling_80))
set_backend_family(default_beam_family(), radiance,
  radiance_material_family(radiance_generic_metal))
set_backend_family(default_column_family(), radiance,
  radiance_material_family(radiance_generic_metal))
set_backend_family(default_door_family(), radiance,
  radiance_material_family(radiance_generic_furniture_50))
set_backend_family(default_panel_family(), radiance,
  radiance_material_family(radiance_glass_material("GenericGlass", gray=0.3)))

# Radiance-specific operations

material_ground = Parameter(radiance_generic_floor_20)

#=
create_ground_plane(shapes, material=material_ground()) =
  if shapes == []
    error("No shapes selected for analysis. Use add-radiance-shape!.")
  else
    let (p0, p1) = bounding_box(union(shapes)),
        (center, ratio) = (quad_center(p0, p1, p2, p3),
                  distance(p0, p4)/distance(p0, p2));
     ratio == 0 ?
      error("Couldn't compute height. Use add-radiance-shape!.") :
      let pts = map(p -> intermediate_loc(center, p, ratio*10), [p0, p1, p2, p3]);
         create_surface_layer(pts, 0, ground_layer(), material)
        end
       end
  end
        w = max(floor_extra_factor()*distance(p0, p1), floor_extra_width())
        with(current_layer,floor_layer()) do
          box(xyz(min(p0.x, p1.x)-w, min(p0.y, p1.y)-w, p0.z-1-floor_distance()),
              xyz(max(p0.x, p1.x)+w, max(p0.y, p1.y)+w, p0.z-0-floor_distance()))
        end
      end
    end

=#

realize(b::Radiance, s::Slab) =
  let base = vz(s.level.height + s.family.coating_thickness - s.family.thickness),
      thickness = vz(s.family.coating_thickness + s.family.thickness),
      path = s.contour,
      modtop = next_modifier(b, realize(b, s.family).top_material),
      modbot = next_modifier(b, realize(b, s.family).bottom_material),
      modside = next_modifier(b, realize(b, s.family).side_material)
    for op in s.openings
      path = subtract_paths(path, op)
    end
    realize_pyramid_fustrum(
      b, modtop, modbot, modside,
      path_vertices(translate(path, base)),
      path_vertices(translate(path, base + thickness)))
  end

#=
#FIXME define the family parameters for beams
realize(b::Radiance, s::Beam) =
    ref(right_cuboid(s.p0, 0.2, 0.2, s.p1, 0))

=#
realize(b::Radiance, s::Panel) =
  let p1 = s.vertices[1],
      p2 = s.vertices[2],
      p3 = s.vertices[3],
      n = vz(s.family.thickness/2, cs_from_o_vx_vy(p1, p2-p1, p3-p1)),
      mod = next_modifier(b, realize(b, s.family).material)
    realize_pyramid_fustrum(
        b, mod, mod, mod,
        map(p -> in_world(p - n), s.vertices),
        map(p -> in_world(p + n), s.vertices))
  end

closed_path_for_height(path, h) =
  let ps = path_vertices(path)
    closed_polygonal_path([ps..., reverse(map(p -> p+vz(h), ps))...])
  end

#=
One important restriction is that Radiance only supports _planar_ polygons.
This forces us to create multiple subpaths.
=#
realize(b::Radiance, w::Wall) =
  let w_base_height = w.bottom_level.height,
      w_height = w.top_level.height - w_base_height,
      r_thickness = r_thickness(w),
      l_thickness = l_thickness(w),
      w_path = translate(w.path, vz(w_base_height)),
      w_paths = subpaths(w_path),
      r_w_paths = subpaths(offset(w_path, -r_thickness)),
      l_w_paths = subpaths(offset(w_path, l_thickness)),
      openings = [w.doors..., w.windows...],
      prevlength = 0,
      modright = next_modifier(b, realize(b, w.family).right_material),
      modleft = next_modifier(b, realize(b, w.family).left_material)
    for (w_seg_path, r_w_path, l_w_path) in zip(w_paths, r_w_paths, l_w_paths)
      let currlength = prevlength + path_length(w_seg_path),
          c_r_w_path = closed_path_for_height(r_w_path, w_height),
          c_l_w_path = closed_path_for_height(l_w_path, w_height)
        realize_pyramid_fustrum(b, modleft, modright, modright, c_l_w_path, c_r_w_path, false)
        openings = filter(openings) do op
          if prevlength <= op.loc.x < currlength ||
             prevlength <= op.loc.x + op.family.width <= currlength # contained (at least, partially)
            let op_height = op.family.height,
                op_at_start = op.loc.x <= prevlength,
                op_at_end = op.loc.x + op.family.width >= currlength,
                op_path = subpath(w_path,
                                  max(prevlength, op.loc.x),
                                  min(currlength, op.loc.x + op.family.width)),
                r_op_path = offset(op_path, -r_thickness),
                l_op_path = offset(op_path, l_thickness),
                fixed_r_op_path =
                  open_polygonal_path([path_start(op_at_start ? r_w_path : r_op_path),
                                       path_end(op_at_end ? r_w_path : r_op_path)]),
                fixed_l_op_path =
                  open_polygonal_path([path_start(op_at_start ? l_w_path : l_op_path),
                                       path_end(op_at_end ? l_w_path : l_op_path)]),
                c_r_op_path = closed_path_for_height(translate(fixed_r_op_path, vz(op.loc.y)), op_height),
                c_l_op_path = closed_path_for_height(translate(fixed_l_op_path, vz(op.loc.y)), op_height),
                idxs = closest_vertices_indexes(path_vertices(c_r_w_path), path_vertices(c_r_op_path))
              realize_pyramid_fustrum(b, modleft, modright, modright, c_r_op_path, c_l_op_path, false)
              c_r_w_path =
                closed_polygonal_path(
                  inject_polygon_vertices_at_indexes(path_vertices(c_r_w_path), path_vertices(c_r_op_path), idxs))
              c_l_w_path =
                closed_polygonal_path(
                  inject_polygon_vertices_at_indexes(path_vertices(c_l_w_path), path_vertices(c_l_op_path), idxs))
              # preserve if not totally contained
              ! (op.loc.x >= prevlength && op.loc.x + op.family.width <= currlength)
            end
          else
            true
          end
        end
        prevlength = currlength
        #realize_polygon(b, modleft, modright, modright, c_l_w_path, false)
        #realize_polygon(b, modleft, modright, modright, c_r_w_path, true)
      end
    end
    void_ref(b)
  end

realize(b::Radiance, w::Window) = nothing
realize(b::Radiance, w::Door) = nothing

backend_wall(b::Radiance, w_path, w_height, l_thickness, r_thickness, family) =
  let w_paths = subpaths(w_path),
      r_w_paths = subpaths(offset(w_path, -r_thickness)),
      l_w_paths = subpaths(offset(w_path, l_thickness)),
      modright = next_modifier(b, realize(b, s.family).material_right),
      modleft = next_modifier(b, realize(b, s.family).material_left)
    for (w_seg_path, r_w_path, l_w_path) in zip(w_paths, r_w_paths, l_w_paths)
      let c_r_w_path = closed_path_for_height(r_w_path, w_height),
          c_l_w_path = closed_path_for_height(l_w_path, w_height)
        realize_pyramid_fustrum(b, modleft, modright, modright, c_l_w_path, c_r_w_path, false)
        # This is surely wrong!
        #realize_polygon(b, family, "wall", c_l_w_path, false)
        #realize_polygon(b, family, "wall", c_r_w_path, true)
      end
    end
    void_ref(b)
  end

#=

Sensors are need on the surfaces that are intended for analysis.
They can be placed in many different ways but that are a few standard ones,
namely those that follow an regular distribution.
=#

analysis_nodes_height = Parameter(0.5)
analysis_nodes_separation_u = Parameter{Union{Missing, Real}}(missing)
analysis_nodes_separation_v = Parameter{Union{Missing, Real}}(missing)
analysis_nodes_separation = Parameter(4.0)

analysis_nodes_u_separation()::Real =
  ismissing(analysis_nodes_separation_u()) ?
    analysis_nodes_separation() :
    analysis_nodes_separation_u()

analysis_nodes_v_separation()::Real =
  ismissing(analysis_nodes_separation_v()) ?
    analysis_nodes_separation() :
    analysis_nodes_separation_v()

#=
sensors_from_surfaces(surfaces=radiance_surfaces(),
 height::Real=analysis_nodes_height(),
 sep_u::Real=analysis_nodes_u_separation(),
 sep_v::Real=analysis_nodes_v_separation()) =
  (map(surface_from -> let_values(nu(nv)(surface_nu_nv(surface_from, sep_u, sep_v))(), let nodes = map_inner_surface_division(pt -> pt && pt+vz(height), surface_from, nu, nv); filter(identity, append*(nodes)) end),
                surfaces))

=#

sensors_locations(s::Slab, b::Radiance) =
  let out_pts = path_vertices(s.countour),
      in_ptss = map(path_vertices, s.openings),
      p_ns = Main.plane_surface_sensor_locations(out_pts, in_ptss)
  end

abstract type Analysis end
abstract type StructuralAnalysis <: Analysis end
abstract type LightingAnalysis <: Analysis end

struct RadianceVisualization <: LightingAnalysis
end

export radiance_visualization
radiance_visualization(b::Radiance=radiance) =
  let path=radiance_simulation_path(),
      radpath = export_geometry(b, path),
      matpath = export_materials(b, path),
      skypath = export_CIE_Overcast_Sky(path), #export_sky(b, path),
      octpath = radiance_oconv(radpath, matpath, skypath)
    radiance_rview(octpath, b.camera, b.target) #Ignoring lens for now
  end

struct DaysimAnalysis <: LightingAnalysis
end

struct RadianceMapAnalysis <: LightingAnalysis
  sensors::Locs
  path::AbstractString
  location::AbstractString
  datetime::DateTime
  sky::AbstractString
  simulation_parameters::AbstractString
  results::Parameter{Any} # HACK: FIX THIS
end

export radiance_map_analysis
radiance_map_analysis(sensors;
                      location="Lisbon",
                      datetime=Dates.now(UTC),
                      sky=sky_for(location, datetime),
                      simulation_parameters="-ab 2 -ad 1000 -as 20 -ar 300 -aa 0.1",
                      path=radiance_simulation_path()) =
  RadianceMapAnalysis(sensors, path, location, datetime, sky, simulation_parameters, Parameter{Any}(nothing))

@defop analyze(a::Analysis)

analyze(a::RadianceMapAnalysis, b::Radiance) =
  let radpath = export_geometry(b, a.path),
      matpath = export_materials(b, a.path),
      skypath = export_sky(b, a.path),
      octpath = path_replace_suffix(a.path, ".oct"),
      ptspath = export_sensors(b, a.path, a.sensors),
      datpath = path_replace_suffix(a.path, ".dat")
    radiance_cmd("""oconv "$(matpath)" "$(skypath)" "$(radpath)" > "$(octpath)" """)
    radiance_cmd("""rtrace -I -h -dp 2048 -ms 0.063 -ds .2 -dt .05 -dc .75 -dr 3 -st .01 -lr 12 -lw .0005 $(a.simulation_parameters) "$(octpath)" < "$(ptspath)" > "$(datpath)" """)
    a.results(read_rtrace_results(a.path))
  end

export_geometry(b::Radiance, path::AbstractString) =
  let radpath = path_replace_suffix(path, ".rad"),
      buf = b.buffer()
    take!(buf)
    for s in b.shapes
      realize(b, s)
    end
    add_ground_plane(b.shapes, buf)
    open(radpath, "w") do out
      write(out, String(take!(buf)))
    end
    radpath
  end

add_ground_plane(shapes, buf) =
  @warn "Not generating ground plane"

used_materials(b::Radiance) =
  let materials = unique(map(f -> realize(s.family, b), b.shapes))
      materials
    end

export_materials(b::Radiance, path::AbstractString) =
  let matpath = path_replace_suffix(path, "_materials.rad")
    open(matpath, "w") do out
      for (k,v) in b.materials
        write_rad_material(out, k)
      end
    end
    matpath
  end

export_sky(b::Radiance, path::AbstractString;
           date::DateTime=DateTime(0, 9, 21),
           latitude::Real=61,
           longitude::Real=150,
           meridian::Real=135,
           sun::Bool=true) =
  let skypath = path_replace_suffix(path, "_sky.rad"),
      d2(i) = lpad(i, 2, '0')
    open(skypath, "w") do out
      let mo = d2(month(date)),
          da = d2(day(date)),
          ho = d2(hour(date)),
          mi = d2(minute(date))
        write(out, "!gensky $(mo) $(da) $(ho):$(mi) $(sun ? "+s" : "-s") -a $(latitude) -o $(longitude) -m $(meridian)\n")
        write(out, extra_sky_rad_contents)
      end
    end
    skypath
  end

export_sky_for_date_location(b::Radiance, path::AbstractString, date::Date, location::String) =
  let (place, latitude, longitude, time_zone, site_elevation, weapath) = get_weather_for_location(path, location)
    export_sky(path,
               date,
               date,
               latitude,
               latitude,
               longitude,
               longitude,
               meridian,
               time_zone)
  end

radiance_sensors = Parameter(Loc[])

export_sensors(path::Path, sensors=radiance_sensors()) =
  let ptspath = path_replace_suffix(path, ".pts")
    open(ptspath, "w") do out
      with(current_layer, sensor_layer()) do
        for p in sensors
          let (wp, wv) = (in_world(p), in_world(vz(1, p)))
            #;Just to see the sensor location
            point(p)
            #    (line p (+z p 3))
            write(out, "$(wp.x),$(wp.y),$(wp.z),$(wv.x),$(wv.y),$(wv.z)\n")
          end
        end
      end
    end
    ptspath
  end


text:
## Defines all hashing functionality
## As usual, check out http://www.lysator.liu.se/~nisse/nettle/nettle.html#Hash-functions

import Base: show
export HashState, HashAlgorithm, HashAlgorithms

# All Hash Algorithms derive from this abstract type
abstract HashAlgorithm

# This is our rather poorly named list of HashAlgorithm types
HashAlgorithms = DataType[]

# This is the user-facing type that is used to actually hash stuff
type HashState{T<:HashAlgorithm}
  ctx::Array{Uint8,1}
end

# This is a mirror of the nettle-meta.h:nettle_hash struct
immutable NettleHash
    name::Ptr{Uint8}
    context_size::Cuint
    digest_size::Cuint
    block_size::Cuint
    init::Ptr{Void}     # nettle_hash_init_func
    update::Ptr{Void}   # nettle_hash_update_func
    digest::Ptr{Void}   # nettle_hash_digest_func
end

# We're going to load in each nettle_hash struct individually, deriving
# HashAlgorithm types off of the names we find, and calculating the output
# and context size from the data members in the C structures
begin
  hash_idx = 1
  while( true )
    nhptr = unsafe_load(cglobal(("nettle_hashes",nettle),Ptr{Ptr{Void}}),hash_idx)
    if nhptr == C_NULL
      break
    end
    nh = unsafe_load(convert(Ptr{NettleHash}, nhptr))

    # Otherwise, we continue on to derive the information from this struct
    name = symbol(uppercase(bytestring(nh.name)))

    # Load in the Ptr{Uint32}'s
    ctx_size = nh.context_size
    dgst_size = nh.digest_size
    block_size = nh.block_size

    # Save the function pointers as well
    fptr_init = nh.init
    fptr_update = nh.update
    fptr_digest = nh.digest

    # First, create the type itself
    @eval immutable $name <: HashAlgorithm; end

    # Next, record all the important information about this hash algorithm
    @eval output_size(::Type{$name}) = $(convert(Int,dgst_size))
    @eval ctx_size(::Type{$name}) = $(convert(Int,ctx_size))
    @eval hash_type(::Type{$name}) = $nhptr

    # Generate the init, update, and digest functions while we're at it!
    # Since we have the function pointers from nh, we'll use those
    @eval function HashState(::Type{$name})
      ctx = Array(Uint8, ctx_size($name))
      ccall($fptr_init,Void,(Ptr{Void},),ctx)
      HashState{$name}(ctx)
    end

    @eval function update!(state::HashState{$name},data)
      ccall($fptr_update,Void,(Ptr{Void},Csize_t,Ptr{Uint8}),state.ctx,sizeof(data),pointer(data))
      state
    end

    @eval function digest!(state::HashState{$name})
      dgst = Array(Uint8,output_size($name))
      ccall($fptr_digest,Void,(Ptr{Void},Uint32,Ptr{Uint8}),state.ctx,sizeof(dgst),pointer(dgst))
      dgst
    end

    # Generate e.g. sha256_hash(string) and sha256_hmac(string)
    name_hash = symbol("$(bytestring(nh.name))_hash")
    @eval $name_hash(string) = digest!(update!(HashState($name), string))

    name_hmac = symbol("$(bytestring(nh.name))_hmac")
    @eval $name_hmac(key, string) = digest!(update!(HMACState($name, key), string))

    # Add this type into the HashAlgorithms group
    @eval push!(HashAlgorithms, $name)

    # Finally, export the type we just created
    for sym in [name, name_hash, name_hmac]
        eval(current_module(), Expr(:toplevel, Expr(:export, sym)))
    end

    hash_idx += 1
  end
end

function show{T<:HashAlgorithm}( io::IO, ::HashState{T} )
  write(io, "$(string(T)) Hash state")
end


text:
module FirstOrderTests

import ..Algopt # ensure we are using the correct Algopt

using Test

using LinearAlgebra

using Algopt.FirstOrder:
Termination,
MaximumGradientDescent, GradientDescent,
ConjugateGradientDescent, CGDProblemState,
search,
descent_step_maxgd, direction_maxgd,
descent_step_gd,
descent_step_cgd!

using Algopt.TestFunctions: Rosenbrock

quad0 = x->x'x
∇quad0 = x->2x

quad3 = x->(x - [3, 3])' * (x - [3, 3])
∇quad3 = x->2(x - [3, 3])

ros = Rosenbrock(a = 1, b = 5)

@testset "Algopt.FirstOrder.MaximumGradientDescent" begin
    @test [-0.6, -0.8] == direction_maxgd([6, 8])
    @test [0, -1] == direction_maxgd([0, 2])

    @test [0, 0] == descent_step_maxgd(quad0, ∇quad0, [3, 4])
    @test [0, 0] == descent_step_maxgd(quad0, ∇quad0, [10, 10])
    @test [3, 3] ≈ descent_step_maxgd(quad3, ∇quad3, [3, 4])
    @test [3, 3] ≈ descent_step_maxgd(quad3, ∇quad3, [10, 10])

    mgd = MaximumGradientDescent()
    x0 = [rand(-10:.1:10), rand(-10:.1:10)]
    @test norm([0, 0] - search(mgd, quad0, ∇quad0, x0)) < cbrt(eps())
    @test norm([0, 0] - search(mgd, quad0, ∇quad0, x0)) < cbrt(eps())
    @test [3, 3] ≈ search(mgd, quad3, ∇quad3, x0)
    @test [3, 3] ≈ search(mgd, quad3, ∇quad3, x0)
    @test norm([1, 1] - search(mgd, ros.f, ros.∇f, x0)) < cbrt(eps())
end

@testset "Algopt.FirstOrder.GradientDescent" begin
    α = 1e-2
    @test [3, 3.98] == descent_step_gd(α, quad3, ∇quad3, [3, 4])
    @test [9.86, 9.86] == descent_step_gd(α, quad3, ∇quad3, [10, 10])

    grd = GradientDescent()
    x0 = [rand(-10:.1:10), rand(-10:.1:10)]
    term = Termination(max_iter = 100_000)
    @test norm([0, 0] - search(grd, quad0, ∇quad0, x0; term = term)) < cbrt(eps())
    @test norm([0, 0] - search(grd, quad0, ∇quad0, x0; term = term)) < cbrt(eps())
    @test norm([3, 3] - search(grd, quad3, ∇quad3, x0; term = term)) < cbrt(eps())
    @test norm([3, 3] - search(grd, quad3, ∇quad3, x0; term = term)) < cbrt(eps())
    @test norm([1, 1] - search(grd, ros.f, ros.∇f, x0; term = term)) < cbrt(eps())
end

@testset "Algopt.FirstOrder.ConjugateGradientDescent" begin
    @test [0, 0] == descent_step_cgd!(CGDProblemState(2), quad0, ∇quad0, [3, 4])
    @test [0, 0] == descent_step_cgd!(CGDProblemState(2), quad0, ∇quad0, [10, 10])
    @test [3, 3] == descent_step_cgd!(CGDProblemState(2), quad3, ∇quad3, [3, 4])
    @test [3, 3] == descent_step_cgd!(CGDProblemState(2), quad3, ∇quad3, [10, 10])

    cgd = ConjugateGradientDescent()
    x0 = [rand(-10:.1:10), rand(-10:.1:10)]
    @test [0, 0] ≈ search(cgd, quad0, ∇quad0, x0)
    @test [0, 0] ≈ search(cgd, quad0, ∇quad0, x0)
    @test [3, 3] ≈ search(cgd, quad3, ∇quad3, x0)
    @test [3, 3] ≈ search(cgd, quad3, ∇quad3, x0)
    @test norm([1, 1] - search(cgd, ros.f, ros.∇f, x0)) < cbrt(eps())
end

end # module


text:
function corner_marker(r::Robot)
  mass = []
  go_to_nord_ost(r, mass) # Занять угол
  mark_corners!(r) # Покрасить углы
  back(r, mass) # Вернуться в начальное положение
end

function go_to_nord_ost(r::Robot, mass::Array)
  while isborder(r, Nord) == false || isborder(r, Ost) == false
    move_if_Possible!(r, mass, Nord)
    move_if_Possible!(r, mass, Ost)
  end
  # Nord Ost ist lage(DE)
end

function moves!(r::Robot, side::HorizonSide)
  while !isborder(r, side)
    move!(r, side)
  end
end

function mark_corners!(r::Robot)
  for side in (Sud, West, Nord, Ost)
    moves!(r, side)
    putmarker!(r)
  end
end

function inverse(side::HorizonSide)
  return HorizonSide(mod(Int(side) + 2, 4))
end

function back(r::Robot, mass::Array)
  while length(mass) > 0
    last_move = mass.pop()
    move!(r, inverse(last_move))
  end
end

function move_if_Possible!(r::Robot, mass::AbstractArray, side::HorizonSide)
  if !isborder(r, side)
      move!(r, side)
      push!(mass, side)
  end
end

max_stars_repo_name: NikitaDestroyer/Basics

max_issues_repo_name: NikitaDestroyer/Basics

max_forks_repo_name: NikitaDestroyer/Basics


text:
# some collection of some painting plus support functions

"""
function ngray(base::Int,digits::Int,value::Int)
Convert a value to a graycode with the given base and digits
"""
function ngray(base::Int,digits::Int,value::Int)
    baseN = zeros(Int,digits)
    gray = zeros(Int,1,digits)

	for i=1:digits
		baseN[i] = value % base
		value = div(value,base)
	end

	shift = 0
	for i=digits:-1:1
		gray[i] = (baseN[i] + shift) % base
		shift = shift + base - gray[i]
	end
	gray
end

"""
function hilbert_curve(c,x,y,lg,i1,i2)
recursive hilbert curve (2D), appends pairs of x,y to c, lg = dimension/length
"""
function hilbert_curve(c,x,y,lg,i1,i2)
    if lg == 1
        append!(c,[x,y])
    else
        lg = lg / 2;
        hilbert_curve(c,x+i1*lg,y+i1*lg,lg,i1,1-i2);
        hilbert_curve(c,x+i2*lg,y+(1-i2)*lg,lg,i1,i2);
        hilbert_curve(c,x+(1-i1)*lg,y+(1-i1)*lg,lg,i1,i2);
        hilbert_curve(c,x+(1-i2)*lg,y+i2*lg,lg,1-i1,i2);
    end
end

function hilbert_colored(surf)

    zscale = 8;
    n1 = 8;
    cr = CairoContext(surf)

    c = Float64[]
    hilbert_curve(c,0,0,64,0,0)

    move_to(cr,0,0)

    translate(cr,zscale/2,zscale/2)
    scale(cr,zscale,zscale)
    set_line_width(cr,zscale/2)
    set_line_cap(cr,Cairo.CAIRO_LINE_CAP_SQUARE)

    for k in zip(collect(1:2:(length(c)-2)),1:(length(1:2:(length(c)-2))))
        move_to(cr,c[k[1]],c[k[1]+1])
        line_to(cr,c[k[1]+2],c[k[1]+3])

        c1 = ngray(n1,3,k[2])
        set_source_rgb(cr,c1[1]/float(n1-1),c1[2]/float(n1-1),c1[3]/float(n1-1))
        stroke(cr)

    end
end

""" function hdraw(s,dim,zscale,linewidth)
draws a hilbert curve with dimension dim (power of 2) and scales the drawing with
zscale.
"""
function hdraw(s,dim,zscale,linewidth)

    cr = CairoContext(s)

    set_source_rgba(cr,0.0,0.0,1.0,0.5)

    set_line_width(cr,linewidth)
    set_line_cap(cr,Cairo.CAIRO_LINE_CAP_SQUARE)
    translate(cr,zscale/2,zscale/2)

    save(cr)
    c = Float64[]

    hilbert_curve(c,0,0,dim,0,0)

    scale(cr,zscale,zscale)

    move_to(cr,0,0)
    for k=1:div(length(c),2)
        line_to(cr,c[(k*2)-1],c[(k*2)])
    end

    stroke(cr)
    restore(cr)
end

"""
function simple_hist(data)
    simple histogram by population count in a Dict
"""
function simple_hist(data)
    # poor man's hist -> pop count in Dict

    pc = Dict()

    for d in data
        if d in keys(pc)
            pc[d] += 1
        else
            pc[d] = 1
        end
    end
    pc
end

"""
function matrix_read(surface)
	paint the input surface into a matrix image of the same size to access
	the pixels.
"""
function matrix_read(surface)
	w = Int(surface.width)
	h = Int(surface.height)
	z = zeros(UInt32,w,h)
	surf = CairoImageSurface(z, Cairo.FORMAT_ARGB32)

    cr = CairoContext(surf)
    set_source_surface(cr,surface,0,0)
    paint(cr)

    surf.data
end


text:
@testset "Pointer" begin
    @testset "constructor" begin
        pointer = WordNet.Pointer("s", 123, 'v', "1234")

        @test pointer.sym    == "s"
        @test pointer.offset == 123
        @test pointer.pos    == 'v'
        @test pointer.source == "12"
        @test pointer.target == "34"
    end

    @testset "is not semantic for non-0" begin
        ptr = WordNet.Pointer("s", 123, 'v', "1234")
        @test WordNet.is_semantic(ptr) == false
    end

    @testset "is semantic for all-0" begin
        ptr = WordNet.Pointer("s", 123, 'v', "0000")
        @test WordNet.is_semantic(ptr) == true
    end
end


text:
#####
##### `Mixer`
#####

struct Mixer{N,T<:NTuple{N,AbstractSynthComponent}} <: AbstractSynthComponent
    inputs::T
end

Mixer(inputs::AbstractSynthComponent...) = Mixer(inputs)

next!(m::Mixer) = sum(next!, m.inputs)

Base.:+(inputs::AbstractSynthComponent...) = Mixer(inputs...)
output_type(::Type{Mixer{N,T}}) where {N,T} = promote_op(+, output_type.(T.parameters)...)

#####
##### `Volume`
#####

struct Volume{S<:AbstractFloat,T<:AbstractSynthComponent} <: AbstractSynthComponent
    volume::S
    input::T
end

next!(v::Volume) = v.volume * next!(v.input)
output_type(::Type{Volume{S,T}}) where {S,T} = promote_op(*, S, output_type(T))

Base.:*(volume::Real, input::AbstractSynthComponent) = Volume(float(volume), input)


text:
GPS_PI = 3.1415926535898

struct TLM_HOW_Data_Struct
    integrity_status_flag::Bool
    TOW::Int64
    alert_flag::Bool
    anti_spoof_flag::Bool
end

struct Subframe_1_Data
    trans_week::Int64
    codeonl2::Int64
    ura::Float64
    svhealth::String
    IODC::String
    l2pcode::Bool
    T_GD::Float64
    t_oc::Int64
    a_f2::Float64
    a_f1::Float64
    a_f0::Float64
end

struct Subframe_2_Data
    IODE::String
    C_rs::Float64
    Δn::Float64
    M_0::Float64
    C_uc::Float64
    e::Float64
    C_us::Float64
    sqrt_A::Float64
    t_oe::Int64
    fit_interval::Bool
    AODO::Int64
end

struct Subframe_3_Data
    C_ic::Float64
    Ω_0::Float64
    C_is::Float64
    i_0::Float64
    C_rc::Float64
    ω::Float64
    Ω_dot::Float64
    IODE::String
    IDOT::Float64
end

struct Subframe_4_Data
end

struct Subframe_5_Data
end

"""
    Decodes words of the first two Words of a subframe
    $(SIGNATURES)

    ´words´: buffer of 60 bits, sliced in 2 Words á 30 Bits

    # Details
    # Decodes the first two Words of the Substring and returns a custom struct for these values
"""
function decode_TLM_HOW(words)
    integrity_status = words[1][23]
    tow = convert(Float64, bin2dec(words[2][1:17]))
    alert_flag = words[2][18]
    if alert_flag
        @warn "Signal URA may be worse than indicated in subframe 1 - Use satellite at own risk!"
    end
    anti_spoof_flag = words[2][19]
    return TLM_HOW_Data_Struct(integrity_status, tow, alert_flag, anti_spoof_flag)
end

"""
    Decodes words of the first subframe
    $(SIGNATURES)

    words: buffer of 300 bits, sliced in 10 Words á 30 Bits

    # Details
    # Decodes the first substring of the Substring and returns a custom struct for these values alongside the computed HOW and TLM Data
"""
function decode_subframe_1(words)
    println("Decoding subframe 1...")

    TLM_HOW = decode_TLM_HOW(words[1:2])
    # * Decoding Word 3
    # Transmission Week
    trans_week = convert(Float64, bin2dec(words[3][1:10]))

    # Codes on L2 Channel
    codeonl2 = convert(Float64, bin2dec(words[3][11:12]))

    if codeonl2 == 3 | 0
        @warn "Code on L2 Channel invalid!"
    end

    # SV Accuracy, user range accuracy
    ura  = bin2dec(words[3][13:16])
    if ura <= 6
        ura = 2^(1 + (ura / 2))
    elseif 6 < ura <= 14
        ura = 2^(ura - 2)
    elseif ura == 15
        @warn "URA unsafe, no accuracy prediction available - use Satellite on own risk!"
        ura = 99999
    end

    # Satellite Health
    svhealth = bitArray2Str(words[3][17:22])
    if words[3][17]
        @warn "Bad LNAV Data, SV-Health critical"
    end

    # Issue of Data Clock
    IODC = bitArray2Str(append!(words[3][23:24], words[8][1:8])) # 2 MSB in Word 2, LSB 8 in Word 8

    # * Decoding Word 4
    # True: LNAV Datastream on PCode commanded OFF
    l2pcode = words[4][1]

    # * Decoding Word 7
    # group time differential
    T_GD = convert(Float64, bin2dec_twoscomp(words[7][17:24])) * 2^-31

    # *Decoding Word 8
    # IODC already in computed in word 3

    # Clock data reference
    t_oc = bin2dec(words[8][9:24]) << 4

    # * Decoding Word 9
    # clock correction parameter a_f2
    a_f2 = bin2dec_twoscomp(words[9][1:8]) * 2^-55

    # clock correction parameter  a_f1
    a_f1 = bin2dec_twoscomp(words[9][9:24]) * 2^-43

    # * Decoding Word 10
    # clock correction parameter a_f0
    a_f0 = bin2dec_twoscomp(words[10][1:22]) * 2^-31

    # * Finish Decoding

    return TLM_HOW, Subframe_1_Data(trans_week, codeonl2, ura, svhealth, IODC, l2pcode, T_GD, t_oc, a_f2, a_f1, a_f0)
end

"""
    Decodes words of the second subframe
    $(SIGNATURES)

    ´words´: buffer of 300 bits, sliced in 10 Words á 30 Bits

    Details:
    # Decodes the second substring of the Substring and returns a custom struct for these values alongside the computed HOW and TLM Data
"""
function decode_subframe_2(words)
    println("Decoding subframe 2...")

    TLM_HOW = decode_TLM_HOW(words[1:2])
    # * Decoding Word 3
    # Issue of ephemeris data
    IODE = bitArray2Str(words[3][1:8])

    # Amplitude of Sine Harmonic Correction Term to Orbit Radius
    C_rs = bin2dec_twoscomp(words[3][9:24]) * 2^-5

    # * Decoding Word 4
    # Mean motion difference from computed value
    Δn = bin2dec_twoscomp(words[4][1:16]) * GPS_PI * 2^-43

    # Mean anomaly at Reference Time (From word 4 and 5)
    M_0 = bin2dec_twoscomp(append!(words[4][17:24], words[5][1:24])) * GPS_PI * 2^-31

    # * Decoding Word 5
    # Mean time anomaly computed in word 4

    # * Decoding Word 6
    # Amplitude of the Cosine Harmonic Correction Term to the Argument Latitude
    C_uc = bin2dec_twoscomp(words[6][1:16]) * 2^-29

    # Eccentricity
    e = bin2dec(append!(words[6][17:24], words[7][1:24])) * 2^-33

    # * Decoding Word 7
    # Eccentricity already computed in word 6

    # * Decoding Word 8
    # Amplitude of the Sine Harmonic Correction Term to the Argument of Latitude
    C_us = bin2dec_twoscomp(words[8][1:16]) * 2^-29

    # Square Root of Semi-Major Axis
    sqrt_A = bin2dec(append!(words[8][17:24], words[9][1:24])) * 2^-19

    # * Decoding Word 9
    # square of A already computed in Word 8

    # * Decoding Word 10
    # Reference Time ephemeris
    t_oe = bin2dec(words[10][1:16]) << 4

    # Curve ftir interval flag - (0: 4 hours | 1: greater than 4 hours)
    fitinterval = words[10][17]

    # AODO Word
    aodo = bin2dec(words[10][18:22])

    # * Finish Decoding
    return TLM_HOW, Subframe_2_Data(IODE, C_rs, Δn, M_0, C_uc, e, C_us, sqrt_A, t_oe, fitinterval, aodo)
end

"""
    Decodes words of the third subframe
    $(SIGNATURES)

    ´words´: buffer of 300 bits, sliced in 10 Words á 30 Bits

    # Details
    # Decodes the third substring of the Substring and returns a custom struct for these values alongside the computed HOW and TLM Data
"""
function decode_subframe_3(words)

    println("Decoding subframe 3...")

    TLM_HOW = decode_TLM_HOW(words[1:2])
    # * Decoding Word 3
    # Amplitude of the Cosine Harmonic Correction to Angle of Inclination
    C_ic = bin2dec_twoscomp(words[3][1:16]) * 2^-29

    # Longitude of Ascending Node of Orbit Plane at Weekly Epoch
    Ω_0 = bin2dec_twoscomp(append!(words[3][17:24], words[4][1:24])) * GPS_PI * 2^-31
    # * Decoding Word 4
    # Omega 0 already in word 3 computed

    # * Decoding Word 5
    # Amplitude of the sine harmonic correction term to angle of Inclination
    C_is = bin2dec_twoscomp(words[5][1:16]) * 2^-29

    # inclination Angle at reference time
    i_0 = bin2dec_twoscomp(append!(words[5][17:24], words[6][1:24])) * GPS_PI * 2^-31

    # * Decoding Word 6
    # i_0 already in Word 5 computed

    # * Decoding Word 7
    # Amplitude of the cosine harmonic correction term to orbit Radius
    C_rc = bin2dec_twoscomp(words[7][1:16]) * 2^-5

    # Argument of Perigee
    ω = bin2dec_twoscomp(append!(words[7][17:24], words[8][1:24])) * GPS_PI * 2^-31

    # * Decoding Word 8
    # Argument of Perigee already computed in word 7

    # * Decoding Word 9
    # Rate of Right Ascension
    Ω_dot = bin2dec_twoscomp(words[9][1:24]) * GPS_PI * 2^-43

    # * Decoding Word 10
    # Issue of Ephemeris Data
    IODE = bitArray2Str(words[10][1:8])

    # Rate of Inclination Angle
    IDOT = bin2dec_twoscomp(words[10][9:22]) * GPS_PI * 2^-43

    # * Finish Decoding
    return TLM_HOW, Subframe_3_Data(C_ic, Ω_0, C_is, i_0, C_rc, ω, Ω_dot, IODE, IDOT)
end

"""
    Decodes words of the fourth subframe
    $(SIGNATURES)

    ´words´: buffer of 300 bits, sliced in 10 Words á 30 Bits

    # Details
    # Does mot do anything at the Moment, except calling the TLM_HOW Decoder and returns its return
"""
function decode_subframe_4(words)
    println("Decoding subframe 4...")

    TLM_HOW = decode_TLM_HOW(words[1:2])
    return TLM_HOW
end

"""
    Decodes words of the fifth subframe
    $(SIGNATURES)

    ´words´: buffer of 300 bits, sliced in 10 Words á 30 Bits

    # Details
    # Does mot do anything at the Moment, except calling the TLM_HOW Decoder and returns its return
"""
function decode_subframe_5(words)
    println("Decoding subframe 5...")

    TLM_HOW = decode_TLM_HOW(words[1:2])
    return TLM_HOW
end

function create_data(
    TLM_HOW_Data::TLM_HOW_Data_Struct,
    subfr_1_data::Subframe_1_Data,
    subfr_2_data::Subframe_2_Data,
    subfr_3_data::Subframe_3_Data)

    data = GPSData(
        TLM_HOW_Data.integrity_status_flag,
        TLM_HOW_Data.TOW,
        TLM_HOW_Data.alert_flag,
        TLM_HOW_Data.anti_spoof_flag,

        subfr_1_data.trans_week,
        subfr_1_data.codeonl2,
        subfr_1_data.ura,
        subfr_1_data.svhealth,
        subfr_1_data.IODC,
        subfr_1_data.l2pcode,
        subfr_1_data.T_GD,
        subfr_1_data.t_oc,
        subfr_1_data.a_f2,
        subfr_1_data.a_f1,
        subfr_1_data.a_f0,

        subfr_2_data.IODE,
        subfr_2_data.C_rs,
        subfr_2_data.Δn,
        subfr_2_data.M_0,
        subfr_2_data.C_uc,
        subfr_2_data.e,
        subfr_2_data.C_us,
        subfr_2_data.sqrt_A,
        subfr_2_data.t_oe,
        subfr_2_data.fit_interval,
        subfr_2_data.AODO,

        subfr_3_data.C_ic,
        subfr_3_data.Ω_0,
        subfr_3_data.C_is,
        subfr_3_data.i_0,
        subfr_3_data.C_rc,
        subfr_3_data.ω,
        subfr_3_data.Ω_dot,
        subfr_3_data.IODE,
        subfr_3_data.IDOT
    )
    return data
end


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

abstract type AbstractScenarioProblems{S <: AbstractScenario} end

abstract type AbstractDecompositionStructure{N} <: AbstractStochasticStructure{N} end

# JuMP #
# ========================== #
function JuMP.objective_function_type(structure::AbstractDecompositionStructure{N}, stage::Integer, scenario_index::Integer) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    stage > 1 || error("There are no scenarios in the first stage.")
    n = num_scenarios(structure, stage)
    1 <= scenario_index <= n || error("Scenario index $scenario_index not in range 1 to $n.")
    return get_from_scenarioproblem(scenarioproblems(structure, stage), scenario_index) do sp, i
        s = fetch(sp).problems[i]
        return jump_function_type(s, MOI.get(backend(s), MOI.ObjectiveFunctionType()))
    end
end

function JuMP.objective_function(structure::AbstractDecompositionStructure{N},
                                 stage::Integer,
                                 scenario_index::Integer,
                                 FunType::Type{<:AbstractJuMPScalar}) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    stage > 1 || error("There are no scenarios in the first stage.")
    n = num_scenarios(structure, stage)
    1 <= scenario_index <= n || error("Scenario index $scenario_index not in range 1 to $n.")
    obj = get_from_scenarioproblem(scenarioproblems(structure, stage), scenario_index, FunType) do sp, i, FunType
        MOIFunType = moi_function_type(FunType)
        subprob = fetch(sp).problems[i]
        obj = MOI.get(subprob, MOI.ObjectiveFunction{MOIFunType}())::MOIFunType
        return obj
    end
    objective = jump_function(structure, stage, scenario_index, obj)::FunType
    return objective
end

function JuMP._moi_optimizer_index(structure::AbstractDecompositionStructure, ci::CI, scenario_index::Integer)
    return JuMP._moi_optimizer_index(scenarioproblems(structure), ci, scenario_index)
end

function DecisionRef(structure::AbstractDecompositionStructure, index::VI, stage::Integer, scenario_index::Integer)
    return DecisionRef(structure.proxy[stage], index)
end
function DecisionRef(structure::AbstractDecompositionStructure, index::VI, at_stage::Integer, stage::Integer, scenario_index::Integer)
    at_stage > 1 || error("There are no scenarios in the first at_stage.")
    n = num_scenarios(structure, at_stage)
    1 <= scenario_index <= n || error("Scenario index $scenario_index not in range 1 to $n.")
    return DecisionRef(structure.proxy[at_stage], index)
end

# Getters #
# ========================== #
function scenarioproblems(structure::AbstractDecompositionStructure{N}, stage::Integer) where N
    1 < stage <= N || error("Stage $stage not in range 2 to $N.")
    stage == 1 && error("Stage 1 does not have scenario problems.")
    N == 2 && (stage == 2 || error("Stage $stage not available in two-stage model."))
    return structure.scenarioproblems[stage-1]
end
function scenarioproblems(structure::AbstractDecompositionStructure{2})
    return scenarioproblems(structure, 2)
end
function scenario_types(structure::AbstractDecompositionStructure{N}) where N
    return ntuple(Val{N-1}()) do s
        scenario_type(scenarioproblems(structure), s + 1)
    end
end
function proxy(structure::AbstractDecompositionStructure{N}, stage::Integer) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    return structure.proxy[stage]
end
function scenario(structure::AbstractDecompositionStructure, stage::Integer, scenario_index::Integer)
    scenario(scenarioproblems(structure, stage), scenario_index)
end
function scenario(structure::AbstractDecompositionStructure{2}, scenario_index::Integer)
    scenario(scenarioproblems(structure, 2), scenario_index)
end
function scenarios(structure::AbstractDecompositionStructure, stage::Integer)
    scenarios(scenarioproblems(structure, stage))
end
function scenarios(structure::AbstractDecompositionStructure{2})
    scenarios(scenarioproblems(structure, 2))
end
function expected(structure::AbstractDecompositionStructure, stage::Integer)
    return expected(scenarioproblems(structure, stage))
end
function expected(structure::AbstractDecompositionStructure{2})
    return expected(scenarioproblems(structure, 2))
end
function scenario_type(structure::AbstractDecompositionStructure, stage::Integer)
    return scenario_type(scenarioproblems(structure, stage))
end
function scenario_type(structure::AbstractDecompositionStructure{2})
    return scenario_type(scenarioproblems(structure, 2))
end
function probability(structure::AbstractDecompositionStructure, stage::Integer, scenario_index::Integer)
    return probability(scenarioproblems(structure, stage), scenario_index)
end
function probability(structure::AbstractDecompositionStructure{2}, scenario_index::Integer)
    return probability(scenarioproblems(structure, 2), scenario_index)
end
function stage_probability(structure::AbstractDecompositionStructure, stage::Integer)
    return probability(scenarioproblems(structure, stage))
end
function stage_probability(structure::AbstractDecompositionStructure{2})
    return probability(scenarioproblems(structure, 2))
end
function subproblem(structure::AbstractDecompositionStructure, stage::Integer, scenario_index::Integer)
    return subproblem(scenarioproblems(structure, stage), scenario_index)
end
function subproblem(structure::AbstractDecompositionStructure{2}, scenario_index::Integer)
    return subproblem(scenarioproblems(structure, 2), scenario_index)
end
function subproblems(structure::AbstractDecompositionStructure, stage::Integer)
    return subproblems(scenarioproblems(structure, stage))
end
function subproblems(structure::AbstractDecompositionStructure{2})
    return subproblems(scenarioproblems(structure, 2))
end
function num_subproblems(structure::AbstractDecompositionStructure, stage::Integer)
    return num_subproblems(scenarioproblems(structure, stage))
end
function num_subproblems(structure::AbstractDecompositionStructure{2})
    return num_subproblems(scenarioproblems(structure, 2))
end
function num_scenarios(structure::AbstractDecompositionStructure, stage::Integer)
    return num_scenarios(scenarioproblems(structure, stage))
end
function num_scenarios(structure::AbstractDecompositionStructure{2})
    return num_scenarios(scenarioproblems(structure, 2))
end
deferred(structure::AbstractDecompositionStructure{N}) where N = deferred(structure, Val(N))
deferred(structure::AbstractDecompositionStructure, ::Val{1}) = deferred_first_stage(structure)
function deferred(structure::AbstractDecompositionStructure, ::Val{N}) where N
    return deferred_stage(structure, N) || deferred(structure, Val(N-1))
end
deferred_first_stage(structure::AbstractDecompositionStructure) = false
function deferred_stage(structure::AbstractDecompositionStructure{N}, stage::Integer) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    stage == 1 && return deferred_first_stage(structure)
    num_subproblems(structure, stage) < num_scenarios(structure, stage)
end
function distributed(structure::AbstractDecompositionStructure{N}, stage::Integer) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    stage == 1 && return false
    return distributed(scenarioproblems(structure, stage))
end

# ========================== #

# Setters
# ========================== #
function update_known_decisions!(structure::AbstractDecompositionStructure, change::DecisionModification, stage::Integer, scenario_index::Integer)
    update_known_decisions!(scenarioproblems(structure, stage), change, scenario_index)
    return nothing
end
function add_scenario!(structure::AbstractDecompositionStructure, stage::Integer, scenario::AbstractScenario)
    add_scenario!(scenarioproblems(structure, stage), scenario)
    return nothing
end
function add_worker_scenario!(structure::AbstractDecompositionStructure, stage::Integer, scenario::AbstractScenario, w::Integer)
    add_scenario!(scenario(structure, stage), scenario, w)
    return nothing
end
function add_scenario!(scenariogenerator::Function, structure::AbstractDecompositionStructure, stage::Integer)
    add_scenario!(scenariogenerator, scenarioproblems(structure, stage))
    return nothing
end
function add_worker_scenario!(scenariogenerator::Function, structure::AbstractDecompositionStructure, stage::Integer, w::Integer)
    add_scenario!(scenariogenerator, scenarioproblems(structure, stage), w)
    return nothing
end
function add_scenarios!(structure::AbstractDecompositionStructure, stage::Integer, scenarios::Vector{<:AbstractScenario})
    add_scenarios!(scenarioproblems(structure, stage), scenarios)
    return nothing
end
function add_worker_scenarios!(structure::AbstractDecompositionStructure, stage::Integer, scenarios::Vector{<:AbstractScenario}, w::Integer)
    add_scenarios!(scenarioproblems(structure, stage), scenarios, w)
    return nothing
end
function add_scenarios!(scenariogenerator::Function, structure::AbstractDecompositionStructure, stage::Integer, n::Integer)
    add_scenarios!(scenariogenerator, scenarioproblems(structure, stage), n)
    return nothing
end
function add_worker_scenarios!(scenariogenerator::Function, structure::AbstractDecompositionStructure, stage::Integer, n::Integer, w::Integer)
    add_scenarios!(scenariogenerator, scenarioproblems(structure, stage), n, w)
    return nothing
end
function sample!(structure::AbstractDecompositionStructure, stage::Integer, sampler::AbstractSampler, n::Integer)
    sample!(scenarioproblems(structure, stage), sampler, n)
    return nothing
end
# ========================== #

# Includes
# ========================== #
include("scenarioproblems.jl")
include("stage_decomposition.jl")
include("scenario_decomposition.jl")