text:
"""
Generate coefficients related to PML boundaries, e.g., damping profiles.
This function outputs a, b, and k variables in eq. 25-26 from Komatitsch 2007 (Geophysics).

	TODO
* At the moment, K=1, and is dummy, unlike in the case of EM when it was useful. Removing K and K_inv might save us memory when using GPUs.
"""
function update_pml!(
    pml::NamedVector{T},
    mgrid, # 1D grid
    flags::Vector{Bool},
    δt::Float64,
    np::Int64,
    velavg::Float64,
    freqpeak::Float64, # dominant frequency in Hz
) where {T<:Data.Array}

    δx = step(mgrid)
    nx = length(mgrid)

    NPOWER = 2.e0
    K_MAX_PML = 1.e0 # from Gedney page 8.11
    ALPHA_MAX_PML = pi * freqpeak # from Festa and Vilotte

    # thickness of the PML layer in meters
    thickness_PML = np * δx

    "reflection coefficient (INRIA report section 6.1) http://hal.inria.fr/docs/00/07/32/19/PDF/RR-3471.pdf"
    Rcoef = 0.001e0

    #! check that NPOWER is okay
    #if(NPOWER < 1) stop "NPOWER must be greater than 1"

    # compute d0 from INRIA report section 6.1 http://hal.inria.fr/docs/00/07/32/19/PDF/RR-3471.pdf
    d0 = -(NPOWER + 1) * (velavg) * log(Rcoef) / (2.e0 * thickness_PML)

    k = zeros(nx)
    k_half = zero(k)
    fill!(k, 1)
    fill!(k_half, 1)
    d = zeros(nx)
    d_half = zero(d)
    alpha = zero(d)
    alpha_half = zero(d)
    a=zero(d)
    b=zero(d)
    a_half=zero(d)
    b_half=zero(d)

    "damping in the X direction"
    "origin of the PML layer (position of right edge minus thickness, in meters)"
    xoriginleft = thickness_PML
    xoriginright = (nx - 1) * δx - thickness_PML

    for ix = 1:nx
        # abscissa of current grid point along the damping profile
        xval = δx * real(ix - 1)

        #---------- left edge
        if (flags[1])

            # define damping profile at the grid points
            abscissa_in_PML = xoriginleft - xval
            if (abscissa_in_PML >= 0.0)
                abscissa_normalized = abscissa_in_PML / thickness_PML
                d[ix] = d0 * abscissa_normalized .^ NPOWER
                # this taken from Gedney page 8.2
                k[ix] = 1.e0 + (K_MAX_PML - 1.e0) * abscissa_normalized .^ NPOWER
                alpha[ix] =
                    ALPHA_MAX_PML * (1.e0 - abscissa_normalized) + 0.01e0 * ALPHA_MAX_PML
            end

            # define damping profile at half the grid points
            abscissa_in_PML = xoriginleft - (xval + δx / 2.e0)
            if (abscissa_in_PML >= 0.0)
                abscissa_normalized = abscissa_in_PML / thickness_PML
                d_half[ix] = d0 * abscissa_normalized .^ NPOWER
                # this taken from Gedney page 8.2
                k_half[ix] = 1.e0 + (K_MAX_PML - 1.e0) * abscissa_normalized .^ NPOWER
                alpha_half[ix] =
                    ALPHA_MAX_PML * (1.e0 - abscissa_normalized) + 0.01e0 * ALPHA_MAX_PML
            end
        end

        #---------- right edge
        if (flags[2])

            # define damping profile at the grid points
            abscissa_in_PML = xval - xoriginright
            if (abscissa_in_PML >= 0.0)
                abscissa_normalized = abscissa_in_PML / thickness_PML
                d[ix] = d0 * abscissa_normalized^NPOWER
                # this taken from Gedney page 8.2
                k[ix] = 1.e0 + (K_MAX_PML - 1.e0) * abscissa_normalized^NPOWER
                alpha[ix] =
                    ALPHA_MAX_PML * (1.e0 - abscissa_normalized) + 0.01e0 * ALPHA_MAX_PML
            end

            # define damping profile at half the grid points
            abscissa_in_PML = xval + δx / 2.e0 - xoriginright
            if (abscissa_in_PML >= 0.0)
                abscissa_normalized = abscissa_in_PML / thickness_PML
                d_half[ix] = d0 * abscissa_normalized^NPOWER
                # this taken from Gedney page 8.2
                k_half[ix] = 1.e0 + (K_MAX_PML - 1.e0) * abscissa_normalized^NPOWER
                alpha_half[ix] =
                    ALPHA_MAX_PML * (1.e0 - abscissa_normalized) + 0.01e0 * ALPHA_MAX_PML
            end

        end
        #
        # just in case, for -5 at the end
        (alpha[ix] < 0.0) ? alpha[ix] = 0.0 : nothing
        (alpha_half[ix] < 0.0) ? alpha_half[ix] = 0.0 : nothing

        # see equation 25 Komatitsch, 2007, get b and a (need k, alpha, and d before this )
        b[ix] = exp(-(d[ix] / k[ix] + alpha[ix]) * δt)
        b_half[ix] = exp(-(d_half[ix] / k_half[ix] + alpha_half[ix]) * δt)

        # this to avoid division by zero outside the PML
        (abs(d[ix]) > 1.e-6) ?
        a[ix] = d[ix] * (b[ix] - 1.e0) / (k[ix] * (d[ix] + k[ix] * alpha[ix])) : nothing
        (abs(d_half[ix]) > 1.e-6) ?
        a_half[ix] =
            d_half[ix] * (b_half[ix] - 1.e0) /
            (k_half[ix] * (d_half[ix] + k_half[ix] * alpha_half[ix])) : nothing

    end
    copyto!(pml[:kI], inv.(k))
    copyto!(pml[:k_halfI], inv.(k_half))
    copyto!(pml[:a],a)
    copyto!(pml[:b],b)
    copyto!(pml[:a_half],a_half)
    copyto!(pml[:b_half],b_half)

end

"""
Generate a NamedArray with PML coefficients for all the dimensions that are then stored in the FDTD structs.
"""
function get_pml(mgrid)
    dnames = dim_names(length(mgrid))
    pnames = [:a, :b, :kI, :a_half, :b_half, :k_halfI]
    return NamedArray(
        [
            NamedArray(
                Data.Array.([
                    zeros(ni),
                    zeros(ni),
                    ones(ni),
                    zeros(ni),
                    zeros(ni),
                    ones(ni),
                ]),
                pnames,
            ) for ni in length.(mgrid)
        ],
        dnames,
    )
end

"""
Just a loop over dims of pml
"""
function update_pml!(
    pml::NamedVector{T},
    mgrid,
    pml_edges::Vector{Symbol},
    args...,
) where {T<:NamedVector}
    for (i, dim) in enumerate(dim_names(length(mgrid)))
        update_pml!(
            pml[dim],
            mgrid[i],
            [
                any(pml_edges .== Symbol(string(dim), "min")),
                any(pml_edges .== Symbol(string(dim), "max")),
            ],
            args...,
        )
    end
end

function update_pml!(pac)
    update_pml!(
        pac.pml,
        pac.exmedium.mgrid,
        pac.pml_edges,
        pac.fc[:dt],
        npml-3,
        Statistics.mean(pac.exmedium.bounds[:vp]),
        pac.fc[:freqpeak],
    )
end

for dimnames in [zip([:1, :2, :3], dim_names(3)), zip([:1, :2], dim_names(2))]
    is = broadcast(x -> Symbol(string("i", x)), getindex.(collect(dimnames), 2))
    ist = Meta.parse(string("(", [string(s, ",") for s in is]..., ")"))
    N = Meta.parse(string(length(is)))

    for (idim, dim) in dimnames
        i = Symbol("i", string(dim))
        ismoff = replace(is, i => :($i + moff))
        isdoff = replace(is, i => :(doff + $i))
        for (fname, fnamenp, idoff) in zip(
            [Symbol("memory", string(dim), "!"), Symbol("memory1", string(dim), "!")],
            [Symbol("memorynp", string(dim), "!"), Symbol("memorynp1", string(dim), "!")],
            [:($i + doff), :($i + doff + 1)],
        )
            @eval @parallel_indices(
                $ist,
                function $fnamenp(memory::Data.Array{$N}, d, a, b, kI, moff, doff)
                    memory[$(ismoff...)] =
                        b[$idoff] * memory[$(ismoff...)] + a[$idoff] * d[$(isdoff...)]
                    d[$(isdoff...)] = d[$(isdoff...)] * kI[$idoff] + memory[$(ismoff...)]
                    return
                end
            )
            @eval function $fname(memory::Data.Array{$N}, d, a, b, kI)
                sm = collect(size(memory))
                setindex!(sm, npml, $idim)
                # first npml points
                @parallel map(x -> (:)(1, x), Tuple(sm)) $fnamenp(memory, d, a, b, kI, 0, 0)
                # last npml points independent of d
                @parallel map(x -> (:)(1, x), Tuple(sm)) $fnamenp(
                    memory,
                    d,
                    a,
                    b,
                    kI,
                    npml,
                    getindex(size(d), $idim) - npml,
                )
            end
        end
    end
end


text:
@testset "96.unique-binary-search-trees.jl" begin
    @test num_of_trees.(1:9) == [1, 2, 5, 14, 42, 132, 429, 1430, 4862]
end

ext: jl


text:
function collect_tags{T <: OSMElement}(osm_elements::Vector{T})
    tags = Set{UTF8String}()
    for n in osm_elements
        if isdefined(n, :tags)
            for k in keys(n.tags)
                push!(tags, k)
            end
        end
    end
    collect(tags)
end

_field(::Type{Node}) = :latlon
_field(::Type{Way}) = :nodes
_field(::Type{Relation}) = :members
_fieldvalue(node::Node) = node.latlon
_fieldvalue(way::Way) = way.nodes
_fieldvalue(relation::Relation) = relation.members
_fieldtype(::Type{Node}) = Tuple{Float64,Float64}
_fieldtype(::Type{Way}) = Vector{Int}
_fieldtype(::Type{Relation}) = Vector{Dict{UTF8String,UTF8String}}

function osm2dataframe{T <: OSMElement}(osm::Vector{T},
                                        colnames::Vector{UTF8String})
    nrows, ncols = length(osm), length(colnames)
    if nrows == 0
        return DataFrames.DataFrame()
    end
    df = DataFrames.DataFrame(vcat([Int,_fieldtype(T)],
                                   fill(UTF8String, ncols)),
                              vcat([:id, _field(T)],
                                   map(Symbol,colnames)),
                              nrows)
    for i in 1:nrows
        element = osm[i]
        df[i, 1] = element.id
        df[i, 2] = _fieldvalue(element)
        for j in 1:ncols
            if isdefined(element, :tags)
                df[i, j+2] = get(element.tags, colnames[j], DataFrames.NA)
            end
        end
    end
    df
end

osm2dataframe{T <: OSMElement}(osm::Vector{T}) =
    osm2dataframe(osm, collect_tags(osm))


text:
# This file is part of Kpax3. License is MIT.

function test_partition_rows_functions()
  for (α, θ) in ((0.4, -0.3), (0.4, 0.0), (0.4, 2.1), (0.0, 2.1), (-2.4, 3))
    ep = Kpax3.EwensPitman(α, θ)

    for i in 1:6
      po = TestPartition.Partition(i)

      pr = 0.0
      qr = 0.0

      for j in 1:po.B
        pr += Kpax3.dPriorRow(po.partition[:, j], ep)
        qr += exp(Kpax3.logdPriorRow(po.partition[:, j], ep))
      end

      @test isapprox(pr, 1.0, atol=ε)
      @test isapprox(qr, 1.0, atol=ε)

      pr = 0.0
      qr = 0.0

      for j in 1:(po.C - 1)
        pr += (po.index[j + 1] - po.index[j]) * Kpax3.dPriorRow(i, po.k[j], po.blocksize[:, j], ep)
        qr += (po.index[j + 1] - po.index[j]) * exp(Kpax3.logdPriorRow(i, po.k[j], po.blocksize[:, j], ep))
      end

      pr += Kpax3.dPriorRow(i, po.k[po.C], po.blocksize[:, po.C], ep)
      qr += exp(Kpax3.logdPriorRow(i, po.k[po.C], po.blocksize[:, po.C], ep))

      @test isapprox(pr, 1.0, atol=ε)
      @test isapprox(qr, 1.0, atol=ε)
    end
  end

  nothing
end

test_partition_rows_functions()


text:
@testset "Hybrid System tests" begin
    test_sys = PSB.build_system(PSB.PSITestSystems, "c_sys14"; add_forecasts = false)

    h_sys = HybridSystem(
        name = "Test H",
        available = true,
        status = true,
        bus = get_component(Bus, test_sys, "Bus 1"),
        active_power = 1.0,
        reactive_power = 1.0,
        thermal_unit = ThermalStandard(nothing),
        electric_load = PowerLoad(nothing),
        storage = GenericBattery(nothing),
        renewable_unit = RenewableDispatch(nothing),
        base_power = 100.0,
        operation_cost = TwoPartCost(nothing),
    )
    add_component!(test_sys, h_sys)

    initial_time = Dates.DateTime("2020-09-01")
    resolution = Dates.Hour(1)
    other_time = initial_time + resolution
    name = "test"
    horizon = 24
    data = Dict(initial_time => ones(horizon), other_time => 5.0 * ones(horizon))
    forecast = Deterministic(name, data, resolution)
    add_time_series!(test_sys, h_sys, forecast)
    ts = get_time_series(Deterministic, h_sys, "test")
    @test isa(ts, Deterministic)
end

@testset "Hybrid System from parsed files" begin
    sys = create_rts_system_with_hybrid_system(add_forecasts = true)
    hybrids = collect(get_components(HybridSystem, sys))
    @test length(hybrids) == 1
    h_sys = hybrids[1]

    electric_load = nothing
    thermal_unit = nothing
    subcomponents = collect(get_subcomponents(h_sys))
    @test length(subcomponents) == 4
    for subcomponent in subcomponents
        @test !PSY.is_attached(subcomponent, sys)
        @test IS.is_attached(subcomponent, sys.data.masked_components)
        if subcomponent isa PowerLoad
            electric_load = subcomponent
        elseif subcomponent isa ThermalStandard
            thermal_unit = subcomponent
        end
    end
    @test electric_load !== nothing
    @test thermal_unit !== nothing

    ts = collect(get_time_series_multiple(h_sys, type = TimeSeriesData))
    @test length(ts) == 0

    @test get_time_series(SingleTimeSeries, electric_load, "max_active_power") isa
          SingleTimeSeries
    @test get_time_series(AbstractDeterministic, electric_load, "max_active_power") isa
          DeterministicSingleTimeSeries

    @test !has_time_series(thermal_unit)
    @test has_time_series(electric_load)
    remove_time_series!(sys, AbstractDeterministic, electric_load, "max_active_power")
    remove_time_series!(sys, SingleTimeSeries, electric_load, "max_active_power")
    @test !has_time_series(electric_load)

    # Can't set the units when the HybridSystem is attached to system.
    @test_throws ArgumentError PSY.set_thermal_unit!(h_sys, thermal_unit)
end

@testset "Hybrid System from unattached subcomponents" begin
    sys = PSB.build_system(PSB.PSITestSystems, "test_RTS_GMLC_sys"; add_forecasts = false)
    thermal_unit = first(get_components(ThermalStandard, sys))
    bus = get_bus(thermal_unit)
    electric_load = first(get_components(PowerLoad, sys))
    storage = first(get_components(GenericBattery, sys))
    renewable_unit = first(get_components(RenewableDispatch, sys))

    for subcomponent in (thermal_unit, electric_load, storage, renewable_unit)
        remove_component!(sys, subcomponent)
    end

    name = "Test H"
    h_sys = HybridSystem(
        name = name,
        available = true,
        status = true,
        bus = bus,
        active_power = 1.0,
        reactive_power = 1.0,
        thermal_unit = thermal_unit,
        electric_load = electric_load,
        storage = storage,
        renewable_unit = renewable_unit,
        base_power = 100.0,
        operation_cost = TwoPartCost(nothing),
    )
    add_component!(sys, h_sys)

    for subcomponent in (thermal_unit, electric_load, storage, renewable_unit)
        @test !PSY.is_attached(subcomponent, sys)
        @test !has_time_series(subcomponent)
        @test IS.is_attached(subcomponent, sys.data.masked_components)
    end

    @test length(IS.get_masked_components(Component, sys.data)) == 4

    # Add time series for a subcomponent.
    initial_time = Dates.DateTime("2020-09-01")
    resolution = Dates.Hour(1)
    other_time = initial_time + resolution
    name = "test"
    horizon = 24
    data = Dict(initial_time => rand(horizon), other_time => 5.0 * ones(horizon))

    forecast = Deterministic(name, data, resolution)
    add_time_series!(sys, thermal_unit, forecast)
    @test get_time_series(Deterministic, thermal_unit, name) isa Deterministic
end


text:
abstract type Setup end

Base.@kwdef struct SphereSetup <: Setup
    name::String
    ncell::Int
    rmin::Float64
    rmax::Float64
    dmin::Float64
    dmax::Float64
    include_in::Bool
    in_ratio::Float64
    ecs_shape::String
    ecs_ratio::Float64
    refinement::Union{Float64,Nothing} = nothing
    D_in::SMatrix{3, 3, Float64}
    D_out::SMatrix{3, 3, Float64}
    D_ecs::SMatrix{3, 3, Float64}
    T₂_in::Float64
    T₂_out::Float64
    T₂_ecs::Float64
    ρ_in::Float64
    ρ_out::Float64
    ρ_ecs::Float64
    κ_in_out::Float64
    κ_out_ecs::Float64
    κ_in::Float64
    κ_out::Float64
    κ_ecs::Float64
end

Base.@kwdef struct CylinderSetup <: Setup
    name::String
    ncell::Int
    rmin::Float64
    rmax::Float64
    dmin::Float64
    dmax::Float64
    height::Float64
    bend::Float64
    twist::Float64
    include_in::Bool
    in_ratio::Float64
    ecs_shape::String
    ecs_ratio::Float64
    refinement::Union{Float64,Nothing} = nothing
    D_in::SMatrix{3, 3, Float64}
    D_out::SMatrix{3, 3, Float64}
    D_ecs::SMatrix{3, 3, Float64}
    T₂_in::Float64
    T₂_out::Float64
    T₂_ecs::Float64
    ρ_in::Float64
    ρ_out::Float64
    ρ_ecs::Float64
    κ_in_out::Float64
    κ_out_ecs::Float64
    κ_in::Float64
    κ_out::Float64
    κ_ecs::Float64
end

Base.@kwdef struct NeuronSetup <: Setup
    name::String
    ncell::Int = 1
    include_in::Bool = false
    in_ratio::Float64 = 0.0
    ecs_shape::String
    ecs_ratio::Float64
    refinement::Union{Float64,Nothing} = nothing
    D_in::SMatrix{3, 3, Float64}
    D_out::SMatrix{3, 3, Float64}
    D_ecs::SMatrix{3, 3, Float64}
    T₂_in::Float64
    T₂_out::Float64
    T₂_ecs::Float64
    ρ_in::Float64
    ρ_out::Float64
    ρ_ecs::Float64
    κ_in_out::Float64
    κ_out_ecs::Float64
    κ_in::Float64
    κ_out::Float64
    κ_ecs::Float64
end


text:
using Catalan
using Base.Test

@test ([5,4,2,2]\[2,1]) == ([5, 4, 2, 2],[2, 1])
@test isrimhook([4,3,2], [2,2,2])
@test !isrimhook([4,3,2], [2,2,1])
@test !isrimhook([4,3,2], [1,1])
λ = [5,4,2,1]
@test partitionsequence(λ) == [1, 0, 1, 0, 1, 1, 0, 1, 0]
@test character(λ, [4,3,2,2,1]) == 0
@test character([1], [1]) == 1


text:
using FastConv
using Test

A = zeros(5,5)
A[3,3] = 1
k = rand(3,3)
Ac = convn(A,k)
@test Ac[3:5,3:5] == k && all(a->a==0, Ac[:,[1,2,6,7]]) && all(a->a==0, Ac[[1,2,6,7], :])

A = zeros(5,6,7)
A[3,4,5] = 1
k = rand(3,3,3)
Ac = convn(A,k)
@test Ac[3:5,4:6,5:7] == k &&
    all(a->a==0, Ac[[1,2,6,7], :, :]) &&
    all(a->a==0, Ac[:, [1,2,3,7,8], :]) &&
    all(a->a==0, Ac[:, :, [1,2,3,4,8,9]])

A = randn(5,5)
k = ones(1,1)
Ac = convn(A,k)
@test Ac == A

A = randn(5,5)
k = zeros(3,3)
Ac = convn(A,k)
@test Ac == zeros(7,7)

A = randn(5,5)
k1 = rand(3,3)
k2 = rand(3,3)
A1_2 = convn(convn(A,k1), k2)
A_12 = fastconv(A, fastconv(k1, k2))
@test A1_2 \approx A_12


text:
########################################################################
#                 REINDEXING & OP BINARIZATION                         #
########################################################################

function reindex(op::Call, st::Dict)
    new_args = [get(st, x, x) for x in op.args]
    new_id = get(st, op.id, op.id)
    return copy_with(op, args=new_args, id=new_id)
end

reindex(op::Assign, st::Dict) = copy_with(op, id=get(st, op.id, op.id),
                                          src_id=get(st, op.src_id, op.src_id))
reindex(op::AbstractOp, st::Dict) = copy_with(op, id=get(st, op.id, op.id))

function reindex_fields!(tape::Tape, st::Dict)
    tape.resultid = get(st, tape.resultid, tape.resultid)
    tape.derivs = Dict(get(st, k, k) => get(st, v, v) for (k, v) in tape.derivs)
    # possibly we also need to reindex .derivs
end

function reindex(tape::Tape, st::Dict)
    new_tape = copy_with(tape, ops=AbstractOp[])
    for op in tape
        push!(new_tape, reindex(op, st))
    end
    return new_tape
end

"""
Transform tape so that all calls to *, /, + and - with multiple arguments were split into
binary form, e.g. transform:

*(a, b, c, d)

to

w1 = a * b
w2 = w1 * c
w3 = w2 * d
"""
function binarize_ops(tape::Tape)
    st = Dict()
    new_tape = copy_with(tape; ops=AbstractOp[])
    for (id, op) in enumerate(tape)
        # id == 11 && break
        # println(id)
        new_id = -1 # to be modified in both branches
        _op = reindex(op, st)   # we need both - original op and reindexed _op
        if op isa Call && op.fn in (+, -, *, /) && length(op.args) > 2
            # record first pair of arguments
            arg1_val, arg2_val = tape[op.args[1]].val, tape[op.args[2]].val
            new_val = op.fn(arg1_val, arg2_val)
            new_id = record!(new_tape, Call, new_val, op.fn, _op.args[1:2])
            # record the rest
            for i=3:length(op.args)
                # calculate new value based on current value and value of original ith argument
                new_val = op.fn(new_val, tape[op.args[i]].val)
                # record new binary op with the previous new_id and new argument id
                new_id = record!(new_tape, Call, new_val, op.fn, [new_id, _op.args[i]])
            end
            # st[id] = new_id
        else
            new_id = length(new_tape) + 1
            push!(new_tape, copy_with(_op; id=new_id))
        end
        if id != new_id
            st[id] = new_id
        end
    end
    reindex_fields!(new_tape, st)
    return new_tape
end

########################################################################
#                       TAPE TRANSFORMATIONS                           #
########################################################################

function find_deps(tape::Tape, ids::Vector{Int}; result=Set{Int}())
    ids = filter(id -> id != -1, ids)
    for id in ids
        # take only calls and skip ops that have already been processed
        if tape[id] isa Call && !in(id, result)
            args = tape[id].args
            find_deps(tape, args; result=result)
            push!(result, args...)
        end
    end
    return result
end

find_deps(tape::Tape, id::Int; result=Set{Int}()) = find_deps(tape, [id]; result=result)

function remove_unused(tape::Tape; keep=[tape.resultid])
    deps = find_deps(tape, keep); push!(deps, keep...)
    st = Dict()
    new_tape = copy_with(tape; ops=AbstractOp[])
    for (id, op) in enumerate(tape)
        if id in deps
            new_id = length(new_tape) + 1
            if id != new_id
                st[id] = new_id
            end
            new_op = reindex(op, st)
            push!(new_tape, new_op)
        end
    end
    reindex_fields!(new_tape, st)
    return new_tape
end

function eliminate_common(tape::Tape)
    new_tape = copy_with(tape, ops=AbstractOp[])
    existing = Dict()
    st = Dict()
    for op in tape
        op = reindex(op, st)
        if isa(op, Call)
            key = (op.fn, op.args...)
            if haskey(existing, key)
                # such expression already exists
                # replace curren one with assignment
                src = existing[key]
                st[op.id] = src
                push!(new_tape, Assign(op.id, src, tape[src].val))
            else
                # new expression - put to `existing` and record to new_tape
                existing[key] = op.id
                push!(new_tape, op)
            end
        else
            push!(new_tape, op)
        end
    end
    return squash_assigned(new_tape)
end

########################################################################
#                       PRE- & POST-PROCESSING                         #
########################################################################

"""
Apply a number of transformations to a tape after tracing and before calculating derivatives
"""
function preprocess(tape::Tape)
    tape = binarize_ops(tape)
    tape = remove_unused(tape)
    return tape
end

"""
Apply a number of transformations after calculating derivatives
"""
function postprocess(tape::Tape)
    tape = eliminate_common(tape)
    return tape
end

########################################################################
#                          HASH & EQUALITY                             #
########################################################################

# note: hash and == don't use operation values and .compiled field
# so different runs of a tape doesn't affect its identity

Base.hash(op::Input, h::UInt) = hash(op.id, h)
Base.hash(op::Constant, h::UInt) = hash(op.id, hash(op.val, h))
Base.hash(op::Assign, h::UInt) = hash(op.id, hash(op.src_id, h))
Base.hash(op::Call, h::UInt) = hash(op.id, hash(op.fn, hash(op.args, h)))

function Base.hash(tape::Tape, h::UInt)
    h = hash(tape.resultid)
    h = hash(tape.derivs, h)
    h = hash(tape.fieldpaths, h)
    h = hash(tape.device, h)
    for op in tape
        h = hash(op, h)
    end
    return h
end

Base.:(==)(op1::Input, op2::Input) = (op1.id == op2.id)
Base.:(==)(op1::Constant, op2::Constant) = (op1.id == op2.id && op1.val == op2.val)
Base.:(==)(op1::Assign, op2::Assign) = (op1.id == op2.id && op1.src_id == op2.src_id)
Base.:(==)(op1::Call, op2::Call) = (op1.id == op2.id && op1.fn == op2.fn && op1.args == op2.args)

function Base.:(==)(tape1::Tape, tape2::Tape)
    return tape1.resultid == tape2.resultid &&
        tape1.derivs == tape2.derivs &&
        tape1.fieldpaths == tape2.fieldpaths &&
        tape1.device == tape2.device &&
        tape1.ops == tape2.ops
end

########################################################################
#                       __NEW__() & __GETPROPERTY__()                  #
########################################################################

"""
Find the variable that was used to construct the specified field of the constructor op
"""
function field_var_from_ctor_op(tape::Tape, ctor_op::Call, fldname::Symbol)
    ctor_typ = ctor_op.fn
    for ctor_def in CONSTRUCTORS
        if ctor_typ == ctor_def[1]
            fld_idx_in_ctor_def = findfirst(isequal(fldname), ctor_def[2:end])
            if fld_idx_in_ctor_def != nothing
                # tape ID of corresponding argument to the constructor
                arg_id = ctor_op.args[fld_idx_in_ctor_def]
                return tape[arg_id]
            end
        end
    end
    return nothing
end

# TODO: make unified API for 2 versions of this function
# 1st - for __new__ constuctor
# 2nd - for primitive constructor created using @ctor macro
function field_var_from_ctor_op(tape::Tape, ctor::Call, getprop_op::Call)
    @assert ctor.fn == __new__
    @assert getprop_op.fn == Base.getproperty
    T = typeof(ctor.val)
    flds = fieldnames(T)
    fld = tape[getprop_op.args[2]].val
    fld_arg_idx = findfirst(x -> x == fld, flds) + 1
    return tape[ctor.args[fld_arg_idx]]
end

"""
Given a tape and getproperty() operation, try to find a variable
that was used to create that field
"""
function find_field_source_var(tape::Tape, getprop_op::Call)
    parent = tape[getprop_op.args[1]]
    if parent isa Call && parent.fn == __new__
        # found constructor for this field, return variable that made up getprop_op's field
        return field_var_from_ctor_op(tape, parent, getprop_op)
    elseif parent isa Call && parent.fn == Base.getproperty
        # nested getproperty, find constructor for the current struct recursively
        ctor = find_field_source_var(tape, parent)
        if ctor != nothing
            return field_var_from_ctor_op(tape, ctor, getprop_op)
        else
            return nothing
        end
    else
        # can't find source field - give up and return nothing
        return nothing
    end
end

"""
Given a tape and getfield() operation, try to find a variable
that was used to create that field.
NOTE: This works only with getfield() on tuples
"""
function find_tuple_field_source_var(tape::Tape, getf_op::Call)
    getf_base_op = tape[getf_op.args[1]]
    if getf_base_op.fn == Base.indexed_iterate
        tuple_op = tape[getf_base_op.args[1]]
        # tuple_idx = tape[ind_it_op.args[2]].val
        tuple_idx = tape[getf_base_op.args[2]].val
    elseif getf_base_op.fn == __tuple__
        tuple_op = getf_base_op
        tuple_idx = tape[getf_op.args[2]].val
    elseif getf_base_op.fn == namedtuple
        tuple_op = tape[getf_base_op.args[2]]
        tuple_fld = tape[getf_op.args[2]].val
        flds = fieldnames(typeof(getf_base_op.val))
        tuple_idx = findfirst(x -> x == tuple_fld, flds)
        # should make recursive?
    else
        throw(AssertionError("Unexpected base op for __getfield__: $(getf_base_op.fn)"))
    end
    # @assert tuple_op isa Call && tuple_op.fn == __tuple__
    src_var = tape[tuple_op.args[tuple_idx]]
    return src_var
end


text:
using Test, PolynomialBases

if !haskey(ENV, "JULIA_PKGEVAL") # sympy is not install on https://github.com/JuliaComputing/NewPkgEval.jl
  import SymPy
  x = SymPy.symbols("x")

  @test 0 == SymPy.simplify( hermite(x, 0) - 1 )
  @test 0 == SymPy.simplify( hermite(x, 1) - ( 2x ) )
  @test 0 == SymPy.simplify( hermite(x, 2) - ( 4x^2 - 2 ) )
  @test 0 == SymPy.simplify( hermite(x, 3) - ( 8x^3 - 12x ) )
  @test 0 == SymPy.simplify( hermite(x, 4) - ( 16x^4 - 48x^2 + 12 ) )

  @test_skip hermite(x, 4)
end
@inferred hermite(10., 4)


text:
using Router, DemosController

route("/") do
  Router.serve_static_file("/welcome.html")
end

route("/demos/outputid/:id::Int", DemosController.demo)

route("/demos/outputidinline/:id::Int") do
  "From inline routes $(@params(:id))"
end

route("/demos", DemosController.index)
route("/dump", DemosController.dumpform, method = POST)


text:
# Create example point data
x = 8 .* rand(100) .- 4
y = 8 .* rand(100) .- 4
z = sin.(x) + cos.(y)
# Tricontour plot
tricont(x, y, z)


text:
# Computes derivatives over the timestep.
using PyPlot
using Limbdark
import Limbdark: Transit_Struct

#include("../../../src/integrate_lightcurve.jl")

# Test it out:

include("../../../test/loglinspace.jl")

t1 = -1.5; t2 = 1.5; nt = 10000; dt = 0.3
t = zeros(nt)
t .= linearspace(t1,t2,nt)
#t = [t[1235]] ; nt =1
# The following compares different tolerances and maxdepths:
r = 0.1; b0 = 0.5; u_n = [0.3,0.3]; nu = length(u_n)
favg0 = zeros(nt,5+nu)
favg1 = zeros(nt,5+nu)
neval1= zeros(Int64,nt)
depthmax1= zeros(Int64,nt)
favg2 = zeros(nt,5+nu)
neval2= zeros(Int64,nt)
depthmax2= zeros(Int64,nt)
favg3 = zeros(nt,5+nu)
favg4 = zeros(nt,5+nu)
favg5 = zeros(nt,5+nu)
favg6 = zeros(nt,5+nu)
favg7 = zeros(nt,5+nu)
neval3= zeros(Int64,nt)
neval4= zeros(Int64,nt)
neval5= zeros(Int64,nt)
neval6= zeros(Int64,nt)
neval7= zeros(Int64,nt)
depthmax3= zeros(Int64,nt)
depthmax4= zeros(Int64,nt)
depthmax5= zeros(Int64,nt)
depthmax6= zeros(Int64,nt)
depthmax7= zeros(Int64,nt)
trans = transit_init(r,b0,u_n,true)
param = [0.0,1.0,b0]   # [t_0,v,b_0]

# Now, carry out speed test:

function compute_lightcurve!(trans::Transit_Struct{T},param::Array{T,1},t::Array{T,1},favg0::Array{T,2},nt::Int64) where {T <: Real}
@inbounds for i=1:nt
  trans.b = sqrt(param[3]^2+(param[2]*(t[i]-param[1]))^2)
  favg0[i,1] =  transit_poly!(trans)-1
  favg0[i,2] =  trans.dfdrb[1]
  favg0[i,3] =  trans.dfdrb[2]/trans.b*param[2]^2*(param[1]-t[i])
  favg0[i,4] =  trans.dfdrb[2]*param[2]/trans.b*(t[i]-param[1])^2
  favg0[i,5] =  trans.dfdrb[2]*param[3]/trans.b
  favg0[i,6:5+trans.n] = trans.dfdu
end
return
end

compute_lightcurve!(trans,param,t,favg0,nt)
timing = zeros(8)
#@time compute_lightcurve!(trans,param,t,favg0,nt)
tstart=time_ns(); compute_lightcurve!(trans,param,t,favg0,nt); timing[1] = (time_ns()-tstart)*1e-9

integrate_lightcurve!(trans,param,t,dt,favg1,nt,1e-4,3,neval1,depthmax1)
#@time integrate_lightcurve!(trans,param,t,dt,favg1,nt,1e-4,8,neval1,depthmax1)
#tic(); integrate_lightcurve!(trans,param,t,dt,favg1,nt,1e-4,8,neval1,depthmax1); timing[2] = toc()
tstart=time_ns(); integrate_lightcurve!(trans,param,t,dt,favg1,nt,1e-4,8,neval1,depthmax1); timing[2] = (time_ns()-tstart)*1e-9
#@time integrate_lightcurve!(trans,param,t,dt,favg2,nt,1e-6,128,neval2,depthmax2)
tstart=time_ns(); integrate_lightcurve!(trans,param,t,dt,favg2,nt,1e-6,128,neval2,depthmax2); timing[3] =  (time_ns()-tstart)*1e-9
#@time integrate_lightcurve!(trans,param,t,dt,favg3,nt,1e-8,128,neval3,depthmax3)
tstart=time_ns(); integrate_lightcurve!(trans,param,t,dt,favg3,nt,1e-8,128,neval3,depthmax3); timing[4] =  (time_ns()-tstart)*1e-9
#@time integrate_lightcurve!(trans,param,t,dt,favg4,nt,1e-10,128,neval4,depthmax4)
tstart=time_ns(); integrate_lightcurve!(trans,param,t,dt,favg4,nt,1e-10,128,neval4,depthmax4); timing[5] =  (time_ns()-tstart)*1e-9
#@time integrate_lightcurve!(trans,param,t,dt,favg5,nt,1e-12,128,neval5,depthmax5)
tstart=time_ns(); integrate_lightcurve!(trans,param,t,dt,favg5,nt,1e-12,128,neval5,depthmax5); timing[6] =  (time_ns()-tstart)*1e-9
#@time integrate_lightcurve!(trans,param,t,dt,favg6,nt,1e-14,128,neval6,depthmax6)
tstart=time_ns(); integrate_lightcurve!(trans,param,t,dt,favg6,nt,1e-14,128,neval6,depthmax6); timing[7] =  (time_ns()-tstart)*1e-9
#@time integrate_lightcurve!(trans,param,t,dt,favg7,nt,1e-16,128,neval7,depthmax7)
tstart=time_ns(); integrate_lightcurve!(trans,param,t,dt,favg7,nt,1e-16,128,neval7,depthmax7); timing[8] =  (time_ns()-tstart)*1e-9

# Now, make some precision plots:

absolute_precision = zeros(7,6)
fractional_precision = zeros(7,6)
depthmax = [maximum(depthmax1),maximum(depthmax2),maximum(depthmax3),maximum(depthmax4),maximum(depthmax5)]
neval = [sum(neval1),sum(neval2),sum(neval3),sum(neval4),sum(neval5),sum(neval6)]
for j=1:7
  absolute_precision[j,1] = maximum(abs.(favg1[:,j]-favg7[:,j]))
  fractional_precision[j,1] = maximum(abs.(favg1[:,j]-favg7[:,j]))/maximum(abs.(favg7[:,j]))
end
for j=1:7
  absolute_precision[j,2] = maximum(abs.(favg2[:,j]-favg7[:,j]))
  fractional_precision[j,2] = maximum(abs.(favg2[:,j]-favg7[:,j]))/maximum(abs.(favg7[:,j]))
end
for j=1:7
  absolute_precision[j,3] = maximum(abs.(favg3[:,j]-favg7[:,j]))
  fractional_precision[j,3] = maximum(abs.(favg3[:,j]-favg7[:,j]))/maximum(abs.(favg7[:,j]))
end
for j=1:7
  absolute_precision[j,4] = maximum(abs.(favg4[:,j]-favg7[:,j]))
  fractional_precision[j,4] = maximum(abs.(favg4[:,j]-favg7[:,j]))/maximum(abs.(favg7[:,j]))
end
for j=1:7
  absolute_precision[j,5] = maximum(abs.(favg5[:,j]-favg7[:,j]))
  fractional_precision[j,5] = maximum(abs.(favg5[:,j]-favg7[:,j]))/maximum(abs.(favg7[:,j]))
end
for j=1:7
  absolute_precision[j,6] = maximum(abs.(favg6[:,j]-favg7[:,j]))
  fractional_precision[j,6] = maximum(abs.(favg6[:,j]-favg7[:,j]))/maximum(abs.(favg7[:,j]))
end

tol = [1e-4,1e-6,1e-8,1e-10,1e-12,1e-14]

clf()
labels = [L"$F-1$",L"$r$",L"$t_0$",L"$v$",L"$b_0$",L"$u_1$",L"$u_2$"]
for j=1:7
#  loglog(tol,absolute_precision[j,:],label=labels[j])
  loglog(neval/nt,absolute_precision[j,:],"o-",label=labels[j])
  #loglog(neval/nt,absolute_precision[j,:],"o",color="k")
#  loglog(tol,fractional_precision[j,:],label=labels[j],linestyle=":")
end
#plot(tol,tol,label="tolerance",linestyle=":")
plot(neval/nt,10*tol,label=L"$10\epsilon_{tol}$",linestyle=":")
legend(loc="lower left")
xlabel("Average number of evaluations per exposure")
ylabel("Absolute precision")
xlim(5, 1200)

twiny()
loglog(timing[2:7]/timing[1],10*tol,label=L"$10\epsilon_{tol}$",linestyle=:dashed)
legend(loc="upper right")
xlim(5, 1200)
xlabel("Time ratio per exposure")

savefig("compare_precision.pdf", bbox_inches="tight")

max_issues_repo_name: rodluger/limbdark


text:
"""
Functions obtanied in:

- https://www.sfu.ca/~ssurjano/optimization.html
- Molga, M., & Smutnicki, C. (2005). Test functions for optimization needs. Test functions for optimization needs, 101, 48.

"""

# N-Dimensional

# ACKLEY FUNCTION
"""
# Description
*Dimensions*: d

The **Ackley function** is widely used for testing optimization algorithms.
In its two-dimensional form, as shown in the plot above, it is characterized
by a nearly flat outer region, and a large hole at the centre.
The function poses a risk for optimization algorithms, particularly hillclimbing algorithms,
to be trapped in one of its many local minima.

Recommended variable values are: **a = 20, b = 0.2 and c = 2π**.

# Input Domain:

The function is usually evaluated on the hypercube **xi ∈ [-32.768, 32.768],
for all i = 1, …, d,** although it may also be restricted to a smaller domain.

# Global Minimum:

**f(X\\*) = 0**, at **X\\* = (0,...,0)**

# Link:

*https://www.sfu.ca/~ssurjano/ackley.html*
"""
function Ackley(X...)
    d = length(X)
    return -20*exp(-0.2*sqrt((1/d)*sum(X.^2))) - exp((1/d)*sum(cos.(2 .*pi.*X))) + 20 + exp(1)
end

# SCHWEFEL’S FUNCTION
"""
# Description
*Dimensions*: d

**Schwefel’s function** is deceptive in that the global minimum is geometrically distant,
over the parameter space, from the next best local minima. Therefore, the search algorithms
are potentially prone to convergence in the wrong direction.

# Input Domain:

Test area is usually restricted to hyphercube **xi ∈ [−500, 500], for all i = 1, . . . , d**.

# Global Minimum:

**f(X\\*) = -418.9829d**, at **X\\* = (420.9687,...,420.9687)**
"""
Schwefel(X...) = sum(@. -X*sin(sqrt(abs(X))))

# GRIEWANK FUNCTION
"""
# Description:
*Dimensions*: d

The **Griewank function** has many widespread local minima, which are regularly distributed. The complexity is shown in the zoomed-in plots.

# Input Domain:
The function is usually evaluated on the hypercube **xi ∈ [-600, 600], for all i = 1, …, d**.

# Global Minimum:
**f(X\\*) = 0**, at **X\\* = (0,...,0)**
"""
function Griewank(X...)
    i = 1:length(X)
    A = sum(X.^2 ./4000)
    B = prod(cos.(X./sqrt.(i)))
    return A - B + 1
end

# ROSENBROCK FUNCTION
"""
*Description*:
# Dimensions: d

The **Rosenbrock function**, also referred to as the Valley or Banana function, is a popular test problem for gradient-based optimization algorithms. It is shown in the plot above in its two-dimensional form.

The function is unimodal, and the global minimum lies in a narrow, parabolic valley.
However, even though this valley is easy to find, convergence to the minimum is difficult (Picheny et al., 2012).

# Input Domain:
The function is usually evaluated on the hypercube **xi ∈ [-5, 10], for all i = 1, …, d**, although it may be restricted
to the hypercube **xi ∈ [-2.048, 2.048], for all i = 1, …, d**.

# Global Minimum:

**f(X\\*) = 0**, at **X\\* = (1,...,1)**

# Link:

https://www.sfu.ca/~ssurjano/rosen.html
"""
function Rosenbrock(X...)
    x1 = X[1:length(X)-1]
    x2 = X[2:length(X)]
    return sum(@. 100*(x2 - x1^2)^2 + (x1 - 1)^2)
end

# MICHALEWICZ FUNCTION
"""
# Description:
*Dimensions*: d

The **Michalewicz function** has d! local minima, and it is multimodal.
The parameter m defines the steepness of they valleys and ridges; a larger m leads to
a more difficult search. The recommended value of m is m = 10.

# Input Domain:
The function is usually evaluated on the hypercube **xi ∈ [0, π], for all i = 1, …, d**.

Global Minima:

at **d = 2**: **f(X\\*) = -1.8013**, at **X\\* = (2.20, 1.57)**
at **d = 5**: **f(X\\*) = -4.687658**
at **d = 10**: **f(X\\*) = -9.66015**

# Link:

*https://www.sfu.ca/~ssurjano/michal.html*
"""
function Michalewicz(X...)
    m = 10
    i = 1:length(X)
    return -sum(@. sin(X)*(sin(i*X^2 / pi)^(2*m)))
end

# 2-Dimensional

# BUKIN FUNCTION N. 6
"""
# Description:
*Dimensions*: 2

The **sixth Bukin function** has many local minima, all of which lie in a ridge.

# Input Domain:

The function is usually evaluated on the rectangle **x1 ∈ [-15, -5], x2 ∈ [-3, 3]**.

# Global Minimum:

**f(X\\*) = 0**, at **X\\* = (-10,1)**

# Link:

*https://www.sfu.ca/~ssurjano/bukin6.html*
"""
Bukin6(x,y) = 100*sqrt(abs(y - 0.01*x^2)) + 0.01*abs(x + 10)

# CROSS IN TRAY FUNCTION
"""
# Description:
*Dimensions*: 2

The **Cross-in-Tray function** has multiple global minima. It is shown here with a smaller domain in the second plot, so that its characteristic "cross" will be visible.

# Input Domain:

The function is usually evaluated on the square **xi ∈ [-10, 10], for all i = 1, 2**.

# Global Minima:

**f(X\\*) = -2.06261**, at **X\\* = (1.3491, -1.3491), (-1.3491, 1.3491), (-1.3491, -1.3491)** and **(1.3491, 1.3491)**

# Link:

*https://www.sfu.ca/~ssurjano/crossit.html*
"""
function Cross_in_tray(x,y)
    s = sqrt(x^2 + y^2)/pi
    return -0.0001*(abs(sin(x)*sin(y)*exp(abs(100-s)))+1)^0.1
end

# DROP WAVE FUNCTION
"""
# Description:
*Dimensions*: 2

The **Drop-Wave function** is multimodal and highly complex. The second plot above shows the function on a smaller input domain, to illustrate its characteristic features.

# Input Domain:

The function is usually evaluated on the square **xi ∈ [-5.12, 5.12], for all i = 1, 2**.

# Global Minimum:

**f(X\\*) = -1**, at **X\\* = (0,0)**

# Link:

*https://www.sfu.ca/~ssurjano/drop.html*
"""
function Drop_wave(X...)
    s = sum(X.^2)
    return -(1 + cos(12*sqrt(s)))/(0.5*(s) + 2)
end

# EGGHOLDER FUNCTION
"""
# Description:
*Dimensions*: 2

The **Eggholder function** is a difficult function to optimize, because of the large number of local minima.

# Input Domain:
The function is usually evaluated on the square **xi ∈ [-512, 512], for all i = 1, 2**.

# Global Minimum:

**f(X\\*) = -959.6407**, at **X\\* = (512,404.2319)**

# Link:

*https://www.sfu.ca/~ssurjano/egg.html*
"""
function Eggholder(x,y)
    A = (y + 47)*sin(sqrt(abs(y + x/2 + 47)))
    B = x*sin(sqrt(abs(x - (y + 47))))
    return -(A + B)
end


text:
include(joinpath(@__DIR__, "../../src/FLS.jl"))

####################
# REFLEXIVE GRAPHS #
####################

reflgraphschema = @acset LabeledGraph begin
    V = 2
    E = 4
    vlabel = [:V, :E]
    elabel = [:id, :refl, :src, :tgt]
    src = [1,1,2,2]
    tgt = [1,2,1,1]
end

reflgraph = FLS(:reflgraph, reflgraphschema, [], [
    (:srcid, [:refl, :src], []),
    (:tgtid, [:refl, :tgt], [])])


text:
"""
A `Robot` represents a mechanism, its state, and other relevant objects.

# Relevant fields in `Robot`
- The arrays `q_lo`, `q_hi`, `v_lo`, `v_hi`, `τ_lo`, `τ_hi` store the lower and upper bounds of the joint positions, velocities, and torques, respectively.
- The number of generalized coordinates, generalized velocities, and actuated joints are stored in `n_q`, `n_v`, `n_τ`, respectively.
- The end-effector frame is stored in `frame_ee`.
"""
struct Robot
    urdfpath::String
    mechanism::Mechanism
    state::MechanismState
    statecache::StateCache
    dynamicsresultcache::DynamicsResultCache
    mvis::MechanismVisualizer

    q_lo::Array{Float64}
    q_hi::Array{Float64}
    v_lo::Array{Float64}
    v_hi::Array{Float64}
    τ_lo::Array{Float64}
    τ_hi::Array{Float64}

    n_q::Int64  # Number of generalized coordinates
    n_v::Int64  # Number of generalized velocities
    n_τ::Int64  # Number of actuated joints

    frame_ee::CartesianFrame3D  # End-effector frame

    @doc """
        Robot(urdfpath, mechanism, frame_ee, mvis)

    Create a new `Robot`.
    """
    function Robot(urdfpath::String,
                   mechanism::Mechanism,
                   frame_ee::CartesianFrame3D,
                   mvis::MechanismVisualizer)
        state = MechanismState(mechanism)
        statecache = StateCache(mechanism)
        dynamicsresultcache = DynamicsResultCache(mechanism)

        q_lo = [lim for joint in joints(mechanism) for lim in map(x -> x.lower, position_bounds(joint))]
        q_hi = [lim for joint in joints(mechanism) for lim in map(x -> x.upper, position_bounds(joint))]
        v_lo = [lim for joint in joints(mechanism) for lim in map(x -> x.lower, velocity_bounds(joint))]
        v_hi = [lim for joint in joints(mechanism) for lim in map(x -> x.upper, velocity_bounds(joint))]
        τ_lo = [lim for joint in joints(mechanism) for lim in map(x -> x.lower,   effort_bounds(joint))]
        τ_hi = [lim for joint in joints(mechanism) for lim in map(x -> x.upper,   effort_bounds(joint))]

        new(
            urdfpath,
            mechanism,
            state,
            statecache,
            dynamicsresultcache,
            mvis,
            q_lo,
            q_hi,
            v_lo,
            v_hi,
            τ_lo,
            τ_hi,
            num_positions(mechanism),
            num_velocities(mechanism),
            num_velocities(mechanism),
            frame_ee
        )
    end
end

"""
    create_robot_kuka_iiwa_14(vis)

Create a new [KUKA LBR iiwa 14](https://www.kuka.com/en-gb/products/robotics-systems/industrial-robots/lbr-iiwa) robot.
"""
function create_robot_kuka_iiwa_14(vis::Visualizer)
    package_path = joinpath(@__DIR__, "..", "iiwa_stack")
    urdfpath = joinpath(@__DIR__, "..", "robots", "iiwa14.urdf")

    mechanism = parse_urdf(urdfpath, remove_fixed_tree_joints=false)
    frame_ee = default_frame(findbody(mechanism, "iiwa_link_pen_tip"))
    remove_fixed_tree_joints!(mechanism)

    urdfvisuals = URDFVisuals(urdfpath, package_path=[package_path])
    mvis = MechanismVisualizer(mechanism, urdfvisuals, vis["robot"])
    # setelement!(mvis, frame_ee)  # Visualize a triad at the end-effector

    Robot(urdfpath, mechanism, frame_ee, mvis)
end

"""
    create_robot_kinova_gen3_lite(vis)

Create a new [Kinova Gen3 lite](https://www.kinovarobotics.com/en/products/gen3-lite-robot) robot.
"""
function create_robot_kinova_gen3_lite(vis::Visualizer)
    package_path = joinpath(@__DIR__, "..", "ros_kortex")
    urdfpath = joinpath(@__DIR__, "..", "robots", "gen3_lite_gen3_lite_2f.urdf")

    mechanism = parse_urdf(urdfpath, remove_fixed_tree_joints=false)
    frame_ee = default_frame(findbody(mechanism, "tool_frame"))
    remove_fixed_tree_joints!(mechanism)

    urdfvisuals = URDFVisuals(urdfpath, package_path=[package_path])
    mvis = MechanismVisualizer(mechanism, urdfvisuals, vis["robot"])
    # setelement!(mvis, frame_ee)  # Visualize a triad at the end-effector

    Robot(urdfpath, mechanism, frame_ee, mvis)
end

export
    Robot,
    create_robot_kuka_iiwa_14,
    create_robot_kinova_gen3_lite


text:
module TestConversions
    using Base.Test
    using DataFrames
    using DataStructures: OrderedDict, SortedDict

    df = DataFrame()
    df[:A] = 1:5
    df[:B] = [:A, :B, :C, :D, :E]
    @test isa(convert(Array, df), Matrix{Any})
    @test convert(Array, df) == convert(Array, convert(NullableArray, df))
    @test isa(convert(Array{Any}, df), Matrix{Any})

    df = DataFrame()
    df[:A] = 1:5
    df[:B] = 1.0:5.0
    # Fails on Julia 0.4 since promote_type(Nullable{Int}, Nullable{Float64}) gives Nullable{T}
    if VERSION >= v"0.5.0-dev"
        @test isa(convert(Array, df), Matrix{Float64})
    end
    @test convert(Array, df) == convert(Array, convert(NullableArray, df))
    @test isa(convert(Array{Any}, df), Matrix{Any})
    @test isa(convert(Array{Float64}, df), Matrix{Float64})

    df = DataFrame()
    df[:A] = 1.0:5.0
    df[:B] = 1.0:5.0
    a = convert(Array, df)
    aa = convert(Array{Any}, df)
    ai = convert(Array{Int}, df)
    @test isa(a, Matrix{Float64})
    @test a == convert(Array, convert(NullableArray, df))
    @test a == convert(Matrix, df)
    @test isa(aa, Matrix{Any})
    @test aa == convert(Matrix{Any}, df)
    @test isa(ai, Matrix{Int})
    @test ai == convert(Matrix{Int}, df)

    df[1,1] = Nullable()
    @test_throws ErrorException convert(Array, df)
    na = convert(NullableArray, df)
    naa = convert(NullableArray{Any}, df)
    nai = convert(NullableArray{Int}, df)
    @test isa(na, NullableMatrix{Float64})
    @test isequal(na, convert(NullableMatrix, df))
    @test isa(naa, NullableMatrix{Any})
    @test isequal(naa, convert(NullableMatrix{Any}, df))
    @test isa(nai, NullableMatrix{Int})
    @test isequal(nai, convert(NullableMatrix{Int}, df))

    a = [1.0,2.0]
    b = [-0.1,3]
    c = [-3.1,7]
    di = Dict("a"=>a, "b"=>b, "c"=>c)

    df = convert(DataFrame,di)
    @test isa(df,DataFrame)
    @test names(df) == Symbol[x for x in sort(collect(keys(di)))]
    @test isequal(df[:a], NullableArray(a))
    @test isequal(df[:b], NullableArray(b))
    @test isequal(df[:c], NullableArray(c))

    od = OrderedDict("c"=>c, "a"=>a, "b"=>b)
    df = convert(DataFrame,od)
    @test isa(df, DataFrame)
    @test names(df) == Symbol[x for x in keys(od)]
    @test isequal(df[:a], NullableArray(a))
    @test isequal(df[:b], NullableArray(b))
    @test isequal(df[:c], NullableArray(c))

    sd = SortedDict("c"=>c, "a"=>a, "b"=>b)
    df = convert(DataFrame,sd)
    @test isa(df, DataFrame)
    @test names(df) == Symbol[x for x in keys(sd)]
    @test isequal(df[:a], NullableArray(a))
    @test isequal(df[:b], NullableArray(b))
    @test isequal(df[:c], NullableArray(c))

    a = [1.0]
    di = Dict("a"=>a, "b"=>b, "c"=>c)
    @test_throws ArgumentError convert(DataFrame,di)

end


text:
using JLD, TimeSeries

include(joinpath(dirname(@__FILE__),"../src/timearray.jl"))

r  = load(joinpath(dirname(@__FILE__),"resp.jld"), "resp")
ta = timearray(r)["Close"]

facts("timearray works on Request object") do

  context("there are 30 rows") do
      @fact  length(ta) --> 30
  end

  context("correct values") do
      @fact  ta[1][1].values[1]  --> 115.01
      @fact  ta[30][1].values[1] --> 109.95
  end
end

facts("Request object inspection check") do

    context("status is 200") do
        @fact  r.status --> 200
    end

    context("data is UInt8") do
        @fact  r.data[1] --> 0x44
    end

# facts("floats and NaNs present") do
#
#   context("NaN fills in missing values slot") do
# #      @pending isnan(sum(md.values[:,1])) --> true
#   end
#
#   context("existing values remain floats") do
# #      @fact  sum(md.values[:,2]) --> 66.86
#   end
#end
end


text:
@testset "execution" begin

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

@eval exec_dummy() = nothing

@testset "cufunction" begin
    @test_throws UndefVarError cufunction(exec_undefined_kernel, ())

    cufunction(dev, exec_dummy, Tuple{})

    # NOTE: other cases are going to be covered by tests below,
    #       as @cuda internally uses cufunction
end

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

@testset "@cuda" begin

@test_throws UndefVarError @cuda exec_undefined_kernel()
@test_throws MethodError @cuda exec_dummy(1)

@testset "compilation params" begin
    @cuda exec_dummy()

    @test_throws CuError @cuda threads=2 maxthreads=1 exec_dummy()
    @cuda threads=2 exec_dummy()
end

@testset "reflection" begin
    CUDAnative.code_lowered(exec_dummy, Tuple{})
    CUDAnative.code_typed(exec_dummy, Tuple{})
    CUDAnative.code_warntype(devnull, exec_dummy, Tuple{})
    CUDAnative.code_llvm(devnull, exec_dummy, Tuple{})
    CUDAnative.code_ptx(devnull, exec_dummy, Tuple{})
    CUDAnative.code_sass(devnull, exec_dummy, Tuple{})

    @device_code_lowered @cuda exec_dummy()
    @device_code_typed @cuda exec_dummy()
    @device_code_warntype io=devnull @cuda exec_dummy()
    @device_code_llvm io=devnull @cuda exec_dummy()
    @device_code_ptx io=devnull @cuda exec_dummy()
    @device_code_sass io=devnull @cuda exec_dummy()

    mktempdir() do dir
        @device_code dir=dir @cuda exec_dummy()
    end

    @test_throws ErrorException @device_code_lowered nothing
end

@testset "shared memory" begin
    @cuda shmem=1 exec_dummy()
end

@testset "streams" begin
    s = CuStream()
    @cuda stream=s exec_dummy()
end

@testset "external kernels" begin
    @eval module KernelModule
        export exec_external_dummy
        exec_external_dummy() = nothing
    end
    import ...KernelModule
    @cuda KernelModule.exec_external_dummy()
    @eval begin
        using ...KernelModule
        @cuda exec_external_dummy()
    end

    @eval module WrapperModule
        using CUDAnative
        @eval exec_dummy() = nothing
        wrapper() = @cuda exec_dummy()
    end
    WrapperModule.wrapper()
end

@testset "calling device function" begin
    @eval @noinline exec_devfun_child(i) = sink(i)
    @eval exec_devfun_parent() = exec_devfun_child(1)

    @cuda exec_devfun_parent()
end

end

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

@testset "argument passing" begin

dims = (16, 16)
len = prod(dims)

@testset "manually allocated" begin
    @eval function exec_pass_ptr(input, output)
        i = (blockIdx().x-1) * blockDim().x + threadIdx().x

        val = unsafe_load(input, i)
        unsafe_store!(output, val, i)
    end

    input = round.(rand(Float32, dims) * 100)
    output = similar(input)

    input_dev = Mem.upload(input)
    output_dev = Mem.alloc(input)

    @cuda threads=len exec_pass_ptr(Base.unsafe_convert(Ptr{Float32}, input_dev),
                                    Base.unsafe_convert(Ptr{Float32}, output_dev))
    Mem.download!(output, output_dev)
    @test input ≈ output
end

@testset "scalar through single-value array" begin
    @eval function exec_pass_scalar(a, x)
        i = (blockIdx().x-1) * blockDim().x + threadIdx().x
        max = gridDim().x * blockDim().x
        if i == max
            val = unsafe_load(a, i)
            unsafe_store!(x, val)
        end
    end

    arr = round.(rand(Float32, dims) * 100)
    val = [0f0]

    arr_dev = Mem.upload(arr)
    val_dev = Mem.upload(val)

    @cuda threads=len exec_pass_scalar(Base.unsafe_convert(Ptr{Float32}, arr_dev),
                                       Base.unsafe_convert(Ptr{Float32}, val_dev))
    @test arr[dims...] ≈ Mem.download(eltype(val), val_dev)[1]
end

@testset "scalar through single-value array, using device function" begin
    @eval @noinline function exec_pass_scalar_devfun(a, x)
        i = (blockIdx().x-1) * blockDim().x + threadIdx().x
        max = gridDim().x * blockDim().x
        if i == max
            val = exec_pass_scalar_devfun_child(a, i)
            unsafe_store!(x, val)
        end
    end
    @eval function exec_pass_scalar_devfun_child(a, i)
        return unsafe_load(a, i)
    end

    arr = round.(rand(Float32, dims) * 100)
    val = [0f0]

    arr_dev = Mem.upload(arr)
    val_dev = Mem.upload(val)

    @cuda threads=len exec_pass_scalar_devfun(Base.unsafe_convert(Ptr{Float32}, arr_dev),
                                              Base.unsafe_convert(Ptr{Float32}, val_dev))
    @test arr[dims...] ≈ Mem.download(eltype(val), val_dev)[1]
end

@testset "tuples" begin
    # issue #7: tuples not passed by pointer

    @eval function exec_pass_tuples(keeps, out)
        if keeps[1]
            unsafe_store!(out, 1)
        else
            unsafe_store!(out, 2)
        end
    end

    keeps = (true,)
    d_out = Mem.alloc(Int)

    @cuda exec_pass_tuples(keeps, Base.unsafe_convert(Ptr{Int}, d_out))
    @test Mem.download(Int, d_out) == [1]
end

@testset "ghost function parameters" begin
    # bug: ghost type function parameters are elided by the compiler

    len = 60
    a = rand(Float32, len)
    b = rand(Float32, len)
    c = similar(a)

    d_a = Mem.upload(a)
    d_b = Mem.upload(b)
    d_c = Mem.alloc(c)

    @eval struct ExecGhost end

    @eval function exec_pass_ghost(ghost, a, b, c)
        i = (blockIdx().x-1) * blockDim().x + threadIdx().x
        unsafe_store!(c, unsafe_load(a,i)+unsafe_load(b,i), i)
    end
    @cuda threads=len exec_pass_ghost(ExecGhost(),
                                      Base.unsafe_convert(Ptr{Float32}, d_a),
                                      Base.unsafe_convert(Ptr{Float32}, d_b),
                                      Base.unsafe_convert(Ptr{Float32}, d_c))
    Mem.download!(c, d_c)
    @test a+b == c

    # bug: ghost type function parameters confused aggregate type rewriting

    @eval function exec_pass_ghost_aggregate(ghost, out, aggregate)
        i = (blockIdx().x-1) * blockDim().x + threadIdx().x
        unsafe_store!(out, aggregate[1], i)
    end
    @cuda threads=len exec_pass_ghost_aggregate(ExecGhost(),
                                                Base.unsafe_convert(Ptr{Float32}, d_c),
                                                (42,))

    Mem.download!(c, d_c)
    @test all(val->val==42, c)
end

@testset "immutables" begin
    # issue #15: immutables not passed by pointer

    @eval function exec_pass_immutables(ptr, b)
        unsafe_store!(ptr, imag(b))
    end

    buf = Mem.upload([0f0])
    x = ComplexF32(2,2)

    @cuda exec_pass_immutables(Base.unsafe_convert(Ptr{Float32}, buf), x)
    @test Mem.download(Float32, buf) == [imag(x)]
end

@testset "automatic recompilation" begin
    buf = Mem.alloc(Int)

    @eval exec_265(ptr) = unsafe_store!(ptr, 1)

    @cuda exec_265(Base.unsafe_convert(Ptr{Int}, buf))
    @test Mem.download(Int, buf) == [1]

    @eval exec_265(ptr) = unsafe_store!(ptr, 2)

    @cuda exec_265(Base.unsafe_convert(Ptr{Int}, buf))
    @test Mem.download(Int, buf) == [2]
end

@testset "non-isbits arguments" begin
    @eval exec_pass_nonbits_unused(T, i) = sink(i)
    @cuda exec_pass_nonbits_unused(Int, 1)

    @eval exec_pass_nonbits_specialized(T, i) = sink(unsafe_trunc(T,i))
    @cuda exec_pass_nonbits_specialized(Int, 1.)

    @eval exec_pass_nonbits_used(i) = sink(unsafe_trunc(Int,i))
    @test_throws CUDAnative.CompilerError @cuda exec_pass_nonbits_used(big"1")
end

@testset "splatting" begin
    @eval function exec_splat(out, a, b)
        unsafe_store!(out, a+b)
    end

    out = [0]
    out_dev = Mem.upload(out)
    out_ptr = Base.unsafe_convert(Ptr{eltype(out)}, out_dev)

    @cuda exec_splat(out_ptr, 1, 2)
    @test Mem.download(eltype(out), out_dev)[1] == 3

    all_splat = (out_ptr, 3, 4)
    @cuda exec_splat(all_splat...)
    @test Mem.download(eltype(out), out_dev)[1] == 7

    partial_splat = (5, 6)
    @cuda exec_splat(out_ptr, partial_splat...)
    @test Mem.download(eltype(out), out_dev)[1] == 11
end

@testset "object invoke" begin
    # this mimics what is generated by closure conversion

    @eval struct exec_object_inner{T} <: Function
        val::T
    end

    @eval function (self::exec_object_inner)(a)
        unsafe_store!(a, self.val)
    end

    @eval function exec_object_outer(a, val)
       inner = exec_object_inner(val)
       @cuda inner(a)
    end

    a = [1.]
    a_dev = Mem.upload(a)

    exec_object_outer(Base.unsafe_convert(Ptr{Float64}, a_dev), 2.)

    @test Mem.download(eltype(a), a_dev) ≈ [2.]
end

@testset "closures" begin
    @eval function exec_closure_outer(a_dev, val)
       function exec_closure_inner(a)
            # captures `val`
            unsafe_store!(a, val)
       end
       @cuda exec_closure_inner(Base.unsafe_convert(Ptr{Float64}, a_dev))
    end

    a = [1.]
    a_dev = Mem.upload(a)

    exec_closure_outer(a_dev, 2.)

    @test Mem.download(eltype(a), a_dev) ≈ [2.]
end

@testset "conversions" begin
    @eval struct Host   end
    @eval struct Device end

    @eval CUDAnative.cudaconvert(a::Host) = Device()

    Base.convert(::Type{Int}, ::Host)   = 1
    Base.convert(::Type{Int}, ::Device) = 2

    out = [0]

    # convert arguments
    out_dev = Mem.upload(out)
    let arg = Host()
        @test Mem.download(eltype(out), out_dev) ≈ [0]
        @eval exec_convert_plain(arg, out) = unsafe_store!(out, convert(Int, arg))
        @cuda exec_convert_plain(arg, Base.unsafe_convert(Ptr{Int}, out_dev))
        @test Mem.download(eltype(out), out_dev) ≈ [2]
    end

    # convert tuples
    out_dev = Mem.upload(out)
    let arg = (Host(),)
        @test Mem.download(eltype(out), out_dev) ≈ [0]
        @eval exec_convert_tuple(arg, out) = unsafe_store!(out, convert(Int, arg[1]))
        @cuda exec_convert_tuple(arg, Base.unsafe_convert(Ptr{Int}, out_dev))
        @test Mem.download(eltype(out), out_dev) ≈ [2]
    end

    # convert named tuples
    out_dev = Mem.upload(out)
    let arg = (a=Host(),)
        @test Mem.download(eltype(out), out_dev) ≈ [0]
        @eval exec_convert_namedtuple(arg, out) = unsafe_store!(out, convert(Int, arg.a))
        @cuda exec_convert_namedtuple(arg, Base.unsafe_convert(Ptr{Int}, out_dev))
        @test Mem.download(eltype(out), out_dev) ≈ [2]
    end

    # don't convert structs
    out_dev = Mem.upload(out)
    @eval struct Nested
        a::Host
    end
    let arg = Nested(Host())
        @test Mem.download(eltype(out), out_dev) ≈ [0]
        @eval exec_convert_namedtuple(arg, out) = unsafe_store!(out, convert(Int, arg.a))
        @cuda exec_convert_namedtuple(arg, Base.unsafe_convert(Ptr{Int}, out_dev))
        @test Mem.download(eltype(out), out_dev) ≈ [1]
    end
end

end

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

end


text:
import ..Topologies
using ..RecursiveApply

function dss_1d!(
    dest,
    src,
    htopology::Topologies.AbstractTopology,
    Nq::Integer,
    Nv::Integer = 1,
)
    idx1 = CartesianIndex(1, 1, 1, 1, 1)
    idx2 = CartesianIndex(Nq, 1, 1, 1, 1)
    for (elem1, face1, elem2, face2, reversed) in
        Topologies.interior_faces(htopology)
        for level in 1:Nv
            @assert face1 == 1 && face2 == 2 && !reversed
            src_slab1 = slab(src, level, elem1)
            src_slab2 = slab(src, level, elem2)
            val = src_slab1[idx1] ⊞ src_slab2[idx2]
            dest_slab1 = slab(dest, level, elem1)
            dest_slab2 = slab(dest, level, elem2)
            dest_slab1[idx1] = dest_slab2[idx2] = val
        end
    end
    return dest
end

"""
    horizontal_dss!(dest, src, topology, Nq)

Apply horizontal direct stiffness summation (DSS) to `src`, storing the result in `dest`.
"""
function dss_2d!(
    dest,
    src,
    topology::Topologies.AbstractTopology,
    Nq::Integer,
    Nv::Integer = 1,
)
    # TODO: generalize to extruded domains by returning a cartesian index?
    # iterate over the interior faces for each element of the mesh
    for (elem1, face1, elem2, face2, reversed) in
        Topologies.interior_faces(topology)
        for level in 1:Nv
            # iterate over non-vertex nodes
            src_slab1 = slab(src, level, elem1)
            src_slab2 = slab(src, level, elem2)
            dest_slab1 = slab(dest, level, elem1)
            dest_slab2 = slab(dest, level, elem2)
            for q in 2:(Nq - 1)
                i1, j1 = Topologies.face_node_index(face1, Nq, q, false)
                i2, j2 = Topologies.face_node_index(face2, Nq, q, reversed)
                value = src_slab1[i1, j1] ⊞ src_slab2[i2, j2]
                dest_slab1[i1, j1] = dest_slab2[i2, j2] = value
            end
        end
    end

    # iterate over all vertices
    for vertex in Topologies.vertices(topology)
        # for each level
        for level in 1:Nv
            # gather: compute sum over shared vertices
            sum_data = mapreduce(⊞, vertex) do (elem, vertex_num)
                i, j = Topologies.vertex_node_index(vertex_num, Nq)
                src_slab = slab(src, level, elem)
                src_slab[i, j]
            end

            # scatter: assign sum to shared vertices
            for (elem, vertex_num) in vertex
                dest_slab = slab(dest, level, elem)
                i, j = Topologies.vertex_node_index(vertex_num, Nq)
                dest_slab[i, j] = sum_data
            end
        end
    end
    return dest
end

function horizontal_dss!(dest, src, space::AbstractSpace)
    if space isa ExtrudedFiniteDifferenceSpace
        Nv = nlevels(space)
        hspace = space.horizontal_space
    else
        Nv = 1
        hspace = space
    end
    htopology = hspace.topology
    Nq = Quadratures.degrees_of_freedom(hspace.quadrature_style)
    if hspace isa SpectralElementSpace1D
        dss_1d!(dest, src, htopology, Nq, Nv)
    elseif hspace isa SpectralElementSpace2D
        dss_2d!(dest, src, htopology, Nq, Nv)
    end
end

horizontal_dss!(data, space::AbstractSpace) =
    horizontal_dss!(data, data, space::AbstractSpace)

weighted_dss!(dest, src, space::AbstractSpace) = horizontal_dss!(
    dest,
    Base.Broadcast.broadcasted(⊠, src, space.dss_weights),
    space,
)

weighted_dss!(dest, src, space::ExtrudedFiniteDifferenceSpace) =
    horizontal_dss!(
        dest,
        Base.Broadcast.broadcasted(⊠, src, space.horizontal_space.dss_weights),
        space,
    )

weighted_dss!(data, space::AbstractSpace) = weighted_dss!(data, data, space)


text:
"""
    TransformedOutput(f, init; tspan::AbstractRange, kw...)

An output that stores the result of some function `f` of the grid/s

# Arguments

- `f`: a function or functor that accepts an `AbstractArray` or `NamedTuple` of
    `AbstractArray` with names matchin `init`. The `AbstractArray` will be a view into
    the grid the same size as the init grids, removing any padding that has been added.
- `init`: initialisation `Array` or `NamedTuple` of `Array`

# Keywords

- `tspan`: `AbstractRange` timespan for the simulation
- `aux`: NamedTuple of arbitrary input data. Use `get(data, Aux(:key), I...)`
    to access from a `Rule` in a type-stable way.
- `mask`: `BitArray` for defining cells that will/will not be run.
- `padval`: padding value for grids with neighborhood rules. The default is `zero(eltype(init))`.
"""
mutable struct TransformedOutput{T,A<:AbstractVector{T},E,F,B} <: Output{T,A}
    frames::A
    running::Bool
    extent::E
    f::F
    buffer::B
end
function TransformedOutput(f::Function, init::Union{NamedTuple,AbstractMatrix}; extent=nothing, kw...)
    extent = extent isa Nothing ? Extent(; init=init, kw...) : extent
    buffer = init isa NamedTuple ? map(zero, init) : zero(init)
    frames = append!([f(init)], map(_ -> f(buffer), tspan(extent)))
    TransformedOutput(frames, false, extent, f, buffer)
end
function TransformedOutput(init; kw...)
    throw(ArgumentError("TransformedOutput must be passed a function and the init grid(s) as arguments"))
end

function storeframe!(o::TransformedOutput, data::AbstractSimData)
    o[frameindex(o, data)] = _store_x(o, grids(data))
end

function _store_x(o::TransformedOutput, grids::NamedTuple)
    # Make a new named tuple of raw arrays without wrappers, copying
    # to the buffer where an OffsetArray was used for padding
    # Often it's faster to copy than use a view when f is sum/mean etc
    nt = map(grids, o.buffer) do g, b
        source(g) isa OffsetArray ? copy!(b, g) : parent(source(g))
    end
    o.f(nt)
end
function _store_x(o::TransformedOutput, grids::NamedTuple{(:_default_,)})
    g = first(grids)
    A = source(g) isa OffsetArray ? copy!(o.buffer, g) : parent(source(g))
    o.f(A)
end

init_output_grids!(o::TransformedOutput, init) = nothing
_initdata!(o::TransformedOutput, init) = nothing


text:
function warmup_cache(cfg::Config,
                      projs_and_mnfsts::AbstractVector{T}) where T <:ProjectFileOrManifestFile
    if isempty(projs_and_mnfsts)
        @warn "No project files or manifest files were provided"
    end
    for proj_or_mnfst in projs_and_mnfsts
        temporary_environment = create_new_temporary_environment(cfg, proj_or_mnfst)
        commands = generate_julia_pkg_commands(
            cfg,
            temporary_environment,
            proj_or_mnfst,
        )
        run_commands(commands)
    end
    return nothing
end


text:
# *********************************************************************************
# REopt, Copyright (c) 2019-2020, Alliance for Sustainable Energy, LLC.
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without modification,
# are permitted provided that the following conditions are met:
#
# Redistributions of source code must retain the above copyright notice, this list
# of conditions and the following disclaimer.
#
# Redistributions in binary form must reproduce the above copyright notice, this
# list of conditions and the following disclaimer in the documentation and/or other
# materials provided with the distribution.
#
# Neither the name of the copyright holder nor the names of its contributors may be
# used to endorse or promote products derived from this software without specific
# prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
# IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
# OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
# OF THE POSSIBILITY OF SUCH DAMAGE.
# *********************************************************************************
"""
    reopt_results(m::JuMP.AbstractModel, p::REoptInputs; _n="")

Create a dictionary of results with string keys for each Scenario structure modeled.
"""
function reopt_results(m::JuMP.AbstractModel, p::REoptInputs; _n="")
	tstart = time()
    d = Dict{String, Any}()
    for b in p.s.storage.types
        if p.s.storage.max_kw[b] > 0
            add_storage_results(m, p, d, b; _n)
        end
    end

    add_electric_tariff_results(m, p, d; _n)
    add_electric_utility_results(m, p, d; _n)
    add_financial_results(m, p, d; _n)
    add_electric_load_results(m, p, d; _n)

	if !isempty(p.techs.pv)
        add_pv_results(m, p, d; _n)
	end

    if "Wind" in p.techs.all
        add_wind_results(m, p, d; _n)
    end

	time_elapsed = time() - tstart
	@info "Base results processing took $(round(time_elapsed, digits=3)) seconds."

	if !isempty(p.techs.gen) && isempty(_n)  # generators not included in multinode model
        tstart = time()
		add_generator_results(m, p, d)
        time_elapsed = time() - tstart
        @info "Generator results processing took $(round(time_elapsed, digits=3)) seconds."
	end

	if !isempty(p.s.electric_utility.outage_durations) && isempty(_n)  # outages not included in multinode model
        tstart = time()
		add_outage_results(m, p, d)
        time_elapsed = time() - tstart
        @info "Outage results processing took $(round(time_elapsed, digits=3)) seconds."
	end

    if !isempty(p.techs.boiler)
        add_existing_boiler_results(m, p, d)
    end

	return d
end

"""
    combine_results(bau::Dict, opt::Dict)

Combine two results dictionaries into one using BAU and optimal scenario results.
"""
function combine_results(p::REoptInputs, bau::Dict, opt::Dict, bau_scenario::BAUScenario)
    bau_outputs = (
        ("Financial", "lcc"),
        ("ElectricTariff", "year_one_energy_cost"),
        ("ElectricTariff", "year_one_demand_cost"),
        ("ElectricTariff", "year_one_fixed_cost"),
        ("ElectricTariff", "year_one_min_charge_adder"),
        ("ElectricTariff", "lifecycle_energy_cost"),
        ("ElectricTariff", "lifecycle_demand_cost"),
        ("ElectricTariff", "lifecycle_fixed_cost"),
        ("ElectricTariff", "lifecycle_min_charge_adder"),
        ("ElectricTariff", "lifecycle_export_benefit"),
        ("ElectricTariff", "year_one_bill"),
        ("ElectricTariff", "year_one_export_benefit"),
        ("ElectricTariff", "year_one_coincident_peak_cost"),
        ("ElectricTariff", "lifecycle_coincident_peak_cost"),
        ("ElectricUtility", "year_one_to_load_series_kw"),
        ("ElectricUtility", "year_one_energy_supplied_kwh"),
        ("PV", "average_annual_energy_produced_kwh"),
        ("PV", "year_one_energy_produced_kwh"),
        ("PV", "lifecycle_om_cost"),
        ("Generator", "fuel_used_gal"),
        ("Generator", "lifecycle_fixed_om_cost"),
        ("Generator", "lifecycle_variable_om_cost"),
        ("Generator", "lifecycle_fuel_cost"),
        ("Generator", "year_one_fuel_cost"),
        ("Generator", "year_one_variable_om_cost"),
        ("Generator", "year_one_fixed_om_cost"),
    )

    for t in bau_outputs
        if t[1] in keys(opt) && t[1] in keys(bau)
            if t[2] in keys(bau[t[1]])
                opt[t[1]][t[2] * "_bau"] = bau[t[1]][t[2]]
            end
        elseif t[1] == "PV" && !isempty(p.techs.pv)
            for pvname in p.techs.pv
                if pvname in keys(opt) && pvname in keys(bau)
                    if t[2] in keys(bau[pvname])
                        opt[pvname][t[2] * "_bau"] = bau[pvname][t[2]]
                    end
                end
            end
        end
    end
    opt["Financial"]["lifecycle_om_costs_bau"] = bau["Financial"]["lifecycle_om_costs_after_tax"]
    opt["Financial"]["npv"] = round(opt["Financial"]["lcc_bau"] - opt["Financial"]["lcc"], digits=2)

    opt["ElectricLoad"]["bau_critical_load_met"] = bau_scenario.outage_outputs.bau_critical_load_met
    opt["ElectricLoad"]["bau_critical_load_met_time_steps"] = bau_scenario.outage_outputs.bau_critical_load_met_time_steps

    return opt
end

max_stars_repo_name: NREL/REoptLite

max_issues_repo_name: NREL/REoptLite

max_forks_repo_name: NREL/REoptLite


text:
module Ryu

include("utils.jl")
include("shortest.jl")
include("fixed.jl")
include("exp.jl")

neededdigits(::Type{Float64}) = 309 + 17
neededdigits(::Type{Float32}) = 39 + 9
neededdigits(::Type{Float16}) = 9 + 5

"""
    Ryu.writeshortest(x, plus=false, space=false, hash=true, precision=-1, expchar=UInt8('e'), padexp=false, decchar=UInt8('.'))
    Ryu.writeshortest(buf::Vector{UInt8}, pos::Int, x, args...)

Convert a float value `x` into its "shortest" decimal string, which can be parsed back to the same value.
This function allows achieving the `%g` printf format.
Note the 2nd method allows passing in a byte buffer and position directly; callers must ensure the buffer has sufficient room to hold the entire decimal string.

Various options for the output format include:
  * `plus`: for positive `x`, prefix decimal string with a `'+'` character
  * `space`: for positive `x`, prefix decimal string with a `' '` character; overridden if `plus=true`
  * `hash`: whether the decimal point should be written, even if no additional digits are needed for precision
  * `precision`: minimum number of significant digits to be included in the decimal string; extra `'0'` characters will be added for padding if necessary
  * `expchar`: character to use exponent component in scientific notation
  * `padexp`: whether two digits should always be written, even for single-digit exponents (e.g. `e+1` becomes `e+01`)
  * `decchar`: decimal point character to be used
"""
function writeshortest(x::T,
        plus::Bool=false,
        space::Bool=false,
        hash::Bool=true,
        precision::Integer=-1,
        expchar::UInt8=UInt8('e'),
        padexp::Bool=false,
        decchar::UInt8=UInt8('.')) where {T <: Base.IEEEFloat}
    buf = Base.StringVector(neededdigits(T))
    pos = writeshortest(buf, 1, x)
    @assert pos - 1 <= length(buf)
    return String(resize!(buf, pos - 1))
end

"""
    Ryu.writefixed(x, plus=false, space=false, hash=true, precision=-1, decchar=UInt8('.'))
    Ryu.writefixed(buf::Vector{UInt8}, pos::Int, x, args...)

Convert a float value `x` into a "fixed" size decimal string.
This function allows achieving the `%f` printf format.
Note the 2nd method allows passing in a byte buffer and position directly; callers must ensure the buffer has sufficient room to hold the entire decimal string.

Various options for the output format include:
  * `plus`: for positive `x`, prefix decimal string with a `'+'` character
  * `space`: for positive `x`, prefix decimal string with a `' '` character; overridden if `plus=true`
  * `hash`: whether the decimal point should be written, even if no additional digits are needed for precision
  * `precision`: minimum number of significant digits to be included in the decimal string; extra `'0'` characters will be added for padding if necessary
  * `decchar`: decimal point character to be used
"""
function writefixed(x::T, precision) where {T <: Base.IEEEFloat}
    buf = Base.StringVector(precision + neededdigits(T))
    pos = writefixed(buf, 1, x, false, false, false, precision)
    @assert pos - 1 <= length(buf)
    return String(resize!(buf, pos - 1))
end

"""
    Ryu.writeexp(x, plus=false, space=false, hash=true, precision=-1, expchar=UInt8('e'), decchar=UInt8('.'))
    Ryu.writeexp(buf::Vector{UInt8}, pos::Int, x, args...)

Convert a float value `x` into a scientific notation decimal string.
This function allows achieving the `%e` printf format.
Note the 2nd method allows passing in a byte buffer and position directly; callers must ensure the buffer has sufficient room to hold the entire decimal string.

Various options for the output format include:
  * `plus`: for positive `x`, prefix decimal string with a `'+'` character
  * `space`: for positive `x`, prefix decimal string with a `' '` character; overridden if `plus=true`
  * `hash`: whether the decimal point should be written, even if no additional digits are needed for precision
  * `precision`: minimum number of significant digits to be included in the decimal string; extra `'0'` characters will be added for padding if necessary
  * `expchar`: character to use exponent component in scientific notation
  * `decchar`: decimal point character to be used
"""
function writeexp(x::T, precision) where {T <: Base.IEEEFloat}
    buf = Base.StringVector(precision + neededdigits(T))
    pos = writeexp(buf, 1, x, false, false, false, precision)
    @assert pos - 1 <= length(buf)
    return String(resize!(buf, pos - 1))
end

function Base.show(io::IO, x::T) where {T <: Base.IEEEFloat}
    if get(io, :compact, false)
        x = round(x, sigdigits=6)
    end
    buf = Base.StringVector(neededdigits(T))
    pos = writeshortest(buf, 1, x)
    @assert pos - 1 <= length(buf)
    write(io, resize!(buf, pos - 1))
    return
end

end # module

max_stars_repo_name: mschauer/julia

max_issues_repo_name: mschauer/julia

max_forks_repo_name: mschauer/julia


text:
export GoalSpec

"Goal specification for a planning problem."
@kwdef struct GoalSpec
    goals::Vector{Term} = Term[] # Goal terms to be satisfied
    metric::Union{Term,Nothing} = nothing # Metric to be minimized
    constraints::Vector{Term} = Term[] # Trajectory constraints
    GoalSpec(goals, metric, constraints) =
        new(flatten_conjs(goals), metric, flatten_conjs(constraints))
end

# Constructors
GoalSpec(goals::Vector{<:Term}) = GoalSpec(goals=goals)
GoalSpec(goal::Term) = GoalSpec(goals=flatten_conjs([goal]))
function GoalSpec(problem::Problem)
    goals = flatten_conjs(problem.goal)
    sign, metric = problem.metric
    if sign > 0 metric = Compound(:-, [metric]) end
    return GoalSpec(goals=goals, metric=metric)
end


text:
module PixelsGT
    include("Krav.jl")
    include("Wignerd.jl")
    include("LKrav.jl")
    include("ImgRots.jl")
    include("ImgAcq.jl")
end


text:
#####
##### General halo filling functions
#####

fill_halo_regions!(::Nothing, args...) = nothing

"""
    fill_halo_regions!(fields, arch)

Fill halo regions for each field in the tuple `fields` according to their boundary
conditions, possibly recursing into `fields` if it is a nested tuple-of-tuples.
"""
function fill_halo_regions!(fields::Union{Tuple,NamedTuple}, arch, args...)
    for field in fields
        fill_halo_regions!(field, arch, args...)
    end
    return nothing
end

fill_halo_regions!(field, arch, args...) =
    fill_halo_regions!(field.data, field.boundary_conditions, arch, field.grid, args...)

"Fill halo regions in x, y, and z for a given field."
function fill_halo_regions!(c::AbstractArray, fieldbcs, arch, grid, args...)
      fill_west_halo!(c, fieldbcs.x.left,   arch, grid, args...)
      fill_east_halo!(c, fieldbcs.x.right,  arch, grid, args...)
     fill_south_halo!(c, fieldbcs.y.left,   arch, grid, args...)
     fill_north_halo!(c, fieldbcs.y.right,  arch, grid, args...)
    fill_bottom_halo!(c, fieldbcs.z.bottom, arch, grid, args...)
       fill_top_halo!(c, fieldbcs.z.top,    arch, grid, args...)
    return nothing
end

#####
##### Halo filling for flux, periodic, and no-penetration boundary conditions.
#####

# For flux boundary conditions we fill halos as for a *no-flux* boundary condition, and add the
# flux divergence associated with the flux boundary condition in a separate step. Note that
# ranges are used to reference the data copied into halos, as this produces views of the correct
# dimension (eg size = (1, Ny, Nz) for the west halos).

  _fill_west_halo!(c, ::FBC, H, N) = @views @. c.parent[1:H, :, :] = c.parent[1+H:1+H,  :, :]
 _fill_south_halo!(c, ::FBC, H, N) = @views @. c.parent[:, 1:H, :] = c.parent[:, 1+H:1+H,  :]
_fill_bottom_halo!(c, ::FBC, H, N) = @views @. c.parent[:, :, 1:H] = c.parent[:, :,  1+H:1+H]

 _fill_east_halo!(c, ::FBC, H, N) = @views @. c.parent[N+H+1:N+2H, :, :] = c.parent[N+H:N+H, :, :]
_fill_north_halo!(c, ::FBC, H, N) = @views @. c.parent[:, N+H+1:N+2H, :] = c.parent[:, N+H:N+H, :]
  _fill_top_halo!(c, ::FBC, H, N) = @views @. c.parent[:, :, N+H+1:N+2H] = c.parent[:, :, N+H:N+H]

# Periodic boundary conditions
  _fill_west_halo!(c, ::PBC, H, N) = @views @. c.parent[1:H, :, :] = c.parent[N+1:N+H, :, :]
 _fill_south_halo!(c, ::PBC, H, N) = @views @. c.parent[:, 1:H, :] = c.parent[:, N+1:N+H, :]
_fill_bottom_halo!(c, ::PBC, H, N) = @views @. c.parent[:, :, 1:H] = c.parent[:, :, N+1:N+H]

 _fill_east_halo!(c, ::PBC, H, N) = @views @. c.parent[N+H+1:N+2H, :, :] = c.parent[1+H:2H, :, :]
_fill_north_halo!(c, ::PBC, H, N) = @views @. c.parent[:, N+H+1:N+2H, :] = c.parent[:, 1+H:2H, :]
  _fill_top_halo!(c, ::PBC, H, N) = @views @. c.parent[:, :, N+H+1:N+2H] = c.parent[:, :, 1+H:2H]

# Recall that, by convention, the first grid point (k=1) in an array with a no-penetration boundary
# condition lies on the boundary, where as the last grid point (k=Nz) lies in the domain.

  _fill_west_halo!(c, ::NPBC, H, N) = @views @. c.parent[1:1+H, :, :] = 0
 _fill_south_halo!(c, ::NPBC, H, N) = @views @. c.parent[:, 1:1+H, :] = 0
_fill_bottom_halo!(c, ::NPBC, H, N) = @views @. c.parent[:, :, 1:1+H] = 0

 _fill_east_halo!(c, ::NPBC, H, N) = @views @. c.parent[N+H+1:N+2H, :, :] = 0
_fill_north_halo!(c, ::NPBC, H, N) = @views @. c.parent[:, N+H+1:N+2H, :] = 0
  _fill_top_halo!(c, ::NPBC, H, N) = @views @. c.parent[:, :, N+H+1:N+2H] = 0

# Generate functions that implement flux, periodic, and no-penetration boundary conditions
sides = (:west, :east, :south, :north, :top, :bottom)
coords = (:x, :x, :y, :y, :z, :z)

for (x, side) in zip(coords, sides)
    outername = Symbol(:fill_, side, :_halo!)
    innername = Symbol(:_fill_, side, :_halo!)
    H = Symbol(:H, x)
    N = Symbol(:N, x)
    @eval begin
        $outername(c, bc::Union{FBC, PBC, NPBC}, arch, grid, args...) =
            $innername(c, bc, grid.$(H), grid.$(N))
    end
end


text:
export minkowski_distance, euclidian_distance, manhattan_distance, chebyshev_distance, braycurtis_distance, canberra_distance, cityblock_distance, mahalanobis_distance,
correlation, cosine_distance, cosine_similarity

using LinearAlgebra
using Statistics: mean

#TODO
# scikitlearn-like kernels ?

# Distance methods

"""
```minkowski_distance(u,v; keywords)```
```minkowski_distance(x::Tuple; keywords)```

Return the Minkowski distance between the given points both of which are given with vector-like structures.

## Arguments
- `u` : first point
- `v` : second point
- `x::Tuple` : tuple as follows: (u,v)

## Keywords
- `p::Int = 2` : p value
- `w::Int` : weights that must be of the same length with u,v
- `rooted = false` : denotes whether or not to execute the final rooting operation
"""
function minkowski_distance(u,v; p = 2, w = nothing, rooted = false)
    _validate_distance_input(u, v, w; p = p, p_is_used = true)
    val = abs.(u .- v)
    val = w != nothing ? val .* w : val
    if p == 1
        val = sum(val)
    elseif p == Inf
        val = maximum(val)
    else
        val = sum(val .^ p)
    end
    return rooted ? val .^ (1 / p) : val
end

minkowski_distance(x::Tuple; p = 2, w = nothing) = minkowski_distance(x...; p = p, w = w)

"""
```euclidian_distance(u,v; keywords)```
```euclidian_distance(x::Tuple; keywords)```

Return the Euclidian distance between the given points both of which are given with vector-like structures.

## Arguments
- `u` : first point
- `v` : second point
- `x::Tuple` : tuple as follows: (u,v)

## Keywords
- `w::Int` : weights that must be of the same length with u,v
- `squared = false` : denotes whether or not to execute the final squaring operation
"""
euclidian_distance(u,v; w = nothing, squared = false) = minkowski_distance(u,v; w = w, p = 2, rooted = squared)
euclidian_distance(x::Tuple;w = nothing, squared = false) = euclidian_distance(x...; w = w, rooted = squared)

"""
```manhattan_distance(u,v; keywords)```
```manhattan_distance(x::Tuple; keywords)```

Return the Manhattan distance between the given points both of which are given with vector-like structures.

## Arguments
- `u` : first point
- `v` : second point
- `x::Tuple` : tuple as follows: (u,v)

## Keywords
- `w::Int` : weights that must be of the same length with u,v
"""
manhattan_distance(u,v; w = nothing) = minkowski_distance(u,v; w = w, p = 1)
manhattan_distance(x::Tuple;w = nothing) = manhattan_distance(x...; w = w)

"""
```chebyshev_distance(u,v; keywords)```
```chebyshev_distance(x::Tuple; keywords)```

Return the Chebyshev distance between the given points both of which are given with vector-like structures.

## Arguments
- `u` : first point
- `v` : second point
- `x::Tuple` : tuple as follows: (u,v)

## Keywords
- `w::Int` : weights that must be of the same length with u,v
"""
chebyshev_distance(u,v; w = nothing) = minkowski_distance(u,v; w = w, p = Inf, rooted=false)
chebyshev_distance(x::Tuple;w = nothing) = chebyshev_distance(x...; w = w)

"""
```braycurtis_distance(u,v; keywords)```
```braycurtis_distance(x::Tuple; keywords)```

Return the Braycurtis distance between the given points both of which are given with vector-like structures.

## Arguments
- `u` : first point
- `v` : second point
- `x::Tuple` : tuple as follows: (u,v)

## Keywords
- `w::Int` : weights that must be of the same length with u,v
"""
function braycurtis_distance(u,v; w = nothing)
    _validate_distance_input(u, v, w)
    difference = abs.(u .- v)
    summation = abs.(u .+ v)
    if w != nothing
        difference = w .* difference
        summation = w .* summation
    end
    sum(difference) / sum(summation)
end

braycurtis_distance(x::Tuple; w = nothing) = braycurtis_distance(x...; w = w)

"""
```canberra_distance(u,v; keywords)```
```canberra_distance(x::Tuple; keywords)```

Return the Canberra distance between the given points both of which are given with vector-like structures.

## Arguments
- `u` : first point
- `v` : second point
- `x::Tuple` : tuple as follows: (u,v)

## Keywords
- `w::Int` : weights that must be of the same length with u,v
"""
function canberra_distance(u,v; w = nothing)
    _validate_distance_input(u, v, w)
    difference = abs.(u .- v)
    result = difference / (abs.(u) .+ abs.(v))
    if w != nothing
        result = result .* w
    end
    return  sum(result[result .!== Inf])
end
canberra_distance(x::Tuple; w = nothing) = canberra_distance(x...; w = w)

"""
```cityblock_distance(u,v; keywords)```
```cityblock_distance(x::Tuple; keywords)```

Return the Cityblock distance between the given points both of which are given with vector-like structures.

## Arguments
- `u` : first point
- `v` : second point
- `x::Tuple` : tuple as follows: (u,v)

## Keywords
- `w::Int` : weights that must be of the same length with u,v
"""
function cityblock_distance(u,v; w = nothing)
    _validate_distance_input(u,v,w)
    result = abs.(u .- v)
    if w != nothing
        result = result .* w
    end
    return sum(result)
end
cityblock_distance(x; w = nothing) = cityblock_distance(x...; w = w)

"""
```mahalanobis_distance(u,v; keywords)```
```mahalanobis_distance(x::Tuple; keywords)```

Return the Cityblock distance between the given points both of which are given with vector-like structures.

## Arguments
- `u` : first point
- `v` : second point
- `x::Tuple` : tuple as follows: (u,v)
- `Vinv::Matrix` : The inverse of the covariance matrix.

"""
function mahalanobis_distance(u,v, Vinv)
    _validate_distance_input(u,v,nothing)
    delta = u .- v
    return sqrt.(dot(dot(delta, convert_1d(Vinv)), delta))
end

mahalanobis_distance(x::Tuple, Vinv) = mahalanobis_distance(x..., Vinv)

##

"""
```correlation(u,v; keywords)```
```correlation(x::Tuple; keywords)```

Return the Cityblock distance between the given points both of which are given with vector-like structures.

## Arguments
- `u` : first point
- `v` : second point
- `x::Tuple` : tuple as follows: (u,v)

## Keywords
- `w::Int` : weights that must be of the same length with u,v
- `centered=true` : If true, u and v will be centered.
"""
function correlation(u,v; w = nothing, centered = true)
    _validate_distance_input(u,v,w)
    w = w != nothing ? w : ones(length(u))
    _u = convert_1d(u)
    _v = convert_1d(v)
    if centered
        umu = mean(_u .* w)
        vmu = mean(_v .* w)
        _u = _u .- umu
        _v = _v .- vmu
    end
    uv = mean(_u .* _v .* w)
    uu = mean( (_u .^ 2) .* w )
    vv = mean( (_v .^ 2) .* w  )
    return abs(1 - uv / sqrt(uu*vv))
end

cosine_distance(u,v;w=nothing) = correlation(u,v; w = w)
cosine_similarity(u,v; w = nothing) = -cosine_distance(u,v; w = w) + 1