text:
# This code is my Julia version of the taut string 1d solver of Barbero and Sra
# The original code of the authors is in https://github.com/albarji/proxTV/blob/master/src/TVL1opt_tautstring.cpp
# Copyright (c) 2013, Álvaro Barbero Jiménez and Suvrit Sra. All rights reserved.

import Base.getindex
import Base.setindex!
import Base.show

mutable struct PointXY
    x::Int
    y::Float64
end

mutable struct SegmentBuffer
    Δx::Vector{Int}
    Δy::Vector{Float64}
    slope::Vector{Float64}
    first::Int
    last::Int
    SegmentBuffer(n::Int) = new(zeros(Int, n), zeros(n), zeros(n), 1, 0)
end

# pretty printer of segment buffers for debugging purposes only
function Base.show(io::IO, sb::SegmentBuffer)
    nsegments = sb.last - sb.first + 1
    print(io, "buffer with ", sb.last - sb.first + 1, " segment(s)")
    (nsegments >= 1) && print(io, "\n",
        join(collect(zip(sb.Δx[sb.first:sb.last], sb.Δy[sb.first:sb.last])), "-"))
end

@inline getindex(sb::SegmentBuffer, i::Int) = (sb.Δx[i], sb.Δy[i], sb.slope[i])

@inline setindex!(sb::SegmentBuffer, val::Tuple{Int, Float64, Float64}, i::Int) =
    ((sb.Δx[i], sb.Δy[i], sb.slope[i]) = val; nothing)

@inline num_segments(sb::SegmentBuffer) = sb.last - sb.first + 1

@inline remove_first_segment!(sb::SegmentBuffer) = (sb.first += 1; nothing)

@inline pop_first_segment!(sb::SegmentBuffer) = (sb.first += 1; sb[sb.first - 1])

@inline first_segment(sb::SegmentBuffer) = sb[sb.first]

@inline last_segment(sb::SegmentBuffer) = sb[sb.last]

@inline first_slope(sb::SegmentBuffer) = sb.slope[sb.first]

@inline last_slope(sb::SegmentBuffer) = sb.slope[sb.last]

@inline pop_last_segment!(sb::SegmentBuffer) = (sb.last -= 1; sb[sb.last + 1])

@inline push_to_last!(sb::SegmentBuffer, Δx::Int, Δy::Float64) =
    (sb.last += 1; sb[sb.last] = (Δx, Δy, Δy / float(Δx)); nothing)

@inline restart_from_state!(sb::SegmentBuffer) =
    (sb.first = sb.last + 1; nothing)
    # (sb.first = sb.last; nothing)
    # (sb.first = sb.last; sb.last -= 1; nothing)

@inline function add_segment_to_minorant!(
    minorant::SegmentBuffer,
    Δx::Int,
    Δy::Float64
)
    δx, δy = Δx, Δy
    if δy < float(δx) * last_slope(minorant) # if not convex
        auxi = num_segments(minorant)
        while true # do-while block
            # Merge last two segments
            dx, dy, s = pop_last_segment!(minorant)
            δx += dx
            δy += dy

            # Check if convex again
            auxi -= 1
            (auxi < 1 || δy >= float(δx) * last_slope(minorant)) && break
        end
    end
    push_to_last!(minorant, δx, δy)
end

@inline function add_segment_to_majorant!(
    majorant::SegmentBuffer,
    Δx::Int,
    Δy::Float64
)
    δx, δy = Δx, Δy
    if δy > float(δx) * last_slope(majorant) # if not concave
        auxi = num_segments(majorant)
        while true # do-while block
            # Merge last two segments
            dx, dy, s = pop_last_segment!(majorant)
            δx += dx
            δy += dy

            # Check if concave again
            auxi -= 1
            (auxi < 1 || δy <= float(δx) * last_slope(majorant)) && break
        end
    end
    push_to_last!(majorant, δx, δy)
end

@inline function new_knot!(
    minorant::SegmentBuffer,
    majorant::SegmentBuffer,
    λ::Float64,
    origin::PointXY,
    last_explored::PointXY
)
    # First segments
    Δxmin, Δymin, slopemin = first_segment(minorant)
    Δxmaj, Δymaj, slopemaj = first_segment(majorant)

    # Shortest segment defines the new knot
    if Δxmin < Δxmaj
        # Remove first segment, rest of the minorant is still valid
        remove_first_segment!(minorant)

        # Majorant is a single segment from new knot to last original point
        restart_from_state!(majorant)
        Δx = last_explored.x - origin.x - Δxmin
        Δy = last_explored.y - origin.y - Δymin - λ
        push_to_last!(majorant, Δx, Δy)

        return Δxmin, Δymin, slopemin
    else # Left-most majorant touching point is the new knot
        remove_first_segment!(majorant)

        # Minorant is a single segment from new knot to last original point
        restart_from_state!(minorant)
        Δx = last_explored.x - origin.x - Δxmaj
        Δy = last_explored.y - origin.y - Δymaj + λ
        push_to_last!(minorant, Δx, Δy)

        return Δxmaj, Δymaj, slopemaj
    end
end

function filter1D!( # Weighted version
    β::AbstractVector{Float64},
    y::AbstractVector{Float64},
    λ::AbstractVector{Float64},
    from::Int,
    to::Int
)
    # problem size
    n = to - from + 1
    # @assert length(y) == length(λ) + 1
    # @assert all([λᵢ >= 0 for λᵢ ∈ λ])

    # @fastmath @inbounds begin # acceleration

    # Initialise memory structures
    majorant = SegmentBuffer(n)
    minorant = SegmentBuffer(n)
    push_to_last!(majorant, 1, y[from] - λ[from])
    push_to_last!(minorant, 1, y[from] + λ[from])

    # Initial point of taut-string and last explored
    origin = PointXY(from - 1, 0.)
    last_explored = PointXY(from, y[from])

    # Iterate along the signal length
    for i = (from + 1):(from + n - 2)
        # Update majorant, minorant and last explored point
        # Take into account difference in lambdas
        Δλ = λ[i] - λ[i - 1]
        add_segment_to_majorant!(majorant, 1, y[i] - Δλ)
        add_segment_to_minorant!(minorant, 1, y[i] + Δλ)
        last_explored.x += 1
        last_explored.y += y[i]

        # Check for slope crossings at the first point
        while first_slope(minorant) < first_slope(majorant)
            # Crossing detected
            Δx, Δy, slope = new_knot!(minorant, majorant, λ[i], origin, last_explored)

            # Write to solution
            @simd for j in 1:Δx
                β[origin.x + j] = slope
            end

            # Update origin
            origin.x += Δx
            origin.y += Δy
        end
    end

    # Update majorant and minorant with last segments
    add_segment_to_majorant!(majorant, 1, y[from + n - 1] + λ[from + n - 2]) # Δλ = 0. - λ[n - 1]
    add_segment_to_minorant!(minorant, 1, y[from + n - 1] - λ[from + n - 2])

    # At this point, because the endpoint of the tube is the same
    # for both majorant and minorant, either the majorant or the minorant
    # is a single straight segment while the other can be larger.
    # The remaining of the taut-string must be the multi-segment component.
    largest_buffer = (num_segments(majorant) > num_segments(minorant) ? majorant : minorant)
    for i in 1:num_segments(largest_buffer)
        # Write to solution
        Δx, Δy, slope = pop_first_segment!(largest_buffer)
        @simd for j in 1:Δx
            β[origin.x + j] = slope
        end
        origin.x += Δx
    end

    # end # end fastmath and inbounds
end

# Unweighted version
function filter1D!(
    β::AbstractVector{Float64},
    y::AbstractVector{Float64},
    λ::Float64,
    from::Int,
    to::Int
)
    filter1D!(β, y, fill(λ, length(y)), from, to)
end

function filter1D!( #
    β::AbstractVector{Float64},
    y::AbstractVector{Float64},
    λ::Union{AbstractVector{Float64}, Float64}
)
    filter1D!(β, y, λ, 1, length(y))
end

# non-inplace version
function filter1D(
    y::AbstractVector{Float64},
    λ::Union{AbstractVector{Float64}, Float64}
)
    # trivial cases
    λ == 0. && (return y)
    β = zeros(length(y))
    filter1D!(β, y, λ)
    β
end


text:
using BeamformerRecipes
using RadioInterferometry
using YAML

"""
    TelInfo(yaml::Dict{Symbol,Any})

Create a `TelInfo` instance based on `yaml` which was loaded from a "telinfo.yml"
file.
"""
function BeamformerRecipes.TelInfo(yaml::Dict{Symbol,Any})
    ti = BeamformerRecipes.TelInfo()
    # Populate the easy fields first
    ti.antenna_position_frame = string(get(yaml, :antenna_position_frame, nothing))
    ti.latitude = dms2d(get(yaml, :latitude, NaN))
    ti.longitude = dms2d(get(yaml, :longitude, NaN))
    ti.altitude = Float64(get(yaml, :altitude, NaN))
    ti.telescope_name = string(get(yaml, :telescope_name, "Unknown"))

    # Now do antenna specific fields
    ants = yaml[:antennas]
    diam = get(yaml, :antenna_diameter, 0.0)
    ti.antenna_names = map(a->string(a[:name]), ants)
    ti.antenna_numbers = map(a->Int(a[:number]), ants)
    ti.antenna_positions = mapreduce(a->convert(Vector{Float64}, a[:position]), hcat, ants)
    ti.antenna_diameters = map(a->Float64(get(a, :diameter, diam)), ants)

    # Validate that antenna positions are 3 dimensional
    @assert size(ti.antenna_positions, 1) == 3 "antenna positions are not 3D"

    ti
end

"""
    TelInfo(yaml::AbstractString)

Create a `TelInfo` instance based on a "telinfo.yml" file named by `yaml`.
"""
function BeamformerRecipes.TelInfo(yaml::AbstractString)
    BeamformerRecipes.TelInfo(YAML.load_file(yaml, dicttype=Dict{Symbol,Any}))
end

"""
    antpos_topo_xyz(telinfo::TelInfo)

Return `telinfo.antenna_positions` from meters in given frame to nanoseconds in
topocentric XYZ frame.
"""
function antpos_topo_xyz(telinfo::TelInfo)
    # If telinfo frame is nothing, assume enu
    inframe = something(telinfo.antenna_position_frame, "enu")
    if inframe == "enu"
        antpos = enu2xyz(telinfo.antenna_positions, deg2rad(telinfo.latitude), 0)
    else
        error("unsupported antenna position frame: $inframe")
    end
end

ext: jl


text:
"""
    composecolors(
        images,
        cmap=["#F00", "#0F0", "#00F"];
        clims,
        stretch,
        contrast,
        bias,
        multiplier
    )

Create a color composite of multiple images by applying `imview` and blending
the results. This function can be used to create RGB composites using any number of channels
(e.g. red, green, blue, and hydrogen alpha) as well as more exotic images like blending
radio and optical data using two different colormaps.

`cmap` should be a list of colorants, named colors (see Colors.jl), or colorschemes (see ColorSchemes.jl).
`clims`, `stretch`, `contrast`, and `bias` are passed on to `imview`. They can be a single value or
a list of different values for each image.

Examples:
```julia
# Basic RGB
composecolors([redimage, greenimage, blueimage])
# Non-linear stretch before blending
composecolors([redimage, greenimage, blueimage], stretch=asinhstretch)
# More than three channels are allowed (H alpha in pink)
composecolors(
    [antred, antgreen, antblue, anthalp],
    ["red", "green", "blue", "maroon1"],
    multiplier=[1,2,1,1]
)
# Can mix
composecolors([radioimage, xrayimage], [:ice, :magma], clims=extrema)
composecolors([radioimage, xrayimage], [:magma, :viridis], clims=[Percent(99), Zscale()])
```
"""
function composecolors(
    images,
    cmap=nothing;
    clims=_default_clims[],
    stretch=_default_stretch[],
    contrast=1.0,
    bias=0.5,
    multiplier=1.0
)
    if isempty(images)
        error("At least one image is required.")
    end
    if !allequal(size.(images))
        error("Images must have the same dimensions to compose them.")
    end
    if length(images) == 3 && isnothing(cmap)
        cmap = ["red", "green", "blue"]
    end
    if length(cmap) < length(images)
        error("Please provide a color channel for each image")
    end

    # Use imview to render each channel to RGBA
    images_rendered = broadcast(images, cmap, clims, stretch, contrast, bias) do image, cmap, clims, stretch, contrast, bias
        imview(image; cmap, clims, stretch, contrast, bias)
    end

    # Now blend, ensuring each color channel never exceeds [0,1]
    combined = mappedarray(images_rendered...) do channels...
        pxblended = sum(channels .* multiplier)
        return typeof(pxblended)(
            clamp(pxblended.r,0,1),
            clamp(pxblended.g,0,1),
            clamp(pxblended.b,0,1),
            clamp(pxblended.alpha,0,1)
        )
    end
    return combined
end
export composechannels


text:
@testset "Simulation Sequence Correct Execution Order" begin
    models_array = [
        DecisionModel(
            MockOperationProblem;
            horizon=48,
            interval=Hour(24),
            steps=2,
            name="DAUC",
        ),
        DecisionModel(
            MockOperationProblem;
            horizon=24,
            interval=Hour(1),
            steps=2 * 24,
            name="HAUC",
        ),
        DecisionModel(
            MockOperationProblem;
            horizon=12,
            interval=Minute(5),
            steps=2 * 24 * 12,
            name="ED",
        ),
    ]
    set_device_model!(
        PSI.get_template(models_array[3]),
        ThermalStandard,
        ThermalBasicDispatch,
    )
    models = SimulationModels(
        models_array,
        EmulationModel(MockEmulationProblem; resolution=Minute(1), name="AGC"),
    )

    test_sequence = SimulationSequence(
        models=models,
        feedforwards=Dict(
            "ED" => [
                SemiContinuousFeedforward(
                    component_type=ThermalStandard,
                    source=OnVariable,
                    affected_values=[ActivePowerVariable],
                ),
            ],
        ),
        ini_cond_chronology=InterProblemChronology(),
    )

    @test !isempty(
        PSI.get_model(PSI.get_template(models_array[3]), ThermalStandard).feedforwards,
    )

    @test length(findall(x -> x == 4, test_sequence.execution_order)) == 24 * 60
    @test length(findall(x -> x == 3, test_sequence.execution_order)) == 24 * 12
    @test length(findall(x -> x == 2, test_sequence.execution_order)) == 24
    @test length(findall(x -> x == 1, test_sequence.execution_order)) == 1

    for model in PSI.get_decision_models(models)
        @test model.internal.simulation_info.sequence_uuid == test_sequence.uuid
    end

    # Test single stage sequence
    test_sequence = SimulationSequence(
        models=SimulationModels(
        # TODO: support passing one model without making a vector
            [DecisionModel(MockOperationProblem; horizon=48, name="DAUC")]),
        ini_cond_chronology=InterProblemChronology(),
    )

    # Disabled temporarily
    # @test isa(test_sequence.ini_cond_chronology, IntraProblemChronology)
    @test test_sequence.execution_order == [1]
end

@testset "Simulation Sequence invalid sequences" begin
    models = SimulationModels([
        DecisionModel(
            MockOperationProblem;
            horizon=48,
            interval=Hour(24),
            steps=2,
            name="DAUC",
        ),
        DecisionModel(
            MockOperationProblem;
            horizon=24,
            interval=Hour(5),
            steps=2 * 24,
            name="HAUC",
        ),
    ])

    @test_throws IS.ConflictingInputsError SimulationSequence(models=models)

    models = SimulationModels([
        DecisionModel(
            MockOperationProblem;
            horizon=2,
            interval=Hour(1),
            steps=2,
            name="DAUC",
        ),
        DecisionModel(
            MockOperationProblem;
            horizon=24,
            interval=Hour(1),
            steps=2 * 24,
            name="HAUC",
        ),
    ])

    @test_throws IS.ConflictingInputsError SimulationSequence(models=models)

    models = SimulationModels([
        DecisionModel(
            MockOperationProblem;
            horizon=24,
            interval=Hour(1),
            steps=2,
            name="DAUC",
        ),
        DecisionModel(
            MockOperationProblem;
            horizon=24,
            interval=Minute(22),
            steps=2 * 24,
            name="HAUC",
        ),
    ])

    @test_throws IS.ConflictingInputsError SimulationSequence(models=models)
end


text:
#  ArgParse
using ArgParse
using spiavalonfpga
using ProgressMeter

function parse_commandline()
s = ArgParseSettings()
@add_arg_table s begin
    #=
    "--opt1"
        help = "an option with an argument"
    "--opt2", "-o"
        help = "another option with an argument"
        arg_type = Int
        default = 0
    "--flag1"
        help = "an option without argument, i.e. a flag"
        action = :store_true
        =#
    "arg1"
        help = "rotator position"
        arg_type = Int
        required = true
end
return parse_args(s)
end

function rotator()
    parsed_args = parse_commandline()
    println("Parsed args:")
    for (arg,val) in parsed_args
        println("  $arg  =>  $val")
    end
    println(parsed_args["arg1"])
    buffer=zeros(UInt8, 1)
    buffer[1] = parsed_args["arg1"]
    spiavalonfpga.SPIManager.writeLEDs(0x0400_0000, buffer)
#    @showprogress for i in 1:5
#    read_back = spiavalonfpga.SPIManager.readbyte(0x0400_0000)
#    println("Read back of motor position = $read_back")
#    end
    while ((position = spiavalonfpga.SPIManager.readbyte(0x0400_0000)) != buffer[1])
        println("Rotation in Progress")
    end
    position = spiavalonfpga.SPIManager.readbyte(0x0400_0000)
    println("Confirming position is $position.")
end

rotator()

max_stars_repo_name: Greg-R/spiavalonfpga

max_issues_repo_name: Greg-R/spiavalonfpga

max_forks_repo_name: Greg-R/spiavalonfpga


text:
using Agents, Random
using GLMakie
using InteractiveDynamics

model, agent_step!, model_step! = Models.social_distancing(isolated = 0.8)

# Various possibilities for colors/sizes/markers
sir_colors(a) = a.status == :S ? "#2b2b33" : a.status == :I ? "#bf2642" : "#338c54"
# sir_sizes(a) = 10*randn()
# sir_sizes(a) = 5*(mod1(a.id, 3)+1)
sir_sizes = 10

# sir_shape(a) = rand(('🐑', '🐺', '🌳'))
# sir_shape(a) = rand(('😹', '🐺', '🌳'))
# sir_shape(a) = rand(('π', '😹', '⚃', '◑', '▼'))
# sir_shape(a) = rand((:diamond, :circle))
sir_shape(a) = a.status == :S ? :circle : a.status == :I ? :diamond : :rect

# function sir_shape(b)
#     φ = atan(b.vel[2], b.vel[1])
#     xs = [(i ∈ (0, 3) ? 2 : 1)*cos(i*2π/3 + φ) for i in 0:3]
#     ys = [(i ∈ (0, 3) ? 2 : 1)*sin(i*2π/3 + φ) for i in 0:3]
#     poly(xs, ys)
# end

fig, abmstepper = abm_plot(model;
    ac = sir_colors, as = sir_sizes, am = sir_shape,
)

display(fig)

# %% same, but interactive
abm_play(model, agent_step!, model_step!;
    ac = sir_colors, as = sir_sizes, am = sir_shape,
)

# %% same, but video
abm_video("socialdist.mp4", model, agent_step!, model_step!;
    ac = sir_colors, as = sir_sizes, am = sir_shape,
    spf = 2, frames = 100,
)

# %% Interactive data plot
model, agent_step!, model_step! = Models.social_distancing(isolated = 0.8)

infected(x) = count(i == :I for i in x)
recovered(x) = count(i == :R for i in x)
adata = [(:status, infected), (:status, recovered)]
alabels = ["I", "R"]
mdata = [nagents]
mlabels = ["N"]

params = Dict(
    :death_rate => 0.02:0.001:1.0,
    :reinfection_probability => 0:0.01:1.0,
    :infection_period => 24:24:2400,
)

when = (model, s) -> s % 50 == 0

fig, adf, mdf = abm_data_exploration(model, agent_step!, model_step!, params;
ac = sir_colors, as = sir_sizes, am = sir_shape,
when = when, mdata = mdata, adata=adata, alabels=alabels, mlabels=mlabels)


text:
mapping_args(ts) = :(@nospecialize(_::Union{$((:(Type{$t}) for t in ts)...)}))

function HLTypeMapping(spec)
    ts = [spec.name]
    has_intermediate_type(spec) && push!(ts, struct_name(spec))
    ex = :(hl_type($(mapping_args(ts))) = $(struct_name(spec, true)))
    HLTypeMapping(spec, ex)
end

function CoreTypeMapping(spec)
    ts = [spec.name, struct_name(spec, true)]
    has_intermediate_type(spec) && push!(ts, struct_name(spec))
    ex = :(core_type($(mapping_args(ts))) = $(spec.name))
    CoreTypeMapping(spec, ex)
end

function IntermediateTypeMapping(spec)
    ts = [struct_name(spec, true), spec.name]
    ex = :(intermediate_type($(mapping_args(ts))) = $(struct_name(spec)))
    IntermediateTypeMapping(spec, ex)
end

max_stars_repo_name: JuliaGPU/VulkanAbstraction

max_issues_repo_name: JuliaGPU/VulkanAbstraction

max_forks_repo_name: JuliaGPU/VulkanAbstraction


text:
struct BoolBridge{T, B<:MOIBC.SetMapBridge{T}, F1, F2, S1, S2} <: MOIBC.SetMapBridge{T, S2, S1, F2, F1}
    con_idx::CI
end

function MOI.supports_constraint(
    bridge::Type{<:BoolBridge{T, B}},
    fct::Type{<:MOI.AbstractFunction},
    set::Type{<:CS.AbstractBoolSet{F1,F2,F1dim,F2dim,S1,S2}}
) where {T, B, F1, F2, F1dim, F2dim, S1, S2}
    # check whether there is at least one constraint which is supported by the inner bridge
    supported = uses_bool_inner(B, fct, set)
    return supported
end

function uses_bool_inner(
    inner_bridge_type,
    fct::Type{<:MOI.AbstractFunction},
    set::Type{<:CS.AbstractBoolSet{F1,F2,F1dim,F2dim,S1,S2}}
) where {F1, F2, F1dim, F2dim, S1, S2}
    return uses_bool_inner(inner_bridge_type, F1, S1) || uses_bool_inner(inner_bridge_type, F2, S2)
end

function uses_bool_inner(
    inner_bridge_type,
    F::Type{<:MOI.AbstractFunction},
    S
)
    return MOI.supports_constraint(inner_bridge_type, F, S)
end

function added_constraint_types(
    bridge::Type{<:BoolBridge{T, B}},
    ::Type{<:CS.AbstractBoolSet{F1,F2,F1dim,F2dim,S1,S2}}
) where {T,B, F1, F2, F1dim, F2dim, S1, S2}
    get_bridged_and_set(bridge,F1,F2,F1dim,F2dim,S1,S2)
end

function MOIBC.concrete_bridge_type(
    ::Type{<:BoolBridge{T,B}},
    G::Type{<:MOI.AbstractFunction},
    ::Type{<:CS.AbstractBoolSet{F1,F2,F1dim,F2dim,S1,S2}}
) where {T,B,F1, F2, F1dim, F2dim, S1, S2}
    return BoolBridge{T,B,F1,F2,S1,S2}
end

function MOIBC.bridge_constraint(bridge::Type{<:BoolBridge{T, B, F1, F2, S1, S2}}, model, fct, set::CS.AbstractBoolSet) where {T, B, F1, F2, S1, S2}
    new_fct = map_function(bridge, fct, set)
    new_set = MOIBC.map_set(bridge, set)
    return BoolBridge{T,B,F1,F2,S1,S2}(MOI.add_constraint(model, new_fct, new_set))
end

function unpack_constraint_types(ct1, ct2)
    new_F1 = ct1[1][1]
    new_F2 = ct2[1][1]
    new_S1 = ct1[1][2]
    new_S2 = ct2[1][2]
    return new_F1, new_F2, new_S1, new_S2
end

function how_supported(
    B::Type{<:MOIBC.SetMapBridge{T}},
    F, S
) where {T}
    supported_by_B = MOI.supports_constraint(B, F, S)
    if !supported_by_B
        return true, nothing
    elseif supported_by_B
        return false, B
    end
end

function get_constraint_types(
    B::Type{<:MOIBC.SetMapBridge{T}},
    F, S
) where {T}
    direct_support, inner_bridge = how_supported(B, F, S)
    if direct_support
        return [(F,S)]
    else
        return MOIB.added_constraint_types(inner_bridge, F, S)
    end
end

function get_mapped_fct(
    B::Type{<:MOIBC.SetMapBridge{T}},
    F, S,
    fct,
) where {T}
    direct_support, inner_bridge = how_supported(B, F, S)
    if direct_support
        return fct
    else
        return MOIBC.map_function(inner_bridge, fct)
    end
end

function get_mapped_set(
    B::Type{<:MOIBC.SetMapBridge{T}},
    F, S,
    set,
) where {T}
    direct_support, inner_bridge = how_supported(B, F, S)
    if direct_support
        return set
    else
        return MOIBC.map_set(inner_bridge, set)
    end
end

function get_bridged_and_set(
    bridge::Type{<:BoolBridge{T, B}},
    F1::Type{<:MOI.AbstractFunction},
    F2::Type{<:MOI.AbstractFunction},
    F1dim,
    F2dim,
    S1::Type{<:AbstractBoolSet{F11,F12,F11dim,F12dim,S11,S12}},
    S2::Type{<:AbstractBoolSet{F21,F22,F21dim,F22dim,S21,S22}}
) where {T, B, F11, F12,F11dim,F12dim, S11, S12, F21, F22, F21dim,F22dim, S21, S22}
    lhs_constraint_types = get_bridged_and_set(bridge, F11, F12, F11dim, F12dim, S11, S12)
    rhs_constraint_types = get_bridged_and_set(bridge, F21, F22, F21dim, F22dim, S21, S22)
    new_F1, new_F2, new_S1, new_S2 = unpack_constraint_types(lhs_constraint_types, rhs_constraint_types)
    return [(MOI.VectorAffineFunction{T},AbstractBoolSet{new_F1,new_F2,F1dim,F2dim,new_S1,new_S2})]
end

function get_bridged_and_set(
    bridge::Type{<:BoolBridge{T, B}},
    F1::Type{<:MOI.AbstractFunction},
    F2::Type{<:MOI.AbstractFunction},
    F1dim,
    F2dim,
    S1::Type{<:AbstractBoolSet{F11,F12,F11dim,F12dim,S11,S12}},
    S2
) where {T, B, F11, F12, F11dim, F12dim, S11, S12}
    lhs_constraint_types = get_bridged_and_set(bridge, F11, F12, F11dim, F12dim, S11, S12)
    rhs_constraint_types = get_constraint_types(B, F2, S2)
    new_F1, new_F2, new_S1, new_S2 = unpack_constraint_types(lhs_constraint_types, rhs_constraint_types)
    return [(MOI.VectorAffineFunction{T},AbstractBoolSet{new_F1,new_F2,F1dim,F2dim,new_S1,new_S2})]
end

function get_bridged_and_set(
    bridge::Type{<:BoolBridge{T, B}},
    F1::Type{<:MOI.AbstractFunction},
    F2::Type{<:MOI.AbstractFunction},
    F1dim,
    F2dim,
    S1,
    S2::Type{<:AbstractBoolSet{F21,F22,F21dim,F22dim,S21,S22}}
) where {T, B, F21, F22, F21dim, F22dim, S21, S22}
    lhs_constraint_types = get_constraint_types(B, F1, S1)
    rhs_constraint_types = get_bridged_and_set(bridge, F21, F22, F21dim, F22dim, S21, S22)
    new_F1, new_F2, new_S1, new_S2 = unpack_constraint_types(lhs_constraint_types, rhs_constraint_types)
    return [(MOI.VectorAffineFunction{T},AbstractBoolSet{new_F1,new_F2,F1dim,F2dim,new_S1,new_S2})]
end

function get_bridged_and_set(
    bridge::Type{<:BoolBridge{T, B}},
    F1::Type{<:MOI.AbstractFunction},
    F2::Type{<:MOI.AbstractFunction},
    F1dim,
    F2dim,
    S1,
    S2,
) where {T, B}
    lhs_constraint_types = get_constraint_types(B, F1, S1)
    rhs_constraint_types = get_constraint_types(B, F2, S2)
    new_F1, new_F2, new_S1, new_S2 = unpack_constraint_types(lhs_constraint_types, rhs_constraint_types)
    return [(MOI.VectorAffineFunction{T}, AbstractBoolSet{new_F1,new_F2,F1dim,F2dim,new_S1,new_S2})]
end

function map_function(
    bridge::Type{<:BoolBridge{T, B}},
    fct,
    set::CS.AbstractBoolSet{F1,F2,F1dim,F2dim,S1,S2}
) where {T,B,F1,F2,F1dim,F2dim,S1<:AbstractBoolSet,S2<:AbstractBoolSet}
    f = MOIU.eachscalar(fct)
    lhs_fct = map_function(bridge, f[1:get_value(F1dim)], set.lhs_set)
    rhs_fct = map_function(bridge, f[end-get_value(F2dim)+1:end], set.rhs_set)
    return MOIU.operate(vcat, T, lhs_fct, rhs_fct)
end

function map_function(
    bridge::Type{<:BoolBridge{T, B}},
    fct,
    set::CS.AbstractBoolSet{F1,F2,F1dim,F2dim,S1,S2}
) where {T,B,F1,F2,F1dim,F2dim,S1<:AbstractBoolSet,S2}
    f = MOIU.eachscalar(fct)
    lhs_fct = map_function(bridge, f[1:get_value(F1dim)], set.lhs_set)
    rhs_fct = get_mapped_fct(B, F2, S2, f[end-get_value(F2dim)+1:end])
    return MOIU.operate(vcat, T, lhs_fct, rhs_fct)
end

function map_function(
    bridge::Type{<:BoolBridge{T, B}},
    fct,
    set::CS.AbstractBoolSet{F1,F2,F1dim,F2dim,S1,S2}
) where {T,B,F1,F2,F1dim,F2dim,S1,S2<:AbstractBoolSet}
    f = MOIU.eachscalar(fct)
    lhs_fct = get_mapped_fct(B, F1, S1, f[1:get_value(F1dim)])
    rhs_fct = map_function(bridge, f[end-get_value(F2dim)+1:end], set.rhs_set)
    return MOIU.operate(vcat, T, lhs_fct, rhs_fct)
end

function map_function(
    bridge::Type{<:BoolBridge{T, B}},
    fct,
    set::CS.AbstractBoolSet{F1,F2,F1dim,F2dim,S1,S2}
) where {T,B,F1,F2,F1dim,F2dim,S1,S2}
    f = MOIU.eachscalar(fct)
    lhs_fct = get_mapped_fct(B, F1, S1, f[1:get_value(F1dim)])
    rhs_fct = get_mapped_fct(B, F2, S2, f[end-get_value(F2dim)+1:end])
    return MOIU.operate(vcat, T, lhs_fct, rhs_fct)
end

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

function MOIBC.map_set(
    bridge::Type{<:BoolBridge{T, B}},
    set::CS.AbstractBoolSet{F1,F2,F1dim,F2dim,S1,S2}
) where {T,B,F1,F2,F1dim,F2dim,S1<:AbstractBoolSet,S2<:AbstractBoolSet}
    lhs_set = MOIBC.map_set(bridge, set.lhs_set)
    rhs_set = MOIBC.map_set(bridge, set.rhs_set)
    return typeof_without_params(set){F1,F2}(lhs_set, rhs_set)
end

function MOIBC.map_set(
    bridge::Type{<:BoolBridge{T, B}},
    set::CS.AbstractBoolSet{F1,F2,F1dim,F2dim,S1,S2}
) where {T,B,F1,F2,F1dim,F2dim,S1<:AbstractBoolSet,S2}
    lhs_set = MOIBC.map_set(bridge, set.lhs_set)
    rhs_set = get_mapped_set(B, F2, S2, set.rhs_set)
    return typeof_without_params(set){F1,F2}(lhs_set, rhs_set)
end

function MOIBC.map_set(
    bridge::Type{<:BoolBridge{T, B}},
    set::CS.AbstractBoolSet{F1,F2,F1dim,F2dim,S1,S2}
) where {T,B,F1,F2,F1dim,F2dim,S1,S2<:AbstractBoolSet}
    lhs_set = get_mapped_set(B, F1, S1, set.lhs_set)
    rhs_set = MOIBC.map_set(bridge, set.rhs_set)
    return typeof_without_params(set){F1,F2}(lhs_set, rhs_set)
end

function MOIBC.map_set(
    bridge::Type{<:BoolBridge{T, B}},
    set::CS.AbstractBoolSet{F1,F2,F1dim,F2dim,S1,S2}
) where {T,B,F1,F2,F1dim,F2dim,S1,S2}
    lhs_set = get_mapped_set(B, F1, S1, set.lhs_set)
    rhs_set = get_mapped_set(B, F2, S2, set.rhs_set)

    return typeof_without_params(set){F1,F2}(lhs_set, rhs_set)
end


text:
function addto!(jm::ModelView, m::EnergyModel, bus::Bus)
    injection = expression(NodalActivePower(), m, bus)
    B, T = AxisArrays.axes(injection)
    @constraint(jm, p_balance[b=B,t=T], injection[b,t] == 0)
end

addcomponent(Bus, :buses, (:B, :T=>:snapshots), joinpath(@__DIR__, "attrs", "buses.csv"))

max_stars_repo_name: coroa/EnergyModels

max_issues_repo_name: coroa/EnergyModels

max_forks_repo_name: coroa/EnergyModels


text:
using ReconTPCF
using Test

#ReconTPCF.loadim("C:/Users/James/.julia/dev/ReconTPCF/test/uniform_circles.png")
path = "C:/Users/James/.julia/dev/ReconTPCF"

begin
    test = loadim("./test/uniform_circles.png")
    dims = size(test)#(200,100)
    maxrng = Int64(round(sqrt(dims[1]^2 + dims[2]^2)+1))
    philen = sum(test)#Int(round(0.3*dims[1]*dims[2]))

    SN = SN_comp(dims, maxrng) #Is SN_comp supposed to be multiplied by 2?

    idx = CartesianIndices(dims)
    wpix = findall(x->x==true, test)
    bpix = setdiff(idx, wpix)
    pix = (wpix, bpix)
    S2N = S2_initialise(pix, maxrng)
    S2 = naninf(S2_finalise(SN, S2N, test))
    #S2[length(S2)÷2:end] .= 0
    C2N = C2_initialise(test, maxrng)
    C2 = naninf(S2_finalise(SN, C2N, test))
end

@testset "ReconTPCF.jl" begin
    # Write your tests here.
#ReconTPCF.loadim("uniform_circles.png")
end


text:
module SexSpecificMigrations

using ForwardPopGenSimulations

# reexport from ForwardPopGenSimulations as users use these to specify model parameters.
export Female,
       Male,
       Mother,
       Father,
       Daughter,
       Son,

       Autosome,
       XChromosome,
       YChromosome,
       Mitochondrion,

       BasicData

export ModelParameters,
       LocusList,
       PopulationSizes,
       MutationRates,
       RecombinationRates,
       ProportionsOfMigrants,
       FitnessOfMigrants,
       FitnessOfForeignTraits,
       TraitTransmissionRates,

       simulate

immutable ForeignTrait end
immutable Migrant end

immutable LocusList
    val::Array{Tuple{ChromosomeType, Int}, 1}
end
Base.call(::Type{LocusList}, val::Tuple{ChromosomeType, Int}...) = LocusList([val...])
Base.length(t::LocusList) = length(t.val)
Base.start(t::LocusList) = start(t.val)
Base.next(t::LocusList, s) = next(t.val, s)
Base.done(t::LocusList, s) = done(t.val, s)

abstract Parameter
abstract DemeAndSexSpecificParameter <: Parameter
abstract TwoGensSexSpecificParameter <: Parameter
abstract SexSpecificParameter <: Parameter
Base.getindex(p::Parameter, i::Integer) = p.data[i]

Base.length(i::DemeAndSexSpecificParameter) = length(i.data)

function Base.call{T<:DemeAndSexSpecificParameter, V<:Number}(
    ::Type{T}, vals::Pair{Int, NTuple{2, Pair{SexType, V}}}...)

    d = Dict([val[1] => Dict(val[2]) for val in vals])
    demes = sort(collect(keys(d)))
    T(tuple([(d[deme][Female], d[deme][Male]) for deme in demes]...))
end

function Base.call{T<:TwoGensSexSpecificParameter, V<:Number}(
    ::Type{T}, vals::Pair{SexType, NTuple{2, Pair{SexType, V}}}...)

    d = Dict([val[1] => Dict(val[2]) for val in vals])
    T(((d[Mother][Daughter], d[Mother][Son]), (d[Father][Daughter], d[Father][Son])))
end

function Base.call{T<:SexSpecificParameter, V<:Number}(
    ::Type{T}, vals::Pair{SexType, V}...)
    d = Dict(vals)
    T((d[Female], d[Male]))
end

immutable PopulationSizes <: DemeAndSexSpecificParameter
    # deme and sex-specific
    data::NTuple{2, NTuple{2, Int}}
    PopulationSizes(d::NTuple{2, NTuple{2, Int}}) = new(d)
end

immutable MutationRates <: SexSpecificParameter
    # sex-specific
    data::NTuple{2, Float64}
    MutationRates(d::NTuple{2, Float64}) = new(d)
end

immutable RecombinationRates <: SexSpecificParameter
    # sex-specific
    data::NTuple{2, Float64}
    RecombinationRates(d::NTuple{2, Float64}) = new(d)
end

immutable ProportionsOfMigrants <: DemeAndSexSpecificParameter
    # deme and sex-specific
    data::NTuple{2, NTuple{2, Float64}}
    ProportionsOfMigrants(d::NTuple{2, NTuple{2, Float64}}) = new(d)
end

abstract Fitness <: DemeAndSexSpecificParameter
immutable FitnessOfMigrants <: Fitness
    # deme and sex-specific
    data::NTuple{2, NTuple{2, Float64}}
    FitnessOfMigrants(d::NTuple{2, NTuple{2, Float64}}) = new(d)
end

immutable FitnessOfForeignTraits <: Fitness
    # deme and sex-specific
    data::NTuple{2, NTuple{2, Float64}}
    FitnessOfForeignTraits(d::NTuple{2, NTuple{2, Float64}}) = new(d)
end

immutable TraitTransmissionRates <: TwoGensSexSpecificParameter
    # sex-specific in parents and offspring
    data::NTuple{2, NTuple{2, Float64}}
    TraitTransmissionRates(d::NTuple{2, NTuple{2, Float64}}) = new(d)
end

function trait(momrate, poprate, momtrait, poptrait, deme)
    if deme == momtrait == poptrait
        deme
    elseif deme != momtrait && deme != poptrait
        rand() < 1.0 - (1.0 - momrate) * (1.0 - poprate) ? momtrait : deme
    elseif deme != momtrait
        rand() < momrate ? momtrait : deme
    else
        rand() < poprate ? poptrait : deme
    end
end

immutable ModelParameters
    chromosomes::LocusList
    population_sizes::PopulationSizes
    mutation_rates::MutationRates
    recombination_rates::RecombinationRates
    proportions_of_migrant::ProportionsOfMigrants
    fitness_of_migrants::FitnessOfMigrants
    fitness_of_foreign_traits::FitnessOfForeignTraits
    transmission_rates_of_trait::TraitTransmissionRates
end

immutable SexSpecificParameters
    deme::Int
    population_size::Int
    number_of_loci::Int
    sex::SexType
    mutation_rate::Float64
    recombination_rate::Array{Float64, 1}
    type_of_locus::Array{ChromosomeType, 1}
    proportion_of_migrant::Float64
    fitness_of_migrant::Float64
    fitness_of_foreign_trait::Float64
    transmission_rate_of_trait::NTuple{2, Float64}
end

function SexSpecificParameters(p::ModelParameters, deme::Integer, sex::SexType)
    sexidx = value(sex)

    chrs = p.chromosomes
    recrates = vcat([[i == 1 ? 0.5 : recrate for i in 1:locus[2]] for locus in chrs]...)
    chrtypes = vcat([[locus[1] for _ in 1:locus[2]] for locus in chrs]...)

    length(recrates) == length(chrtypes) || error()

    nloci = length(chrtypes)

    SexSpecificParameters(
        deme,
        p.population_sizes[deme][sexidx],
        nloci,
        sex,
        p.mutation_rates[sexidx],
        recrates,
        chrtypes,
        p.proportions_of_migrant[deme][sexidx],
        p.fitness_of_migrants[deme][sexidx],
        p.fitness_of_foreign_traits[deme][sexidx],
        p.transmission_rates_of_trait[sexidx]
    )
end

type SexSpecificPopulation
    data::Array{Int, 2} # Chromosome is along row, and locus is along column.
    trait::Array{Int, 1}
    params::SexSpecificParameters
end

function SexSpecificPopulation(params::ModelParameters, deme::Int, sex::SexType)
    ps = SexSpecificParameters(params, deme, sex)
    nloci = length(ps.type_of_locus)
    popsize = ps.population_size
    SexSpecificPopulation(zeros(Int, 2, nloci * popsize), fill(deme, popsize), ps)
end

get(p::SexSpecificPopulation, s::Symbol) = getfield(p.params, s)
get(p::SexSpecificPopulation, s::Symbol, sex::SexType) = getfield(p.params, s)[value(sex)]
get(p::SexSpecificPopulation, s::Symbol, i::Integer) = getfield(p.params, s)[i]
Base.setindex!(p::SexSpecificPopulation, val::Int, ::Type{Val{:trait}}, i::Integer) = p.trait[i] = val
Base.getindex(p::SexSpecificPopulation, ::Type{Val{:trait}}, i::Integer) = p.trait[i]
Base.length(p::SexSpecificPopulation) = p.params.population_size
Base.eachindex(p::SexSpecificPopulation) = eachindex(p.data)
Base.getindex(p::SexSpecificPopulation, i::Integer) = p.data[i]
Base.getindex(p::SexSpecificPopulation, c::Integer, l::Integer) = p.data[c, l]
Base.setindex!(p::SexSpecificPopulation, val::Int, i::Integer) = p.data[i] = val
Base.setindex!(p::SexSpecificPopulation, val::Int, c::Integer, l::Integer) = p.data[c, l] = val

type Population
    data::NTuple{2, SexSpecificPopulation}
end
Base.length(p::Population) = length(p.data)
Base.getindex(p::Population, i::Integer) = p.data[i]
Base.start(p::Population) = 0
Base.next(p::Population, state) = (state += 1; (p.data[state], state))
Base.done(p::Population, state) = length(p.data) == state
Base.eachindex(p::Population) = 1:length(p.data)

function createpopulations(params::ModelParameters)
    ndemes = length(params.population_sizes)
    pops = Array{Population}(ndemes)
    sexes = (Female, Male)
    for i = eachindex(pops)

        spops = Array{SexSpecificPopulation}(2)
        for (j, sex) = enumerate(sexes)
            spops[j] = SexSpecificPopulation(params, i, sex)
        end

        pops[i] = Population(tuple(spops...))
    end
    pops
end

function setinitialgenotype2!(core, p, i)
    for j = 1:2
        p[j,i] = insert!(db(core), GeneRecord(0, nextstate!(core)))
    end
    nothing
end

function setinitialgenotype1!(core, p, i, j)
    p[j,i] = insert!(db(core), GeneRecord(0, nextstate!(core)))
    nothing
end

function setinitialgenotype!(
    core::BasicData,
    p::SexSpecificPopulation,
    ind,
    locus,
    ::Type{Val{Female}},
    ::Union(Type{Val{Autosome}}, Type{Val{XChromosome}}))

    setinitialgenotype2!(core, p, toposition(p, ind, locus))
end

function setinitialgenotype!(
    core::BasicData,
    p::SexSpecificPopulation,
    ind,
    locus,
    ::Type{Val{Female}},
    ::Type{Val{YChromosome}})

    nothing
end

function setinitialgenotype!(
    core::BasicData,
    p::SexSpecificPopulation,
    ind,
    locus,
    ::Type{Val{Female}},
    ::Type{Val{Mitochondrion}})

    setinitialgenotype1!(core, p, toposition(p, ind, locus), 1)
end

function setinitialgenotype!(
    core::BasicData,
    p::SexSpecificPopulation,
    ind,
    locus,
    ::Type{Val{Male}},
    ::Type{Val{Autosome}})

    setinitialgenotype2!(core, p, toposition(p, ind, locus))
end

function setinitialgenotype!(
    core::BasicData,
    p::SexSpecificPopulation,
    ind,
    locus,
    ::Type{Val{Male}},
    ::Type{Val{XChromosome}})

    setinitialgenotype1!(core, p, toposition(p, ind, locus), 1)
end

function setinitialgenotype!(
    core::BasicData,
    p::SexSpecificPopulation,
    ind,
    locus,
    ::Type{Val{Male}},
    ::Type{Val{YChromosome}})

    setinitialgenotype1!(core, p, toposition(p, ind, locus), 2)
end

function setinitialgenotype!(
    core::BasicData,
    p::SexSpecificPopulation,
    ind,
    locus,
    ::Type{Val{Male}},
    ::Type{Val{Mitochondrion}})

    nothing
end

function initialize!(core, pops)
    for pop in pops
        for sexpop in pop
            nl = get(sexpop, :number_of_loci)
            sex = get(sexpop, :sex)
            popsize = length(sexpop)
            for j = 1:popsize, i in 1:nl
                chr = get(sexpop, :type_of_locus, i)
                setinitialgenotype!(core, sexpop, j, i, Val{sex}, Val{chr})
            end
        end
    end
end

function reinitialize!(oldcore, pops)
    oldgdb = db(oldcore)
    core = BasicData()
    gdb = db(core)
    # Preemptively map 0 to 0 as this value has a special meaning.
    smap = Dict{Int, Int}(0 => 0)
    smax = 1
    for pop in pops
        for sexpop in pop
            for idx in eachindex(sexpop.data)
                sexpop[idx] == 0 && continue
                state = oldgdb[sexpop[idx]].state
                if !haskey(smap, state)
                    smap[state] = smax
                    smax += 1
                end
                sexpop[idx] = insert!(gdb, GeneRecord(0, smap[state]))
            end
        end
    end
    core.state = maximum(keys(db(core)))
    core
end

function listgenes!(gids::AbstractArray, pops, locus, idx)
    for pop in pops
        for sexpop in pop
            chrtype = Val{get(sexpop, :type_of_locus, locus)}
            sex = get(sexpop, :sex)
            idx = listgenes!(gids, sexpop, locus, idx, chrtype, Val{sex})
        end
    end
    idx
end

listgenes!(gids::AbstractArray, pops, locus) = listgenes!(gids, pops, locus, 1)

listgenes!(gids::AbstractArray, sexpop, locus, idx, chrtype, sex) = idx

function listgenes!(gids::AbstractArray, sexpop, locus, idx, ::Type{Val{Autosome}}, sex)

    psize = length(sexpop)
    offset = get(sexpop, :number_of_loci)
    for i = 1:psize, chr = 1:2
        gids[idx] = sexpop[chr, (i - 1) * offset + locus]
        idx += 1
    end
    idx
end

function listgenes!(gids::AbstractArray, sexpop, locus, idx,
    ::Type{Val{XChromosome}}, ::Type{Val{Female}})

    psize = length(sexpop)
    offset = get(sexpop, :number_of_loci)
    for i = 1:psize, chr = 1:2
        gids[idx] = sexpop[chr, (i - 1) * offset + locus]
        idx += 1
    end
    idx
end

function listgenes!(gids::AbstractArray, sexpop, locus, idx,
    ::Type{Val{XChromosome}}, ::Type{Val{Male}})

    psize = length(sexpop)
    offset = get(sexpop, :number_of_loci)
    for i = 1:psize
        gids[idx] = sexpop[1, (i - 1) * offset + locus]
        idx += 1
    end
    idx
end

function listgenes!(gids::AbstractArray, sexpop, locus, idx,
    ::Type{Val{YChromosome}}, ::Type{Val{Male}})

    psize = length(sexpop)
    offset = get(sexpop, :number_of_loci)
    for i = 1:psize
        gids[idx] = sexpop[2, (i - 1) * offset + locus]
        idx += 1
    end
    idx
end

function listgenes!(gids::AbstractArray, sexpop, locus, idx,
    ::Type{Val{Mitochondrion}}, ::Type{Val{Female}})

    psize = length(sexpop)
    offset = get(sexpop, :number_of_loci)
    for i = 1:psize
        gids[idx] = sexpop[1, (i - 1) * offset + locus]
        idx += 1
    end
    idx
end

function reproduce!(core, parpop, ppos, chpop, cpos)
    # Maternally derived genes are always stored on the first row, whereas
    # paternally derived genes are on the second row.
    # This function deals with a single offspring, not the entire population of offspring.
    csex = get(chpop, :sex)
    cdeme = get(chpop, :deme)
    pdeme = get(parpop[1], :deme) # both parents are from the same deme.
    for (sexpop, pos) in zip(parpop, ppos)
        psex = get(sexpop, :sex)
        nl = get(sexpop, :number_of_loci)
        cchr = value(psex)
        chr = rand(1:2)
        mut = get(sexpop, :mutation_rate)
        for i = 1:nl
            chrtype = get(sexpop, :type_of_locus, i)
            rrate = get(sexpop, :recombination_rate, i)

            chr, passgene = selectchr(chr, rrate, Val{chrtype}, Val{psex}, Val{csex})
            passgene || continue

            pos = toposition(chpop, cpos, i)
            if rand() < mut
                state = nextstate!(core)
            else
                state = -1
            end
            pgene = sexpop[chr, pos]
            geneid = transmit!(db(core), time(core), pgene, state=state, src=pdeme, dest=cdeme)
            chpop[cchr, pos] = geneid
        end
    end
end

selectchr(chr, rate, ::Type{Val{Autosome}}, psex, csex) = pickchr(chr, rate), true
selectchr(chr, rate, ::Type{Val{XChromosome}}, ::Type{Val{Female}}, csex) =
    pickchr(chr, rate), true
selectchr(chr, rate, ::Type{Val{XChromosome}}, ::Type{Val{Male}}, ::Type{Val{Female}}) = 1, true
selectchr(chr, rate, ::Type{Val{YChromosome}}, ::Type{Val{Male}}, ::Type{Val{Male}}) = 2, true
selectchr(chr, rate, ::Type{Val{Mitochondrion}}, ::Type{Val{Female}}, ::Type{Val{Female}}) =
    1, true
selectchr(chr, rate, chrtype, psex, csex) = chr, false

pickchr(chr, rate) = rand() < rate ? 3 - chr : chr

toposition(pop, ind, locus) = get(pop, :number_of_loci) * (ind - 1) + locus

function evolve!(
    core::BasicData,
    parpops::Array{Population, 1},
    termon::Int,
    tclean::Int)

    info("evolution started on ", now(), ".")

    ndemes = length(parpops) # basically fixed at 2
    nloci = get(parpops[1][1], :number_of_loci)
    loci = get(parpops[1][1], :type_of_locus)

    chpops = deepcopy(parpops)

    ncoals = [0 for _ = 1:nloci]
    lastcoals = [0 for _ = 1:nloci]

    # storage of gene ids. Autosome
    gidstore = Array{Array{Int, 1}, 1}(4)
    # Autosome
    gidstore[1] = zeros(
        Int,
        2 * sum((vcat([[length(pop) for pop in parpop] for parpop in parpops]...)))
    )
    # XChromosome
    gidstore[2] = zeros(
        Int,
        2 * sum([length(parpop[1]) for parpop in parpops]) +
        sum([length(parpop[2]) for parpop in parpops])
    )
    # YChromosome
    gidstore[3] = zeros(
        Int,
        sum([length(parpop[2]) for parpop in parpops])
    )
    # Mitochondrion
    gidstore[4] = zeros(
        Int,
        sum([length(parpop[1]) for parpop in parpops])
    )

    nnewids = sum(map(x -> length(gidstore[value(x)]), loci))

    gen = 0

    # main loop of a simulation
    for gen = core
        for pop in chpops
            for sexpop in pop
                # cdeme and csex indicate the location and sex of offspring (child).
                cdeme = get(sexpop, :deme)
                csex = get(sexpop, :sex)
                # i = toposition(sexpop)
                propmig = get(sexpop, :proportion_of_migrant)
                migfit = get(sexpop, :fitness_of_migrant)
                for org = 1:length(sexpop)
                    while true
                        # decide if an offspring is a native or migrant.
                        pdeme = cdeme
                        if rand() < propmig
                            rand() < migfit || continue
                            pdeme = 3 - cdeme
                        end

                        # If cdeme == pdeme, an offspring is born in pdeme and stay there
                        # for reproduction.
                        # If cdeme != pdeme, an offspring is born in pdeme, and it migrates
                        # to cdeme before reproduction.
                        # Because all offspring are born in pdeme, it is OK to have a parent
                        # to have offspring in both demes.

                        pp = parpops[pdeme]

                        # sample parents
                        mom = rand(1:length(pp[1]))
                        pop = rand(1:length(pp[2]))

                        # sample the trait of a offspring
                        momtrait = pp[1][Val{:trait}, mom]
                        poptrait = pp[2][Val{:trait}, pop]
                        ctrait = trait(
                            get(pp[1], :transmission_rate_of_trait, csex),
                            get(pp[2], :transmission_rate_of_trait, csex),
                            momtrait, poptrait, pdeme)

                        if pdeme != cdeme && ctrait != cdeme
                            tfit = get(sexpop, :fitness_of_foreign_trait)
                            rand() < tfit || continue
                        end

                        reproduce!(core, pp, (mom, pop), sexpop, org)

                        sexpop[Val{:trait}, org] = ctrait
                        break
                    end
                end
            end
        end

        parpops, chpops = chpops, parpops

        for (idx, locus) in enumerate(loci)
            lpos = value(locus)
            listgenes!(gidstore[lpos], parpops, idx)
            anc = mrca(db(core), gidstore[lpos])
            if anc.epoch > lastcoals[idx]
                ncoals[idx] += 1
                lastcoals[idx] = gen
                info("turn over ", ncoals[lpos], " of locus ", idx, " on generation ", gen)
            end
        end
        if termon == minimum(ncoals)
            break
        end
        gen % tclean == 0 && clean!(db(core), nnewids * gen + 1, nnewids * (gen + 1))
    end
    clean!(db(core), nnewids * gen + 1, nnewids * (gen + 1))

    if termon == minimum(ncoals)
        info("evolution terminated by number of turn-overs in gen ", gen, " on ", now(), ".")
    else
        info("evolution terminated by reaching max gen ", gen, " on ", now(), ".")
    end
    parpops, core

end

function simulate(params::ModelParameters, burnin::Int, t::Int, termon::Int, tclean::Int, rep=1)
    info("process started on ", now(), ".")
    core = BasicData()
    pops = createpopulations(params)

    initialize!(core, pops)
    settmax!(core, burnin)
    pops, core  = evolve!(core, pops, -1, tclean)

    datastore = Array{Tuple{typeof(pops), typeof(db(core)), typeof(time(core))}}(0)

    for i = 1:rep
        pops = deepcopy(pops) # clone population that to avoid overwrite.
        core = reinitialize!(core, pops)
        settmax!(core, t)
        pops, core = evolve!(core, pops, termon, tclean)
        push!(datastore, (pops, db(core), time(core)))
    end
    info("process terminated on ", now(), ".")
    datastore
end

end # module


text:
const defaultValences = (
    1,
    0,
    1,
    2,
    3,
    4,
    5,
    -2,
    1,
    0,
    1,
    2,
    3,
    4,
    5,
    6,
    5,
    0,
    1,
    2,
    3,
    4,
    5,
    2,
    2,
    2,
    2,
    2,
    2,
    2,
    3,
    4,
    3,
    6,
    5,
    0,
    1,
    2,
    3,
    4,
    5,
    6,
    2,
    4,
    4,
    2,
    1,
    2,
    3,
    2,
    3,
    4,
    5,
    0,
    1,
    2,
    3,
    4,
    3,
    3,
    3,
    3,
    3,
    3,
    3,
    3,
    3,
    3,
    3,
    3,
    3,
    4,
    5,
    6,
    4,
    4,
    4,
    4,
    3,
    2,
    1,
    2,
    3,
    4,
    5,
    0,
    1,
    2,
    3,
    4,
    5,
    4,
    4,
    4,
)

"""
    asoxide(elm::Element, valences = NeXLCore.defaultValences)

Compute the oxidized form of the specified element using the valences provided in `val`.  By default,
`val = NeXLCore.defaultValences`, a typical set of valences.
"""
function asoxide(
    elm::Element;
    valences = NeXLCore.defaultValences,
    name = nothing,
    atomicweights::Dict{Element,<:AbstractFloat} = Dict{Element,Float64}(),
)
    function buildoxidefraction(elm, val)
        den = gcd(val[z(elm)], -val[z(n"O")])
        return (n"O" => val[z(elm)] ÷ den, elm => -val[z(n"O")] ÷ den)
    end
    function buildoxidename(elm, val)::String
        nnn(n) = n >=10 ? nnn(n÷10)*nnn(n%10) : ("₀","₁","₂","₃","₄","₅","₆","₇","₈","₉" )[n%10+1]
        nn(n) = n==1 ? "" : nnn(n)
        den = gcd(val[z(elm)], -val[z(n"O")])
        ne, no = -val[z(n"O")] ÷ den, val[z(elm)] ÷ den
        return "$(symbol(elm))$(nn(ne))O$(nn(no))"
    end
    name = something(name, buildoxidename(elm, valences))
    return atomicfraction(
        name,
        buildoxidefraction(elm, valences)...,
        atomicweights = atomicweights,
    )
end

"""
    asoxide(elms::Pair{Element, <:AbstractFloat}...; valences = NeXLCore.defaultValences, atomicweights::Dict{Element,<:AbstractFloat} = Dict{Element,Float64}())
    asoxide(elms::Dict{Element, <:AbstractFloat}...; valences = NeXLCore.defaultValences, atomicweights::Dict{Element,<:AbstractFloat} = Dict{Element,Float64}())

Providing the mass-fraction of the consituent elements in `elms`, compute the corresponding amounts of the oxide forms
of the elements.  This can be used to answer the question: If I measure this amount of these elements, what mass fraction of
the oxide-forms of each element does this correspond to?  The example below demonstrates that Albite is 68.74% SiO₂ by mass.
By default, `val = NeXLCore.valences`, a typical set of valences.  See also `obystoichiometry(...)`

Example:

    julia> mat"NaAlSi3O8"
      NaAlSi3O8[Al=0.1029,Na=0.0877,Si=0.3213,O=0.4881]
    julia> asoxide(n"Al"=>0.1029, n"Na"=>0.0877, n"Si"=>0.3213)
      Dict{Material, Float64} with 3 entries:
      SiO₂[O=0.5326,Si=0.4674]  => 0.687366
      Al₂O₃[Al=0.5293,O=0.4707] => 0.194424
      Na₂O[Na=0.7419,O=0.2581]  => 0.118216
    julia> sum(asoxide(n"Al"=>0.1029, n"Na"=>0.0877, n"Si"=>0.3213), name="Albite")
      Albite[Al=0.1029,Na=0.0877,O=0.4881,Si=0.3213]
    julia> asoxide(filter(kv->kv[1]!=n"O", massfraction(mat"NaAlSi3O8")))
      Dict{Material, Float64} with 3 entries:
        SiO₂[O=0.5326,Si=0.4674]  => 0.687401
        Al₂O₃[Al=0.5293,O=0.4707] => 0.194418
        Na₂O[Na=0.7419,O=0.2581]  => 0.118181

"""
function asoxide(
    elms::Pair{Element,<:AbstractFloat}...;
    valences = NeXLCore.defaultValences,
    atomicweights::Dict{Element,<:AbstractFloat} = Dict{Element,Float64}()
) :: Dict{Material, AbstractFloat}
    return  Dict{Material, AbstractFloat}(
        map(elms) do (elm, f)
            @assert elm ≠ n"O" "Don't include the element O in the `elms` argument to asoxide(...)."
            ox = asoxide(elm, valences=valences, atomicweights=atomicweights)
            ox=>(f/ox[elm])
        end
    )
end
function asoxide(
    elms::Dict{Element, <:AbstractFloat};
    valences = NeXLCore.defaultValences,
    name = nothing,
    atomicweights::Dict{Element,<:AbstractFloat} = Dict{Element,Float64}()
)
    return  Dict{Material, AbstractFloat}(
        map(collect(keys(elms))) do elm
            @assert elm ≠ n"O" "Don't include the element O in the `elms` argument to asoxide(...)."
            ox = asoxide(elm, valences=valences, atomicweights=atomicweights)
            ox=>(elms[elm]/ox[elm])
        end
    )
end

"""
    obystoichiometry(elms::Pair{Element, <:AbstractFloat}..., valences = NeXLCore.defaultValences)
    obystoichiometry(elms::Dict{Element, <:AbstractFloat}; valences = NeXLCore.defaultValences)

Compute O-by-stoichiometry from the provided mass fractions of elements.

Example:

    obystoichiometry(n"Mg"=>0.1099, n"Al"=>0.0443, n"Si"=>0.1941, n"Ca"=>0.1034, n"Fe"=>0.0756)
    0.39582340257233467
"""
function obystoichiometry(elms::Dict{Element,<:AbstractFloat}; valences = NeXLCore.defaultValences)
    sum(elms) do (elm, f)
        elm ≠ n"O" ?  f * (-valences[z(elm)] * a(n"O")) / (a(elm) * valences[z(n"O")]) : 0.0
    end
end
function obystoichiometry(elms::Pair{Element,<:AbstractFloat}...; valences = NeXLCore.defaultValences)
    obystoichiometry(Dict(elms...); valences = valences)
end


text:
module TestOneDimensionalRanges

using Test
using MLJBase

mutable struct DummyModel <: Deterministic
    K::Int
    metric::Float64
    kernel::Char
end

dummy_model = DummyModel(4, 9.5, 'k')

mutable struct SuperModel <: Deterministic
    lambda::Float64
    model1::DummyModel
    model2::DummyModel
end

dummy1 = DummyModel(1, 9.5, 'k')
dummy2 = DummyModel(2, 9.5, 'k')
super_model = SuperModel(0.5, dummy1, dummy2)

@testset "constructors" begin
    @test_throws ArgumentError range(dummy_model, :K, lower=Inf,
                                      origin=1, unit=1)
    @test_throws ArgumentError range(dummy_model, :K, upper=-Inf,
                                      origin=1, unit=1)

    @test_throws DomainError range(dummy_model, :K, lower=1)
    @test_throws DomainError range(dummy_model, :K, lower=1, upper=Inf)
    @test_throws DomainError range(dummy_model, :K, upper=1)
    @test_throws DomainError range(dummy_model, :K, upper=1, lower=-Inf)

    @test_throws DomainError range(dummy_model, :K, lower=1, origin=2)
    @test_throws DomainError range(dummy_model, :K, lower=1, upper=Inf,
                                      origin=2)
    @test_throws DomainError range(dummy_model, :K, upper=1, origin=2)
    @test_throws DomainError range(dummy_model, :K, upper=1, lower=-Inf,
                                      origin=2)

    @test_throws DomainError range(dummy_model, :K,
                                      lower=3, unit=0, origin=4)
    @test_throws DomainError range(dummy_model, :K,
                                      lower=3, unit=1, origin=2)

    @test_throws DomainError range(dummy_model, :K, origin=2)
    @test_throws DomainError range(dummy_model, :K, unit=1)
    @test_throws DomainError range(dummy_model, :K)

    @test_throws ArgumentError range(dummy_model, :kernel)

    @test_throws ArgumentError range(dummy_model, :K, values=['c', 'd'])
    @test_throws ArgumentError range(Int, :K, values=['c', 'd'])

    @test range(dummy_model, :K, values=[1, 7]) ==
        range(Int, :K, values=[1, 7])

    z1 = range(dummy_model, :K, lower=1, upper=10)
    @test z1.origin == 5.5
    @test z1.unit == 4.5
    @test z1.scale == :linear

    z2 = range(dummy_model, :K, lower=10, origin=10^6, unit=10^5)
    @test z2.origin == 10^6
    @test z2.unit == 10^5
    @test z2.upper  == Inf
    @test z2.scale == :log

    z3 = range(dummy_model, :K, upper=-10, origin=-10^6, unit=10^5)
    @test z3.origin == -10^6
    @test z3.unit == 10^5
    @test z3.lower  == -Inf
    @test z3.scale == :logminus

    z4 = range(super_model, :lambda, lower=1, upper=10)
    @test z4.origin == 5.5
    @test z4.unit == 4.5
    @test z4.scale == :linear

    z5 = range(dummy_model, :K, origin=10, unit=20)
    @test z5.scale == :linear

    p1 = range(dummy_model, :K, lower=1, upper=10, scale=:log10)
    p2 = range(dummy_model, :kernel, values=['c', 'd'])
    p3 = range(super_model, :lambda, lower=0.1, upper=1, scale=:log2)
    p4 = range(dummy_model, :K, lower=1, upper=3, scale=x->2x)

    # test that you can replace model with type:
    @test z1 == range(Int, :K, lower=1, upper=10)
    @test z4 == range(Float64, :lambda, lower=1, upper=10)
    @test p2 == range(Char, :kernel, values=['c', 'd'])
end

@testset "range constructors for nested parameters" begin
    p1 = range(dummy_model, :K, lower=1, upper=10, scale=:log10)
    q1 = range(super_model, :(model1.K) , lower=1, upper=10, scale=:log10)
    @test iterator(q1, 5) == iterator(p1, 5)
    q2 = range
end

end
true


text:
using HDF5
using Statistics
function plot_acquisition_benchmark(f)
    # error = read(f, "mean_squared_error")
    # error = read(f, "rmse")
    error = read(f, "r2")
    mean_rmse = mean(error, dims = 4) # average over stripes
    std_rmse = std(error, dims = 4)
    # unpack
    μ = @view mean_rmse[:, 1, :] # mean rmse
    σ = @view std_rmse[:, 1, :]
    dμ = @view mean_rmse[:, 2, :] # mean rmse of derivative
    dσ = @view std_rmse[:, 2, :]

    nsample, npolicies = read(f, "nsample"), read(f, "npolicies")
    μ = reshape(μ, nsample, npolicies)
    dμ = reshape(dμ, nsample, npolicies)
    σ = reshape(σ, nsample, npolicies)
    # for each kernel, use length scale with best performance
    # lind = zeros(Int, npolicies)
    # for i in 1:npolicies
    #     A = @view dμ[end, i, :]
    #     # lind[i] = argmin(A)
    #     lind[i] = argmax(A)
    # end
    A = zeros(nsample, npolicies)
    B = zeros(nsample, npolicies)
    for i in 1:npolicies
        @. A[:, i] = dμ[:, i] #, lind[i]]
        @. B[:, i] = dσ[:, i] #, lind[i]]
    end
    n = size(error, 5)
    confidence = maximum(B) / sqrt(n)
    println("confidence = $confidence")
    names = read(f, "policy names")
    plot(xlabel = "# of points", ylabel = "R2") #, yscale = :log10)
    plot!(A, label = permutedims(names), legend = :bottomleft, linewidth = 2) #, linestyle = :auto)
    gui()
end

path = "SARA/NatCom2020/results/"
# file = "inner_acquisition_benchmark_iSARA_short.h5"
# file = "inner_acquisition_benchmark_EQ_short.h5"
file = "inner_acquisition_benchmark_EQ_long.h5"

f = h5open(path * file)

# plotly()
# plotlyjs()
using Plots
pyplot()
# Plots.scalefontsizes(1.5)
plot_acquisition_benchmark(f)
# savefig("inner_kernel_plot.pdf")


text:
using HTTP
using JSON

# will later get secrets from config.json
api_dev_key = ""
api_user_name = ""
api_user_password = ""

# 0=public 1=unlisted 2=private
api_paste_private 	  = "1"
#api_paste_expire_date	= "N" # never expires
api_paste_expire_date	= "10M" # expires in 10 minutes

# "text" applies no formatting or highlighting. See https://pastebin.com/api
# for other options
api_paste_format 		  = "text"

# if an invalid or expired api_user_key is used, an error will spawn. If no
# api_user_key is used, a guest paste will be created
api_user_key 		  	  = ""

url 				= "https:#pastebin.com/api/api_post.php"

# paste name
api_paste_name			= "justmyfilename.php"

# paste content
api_paste_code 	  	= "just some random text you :)"

# get secrets from config.json
function load_config()
  j = nothing
  open("config.json", "r") do io
    j = JSON.parse(io)
  end
  j
end

function login(dev_key, username, password)
  login_url = "https://pastebin.com/api/api_login.php"

  request_body = "api_dev_key=$dev_key"
  request_body *= "&api_user_name=$username"
  request_body *= "&api_user_password=$password"

  println("sending this request body: ", request_body)

  HTTP.request("POST", login_url, [], request_body)
end

function main()
  config = load_config()

  api_dev_key = config["api_dev_key"]
  api_user_name = config["api_user_name"]
  api_user_password = config["api_user_password"]

  response = login(api_dev_key, api_user_name, api_user_password)
  println("this is the body: ", String(response.body))
end


text:
export MtlComputePipelineState

const MTLComputePipelineState = Ptr{MtComputePipelineState}

"""
    MtlComputePipelineState(d::MtlDevice, f::MtlFunction)

Create an object that stores information about the execution parameters of a MtlFunction.
"""
mutable struct MtlComputePipelineState
    handle::MTLComputePipelineState
    device::MtlDevice
end

Base.unsafe_convert(::Type{MTLComputePipelineState}, q::MtlComputePipelineState) = q.handle

Base.:(==)(a::MtlComputePipelineState, b::MtlComputePipelineState) = a.handle == b.handle
Base.hash(q::MtlComputePipelineState, h::UInt) = hash(q.handle, h)

function unsafe_destroy!(cce::MtlComputePipelineState)
    mtRelease(cce.handle)
end

function MtlComputePipelineState(d::MtlDevice, f::MtlFunction)
    handle = @mtlthrows _errptr mtNewComputePipelineStateWithFunction(d, f, _errptr)

    obj = MtlComputePipelineState(handle, d)
    finalizer(unsafe_destroy!, obj)
    return obj
end

## properties

Base.propertynames(o::MtlComputePipelineState) = (
    # identification
    :device, :label,
    # threadgroup attributes
    :maxTotalThreadsPerThreadgroup, :threadExecutionWidth, :staticThreadgroupMemoryLength,
    #other
    #=supportIndirectCommandBuffers=#
)

function Base.getproperty(o::MtlComputePipelineState, f::Symbol)
    if f === :device
        return MtlDevice(mtComputePipelineDevice(o))
    elseif f === :label
        ptr = mtComputePipelineLabel(o)
        ptr == C_NULL ? nothing : unsafe_string(ptr)
    elseif f === :maxTotalThreadsPerThreadgroup
        return Int(mtComputePipelineMaxTotalThreadsPerThreadgroup(o))
    elseif f === :threadExecutionWidth
        return Int(mtComputePipelineThreadExecutionWidth(o))
    elseif f === :staticThreadgroupMemoryLength
        return Int(mtComputePipelineStaticThreadgroupMemoryLength(o))
    else
        getfield(o, f)
    end
end


text:
const T0 = 288.15
const gamma = -6.5e-3
const P0 = 101325
const gravity = 9.8076
const Rd = 287.053

using Interpolations
export interpolation2pressure, interpolate_height2pressure, interpolate_pressure2height, pressure2height, height2pressure, hello, Line

"""
interpolation2pressure(A::ClimArray, pressure::ClimArray, pressure_levels::Vector; heightname=Hei(), extrapolation_bc=NaN )
Return a ClimArray where the vertical coordinate is pressure. Pressure values need to be ascending or descending in order.

`A`: `ClimArray` with some arbitrary height coordinate (height above ground, model levels, etc...) with the variable that should be interpolated
`pressure`: `ClimArray` with pressure values in Pascal and the same coordinates as `A`.
`pressure_levels`: Vector which contains the pressure levels in Pascal that A should be interpolated on.
`vertical_coord` : Name of the vertical coordinate, which is usually `Hei` but can be model levels or something else that does not represent physical height.
`extrapolation_bc`: Extrapolation is set to `NaN` by default, but can be any value or linear extrapolation with `extrapolation_bc = Line()` For other extrapolation methods use the `Interpolations` package.
`descending` : specifies whether pressure is ordered descendingly or ascendingly along the vertical coordinate. Descending by default.
"""
function interpolation2pressure(A::ClimArray, pressure::ClimArray, pressure_levels::Vector; vertical_coord=Hei(), extrapolation_bc=NaN,descending=true )

    dims(A) == dims(pressure) || error("Pressure and Variable Array do not have the same dimensions")
    hasdim(A, vertical_coord) || error("Vertical dimension " + string(vertical_coord) +  " not found in array")

    # The interplation requires ascending pressure values:
    if descending == true
        A = reverse(A,dims=vertical_coord)
        pressure = reverse(pressure,dims=vertical_coord)
    end

    # construct output Array:
    pre = Pre(pressure_levels; metadata = Dict())
    dims_no_height = otherdims(A, vertical_coord)
    out_dims = (dims_no_height...,pre)
    dimension_lengths = length.(out_dims)

    int_array = ClimArray(zeros(eltype(A), Tuple(dimension_lengths)), out_dims ; name = A.name, attrib = A.attrib)
    for i in otheridxs(A, vertical_coord)
        itp = LinearInterpolation(pressure[i],A[i],extrapolation_bc=extrapolation_bc)
        int_array[i] = itp(pressure_levels)

    end

    return int_array
end

"""
    pressure2height(height)
Converts pressure [Pa] to height above ground [m].
The calculation is done under the assumption of hydrostatic equilibrium with a constant lapse rate of -6.5 K/km, starting from a surface temperature T0=288.15 K.
See also equations 39 and 40 in: Berberan-Santos, M. N., Bodunov, E. N., & Pogliani, L. (1997). On the barometric formula. American Journal of Physics, 65(5), 404-412.
"""
function pressure2height(pressure)
    return (T0/gamma) * ( (pressure/P0)^( - gamma*Rd/gravity ) - 1.0 )
end

"""
    height2pressure(height)
Converts height above ground [m] to pressure [Pa].
The calculation is done under the assumption of hydrostatic equilibrium with a constant lapse rate of -6.5 K/km, starting from a surface temperature T0=288.15 K.
See also equations 39 and 40 in: Berberan-Santos, M. N., Bodunov, E. N., & Pogliani, L. (1997). On the barometric formula. American Journal of Physics, 65(5), 404-412.
"""
function height2pressure(height)
    return P0 * exp(- gravity/(Rd * gamma ) * log(1.0 + (height*gamma)/T0 ))
end

"""
    interpolate_height2pressure(A::ClimArray, pressure_levels::Vector; extrapolation_bc=NaN )
Return a `ClimArray` where the vertical coordinate is pressure. Pressure levels are calculated with the `height2pressure` function, based on hydrostatic equilibrium, and then interpolated to the given levels.

`A`: `ClimArray` with height above ground [m] as height coordinate.
`pressure_levels`: Vector which contains the pressure levels in Pascal that `A` should be interpolated on.
`extrapolation_bc`: Extrapolation is set to `NaN` by default, but can be any value or linear extrapolation with `extrapolation_bc = Line()`. For other extrapolation methods use the `Interpolations` package.
"""
function interpolate_height2pressure(A::ClimArray,pressure_levels::Vector; extrapolation_bc=NaN)

    hasdim(A, Hei) || error("Hei() dimension not found in array")

    # The interplation requires ascending coordinates:
    # Pressure has to be ascending, therefore height has to be descending
    if issorted(dims(A,Hei)) ==  true
        A = reverse(A; dims = Hei)
    else
        error("Height is not an ascending coordinate.")
    end

    pressure = height2pressure.(dims(A,Hei()).val)

    # construct output Array:
    pre = Pre(pressure_levels; metadata = Dict())
    dims_no_height = otherdims(A, Hei())
    out_dims = (dims_no_height...,pre)
    dimension_lengths = length.(out_dims)

    int_array = ClimArray(zeros(eltype(A), Tuple(dimension_lengths)), out_dims ; name = A.name, attrib = A.attrib)

    for i in otheridxs(A, Hei())

        itp = LinearInterpolation(pressure,A[i],extrapolation_bc=extrapolation_bc)
        int_array[i] = itp(pressure_levels)

    end
    return int_array

end

"""
    interpolate_pressure2height(A::ClimArray,heights::Vector; extrapolation_bc=NaN)
Return a `ClimArray` where the vertical coordinate is height above ground. Height above ground is calculated with the `pressure2height` function, based on hydrostatic equilibrium, and then interpolated to the given heights.

`A`: ClimArray with pressure [Pa] as height coordinate.
`heights`: Vector which contains the heights that `A` should be interpolated on in meters above ground.
`extrapolation_bc`: Extrapolation is set to `NaN` by default, but can be any value or linear extrapolation with `extrapolation_bc = Line()`. For other extrapolation methods use the `Interpolations` package.
"""
function interpolate_pressure2height(A::ClimArray,heights::Vector; extrapolation_bc=NaN)

    hasdim(A, Pre) || error("Pre() dimension not found in array")

    height = pressure2height.(dims(A,Pre()).val)

    # construct output Array:
    hei = Hei(heights; metadata = Dict())
    dims_no_pre = otherdims(A, Pre())
    out_dims = (dims_no_pre...,hei)
    dimension_lengths = length.(out_dims)

    int_array = ClimArray(zeros(eltype(A), Tuple(dimension_lengths)), out_dims ; name = A.name, attrib = A.attrib)

    for i in otheridxs(A, Pre())

        itp = LinearInterpolation(height,A[i],extrapolation_bc=extrapolation_bc)
        int_array[i] = itp(heights)

    end
    return int_array

end


text:
# Copyright 2020 Digital Domain 3.0
#
# Licensed under the Apache License, Version 2.0 (the "Apache License")
# with the following modification; you may not use this file except in
# compliance with the Apache License and the following modification to it:
# Section 6. Trademarks. is deleted and replaced with:
#
# 6. Trademarks. This License does not grant permission to use the trade
#    names, trademarks, service marks, or product names of the Licensor
#    and its affiliates, except as required to comply with Section 4(c) of
#    the License and to reproduce the content of the NOTICE file.
#
# You may obtain a copy of the Apache License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the Apache License with the above modification is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the Apache License for the specific
# language governing permissions and limitations under the Apache License.

module DiscreteDifferentialGeometry

#!me check that the signs of the cell differentials satisfy 𝑑_cellⁿ⁻ᵏ == (-1)ᵏ(𝑑_primalᵏ⁻¹)ᵀ
#!me might not have the sign right

"""
Ωᵢ : k-form (k == i) domain of the mesh.  i.e. Ω₀ is on verts Ω₁ is on edges, Ω₂ is on faces
Ω⋆ᵢ : k-form (k == i) dual domain of the mesh.  i.e. Ω⋆₀ is on dual verts Ω₁ is on dual edges, Ω₂ is on dual faces
⋆ᵢ : hodge star operator
𝑑ᵢ : differential operator

   Ω₀---𝑑₀-→Ω₁----𝑑₁-→Ω₂
   ↑        ↑         ↑
 ⋆₀|⋆₀⁻   ⋆₁|⋆₁⁻    ⋆₂|⋆₂⁻
   ↓        ↓         ↓
   Ω⋆₂←-𝑑₁--Ω⋆₁←-𝑑₀---Ω⋆₀

"""

using Multivectors
using LinearAlgebra, SparseArrays
using HalfEdges
using HalfEdges: incidence, area, nvertices, faces, isboundary, head, edge, next, opposite
import Multivectors: ⋆

include("KForms.jl")

export
DiscreteDifferentialOperator,
PrimalKForm,
CellKForm,
array,
δ, codifferential,
Δ, laplace_beltrami,
Δcell, laplace_beltrami_cell,
𝑑, differential,
⋆, star,
poisson,
constrain_system,
constrain_rhs,
interpolate,
CircumcentricDual,
BarycentricDual

struct Curve
  indices::Vector{Int}
  P
end

partial = (f::Function,y...)->(z...)->f(y...,z...)
∂(f::Function,y...) = partial(f,y...)

abstract type DifferentialOperator end
abstract type DiscreteDifferentialOperator{N, A} <: DifferentialOperator end

struct DiscreteDifferential{N, A} <: DiscreteDifferentialOperator{N, A}
  D::A
end

struct DiscreteCellDifferential{N, A} <: DiscreteDifferentialOperator{N, A}
  D::A
end

struct DiscreteHodgeStar{N,A} <: DiscreteDifferentialOperator{N, A}
  D::A
end

struct DiscreteHodgeStarInv{N,A} <: DiscreteDifferentialOperator{N, A}
  D::A
end

struct ZeroOperator{N,A} <: DiscreteDifferentialOperator{N,A}
  D::A
end

"""
    DiscreteKForm{N,A}

A discrete k-form of grade N, represented by an indexible container A ( usually a Vector )
A 0-form (N==0) is a k-form with a value per vertex
A 1-form (N==1) is a k-form with a value per edge
A 2-form (N==2) is a k-form with a value per face
"""
struct PrimalKForm{N,A} <: DiscreteDifferentialOperator{N,A}
  D::A
end

PrimalKForm(K::A, N) where A = PrimalKForm{N,A}(K)
PrimalKForm{N}(K::A) where {A<:Diagonal, N} = PrimalKForm{N,A}(K)
PrimalKForm{N}(K::AbstractVector) where {N} = PrimalKForm{N}(Diagonal(K))
PrimalKForm(K::V, N) where V<:AbstractVector = PrimalKForm{N}(Diagonal(K))

"""
    DiscreteCellKForm{N,A}

A discrete k-form of grade N, represented by an indexible container A ( usually a Vector )
A 0-form (N==0) is a k-form with a value per dual face ( point at circumcentre, on cell boundary )
A 1-form (N==1) is a k-form with a value per dual edge ( boundary of one ring cell )
A 2-form (N==2) is a k-form with a value per dual vertex ( one ring cell around vertex )

"""
struct CellKForm{N,A} <: DiscreteDifferentialOperator{N,A}
  D::A
end

CellKForm(K::A, k) where A = CellKForm{k,A}(K)
CellKForm{N}(K::A) where {A<:Diagonal, N} = CellKForm{N,A}(K)
CellKForm{N}(K::AbstractVector) where {N} = CellKForm{N}(Diagonal(K))
CellKForm(K::V, k) where V<:AbstractVector = CellKForm{k}(Diagonal(K))

array(D::DiscreteDifferentialOperator) = D.D
array(D::PrimalKForm) = D.D
array(D::CellKForm) = D.D
array(P::T) where {V<:AbstractArray,T<:AbstractArray{V}} = reduce(hcat, P)'

isdual( ::DiscreteHodgeStar ) = true
isdual( ::DiscreteHodgeStarInv ) = false
isdual( ::DiscreteCellDifferential ) = true
isdual( ::DiscreteDifferential ) = false
isdual( ::CellKForm ) = true
isdual( ::PrimalKForm ) = false

Base.:*( 𝑑̂n::D, A::T ) where {D<:DiscreteDifferentialOperator, T<:AbstractArray} = 𝑑̂n.D*A
Base.:*( A::T, 𝑑̂n::D) where {D<:DiscreteDifferentialOperator, T<:AbstractArray} = A*𝑑̂n.D
Base.iszero( 𝑑̂n::D ) where {D<:DiscreteDifferentialOperator, T<:AbstractArray} = iszero(𝑑̂n.D)
Base.:*( 𝑑̂n::D, s::T ) where {D<:DiscreteDifferentialOperator, T<:Real} = D(𝑑̂n.D*s)
Base.:*( s::T, 𝑑̂n::D ) where {D<:DiscreteDifferentialOperator, T<:Real} = D(s*𝑑̂n.D)
Base.:/( 𝑑̂n::D, s::T ) where {D<:DiscreteDifferentialOperator, T<:Real} = D(𝑑̂n.D/s)

Base.:+( n::Union{D, ZeroOperator}, m::Union{D, ZeroOperator} ) where {D<:DiscreteDifferentialOperator} = D(n.D+m.D)
Base.:-( n::Union{D, ZeroOperator}, m::Union{D, ZeroOperator} ) where {D<:DiscreteDifferentialOperator} = D(n.D-m.D)

struct MultiDifferential <: DifferentialOperator
  MD::Vector{DifferentialOperator}
end

array(D::MultiDifferential) = sum(D.MD)

Base.:+(a::DifferentialOperator, b::DifferentialOperator) = MultiDifferential( [a,b] )
Base.:-(a::DifferentialOperator, b::DifferentialOperator) = MultiDifferential( [a,-b] )
Base.:-(a::MultiDifferential) = MultiDifferential(-1.0.*(a.MD))
Base.:+(a::MultiDifferential, b::DifferentialOperator) = MultiDifferential( vcat(a.MD, b) )
Base.:+(a::DifferentialOperator, b::MultiDifferential) = MultiDifferential( vcat(b.MD, b) )
Base.:+(a::MultiDifferential, b::MultiDifferential) = MultiDifferential( vcat(a.MD, b.MD) )

Base.:(==)( a::DifferentialOperator, b::DifferentialOperator ) = array(a) == array(b)

"""the cotangent of the angle between two vectors"""
cotan( v1::T, v2::T ) where {T<:AbstractArray} = (v1⋅v2)/norm(v1×v2)

"""
  across( topo, P, heh )

the triangle vertex across from the halfedge.
"""
function across( topo::Topology, P, heh::HalfEdgeHandle ) where {PT}
  if isboundary(topo, heh)
    nothing
  else
    @inbounds P[head(topo, next(topo, heh))]
  end
end

"""
vector cotangent of angle at Pₐ, ccw from P₁ to P₂
"""
function cotanweight( P₁::T, P₂::T, Pₐ::T ) where T
  cotan( P₁-Pₐ, P₂-Pₐ )
end

cotanweight( P₁::T, P₂::T, P₃::Nothing ) where T =  zero(eltype(T))

"""
  cotanweight( topo, P, h )

calculate cotangent weight across an edge.
0.5*(cot(α) + cot(Β)), where α,Β are the angles of edges opposite our halfedge
"""
function cotanweight( topo::Topology, P, he::HalfEdgeHandle )
  αᵢ, αⱼ = edge(topo, he)
  Pαᵢ = P[αᵢ]
  Pαⱼ = P[αⱼ]
  cotα = cotanweight(Pαᵢ, Pαⱼ, across(topo, P, he))
  cotΒ = cotanweight(Pαⱼ, Pαᵢ, across(topo, P, opposite(topo, he)))
  0.5*(cotα + cotΒ)
end

"""
Barycentric dual vertex area. One third of the summed areas of triangles incident to given vertex
"""
function vertexarea( mesh, P, i::VertexHandle )
  ring = OneRing(mesh,i)
  iring = Iterators.filter(heh->!isboundary(mesh,heh),ring)
  mapreduce(heh->area(mesh,P,heh)/3.0,+,iring)
end

##===== Operators implementations =======##

#======= MESHES =======#

struct CircumcentricDual end
struct BarycentricDual end

⋆(topo::Topology, P, D::DO, ::DT) where {DO <: DifferentialOperator, DT<:Union{CircumcentricDual, BarycentricDual}} = ⋆(topo, P, D)

# From Primal 0

""" ⋆(topo, P)  Ω₀->⋆₀->Ω⋆₂.  A map to dual 2-forms (dual vertices) """
function ⋆(topo::Topology, P, ::CircumcentricDual )
  aa = map(VertexHandle.(1:nvertices(topo))) do vh
    reduce(Iterators.filter(!∂(isboundary, topo), (OneRing(topo, vh))); init=0.0) do acc, h
      i,j = edge(topo, h)
      k = head(topo, next(topo, h))
      Pi, Pj, Pk = P[[i, j, k]]
      Pki = Pi-Pk; Pij = P[j]-P[i]
      acc+0.125*((Pki⋅Pki)*cotanweight( Pk, Pi, Pj) +
           (Pij⋅Pij)*cotanweight( Pi, Pj, Pk))
    end
  end |> Diagonal
  DiscreteHodgeStar{0,typeof(aa)}(aa)
end

""" ⋆(topo, P)  Ω₀->⋆₀->Ω⋆₂.  A map to dual 2-forms (dual vertices) """
function ⋆(topo::Topology, P, ::BarycentricDual )
  aa = map(VertexHandle.(1:nvertices(topo))) do vh
    vertexarea( topo, P, vh )
  end |> Diagonal
  DiscreteHodgeStar{0,typeof(aa)}(aa)
end

""" ⋆(topo, P)   Ω₀->⋆₀->Ω⋆₂.  A map to dual 2-forms (dual vertices) """
⋆(topo::Topology, P) = ⋆(topo, P, CircumcentricDual())

const Omega0{A} = Union{DiscreteHodgeStarInv{0,A}, PrimalKForm{0,A}}

""" Ω₀->𝑑₀->Ω₁ """
function 𝑑( topo::Topology )
  𝑑0 = incidence( topo, VertexHandle, EdgeHandle, true )
  DiscreteDifferential{0, typeof(𝑑0)}(𝑑0)
end

""" ⋆₀⁻¹->𝑑₀->Ω₁ """
function 𝑑( topo::Topology, D::Omega0 )
  𝑑0star = incidence( topo, VertexHandle, EdgeHandle, true )*array(D)
  DiscreteDifferential{0, typeof(𝑑0star)}(𝑑0star)
end

"""
        Ω₀->⋆₀->Ω⋆₂
  ⋆₀⁻¹->Ω₀->⋆₀->Ω⋆₂

A map to 2-forms (vertex cell).
"""
function ⋆(topo::Topology, P, D::Omega0{A}) where A
  aa = (⋆(topo, P)).D*D.D
  DiscreteHodgeStar{0,typeof(aa)}(aa)
end

# From Primal 1

const Omega1{A} = Union{DiscreteDifferential{0,A}, DiscreteHodgeStarInv{1,A}, PrimalKForm{1,A}}

"""
    𝑑( topo, Ω₁ )

     Ω₁->𝑑₁->Ω₂
 𝑑₀->Ω₁->𝑑₁->Ω₂
⋆⁻¹->Ω₁->𝑑₁->Ω₂

A map to 2-forms (faces)
"""
function 𝑑( topo::Topology, D::Omega1{A} ) where A
  𝑑1 = incidence( topo, EdgeHandle, FaceHandle, true )*D.D
  DiscreteDifferential{1, typeof(𝑑1)}(𝑑1)
end

"""
    ⋆( topo, Ω₁ )

     Ω₁->⋆₁->Ω⋆₁
 𝑑₀->Ω₁->⋆₁->Ω⋆₁
⋆⁻₁->Ω₁->⋆₁->Ω⋆₁

A map to 1-forms (cell edges)
"""
function ⋆( topo::Topology, P, D::Omega1{A} ) where {A}
  stard0 = (map(∂(cotanweight, topo, P), UniqueHalfEdges(topo)) |> Diagonal)*D.D
  DiscreteHodgeStar{1, typeof(stard0)}(stard0)
end

# From Primal 2

const Omega2{A} = Union{DiscreteDifferential{1,A}, DiscreteHodgeStarInv{2,A}, PrimalKForm{2,A}}

"""
    ⋆( topo, Ω₂ )

     Ω₂->⋆₂->Ω⋆₀
 𝑑₁->Ω₂->⋆₂->Ω⋆₀
⋆⁻₂->Ω₂->⋆₂->Ω⋆₀

A map to 0-forms (cell centred vertex)
"""
function ⋆( topo::Topology, P, D::Omega2{A} ) where {A}
  stard1 = ((inv∘∂(area, topo, P)).((((x->Polygon(topo, x)).(faces(topo))))) |> Diagonal)*array(D)
  DiscreteHodgeStar{2, typeof(stard1)}(stard1)
end

# From Dual 2

const OmegaStar2{A} = Union{DiscreteHodgeStar{0,A}, DiscreteCellDifferential{1,A}, CellKForm{2,A}}

""" ->⋆₀->⋆₀⁻¹->Ω₀ """
function ⋆(topo::Topology, P, D::OmegaStar2{A}) where A
  aa = inv(⋆(topo, P).D)*D.D
  DiscreteHodgeStarInv{0,typeof(aa)}(aa)
end

#==
""" ->𝑑₁ᵀ->⋆₀⁻¹->Ω₀ """
function ⋆( topo::Topology, P, 𝑑̂1t::D1, dual_area_type = CircumcentricDual()) where {A, D1<:DiscreteDifferentialT{1, A}}
  stard1t = (⋆(topo, P, dual_area_type) |> inv∘array)*𝑑̂1t.D
  DiscreteHodgeStarInv{0, typeof(stard1t)}(stard1t)
end
==#

# From Dual 1

const OmegaStar1{A} = Union{DiscreteCellDifferential{0,A}, DiscreteHodgeStar{1,A}, CellKForm{1,A}}

"""
    𝑑( topo, Ω⋆₁ )

    Ω⋆₁->𝑑₁->Ω⋆₂
𝑑₀->Ω⋆₁->𝑑₁->Ω⋆₂
⋆₁->Ω⋆₁->𝑑₁->Ω⋆₂
"""
function 𝑑( topo::Topology, D::OmegaStar1{A} ) where A
  𝑑1t = adjoint(incidence( topo, VertexHandle, EdgeHandle, true ))*D.D
  DiscreteCellDifferential{1, typeof(𝑑1t)}(𝑑1t)
end

"""
    ⋆( topo, Ω⋆₁ )

    Ω⋆₁->⋆₁⁻¹->Ω₁
𝑑₀->Ω⋆₁->⋆₁⁻¹->Ω₁
⋆₁->Ω⋆₁->⋆₁⁻¹->Ω₁
"""
function ⋆( topo, P, D::OmegaStar1 )
  starstar = (map(∂(cotanweight, topo, P), UniqueHalfEdges(topo)) |> inv∘Diagonal)*D.D
  DiscreteHodgeStarInv{1, typeof(starstar)}(starstar)
end

# From Dual 0
const OmegaStar0{A} = Union{DiscreteHodgeStar{2,A}, CellKForm{0,A}}

""" ->⋆₂->⋆₂⁻¹->Ω₂ """
function ⋆( topo::Topology, P, D::OmegaStar0 )
  starstar = ((∂(area, topo, P)).((((x->Polygon(topo, x)).(faces(topo))))) |> Diagonal)*array(D)
  DiscreteHodgeStarInv{2, typeof(starstar)}(starstar)
end

"""
    𝑑( topo, s )

Ω₂⋆->𝑑₀->Ω⋆₁.  A map to dual edges
"""
function 𝑑( topo::Topology, D::OmegaStar0{A}) where {A}
  dstar = adjoint(incidence( topo, EdgeHandle, FaceHandle, true ))*array(D)
  DiscreteCellDifferential{0, typeof(dstar)}(dstar)
end
#==
"""
    𝑑( topo, s )

Ω₂⋆->𝑑₀->Ω⋆₁.  A map to dual edges
"""
function 𝑑( topo::Topology, s::H2 ) where {A, H2<:DiscreteHodgeStar{2, A}}
  d0tstar = adjoint(incidence( topo, EdgeHandle, FaceHandle, true ))*array(s)
  DiscreteCellDifferential{0, typeof(d0tstar)}(d0tstar)
end
==#

#==
""" ->⋆₀⁻¹->⋆₀->Ω₀ """
function ⋆(topo::Topology, P, D::DiscreteHodgeStarInv{0,A}) where A
  aa = (⋆(topo, P)).D*D.D
  DiscreteHodgeStar{0,typeof(aa)}(aa)
end
==#

zero_operator(::Type{A}) where A = (𝑑z = zero(eltype(A))*I; ZeroOperator{0, typeof(𝑑z)}(𝑑z))

""" 𝑑(topo, Z) differential of zero is zero """
𝑑( topo::Topology, Z::ZeroOperator{N,A} ) where {N,A} = zero_operator(A)

""" 𝑑(topo, 𝑑₀ᵗ) differential of 𝑑₀ᵗ which is a map to dual vertices is zero, since dual simplex of degree 3 doesn't exist on 𝑅² manifold """
𝑑( topo::Topology, D::DiscreteCellDifferential{1,A} ) where A = zero_operator(A)

""" 𝑑(topo, 𝑑₁) differential of 𝑑₁ which is a map to faces is zero, since simplex of degree 3 doesn't exist on 𝑅² manifold """
𝑑( topo::Topology, D::DiscreteDifferential{1,A} ) where A = zero_operator(A)

""" 𝑑(topo, ⋆₀) differential of ⋆₀ which is a map to dual vertices is zero, since dual simplex of degree 3 doesn't exist on 𝑅² manifold """
𝑑( topo::Topology, H::DiscreteHodgeStar{0,A} ) where A = zero_operator(A)

""" 𝑑(topo, ⋆₂⁻¹) differential of ⋆₂⁻¹ which is a map to faces is zero, since simplex of degree 3 doesn't exist on 𝑅² manifold """
𝑑( topo::Topology, H::DiscreteHodgeStarInv{2,A} ) where A = zero_operator(A)

⋆( topo::Topology, P, Z::ZeroOperator ) = Z

"""
    δ(topo, P)

Discrete coderivative
"""
δ(topo::Topology, P) = ⋆(topo, P, 𝑑(topo, ⋆(topo, P)))

δ(topo::Topology, P, D::DO) where DO<:DifferentialOperator = ⋆(topo, P, 𝑑(topo, ⋆(topo, P, D)))

δ(topo::Topology, P, bcd::BarycentricDual) =
  ⋆(topo, P, 𝑑(topo, ⋆(topo, P, bcd)), bcd)

δ(topo::Topology, P, D::DO, bcd::BarycentricDual) where DO<:DifferentialOperator =
  ⋆(topo, P, 𝑑(topo, ⋆(topo, P, D, bcd)), bcd)

"""
    Δ(topo, P)

Discrete Laplace-Beltrami operator
"""
Δ(topo::Topology, P) = δ(topo, P, 𝑑(topo)) + 𝑑(topo, δ(topo, P))
Δ(topo::Topology, P, bcd::BarycentricDual) = δ(topo, P, 𝑑(topo), bcd) + 𝑑(topo, δ(topo, P, bcd))

"""
    Δcell(topo, P)

The symmetric Laplace-Beltrami operator.  ⋆Δcell = Δ.
The cell Δ differs from the (strong) Δ in that cell version maps to the dual simplex.

i.e. Δcell(topo, P) -> Ω⋆₀ and Δ(topo, P) -> Ω₀

The rhs must also be on dual. You must apply ⋆ to the rhs.
Solve: Δ(topo, P) = b  -----> Solve: Δcell(topo, P) = ⋆(topo,P)*b

Useful for constructing symmetric operators.
"""
Δcell(topo, P) = 𝑑(topo, ⋆(topo, P, 𝑑(topo)))
# Note: we have dropped the term corresponding to ⋆𝑑δ(topo, P) as it is 0 when starting from vertices

# Optimized Laplace-Beltrami operator given face and vertex lists
"""
local Laplace matrix for a triangle
"""
function Δ( i::PT, j::PT, k::PT ) where { T, PT<:AbstractVector{T} }
  ij = j-i
  jk = k-j
  ki = i-k
  a = cotan(ij,-ki)*T(0.5)
  b = cotan(jk,-ij)*T(0.5)
  c = cotan(ki,-jk)*T(0.5)
  Matrix3( [b+c -c -b; -c c+a -a;-b -a a+b] )
end

"""
    Δcell(F, P)

Laplace-Beltrami operator built from list of triangle indices and points.  Optimized for performance.

"""
function Δcell( F::Vector{T}, P::Vector{PT} ) where {T<:Union{AbstractArray,NTuple{3}}, FT, PT<:AbstractVector{FT}}
  n = length(F)*9
  indexI = Vector{Int}(undef,n)
  indexJ = Vector{Int}(undef,n)
  element = Vector{FT}(undef,n)
  s = 1
  @inbounds for ijk in F
    i,j,k = ijk
    Pᵢ = P[i]
    Pⱼ = P[j]
    Pk = P[k]
    ij = Pⱼ-Pᵢ
    jk = Pk-Pⱼ
    ki = Pᵢ-Pk
    a = cotan(ij,-ki)*FT(0.5)
    b = cotan(jk,-ij)*FT(0.5)
    c = cotan(ki,-jk)*FT(0.5)
    #L[i,i] += b+c
    indexI[s] = i; indexJ[s] = i; element[s] = b+c; s+=1
    #L[i,j] -= c
    indexI[s] = i; indexJ[s] = j; element[s] = -c; s+=1
    #L[i,k] -= b
    indexI[s] = i; indexJ[s] = k; element[s] = -b; s+=1
    #L[j,i] -= c
    indexI[s] = j; indexJ[s] = i; element[s] = -c; s+=1
    #L[j,j] += c+a
    indexI[s] = j; indexJ[s] = j; element[s] = c+a; s+=1
    #L[j,k] -= a
    indexI[s] = j; indexJ[s] = k; element[s] = -a; s+=1
    #L[k,i] -= b
    indexI[s] = k; indexJ[s] = i; element[s] = -b; s+=1
    #L[k,j] -= a
    indexI[s] = k; indexJ[s] = j; element[s] = -a; s+=1
    #L[k,k] += a+b
    indexI[s] = k; indexJ[s] = k; element[s] = a+b; s+=1
  end
  sparse(indexI, indexJ, element)
end

Δcell( F::Vector{IT}, P ) where {IT<:Integer} = Δ(Iterators.partition( F, 3 ) |> collect, P)

#======= SHARP, FLAT =====#

#♯₀(topo, P) = ⋆(topo, P, ⋆(topo, P))
#

#======= ALIASES =========#

const codifferential = δ
const laplace_beltrami = Δ
const laplace_beltrami_cell = Δcell
const differential = 𝑑
const star = ⋆

#======= CURVES =======#

"""
    𝑑(curve)

Ω₀->𝑑₀->Ω₁
map from vertex to edges
"""
function 𝑑( c::Curve )
  V = length(c.P)
  E = length(c.indices)-1
  T = eltype(c.P)
  𝑑0 = sparse( reduce(vcat, [[i,i] for i in 1:E]),
               reduce( vcat, [ [ij...] for ij in zip(c.indices[1:end-1], c.indices[2:end])]),
               reduce(vcat, [[T(1),-T(1)] for i in 1:E]),
               E, V)
  DiscreteDifferential{0, typeof(𝑑0)}(𝑑0)
end

#======= VECTOR FIELDS ======#
# some useful details here: https://www.math.arizona.edu/~agillette/research/decNotes.pdf
# Discrete Vector Fields are valued at either vertices or dual triangles.  i.e. 0-simplices
# The dual 0-simplex is the circumcentre of a primal triangle ( 2-simplex ). i.e. a dual triangle
# This is not really explained in the CMU course at the point I've read up to, but is used in their code
# Also geodesic heat method uses it for computing the gradient.
# So heat values are stored on primal vertices and then interpolated to dual 0-simplex... hmmm
# could you instead compute heat values on dual 0-simplex and then apply differential?
# actually circumetric location is used because the tangent plane is well defined on a face.
# ϕᵢ𝑑ϕⱼ - ϕⱼ𝑑ϕᵢ

baryc(p::K, (b, c)::KK, A⁻¹::KB) where {K<:KVector, KK<:Tuple{KVector,KVector}, KB<:Union{KVector, Blade}} =
  grade(((c-b)∧(p-b))*A⁻¹, 0)

#==
ϕ₁(p::K, a::K,b::K,c::K, A⁻¹::K) where K<:KVector = ϕᵢ(p, (b,c), A⁻¹)
ϕ₂(p::K, a::K,b::K,c::K, A⁻¹::K) where K<:KVector = ϕᵢ(p, (c,a), A⁻¹)
ϕ₃(p::K, a::K,b::K,c::K, A⁻¹::K) where K<:KVector = ϕᵢ(p, (a,b), A⁻¹)
==#

function barycentric_coords( p::K ) where {K<:KVector}
  A = (b-a)∧(c-a)
  ((c-b)∧(p-b)/A, (a-c)∧(p-c)/A, (b-a)∧(p-a)/A)  # a tuple of coords
end

function ϕᵢ(edge_bc::KK, A⁻¹::KB) where {K<:KVector, KK<:Tuple{KVector,KVector}, KB<:Union{KVector, Blade}}
  A₁₂, A₁₃, A₂₃ = coords(A⁻¹)
  b,c = edge_bc
  e₁, e₂, e₃ = basis_1blades(b[1])
  cbx, cby, cbz = coords(c-b)
  ∂ϕᵢ∂x =  cby*(A₁₂) + cbz*(A₁₃)
  ∂ϕᵢ∂y =  cbz*(A₂₃) - cbx*(A₁₂)
  ∂ϕᵢ∂z = -cbx*(A₁₃) - cby*(A₂₃)

  ZForm{3}( (x, y, z)->baryc(e₁(x)+e₂(y)+e₃(z), edge_bc, A⁻¹),
           [(x, y, z)->∂ϕᵢ∂x,
            (x, y, z)->∂ϕᵢ∂y,
            (x, y, z)->∂ϕᵢ∂z])
end

function ϕᵢ(a::KVector, b::KVector, c::KVector)
  A = (b-a)∧(c-a)
  A⁻¹ = sortbasis(inv(A)) # ordered e₁₂, e₁₃, e₂₃
  (ϕᵢ((b, c), A⁻¹), ϕᵢ((c, a), A⁻¹), ϕᵢ((a, b), A⁻¹))
end

𝑑ϕ( ϕ, p ) = map(∂ϕ->∂ϕ(coords(p)...), ϕ.gradient) .* basis_1blades(p)

"""
    interpolate( (a,b,c), edgevals, p)

Interpolate the discrete values on the edges of a triangle with points a,b,c to a vector at location p.

edgevals: values on edges should be ordered so the value at index i corresponds to the edge across from vertex i.
"""
function interpolate( (a,b,c)::Tri, edgeval::V, p::X ) where
  {Tri<:Union{Tuple{KVector,KVector,KVector}, AbstractVector{<:KVector}}, V<:Union{AbstractArray, Tuple}, X<:KVector}

  ϕ₁, ϕ₂, ϕ₃ = ϕᵢ(a, b, c)
  xyz = coords(p)

  wab = ϕ₁.f(xyz...)*𝑑ϕ(ϕ₂, p) - ϕ₂.f(xyz...)*𝑑ϕ(ϕ₁, p)
  wbc = ϕ₂.f(xyz...)*𝑑ϕ(ϕ₃, p) - ϕ₃.f(xyz...)*𝑑ϕ(ϕ₂, p)
  wca = ϕ₃.f(xyz...)*𝑑ϕ(ϕ₁, p) - ϕ₁.f(xyz...)*𝑑ϕ(ϕ₃, p)

  v = edgeval[3]*wab + edgeval[1]*wbc + edgeval[2]*wca
  KVector(v)
end

interpolate(a,b,c,eva,evb,evc,p) = interpolate((a,b,c),(eva,evb,evc),p)

#======= BOUNDARY ======#

"""
    constrain( L, pin )

Add Diriclete boundary conditions to Linear Operator L.
return tuple with reduced L and columns of L representing constrained dofs
"""
function constrain_system( L_in::S, pin ) where {S<:AbstractArray}
  L = copy(L_in)
  I_ = sparse(one(eltype(S))*I,size(L)...)
  L_constrained = L[:,pin]
  L[:,pin] = I_[:,pin]
  L[pin,:] = I_[pin,:]

  (L,L_constrained)
end

"""
  add boundary conditions to rhs of linear system
"""
function constrain_rhs(rhs, u, pin, L_constrained)
  b = rhs - L_constrained*u[pin,:]

  for ic in pin
    b[ic,:] = u[ic,:]
  end
  b
end

"""
    poisson(L, b, u_boundary, boundary_index)

Solve the Poisson problem Lu = b, where L is a Laplace-beltrami operator with Dirichlet boundary specified by u_boundary at locations boundary_index.
"""
function poisson(L, b, u_boundary, boundary_index)
  (Lₒ, Lpin) = constrain_system(L, boundary_index)

  rhs = constrain_rhs(b, u_boundary, boundary_index, Lpin)

  Lₒ\rhs
end

"""
    apply(𝑑̂ₐ, α, P, k)

Apply discrete differential operator 𝑑̂ₐ to k-form α at points P with k-vector k.
"""
function apply(𝑑̂ₐ::D, α::KF, P::A, k::K) where{D<:DiscreteDifferentialOperator,
                                        F<:ZForm, KF<:KVector{F},
                                        T, V<:AbstractVector{T}, A<:AbstractVector{V},
                                        K<:KVector{T}}
  ϕₐ = map(P) do pᵢ
    apply(α, pᵢ, k)
  end
  apply(𝑑̂ₐ, ϕₐ)
end

"""
    apply(𝑑̂ₐ, ϕₐ)

Apply discrete differential operator against ϕₐ samples taken at k-simplicies.
"""
apply(𝑑̂ₐ::D, ϕₐ::T) where {D<:DiscreteDifferentialOperator, T<:AbstractVector} = 𝑑̂ₐ.D*ϕₐ

# . It is helpful, as a matter of notation first, to consider differentiation as an abstract operation that accepts a function and returns another function

#==
𝑑̂(D::A) where A<:AbstractArray = (ϕ::Fu)->DiscreteDifferentialOperator()
⋆(DiscreteDifferentialOperator

abstract type Simplex end
struct Vertex <: Simplex end
struct Edge <: Simplex end
struct Face <: Simplex end

simplex_trait(::Type{VertexHandle}) = Vertex()
simplex_trait(::Type{HalfEdgeHandle}) = Edge()
simplex_trait(::Type{EdgeHandle}) = Edge()
simplex_trait(::Type{FaceHandle}) = Face()

𝑑( topo::Topology, ::Type{VertexHandle} ) = 𝑑̂( incidence(topo, VertexHandle, EdgeHandle, true) )
𝑑( topo::Topology, ::Type{EdgeHandle} ) = 𝑑̂( incidence(topo, EdgeHandle, FaceHandle, true) )
𝑑( topo::Topology, ::Type{FaceHandle} ) = 𝑑̂( [0.0] )

apply(𝑑̂₀(ϕ), P_flat_centred)

𝑑(::Vertex, mesh::T) where T =
𝑑(::Edge, mesh::T) where T = (u, v, w)->edge(mesh, u)

function ⋆̂(::Edge, mesh::T, h::H) where {T}

end

⋆̂(mesh::T, h::H) where {T,H} = ⋆̂(simplex_trait(H), mesh, h)
==#

end # module


text:
module SegwayRobot3

    include("src3/segway.jl")
    #export Segway_mutable, Segway_unmutable

    include("src3/dynamics_type.jl")
    include("src3/dynamics_vec.jl")
    #export segway_rk4, segway_rk4!, dxdt_segway

    const params_segway = [-1.23, -1.20, -3.70, -7.51, -14.7, -10.2, -11.9, -3.43]
    const s_init  = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
    const u_range = [1.0, 1.0]

    include("src3/dynmodel.jl")

    cost_func = segway_cost # set standard Environment cost-function

    function rk4_func(xt::AbstractVector, ut::AbstractVector, dt::Number)

        #run the simulation
        xt1 = dynmodel(xt, ut, dt, params_segway)

        return xt1
    end

    export rk4_func, cost_func, s_init, u_range
end # module SegwayRobot3

max_stars_repo_name: NikonPic/REINFORCEMENT-JULIA

max_issues_repo_name: NikonPic/REINFORCEMENT-JULIA

max_forks_repo_name: NikonPic/REINFORCEMENT-JULIA