text:
"""
    peakwidths(peaks, x, proms;
        strict=true,
        relheight=0.5,
        minwidth=nothing,
        maxwidth=nothing
    ) -> (peaks, widths, leftedge, rightedge)

Calculate the widths of `peaks` in `x` at a reference level based on `proms` and
`relheight`, removing peaks with widths less than `minwidth` or greater than `maxwidth`, if
either are given. Returns the modified peaks, widths, and the left and right edges at the
reference level.

Peak width is the distance between the signal crossing a reference level before and after
the peak. Signal crossings are linearly interpolated between indices. The reference level is
the difference between the peak height and `relheight` times the peak prominence. Width
cannot be calculated for a `NaN` or `missing` prominence.

The width for a peak with a gap in the signal (e.g. `NaN`, `missing`) at the reference level
will match the value/type of the signal gap if `strict == true`. For `strict ==
false`, the signal crossing will be linearly interpolated between the edges of the gap.

See also: [`peakprom`](@ref), [`findminima`](@ref), [`findmaxima`](@ref)

# Examples
```jldoctest
julia> x = [0,1,0,-1.];

julia> xpks = argmaxima(x)
1-element Vector{Int64}:
 2

julia> peakwidths(xpks, x, [1])
([2], [1.0], [1.5], [2.5])

julia> x[3] = NaN;

julia> peakwidths(xpks, x, [1])
([2], [NaN], [1.5], [NaN])

julia> peakwidths(xpks, x, [1]; strict=false)
([2], [1.0], [1.5], [2.5])
```
"""
function peakwidths(
    peaks::AbstractVector{Int}, x::AbstractVector, proms::AbstractVector;
    strict=true, relheight=0.5, minwidth=nothing, maxwidth=nothing,
)
    if !isnothing(minwidth) || !isnothing(maxwidth)
        _peaks = copy(peaks)
    else
        # peaks will not be modified
        _peaks = peaks
    end
    peakwidths!(_peaks, x, proms; strict=strict, relheight=relheight,
        minwidth=minwidth, maxwidth=maxwidth)
end

"""
    peakwidths!(peaks, x, proms;
        strict=true,
        relheight=0.5,
        minwidth=nothing,
        maxwidth=nothing
    ) -> (peaks, widths, leftedge, rightedge)

Calculate the widths of `peaks` in `x` at a reference level based on `proms` and
`relheight`, removing peaks with widths less than `minwidth` or greater than `maxwidth`, if
either are given. Returns the modified peaks, widths, and the left and right edges at the
reference level.

See also: [`peakwidths`](@ref), [`peakproms`](@ref), [`findminima`](@ref), [`findmaxima`](@ref)
"""
function peakwidths!(
    peaks::AbstractVector{Int}, x::AbstractVector{T}, proms::AbstractVector{U};
    strict=true, relheight=0.5, minwidth=nothing, maxwidth=nothing,
) where {T, U}
    if !isnothing(minwidth) && !isnothing(maxwidth)
        minwidth < maxwidth || throw(ArgumentError("maxwidth must be greater than minwidth"))
    end
    all(∈(eachindex(x)), peaks) ||
        throw(ArgumentError("peaks contains invalid indices to x"))

    # if peaks was calculated with strict=false, first(peaks) could be minima at firstindex
    fp = length(peaks) > 1 ? peaks[2] : first(peaks)
    if fp > 1 && ((x[fp] < x[fp-1]) === true)
        pktype = :minima
    else
        pktype = :maxima
    end
    cmp = pktype === :maxima ? (≤) : (≥)
    op = pktype === :maxima ? (-) : (+)

    V1 = promote_type(T,U)
    _bad = Missing <: V1 ? missing : float(Int)(NaN)

    V = promote_type(V1, typeof(_bad))
    ledge = similar(proms, V)
    redge = similar(proms, V)

    if strict
        lst, fst = _bad, _bad
    else
        lst = lastindex(x)
        fst = firstindex(x)
    end

    for i in eachindex(peaks, ledge, redge)
        prom = proms[i]
        if ismissing(prom) || isnan(prom)
            redge[i] = _bad
            ledge[i] = _bad
        else
            ht = op(x[peaks[i]], relheight*proms[i])
            lo = findprev(v -> !ismissing(v) && cmp(v,ht), x, peaks[i])
            up = findnext(v -> !ismissing(v) && cmp(v,ht), x, peaks[i])

            if !strict
                if !isnothing(lo)
                    lo1 = findnext(v -> !ismissing(v) && cmp(ht,v), x, lo+1)
                    lo += (ht - x[lo])/(x[lo1] - x[lo])*(lo1 - lo)
                end
                if !isnothing(up)
                    up1 = findprev(v -> !ismissing(v) && cmp(ht,v), x, up-1)
                    up -= (ht - x[up])/(x[up1] - x[up])*(up - up1)
                end
            else
                !isnothing(lo) && (lo += (ht - x[lo])/(x[lo+1] - x[lo]))
                !isnothing(up) && (up -= (ht - x[up])/(x[up-1] - x[up]))
            end

            redge[i] = something(up, lst)
            ledge[i] = something(lo, fst)
        end
    end

    widths::Vector{V} = redge - ledge

    if !isnothing(minwidth) || !isnothing(maxwidth)
        lo = something(minwidth, zero(eltype(widths)))
        up = something(maxwidth, typemax(Base.nonmissingtype(eltype(widths))))
        matched = findall(x -> !ismissing(x) && !(lo ≤ x ≤ up), widths)
        deleteat!(peaks, matched)
        deleteat!(ledge, matched)
        deleteat!(redge, matched)
        deleteat!(widths, matched)
    end

    return peaks, widths, ledge, redge
end

max_stars_repo_name: halleysfifthinc/Peaks

max_issues_repo_name: halleysfifthinc/Peaks

max_forks_repo_name: halleysfifthinc/Peaks


text:
using Turing, Test

ProjDir = mktempdir()

@testset "serialization read and write test" begin
    @model gdemo(x) = begin
        s ~ InverseGamma(2,3)
        m ~ Normal(0, sqrt(s))
        for i in eachindex(x)
            x[i] ~ Normal(m, sqrt(s))
        end
    end

    model = gdemo([1.5, 2.0])
    sampler = HMC(0.01, 5)
    chn1 = sample(model, sampler, 500, save_state=true)

    write(joinpath(ProjDir, "chn1.jls"), chn1)
    chn2 = read(joinpath(ProjDir, "chn1.jls"), MCMCChains.Chains)

    open(joinpath(ProjDir, "chn1.txt"), "w") do io
        describe(io, chn1);
    end

    open(joinpath(ProjDir, "chn2.txt"), "w") do io
        describe(io, chn2);
    end

    @test open(f->read(f, String), joinpath(ProjDir, "chn1.txt")) ==
        open(f->read(f, String), joinpath(ProjDir, "chn2.txt"))

    # Test whether sampler state made it through serialization/deserialization.
    chn3 = Turing.Inference.resume(chn2, 100)
    @test range(chn3) == 1:1:100
end

rm(ProjDir, force=true, recursive=true)

ext: jl


text:
for i in 1:5

    data = 10 * randn(440, 1) + 20 * sin(2 * pi * 0.01 * collect(1:440))
    cleaned = data .- Private.savitsky_golay(vec(data), 31, 2)

    @test mean(data[1:50]) > 6
    @test mean(cleaned[1:50]) < 1.5

    @test mean(data[100:150]) > 6
    @test mean(cleaned[100:150]) < 1.5

    # display(plot(data, lab = "raw"))
    # display(plot(data .- Private.savitsky_golay(vec(data), 11, 2), lab = "detrended", c = :red))
end

fname = joinpath(dirname(@__FILE__), "../data", "test_Hz19.5-testing.bdf")

a = read_SSR(fname)
a = detrend(a)


text:
struct UnrollSymbols
    u₁loopsym::Symbol
    u₂loopsym::Symbol
    vectorized::Symbol
end
struct UnrollArgs{T <: Union{Nothing,Int}}
    u₁::Int
    u₁loopsym::Symbol
    u₂loopsym::Symbol
    vectorized::Symbol
    u₂max::Int
    suffix::T
end
function UnrollArgs(u₁::Int, unrollsyms::UnrollSymbols, u₂max::Int, suffix)
    @unpack u₁loopsym, u₂loopsym, vectorized = unrollsyms
    UnrollArgs(u₁, u₁loopsym, u₂loopsym, vectorized, u₂max, suffix)
end
function UnrollArgs(ua::UnrollArgs, u::Int)
    @unpack u₁loopsym, u₂loopsym, vectorized, u₂max, suffix = ua
    UnrollArgs(u, u₁loopsym, u₂loopsym, vectorized, u₂max, suffix)
end
# UnrollSymbols(ua::UnrollArgs) = UnrollSymbols(ua.u₁loopsym, ua.u₂loopsym, ua.vectorized)

# isfirst(ua::UnrollArgs{Nothing}) = iszero(ua.u₁)
# isfirst(ua::UnrollArgs{Int}) = iszero(ua.u₁) & iszero(ua.suffix)

struct UnrollSpecification
    u₁loopnum::Int
    u₂loopnum::Int
    vectorizedloopnum::Int
    u₁::Int
    u₂::Int
end
# UnrollSpecification(ls::LoopSet, u₁loop::Loop, vectorized::Symbol, u₁, u₂) = UnrollSpecification(ls, u₁loop.itersymbol, vectorized, u₁, u₂)
function UnrollSpecification(us::UnrollSpecification, u₁, u₂)
    @unpack u₁loopnum, u₂loopnum, vectorizedloopnum = us
    UnrollSpecification(u₁loopnum, u₂loopnum, vectorizedloopnum, u₁, u₂)
end
# function UnrollSpecification(us::UnrollSpecification; u₁ = us.u₁, u₂ = us.u₂)
#     @unpack u₁loopnum, u₂loopnum, vectorizedloopnum = us
#     UnrollSpecification(u₁loopnum, u₂loopnum, vectorizedloopnum, u₁, u₂)
# end
isunrolled1(us::UnrollSpecification, n::Int) = us.u₁loopnum == n
isunrolled2(us::UnrollSpecification, n::Int) = !isunrolled1(us, n) && us.u₂loopnum == n
isvectorized(us::UnrollSpecification, n::Int) = us.vectorizedloopnum == n
function unrollfactor(us::UnrollSpecification, n::Int)
    @unpack u₁loopnum, u₂loopnum, u₁, u₂ = us
    (u₁loopnum == n) ? u₁ : ((u₂loopnum == n) ? u₂ : 1)
end

struct Loop
    itersymbol::Symbol
    starthint::Int
    stophint::Int
    startsym::Symbol
    stopsym::Symbol
    startexact::Bool
    stopexact::Bool
end
startstopint(s::Int, start) = s
startstopint(s::Symbol, start) = start ? 1 : 1024
startstopsym(s::Int) = Symbol("##UNDEFINED##")
startstopsym(s::Symbol) = s
function Loop(itersymbol::Symbol, start::Union{Int,Symbol}, stop::Union{Int,Symbol})
    Loop(
        itersymbol, startstopint(start,true), startstopint(stop,false),
        startstopsym(start), startstopsym(stop), start isa Int, stop isa Int
    )
end
Base.length(loop::Loop) = 1 + loop.stophint - loop.starthint
isstaticloop(loop::Loop) = loop.startexact & loop.stopexact

function startloop(loop::Loop, itersymbol)
    startexact = loop.startexact
    if startexact
        Expr(:(=), itersymbol, loop.starthint - 1)
    else
        Expr(:(=), itersymbol, Expr(:call, lv(:staticm1), Expr(:call, lv(:unwrap), loop.startsym)))
    end
end
addexpr(ex, incr) = Expr(:call, lv(:vadd), ex, incr)
function addexpr(ex, incr::Number)
    if iszero(incr)
        incr
    elseif incr > 0
        Expr(:call, lv(:vadd), ex, incr)
    else
        Expr(:call, lv(:vsub), ex, -incr)
    end
end
addexpr(ex::Number, incr::Number) = ex + incr
subexpr(ex, incr) = Expr(:call, lv(:vsub), ex, incr)
subexpr(ex::Number, incr::Number) = ex - incr
subexpr(ex, incr::Number) = addexpr(ex,  -incr)

staticmulincr(ptr, incr) = Expr(:call, lv(:staticmul), Expr(:call, :eltype, ptr), incr)
callpointerforcomparison(sym) = Expr(:call, lv(:pointerforcomparison), sym)
function vec_looprange(loopmax, UF::Int, mangledname::Symbol, ptrcomp::Bool)
    if ptrcomp
        vec_looprange(loopmax, UF, callpointerforcomparison(mangledname), staticmulincr(mangledname, VECTORWIDTHSYMBOL))
    else
        vec_looprange(loopmax, UF, mangledname, VECTORWIDTHSYMBOL)
    end
end
function vec_looprange(loopmax, UF::Int, mangledname, W)
    incr = if isone(UF)
        Expr(:call, lv(:valsub), W, 1)
    else
        Expr(:call, lv(:valmulsub), W, UF, 1)
    end
    compexpr = subexpr(loopmax, incr)
    Expr(:call, :<, mangledname, compexpr)
end

function looprange(stopcon, incr::Int, mangledname)
    if iszero(incr)
        Expr(:call, :≤, mangledname, stopcon)
    else
        Expr(:call, :≤, mangledname, subexpr(stopcon, incr))
    end
end
function looprange(loop::Loop, incr::Int, mangledname)
    loop.stopexact ? looprange(loop.stophint, incr, mangledname) : looprange(loop.stopsym, incr, mangledname)
end
function terminatecondition(
    loop::Loop, us::UnrollSpecification, n::Int, mangledname::Symbol, inclmask::Bool, UF::Int = unrollfactor(us, n)
)
    if !isvectorized(us, n)
        looprange(loop, UF, mangledname)
    elseif inclmask
        looprange(loop, 1, mangledname)
    elseif loop.stopexact
        vec_looprange(loop.stophint, UF, mangledname, false) # may not be u₂loop
    else
        vec_looprange(loop.stopsym, UF, mangledname, false) # may not be u₂loop
    end
end
function incrementloopcounter(us::UnrollSpecification, n::Int, mangledname::Symbol, UF::Int = unrollfactor(us, n))
    if isvectorized(us, n)
        if isone(UF)
            Expr(:(=), mangledname, Expr(:call, lv(:valadd), VECTORWIDTHSYMBOL, mangledname))
        else
            Expr(:(=), mangledname, Expr(:call, lv(:valmuladd), VECTORWIDTHSYMBOL, UF, mangledname))
        end
    else
        Expr(:(=), mangledname, Expr(:call, lv(:vadd), mangledname, UF))
    end
end
function incrementloopcounter!(q, us::UnrollSpecification, n::Int, UF::Int = unrollfactor(us, n))
    if isvectorized(us, n)
        if isone(UF)
            push!(q.args, Expr(:call, lv(:Static), VECTORWIDTHSYMBOL))
        else
            push!(q.args, Expr(:call, lv(:valmul), VECTORWIDTHSYMBOL, Expr(:call, Expr(:curly, lv(:Static), UF))))
        end
    else
        push!(q.args, Expr(:call, Expr(:curly, lv(:Static), UF)))
    end
end
function looplengthexpr(loop::Loop)
    if loop.stopexact
        if loop.startexact
            length(loop)
        else
            Expr(:call, lv(:vsub), loop.stophint + 1, loop.startsym)
        end
    elseif loop.startexact
        if isone(loop.starthint)
            loop.stopsym
        else
            Expr(:call, lv(:vsub), loop.stopsym, loop.starthint - 1)
        end
    else
        Expr(:call, lv(:vsub), loop.stopsym, Expr(:call, lv(:staticm1), loop.startsym))
    end
end
# use_expect() = false
use_expect() = true
function looplengthexpr(loop, n)
    le = looplengthexpr(loop)
    # if false && use_expect() && isone(n) && !isstaticloop(loop)
    #     le = expect(le)
    #     push!(le.args, Expr(:call, Expr(:curly, :Val, length(loop))))
    # end
    le
end

# load/compute/store × isunrolled × istiled × pre/post loop × Loop number
struct LoopOrder <: AbstractArray{Vector{Operation},5}
    oporder::Vector{Vector{Operation}}
    loopnames::Vector{Symbol}
    bestorder::Vector{Symbol}
end
# function LoopOrder(N::Int)
#     LoopOrder(
#         [ Operation[] for _ ∈ 1:8N ],
#         Vector{Symbol}(undef, N), Vector{Symbol}(undef, N)
#     )
# end
LoopOrder() = LoopOrder(Vector{Operation}[],Symbol[],Symbol[])
Base.empty!(lo::LoopOrder) = foreach(empty!, lo.oporder)
function Base.resize!(lo::LoopOrder, N::Int)
    Nold = length(lo.loopnames)
    resize!(lo.oporder, 8N)
    for n ∈ 8Nold+1:8N
        lo.oporder[n] = Operation[]
    end
    resize!(lo.loopnames, N)
    resize!(lo.bestorder, N)
    lo
end
Base.size(lo::LoopOrder) = (2,2,2,length(lo.loopnames))
Base.@propagate_inbounds Base.getindex(lo::LoopOrder, i::Int) = lo.oporder[i]
Base.@propagate_inbounds Base.getindex(lo::LoopOrder, i...) = lo.oporder[LinearIndices(size(lo))[i...]]

@enum NumberType::Int8 HardInt HardFloat IntOrFloat INVALID

struct LoopStartStopManager
    terminators::Vector{Int}
    incrementedptrs::Vector{Vector{ArrayReferenceMeta}}
    uniquearrayrefs::Vector{ArrayReferenceMeta}
end
# Must make it easy to iterate
# outer_reductions is a vector of indices (within operation vectors) of the reduction operation, eg the vmuladd op in a dot product
# O(N) search is faster at small sizes
struct LoopSet
    loopsymbols::Vector{Symbol}
    loopsymbol_offsets::Vector{Int}  # symbol loopsymbols[i] corresponds to loops[lso[i]+1:lso[i+1]] (CartesianIndex handling)
    loops::Vector{Loop}
    opdict::Dict{Symbol,Operation}
    operations::Vector{Operation} # Split them to make it easier to iterate over just a subset
    operation_offsets::Vector{Int}
    outer_reductions::Vector{Int} # IDs of reduction operations that need to be reduced at end.
    loop_order::LoopOrder
    preamble::Expr
    prepreamble::Expr # performs extractions that must be performed first, and don't need further registering
    preamble_symsym::Vector{Tuple{Int,Symbol}}
    preamble_symint::Vector{Tuple{Int,Int}}
    preamble_symfloat::Vector{Tuple{Int,Float64}}
    preamble_zeros::Vector{Tuple{Int,NumberType}}
    preamble_funcofeltypes::Vector{Tuple{Int,Symbol}}
    includedarrays::Vector{Symbol}
    includedactualarrays::Vector{Symbol}
    syms_aliasing_refs::Vector{Symbol}
    refs_aliasing_syms::Vector{ArrayReferenceMeta}
    cost_vec::Matrix{Float64}
    reg_pres::Matrix{Float64}
    included_vars::Vector{Bool}
    place_after_loop::Vector{Bool}
    unrollspecification::Base.RefValue{UnrollSpecification}
    loadelimination::Base.RefValue{Bool}
    lssm::Base.RefValue{LoopStartStopManager}
    vector_width::Base.RefValue{Int}
    isbroadcast::Base.RefValue{Bool}
    mod::Symbol
end

function cost_vec_buf(ls::LoopSet)
    cv = @view(ls.cost_vec[:,2])
    @inbounds for i ∈ 1:4
        cv[i] = 0.0
    end
    cv
end
function reg_pres_buf(ls::LoopSet)
    ps = @view(ls.reg_pres[:,2])
    @inbounds for i ∈ 1:5
        ps[i] = 0
    end
    ps
end
function save_tilecost!(ls::LoopSet)
    @inbounds for i ∈ 1:4
        ls.cost_vec[i,1] = ls.cost_vec[i,2]
        ls.reg_pres[i,1] = ls.reg_pres[i,2]
    end
    ls.reg_pres[5,1] = ls.reg_pres[5,2]
end

pushprepreamble!(ls::LoopSet, ex) = push!(ls.prepreamble.args, ex)
function pushpreamble!(ls::LoopSet, op::Operation, v::Symbol)
    if v !== mangledvar(op)
        push!(ls.preamble_symsym, (identifier(op),v))
    end
    nothing
end
function pushpreamble!(ls::LoopSet, op::Operation, v::Number)
    typ = v isa Integer ? HardInt : HardFloat
    id = identifier(op)
    if iszero(v)
        push!(ls.preamble_zeros, (id, typ))
    elseif isone(v)
        push!(ls.preamble_funcofeltypes, (id, :one))
    elseif v isa Integer
        push!(ls.preamble_symint, (id, convert(Int,v)))
    else
        push!(ls.preamble_symfloat, (id, convert(Float64,v)))
    end
end
pushpreamble!(ls::LoopSet, ex::Expr) = push!(ls.preamble.args, ex)
# function pushpreamble!(ls::LoopSet, op::Operation, RHS::Expr)
#     c = gensym(:licmconst)
#     if RHS.head === :call && first(RHS.args) === :zero
#         push!(ls.preamble_zeros, (identifier(op), IntOrFloat))
#     elseif RHS.head === :call && first(RHS.args) === :one
#         push!(ls.preamble_funcofeltypes, (identifier(op), :one))
#     else
#         pushpreamble!(ls, Expr(:(=), c, RHS))
#         pushpreamble!(ls, op, c)
#     end
#     nothing
# end
function zerotype(ls::LoopSet, op::Operation)
    opid = identifier(op)
    for (id,typ) ∈ ls.preamble_zeros
        id == opid && return typ
    end
    INVALID
end

includesarray(ls::LoopSet, array::Symbol) = array ∈ ls.includedarrays

function LoopSet(mod::Symbol)
    LoopSet(
        Symbol[], [0], Loop[],
        Dict{Symbol,Operation}(),
        Operation[], [0],
        Int[],
        LoopOrder(),
        Expr(:block),Expr(:block),
        Tuple{Int,Symbol}[],
        Tuple{Int,Int}[],
        Tuple{Int,Float64}[],
        Int[],Int[],
        Symbol[], Symbol[], Symbol[],
        ArrayReferenceMeta[],
        Matrix{Float64}(undef, 4, 2),
        Matrix{Float64}(undef, 5, 2),
        Bool[], Bool[], Ref{UnrollSpecification}(),
        Ref(false), Ref{LoopStartStopManager}(),
        Ref(0), Ref(false), mod
    )
end

cacheunrolled!(ls::LoopSet, u₁loop, u₂loop, vectorized) = foreach(op -> setunrolled!(op, u₁loop, u₂loop, vectorized), operations(ls))

num_loops(ls::LoopSet) = length(ls.loops)
function oporder(ls::LoopSet)
    N = length(ls.loop_order.loopnames)
    reshape(ls.loop_order.oporder, (2,2,2,N))
end
names(ls::LoopSet) = ls.loop_order.loopnames
reversenames(ls::LoopSet) = ls.loop_order.bestorder
function getloopid_or_nothing(ls::LoopSet, s::Symbol)
    for (loopnum,sym) ∈ enumerate(ls.loopsymbols)
        s === sym && return loopnum
    end
end
getloopid(ls::LoopSet, s::Symbol) = getloopid_or_nothing(ls, s)::Int
getloop(ls::LoopSet, s::Symbol) = ls.loops[getloopid(ls, s)]
# getloop(ls::LoopSet, i::Integer) = ls.loops[i]
getloopsym(ls::LoopSet, i::Integer) = ls.loopsymbols[i]
Base.length(ls::LoopSet, s::Symbol) = length(getloop(ls, s))

# isstaticloop(ls::LoopSet, s::Symbol) = isstaticloop(getloop(ls,s))
# looprangehint(ls::LoopSet, s::Symbol) = length(getloop(ls, s))
# looprangesym(ls::LoopSet, s::Symbol) = getloop(ls, s).rangesym

"""
getop only works while construction a LoopSet object. You cannot use it while lowering.
"""
getop(ls::LoopSet, var::Number, elementbytes) = add_constant!(ls, var, elementbytes)
function getop(ls::LoopSet, var::Symbol, elementbytes::Int)
    get!(ls.opdict, var) do
        add_constant!(ls, var, elementbytes)
    end
end
function getop(ls::LoopSet, var::Symbol, deps, elementbytes::Int)
    get!(ls.opdict, var) do
        add_constant!(ls, var, deps, gensym(:constant), elementbytes)
    end
end
getop(ls::LoopSet, i::Int) = ls.operations[i]

# """
# Returns an operation with the same name as `s`.
# """
# function getoperation(ls::LoopSet, s::Symbol)
#     for op ∈ Iterators.Reverse(operations(ls))
#         name(op) === s && return op
#     end
#     throw("Symbol $s not found among operations(ls).")
# end

function Operation(
    ls::LoopSet, variable, elementbytes, instruction,
    node_type, dependencies, reduced_deps, parents, ref = NOTAREFERENCE
)
    Operation(
        length(operations(ls)), variable, elementbytes, instruction,
        node_type, dependencies, reduced_deps, parents, ref
    )
end
function Operation(ls::LoopSet, variable, elementbytes, instr, optype, mpref::ArrayReferenceMetaPosition)
    Operation(length(operations(ls)), variable, elementbytes, instr, optype, mpref)
end

operations(ls::LoopSet) = ls.operations
function pushop!(ls::LoopSet, op::Operation, var::Symbol = name(op))
    for opp ∈ operations(ls)
        if matches(op, opp)
            ls.opdict[var] = opp
            return opp
        end
    end
    push!(ls.operations, op)
    ls.opdict[var] = op
    op
end
function add_block!(ls::LoopSet, ex::Expr, elementbytes::Int, position::Int)
    for x ∈ ex.args
        x isa Expr || continue # be that general?
        x.head === :inbounds && continue
        push!(ls, x, elementbytes, position)
    end
end
function maybestatic!(expr::Expr)
    if expr.head === :call
        f = first(expr.args)
        if f === :length
            expr.args[1] = lv(:maybestaticlength)
        elseif f === :size && length(expr.args) == 3
            i = expr.args[3]
            if i isa Integer
                expr.args[1] = lv(:maybestaticsize)
                expr.args[3] = Expr(:call, Expr(:curly, :Val, convert(Int, i)))
            end
        end
    end
    expr
end
function add_loop_bound!(ls::LoopSet, itersym::Symbol, bound, upper::Bool = true)
    (bound isa Symbol && upper) && return bound
    bound isa Expr && maybestatic!(bound)
    N = gensym(string(itersym) * (upper ? "_loop_upper_bound" : "_loop_lower_bound"))
    pushprepreamble!(ls, Expr(:(=), N, bound))
    N
end

"""
This function creates a loop, while switching from 1 to 0 based indices
"""
function register_single_loop!(ls::LoopSet, looprange::Expr)
    itersym = (looprange.args[1])::Symbol
    r = looprange.args[2]
    if isexpr(r, :call)
        f = first(r.args)
        loop::Loop = if f === :(:)
            lower = r.args[2]
            upper = r.args[3]
            lii::Bool = lower isa Integer
            liiv::Int = lii ? convert(Int, lower) : 1
            uii::Bool = upper isa Integer
            if lii & uii # both are integers
                Loop(itersym, liiv, convert(Int, upper))
            elseif lii # only lower bound is an integer
                if upper isa Symbol
                    Loop(itersym, liiv, upper)
                elseif upper isa Expr
                    Loop(itersym, liiv, add_loop_bound!(ls, itersym, upper, true))
                else
                    Loop(itersym, liiv, add_loop_bound!(ls, itersym, upper, true))
                end
            elseif uii # only upper bound is an integer
                uiiv = convert(Int, upper)
                Loop(itersym, add_loop_bound!(ls, itersym, lower, false), uiiv)
            else # neither are integers
                L = add_loop_bound!(ls, itersym, lower, false)
                U = add_loop_bound!(ls, itersym, upper, true)
                Loop(itersym, L, U)
            end
        elseif f === :OneTo || isscopedname(f, :Base, :OneTo)
            otN = r.args[2]
            if otN isa Integer
                Loop(itersym, 1, otN)
            else
                otN isa Expr && maybestatic!(otN)
                N = gensym("loop" * string(itersym))
                pushprepreamble!(ls, Expr(:(=), N, otN))
                Loop(itersym, 1, N)
            end
        else
            N = gensym("loop" * string(itersym))
            pushprepreamble!(ls, Expr(:(=), N, Expr(:call, lv(:maybestaticrange), r)))
            L = add_loop_bound!(ls, itersym, Expr(:call, lv(:maybestaticfirst), N), false)
            U = add_loop_bound!(ls, itersym, Expr(:call, lv(:maybestaticlast), N), true)
            Loop(itersym, L, U)
        end
    elseif isa(r, Symbol)
        # Treat similar to `eachindex`
        N = gensym("loop" * string(itersym))
        pushprepreamble!(ls, Expr(:(=), N, Expr(:call, lv(:maybestaticrange), r)))
        L = add_loop_bound!(ls, itersym, Expr(:call, lv(:maybestaticfirst), N), false)
        U = add_loop_bound!(ls, itersym, Expr(:call, lv(:maybestaticlast), N), true)
        loop = Loop(itersym, L, U)
    else
        throw("Unrecognized loop range type: $r.")
    end
    add_loop!(ls, loop, itersym)
    nothing
end
function register_loop!(ls::LoopSet, looprange::Expr)
    if looprange.head === :block # multiple loops
        for lr ∈ looprange.args
            register_single_loop!(ls, lr)
        end
    else
        @assert looprange.head === :(=)
        register_single_loop!(ls, looprange)
    end
end
function add_loop!(ls::LoopSet, q::Expr, elementbytes::Int)
    register_loop!(ls, q.args[1])
    body = q.args[2]
    position = length(ls.loopsymbols)
    if body.head === :block
        add_block!(ls, body, elementbytes, position)
    else
        push!(ls, q, elementbytes, position)
    end
end
function add_loop!(ls::LoopSet, loop::Loop, itersym::Symbol = loop.itersymbol)
    push!(ls.loopsymbols, itersym)
    push!(ls.loops, loop)
    nothing
end

# function instruction(x)
#     x isa Symbol ? x : last(x.args).value
# end
# instruction(ls::LoopSet, f::Symbol) = instruction!(ls, f)
function instruction!(ls::LoopSet, x::Expr)
    x isa Symbol && return x
    instr = last(x.args).value
    if instr ∉ keys(COST)
        instr = gensym(:f)
        pushpreamble!(ls, Expr(:(=), instr, x))
    end
    Instruction(Symbol(""), instr)
end
instruction!(ls::LoopSet, x::Symbol) = instruction(x)

function maybe_const_compute!(ls::LoopSet, LHS::Symbol, op::Operation, elementbytes::Int, position::Int)
    # return op
    if iscompute(op) && iszero(length(loopdependencies(op)))
        ls.opdict[LHS] = add_constant!(ls, LHS, ls.loopsymbols[1:position], gensym(instruction(op).instr), elementbytes, :numericconstant)
    else
        op
    end
end
function strip_op_linenumber_nodes(q::Expr)
    filtered = filter(x -> !isa(x, LineNumberNode), q.args)
    if VERSION ≥ v"1.4"
        only(filtered)
    else
        @assert isone(length(filtered))
        first(filtered)
    end
end

function add_operation!(ls::LoopSet, LHS::Symbol, RHS::Symbol, elementbytes::Int, position::Int)
    add_constant!(ls, RHS, ls.loopsymbols[1:position], LHS, elementbytes)
end
function add_operation!(
    ls::LoopSet, LHS::Symbol, RHS::Expr, elementbytes::Int, position::Int
)
    if RHS.head === :ref
        add_load_ref!(ls, LHS, RHS, elementbytes)
    elseif RHS.head === :call
        f = first(RHS.args)
        if f === :getindex
            add_load_getindex!(ls, LHS, RHS, elementbytes)
        elseif f === :zero || f === :one
            c = gensym(f)
            op = add_constant!(ls, c, ls.loopsymbols[1:position], LHS, elementbytes, :numericconstant)
            if f === :zero
                push!(ls.preamble_zeros, (identifier(op), IntOrFloat))
            else
                push!(ls.preamble_funcofeltypes, (identifier(op), :one))
            end
            op
        else
            # maybe_const_compute!(ls, add_compute!(ls, LHS, RHS, elementbytes, position), elementbytes, position)
            add_compute!(ls, LHS, RHS, elementbytes, position)
        end
    elseif RHS.head === :if
        add_if!(ls, LHS, RHS, elementbytes, position)
    elseif RHS.head === :block
        add_operation!(ls, LHS, strip_op_linenumber_nodes(RHS), elementbytes, position)
    else
        println(RHS)
        throw("Expression not recognized.")
    end
end
add_operation!(ls::LoopSet, RHS::Expr, elementbytes::Int, position::Int) = add_operation!(ls, gensym(:LHS), RHS, elementbytes, position)
function add_operation!(
    ls::LoopSet, LHS_sym::Symbol, RHS::Expr, LHS_ref::ArrayReferenceMetaPosition, elementbytes::Int, position::Int
)
    if RHS.head === :ref# || (RHS.head === :call && first(RHS.args) === :getindex)
        array, rawindices = ref_from_expr!(ls, RHS)
        RHS_ref = array_reference_meta!(ls, array, rawindices, elementbytes, gensym(LHS_sym))
        op = add_load!(ls, RHS_ref, elementbytes)
        iop = add_compute!(ls, LHS_sym, :identity, [op], elementbytes)
        # pushfirst!(LHS_ref.parents, iop)
    elseif RHS.head === :call
        f = first(RHS.args)
        if f === :getindex
            add_load!(ls, LHS_sym, LHS_ref, elementbytes)
        elseif f === :zero || f === :one
            c = gensym(f)
            op = add_constant!(ls, c, ls.loopsymbols[1:position], LHS_sym, elementbytes, :numericconstant)
            # op = add_constant!(ls, c, Symbol[], LHS_sym, elementbytes, :numericconstant)
            if f === :zero
                push!(ls.preamble_zeros, (identifier(op), IntOrFloat))
            else
                push!(ls.preamble_funcofeltypes, (identifier(op), :one))
            end
            op
        else
            add_compute!(ls, LHS_sym, RHS, elementbytes, position, LHS_ref)
        end
    elseif RHS.head === :if
        add_if!(ls, LHS_sym, RHS, elementbytes, position, LHS_ref)
    elseif RHS.head === :block
        add_operation!(ls, LHS, strip_op_linenumber_nodes(RHS), elementbytes, position)
    else
        println(RHS)
        throw("Expression not recognized.")
    end
end

function prepare_rhs_for_storage!(ls::LoopSet, RHS::Symbol, array, rawindices, elementbytes::Int, position::Int)
    add_store!(ls, RHS, array, rawindices, elementbytes)
end
function prepare_rhs_for_storage!(ls::LoopSet, RHS::Expr, array, rawindices, elementbytes::Int, position::Int)
    mpref = array_reference_meta!(ls, array, rawindices, elementbytes)
    cachedparents = copy(mpref.parents)
    ref = mpref.mref.ref
    # id = findfirst(r -> r == ref, ls.refs_aliasing_syms)
    # lrhs = id === nothing ? gensym(:RHS) : ls.syms_aliasing_refs[id]
    lrhs = gensym(:RHS)
    mpref.varname = lrhs
    add_operation!(ls, lrhs, RHS, mpref, elementbytes, position)
    mpref.parents = cachedparents
    add_store!(ls, mpref, elementbytes)
end

function Base.push!(ls::LoopSet, ex::Expr, elementbytes::Int, position::Int)
    if ex.head === :call
        finex = first(ex.args)::Symbol
        if finex === :setindex!
            array, rawindices = ref_from_setindex!(ls, ex)
            prepare_rhs_for_storage!(ls, ex.args[3], array, rawindices, elementbytes, position)
        else
            throw("Function $finex not recognized.")
        end
    elseif ex.head === :(=)
        LHS = ex.args[1]
        RHS = ex.args[2]
        if LHS isa Symbol
            if RHS isa Expr
                maybe_const_compute!(ls, LHS, add_operation!(ls, LHS, RHS, elementbytes, position), elementbytes, position)
            else
                add_constant!(ls, RHS, ls.loopsymbols[1:position], LHS, elementbytes)
            end
        elseif LHS isa Expr
            if LHS.head === :ref
                if RHS isa Symbol
                    add_store_ref!(ls, RHS, LHS, elementbytes)
                elseif RHS isa Expr
                    # need to check if LHS appears in RHS
                    # assign RHS to lrhs
                    array, rawindices = ref_from_expr!(ls, LHS)
                    prepare_rhs_for_storage!(ls, RHS, array, rawindices, elementbytes, position)
                else
                    add_store_ref!(ls, RHS, LHS, elementbytes)
                end
            elseif LHS.head === :tuple
                @assert length(LHS.args) ≤ 9 "Functions returning more than 9 values aren't currently supported."
                lhstemp = gensym(:lhstuple)
                vparents = Operation[maybe_const_compute!(ls, lhstemp, add_operation!(ls, lhstemp, RHS, elementbytes, position), elementbytes, position)]
                for i ∈ eachindex(LHS.args)
                    f = (:first,:second,:third,:fourth,:fifth,:sixth,:seventh,:eighth,:ninth)[i]
                    lhsi = LHS.args[i]
                    if lhsi isa Symbol
                        add_compute!(ls, lhsi, f, vparents, elementbytes)
                    elseif lhsi isa Expr && lhsi.head === :ref
                        tempunpacksym = gensym(:tempunpack)
                        add_compute!(ls, tempunpacksym, f, vparents, elementbytes)
                        add_store_ref!(ls, tempunpacksym, lhsi, elementbytes)
                    else
                        println(lhsi)
                        throw("Unpacking the above expression in the left hand side was not understood/supported.")
                    end
                end
                first(vparents)
            else
                println(LHS)
                throw("LHS not understood; only `:ref`s and `:tuple`s are currently supported.")
            end
        else
            println(LHS)
            throw("LHS not understood.")
        end
    elseif ex.head === :block
        add_block!(ls, ex, elementbytes, position)
    elseif ex.head === :for
        add_loop!(ls, ex, elementbytes)
    elseif ex.head === :&&
        add_andblock!(ls, ex, elementbytes, position)
    elseif ex.head === :||
        add_orblock!(ls, ex, elementbytes, position)
    elseif ex.head === :local # Handle locals introduced by `@inbounds`; using `local` with `@avx` is not recomended (nor is `@inbounds`; which applies automatically regardless)
        @assert length(ex.args) == 1 # TODO replace assert + first with "only" once support for Julia < 1.4 is dropped
        localbody = first(ex.args)
        @assert localbody.head === :(=)
        @assert length(localbody.args) == 2
        LHS = (localbody.args[1])::Symbol
        RHS = push!(ls, (localbody.args[2]), elementbytes, position)
        if isstore(RHS)
            RHS
        else
            add_compute!(ls, LHS, :identity, [RHS], elementbytes)
        end
    else
        println(ex)
        throw("Don't know how to handle expression.")
    end
end

function UnrollSpecification(ls::LoopSet, u₁loop::Symbol, u₂loop::Symbol, vectorized::Symbol, u₁, u₂)
    order = names(ls)
    nu₁ = findfirst(isequal(u₁loop), order)::Int
    nu₂ = u₂ == -1 ? nu₁ : findfirst(isequal(u₂loop), order)::Int
    nv = findfirst(isequal(vectorized), order)::Int
    UnrollSpecification(nu₁, nu₂, nv, u₁, u₂)
end

"""
  looplengthprod(ls::LoopSet)

Convert to `Float64` for the sake of non-64 bit platforms.
"""
looplengthprod(ls::LoopSet) = prod(Float64 ∘ length, ls.loops)

function looplength(ls::LoopSet, s::Symbol)
    # search_tree(parents(operations(ls)[i]), name(op)) && return true
    id = getloopid_or_nothing(ls, s)
    if isnothing(id)
        l = 0.0
        # TODO: we could double count a loop.
        for op ∈ operations(ls)
            name(op) === s || continue
            for opp ∈ parents(op)
                if isloopvalue(opp)
                    oppname = first(loopdependencies(opp))
                    l += looplength(ls, oppname)
                elseif iscompute(opp)
                    oppname = name(opp)
                    l += looplength(ls, oppname)
                    # TODO elseif isconstant(opp)
                end
            end
            l += 1 - length(parents(op))
        end
        l
    else
        Float64(length(ls.loops[id]))
    end
end

# function getunrolled(ls::LoopSet)
#     order = names(ls)
#     us = ls.unrollspecification[]
#     @unpack u₁loopnum, u₂loopnum = us
#     order[u₁loopnum], order[u₂loopnum]
# end


text:
module SparseJacobian

using DelimitedFiles
import ForwardDiff
import FLOWFarm
using LinearAlgebra
using SparseArrays, SparseDiffTools
using FLOWMath: gradient
using Colors
using Traceur

## Derivative calculation
function sparseAEPJacobian!(x::Vector{Float64},params,deriv)
    deriv.Jacobian .= 0
    for i = 1:length(params.windresource.wind_directions)
        # calculate pattern
        calculateSparsityPattern!(x,deriv.pattern,params,i,deriv.sorted_index)

        if !isequal(deriv.pattern,deriv.patterns[:,:,i])
            deriv.patterns[:,:,i] .= deriv.pattern
            allocateNewCache!(x,deriv,i)
        end

        # calculate state Jacobian
        calculateSparseJacobian!(x,deriv,i)

        # sum jacobian columns and add to other states
        deriv.Jacobian .= deriv.Jacobian .+ transpose(sum(deriv.jac,dims = 1) .* params.windresource.wind_probabilities[i])
    end
    deriv.Jacobian .= deriv.Jacobian .* params.obj_scale .* 365.25 .* 24
end

function calculateSparseJacobian!(x::Vector{Float64},deriv,currentState::Int64)
    deriv.jac .= sparse(deriv.patterns[:,:,currentState])
    if iszero(deriv.jac)
        return nothing
    end

    forwarddiff_color_jacobian!(deriv.jac,
                            deriv.functions[currentState],
                            x,
                            deriv.caches[currentState])
end

function calculateSparsityPattern!(x::Vector{Float64},pattern::Matrix{Float64},params,currentState::Int64,sorted_index::Vector{Int64})
    n = Int(length(x)/2)
    rot_x, rot_y = ff.rotate_to_wind_direction(x[1:n],x[n+1:end],params.windresource.wind_directions[currentState])::Tuple{Vector{Float64}, Vector{Float64}}
    rot_x .= rot_x ./ params.rotor_diameter
    rot_y .= rot_y ./ params.rotor_diameter
    pattern .= 0
    sorted_index .= sortperm(rot_x)::Vector{Int64}
    d = 3.4 #3.4 #4.5

    # determine if current turbine is affected by the other turbines
    for i = 1:Int(length(x)/2)
        xCurrent = rot_x[i]
        yCurrent = rot_y[i]
        for j = 1:Int(length(x)/2)
            jt = sorted_index[j]
            xOther = rot_x[jt]
            xdif = xCurrent - xOther
            if i == jt
                pattern[i,jt+n] = 1.0
                pattern[i,jt] = 1.0
                continue
            elseif xdif < 0
                break
            end

            yOther = rot_y[jt]
            if yOther < yCurrent
                yOther = 2*yCurrent - yOther
            end

            ydif = yOther - yCurrent

            if ydif < d
                pattern[i,jt+n] = 1.0
                pattern[i,jt] = 1.0
                continue
            elseif ydif > (d + xdif*tand(15)) #15
                continue
            elseif ydif == 0
                pattern[i,jt] = 1.0
            else
                pattern[i,jt+n] = 1.0
                pattern[i,jt] = 1.0
            end
        end
    end
end

mutable struct deriv_struct{MF, SMFI, VF, AF3, AA, VI}
    pattern::MF
    jac::SMFI
    Jacobian::VF
    patterns::AF3
    caches::AA
    sorted_index::VI
    functions::AA
end

function allocateNewCache!(x::Vector{Float64},deriv,currentState::Int64)
    # global newCaches = newCaches+1
    deriv.jac .= dropzeros(sparse(deriv.patterns[:,:,currentState]))
    if iszero(deriv.jac)
        colors = collect(Int64,1:length(x))
    else
        colors = matrix_colors(deriv.jac)
    end
    deriv.caches[currentState] = ForwardColorJacCache(deriv.functions[currentState],x,
                              dx = similar(x[1:Int(length(x)/2)]),
                              colorvec=colors,
                              sparsity = deriv.jac)
end

function allocateContainers(x::Vector{Float64},params)
    n = Int(length(x)/2)
    numStates = length(params.windresource.wind_directions)
    pattern = zeros(n,n*2)
    jac = spzeros(n,n*2)
    Jacobian = zeros(n*2)
    patterns = zeros(n,n*2,numStates)
    caches = Array{Any,1}(undef,numStates)
    functions = Array{Any,1}(undef,numStates)
    colors = zeros(Int,n*2)
    index = zeros(Int,n)

    for i = 1:numStates
        calculateSparsityPattern!(x,pattern,params,i,index)
        patterns[:,:,i] .= pattern
        jac .= dropzeros(sparse(patterns[:,:,i]))
        if iszero(jac)
            colors = collect(Int64,1:length(colors))
        else
            colors = matrix_colors(jac)
        end
        f = (wt_power,x) -> p_wrapper!(wt_power,x,params,i)
        functions[i] = f
        caches[i] = ForwardColorJacCache(f,x,
                              dx = similar(x[1:n]),
                              colorvec=colors,
                              sparsity = jac)
    end
    deriv = deriv_struct(pattern,jac,Jacobian,patterns,caches,index,functions)
    return deriv
end

## Optimization calculations
function opt!(g,df,dg,x,params)
    # calculate spacing constraint value and jacobian
    spacing_con = spacing_wrapper(x)

    # calculate boundary constraint and jacobian
    boundary_con = boundary_wrapper(x)

    # combine constaint values and jacobians into overall constaint value and jacobian arrays
    c = [spacing_con; boundary_con]
    g[:] = c[:]

    # calculate the objective function and jacobian (negative sign in order to maximize AEP)
    AEP = -aep_wrapper(x)[1]

    df[:] = -ForwardDiff.jacobian(aep_wrapper,x)
    ds_dx = ForwardDiff.jacobian(spacing_wrapper, x)
    db_dx = ForwardDiff.jacobian(boundary_wrapper, x)
    dcdx = [ds_dx; db_dx]
    dg[:] = dcdx[:]

    return AEP
end

function opt_fast!(g,df,dg,x,params,deriv)
    # calculate spacing constraint value and jacobian
    spacing_con = spacing_wrapper(x)

    # calculate boundary constraint and jacobian
    boundary_con = boundary_wrapper(x)

    # combine constaint values and jacobians into overall constaint value and jacobian arrays
    c = [spacing_con; boundary_con]
    g[:] = c[:]

    # calculate the objective function and jacobian (negative sign in order to maximize AEP)
    AEP = -aep_wrapper(x)[1]

    sparseAEPJacobian!(x,params,deriv)

    df[:] = -deriv.Jacobian
    ds_dx = ForwardDiff.jacobian(spacing_wrapper, x)
    db_dx = ForwardDiff.jacobian(boundary_wrapper, x)
    dcdx = [ds_dx; db_dx]
    dg[:] = dcdx[:]

    return AEP
end

## Farm setup
struct params_struct{MS, AF, F, AI, ACTM, WR, APM, SPAF, SPAF2D}
    model_set::MS
    rotor_points_y::AF
    rotor_points_z::AF
    turbine_z::AF
    rotor_diameter::AF
    boundary_center::AF
    boundary_radius::F
    obj_scale::F
    hub_height::AF
    turbine_yaw::AF
    ct_models::ACTM
    generator_efficiency::AF
    cut_in_speed::AF
    cut_out_speed::AF
    rated_speed::AF
    rated_power::AF
    windresource::WR
    power_models::APM
    i::AI
    dx::SPAF
    jac::SPAF2D
    color::SPAF
end

# set up boundary constraint wrapper function
function boundary_wrapper(x, params)

    # get number of turbines
    nturbines = Int(length(x)/2)

    # extract x and y locations of turbines from design variables vector
    turbine_x = x[1:nturbines]
    turbine_y = x[nturbines+1:end]

    # get and return boundary distances
    return ff.circle_boundary(params.boundary_center, params.boundary_radius, turbine_x, turbine_y) .* 1E-7
end

# set up spacing constraint wrapper function
function spacing_wrapper(x, params)

    # get number of turbines
    nturbines = Int(length(x)/2)

    # extract x and y locations of turbines from design variables vector
    turbine_x = x[1:nturbines]
    turbine_y = x[nturbines+1:end]

    # get and return spacing distances
    return (2.0*params.rotor_diameter[1] .- ff.turbine_spacing(turbine_x,turbine_y)) .* 1E-4
end

# set up objective wrapper function
function aep_wrapper(x, params)

    # get number of turbines
    nturbines = Int(length(x)/2)

    # extract x and y locations of turbines from design variables vector
    turbine_x = x[1:nturbines]
    turbine_y = x[nturbines+1:end]

    # calculate AEP
    obj_scale = params.obj_scale
    AEP = obj_scale*ff.calculate_aep(turbine_x, turbine_y, params.turbine_z, params.rotor_diameter,
                params.hub_height, params.turbine_yaw, params.ct_models, params.generator_efficiency, params.cut_in_speed,
                params.cut_out_speed, params.rated_speed, params.rated_power, params.windresource, params.power_models, params.model_set,
                rotor_sample_points_y=params.rotor_points_y,rotor_sample_points_z=params.rotor_points_z)

    # return the objective as an array
    return [AEP]
end

function p_wrapper(x,params,state)
    n = Int(length(x)/2)
    turbine_x = x[1:n]
    turbine_y = x[n+1:end]

    rot_x, rot_y = ff.rotate_to_wind_direction(turbine_x, turbine_y, params.windresource.wind_directions[state])

    sorted_turbine_index = sortperm(rot_x)
    turbine_velocities = ff.turbine_velocities_one_direction(rot_x, rot_y, params.turbine_z, params.rotor_diameter, params.hub_height, params.turbine_yaw,
                            sorted_turbine_index, params.ct_models, params.rotor_points_y, params.rotor_points_z, params.windresource,
                            params.model_set, wind_farm_state_id=state, velocity_only=true)

    wt_power = ff.turbine_powers_one_direction(params.generator_efficiency, params.cut_in_speed, params.cut_out_speed, params.rated_speed,
                            params.rated_power, params.rotor_diameter, turbine_velocities, params.turbine_yaw, params.windresource.air_density, params.power_models)
end

function p_wrapper_scaled(x,params,state)
    n = Int(length(x)/2)
    turbine_x = x[1:n]
    turbine_y = x[n+1:end]

    rot_x, rot_y = ff.rotate_to_wind_direction(turbine_x, turbine_y, params.windresource.wind_directions[state])

    sorted_turbine_index = sortperm(rot_x)
    turbine_velocities = ff.turbine_velocities_one_direction(rot_x, rot_y, params.turbine_z, params.rotor_diameter, params.hub_height, params.turbine_yaw,
                            sorted_turbine_index, params.ct_models, params.rotor_points_y, params.rotor_points_z, params.windresource,
                            params.model_set, wind_farm_state_id=state, velocity_only=true)

    wt_power = ff.turbine_powers_one_direction(params.generator_efficiency, params.cut_in_speed, params.cut_out_speed, params.rated_speed,
                            params.rated_power, params.rotor_diameter, turbine_velocities, params.turbine_yaw, params.windresource.air_density, params.power_models)

    wt_power = wt_power * params.obj_scale .* 365.25 .* 24
end

function p_wrapper!(wt_power,x,params,currentState)
    n = Int(length(x)/2)
    turbine_x = x[1:n]
    turbine_y = x[n+1:end]

    rot_x, rot_y = ff.rotate_to_wind_direction(turbine_x, turbine_y, params.windresource.wind_directions[currentState])
    sorted_turbine_index = sortperm(rot_x)
    turbine_velocities = ff.turbine_velocities_one_direction(rot_x, rot_y, params.turbine_z, params.rotor_diameter, params.hub_height, params.turbine_yaw,
                            sorted_turbine_index, params.ct_models, params.rotor_points_y, params.rotor_points_z, params.windresource,
                            params.model_set, wind_farm_state_id=currentState, velocity_only=true)

    wt_power .= ff.turbine_powers_one_direction(params.generator_efficiency, params.cut_in_speed, params.cut_out_speed, params.rated_speed,
                            params.rated_power, params.rotor_diameter, turbine_velocities, params.turbine_yaw, params.windresource.air_density, params.power_models)

    nothing
end

## dirs is number of wind directions 12, 72, 1
## num is number of turbines 38 is standard
## angle is angle of the farm, 0 is default
## returns turbine_x,turbine_y,params
function loadRoundFarm(dirs,num,angle,radius=1225.8227848101264)
    cd("/Users/benjaminvarela/.julia/dev/FLOWFarm/test")

    # scale objective to be between 0.1 and 1
    if num >= 19
        obj_scale = 1E-18
    elseif num == 7
        obj_scale = 1E-17
    elseif num == 1
        obj_scale = 1.0
    else
        obj_scale = 1E-13
    end

    # set initial turbine x and y locations
    diam = 80.0
    data = readdlm("inputfiles/layout_38turb_round.txt",  ' ', skipstart=1)
    turbine_x = data[:, 1].*diam
    nturbines = length(turbine_x)
    nturbines = num
    turbine_y = data[:, 2].*diam
    turbine_z = zeros(nturbines)

    turbine_x = turbine_x .- turbine_x[1]
    turbine_y = turbine_y .- turbine_y[1]

    angle = 270 + angle

    turbine_x, turbine_y = ff.rotate_to_wind_direction(turbine_x, turbine_y, angle*pi/180)

    turbine_z = zeros(nturbines)

    # set turbine base heights
    turbine_z = zeros(nturbines) .+ 0.0

    # set turbine yaw values
    turbine_yaw = zeros(nturbines)

    # set turbine design parameters
    rotor_diameter = zeros(nturbines) .+ diam # m
    hub_height = zeros(nturbines) .+ 70.0   # m
    cut_in_speed = zeros(nturbines) .+ 0.  # m/s 4.0
    cut_out_speed = zeros(nturbines) .+ 25.  # m/s
    rated_speed = zeros(nturbines) .+ 16.  # m/s
    rated_power = zeros(nturbines) .+ 2.0E6  # W
    generator_efficiency = zeros(nturbines) .+ 0.944

    # rotor swept area sample points (normalized by rotor radius)
    rotor_points_y = [0.0]
    rotor_points_z = [0.0]

    # set flow parameters
    if dirs == 12
        data = readdlm("inputfiles/windrose_nantucket_12dir.txt",  ' ', skipstart=1)
        winddirections = data[:, 1].*pi/180.0
        windspeeds = data[:,2]
        windprobabilities = data[:, 3]
        nstates = length(windspeeds)
    elseif dirs == 72
        data = readdlm("inputfiles/windrose_amalia_72dirs.txt",  ' ', skipstart=1)
        winddirections = data[:, 1].*pi/180.0
        windspeeds = data[:,2]
        windprobabilities = data[:, 3]
        nstates = length(windspeeds)
    else
        winddirections = [270.0*pi/180]
        windspeeds = [8.0]
        windprobabilities = [1.0]
        nstates = length(windspeeds)
    end

    air_density = 1.1716  # kg/m^3
    ambient_ti = 0.077
    shearexponent = 0.15
    ambient_tis = zeros(nstates) .+ ambient_ti
    measurementheight = zeros(nstates) .+ hub_height[1]

    # load power curve
    powerdata = readdlm("inputfiles/niayifar_vestas_v80_power_curve_observed.txt",  ',', skipstart=1)
    velpoints = powerdata[:,1]
    powerpoints = powerdata[:,2]*1E6

    # initialize power model
    power_model = ff.PowerModelCpPoints(velpoints, powerpoints)
    power_models = Vector{typeof(power_model)}(undef, nturbines)
    for i = 1:nturbines
        power_models[i] = power_model
    end

    # load thrust curve
    ctdata = readdlm("inputfiles/predicted_ct_vestas_v80_niayifar2016.txt",  ',', skipstart=1)
    velpoints = ctdata[:,1]
    ctpoints = ctdata[:,2]

    # initialize thurst model
    ct_model = ff.ThrustModelCtPoints(velpoints, ctpoints)
    ct_models = Vector{typeof(ct_model)}(undef, nturbines)
    for i = 1:nturbines
        ct_models[i] = ct_model
    end

    # initialize wind shear model
    wind_shear_model = ff.PowerLawWindShear(shearexponent)

    # get sorted indecies
    # sorted_turbine_index = sortperm(turbine_x)

    # initialize the wind resource definition
    windresource = ff.DiscretizedWindResource(winddirections, windspeeds, windprobabilities, measurementheight, air_density, ambient_tis, wind_shear_model)

    # set up wake and related models
    wakedeficitmodel = ff.GaussYawVariableSpread()
    wakedeficitmodel.wec_factor[1] = 1.0

    wakedeflectionmodel = ff.GaussYawVariableSpreadDeflection()
    wakecombinationmodel = ff.LinearLocalVelocitySuperposition()
    localtimodel = ff.LocalTIModelNoLocalTI()

    # initialize model set
    model_set = ff.WindFarmModelSet(wakedeficitmodel, wakedeflectionmodel, wakecombinationmodel, localtimodel)

    # initialize sparse partial derivatives array for turbine velocity
    dx = spzeros(nturbines)

    # set wind farm boundary parameters
    boundary_center = [0.0,0.0]
    boundary_radius = radius

    global params = params_struct(model_set, rotor_points_y, rotor_points_z, turbine_z,
    rotor_diameter, boundary_center, boundary_radius, obj_scale, hub_height, turbine_yaw,
    ct_models, generator_efficiency, cut_in_speed, cut_out_speed, rated_speed, rated_power,
    windresource, power_models, [0], dx, [spzeros(nturbines, 2*nturbines)], spzeros(2*nturbines));

    return turbine_x,turbine_y,params
end

## function to be build different sized farms
## num is rings of turbines
## dirs is number of wind directions 12, 72, 1
## angle is angle of the farm, 0 is default
## returns turbine_x,turbine_y,params
function buildRoundFarm(dirs,rings;angle=0)
    x,y = roundFarms(rings)
    xTemp,yTemp,params = loadRoundFarm(dirs,length(x),angle,maximum(x))

    return x,y,params
end

function roundFarms(rings)
    dr = 5.107594936708861
    diameter = 80
    x = zeros(1)
    y = zeros(1)
    if rings == 0
        return x .*diameter,y .*diameter
    end

    for i = 1:rings
        r = dr * i
        angle,numTurbines = calculateAngle(r,dr)
        for i = 1:numTurbines
            x = [x;r*cos(angle*(i-1))]
            y = [y;r*sin(angle*(i-1))]
        end
    end

    return x .*diameter,y .*diameter
end

function calculateAngle(r,dr)
    angle = 0
    dr = floor(dr)
    numTurbines = 25
    first = zeros(2)
    second = zeros(2)
    while true
        angle = 2*pi/numTurbines
        first = [r;0]
        second = [r*cos(angle);r*sin(angle)]
        diff = first .- second
        d = hypot(diff[1],diff[2])
        if d < dr
            numTurbines = numTurbines - 1
        else
            angleTemp = 2*pi/(numTurbines+1)
            first = [r;0]
            second = [r*cos(angleTemp);r*sin(angleTemp)]
            diff = first .- second
            dCheck = hypot(diff[1],diff[2])
            if dCheck < dr
                break
            else
                numTurbines = Int(numTurbines * 2)
            end
        end
    end
    return angle,numTurbines
end

spacing_wrapper(x) = spacing_wrapper(x, params)
aep_wrapper(x) = aep_wrapper(x, params)
boundary_wrapper(x) = boundary_wrapper(x, params)
opt!(g,df,dg,x) = opt!(g,df,dg,x,params)

const ff = FLOWFarm

### single pattern work

function calculateSingleSparsityPattern2!(x::Vector{Float64},pattern::Matrix{Float64},params,currentState::Int64,sorted_index::Vector{Int64})
    n = Int(length(x)/2)
    rot_x, rot_y = ff.rotate_to_wind_direction(x[1:n],x[n+1:end],params.windresource.wind_directions[currentState])::Tuple{Vector{Float64}, Vector{Float64}}
    rot_x .= rot_x ./ params.rotor_diameter
    rot_y .= rot_y ./ params.rotor_diameter
    pattern .= 0
    sorted_index .= sortperm(rot_x)::Vector{Int64}
    d = 8 #72 dirs 19 turbines even 20 doesnt work, 72 63 doesnt work
    # dx = 30

    # determine if current turbine is affected by the other turbines
    for i = 1:Int(length(x)/2)
        xCurrent = rot_x[i]
        yCurrent = rot_y[i]
        for j = 1:Int(length(x)/2)
            jt = sorted_index[j]
            xOther = rot_x[jt]
            xdif = xCurrent - xOther
            if i == jt
                pattern[i,jt+n] = 1.0
                pattern[i,jt] = 1.0
                continue
            elseif xdif < 0
                break
            # elseif xdif > dx
            #     continue
            end

            yOther = rot_y[jt]
            ydif = yOther - yCurrent

            if abs(ydif) <= d
                pattern[i,jt+n] = 1.0
                pattern[i,jt] = 1.0
                continue
            end
        end
        temp = pattern[i,:]
        temp = temp[temp .> 0.0]
        if length(temp) == 2
            pattern[i,:] .= 0
        end
    end
end

function allocateContainersSingle(x::Vector{Float64},params)
    n = Int(length(x)/2)
    numStates = length(params.windresource.wind_directions)
    pattern = zeros(n,n*2)
    jac = spzeros(n,n*2)
    Jacobian = zeros(n*2)
    patterns = zeros(n,n*2,numStates)
    caches = Array{Any,1}(undef,numStates)
    functions = Array{Any,1}(undef,numStates)
    colors = zeros(Int,n*2)
    index = zeros(Int,n)

    for i = 1:numStates
        calculateSingleSparsityPattern2!(x,pattern,params,i,index)
        patterns[:,:,i] .= pattern
        jac .= dropzeros(sparse(patterns[:,:,i]))
        if iszero(jac)
            colors = collect(Int64,1:length(colors))
        else
            colors = matrix_colors(jac)
        end
        f = (wt_power,x) -> p_wrapper!(wt_power,x,params,i)
        functions[i] = f
        caches[i] = ForwardColorJacCache(f,x,
                              dx = similar(x[1:n]),
                              colorvec=colors,
                              sparsity = jac)
    end
    deriv = deriv_struct(pattern,jac,Jacobian,patterns,caches,index,functions)
    return deriv
end

function opt_fast_single!(g,df,dg,x,params,deriv)
    # calculate spacing constraint value and jacobian
    spacing_con = spacing_wrapper(x,params)

    # calculate boundary constraint and jacobian
    boundary_con = boundary_wrapper(x,params)

    # combine constaint values and jacobians into overall constaint value and jacobian arrays
    c = [spacing_con; boundary_con]
    g[:] = c[:]

    # calculate the objective function and jacobian (negative sign in order to maximize AEP)
    AEP = -aep_wrapper(x,params)[1]

    sparseAEPJacobianSinglePattern!(x,params,deriv)

    df[:] = -deriv.Jacobian
    ds_dx = ForwardDiff.jacobian(spacing_wrapper, x)
    db_dx = ForwardDiff.jacobian(boundary_wrapper, x)
    dcdx = [ds_dx; db_dx]
    dg[:] = dcdx[:]

    return AEP
end

function sparseAEPJacobianSinglePattern!(x::Vector{Float64},params,deriv)
    deriv.Jacobian .= 0
    for i = 1:length(params.windresource.wind_directions)
        # calculate state Jacobian
        calculateSparseJacobian!(x,deriv,i)

        # sum jacobian columns and add to other states
        deriv.Jacobian .= deriv.Jacobian .+ transpose(sum(deriv.jac,dims = 1) .* params.windresource.wind_probabilities[i])
    end
    deriv.Jacobian .= deriv.Jacobian .* params.obj_scale .* 365.25 .* 24
end

### update pattern every # function cals
function sparseAEPJacobianUpdateFunctionCalls!(x::Vector{Float64},params,deriv,numCalls)
    deriv.Jacobian .= 0
    maxCalls = 20 #200
    numCalls += 1
    for i = 1:length(params.windresource.wind_directions)
        # calculate pattern
        calculateSingleSparsityPattern2!(x,deriv.pattern,params,i,deriv.sorted_index)

        if numCalls > maxCalls
            deriv.patterns[:,:,i] .= deriv.pattern
            allocateNewCache!(x,deriv,i)
        end

        # calculate state Jacobian
        calculateSparseJacobian!(x,deriv,i)

        # sum jacobian columns and add to other states
        deriv.Jacobian .= deriv.Jacobian .+ transpose(sum(deriv.jac,dims = 1) .* params.windresource.wind_probabilities[i])
    end
    if numCalls > maxCalls
        numCalls = 0
    end
    deriv.Jacobian .= deriv.Jacobian .* params.obj_scale .* 365.25 .* 24
end

function opt_fast_update!(g,df,dg,x,params,deriv,numCalls)
    # calculate spacing constraint value and jacobian
    spacing_con = spacing_wrapper(x,params)

    # calculate boundary constraint and jacobian
    boundary_con = boundary_wrapper(x,params)

    # combine constaint values and jacobians into overall constaint value and jacobian arrays
    c = [spacing_con; boundary_con]
    g[:] = c[:]

    # calculate the objective function and jacobian (negative sign in order to maximize AEP)
    AEP = -aep_wrapper(x,params)[1]

    sparseAEPJacobianUpdateFunctionCalls!(x,params,deriv,numCalls)

    df[:] = -deriv.Jacobian
    ds_dx = ForwardDiff.jacobian(spacing_wrapper, x)
    db_dx = ForwardDiff.jacobian(boundary_wrapper, x)
    dcdx = [ds_dx; db_dx]
    dg[:] = dcdx[:]

    return AEP
end

## update pattern by jacobian cutoff
function sparseJacobianFromCutoff!(x::Vector{Float64},params,deriv,numCalls,cut)
    deriv.Jacobian .= 0
    maxCalls = 0
    numCalls += 1
    for i = 1:length(params.windresource.wind_directions)
        # calculate pattern
        if numCalls > maxCalls
            calculateSparsityPatternCutoff!(x,deriv.pattern,params,i,cut)
            deriv.patterns[:,:,i] .= deriv.pattern
            allocateNewCache!(x,deriv,i)
        end

        # calculate state Jacobian
        calculateSparseJacobian!(x,deriv,i)

        # sum jacobian columns and add to other states
        deriv.Jacobian .= deriv.Jacobian .+ transpose(sum(deriv.jac,dims = 1) .* params.windresource.wind_probabilities[i])
    end
    if numCalls > maxCalls
        numCalls = 0
    end
    deriv.jac .= deriv.jac .* params.obj_scale .* 365.25 .* 24
    deriv.Jacobian .= deriv.Jacobian .* params.obj_scale .* 365.25 .* 24
end

function sparseJacobianFromColors!(x::Vector{Float64},params,deriv,numCalls,numColors)
    deriv.Jacobian .= 0
    maxCalls = 0
    numCalls += 1
    for i = 1:length(params.windresource.wind_directions)
        # calculate pattern
        if numCalls > maxCalls
            calculateSparsityPatternColors!(x,deriv.pattern,params,i,numColors)
            deriv.patterns[:,:,i] .= deriv.pattern
            allocateNewCache!(x,deriv,i)
        end

        # calculate state Jacobian
        calculateSparseJacobian!(x,deriv,i)

        # sum jacobian columns and add to other states
        deriv.Jacobian .= deriv.Jacobian .+ transpose(sum(deriv.jac,dims = 1) .* params.windresource.wind_probabilities[i])
    end
    if numCalls > maxCalls
        numCalls = 0
    end
    deriv.jac .= deriv.jac .* params.obj_scale .* 365.25 .* 24
    deriv.Jacobian .= deriv.Jacobian .* params.obj_scale .* 365.25 .* 24
end

function calculateSparsityPatternColors!(x::Vector{Float64},pattern::Matrix{Float64},params,currentState::Int64,numColors)
    p_wrapper(x) = SparseJacobian.p_wrapper(x,params,currentState)
    dense = abs.(ForwardDiff.jacobian(p_wrapper,x))
    maxJac = maximum(dense)
    if numColors == length(x)
        pattern .= dense
        return
    end

    function reduceColors(cut)
        if cut == 1.0
            return -numColors
        elseif cut == 0.0
            return length(x) - numColors
        end
        pattern .= dense
        pattern[pattern .<= (maxJac*cut)] .= 0
        if maximum(pattern) > 0
            return maximum(matrix_colors(sparse(pattern))) - numColors
        end
        return -numColors
    end

    left = 0.0
    right = 1.0
    mid = 0.5
    for i = 1:500
        num = reduceColors(mid)
        if num < 0
            right = mid
        elseif num > 0
            left = mid
        else
            right = mid
            # println("Bracket: ",i)
            break
        end
        mid = (left+right)/2
    end

    for k = 1:500
        mid = (left+right)/2
        num = reduceColors(mid)
        if num < 1
            right = mid
        else
            # println("Refine 1: ",k)
            break
        end
        mid = (left+right)/2
    end

    for k = 1:500
        mid = (left+right)/2
        num = reduceColors(mid)
        if num > 0
            left = mid
        else
            # println("Refine 2: ",k)
            break
        end
        mid = (left+right)/2
    end

    start = left
    pattern .= copy(dense)
    pattern[pattern .< (maxJac*(start))] .= 0
    # println("Colors: ",maximum(matrix_colors(sparse(pattern))))

    i = 0
    while true
        if maximum(pattern) == 0
            break
        end
        if maximum(matrix_colors(sparse(pattern))) <= numColors
            break
        end
        i += 1
        pattern[pattern .== 0] .= Inf
        pattern[findmin(pattern)[2]] = 0
        pattern[pattern .== Inf] .= 0
        if i > 200
            break
        end
    end
    # println("Remove: ",i)
end

function calculateSparsityPatternCutoff!(x::Vector{Float64},pattern::Matrix{Float64},params,currentState::Int64,cut)
    p_wrapper(x) = SparseJacobian.p_wrapper(x,params,currentState)
    dense = log10.(abs.(ForwardDiff.jacobian(p_wrapper,x)) .+ 1E-200)
    maxJac = maximum(dense)
    pattern .= copy(dense)
    pattern[pattern .< (maxJac-cut)] .= 0
end

function allocateContainersCutoff(x::Vector{Float64},params,cut)
    n = Int(length(x)/2)
    numStates = length(params.windresource.wind_directions)
    pattern = zeros(n,n*2)
    jac = spzeros(n,n*2)
    Jacobian = zeros(n*2)
    patterns = zeros(n,n*2,numStates)
    caches = Array{Any,1}(undef,numStates)
    functions = Array{Any,1}(undef,numStates)
    colors = zeros(Int,n*2)
    index = zeros(Int,n)

    for i = 1:numStates
        calculateSparsityPatternCutoff!(x,pattern,params,i,cut)
        patterns[:,:,i] .= pattern
        jac .= dropzeros(sparse(patterns[:,:,i]))
        if iszero(jac)
            colors = collect(Int64,1:length(colors))
        else
            colors = matrix_colors(jac)
        end
        f = (wt_power,x) -> p_wrapper!(wt_power,x,params,i)
        functions[i] = f
        caches[i] = ForwardColorJacCache(f,x,
                              dx = similar(x[1:n]),
                              colorvec=colors,
                              sparsity = jac)
    end
    deriv = deriv_struct(pattern,jac,Jacobian,patterns,caches,index,functions)
    return deriv
end

function allocateContainersColors(x::Vector{Float64},params,numColors)
    n = Int(length(x)/2)
    numStates = length(params.windresource.wind_directions)
    pattern = zeros(n,n*2)
    jac = spzeros(n,n*2)
    Jacobian = zeros(n*2)
    patterns = zeros(n,n*2,numStates)
    caches = Array{Any,1}(undef,numStates)
    functions = Array{Any,1}(undef,numStates)
    colors = zeros(Int,n*2)
    index = zeros(Int,n)

    for i = 1:numStates
        calculateSparsityPatternColors!(x,pattern,params,i,numColors)
        patterns[:,:,i] .= pattern
        jac .= dropzeros(sparse(patterns[:,:,i]))
        if iszero(jac)
            colors = collect(Int64,1:length(colors))
        else
            colors = matrix_colors(jac)
        end
        f = (wt_power,x) -> p_wrapper!(wt_power,x,params,i)
        functions[i] = f
        caches[i] = ForwardColorJacCache(f,x,
                              dx = similar(x[1:n]),
                              colorvec=colors,
                              sparsity = jac)
    end
    deriv = deriv_struct(pattern,jac,Jacobian,patterns,caches,index,functions)
    return deriv
end

function opt_colors!(g,df,dg,x,params,deriv,numCalls)
    # calculate spacing constraint value and jacobian
    spacing_con = spacing_wrapper(x,params)

    # calculate boundary constraint and jacobian
    boundary_con = boundary_wrapper(x,params)

    # combine constaint values and jacobians into overall constaint value and jacobian arrays
    c = [spacing_con; boundary_con]
    g[:] = c[:]

    # calculate the objective function and jacobian (negative sign in order to maximize AEP)
    AEP = -aep_wrapper(x,params)[1]

    colors = Int64(length(x)/2)
    sparseJacobianFromCutoff!(x,params,deriv,numCalls,colors)

    df[:] = -deriv.Jacobian
    ds_dx = ForwardDiff.jacobian(spacing_wrapper, x)
    db_dx = ForwardDiff.jacobian(boundary_wrapper, x)
    dcdx = [ds_dx; db_dx]
    dg[:] = dcdx[:]

    return AEP
end

end


text:
using SoapySDR
using Documenter

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

makedocs(;
    modules=[SoapySDR],
    authors="JuliaTelecom and contributors",
    repo="https://github.com/JuliaTelecom/SoapySDR.jl/blob/{commit}{path}#{line}",
    sitename="SoapySDR.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://JuliaTelecom.github.io/SoapySDR.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Tutorial" => "tutorial.md",
        "High Level API" => "highlevel.md",
        "Low Level API" => "lowlevel.md",
        "Troubleshooting" => "troubleshooting.md",
    ],
)

deploydocs(;
    repo="github.com/JuliaTelecom/SoapySDR.jl",
    devbranch = "main"
)


text:
# fibonacci seq - LaTeX
#JL Plots

Pkg.add("LaTeXStrings")
Pkg.add("Measures")
using LaTeXStrings, Measures
gr(bg = :white)

fib = zeros(12);
for i = 1:12
    fib[i] = (((1+sqrt(5))/2)^i - ((1-sqrt(5))/2)^i)/sqrt(5);
end

Plots.plot(fib,
    marker=:circle,
    linewidth= 2,
    xlabel=L"n",
    ylabel=L"F_n",
    annotation=(5, 100, L"F_n = \frac{1}{\sqrt{5}} \left[\left( \frac{1+\sqrt{5}}{2} \right)^n - \left( \frac{1-\sqrt{5}}{2} \right)^n \right]"))

max_stars_repo_name: Ellon-M/visualizations

max_issues_repo_name: Ellon-M/visualizations

max_forks_repo_name: Ellon-M/visualizations


text:
using Test;

using AdalmPluto;

@testset "libIIO/buffer.jl" begin
    toggleNoAssertions(true);

    try
        global uri          = AdalmPluto.scan("usb", 1, false);
        global ctx          = AdalmPluto.createContext(uri);
        global tx, rx       = AdalmPluto.findTRXDevices(ctx);
        global txc, rxc     = AdalmPluto.findTRXChannels(ctx);
        global txloc, rxloc = AdalmPluto.findLoChannels(ctx);
        global txi, txq     = AdalmPluto.findIQChannels(tx, "voltage0", "voltage1", true);
        global rxi, rxq     = AdalmPluto.findIQChannels(rx, "voltage0", "voltage1", false);

        C_iio_channel_enable(txi);
        C_iio_channel_enable(txq);
        C_iio_channel_enable(rxi);
        C_iio_channel_enable(rxq);
    catch
        @error "Could not get a valid uri or context or channels, skipping the remaining tests.";
        C_iio_context_destroy(ctx);
        return;
    end

    #  C_iio_device_create_buffer (device::Ptr, …)
    @test (global buf = C_iio_device_create_buffer(rx, UInt(1024), false)) != C_NULL;
    @testset "Buffer properties" begin
        #  C_iio_buffer_get_device (buffer::Ptr, …)
        @test C_iio_buffer_get_device(buf) == rx;
        #  C_iio_buffer_get_poll_fd (buffer::Ptr, …)
        @test_skip C_iio_buffer_get_poll_fd(buf) > 0;                   # The function returns -38: Function not implemented
        #  C_iio_buffer_set_blocking_mode (buffer::Ptr, …)
        @test_skip C_iio_buffer_set_blocking_mode(buf, false) == 0;     # The function returns -38: Function not implemented
        @test_skip C_iio_buffer_set_blocking_mode(buf, true) == 0;      # The function returns -38: Function not implemented
        # By default devices are blocking

        saved_data_ptr = C_iio_buffer_get_data(buf);
        test_data = b"Test data\0";
        #  C_iio_buffer_set_data (buffer::Ptr, …)
        @test C_iio_buffer_set_data(buf, Ptr{Cvoid}(pointer(test_data))) === nothing;
        #  C_iio_buffer_get_data (buffer::Ptr, …)
        @test C_iio_buffer_get_data(buf) == Ptr{Cvoid}(pointer(test_data));
        C_iio_buffer_set_data(buf, saved_data_ptr);

        #  C_iio_buffer_start (buffer::Ptr, …)
        @test C_iio_buffer_start(buf) != C_NULL;
        #  C_iio_buffer_first (buffer::Ptr, …)
        @test C_iio_buffer_first(buf, rxi) != C_NULL;
        #  C_iio_buffer_step (buffer::Ptr, …)
        @test C_iio_buffer_step(buf) > 0;
        #  C_iio_buffer_end (buffer::Ptr, …)
        @test C_iio_buffer_end(buf) != C_NULL;
    end

    @testset "Buffer read/write" begin
        #  C_iio_buffer_refill (buffer::Ptr, …)
        @test C_iio_buffer_refill(buf) > 0;

        #  C_iio_channel_read! (chn::Ptr, …)
        smp = zeros(UInt8, 512 * C_iio_device_get_sample_size(rx));
        @test C_iio_channel_read!(rxi, buf, smp) == 512 * C_iio_device_get_sample_size(rx);

        # Testing those implies you can listen to what the radio emits
        #  C_iio_buffer_push (buffer::Ptr, …)
        @test_skip C_iio_buffer_push(buf) > 0;
        #  C_iio_buffer_push_partial (buffer::Ptr, …)
        @test_skip C_iio_buffer_push_partial(buf, samples_count) > 0;

        #  C_iio_buffer_cancel (buffer::Ptr, …)
        @test_skip C_iio_buffer_cancel(buf) === nothing;
    end

    @testset "Placeholders" begin
        #  C_iio_buffer_foreach_sample
        @test_skip C_iio_buffer_foreach_sample() === nothing;
    end

    C_iio_channel_disable(txi);
    C_iio_channel_disable(txq);
    C_iio_channel_disable(rxi);
    C_iio_channel_disable(rxq);

    #  C_iio_buffer_destroy (buffer::Ptr, …)
    @test C_iio_buffer_destroy(buf) === nothing;

    C_iio_context_destroy(ctx);
end


text:
# This script is licensed under the MIT License (MIT).

if !Sys.isunix() || Sys.isapple()
    for envvar_name in ["JUPYTER", "PYTHON"]
        if !haskey(ENV, envvar_name)
            @error "Please set environment variable $envvar_name, run \"julia system_check.jl\" for details."
            exit(1)
        end
    end
end

using Pkg

Pkg.add("ArraysOfArrays")
Pkg.add("BufferedStreams")
Pkg.add("CoordinateTransformations")
Pkg.add("CSV")
Pkg.add("DataFrames")
Pkg.add("Distributions")
Pkg.add("DSP")
Pkg.add("ElasticArrays")
Pkg.add("FFTW")
Pkg.add("FillArrays")
Pkg.add("GPUArrays")
Pkg.add("HDF5")
Pkg.add("IJulia")
Pkg.add("Interact")
Pkg.add("LaTeXStrings")
Pkg.add("LsqFit")
Pkg.add("Measurements")
Pkg.add("Optim")
Pkg.add("ParallelProcessingTools")
Pkg.add("Plots")
Pkg.add("PyCall")
Pkg.add("PyPlot")
Pkg.add("RadiationSpectra")
Pkg.add("RecipesBase")
Pkg.add("SolidStateDetectors")
Pkg.add("SpecialFunctions")
Pkg.add("SplitApplyCombine")
Pkg.add("StaticArrays")
Pkg.add("StatsBase")
Pkg.add("Tables")
Pkg.add("TypedTables")
Pkg.add("Unitful")
Pkg.add("UnsafeArrays")
Pkg.add("Weave")

Pkg.add(PackageSpec(url="https://github.com/JuliaHEP/RadiationDetectorSignals.jl.git", rev="master"))
Pkg.add(PackageSpec(url="https://github.com/JuliaHEP/RadiationDetectorDSP.jl.git", rev="master"))

Pkg.add(PackageSpec(url="https://github.com/legend-exp/LegendDataTypes.jl.git", rev="master"))
Pkg.add(PackageSpec(url="https://github.com/legend-exp/LegendDSP.jl.git", rev="master"))
Pkg.add(PackageSpec(url="https://github.com/legend-exp/LegendHDF5IO.jl.git", rev="master"))
Pkg.add(PackageSpec(url="https://github.com/legend-exp/LegendTextIO.jl.git", rev="master"))

pkg"""precompile"""


text:
using Distributed

println("nprocs = ", nprocs())

@everywhere using BulletsMaggots
@everywhere build_table(false)

function jonah_parallel(reps::Int = 1000, full_list_limit::Int = 0)

    println(
        "Experiment with $reps repetitions and full list used for $full_list_limit of the first moves",
    )

    total = @distributed (+) for i = 1:reps
        jonah_solver(guess_any(), false, full_list_limit)
    end
    return total / reps
end

function bm_parallel(reps::Int = 1000)
    total = @distributed (+) for i = 1:reps
        bm_solver(guess_any(), false)
    end
    return total / reps
end


text:
using DSGE, JLD2, DataStructures
import Test: @test, @testset
import DSGE: @test_matrix_approx_eq

m = KrusellSmith()

# Load pre-reduced canonical form matrices
g1 = load("../../test_outputs/solve/canonical.jld2", "g1")
g0 = load("../../test_outputs/solve/canonical.jld2", "g0")
psi = load("../../test_outputs/solve/canonical.jld2", "psi")
pi = load("../../test_outputs/solve/canonical.jld2", "pi")
c = load("../../test_outputs/solve/canonical.jld2", "c")

# Load reduced canonical form matrices
g1_kry = load("../../test_outputs/solve/canonical_kry.jld2", "g1")
g0_kry = load("../../test_outputs/solve/canonical_kry.jld2", "g0")
psi_kry = load("../../test_outputs/solve/canonical_kry.jld2", "psi")
pi_kry = load("../../test_outputs/solve/canonical_kry.jld2", "pi")
c_kry = load("../../test_outputs/solve/canonical_kry.jld2", "c")

# Load Krylov + spline reduced canonical form matrices
g1_spl = load("../../test_outputs/solve/canonical_spl.jld2", "g1")
g0_spl = load("../../test_outputs/solve/canonical_spl.jld2", "g0")
psi_spl = load("../../test_outputs/solve/canonical_spl.jld2", "psi")
pi_spl = load("../../test_outputs/solve/canonical_spl.jld2", "pi")
c_spl = load("../../test_outputs/solve/canonical_spl.jld2", "c")

# Update model
update!(m.settings[:I], Setting(:I, 100))
update!(m.settings[:amin], Setting(:amin, 0.))
update!(m.settings[:amax], Setting(:amax, 100.))
update!(m.settings[:a], Setting(:a, collect(linspace(0., 100., 100))))
update!(m.settings[:da], Setting(:da, 100./(100 - 1)))
update!(m.settings[:aaa], Setting(:aaa, reshape(get_setting(m, :aa), 2 * 100, 1)))
update!(m.settings[:J], Setting(:J, 2))
update!(m.settings[:z], Setting(:z, [0., 1.]))
update!(m.settings[:dz], Setting(:dz, 1.))
update!(m.settings[:zz], Setting(:zz, ones(100, 1) * get_setting(m, :z)'))
update!(m.settings[:zzz], Setting(:zzz, reshape(get_setting(m, :zz), 2 * 100, 1)))
update!(m.settings[:n_jump_vars], Setting(:n_jump_vars, 200))
update!(m.settings[:n_state_vars], Setting(:n_state_vars, 199 + 1))
update!(m.settings[:n_state_vars_unreduce], Setting(:n_state_vars_unreduce, 0))
update!(m.settings[:reduce_state_vars], Setting(:reduce_state_vars, true))
update!(m.settings[:reduce_v], Setting(:reduce_v, true))
update!(m.settings[:krylov_dim], Setting(:krylov_dim, 5))
update!(m.settings[:n_knots], Setting(:n_knots, 12))
update!(m.settings[:c_power], Setting(:c_power, 7))
update!(m.settings[:knots_dict], Setting(:knots_dict, Dict(1 => collect(linspace(get_setting(m, :amin), get_setting(m, :amax), 12 - 1)))))
update!(m.settings[:spline_grid], Setting(:spline_grid, get_setting(m, :a)))
update!(m.settings[:n_prior], Setting(:n_prior, 1))
update!(m.settings[:n_post], Setting(:n_post, 2))
update!(m.settings[:F], Setting(:F, identity))

# Test non-sparse methods
Γ0, Γ1, Ψ, Π, C, ~, ~ = krylov_reduction(m, full(g0), full(g1), full(psi), full(pi), full(c))
@testset "Krylov Reduction (non-sparse)" begin
    @test @test_matrix_approx_eq full(g1_kry) Γ1
    @test @test_matrix_approx_eq full(g0_kry) Γ0
    @test vec(full(psi_kry)) ≈ vec(full(Ψ)) # allow machine error
    @test vec(full(pi_kry)) ≈ vec(full(Π))
    @test vec(full(c_kry)) ≈ vec(full(C))
end

Γ0, Γ1, Ψ, Π, C, ~, ~ = valuef_reduction(m, full(g0_kry), full(g1_kry), full(psi_kry), full(pi_kry), full(c_kry))
@testset "Value Function Reduction via Spline Basis (non-sparse)" begin
    @test @test_matrix_approx_eq full(g1_spl) Γ1
    @test @test_matrix_approx_eq full(g0_spl) Γ0
    @test vec(full(psi_spl)) ≈ vec(full(Ψ))
    @test vec(full(pi_spl)) ≈ vec(full(Π))
    @test vec(full(c_spl)) ≈ vec(full(C))
end

# Test sparse methods
Γ0, Γ1, Ψ, Π, C, ~, ~ = krylov_reduction(m, sparse(g0), sparse(g1), sparse(psi), sparse(pi), sparse(c))
@testset "Krylov Reduction (sparse)" begin
    @test typeof(Γ1) == SparseMatrixCSC{Float64,Int64}
    @test typeof(Γ0) == SparseMatrixCSC{Float64,Int64}
    @test typeof(Ψ) == SparseMatrixCSC{Float64,Int64}
    @test typeof(Π) == SparseMatrixCSC{Float64,Int64}
    @test typeof(C) == SparseMatrixCSC{Float64,Int64}
    @test @test_matrix_approx_eq full(g1_kry) full(Γ1)
    @test @test_matrix_approx_eq full(g0_kry) full(Γ0)
    @test vec(full(psi_kry)) ≈ vec(full(Ψ)) # allow machine error
    @test vec(full(pi_kry)) ≈ vec(full(Π))
    @test vec(full(c_kry)) ≈ vec(full(C))
end

Γ0, Γ1, Ψ, Π, C, ~, ~ = valuef_reduction(m, sparse(g0_kry), sparse(g1_kry), sparse(psi_kry),
                                         sparse(pi_kry), sparse(c_kry))
@testset "Value Function Reduction via Spline Basis (sparse)" begin
    @test typeof(Γ1) == SparseMatrixCSC{Float64,Int64}
    @test typeof(Γ0) == SparseMatrixCSC{Float64,Int64}
    @test typeof(Ψ) == SparseMatrixCSC{Float64,Int64}
    @test typeof(Π) == SparseMatrixCSC{Float64,Int64}
    @test typeof(C) == SparseMatrixCSC{Float64,Int64}
    @test @test_matrix_approx_eq full(g1_spl) full(Γ1)
    @test @test_matrix_approx_eq full(g0_spl) full(Γ0)
    @test vec(full(psi_spl)) ≈ vec(full(Ψ))
    @test vec(full(pi_spl)) ≈ vec(full(Π))
    @test vec(full(c_spl)) ≈ vec(full(C))
end


text:
module Connectivity
    include("DataTestTypes.jl") # How do I make it obvious that block comes from DataTestTypes
    # using .DataTestTypes: Block
    import .DataTestTypes: Block
    export add_blocks, Block

    function add_blocks(block1::Block, block2::Block)
        return block1.a + block2.b
    end

end

max_stars_repo_name: pjuangph/JuliaLearning

max_issues_repo_name: pjuangph/JuliaLearning

max_forks_repo_name: pjuangph/JuliaLearning


text:
using Plots;
using Random;
using LinearAlgebra;
using Statistics;
pyplot();

#===============================================#
#                                               #
#                Tools functions                #
#                                               #
#===============================================#

function generate_dataset(sample_size::Int, feature_dim::Int, num_classes::Int)
    CORE_BIAS = 5
    # random get
    Y = rand(1:num_classes,sample_size)
    # make features
    X = randn(sample_size,feature_dim) .+ (CORE_BIAS .* Y)
    # make label
    Y = Int.([y == class for y in Y, class = 1:num_classes])
    return  X, Y
end

function rand_center(data,k)
    n = size(data,2)
    centroids = zeros(k,n)

    for i in 1:n
        dmin, dmax = minimum(data[1:end,i]), maximum(data[1:end,i])
        centroids[:,i] = dmin .+ (dmax - dmin) * rand(k)
    end
    return centroids
end

function draw(X,label,centroids)
    # draw the line and points
    color_list = [:red, :green, :blue, :yellow, :pink]
    colors = [color_list[l] for l in label]
    scatter(X[:,1], X[:,2],  markercolor = colors)
    scatter!(centroids[:,1],centroids[:,2],markercolor=color_list,marker=10)
end

#===============================================#
#                                               #
#                     K-Means                   #
#                                               #
#===============================================#

distance(P1,P2) = sqrt(sum((P1 - P2).^2))

function kmeans(X,NUM_CLASSES)
    centroids = rand_center(X,NUM_CLASSES)      # center points
    labels = zeros(Int,size(X,1))               # labels
    assignment = zeros(size(X,1))               # 方差

    while true
        old_labels = deepcopy(labels)
        for i in 1:size(X,1)
            min_dist, min_index = Inf, -1
            for j in 1:NUM_CLASSES
                dist = distance(X[i,:],centroids[j,:])
                if dist < min_dist
                    min_dist, min_index = dist, j
                    labels[i] = j
                    assignment[i] = dist^2
                end
            end
        end

        for m in 1:NUM_CLASSES
            M_class = X[[i for (i,value) in enumerate(labels) if value == m],:]
            centroids[m,:] = mean(M_class,dims = 1)
        end
        if old_labels == labels
            break
        end
    end
    return centroids,labels
end

SAMPLE_SIZE = 50
FEATURE_DIM = 2
NUM_CLASSES = 2
X,Y = generate_dataset(SAMPLE_SIZE,FEATURE_DIM,NUM_CLASSES)
centroids,labels = kmeans(X,NUM_CLASSES)
draw(X,labels,centroids)

#===============================================#
#                                               #
#                       EM                      #
#                                               #
#===============================================#

function Init_params(Data, K)
    N, D = size(Data)            # 元素数目，元素特性数目
    πs = zeros(1, K)             # k类GMM所占权重（influence factor）
    μs = rand_center(Data,K)     # k*D, 即k类的中心点
    Σs = zeros(K, D, D)          # k类GMM的协方差矩阵，每个是D*D的

    Dist = repeat(sum(Data .^ 2, dims = 2), outer = (1,K)) + repeat(sum(μs' .^ 2, dims = 1), outer = (N,1)) - 2Data * μs'
    labels = mapslices(argmin, Dist, dims = 2)

    for k in 1:K
        Xk = Data[[index for (index,value) in enumerate(labels) if value == k],:]
        πs[k] = size(Xk,1) / N
        Σs[k,:,:] = cov(Xk)
    end
    return μs, πs, Σs
end

function Gaussian(Data,μ,Σ)
    D = size(Data,1)          # 元素数目，元素特性数目
    detΣ = det(Σ)
    if detΣ == 0              # 以防行列式为0
        Σ = Σ + Matrix{Float64}(I, D, D) .* 0.01
        detΣ = det(Σ)
    end
    dist = (Data - μ)
    return (2π)^-0.5D * detΣ^-0.5 * exp(-0.5 * dist' * Σ^-1 * dist)
end

function γ(Data,πs,μs,Σs,K)
    N, D = size(Data)         # 元素数目，元素特性数目
    γs = zeros(N,K)
    for i = 1:N
        res = [πs[k] * Gaussian(Data[i,:],μs[k,:],Σs[k,:,:]) for k in 1:K]
        γs[[i],:] = res ./ sum(res)
    end
    return γs
end

function GMM(Data,K)
    N, D = size(Data)
    μs, πs, Σs = Init_params(Data, K)

    losses = Float64[]
    γs = zeros(N,K)
    for i in 1:20
        # E
        γs = γ(Data,πs,μs,Σs,K)

        # M
        Nk = sum(γs, dims = 1)
        πs = Nk ./ sum(Nk)

        for k in 1:K
            μs[[k],:] = Nk[k]^-1 .* sum(γs[:,k] .* Data, dims = 1)

            dist = γs[:,[k]].^0.5 .* (Data .- μs[[k],:])
            Σs[k,:,:] = Nk[k]^-1 .* dist' * dist
        end

        loglikelyhood = sum([log(sum([πs[k] * Gaussian(Data[n,:],μs[k,:],Σs[k,:,:]) for k in 1:K])) for n in 1:N])
        push!(losses,-loglikelyhood)
    end

    return μs, mapslices(argmax,γs,dims = 2), losses
end

SAMPLE_SIZE = 100
FEATURE_DIM = 2
NUM_CLASSES = 3
X,Y = generate_dataset(SAMPLE_SIZE,FEATURE_DIM,NUM_CLASSES)
centroids,labels,losses = GMM(X,NUM_CLASSES)
draw(X,labels,centroids)
plot(1:length(losses),losses)

#===============================================#
#                                               #
#                       UCI                     #
#                                               #
#===============================================#

using DataFrames, CSV;
function toLabel(labels)
    classes = unique(labels)
    function λ(label)
        for (index, value) in enumerate(classes)
            if label == value
                return index
            end
        end
    end
    return Int.([y == class for y in λ.(labels), class = 1:size(classes,1)]),classes
end

datasets = CSV.read("data/iris-cut.csv");
labels = Array{String}(copy(datasets[:Species]))

X = Matrix{Float64}(datasets[:,1:2])
Y,classes = toLabel(labels)
real_label = (Y[:,1] .+ 1)

centroids,labels,losses = GMM(X,2)
draw(X,labels,centroids)
plot(1:length(losses),losses)
sum(Int.(real_label .== labels))

max_stars_repo_name: AlongWY/HIT-WorkShop

max_issues_repo_name: AlongWY/HIT-WorkShop

max_forks_repo_name: AlongWY/HIT-WorkShop


text:
@testset "Testing interp.jl" begin

    # uniform interpolation
    breaks = range(-3, stop=3, length=100)
    vals = exp.(breaks)
    vals2 = [vals sin.(breaks)]

    li = interp(breaks, vals)
    li2 = LinInterp(breaks, vals)

    li_mat = interp(breaks, vals2)
    li_mat2 = LinInterp(breaks, vals2)

    # test constructor
    @test li == li2
    @test li_mat == li_mat2

    # make sure evaluation is inferrable
    for T in (Float64, Float32, Float16, Int64, Int32, Int16)
        @test begin
            @inferred li(one(T))
            true
        end
        @test begin
            @inferred li_mat(one(T))
            true
        end
    end

    # on grid is exact
    for i in 1:length(breaks)
        @test abs(li(breaks[i]) - vals[i]) < 1e-15
        @test all(abs.(li_mat(breaks[i]) - vals2[i, :] .< 1e-15))
    end

    # off grid is close
    for x in range(-3, stop=3, length=300)
        @test abs(li(x) - exp(x)) < 1e-2
        @test all(abs.(li_mat(x) .- [exp(x), sin(x)]) .< 1e-2)
        @test li(x) ≈ li_mat(x, 1)
    end

    # test errors for col spec for li_mat being wrong
    @test_throws BoundsError li_mat(0.5, 0)
    @test_throws BoundsError li_mat(0.5, 3)
    @test_throws BoundsError li_mat(0.5, [0, 1])
    @test_throws BoundsError li_mat(0.5, [2, 3])

    # non-uniform
    breaks = cumsum(0.1 .* rand(20))
    vals = 0.1 .* map(sin, breaks)
    li = interp(breaks, vals)
    li_mat = interp(breaks, [vals vals .+ 1])

    # on grid is exact
    for i in 1:length(breaks)
        @test abs(li(breaks[i]) - vals[i]) < 1e-15
        @test all(abs.(li_mat(breaks[i]) .- [vals[i], vals[i]+1]) .< 1e-15)
    end

    # off grid is close
    for x in range(minimum(breaks), stop = maximum(breaks), length = 30)
        @test abs(li(x) - 0.1*sin(x)) < 1e-2
        @test all(abs.(li_mat(x) - [0.1*sin(x), 0.1*sin(x)+1]) .< 1e-2)

    end

    # un-sorted works for `interp` function, but not `LinInterp`
    breaks = rand(10)
    vals = map(sin, breaks)

    @inferred interp(breaks, vals)
    @inferred interp(breaks, [vals vals .+ 1])
    @test_throws ArgumentError LinInterp(breaks, vals)
    @test_throws ArgumentError LinInterp(breaks, [vals vals .+ 1])

    # dimension mismatch
    breaks = cumsum(rand(10))
    vals = rand(8)

    @test_throws DimensionMismatch interp(breaks, vals)
    @test_throws DimensionMismatch LinInterp(breaks, vals)

    @test_throws DimensionMismatch interp(breaks, [vals vals .+ 1])
    @test_throws DimensionMismatch LinInterp(breaks, [vals vals .+ 1])

end  # @testset


text:
module VectorWithEntryRemovedTests

using Test
using Gridap.Arrays

a = collect(1:10)
for i=1:length(a)
   b    = VectorWithEntryRemoved(a,i)
   test_array(vcat(a[1:i-1],a[i+1:end]),b)
end

for i=1:length(a)
  b    = VectorWithEntryRemoved(a,i)
  b[rand(1:length(a)-1)]=rand(Int)
  test_array(vcat(a[1:i-1],a[i+1:end]),b)
end

end # module


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct GenericBattery <: Storage
        name::String
        available::Bool
        bus::Bus
        primemover::PrimeMovers
        energy::Float64
        capacity::Min_Max
        rating::Float64
        activepower::Float64
        inputactivepowerlimits::Min_Max
        outputactivepowerlimits::Min_Max
        efficiency::NamedTuple{(:in, :out), Tuple{Float64, Float64}}
        reactivepower::Float64
        reactivepowerlimits::Union{Nothing, Min_Max}
        _forecasts::InfrastructureSystems.Forecasts
        internal::InfrastructureSystemsInternal
    end

Data structure for a generic battery

# Arguments
- `name::String`
- `available::Bool`
- `bus::Bus`
- `primemover::PrimeMovers`: PrimeMover Technology according to EIA 923
- `energy::Float64`: State of Charge of the Battery p.u.-hr
- `capacity::Min_Max`: Maximum and Minimum storage capacity in p.u.-hr
- `rating::Float64`
- `activepower::Float64`
- `inputactivepowerlimits::Min_Max`
- `outputactivepowerlimits::Min_Max`
- `efficiency::NamedTuple{(:in, :out), Tuple{Float64, Float64}}`
- `reactivepower::Float64`
- `reactivepowerlimits::Union{Nothing, Min_Max}`
- `_forecasts::InfrastructureSystems.Forecasts`
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct GenericBattery <: Storage
    name::String
    available::Bool
    bus::Bus
    "PrimeMover Technology according to EIA 923"
    primemover::PrimeMovers
    "State of Charge of the Battery p.u.-hr"
    energy::Float64
    "Maximum and Minimum storage capacity in p.u.-hr"
    capacity::Min_Max
    rating::Float64
    activepower::Float64
    inputactivepowerlimits::Min_Max
    outputactivepowerlimits::Min_Max
    efficiency::NamedTuple{(:in, :out), Tuple{Float64, Float64}}
    reactivepower::Float64
    reactivepowerlimits::Union{Nothing, Min_Max}
    _forecasts::InfrastructureSystems.Forecasts
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function GenericBattery(name, available, bus, primemover, energy, capacity, rating, activepower, inputactivepowerlimits, outputactivepowerlimits, efficiency, reactivepower, reactivepowerlimits, _forecasts=InfrastructureSystems.Forecasts(), )
    GenericBattery(name, available, bus, primemover, energy, capacity, rating, activepower, inputactivepowerlimits, outputactivepowerlimits, efficiency, reactivepower, reactivepowerlimits, _forecasts, InfrastructureSystemsInternal())
end

function GenericBattery(; name, available, bus, primemover, energy, capacity, rating, activepower, inputactivepowerlimits, outputactivepowerlimits, efficiency, reactivepower, reactivepowerlimits, _forecasts=InfrastructureSystems.Forecasts(), )
    GenericBattery(name, available, bus, primemover, energy, capacity, rating, activepower, inputactivepowerlimits, outputactivepowerlimits, efficiency, reactivepower, reactivepowerlimits, _forecasts, )
end

# Constructor for demo purposes; non-functional.

function GenericBattery(::Nothing)
    GenericBattery(;
        name="init",
        available=false,
        bus=Bus(nothing),
        primemover=BA::PrimeMovers,
        energy=0.0,
        capacity=(min=0.0, max=0.0),
        rating=0.0,
        activepower=0.0,
        inputactivepowerlimits=(min=0.0, max=0.0),
        outputactivepowerlimits=(min=0.0, max=0.0),
        efficiency=(in=0.0, out=0.0),
        reactivepower=0.0,
        reactivepowerlimits=(min=0.0, max=0.0),
        _forecasts=InfrastructureSystems.Forecasts(),
    )
end

"""Get GenericBattery name."""
get_name(value::GenericBattery) = value.name
"""Get GenericBattery available."""
get_available(value::GenericBattery) = value.available
"""Get GenericBattery bus."""
get_bus(value::GenericBattery) = value.bus
"""Get GenericBattery primemover."""
get_primemover(value::GenericBattery) = value.primemover
"""Get GenericBattery energy."""
get_energy(value::GenericBattery) = value.energy
"""Get GenericBattery capacity."""
get_capacity(value::GenericBattery) = value.capacity
"""Get GenericBattery rating."""
get_rating(value::GenericBattery) = value.rating
"""Get GenericBattery activepower."""
get_activepower(value::GenericBattery) = value.activepower
"""Get GenericBattery inputactivepowerlimits."""
get_inputactivepowerlimits(value::GenericBattery) = value.inputactivepowerlimits
"""Get GenericBattery outputactivepowerlimits."""
get_outputactivepowerlimits(value::GenericBattery) = value.outputactivepowerlimits
"""Get GenericBattery efficiency."""
get_efficiency(value::GenericBattery) = value.efficiency
"""Get GenericBattery reactivepower."""
get_reactivepower(value::GenericBattery) = value.reactivepower
"""Get GenericBattery reactivepowerlimits."""
get_reactivepowerlimits(value::GenericBattery) = value.reactivepowerlimits
"""Get GenericBattery _forecasts."""
get__forecasts(value::GenericBattery) = value._forecasts
"""Get GenericBattery internal."""
get_internal(value::GenericBattery) = value.internal