text:
using PyPlot
using LinearAlgebra
using PolynomialBasis
using ImplicitDomainQuadrature
using Revise
using CutCell
include("../../../test/useful_routines.jl")

function bulk_modulus(l, m)
    return l + 2m / 3
end

function lame_lambda(k, m)
    return k - 2m / 3
end

function onleftboundary(x, L, W)
    return x[1] ≈ 0.0
end

function onbottomboundary(x, L, W)
    return x[2] ≈ 0.0
end

function onrightboundary(x, L, W)
    return x[1] ≈ L
end

function ontopboundary(x, L, W)
    return x[2] ≈ W
end

function solve_for_displacement(
    basis,
    cellquads,
    interfacequads,
    facequads,
    cutmesh,
    levelset,
    levelsetcoeffs,
    stiffness,
    transfstress,
    penalty,
    applydisplacement
)

    L,W = CutCell.widths(cutmesh)

    bilinearforms = CutCell.BilinearForms(basis, cellquads, stiffness, cutmesh)
    interfacecondition = CutCell.coherent_interface_condition(
        basis,
        interfacequads,
        stiffness,
        cutmesh,
        penalty,
    )

    leftbc = CutCell.DisplacementComponentCondition(
        x -> 0.0,
        basis,
        facequads,
        stiffness,
        cutmesh,
        x->onleftboundary(x,L,W),
        [1.0, 0.0],
        penalty,
    )
    bottombc = CutCell.DisplacementComponentCondition(
        x -> 0.0,
        basis,
        facequads,
        stiffness,
        cutmesh,
        x->onbottomboundary(x,L,W),
        [0.0, 1.0],
        penalty,
    )
    rightbc = CutCell.DisplacementComponentCondition(
        x -> applydisplacement,
        basis,
        facequads,
        stiffness,
        cutmesh,
        x->onrightboundary(x,L,W),
        [1.0, 0.0],
        penalty,
    )

    sysmatrix = CutCell.SystemMatrix()
    sysrhs = CutCell.SystemRHS()

    CutCell.assemble_bilinear_form!(sysmatrix, bilinearforms, cutmesh)
    CutCell.assemble_interface_condition!(sysmatrix, interfacecondition, cutmesh)
    CutCell.assemble_bulk_transformation_linear_form!(
        sysrhs,
        transfstress,
        basis,
        cellquads,
        cutmesh,
    )

    CutCell.assemble_coherent_interface_transformation_rhs!(
        sysrhs,
        transfstress,
        basis,
        interfacequads,
        cutmesh,
    )

    CutCell.assemble_penalty_displacement_bc!(sysmatrix,sysrhs,leftbc,cutmesh)
    CutCell.assemble_penalty_displacement_component_transformation_rhs!(
        sysrhs,
        transfstress,
        basis,
        facequads,
        cutmesh,
        x -> onleftboundary(x, L, W),
        [1.0,0.0]
    )

    CutCell.assemble_penalty_displacement_bc!(sysmatrix,sysrhs,bottombc,cutmesh)
    CutCell.assemble_penalty_displacement_component_transformation_rhs!(
        sysrhs,
        transfstress,
        basis,
        facequads,
        cutmesh,
        x -> onbottomboundary(x, L, W),
        [0.0,1.0]
    )

    CutCell.assemble_penalty_displacement_bc!(sysmatrix,sysrhs,rightbc,cutmesh)
    CutCell.assemble_penalty_displacement_component_transformation_rhs!(
        sysrhs,
        transfstress,
        basis,
        facequads,
        cutmesh,
        x -> onrightboundary(x, L, W),
        [1.0,0.0]
    )

    matrix = CutCell.make_sparse(sysmatrix, cutmesh)
    rhs = CutCell.rhs(sysrhs, cutmesh)

    nodaldisplacement = matrix \ rhs

    return nodaldisplacement
end

function spatial_closest_points(refclosestpoints, refclosestcellids, mesh)
    dim, npts = size(refclosestpoints)
    spclosestpoints = zeros(dim, npts)
    for i = 1:npts
        cellmap = CutCell.cell_map(mesh, refclosestcellids[i])
        spclosestpoints[:, i] .= cellmap(refclosestpoints[:, i])
    end
    return spclosestpoints
end

function normal_stress_component(stressvector, normal)
    snn =
        normal[1] * stressvector[1] * normal[1] +
        2.0 * normal[1] * stressvector[3] * normal[2] +
        normal[2] * stressvector[2] * normal[2]
end

function normal_stress_component_over_points(stresses,normals)
    numstresscomponents,npts = size(stresses)
    @assert size(normals) == (2,npts)

    normalstresscomp = zeros(npts)
    for i = 1:npts
        normalstresscomp[i] = normal_stress_component(stresses[:,i],normals[:,i])
    end
    return normalstresscomp
end

function stress_inner_product_over_points(stresses)
    return stresses[1,:].^2 + stresses[2,:].^2 + 2.0*stresses[3,:].^2 + stresses[4,:].^2
end

K1, K2 = 247.0e9, 162.0e9    # Pa
mu1, mu2 = 126.0e9, 126.0e9   # Pa
lambda1 = lame_lambda(K1, mu1)
lambda2 = lame_lambda(K2, mu2)
stiffness = CutCell.HookeStiffness(lambda1, mu1, lambda2, mu2)

rho1 = 3.93e3   # Kg/m^3
rho2 = 3.93e3   # Kg/m^3

V1 = inv(rho1)
V2 = inv(rho2)
theta0 = 0.0
transfstress = CutCell.plane_strain_transformation_stress(lambda1, mu1, theta0)
ΔG = -1022.0    # J/mol
molarmass = 0.147   # Kg/mol
diffG0 = ΔG / molarmass # J/Kg

width = 1.0e-3
penaltyfactor = 1e2
applydisplacement = 0.01*width

polyorder = 1
numqp = required_quadrature_order(polyorder) + 2
nelmts = 11
interfacex0 = [0.8width,0.0]
interfaceangle = 30.0
interfacenormal = [cosd(interfaceangle),sind(interfaceangle)]

dx = width / nelmts
meanmoduli = 0.5 * (lambda1 + lambda2 + mu1 + mu2)
penalty = penaltyfactor / dx * meanmoduli

basis = TensorProductBasis(2, polyorder)
mesh = CutCell.Mesh([0.0, 0.0], [width, width], [nelmts, nelmts], basis)
levelset = InterpolatingPolynomial(1, basis)
levelsetcoeffs =
    CutCell.levelset_coefficients(x -> plane_distance_function(x,interfacenormal,interfacex0), mesh)

cutmesh = CutCell.CutMesh(levelset, levelsetcoeffs, mesh)
cellquads = CutCell.CellQuadratures(levelset, levelsetcoeffs, cutmesh, numqp)
interfacequads = CutCell.InterfaceQuadratures(levelset, levelsetcoeffs, cutmesh, numqp)
facequads = CutCell.FaceQuadratures(levelset, levelsetcoeffs, cutmesh, numqp)

nodaldisplacement = solve_for_displacement(
    basis,
    cellquads,
    interfacequads,
    facequads,
    cutmesh,
    levelset,
    levelsetcoeffs,
    stiffness,
    transfstress,
    penalty,
    applydisplacement
)

refseedpoints, spatialseedpoints, seedcellids =
    CutCell.seed_zero_levelset_with_interfacequads(interfacequads, cutmesh)
nodalcoordinates = CutCell.nodal_coordinates(cutmesh)

tol = 1e-8
boundingradius = 3.0
refclosestpoints, refclosestcellids, refgradients =
    CutCell.closest_reference_points_on_zero_levelset(
        nodalcoordinates,
        refseedpoints,
        spatialseedpoints,
        seedcellids,
        levelset,
        levelsetcoeffs,
        mesh,
        tol,
        boundingradius,
    )
spclosestpoints = spatial_closest_points(refclosestpoints,refclosestcellids,cutmesh)

sortidx = sortperm(spclosestpoints[2,:])

spclosestpoints = spclosestpoints[:,sortidx]
refclosestpoints = refclosestpoints[:,sortidx]
refclosestcellids = refclosestcellids[sortidx]
refgradients = refgradients[:,sortidx]

invjac = CutCell.inverse_jacobian(cutmesh)
normals = diagm(invjac) * refgradients
CutCell.normalize_normals!(normals)

product_stress_at_cp = CutCell.product_stress_at_reference_points(
    refclosestpoints,
    refclosestcellids,
    basis,
    stiffness,
    transfstress,
    theta0,
    nodaldisplacement,
    cutmesh,
)
parent_stress_at_cp = CutCell.parent_stress_at_reference_points(
    refclosestpoints,
    refclosestcellids,
    basis,
    stiffness,
    nodaldisplacement,
    cutmesh,
)

productpressure = CutCell.pressure_at_points(product_stress_at_cp)
parentpressure = CutCell.pressure_at_points(parent_stress_at_cp)

productdevstress = CutCell.deviatoric_stress_at_points(product_stress_at_cp,productpressure)
parentdevstress = CutCell.deviatoric_stress_at_points(parent_stress_at_cp,parentpressure)

productdevnorm = stress_inner_product_over_points(productdevstress)
parentdevnorm = stress_inner_product_over_points(parentdevstress)

productnormaldevstress = normal_stress_component_over_points(productdevstress,normals)
parentnormaldevstress = normal_stress_component_over_points(parentdevstress,normals)

productspecificvolume = V1*(1.0 .- productpressure/K1)
parentspecificvolume = V2*(1.0 .- parentpressure/K2)

productp1 = productspecificvolume .* productpressure
productp2 = V1/(2K1)*(productpressure.^2)
productp3 = -(productspecificvolume .* productnormaldevstress)
productp4 = V1/(4mu1)*productdevnorm
productpotential = productp1+productp2+productp3+productp4

parentp1 = parentspecificvolume .* parentpressure
parentp2 = V2/(2K2)*(parentpressure.^2)
parentp3 = -(parentspecificvolume .* parentnormaldevstress)
parentp4 = V2/(4mu2)*parentdevnorm
parentpotential = parentp1+parentp2+parentp3+parentp4

potentialdifference = productpotential - parentpotential

folderpath = "examples/transformation-strain/plane-interface/"
cpycoords = spclosestpoints[2,:]

# fig,ax = PyPlot.subplots()
# ax.plot(cpycoords,parentp1,label="term1")
# ax.plot(cpycoords,parentp2,label="term2")
# ax.plot(cpycoords,parentp3,label="term3")
# ax.plot(cpycoords,parentp4,label="term4")
# ax.plot(cpycoords,parentpotential,label="potential")
# ax.set_xlabel("Height (mm)")
# ax.set_ylabel("Energy density (J/Kg)")
# ax.set_title("Potential on parent side of interface")
# ax.grid()
# ax.legend()
# fig
# fig.tight_layout()
# fig.savefig(folderpath*"parent-potential.png")

# fig,ax = PyPlot.subplots()
# ax.plot(cpycoords,productp1,label="term1")
# ax.plot(cpycoords,productp2,label="term2")
# ax.plot(cpycoords,productp3,label="term3")
# ax.plot(cpycoords,productp4,label="term4")
# ax.plot(cpycoords,productpotential,label="potential")
# ax.set_xlabel("Height (mm)")
# ax.set_ylabel("Energy density (J/Kg)")
# ax.set_title("Potential on parent side of interface")
# ax.grid()
# ax.legend()
# fig
# fig.tight_layout()
# fig.savefig(folderpath*"parent-potential.png")

fig,ax = PyPlot.subplots()
ax.plot(cpycoords,potentialdifference)
ax.grid()
ax.set_xlabel("Height (mm)")
ax.set_ylabel("Potential difference (J/Kg)")
ax.set_title("Potential difference along interface")
# ax.set_ylim(-1,1)
fig.tight_layout()
fig
# fig.savefig(folderpath*"potential-difference.png")


text:
"""
    to_cartesian(N, args::Index...)

Convert `args` to a tuple of `CartesianIndex{N}`. Each item in `args` is either `Int`, `NTuple{N,
Int}` or `CartesianIndex{N}`; if it is `Int` it will be repeated, i.e., `ntuple(_->i, N)`.

# Examples

```jldoctest; setup=:(using BlockMatching: to_cartesian)
julia> to_cartesian(2, 3, 1, 4)
(CartesianIndex(3, 3), CartesianIndex(1, 1), CartesianIndex(4, 4))

julia> to_cartesian(2, CartesianIndex(1, 3), 1, (4, 4))
(CartesianIndex(1, 3), CartesianIndex(1, 1), CartesianIndex(4, 4))
```
"""
to_cartesian(N::Int, args...) = to_cartesian(Val(N), args...)
to_cartesian(::Val{N}, args::Union{Int, NTuple{N, Int}, CartesianIndex{N}}...) where N = map(i->_to_cartesian(Val(N), i), args)
_to_cartesian(::Val{N}, i::Int) where N = __to_cartesian(ntuple(_->i, N))
_to_cartesian(::Val{N}, t::NTuple{N, Int}) where N = __to_cartesian(t)
_to_cartesian(::Val{N}, t::CartesianIndex{N}) where N = t
__to_cartesian(t::NTuple{N, Int}) where N = CartesianIndex(t)

function size_check(ref, frame)
    Base.require_one_based_indexing(ref)
    Base.require_one_based_indexing(frame)
    any(size(ref) .< size(frame)) && throw(ArgumentError("`size(ref) = $(size(ref))` should be larger than `size(frame) = $(size(frame))`."))
end

###
# Min/Max heap-based Top-K algorithm
#
# Modified from https://github.com/JuliaCollections/DataStructures.jl/blob/0230564f0c88a3f65815cb98a8fe78c2658d8db8/src/heaps/arrays_as_heaps.jl
#
# No longer used but kept here in case we still need it in the future.
###
heapleft(i::Integer) = 2i
heapright(i::Integer) = 2i + 1
heapparent(i::Integer) = div(i, 2)

function percolate_down!(vals, inds, i, o::Base.Ordering)
    x, p, len = vals[i], inds[i], length(vals)
    @inbounds while (l = heapleft(i)) <= len
        r = heapright(i)
        j = r > len || Base.Order.lt(o, vals[l], vals[r]) ? l : r
        if Base.Order.lt(o, vals[j], x)
            vals[i] = vals[j]
            inds[i] = inds[j]
            i = j
        else
            break
        end
    end
    vals[i] = x
    inds[i] = p
    return vals, inds
end

function heapify!(vals, inds, o::Base.Ordering)
    for i in heapparent(length(vals)):-1:1
        percolate_down!(vals, inds, i, o)
    end
    return vals, inds
end

function topk!(vals, inds, data, o::Base.Ordering)
    @assert length(vals) == length(inds) <= length(data)
    @inbounds for p = 1:length(data)
        x = data[p]
        if Base.Order.lt(o, first(vals), x)
            vals[1] = x
            inds[1] = p
            heapify!(vals, inds, o)
        end
    end
    return vals, inds
end

"""
    topk(data, K, o=Forward) -> (vals, inds)

Given input `data`, find the largest `K` numbers `vals` and their associated indices `inds`.
If `o=Base.Order.Reverse`, then it's the lowest `K` numbers.

```jldoctest; setup=:(using BlockMatching, Random; Random.seed!(1234))
julia> data = rand(1:20, 20);

julia> vals, inds = BlockMatching.topk(data, 5)
([16, 16, 19, 19, 19], [16, 14, 7, 3, 17])

julia> data[inds] == vals
true
```
"""
function topk(data::AbstractArray{T}, K::Integer, o::Base.Ordering=Base.Forward) where T
    mn, mx = typemin(T), typemax(T)
    fillvalue = Base.Order.lt(o, mn, mx) ? mn : mx
    vals = fill(fillvalue, K)
    inds = fill(0, K)
    return topk!(vals, inds, data, o)
end

ext: jl


text:
using LinearAlgebra
using Random
using Statistics

include("../utils/simula01.jl")
include("../ACQR/ACQR_simula.jl")
include("../utils/byrow.jl")

include("ACQRS_sub.jl")


text:
using .Hidato       # Note that the . here means to look locally for the module rather than in the libraries

const hopid = """
 . 0 0 . 0 0 .
 0 0 0 0 0 0 0
 0 0 0 0 0 0 0
 . 0 0 0 0 0 .
 . . 0 0 0 . .
 . . . 0 . . . """

const hopidomoves = [[-3, 0], [0, -3], [-2, -2], [-2, 2], [2, -2], [0, 3], [3, 0], [2, 2]]

board, maxmoves, fixed, starts = hidatoconfigure(hopid)
printboard(board, " 0", "  ")
hidatosolve(board, maxmoves, hopidomoves, fixed, starts[1][1], starts[1][2], 1)
printboard(board)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using PlutoCLI
using Test

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


text:
using Makie
using StaticArrays
using Printf

room_count = [0]  # global
v = Vector{Int64}()  # global

struct NODE
    dir_pair::String
    x::Union{Nothing,Float64}
    y::Union{Nothing,Float64}
    function NODE(; x=nothing, y=nothing, dir_pair::String)
        return new(dir_pair, x, y)
    end
end
is_x(n::NODE) = !isnothing(n.x)
is_y(n::NODE) = !isnothing(n.y)

###

mutable struct Attached_Type
    room
    dir::String
    Δ_pos::Float64
    door_pos_in_child::Float64
    function Attached_Type(room, dir::String, Δ_pos::Float64, door_pos_in_child::Float64)
        return new(room, dir, Δ_pos, door_pos_in_child)
    end
end

###

mutable struct Room
    id::Int64
    x::Float64
    y::Float64
    refine_method
    name::String
    attached::Vector{Attached_Type}
    nodes::Vector{NODE}
    scene_element
    is_active::Bool
    global_bounds::NTuple{4,Float64}
    is_pos_x::Bool
    is_pos_y::Bool
    function Room(x, y, refine_method, name::String)
        attached = Vector{Attached_Type}()
        nodes = Vector{NODE}()
        room_count[1] += 1
        return new(room_count[1], x, y, refine_method, name, attached, nodes, nothing, false)
    end
end

get_Δx(r::Room) = r.global_bounds[2] -  r.global_bounds[1]
get_Δy(r::Room) = r.global_bounds[4] -  r.global_bounds[3]
get_x0(r::Room) = r.global_bounds[1]
get_x1(r::Room) = r.global_bounds[2]
get_y0(r::Room) = r.global_bounds[3]
get_y1(r::Room) = r.global_bounds[4]

get_nx(r::Room) = get_nx(r.refine_method)
get_ny(r::Room) = get_ny(r.refine_method)

function search_for_x_constraint(r::Room)
    nodes = r.nodes
    for node_k = nodes
        is_x(node_k) && (return node_k.x)
    end
    return nothing
end
function search_for_y_constraint(r::Room)
    nodes = r.nodes
    for node_k = nodes
        is_y(node_k) && (return node_k.y)
    end
    return nothing
end

is_con_y(r::Room) = !isnothing(search_for_y_constraint(r))
is_con_x(r::Room) = !isnothing(search_for_x_constraint(r))

###

abstract type RefineType end

mutable struct RefineHoriz <: RefineType
    n::Int64
    function RefineHoriz()
        return new(1)
    end
end
mutable struct RefineVert <: RefineType
    n::Int64
    function RefineVert()
        return new(1)
    end
end
mutable struct RefineRectUR <: RefineType
    nx::Int64
    ny::Int64
    function RefineRectUR()
        return new(1, 1)
    end
end
mutable struct RefineRectDR <: RefineType
    nx::Int64
    ny::Int64
    function RefineRectDR()
        return new(1, 1)
    end
end
mutable struct RefineTri <: RefineType end

refine_x!(r::Union{RefineRectUR,RefineRectDR}) = (r.nx += 1; return nothing)
refine_y!(r::Union{RefineRectUR,RefineRectDR}) = (r.ny += 1; return nothing)
refine_x!(r::Union{RefineHoriz,RefineVert})    = (r.n += 1;  return nothing)
refine_y!(r::Union{RefineHoriz,RefineVert})    = (r.n += 1;  return nothing)

get_nx(r::Union{RefineRectUR,RefineRectDR}) = r.nx
get_ny(r::Union{RefineRectUR,RefineRectDR}) = r.ny

get_nx(r::RefineHoriz) = r.n
get_ny(r::RefineHoriz) = 1

get_nx(r::RefineVert) = 1
get_ny(r::RefineVert) = r.n

###

function add_attached!(r::Room, r_new::Room, dir::String, offset::Float64, opening::Float64)
    a = Attached_Type(r_new, dir, offset, opening)
    push!(r.attached, a)
end

plot_room!(r::Room) = plot_room!(r, 0.0, 0.0)

function plot_room!(r::Room, x0::Float64, y0::Float64)
    x1 = r.x + x0
    y1 = r.y + y0
    x = [x0, x1, x1, x0, x0]
    y = [y0, y0, y1, y1, y0]
    plot!(x, y, color="black")
    r.global_bounds = (x0, x1, y0, y1)
    for attached_k = r.attached
        x0_k, y0_k = calc_offset(r, attached_k, x0, y0)
        plot_room!(attached_k.room, x0_k, y0_k)
    end
end

function calc_offset(r::Room, a::Attached_Type, x0::Float64, y0::Float64)
    if a.dir == "r"
        x0 += r.x
        y0 += a.Δ_pos
        door_pos = y0 + a.door_pos_in_child
        push!(r.nodes,      NODE(y=door_pos, dir_pair="r"))
        push!(a.room.nodes, NODE(y=door_pos, dir_pair="l"))
    elseif a.dir == "u"
        y0 += r.y
        x0 += a.Δ_pos
        door_pos = x0 + a.door_pos_in_child
        push!(r.nodes,      NODE(x=door_pos, dir_pair="u"))
        push!(a.room.nodes, NODE(x=door_pos, dir_pair="d"))
    elseif a.dir == "d"
        y0 -= a.room.y
        x0 += a.Δ_pos
        door_pos = x0 + a.door_pos_in_child
        push!(r.nodes,      NODE(x=door_pos, dir_pair="d"))
        push!(a.room.nodes, NODE(x=door_pos, dir_pair="u"))
    elseif a.dir == "l"
        x0 -= a.room.x
        y0 += a.Δ_pos
        door_pos = y0 + a.door_pos_in_child
        push!(r.nodes,      NODE(y=door_pos, dir_pair="l"))
        push!(a.room.nodes, NODE(y=door_pos, dir_pair="r"))
    else
        error("not implemented error")
    end

    return x0, y0
end

function do_start_smart(v_con::Nothing, r::Room, n::Int64, is_x::Bool)
    V0 = ifelse(is_x, get_x0, get_y0)(r)
    V1 = ifelse(is_x, get_x1, get_y1)(r)
    v_space = (V1 - V0) / n
    v_start = V0 + v_space / 2
    return v_start, v_space, true
end

function find_con_side(r::Room, v_con::Float64, is_x::Bool)
    v_con = ifelse(is_x, search_for_x_constraint, search_for_y_constraint)(r)
    V0 = ifelse(is_x, get_x0, get_y0)(r)
    V1 = ifelse(is_x, get_x1, get_y1)(r)
    return abs(v_con - V0) > abs(v_con - V1)
end

function do_start_smart(v_con::Float64, r::Room, n::Int64, is_x::Bool)
    (n == 1) && (return v_con, 0.0, true)
    V0 = ifelse(is_x, get_x0, get_y0)(r)
    V1 = ifelse(is_x, get_x1, get_y1)(r)
    is_pos = abs(v_con - V0) < abs(v_con - V1)
    v0 = v_con
    v1 = V1 + (V0 - v0)
    v_space = (v1 - v0) / (n - 1)
    if !is_pos
        v_space *= -1.0
        v0 = v1
    end
    return v0, v_space, is_pos
end

function do_start_smart(r::Room, is_x::Bool)
    n = ifelse(is_x, get_nx, get_ny)(r)
    v_con = ifelse(is_x, search_for_x_constraint, search_for_y_constraint)(r)
    v_start, v_space, is_pos = do_start_smart(v_con, r, n, is_x)
    return v_start, v_space, n, is_pos
end

function plot_node(r::Room, ::RefineTri)
    n1 = r.nodes[1]
    n2 = r.nodes[2]

    is_y(n1) || error("not implemented")
    is_y(n2) || error("not implemented")
    (length(r.nodes) == 2) || error("not implemented")

    X0 = r.global_bounds[1]
    X1 = r.global_bounds[2]
    ΔX = X1 - X0

    x0 = X0 + ΔX * 0.25
    x1 = X0 + ΔX * 0.75

    y0 = n1.y
    y1 = n2.y
    y1_below = y1 - 2*(y1 - y0)

    x = [x0, x1, x1]
    y = [y0, y1, y1_below]

    return x, y, true, true
end

function plot_node(r::Room, ::Union{RefineRectDR,RefineVert,RefineHoriz,RefineRectUR})
    start_x, space_x, nx, is_pos_x = do_start_smart(r, true)
    start_y, space_y, ny, is_pos_y = do_start_smart(r, false)
    x = Vector{Float64}()
    y = Vector{Float64}()
    for ix = 0:(nx-1)
        for iy = 0:(ny-1)
            push!(x, start_x + space_x * ix)
            push!(y, start_y + space_y * iy)
        end
    end
    return x, y, is_pos_x, is_pos_y
end

function plot_node!(r::Room)
    x, y, is_pos_x, is_pos_y = plot_node(r, r.refine_method)
    r.is_pos_x = is_pos_x
    r.is_pos_y = is_pos_y

    # prevent remove points from previous iterations if necessary
    index_in_v = findall(v .== r.id)
    if length(index_in_v) != 0
        index = only(index_in_v)
        delete!(scene, scene[end - index + 1])
        deleteat!(v, index)
    end
    scatter!(x, y)
    pushfirst!(v, r.id)

    for attached_k = r.attached # recurse
        is_pos_x, is_pos_y = plot_node!(attached_k.room)
    end

    return is_pos_x, is_pos_y
end

###

struct Graph
    points::Vector{SVector{2,Float64}}
    edges::Vector{SVector{2,Int64}}
    function Graph()
        points = Vector{SVector{2,Float64}}()
        edges = Vector{SVector{2,Int64}}()
        return new(points, edges)
    end
end

function add_graph_room!(g::Graph, r::Room, ::RefineTri, ee_parent)
    x, y = plot_node(r, RefineTri())

    len_orig = length(g.points)

    for k = 1:3
        push!(g.points, SVector{2,Float64}(x[k], y[k]))
    end

    for k = 1:3
        push!(g.edges, SVector{2,Int64}(k, rem(k, 3) + 1) .+ len_orig)
    end

    return ExtremeInd(1, 1, 3, 2, len_orig)
end

function add_graph_room!(g::Graph, r::Room, ::Union{RefineRectDR,RefineVert,RefineHoriz,RefineRectUR}, ee_parent)
    ind_2_node(ind_y::Int64, ind_x::Int64) = ind_y + (ind_x - 1) * ny

    start_x, space_x, nx, is_reg_x = do_start_smart(r, true)
    start_y, space_y, ny, is_reg_y = do_start_smart(r, false)
    point_matrix = Matrix{SVector{2,Float64}}(undef, ny, nx)
    edge_vec = Vector{SVector{2,Int64}}()
    for ix = 1:nx
        for iy = 1:ny
            x = start_x + space_x * (ix - 1)
            y = start_y + space_y * (iy - 1)
            point_matrix[iy, ix] = SVector{2,Float64}(x, y)
            index = ind_2_node(iy, ix)
            if ix != nx
                index_right = ind_2_node(iy, ix + 1)
                push!(edge_vec, SVector{2,Int64}(index, index_right))
            end
            if iy != ny
                index_down = ind_2_node(iy + 1, ix)
                push!(edge_vec, SVector{2,Int64}(index, index_down))
            end
        end
    end

    len_orig = length(g.points)

    for point_k = point_matrix
        # point_k = SVector{2,Float64}(point_k[1] + randn() * 0.15, point_k[2] + randn() * 0.15)
        push!(g.points, point_k)
    end

    for edge_k = edge_vec
        push!(g.edges, edge_k .+ len_orig)
    end

    return ExtremeInd(ind_2_node(ny, 1), ind_2_node(1, 1), ind_2_node(ny, nx), ind_2_node(1, nx), len_orig)
end

struct ExtremeInd
    i_ul::Int64
    i_dl::Int64
    i_ur::Int64
    i_dr::Int64
    function ExtremeInd(i_ul, i_dl, i_ur, i_dr, len_orig)
        return new(i_ul + len_orig, i_dl + len_orig, i_ur + len_orig, i_dr + len_orig)
    end
end

function add_graph!(g::Graph, r::Room, a=nothing)
    ee_parent = add_graph_room!(g, r, r.refine_method, a)

    for attached_k = r.attached
        ee_child = add_graph!(g, attached_k.room, attached_k)

        if attached_k.dir == "r"
            if r.is_pos_y
                i_parent = ee_parent.i_dr
            else
                i_parent = ee_parent.i_ur
            end
            i_child = ee_child.i_ul
            push!(g.edges, SVector{2,Int64}(i_parent, i_child))
        elseif attached_k.dir == "u"
            v_con = search_for_x_constraint(attached_k.room)
            i_child = ifelse(find_con_side(attached_k.room, v_con, true), ee_child.i_dr, ee_child.i_dl)
            i_parent = ifelse(find_con_side(r, v_con, true), ee_parent.i_ur, ee_parent.i_ul)
            push!(g.edges, SVector{2,Int64}(i_parent, i_child))
        elseif attached_k.dir == "d"
            v_con = search_for_x_constraint(attached_k.room)
            i_child = ifelse(find_con_side(attached_k.room, v_con, true), ee_child.i_ur, ee_child.i_ul)
            i_parent = ifelse(find_con_side(r, v_con, true), ee_parent.i_dr, ee_parent.i_dl)
            push!(g.edges, SVector{2,Int64}(i_parent, i_child))
        elseif attached_k.dir == "l"
            v_con = search_for_y_constraint(attached_k.room)
            i_child = ifelse(find_con_side(attached_k.room, v_con, false), ee_child.i_ur, ee_child.i_dr)
            i_parent = ifelse(find_con_side(r, v_con, false), ee_parent.i_ul, ee_parent.i_dl)
            push!(g.edges, SVector{2,Int64}(i_parent, i_child))
        end

    end

    return ee_parent
end

function print_graph(g::Graph)
    println("===================")
    println()

    dx = -6.95
    dy = 31.8

    # translate graph to my coordinates
    i = 1
    for point_k = g.points
        s0 = @sprintf("%d", i)
        s1 = @sprintf("%0.5f", point_k[1] + dx)
        s2 = @sprintf("%0.5f", point_k[2] + dy)
        println(s0," ", s1, " ", s2)
        i += 1
    end
    println()
    j = 1
    for edge_k = g.edges
        println(j," ",edge_k[1], " ", edge_k[2])
        j += 1
    end
    println()
end

function plot_graph(g::Graph)

    # # prevent remove points from previous iterations if necessary
    # index_in_v = findall(v .== -9999)
    # if length(index_in_v) != 0
    #     index = only(index_in_v)
    #     delete!(scene, scene[end - index + 1])
    #     deleteat!(v, index)
    # end

    for edge_k = g.edges
        p1 = g.points[edge_k[1]]
        p2 = g.points[edge_k[2]]
        plot!([p1[1], p2[1]], [p1[2], p2[2]])
    end

    # pushfirst!(v, -9999)
end

###

function create_layout_gazebo(is_debug::Bool=false)
    width = 3.9
    wid2 = 1.95

    # cafe
    r1 = Room(10.0, 17.0, RefineRectUR(), "room_1")

    is_debug && (r1.refine_method.nx = 3)
    is_debug && (r1.refine_method.ny = 5)

    # H3
    r2 = Room(10.0, 8.3, RefineTri(), "room_2")
    add_attached!(r1, r2, "r", 10.5, 3.3)

    # H2 - Stairs
    r3 = Room(12.0, width, RefineHoriz(), "room_3")
    add_attached!(r2, r3, "r", 0.0, wid2)

    is_debug && (r3.refine_method.n = 4)

    r_last_hall = r3

    for k = 1:4

        hx = ifelse(k == 4, 9.4, 9.3)
        r4 = Room(hx, width, RefineHoriz(), "room_4")
        add_attached!(r_last_hall, r4, "r", 0.0, wid2)

        r5 = Room(7.9, 5.5, RefineRectDR(), "room_5")
        add_attached!(r4, r5, "u", 0.0, 7.45)

        r_last_hall = r4

        is_debug && (r_last_hall.refine_method.n = 2)
        is_debug && (r5.refine_method.nx = 2)
        is_debug && (r5.refine_method.ny = 2)
    end

    width = 3.4
    for k = 1:3

        h = [4.1, 3.6, 4.8]
        wd = 1.5

        # H1
        r6 = Room(width, h[k], RefineVert(), "room_6")

        Δx = ifelse(k == 1, 6.0, 0.0)
        add_attached!(r_last_hall, r6, "d", Δx, wd)

        is_debug && (r6.refine_method.n = 2)

        # rooms A-B-C
        r7 = Room(5.83, 3.0, RefineHoriz(), "room_7")

        dh = ifelse(k==3, 1.2, 0.0)
        add_attached!(r6, r7, "l", dh, 1.5)

        is_debug && (r7.refine_method.n = 3)

        r_last_hall = r6
    end

    w8 = 17.0
    r8 = Room(w8, 31.0, RefineRectUR(), "room_8")
    add_attached!(r_last_hall, r8, "d", - 13.6, 15.8)

    is_debug && (r8.refine_method.nx = 3)
    is_debug && (r8.refine_method.ny = 5)

    return r1
end

function create_layout_bloat(is_debug::Bool=false)

    width = 4.29
    wid2 = width/2

    # cafe dimensions: (12.65, 23.2)
    r1 = Room(12.65, 23.2, RefineRectUR(), "room_1")

    is_debug && (r1.refine_method.nx = 3)
    is_debug && (r1.refine_method.ny = 5)

    # Hall 3 dimensions (11.0, 9.13)
    r2 = Room(11.0, 9.13, RefineTri(), "room_2")
    # Attach cafe-h3:  16.05, 2.4
    add_attached!(r1, r2, "r", 16.05, wid2)

    # H2.S dimensions (8.58, 4.29)
    r3 = Room(8.58, width, RefineHoriz(), "room_3")
    # H2.S attach  (0.0, 2.145)
    add_attached!(r2, r3, "r", 0.0, wid2)

    is_debug && (r3.refine_method.n = 4)

    r_last_hall = r3

    # rooms G-F-E-D (11.99, 4.29) (10.45, 4.29)
    wh2 = 4.29
    wh2h = wh2/2
    for k = 1:4

        # H2.D dimensions <> H2.G, H2.F, H2.E dimensions
        hx = ifelse(k == 4, 10.45, 11.99)
        r4 = Room(hx, wh2, RefineHoriz(), "room_4")
        # Attach classrooms: (0, 2.145)
        add_attached!(r_last_hall, r4, "r", 0.0, wh2h)

        # G-D rooms dimensions (10.45, 6.05)
        r5 = Room(10.45, 6.05, RefineRectDR(), "room_5")
        # G-D attach doors (0.0, 9.2)
        add_attached!(r4, r5, "u", 0.0, 9.2)

        r_last_hall = r4

        is_debug && (r_last_hall.refine_method.n = 2)
        is_debug && (r5.refine_method.nx = 2)
        is_debug && (r5.refine_method.ny = 2)
    end

    wh1 = 3.74
    wh1h = wh1/2

    for k = 1:3

        # H1C (3.74, 4.51) , H1B (3.74, 3.96), H1A (3.74, 5.28)
        h = [4.51, 3.96, 5.28]
        r6 = Room(width, h[k], RefineVert(), "room_6")

        # Attach H1-H2 (6.71, 1.87)
        Δx = ifelse(k == 1, 6.71, 0.0)
        add_attached!(r_last_hall, r6, "d", Δx, wh1h)

        is_debug && (r6.refine_method.n = 2)

        # Rooms A-B-C dimensions (6.41, 3.3)
        r7 = Room(6.41, 3.3, RefineHoriz(), "room_7")

        # H1B-C attach (0.0, 1.8), H1-A attach (1.32, 1.8)
        dh = ifelse(k==3, 1.32, 0.0)
        add_attached!(r6, r7, "l", dh, 1.8)

        is_debug && (r7.refine_method.n = 3)

        r_last_hall = r6
    end

    w8 = 18.7
    # Gym dimensions (18.7, 34.1)
    r8 = Room(w8, 34.1, RefineRectUR(), "room_8")
    # Attach Gym (-14.96, 16.83)
    add_attached!(r_last_hall, r8, "d", - 14.96, 16.83)

    is_debug && (r8.refine_method.nx = 3)
    is_debug && (r8.refine_method.ny = 5)

    return r1
end

scene = Scene()
r1 = create_layout_bloat()
plot_room!(r1)
plot_node!(r1)

### Event-based stuff

function print_active_room(r::Room)
    if r.is_active
        println(r.name, " is active!")
    end

    for attached_k = r.attached
        print_active_room(attached_k.room)
    end
end

function set_all_room_inactive(r::Room)
    r.is_active = false
    for attached_k = r.attached
        set_all_room_inactive(attached_k.room)
    end
end

function set_active_room(r::Room, x::Float32, y::Float32)
    if (r.global_bounds[1] < x < r.global_bounds[2])
        if (r.global_bounds[3] < y < r.global_bounds[4])
            r.is_active = true
        end
    end
    for attached_k = r.attached
        set_active_room(attached_k.room, x, y)
    end
end

function refine_active!(r::Room, fun_name::Function)
    r.is_active && fun_name(r.refine_method)
    for attached_k = r.attached
        refine_active!(attached_k.room, fun_name)
    end
end

display(scene)

xlims!(-40, 75.0)
ylims!(-40, 40.0)

if true
    dir = lift(scene.events.keyboardbuttons) do but
        is_lr = ispressed(but, Keyboard.left) || ispressed(but, Keyboard.right)
        is_ud = ispressed(but, Keyboard.up) || ispressed(but, Keyboard.down)

        if is_lr
            println("in keyboard event loop x")
            refine_active!(r1, refine_x!)
            plot_node!(r1)
        end

        if is_ud
            println("in keyboard event loop y")
            refine_active!(r1, refine_y!)
            plot_node!(r1)
        end

        if ispressed(but, Keyboard.space)
            g = Graph()
            add_graph!(g, r1)
            plot_graph(g)
            print_graph(g)
        end
        return 0
    end

    pos = lift(scene.events.mouseposition) do mpos
        # println("HIT")

        x0, x1 = mouseposition(scene)

        set_all_room_inactive(r1)
        set_active_room(r1, x0, x1)
        return 0
    end
end

println()
println()
println("DONE!!!")


text:
using Base.Threads
import Base.Threads.@spawn

"""Merge two `CellList{d}`s with same dimension `d`.

# Examples
```julia
n, d, r = 100, 2, 0.01
p = (rand(n, d), rand(n, d))
merge(CellList(p[1], r), CellList(p[2], r))
```
"""
function Base.merge(c1::CellList{d}, c2::CellList{d}) where d
    CellList{d}(merge(vcat, c1.data, c2.data))
end

"""Multithreaded CellList constructor.

# Examples
```julia
n, d, r = 100, 2, 0.01
p = rand(n, d)
c = CellList(p, r, Val(:threads))
```
"""
function CellList(p::Array{T, 2}, r::T, ::Val{:threads}) where T <: AbstractFloat
    t = nthreads()
    n, d = size(p)
    cs = cumsum(fill(fld(n, t), t-1))
    parts = zip([0; cs], [cs; n])
    res = Array{CellList{d}, 1}(undef, t)
    @threads for (i, (a, b)) in collect(enumerate(parts))
        res[i] = CellList(p[(a+1):b, :], r; offset=a)
    end
    return reduce(merge, res)
end

function brute_force_cell!(pts, cell, is, p, r, data, offsets)
    # Pairs of points within the cell
    brute_force!(pts[threadid()], is, p, r)
    # Pairs of points with non-empty neighboring cells
    for offset in offsets
        neigh_cell = cell + offset
        if haskey(data, neigh_cell)
            @inbounds js = data[neigh_cell]
            brute_force!(pts[threadid()], is, js, p, r)
        end
    end
end

"""Multithreaded near neighbors.

# Examples
```julia-repl
julia> near_neighbors(c, p, r, Val(:threads))
[(1, 4), (3, 11), ...]
```
"""
function near_neighbors(c::CellList{d}, p::Array{T, 2}, r::T, ::Val{:threads}) where d where T <: AbstractFloat
    offsets = neighbors(d)
    pts = [Vector{Tuple{Int, Int}}() for _ in 1:nthreads()]
    data = collect(c.data)
    perm = sortperm(@. length(getfield(data, 2)))
    # Iterate over non-empty cells in decreasing order of points per cell
    @sync for (cell, is) in reverse(data[perm])
        @spawn brute_force_cell!(pts, cell, is, p, r, c.data, offsets)
    end
    return reduce(vcat, pts)
end


text:
export ISMRMRDFile, saveasISMRMRDFile

include("HDF5Types.jl")

struct ISMRMRDFile
  filename::String
end

"""
    RawAcquisitionData(f::ISMRMRDFile, dataset="dataset")

reads the `ISMRMRDFile` f and stores the result in a `RawAcquisitionDataObject`
"""
function RawAcquisitionData(f::ISMRMRDFile, dataset="dataset")

  h5open(f.filename) do h
    headerStr = read(h["/$(dataset)/xml"])
    params = GeneralParameters(headerStr[1])

    d = read(h["/$(dataset)/data"])
    dtypehead = h5t_get_member_type(datatype(h["/$(dataset)/data"]),0)

    M = length(d)

    profiles = Profile[]

    for m=1:M
      head = read_header(d[m].data[1], dtypehead)
      D = Int(head.trajectory_dimensions)
      chan = Int(head.active_channels)

      N = reinterpret(Int64, d[m].data[2][1:8])[1]
      if N > 0
        ptr = reinterpret(Int64, d[m].data[2][9:end])[1]
        traj = unsafe_wrap(Array,Ptr{Float32}(ptr),(D,div(N,D)))
      else
        traj = Matrix{Float32}(undef,0,0)
      end

      N = reinterpret(Int64, d[m].data[3][1:8])[1]
      if N > 0
        ptr = reinterpret(Int64, d[m].data[3][9:end])[1]
        U = unsafe_wrap(Array,Ptr{ComplexF32}(ptr),(div(N,2),))
        dat = reshape(U,div(N,2*chan),chan)
      else
        dat = Matrix{ComplexF32}(undef,0,0)
      end

      push!(profiles, Profile(head,traj,dat) )
    end
    return RawAcquisitionData(params, profiles)
  end
end

"""
    AcquisitionData(f::ISMRMRDFile, dataset="dataset")

reads the `ISMRMRDFile` f and stores the result in an `AcquisitionDataObject`
"""
function AcquisitionData(f::ISMRMRDFile, dataset="dataset")
  return AcquisitionData(RawAcquisitionData(f,dataset))
end

function read_header(header, dtypehead)
  # very important to use the offset here, since the offset varies from
  # file to file

  buf = IOBuffer(header)

  seek(buf, h5t_get_member_offset(dtypehead, 0))
  version = read(buf,Int16)
  seek(buf, h5t_get_member_offset(dtypehead, 1))
  flags = read(buf,Int64)
  seek(buf, h5t_get_member_offset(dtypehead, 2))
  measurement_uid = read(buf,Int32)
  seek(buf, h5t_get_member_offset(dtypehead, 3))
  scan_counter = read(buf,Int32)
  seek(buf, h5t_get_member_offset(dtypehead, 4))
  acquisition_time_stamp = read(buf,Int32)
  seek(buf, h5t_get_member_offset(dtypehead, 5))
  physiology_time_stamp = ntuple(i->read(buf,Int32),3)
  seek(buf, h5t_get_member_offset(dtypehead, 6))
  number_of_samples = read(buf,Int16)
  seek(buf, h5t_get_member_offset(dtypehead, 7))
  available_channels = read(buf,Int16)
  seek(buf, h5t_get_member_offset(dtypehead, 8))
  active_channels = read(buf,Int16)
  seek(buf, h5t_get_member_offset(dtypehead, 9))
  channel_mask = ntuple(i->read(buf,Int64),16)
  seek(buf, h5t_get_member_offset(dtypehead, 10))
  discard_pre = read(buf,Int16)
  seek(buf, h5t_get_member_offset(dtypehead, 11))
  discard_post = read(buf,Int16)
  seek(buf, h5t_get_member_offset(dtypehead, 12))
  center_sample = read(buf,Int16)
  seek(buf, h5t_get_member_offset(dtypehead, 13))
  encoding_space_ref = read(buf,Int16)
  seek(buf, h5t_get_member_offset(dtypehead, 14))
  trajectory_dimensions = read(buf,Int16)
  seek(buf, h5t_get_member_offset(dtypehead, 15))
  sample_time_us = read(buf,Float32)
  seek(buf, h5t_get_member_offset(dtypehead, 16))
  position = ntuple(i->read(buf,Float32),3)
  seek(buf, h5t_get_member_offset(dtypehead, 17))
  read_dir = ntuple(i->read(buf,Float32),3)
  seek(buf, h5t_get_member_offset(dtypehead, 18))
  phase_dir = ntuple(i->read(buf,Float32),3)
  seek(buf, h5t_get_member_offset(dtypehead, 19))
  slice_dir = ntuple(i->read(buf,Float32),3)
  seek(buf, h5t_get_member_offset(dtypehead, 20))
  patient_table_position = ntuple(i->read(buf,Float32),3)
  ###
  dtype_enc = h5t_get_member_type(dtypehead,21)
  off = h5t_get_member_offset(dtypehead, 21)
  seek(buf, h5t_get_member_offset(dtype_enc, 0) + off)
  kspace_encode_step_1 = read(buf,Int16)
  seek(buf, h5t_get_member_offset(dtype_enc, 1) + off)
  kspace_encode_step_2 = read(buf,Int16)
  seek(buf, h5t_get_member_offset(dtype_enc, 2) + off)
  average = read(buf,Int16)
  seek(buf, h5t_get_member_offset(dtype_enc, 3) + off)
  slice = read(buf,Int16)
  seek(buf, h5t_get_member_offset(dtype_enc, 4) + off)
  contrast = read(buf,Int16)
  seek(buf, h5t_get_member_offset(dtype_enc, 5) + off)
  phase = read(buf,Int16)
  seek(buf, h5t_get_member_offset(dtype_enc, 6) + off)
  repetition = read(buf,Int16)
  seek(buf, h5t_get_member_offset(dtype_enc, 7) + off)
  set = read(buf,Int16)
  seek(buf, h5t_get_member_offset(dtype_enc, 8) + off)
  segment = read(buf,Int16)
  seek(buf, h5t_get_member_offset(dtype_enc, 9) + off)
  user = ntuple(i->read(buf,Int16),8)
  ###
  idx = EncodingCounters(kspace_encode_step_1,kspace_encode_step_2,
                      average,slice,contrast,phase,repetition,set,segment,user)

  seek(buf, h5t_get_member_offset(dtypehead, 22))
  user_int = ntuple(i->read(buf,Int32),8)
  seek(buf, h5t_get_member_offset(dtypehead, 23))
  user_float = ntuple(i->read(buf,Float32),8)

  head = AcquisitionHeader( version, flags, measurement_uid,
    scan_counter, acquisition_time_stamp, physiology_time_stamp, number_of_samples,
    available_channels, active_channels, channel_mask, discard_pre, discard_post,
    center_sample, encoding_space_ref, trajectory_dimensions,sample_time_us,
    position, read_dir, phase_dir, slice_dir, patient_table_position, idx, user_int,
    user_float
    )

  return head
end

function FileIO.save(f::ISMRMRDFile, acq::RawAcquisitionData, dataset="dataset")
  h5open(f.filename, "w") do file
    headerStr = GeneralParametersToXML(acq.params)
    write(file, "/$(dataset)/xml", [headerStr]) # [] ensures that we store as var string
    writeProfiles(file, "/$(dataset)/data", acq.profiles)
  end
end


text:
"""
    Economics(
        E0, γ, β, ρ, Finf,
        mitigate_cost, remove_cost, geoeng_cost, adapt_cost,
        mitigate_init, remove_init, geoeng_init, adapt_init,
        baseline_emissions,
        extra_CO₂
    )

Create data structure for economic input parameters for `ClimateModel` struct,
including a baseline emissions scenario.

### Arguments
- `E0::Array{Float64,1}`: Gross World Product timeseries [10^12 USD / year]
- `γ::Float64`: economic growth rate [fraction]
- `β::Float64`: climate damage parameter [% GWP / (°C)^2].
- `ρ::Float64`: typically denoted ρ in economic references [fraction].
- `Finf::Float64`: maximum SRM forcing
- `[control]_cost::Float64`: scaling cost of full control deployment [10^12 USD / year OR % of GWP].
- `[control]_init::Float64`: fixed initial condition for control deployment [10^12 USD / year].
- `baseline_emissions::Array{Float64,1}`: prescribed baseline CO₂ equivalent emissions [ppm / yr].
- `extra_CO₂::Array{Float64,1}`: optional additional CO₂ input that is used only for Social Cost of Carbon calculations [ppm].

See also: [`ClimateModel`](@ref), [`baseline_emissions`](@ref), [`GWP`](@ref).

"""
mutable struct Economics
    E0::Float64
    γ::Float64
    β::Float64
    ρ::Float64
    Finf::Float64

    mitigate_cost::Float64
    remove_cost::Float64
    geoeng_cost::Float64
    adapt_cost::Float64

    mitigate_init
    remove_init
    geoeng_init
    adapt_init

    baseline_emissions::Array{Float64,1}
    extra_CO₂::Array{Float64,1}
end

# constructor without `extra_CO₂`, setting it to zero
function Economics(E0::Float64,
    γ::Float64,
    β::Float64,
    ρ::Float64,
    Finf::Float64,
    mitigate_cost::Float64,
    remove_cost::Float64,
    geoeng_cost::Float64,
    adapt_cost::Float64,
    mitigate_init,
    remove_init,
    geoeng_init,
    adapt_init,
    baseline_emissions::Array{Float64,1})
    return Economics(
        E0,
        γ,
        β,
        ρ,
        Finf,
        mitigate_cost,
        remove_cost,
        geoeng_cost,
        adapt_cost,
        mitigate_init,
        remove_init,
        geoeng_init,
        adapt_init,
        baseline_emissions,
        zeros(size(baseline_emissions))
    )
end


text:
using Combinatorics

part1(input) = find_best_distance(input, min)
part2(input) = find_best_distance(input, max)

function find_best_distance(input, order)
    distances = parse_line.(eachline(input))
    places = unique(foldl(vcat, first.(distances)))
    n = length(places)
    ids = Dict(place => i for (i, place) in enumerate(places))
    D = fill(0, n, n)
    for ((place1, place2), distance) in distances
        i, j = ids[place1], ids[place2]
        D[i, j] = distance
        D[j, i] = distance
    end
    return foldl(order, (route_length(D, route) for route in permutations(1:n)))
end

function parse_line(line)
    parts = split(line)
    return parts[[1, 3]], parse(Int, parts[5])
end

route_length(D, route) = sum(D[route[i - 1], route[i]] for i in 2:length(route))


text:
##########
# Parser #
##########

mutable struct Parser
    l::Lexer
    current_token::Token
    Parser(str) = new(Lexer(str))
end

current_token(ps::Parser) = ps.current_token
next_token!(ps::Parser) = return (ps.current_token = gettok(ps.l))

# Operator precedence
const BinopPrecedence = Dict{Kinds.Kind, Int}()
BinopPrecedence[Kinds.EQUAL]   = 2
BinopPrecedence[Kinds.LESS]    = 10
BinopPrecedence[Kinds.GREATER] = 10
BinopPrecedence[Kinds.PLUS]    = 20
BinopPrecedence[Kinds.MINUS]   = 20
BinopPrecedence[Kinds.STAR]    = 40
BinopPrecedence[Kinds.SLASH]   = 40

function operator_precedence(ps)
    v = current_token(ps).kind
    return (v in keys(BinopPrecedence)) ? BinopPrecedence[v] : -1
end

#############
# AST Nodes #
#############

abstract type ExprAST end

struct NumberExprAST <: ExprAST
    val::Float64
end
Base.show(io::IO, expr::NumberExprAST) = print(io, expr.val)

struct VariableExprAST <: ExprAST
    name::String
end
Base.show(io::IO, expr::VariableExprAST) = print(io, expr.name)

struct BinaryExprAST <: ExprAST
    op::Kinds.Kind
    lhs::ExprAST
    rhs::ExprAST
end
Base.show(io::IO, expr::BinaryExprAST) = print(io, expr.op, "(", expr.lhs, ", ", expr.rhs, ")")

struct CallExprAST <: ExprAST
    callee::String
    args::Vector{ExprAST}
end

struct IfExprAST <: ExprAST
    cond::ExprAST
    then::ExprAST
    elsee::ExprAST
end

struct ForExprAST <: ExprAST
    varname::String
    start::ExprAST
    endd::ExprAST
    step::ExprAST
    body::ExprAST
end

struct VarExprAST <: ExprAST
    varnames::Vector{Tuple{String, ExprAST}}
end

struct BlockExprAST <: ExprAST
    exprs::Vector{ExprAST}
end

struct PrototypeAST
    name::String
    args::Vector{String}
end

struct FunctionAST
    proto::PrototypeAST
    body::ExprAST
end

#####################
# Parse Expressions #
#####################

function ParseNumberExpr(ps::Parser)::NumberExprAST
    result = NumberExprAST(Base.parse(Float64, current_token(ps).val))
    return result
end

function ParseIdentifierExpr(ps::Parser)::Union{ VariableExprAST, CallExprAST}
    idname = current_token(ps).val
    next_token!(ps) # eat the idname

    if current_token(ps).kind != Kinds.LPAR
        return VariableExprAST(idname)
    end

    next_token!(ps) # eat '('
    args = ExprAST[]
    first = true
    while true
        if current_token(ps).kind == Kinds.RPAR
            break
        end
        if !first
            if current_token(ps).kind != Kinds.COMMA
                error("Expected ')' or ',' in argument list, got $(current_token(ps))")
            end
            next_token!(ps) # eat the ','
        end
        first = false
        push!(args, ParseExpression(ps))
    end
    next_token!(ps) # eat ')'
    return CallExprAST(idname, args)
end

function ParseIfExpr(ps)::IfExprAST
    # if
    next_token!(ps) # eat 'if'
    cond = ParseExpression(ps)

    # then
    if current_token(ps).kind != Kinds.THEN
        error("expected 'then'")
    end
    next_token!(ps) # eat 'then'
    then = ParseExpression(ps)

    # else
    if current_token(ps).kind != Kinds.ELSE
        error("expected 'else'")
    end
    next_token!(ps)
    elsee = ParseExpression(ps)

    return IfExprAST(cond, then, elsee)
end

function ParsePrototype(ps)::PrototypeAST
    if current_token(ps).kind != Kinds.IDENTIFIER
        error("Expected function name in prototype, got $(current_token(ps))")
    end

    func_name = current_token(ps).val
    tok = next_token!(ps) # eat identifier

    if tok.kind != Kinds.LPAR
        error("Expected '(' in prototype")
    end

    argnames = String[]
    while (next_token!(ps).kind == Kinds.IDENTIFIER)
        push!(argnames, current_token(ps).val)
    end

    if current_token(ps).kind != Kinds.RPAR
        error("Expected ')' in prototype, got $(current_token(ps))")
    end

    next_token!(ps)

    return PrototypeAST(func_name, argnames)
end

function ParseDefinition(ps)::FunctionAST
    next_token!(ps) # eat def
    proto = ParsePrototype(ps)
    E = ParseExpression(ps)
    return FunctionAST(proto, E)
end

function ParseParenExpr(ps::Parser)::ExprAST
    next_token!(ps) # eat '('
    V = ParseExpression(ps)
    if current_token(ps).kind != Kinds.RPAR
        error("expected ')'")
    end
    next_token!(ps) # eat ')'
    return V
end

function ParseBinOpRHS(ps, ExprPrec::Int, LHS::ExprAST)::ExprAST
    while true
        tokprec = operator_precedence(ps)
        if tokprec < ExprPrec
            return LHS
        end

        bin_op = current_token(ps)
        next_token!(ps) # eat binary token

        RHS = ParsePrimary(ps)
        nextprec = operator_precedence(ps)
        if tokprec < nextprec
            RHS = ParseBinOpRHS(ps, tokprec + 1, RHS)
        end

        LHS = BinaryExprAST(bin_op.kind, LHS, RHS)
    end
end

function ParseForExpr(ps)::ForExprAST
    next_token!(ps) # eat 'for'

    if current_token(ps).kind != Kinds.IDENTIFIER
        error("expected identifier after for")
    end

    idname = current_token(ps).val
    next_token!(ps) # eat identifier

    if current_token(ps).kind != Kinds.EQUAL
        error("expected `=` after identifier in for expression")
    end
    next_token!(ps) # eat =

    start = ParseExpression(ps)
    if current_token(ps).kind != Kinds.COMMA
        error("expected `,` after for start value")
    end
    next_token!(ps) # eat ,

    endd = ParseExpression(ps)

    # TODO: make optional
    if current_token(ps).kind != Kinds.COMMA
        error("expected ',' after for end value")
    end
    next_token!(ps) # eat ,
    step = ParseExpression(ps)

    if current_token(ps).kind != Kinds.IN
        error("expected 'in' after for")
    end
    next_token!(ps) # eat in

    body = ParseExpression(ps)

    return ForExprAST(idname, start, endd, step, body)
end

function ParseExtern(ps)::PrototypeAST
    next_token!(ps) # eat 'extern'
    return ParsePrototype(ps)
end

function ParseVarExpr(ps)
    next_token!(ps) #eat the var
    varnames = Tuple{String, ExprAST}[]
    if current_token(ps).kind != Kinds.IDENTIFIER
        error("expected identifier after var")
    end

    while true
        name = current_token(ps).val
        next_token!(ps) # eat the identifier
        # TODO: Optional initializer
        if current_token(ps).kind != Kinds.EQUAL
            error("expected equal after var identifier")
        end
        next_token!(ps) # eat the '='
        init = ParseExpression(ps)

        push!(varnames, (name, init))

        current_token(ps).kind != Kinds.COMMA && break
        next_token!(ps) # eat the '.'
        if current_token(ps).kind != Kinds.IDENTIFIER
            error("expected identifier list after var")
        end
    end

    return VarExprAST(varnames)
end

function ParseTopLevelExpr(ps)::FunctionAST
    E = ParseExpression(ps)
    proto = PrototypeAST("__anon_expr", String[])
    return FunctionAST(proto, E)
end

@noinline function ParseExpression(ps)::ExprAST
    LHS = ParsePrimary(ps)
    return ParseBinOpRHS(ps, 0, LHS)
end

function ParseBlockExpr(ps)::BlockExprAST
    next_token!(ps) # eat the '{'
    exprs = ExprAST[]
    while true
        if current_token(ps).kind == Kinds.RBRACE
            next_token!(ps) # eat the '}'
            break
        end
        push!(exprs, ParseExpression(ps))
    end
    return BlockExprAST(exprs)
end

function ParsePrimary(ps)::ExprAST
    curtok = current_token(ps)
    if curtok.kind == Kinds.IDENTIFIER
        return ParseIdentifierExpr(ps)
    elseif curtok.kind == Kinds.NUMBER
        return ParseNumberExpr(ps)
    elseif curtok.kind == Kinds.LPAR
        return ParseParenExpr(ps)
    elseif curtok.kind == Kinds.IF
        return ParseIfExpr(ps)
    elseif curtok.kind == Kinds.FOR
        return ParseForExpr(ps)
    elseif curtok.kind == Kinds.VAR
        return ParseVarExpr(ps)
    elseif curtok.kind == Kinds.LBRACE
        return ParseBlockExpr(ps)
    else
        error("unexpected token: $curtok")
    end
end


text:
"""
    @floop begin
        s₁ = initialization of s₁
        s₂  # pre-initialized variable
        ...
        for x in xs, ...
            ...
        end
    end

`@floop begin ... end` expects a (possibly empty) series of
assignments or variable declaration (as in `s₂` above) followed by a
`for` loop.

When there is no induction variables, `begin ... end` can be omitted:

    @floop for x in xs, ...
        ...
    end

Use [`@reduce`](@ref) for parallel execution:

    @floop for x in xs, ...
        ...
        @reduce ...
    end

`@floop` can also take an `executor` argument (which should be an instance of
executor such as `SequentialEx`, `ThreadedEx` and `DistributedEx`) or a
`nothing` (indicating an appropriate parallel executor):

    @floop executor for x in xs, ...
        ...
        @reduce ...
    end

See the module docstring of [`FLoops`](@ref) for examples.
"""
macro floop(ex)
    ctx = MacroContext(__source__, __module__)
    ex, simd = remove_at_simd(__module__, ex)
    exx = macroexpand(__module__, ex)
    is_parallel(exx) && return esc(floop_parallel(ctx, exx, simd))
    esc(floop(exx, simd))
end

macro floop(executor, ex)
    ctx = MacroContext(__source__, __module__)
    ex, simd = remove_at_simd(__module__, ex)
    exx = macroexpand(__module__, ex)
    esc(floop_parallel(ctx, exx, simd, executor))
end

struct Return{T}
    value::T
end

struct Goto{label,T}
    acc::T
end

Goto{label}(acc::T) where {label,T} = Goto{label,T}(acc)
gotoexpr(label::Symbol) = :($Goto{$(QuoteNode(label))})

function floop(ex, simd)
    pre = post = Union{}[]
    ansvar = :_
    if isexpr(ex, :for)
        loops, body = _for_loop(ex)
    elseif isexpr(ex, :block)
        args = flattenblockargs(ex)
        i, ansvar, loops, body =
            find_first_for_loop(args) |> ifnothing() do
                throw(ArgumentError("Unsupported expression:\n$ex"))
            end
        pre = args[1:i-1]
        post = args[i+1:end]
        if find_first_for_loop(post) !== nothing
            throw(ArgumentError("Multiple top-level `for` loop found in:\n$ex"))
        end
    else
        throw(ArgumentError("Unsupported expression:\n$ex"))
    end
    if ansvar !== :_
        post = vcat(post, ansvar)
    end
    pre = vcat(pre, something(EXTRA_STATE_VARIABLES[], Union{}[]))

    init_vars = mapcat(assigned_vars, pre)
    foldlex = @match loops begin
        Expr(:block, loop_axes...) => begin
            rf_arg, coll = transform_multi_loop(loop_axes)
            asfoldl(rf_arg, coll, body, init_vars, simd)
        end
        Expr(:(=), rf_arg, coll) => begin
            asfoldl(rf_arg, coll, body, init_vars, simd)
        end
    end
    return Expr(:block, pre..., :($ansvar = $foldlex), post...)
end

function _for_loop(ex)
    if isexpr(ex, :for) && length(ex.args) == 2
        return ex.args
    end
    throw(ArgumentError("Malformed `for` block:\n$ex"))
end

find_first_for_loop(args) =
    firstsomething(enumerate(args)) do (i, x)
        @match x begin
            Expr(:for, loops, body) => return (i, :_, loops, body)
            Expr(:local, Expr(:(=), ansvar, Expr(:for, loops, body))) =>
                return (i, ansvar, loops, body)
            _ => nothing
        end
    end

function gotos_for(external_labels::Vector{Symbol}, unpack_state::Expr, acc::Symbol)
    gotos = map(external_labels) do label
        quote
            if $acc isa $(gotoexpr(label))
                let $acc = $acc.acc
                    $unpack_state
                end
                $Base.@goto $label
            end
        end
    end
    return gotos
end

function asfoldl(rf_arg, coll, body0, state_vars, simd)
    @assert simd in (false, true, :ivdep)
    body, info = transform_loop_body(body0, state_vars)
    @gensym step acc xf foldable
    pack_state = info.pack_state
    unpack_state = :(($(state_vars...),) = $acc)
    gotos = gotos_for(info.external_labels, unpack_state, acc)
    state_declarations = [:(local $v) for v in state_vars]
    quote
        @inline function $step($acc, $rf_arg)
            $(state_declarations...)
            $unpack_state
            $body
            return $pack_state
        end
        $xf, $foldable = $extract_transducer($coll)
        $acc = $foldl($step, $xf, $foldable; init = $pack_state, simd = $(Val(simd)))
        $acc isa $Return && return $acc.value
        $(gotos...)
        $unpack_state
        nothing
    end
end
# For now, using (private interface) `extract_transducer` as done in:
# https://github.com/tkf/ThreadsX.jl/pull/106.  But this should be
# done automatically in Transducers.jl.

function transform_loop_body(body, state_vars)
    external_labels::Vector{Symbol} = setdiff(gotos_in(body), labels_in(body))
    # state_vars = extract_state_vars(body)
    pack_state = :(($(state_vars...),))
    info = (
        state_vars = state_vars,
        pack_state = pack_state,
        external_labels = external_labels,
        nested_for = false,
    )
    return as_rf_body(body, info), info
end

function as_rf_body(body, info)
    is_function(body) && return body
    @match body begin
        # Do not transform `break` in inside other `for` loops:
        Expr(:for, nloop, nbody) =>
            Expr(:for, nloop, as_rf_body(nbody, @set info.nested_for = true))

        Expr(:continue) => info.nested_for ? body : :(return $(info.pack_state))
        Expr(:break) => info.nested_for ? body : :(return $reduced($(info.pack_state)))
        Expr(:return) => :(return $reduced($Return(nothing)))
        Expr(:return, value) => :(return $reduced($Return($value)))
        Expr(:symbolicgoto, label) => if label in info.external_labels
            :(return $reduced($(gotoexpr(label))($(info.pack_state))))
        else
            body
        end

        Expr(:let, let_bindings_, let_body) => begin
            let_bindings = @match let_bindings_ begin
                Expr(:block, args...) => collect(args)
                b => [b]
            end
            let_vars = map(let_bindings) do x
                @match x begin
                    _::Symbol => x
                    Expr(:(=), x, _) => x
                end
            end
            frozen_vars = intersect(info.state_vars, let_vars)
            gensym_vars = map(gensym, frozen_vars)
            d = Dict(zip(frozen_vars, gensym_vars))
            pack_state_args = @match info.pack_state begin
                :(($(args...),)) => begin
                    map(a -> get(d, a, a), args)
                end
            end
            info = @set info.state_vars = setdiff(info.state_vars, frozen_vars)
            info = @set info.pack_state = :(($(pack_state_args...),))
            all_bindings = map(gensym_vars, frozen_vars) do g, f
                :($g = $f)
            end
            append!(all_bindings, let_bindings)
            Expr(:let, Expr(:block, all_bindings...), as_rf_body(let_body, info))
        end

        Expr(head, args...) => Expr(head, [as_rf_body(a, info) for a in args]...)

        _ => body
    end
end

assigned_vars(x::Symbol) = [x]
assigned_vars(::Any) = Symbol[]
function assigned_vars(ex::Expr)
    @match ex begin
        Expr(:(=), lhs, rhs) => vcat(vars_in(lhs), assigned_vars(rhs))
        Expr(:tuple, lhs1..., Expr(:(=), lhs2, rhs)) =>
            vcat(mapfoldl(vars_in, vcat, lhs1), vars_in(lhs2), assigned_vars(rhs))
        # Is it better to be permissive here?
        _ => Symbol[]
    end
end

vars_in(x::Symbol) = [x]
function vars_in(ex)
    @match ex begin
        Expr(:tuple, vars...) => vars
        _ => Symbol[]
    end
end

# extract_state_vars(body) =
#     mutated_global_vars(solve!(simplify_ex(deepcopy(body))))

# function mutated_global_vars(ex)
#     return @match ex begin
#         Expr(:(=), lhs, rhs) => begin
#             if lhs isa Var && lhs.is_global
#                 append!([lhs.name], mutated_global_vars(rhs))
#             else
#                 mutated_global_vars(rhs)
#             end
#         end
#         Expr(_, args...) => mapreduce(mutated_global_vars, vcat, args)
#         _ => Union{}[]
#     end
# end

gotos_in(_) = Symbol[]
function gotos_in(ex::Expr)::Vector{Symbol}
    @match ex begin
        Expr(:symbolicgoto, label) => [label]
        Expr(_, args...) => mapcat(gotos_in, args)
    end
end

labels_in(_) = Symbol[]
function labels_in(ex::Expr)::Vector{Symbol}
    @match ex begin
        Expr(:symboliclabel, label) => [label]
        Expr(_, args...) => mapcat(labels_in, args)
    end
end

unbound_rhs(@nospecialize _) = Symbol[]
unbound_rhs(x::Symbol) = [x]
unbound_rhs(ex::Expr) = _global_rhs(solve!(simplify_ex(deepcopy(ex))))
_global_rhs(_) = Symbol[]
_global_rhs(x::Var) = x.is_global ? [x.name] : Symbol[]
function _global_rhs(ex::Expr)
    @match ex begin
        Expr(:(=), _, rhs) => _global_rhs(rhs)
        Expr(:tuple, _..., Expr(:(=), _, rhs)) => _global_rhs(rhs)
        Expr(_, args...) => mapfoldl(_global_rhs, append!, args; init = Symbol[])
    end
end

function transform_multi_loop(loop_axes)
    all_loop_vars = Symbol[]
    loop_vars = []
    loop_collections = []
    is_triangular = false
    for ex in loop_axes
        @match ex begin
            Expr(:(=), var, axis) => begin
                is_triangular |= !isempty(intersect(all_loop_vars, unbound_rhs(axis)))
                append!(all_loop_vars, vars_in(var))
                push!(loop_vars, var)
                push!(loop_collections, axis)
            end
        end
    end
    rf_arg = :(($(reverse(loop_vars)...),))
    if is_triangular
        coll = foldr(
            collect(zip(loop_vars, loop_collections));  # `collect` for Julia 1.0
            init = :($(reverse(loop_vars)...),),
        ) do (v, axis), coll
            f = gensym("loop_axis_$v")
            Expr(:call, |>, axis, quote
                $f($v) = $coll  # should it be inlined?
                $Map($f)
            end)
        end
        for _ in 2:length(loop_collections)
            coll = Expr(:call, |>, coll, :($Cat()))
        end
        return (rf_arg, coll)
    else
        return (rf_arg, :($Iterators.product($(reverse(loop_collections)...))))
    end
end
#
# Notes on triangular loops: For example, a `for` loop specification
# `for x in xs, y in f(x), z in g(x, y)` is transformed to
#
#     xs |> Map() do x
#         f(x) |> Map() do y
#             g(x, y) |> Map() do z
#                 (z, y, x)
#             end
#         end
#     end |> Cat() |> Cat()
#
# `Cat`s are composed outside so that `foldxt` can (in principle)
# rewrite the transducer to use `TCat` to parallelize inner loops.

function resolvesymbol(m, e)
    @match e begin
        Expr(:., a, QuoteNode(b::Symbol)) => begin
            x = resolvesymbol(m, a)
            x isa Module ? getfield(x, b) : nothing
        end
        (a::Symbol) => getfield(m, a)
        _ => nothing
    end
end

"""
    remove_at_simd(__module__, ex) -> (ex′, simd)

Return a tuple with following items:

* `ex′`: `ex` without the first top-level macrocall to `Base.@simd` removed.
* `simd`: `true` if `@simd for ... end` is found. `:ivdep` if `@simd
  ivdep for ... end` is found.  `false` otherwise.

Macros that happened to have the name `@simd` but not identical to
`Base.@simd` are ignored.
"""
function remove_at_simd(__module__, ex)
    ans = @match ex begin
        Expr(:macrocall, mcr, ::LineNumberNode, options..., loop) => begin
            if (
                isexpr(loop, :for) &&
                resolvesymbol(__module__, mcr) === getfield(Base, Symbol("@simd"))
            )
                if options == [:ivdep]
                    return loop, :ivdep
                elseif options == []
                    return loop, true
                end
            end
        end
        Expr(:block, args...) => begin
            for (i, x) in enumerate(args)
                y, simd = remove_at_simd(__module__, x)
                if simd !== false
                    return Expr(:block, args[1:i-1]..., y, args[i+1:end]...), simd
                end
            end
        end
        _ => nothing
    end
    ans === nothing || return ans
    return ex, false
end