text:
include("fixtures.jl")

@testset "Basic Example 1" begin
    xin = 1.0:1.0:7.0
    yin = [1.,2.,3.,4.,5.,6.,7.]
    xout = [3.2, 5.2]
    @test nuresample(xin, yin, xout, required_points_per_slice=1) == [3.5, 5.5]
end

@testset "Basic Example 2" begin
    xin = 1.0:1.0:9.0
    # xin=[1, 2, 3, 4, 5, 6, 7, 8, 9]
    yin = [0, 8, 2, 7, 4, 9, 7, 8, 4]
    xout = [3.2, 6] # Δ = 2.8, Δ/2 = 1.4
    # first point: from 1.8 to 4.6 -> 2 3 | 4
    # second points: from 4.6 to 7.4 -> 5 |6 7
    # (middle one should be counted to the higher interval)
    @test nuresample(xin, yin, xout, required_points_per_slice=1) == [
        mean(yin[2:4]),
        mean(yin[5:7])
        ]
end

@testset "Example that Should not Cause NaN Result" begin
    yin = [1.0, 212.47605778632067, 118.55093948572791, 273.88610509333614,
    307.67343485176576, 257.1875847477985, 357.4538672443729, 484.87523154603514,
    400.1796499550531, 398.29054790171745, 317.0589375578311, 501.18723362727246,
    608.6533266004108, 428.7113682873226, 869.7311001362467]
    xin = range(0.0, step=0.18310546955466883, length=length(yin))
    xout = [1.0, 1.0599939505656388]
    output = nuresample(xin, yin, xout)
    @test all(.!isnan.(output))
end

@testset "Balanced Average for Asymmetric Windows" begin
    xout = [2.9, 3.8, 7.5] # middle point highly asymmetric
    # yin=[9, 7, 8, 7, 1, 2, 4, 1, 9, 8, 5, 3]
    #      1--2-*3-*4--5--6--7-*8--9--10-11-12  (output points as "*")
    #            ll*rrrrrrrrrrr                 (point under test)

    # without upsampling
    output = nuresample(xin, yin, xout, rect_window(1.),
        required_points_per_slice=1)
    @test output[2] == mean(yin[3:7])

    # now with upsampling
    output_asymmetric_upsampling = nuresample(xin, yin, xout, rect_window(1.),
        required_points_per_slice=2, upsampling_function=tri_window(1.))
    # expected: left slice only contains 1 point, but 2 are required
    # so both slices should be upsampled with a step of 0.9/2 = 0.45
    # upsampling positions therefore are: [3.125, 3.575] (left) and
    # [4.025, 4.475, 4.925, 5.375, 5.825, 6.275, 6.725, 7.175] (right)
    # the mean of all those points, linearly interpolated, should be the result
    itp = interpolate(yin, BSpline(Linear()))
    expected_positions = [
        3.125, 3.575, 4.025, 4.475, 4.925, 5.375, 5.825, 6.275, 6.725, 7.175]
    @test mean(itp(expected_positions)) ≈ output_asymmetric_upsampling[2]

    # relative change from enabling upsampling should be smaller than 8.2 %
    # (value that is just achieved with implementation at time of writing test)
    rel_change = abs(output_asymmetric_upsampling[2] - output[2])/output[2]
    @test rel_change < .082

    # compare with large oversampling
    big_oversampling = nuresample(xin, yin, xout, rect_window(1.),
    required_points_per_slice=9, upsampling_function=tri_window(1.))[2]
    # relative difference between 2x and 9x oversampling should be smaller 0.96 %
    # (value that is just achieved with implementation at time of writing test)
    rel_change = abs(output_asymmetric_upsampling[2] - big_oversampling)/
        big_oversampling
    @test rel_change < .0096
end

@testset "Large Oversampling Value Should Not Error" begin
    xout = [2.9, 3.8, 7.5] # middle point highly asymmetric

    result8 = nuresample(xin, yin, xout, rect_window(1.),
    required_points_per_slice=8, upsampling_function=tri_window(1.))[2]

    result32 = nuresample(xin, yin, xout, rect_window(1.),
    required_points_per_slice=32, upsampling_function=tri_window(1.))[2]

    # results should be very comparable
    @test abs(result8 - result32)/result32 < .001
end

@testset "Output Value Does Not Jump when Output Point Moves over Input Point" begin
    # output point left of index 6
    output_left = nuresample(xin, yin, [3.9, 5.9, 7.9], rect_window(1.),
        required_points_per_slice=1)
    # output point right of index 6
    output_right = nuresample(xin, yin, [3.9, 6.1, 7.9], rect_window(1.),
        required_points_per_slice=1)
    @test output_left[2] == output_right[2]
end

@testset "Error When not Enough Points at Beginning or End" begin
    # Just enough points at beginning
    nuresample(xin, yin, [2., 3.9999999999999], rect_window(1.),
        required_points_per_slice=1)

    # not enough points (the point at index -1 would be needed but does not exist)
    @test_throws Exception nuresample(xin, yin, [2., 4.], rect_window(1.),
        required_points_per_slice=1)

    # Just enough points at end
    nuresample(xin, yin, [9.00000000000001, 11.], rect_window(1.),
        required_points_per_slice=1)

    # not enough points (the point at index 13 would be needed but does not exist)
    @test_throws Exception nuresample(xin, yin, [9., 11.], rect_window(1.),
        required_points_per_slice=1)
end

@testset "xout not Monotonically Increasing Throws Error" begin
    # not monotonic
    @test_throws(
        ArgumentError("xout must be increasing everywhere. Violated between index 2 and 3: 6.5 >= 6.0"),
        nuresample(xin, yin, [4.5, 6.5, 6., 8.5]))

    # monotonic, but not strict
    @test_throws(
        ArgumentError("xout must be increasing everywhere. Violated between index 2 and 3: 6.5 >= 6.5"),
        nuresample(xin, yin, [4.5, 6.5, 6.5, 8.5]))

    # not monotonic at beginning
    @test_throws(
        ArgumentError("xout must be increasing everywhere. Violated between index 1 and 2: 4.5 >= 3.5"),
        nuresample(xin, yin, [4.5, 3.5, 6., 8.5]))

    # not strict monotonic at beginning
    @test_throws(
        ArgumentError("xout must be increasing everywhere. Violated between index 1 and 2: 4.5 >= 4.5"),
        nuresample(xin, yin, [4.5, 4.5, 6., 8.5]))

    # not monotonic at end
    @test_throws(
        ArgumentError("xout must be increasing everywhere. Violated between index 3 and 4: 6.0 >= 5.5"),
        nuresample(xin, yin, [4.5, 5., 6., 5.5]))

    # not strict monotonic at end
    @test_throws(
        ArgumentError("xout must be increasing everywhere. Violated between index 3 and 4: 6.0 >= 6.0"),
        nuresample(xin, yin, [4.5, 5., 6., 6.]))
end

@testset "Required Input Point per Slice" begin
    # The default for required_points_per_slice should be
    # required_points_per_unit = 4
    # required_points_per_slice = round(required_points_per_unit * smoothing_kernel.width)

    # So, with the default kernel of Rect(.5), at least 2 points have to be in the slice.

    # 1 2 3 4 5 6 7 8 9 10 11 12
    #    *   | *   |
    # (.9:4.1) (4.1:7.3)
    # 2    2|1   2 -> no upsampling on first point, but on second point
    result = nuresample(xin, yin, [2.5, 5.7], upsampling_function=tri_window(1.))
    @test result[1] == mean(yin[1:4])
    itp = interpolate(yin, BSpline(Linear()))
    expected_interpolation_positions = [
        4.1+1.6/2*.5, 4.1+1.6/2*1.5, 5.7+.5*1.6/2, 5.7+1.5*1.6/2
    ]
    @test result[2] == mean(itp(expected_interpolation_positions))

    # When using a larger kernel, more points need to be in the slice
    # e.g. Rect(.75) -> 3 points have to be in the slice
    # 1 2 3 4 5 6 7 8 9 10 11 12
    #      *     *
    # (1.1:5.9) (4.3:9.1)
    # 2    2|2   3 -> upsampling on both points
    result = nuresample(xin, yin, [3.5, 6.7], rect_window(.75),
        upsampling_function=tri_window(1.))
    itp = interpolate(yin, BSpline(Linear()))
    s = 2.4/3 # step
    expected_interpolation_positions1 = [
        1.1+s*.5, 1.1+s*1.5, 1.1+s*2.5, 3.5+s*.5, 3.5+s*1.5, 3.5+s*2.5
    ]
    @test result[1] == mean(itp(expected_interpolation_positions1))
    expected_interpolation_positions2 = [
        4.3+s*.5, 4.3+s*1.5, 4.3+s*2.5, 6.7+s*.5, 6.7+s*1.5, 6.7+s*2.5
    ]
    @test result[2] ≈ mean(itp(expected_interpolation_positions2))

    # with "disabled" interpolation
    no_interpolation = nuresample(xin, yin, [2.5, 5.7], required_points_per_slice = 1)
    @test no_interpolation[1] == mean(yin[1:4])
    @test no_interpolation[2] == mean(yin[5:7])
end

@testset "Strange required_points_per_slice Values" begin
    # disallow values smaller than 0, because slices without values in them do
    # not contribute anything, so the value of a point might become undefined
    # without upsampling, which then is never triggered. If upsampling should
    # never happen, use nearest neighbor interpolation from a different package.
    @test_throws AssertionError nuresample(xin, yin, [4.5, 7.7], required_points_per_slice=0)
    @test_throws AssertionError nuresample(xin, yin, [4.5, 7.7], required_points_per_slice=-1)
    # disallow float values, because they result in discontinuous values with
    # the upsampling algorithm, i.e. the outputs values with 4.5 are very
    # different from 4 or 5
    @test_throws TypeError nuresample(xin, yin, [4.5, 7.7], required_points_per_slice=4.5)
end

@testset "Invalid Window Widths" begin
    # Smoothing Kernel Width 0 => just use the upsampling kernel as smoothing kernel instead
    @test_throws AssertionError nuresample(xin, yin, [2.5, 4.9, 5.0, 6.4],
    rect_window(0),
    upsampling_function=rect_window(.5),
    required_points_per_slice=1)

    # Upsampling Kernel Width < 1 => not guaranteed to always have a point to each side
    @test_throws AssertionError nuresample(xin, yin, [2.5, 4.9, 5.0, 6.4],
    upsampling_function=rect_window(1 - eps(1.)),
    required_points_per_slice=1)
end

@testset "Conversion from UnitRange to StepRangeLen" begin
    @test nuresample(xin, yin, [4.2, 5.5]) == nuresample(1:12, yin, [4.2, 5.5])
    @test_throws MethodError nuresample(xin, yin) # check this dispatch isn't recursive
    @test_throws MethodError nuresample(xin, yin, [2.2, 3.1]; abc=12) # check this dispatch isn't recursive
end

@testset "StackOverflowError with wrong argument types" begin
    @test_broken nuresample(xin, yin, [3.2, 5.2], 1.3)
end


text:
@inline function generate_bc_eqs!(bceqs, s, boundaryvalfuncs, interiormap, boundary::AbstractTruncatingBoundary)
    args = params(depvar(boundary.u, s), s)
    indexmap = Dict([args[i]=>i for i in 1:length(args)])
    push!(bceqs, generate_bc_eqs(s, boundaryvalfuncs, boundary, interiormap, indexmap))
end

function generate_bc_eqs!(bceqs, s::DiscreteSpace{N}, boundaryvalfuncs, interiormap, boundary::PeriodicBoundary) where N
    # depvarbcmaps will dictate what to replace the variable terms with in the bcs
    # replace u(t,0) with u₁, etc
    u_, x_ = getvars(boundary)
    j = x2i(s, depvar(u_, s), x_)
    # * Assume that the periodic BC is of the simple form u(t,0) ~ u(t,1)
    Ioffset = unitindex(N, j)*(length(s, x_) - 1)
    local disc
    for u in s.ū
        if isequal(operation(u), operation(u_))
            disc =  s.discvars[u]
            break
        end
    end

    push!(bceqs, vec(map(edge(s, boundary, interiormap)) do II
        disc[II] ~ disc[II + Ioffset]
    end))

end

function generate_boundary_val_funcs(s, depvars, boundarymap, indexmap, derivweights)
    return mapreduce(vcat, values(boundarymap)) do boundaries
        map(reduce(vcat, values(boundaries))) do b
            if b isa PeriodicBoundary
                II -> []
            else
                II -> boundary_value_maps(II, s, b, derivweights, indexmap)
            end
        end
    end
end

function boundary_value_maps(II, s::DiscreteSpace{N,M,G}, boundary, derivweights, indexmap) where {N,M,G<:EdgeAlignedGrid}
    u_, x_ = getvars(boundary)

    ufunc(v, I, x) = s.discvars[v][I]

    # depvarbcmaps will dictate what to replace the variable terms with in the bcs
    # replace u(t,0) with u₁, etc

    u = depvar(u_, s)
    args = params(u, s)
    j = findfirst(isequal(x_), args)

    # We need to construct a new index in case the value at the boundary appears in an equation one dimension lower
    is = [II[indexmap[x]] for x in filter(!isequal(x_), args)]

    is = [is[1:j-1]..., idx(boundary, s), is[j:end]...]
    II = CartesianIndex(is...)

    # Shift depending on the boundary
    shift(::LowerBoundary) = zero(II)
    shift(::UpperBoundary) = unitindex(N, j)

    depvarderivbcmaps = [(Differential(x_)^d)(u_) => half_offset_centered_difference(derivweights.halfoffsetmap[1][Differential(x_)^d], II-shift(boundary), s, isperiodic(boundary), (j,x_), u, ufunc) for d in derivweights.orders[x_]]

    depvarbcmaps = [u_ => half_offset_centered_difference(derivweights.interpmap[x_], II-shift(boundary), s, isperiodic(boundary), (j,x_), u, ufunc)]

    return vcat(depvarderivbcmaps, depvarbcmaps)
end

function boundary_value_maps(II, s::DiscreteSpace{N,M,G}, boundary, derivweights, indexmap) where {N,M,G<:CenterAlignedGrid}
    u_, x_ = getvars(boundary)
    ufunc(v, I, x) = s.discvars[v][I]

    depvarderivbcmaps = []
    depvarbcmaps = []

    # * Assume that the BC is in terms of an explicit expression, not containing references to variables other than u_ at the boundary
    u = depvar(u_, s)
    args = params(u, s)
    j = findfirst(isequal(x_), args)

    # We need to construct a new index in case the value at the boundary appears in an equation one dimension lower
    is = [II[indexmap[x]] for x in filter(!isequal(x_), args)]

    is = [is[1:j-1]..., idx(boundary, s), is[j:end]...]
    II = CartesianIndex(is...)

    depvarderivbcmaps = [(Differential(x_)^d)(u_) => central_difference(derivweights.map[Differential(x_)^d], II, s,isperiodic(boundary), (x2i(s, u, x_), x_), u, ufunc) for d in derivweights.orders[x_]]
    # ? Does this need to be done for all variables at the boundary?
    depvarbcmaps = [u_ => s.discvars[u][II]]

    return vcat(depvarderivbcmaps, depvarbcmaps)
end

function generate_bc_eqs(s::DiscreteSpace{N,M,G}, boundaryvalfuncs, boundary::AbstractTruncatingBoundary, interiormap, indexmap) where {N, M, G}
    bc = boundary.eq
    return vec(map(edge(s, boundary, interiormap)) do II
        boundaryvalrules = mapreduce(f -> f(II), vcat, boundaryvalfuncs)
        vmaps = varmaps(s, boundary.depvars, II, indexmap)
        varrules = axiesvals(s, depvar(boundary.u, s), boundary.x, II)
        rules = vcat(boundaryvalrules, vmaps, varrules)

        substitute(bc.lhs, rules) ~ substitute(bc.rhs, rules)
    end)
end

#TODO: Benchmark and optimize this

@inline function generate_corner_eqs!(bceqs, s, interiormap, N, u)
    interior = interiormap.I[interiormap.pde[u]]
    sd(i, j) = selectdim(interior, j, i)
    domain = setdiff(s.Igrid[u], interior)
    II1 = unitindices(N)
    for j in 1:N
        I1 = II1[j]
        edge = sd(1, j)
        offset = edge[1][j]-1
        for k in 1:offset
            setdiff!(domain, vec(copy(edge).-[I1*k]))
        end
        edge = sd(size(interior, j), j)
        offset = size(s.discvars[u], j) - size(interior, j)
        for k in 1:offset
            setdiff!(domain, vec(copy(edge).+[I1*k]))
        end
    end
    push!(bceqs, s.discvars[u][domain] .~ 0)
end

"""
Create a vector containing indices of the corners of the domain.
"""
@inline function findcorners(s::DiscreteSpace, lower, upper, u) where {M}
    args = remove(arguments(u), s.time)
    if any(lower.==0) && any(upper.==0)
        return CartesianIndex{2}[]
    end
    return reduce(vcat, vec.(map(0 : 3) do n
        dig = digits(n, base = 2, pad = 2)
        CartesianIndices(Tuple(map(enumerate(dig)) do (i, b)
            x = args[i]
            if b == 1
                1:lower[i]
            elseif b == 0
                length(s, x)-upper[i]+1:length(s, x)
            end
        end))
    end))
end

@inline function generate_corner_eqs!(bceqs, s, interiormap, pde)
    u = interiormap.var[pde]
    N = ndims(u, s)
    if N <= 1
        return
    elseif N == 2
        Icorners = findcorners(s, interiormap.lower[pde], interiormap.upper[pde], u)
        push!(bceqs, s.discvars[u][Icorners] .~ 0)
    else
        generate_corner_eqs!(bceqs, s, interiormap, N, u)
    end
end

ext: jl


text:
# TODO: Add methods for other matrix-creation functions, such as zeros() and rand().

StaticArrays.SOneTo(::StaticNumbers.StaticOneTo{L}) where {L} = StaticArrays.SOneTo{L}()
StaticNumbers.StaticOneTo(::StaticArrays.SOneTo{L}) where {L} = StaticNumbers.StaticOneTo{L}()
StaticNumbers.staticlength(r::StaticArrays.SOneTo) = StaticNumbers.StaticOneTo(r)
StaticNumbers.static(r::StaticArrays.SOneTo) = StaticNumbers.StaticOneTo(r)

# TODO: Maybe just SOneTo from StaticArrays instead of defining StaticOneTo?

@inline StaticArrays.SUnitRange(r::StaticNumbers.LengthUnitRange) = StaticArrays.SUnitRange{Int(first(r)), Int(length(r))}()
StaticNumbers.LengthUnitRange(::StaticArrays.SUnitRange{Start, L}) where {Start, L} = StaticNumbers.LengthUnitRange(StaticNumbers.static(Start-1), StaticNumbers.static(L))
StaticNumbers.staticlength(r::StaticArrays.SUnitRange) = StaticNumbers.LengthUnitRange(r)
StaticNumbers.static(r::StaticArrays.SUnitRange) = StaticNumbers.LengthUnitRange(r)

StaticArrays.Size(::LengthRange{T,Z,S,StaticInteger{L}}) where {T,Z,S,L} = StaticArrays.Size(Int(L))
StaticArrays.Size(s::Tuple{StaticInteger, Vararg{StaticInteger}}) = StaticArrays.Size(Int.(s))
StaticArrays.Size(s::StaticInteger, ss::StaticInteger...) = StaticArrays.Size((s, ss...))
StaticNumbers.static(::StaticArrays.Size{S}) where {S} = StaticNumbers.static.(S)

for AT in (Array, StaticArrays.MArray), RT in (LengthStepRange{<:Integer,Z,S,<:StaticInteger} where {Z,S}, LengthUnitRange{<:Integer,Z,<:StaticInteger} where {Z})
    Base.getindex(A::AT, r::RT) = StaticArrays.MVector(ntuple(i -> A[r[i]], length(r)))
end

for AT in (StaticArrays.SArray,), RT in (LengthStepRange{<:Integer,Z,S,<:StaticInteger} where {Z,S}, LengthUnitRange{<:Integer,Z,<:StaticInteger} where {Z})
    Base.getindex(A::AT, r::RT) = StaticArrays.SVector(ntuple(i -> A[r[i]], length(r)))
end

# Allow "end" to become static when indexing into static arrays
@inline StaticNumbers.maybe_static(::typeof(lastindex), A::StaticArrays.StaticArray) = StaticNumbers.static(lastindex(A))
@inline StaticNumbers.maybe_static(::typeof(lastindex), A::StaticArrays.StaticArray, d::StaticInteger) = StaticNumbers.static(lastindex(A, i))

@inline StaticNumbers.maybe_static(::typeof(size), A::StaticArrays.StaticArray) = map(static, size(A))
@inline StaticNumbers.maybe_static(::typeof(size), A::StaticArrays.StaticArray, d::Integer) = static(size(A, d))

@inline (::Type{T})(g::Base.Generator{<:LengthRange}) where {T<:StaticVector} = T(Tuple(g))
@inline (::Type{T})(iter::StaticNumbers.LengthRange) where {T<:StaticVector} = T(Tuple(iter))

@inline function (::Type{T})(g::Base.Generator{<:Base.Iterators.ProductIterator{<:Tuple{StaticNumbers.StaticLengthRange, Vararg{StaticNumbers.StaticLengthRange}}}}) where {T<:StaticArrays.StaticArray}
    sz = StaticArrays.Size(size(g))
    data = Tuple(g)
    ST = similar_type(T, eltype(data), sz)
    ST <: T || error("Generator yields the wrong type of static array.")
    ST(data)
end

Base.@propagate_inbounds function Base.convert(::Type{SA}, a::SubArray{T,N,<:AbstractArray{T},<:Tuple{StaticNumbers.StaticLengthRange, Vararg{StaticNumbers.StaticLengthRange}}}) where {SA <: StaticArrays.StaticArray, T, N}
    Z = similar_type(SA, T, StaticArrays.Size(size(a)))
    Z <: SA || error("Input array yields the wrong type of static array.")
    Z(Tuple(a))
end

# # TODO: Interface to StaticArrays for multi-dimensional indexing


text:
using ChaosTools
using Test, StaticArrays, OrdinaryDiffEq

println("\nTesting orbit diagrams...")

@testset "Orbit Diagram" begin
    @testset "Discrete 1D" begin
        ds = Systems.logistic()
        i = 1
        parameter = 1
        pvalues = 2:0.01:4
        ics = [rand() for m in 1:10]
        n = 50
        Ttr = 5000
        output = orbitdiagram(ds, i, parameter, pvalues; n = n, Ttr = Ttr)
        @test output[1][1] ≈ 0.5
        @test output[2][1] == 0.5024875621890547
        @test output[end][1] != output[end][2] != output[end][3]
    end

    @testset "Discrete 2D" begin
        ds = Systems.standardmap()
        i = 2
        parameter = 1
        pvalues = 0:0.005:2
        ics = [0.001rand(2) for m in 1:10]
        n = 50
        Ttr = 5000
        output = orbitdiagram(ds, i, parameter, pvalues; n = n, Ttr = Ttr)
        @test length(output[1]) == n
    end

    @testset "Discrete IIP" begin
        ds = Systems.henon_iip()
        i = 2
        parameter = 1
        pvalues = [1.4, 1.5]
        ics = [rand(2), rand(2)]
        n = 50
        Ttr = 5000
        output = orbitdiagram(ds, i, parameter, pvalues; n = n, Ttr = Ttr)
        @test length(output[1]) == n
    end
end

@testset "Poincare SOS" begin
    @testset "Henon Helies" begin
        ds = Systems.henonheiles([0, .483000, .278980390, 0] )
        psos = poincaresos(ds, (2, 0.0), 1000.0)
        xcross = psos[:, 2]
        @test length(xcross) > 1
        for x in xcross
            @test abs(x) < 1e-3
        end

        # @inline Vhh(q1, q2) = 1//2 * (q1^2 + q2^2 + 2q1^2 * q2 - 2//3 * q2^3)
        # @inline Thh(p1, p2) = 1//2 * (p1^2 + p2^2)
        # @inline Hhh(q1, q2, p1, p2) = Thh(p1, p2) + Vhh(q1, q2)
        # @inline Hhh(u::AbstractVector) = Hhh(u...)
        #
        # E = [Hhh(p) for p in psos]
        #
        # @test std(E) < 1e-10
        # @test maximum(@. E - E[1]) < 1e-10
    end
    @testset "Gissinger crazy plane" begin
        gis = Systems.gissinger([2.32865, 2.02514, 1.98312])
        # Define appropriate hyperplane for gissinger system
        ν = 0.1
        Γ = 0.9 # default parameters of the system

        # I want hyperperplane defined by these two points:
        Np(μ) = SVector{3}(sqrt(ν + Γ*sqrt(ν/μ)), -sqrt(μ + Γ*sqrt(μ/ν)), -sqrt(μ*ν))
        Nm(μ) = SVector{3}(-sqrt(ν + Γ*sqrt(ν/μ)), sqrt(μ + Γ*sqrt(μ/ν)), -sqrt(μ*ν))

        # Create hyperplane using normal vector to vector connecting points:
        gis_plane(μ) = (d = (Np(μ) - Nm(μ)); [d[2], -d[1], 0, 0])

        μ = 0.12
        set_parameter!(gis, 1, μ)
        psos = poincaresos(gis, gis_plane(μ), 10000.0, Ttr = 200.0, direction = -1)
        @test length(psos) > 1
        @test generalized_dim(2, psos) < 1
    end
    @testset "beginning on the plane" begin

        @inbounds @inline function ż(z, p, t)
            A, B, D = p
            p₀, p₂ = z[1:2]
            q₀, q₂ = z[3:4]

            return SVector{4}(
                -A * q₀ - 3 * B / √2 * (q₂^2 - q₀^2) - D * q₀ * (q₀^2 + q₂^2),
                -q₂ * (A + 3 * √2 * B * q₀ + D * (q₀^2 + q₂^2)),
                A * p₀,
                A * p₂
            )
        end
        z0 = SVector{4}(2.3499921651423565, -9.801351029039825, 0.0, -2.7230316965872268)

        ds = ContinuousDynamicalSystem(ż, z0, (A=1, B=0.55, D=0.4))
        psos = poincaresos(ds, (3, 0), 10., direction=1, idxs=[2,4])
        @test length(psos) > 0
        @test size(psos)[2] == 2
    end
end

@testset "Produce OD" begin
  @testset "Shinriki" begin
    ds = Systems.shinriki([-2, 0, 0.2])

    pvalues = range(19,stop=22,length=11)
    parameter = 1
    i = 1
    j = 2
    tf = 200.0

    de = Dict(:abstol=>1e-9, :reltol => 1e-9)
    output = produce_orbitdiagram(ds, (j, 0.0), i, parameter, pvalues; tfinal = tf,
    Ttr = 100.0, printparams = false, de...)
    @test length(output) == length(pvalues)

    v = round.(output[1], digits = 4)
    s = collect(Set(v))
    @test length(s) == 1
    @test s[1] == -0.856

    v = round.(output[4], digits = 4)
    s = Set(v)
    @test length(s) == 2
    @test s == Set([-0.376, -1.2887])

  end
end

# # This test offers nothing to the Suite (does pass though)
# @testset "Stroboscopic" begin
#   ds = Systems.duffing(β = -1, ω = 1, f = 0.3)
#   a = trajectory(ds, 100000.0, dt = 2π)
#   D = information_dim(a)
#   @test 1.25 < D < 1.5
# end


text:
using Revise
using AdFem
using PyPlot
using ADCME

m = 100
n = 100
h = 1/n
Δt = 0.1
ρ = 1

function step3(u0, dp)
    A = constant(compute_fem_mass_matrix(m, n, h))
    bd = bcnode("all", m, n, h)

    grad_dp_fem = - compute_interaction_term(dp, m, n, h)
    b = A * u0 - Δt/ρ * grad_dp_fem

    A, _ = fem_impose_Dirichlet_boundary_condition(A, bd, m, n, h) # 2D boundary condition
    bd = [bd; (m+1)*(n+1) .+ bd]
    b = scatter_update(b, bd, zeros(length(bd)))
    sol = A\b
    return sol
end

xy = fvm_nodes(m, n, h)
x, y = xy[:,1], xy[:,2]
dp = @. x^2*(1-x)^2*y^2*(1 -y)^2

xy = fem_nodes(m, n, h)
x, y = xy[:,1], xy[:,2]
u0 = [@. (x*(1-x)*y*(1-y)); @. (x*(1-x)*y*(1-y))]

grad_dp_x = x.^2 .*y.^2 .*(1 .- y).^2 .*(2*x .- 2) + 2*x.*y.^2 .*(1 .- x).^2 .*(1 .- y).^2
grad_dp_y = x.^2 .*y.^2 .*(1 .- x).^2 .*(2*y .- 2) + 2*x.^2 .*y.*(1 .- x).^2 .*(1 .- y).^2

# u1 = [@. (x*(1-x)*y*(1-y)- Δt * (1-2x)*y*(1-y) ); @. (x*(1-x)*y*(1-y) - Δt * (1-2y)*x*(1-x) )]
u1_x = u0[1:(m+1)*(n+1)]- Δt / ρ * grad_dp_x
u1_y = u0[(m+1)*(n+1)+1:end] - Δt / ρ *grad_dp_y

# A = constant(compute_fem_mass_matrix(m, n, h))
# bd = bcnode("all", m, n, h)

# grad_dp_fem = - compute_interaction_term(dp, m, n, h)
# b = A * u0 - Δt/ρ * grad_dp_fem

# A, _ = fem_impose_Dirichlet_boundary_condition(A, bd, m, n, h) # 2D boundary condition
# bd = [bd; (m+1)*(n+1) .+ bd]
# b = scatter_update(b, bd, zeros(length(bd)))

# sol = A\b

sol = step3(u0, dp)

sess = Session(); init(sess)
sol = run(sess, sol)
# u1 = run(sess, u1)

# exact = u1(x, y)

figure(figsize=(10,4))
subplot(121)
visualize_scalar_on_fem_points(sol[1:(m+1)*(n+1)]-u1_x, m, n, h)
subplot(122)
visualize_scalar_on_fem_points(u1_x, m, n, h)

figure(figsize=(10,4))
subplot(121)
visualize_scalar_on_fem_points(sol[(m+1)*(n+1)+1:end]-u1_y, m, n, h)
subplot(122)
visualize_scalar_on_fem_points(u1_y, m, n, h)


text:
s = """[i for i in 1:2]
"""
msgs = lintstr(s)
@test isempty(msgs)

s = if VERSION ≥ v"0.5.0-"
    """
    Dict(j => j*j for j in 1:2)
    """
else
    """
    [j => j*j for j in 1:2]
    """
end
msgs = lintstr(s)
@test isempty(msgs)

s = if VERSION ≥ v"0.5.0-"
    """
    Dict{Int,Int}(y2 => y2*y2 for y2 in 1:2)
    """
else
    """
    (Int=>Int)[y2 => y2*y2 for y2 in 1:2]
    """
end
msgs = lintstr(s)
@test isempty(msgs)


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta

const GLMRegressor_PostEvaluationTransform = (;[
    Symbol("NoTransform") => Int32(0),
    Symbol("Logit") => Int32(1),
    Symbol("Probit") => Int32(2),
]...)

mutable struct GLMRegressor_DoubleArray <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function GLMRegressor_DoubleArray(; kwargs...)
        obj = new(meta(GLMRegressor_DoubleArray), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct GLMRegressor_DoubleArray
const __meta_GLMRegressor_DoubleArray = Ref{ProtoMeta}()
function meta(::Type{GLMRegressor_DoubleArray})
    ProtoBuf.metalock() do
        if !isassigned(__meta_GLMRegressor_DoubleArray)
            __meta_GLMRegressor_DoubleArray[] = target = ProtoMeta(GLMRegressor_DoubleArray)
            pack = Symbol[:value]
            allflds = Pair{Symbol,Union{Type,String}}[:value => Base.Vector{Float64}]
            meta(target, GLMRegressor_DoubleArray, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, pack, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_GLMRegressor_DoubleArray[]
    end
end
function Base.getproperty(obj::GLMRegressor_DoubleArray, name::Symbol)
    if name === :value
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{Float64}
    else
        getfield(obj, name)
    end
end

mutable struct GLMRegressor <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function GLMRegressor(; kwargs...)
        obj = new(meta(GLMRegressor), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct GLMRegressor
const __meta_GLMRegressor = Ref{ProtoMeta}()
function meta(::Type{GLMRegressor})
    ProtoBuf.metalock() do
        if !isassigned(__meta_GLMRegressor)
            __meta_GLMRegressor[] = target = ProtoMeta(GLMRegressor)
            pack = Symbol[:offset]
            allflds = Pair{Symbol,Union{Type,String}}[:weights => Base.Vector{GLMRegressor_DoubleArray}, :offset => Base.Vector{Float64}, :postEvaluationTransform => Int32]
            meta(target, GLMRegressor, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, pack, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_GLMRegressor[]
    end
end
function Base.getproperty(obj::GLMRegressor, name::Symbol)
    if name === :weights
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{GLMRegressor_DoubleArray}
    elseif name === :offset
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{Float64}
    elseif name === :postEvaluationTransform
        return (obj.__protobuf_jl_internal_values[name])::Int32
    else
        getfield(obj, name)
    end
end

export GLMRegressor_PostEvaluationTransform, GLMRegressor_DoubleArray, GLMRegressor


text:
struct Solver
    path::String
    args::String
    source::String
end

function run(solver::Solver, folder::String)
    binary = solver.path
    args = solver.args
    source = solver.source

    old_pwd = pwd()
    cd(folder)

    out = joinpath(folder, basename(binary) * "UncertaintyQuantification.out")
    err = joinpath(folder, basename(binary) * "UncertaintyQuantification.err")

    if !isempty(args)
        Base.run(pipeline(`$binary $args $source`; stdout=out, stderr=err))
    else
        Base.run(pipeline(`$binary $source`; stdout=out, stderr=err))
    end

    return cd(old_pwd)
end

max_stars_repo_name: FriesischScott/julia-uq

max_issues_repo_name: FriesischScott/julia-uq

max_forks_repo_name: FriesischScott/julia-uq


text:
generic_sphere(ref_dirs) = ref_dirs ./ ([norm(v) for v in ref_dirs])

function DTLZ_hypersphere!(fx, x; α = 1)
    m = length(fx)
    for i = 1:m
        if i < m
            fx[i] *= prod(cos.((π*0.5) * x[1:m-i] .^ α ))
        end
        if i > 1
            fx[i] *= sin((π*0.5) * x[1+m - i] .^ α)
        end
    end

    return fx
end

function DTLZ_hyperplane!(fx, x)
    m = length(fx)
    for i in 1:m
        fx[i] *= prod(x[1:m-i])
        if i > 1
            fx[i] *= 1 - x[1+m - i]
        end
    end
    fx
end

function DTLZ_g1(x, m)
    y = view(x, m:length(x)) .- 0.5
    return 100*( length(y) + sum( y.^2 - cos.(20π*y) ))
end

function DTLZ_g2(x, m)
    return sum( (view(x, m:length(x)) .- 0.5).^2 )
end

function DTLZ1_f(x, m = 3)
    g = DTLZ_g1(x, m)
    D = length(x)

    fx = fill(0.5*(1 + g), m)
    DTLZ_hyperplane!(fx, x)

    return fx, [0.0], [0.0]
end

"""
    DTLZ1(m = 3, ref_dirs = gen_ref_dirs(m, 12))

DTLZ1 returns `(f::function, bounds::Matrix{Float64}, pareto_set::Array{xFgh_indiv})`
where `f` is the objective function and `pareto_set` is an array with optimal Pareto solutions
with `n_solutions`.

### Parameters
- `m` number of objective functions
- `ref_dirs` number of Pareto solutions (default: Das and Dennis' method).

Main properties:
- convex
- multifrontal
"""
function DTLZ1(m=3, ref_dirs = gen_ref_dirs(m, 12))

    D = 10 + m - 1

    bounds = Array([zeros(D) ones(D)]')

    pf = 0.5ref_dirs
    pareto_set = [ generateChild(zeros(0), (fx, [0.0], [0.0])) for fx in pf ]

    return DTLZ1_f, bounds, pareto_set
end

function DTLZ2_f(x, m = 3)
    g = DTLZ_g2(x, m)
    fx = fill(1.0 + g, m)
    DTLZ_hypersphere!(fx, x)

    return fx, [0.0], [0.0]
end

"""
    DTLZ2(m = 3, ref_dirs = gen_ref_dirs(m, 12))

DTLZ2 returns `(f::function, bounds::Matrix{Float64}, pareto_set::Array{xFgh_indiv})`
where `f` is the objective function and `pareto_set` is an array with optimal Pareto solutions
with `n_solutions`.

### Parameters
- `m` number of objective functions
- `ref_dirs` number of Pareto solutions (default: Das and Dennis' method).

Main properties:
- nonconvex
- unifrontal
"""
function DTLZ2(m=3, ref_dirs = gen_ref_dirs(m, 12))

    D = 10 + m - 1

    bounds = Array([zeros(D) ones(D)]')

    pf = generic_sphere(ref_dirs)
    pareto_set = [ generateChild(zeros(0), (fx, [0.0], [0.0])) for fx in pf ]

    return DTLZ2_f, bounds, pareto_set
end

function DTLZ3_f(x,m=3)
    g = DTLZ_g1(x, m)
    fx = fill(1 + g, m)
    DTLZ_hypersphere!(fx, x)

    return fx, [0.0], [0.0]
end

"""
    DTLZ3(m = 3, ref_dirs = gen_ref_dirs(m, 12))

DTLZ3 returns `(f::function, bounds::Matrix{Float64}, pareto_set::Array{xFgh_indiv})`
where `f` is the objective function and `pareto_set` is an array with optimal Pareto solutions
with `n_solutions`.

### Parameters
- `m` number of objective functions
- `ref_dirs` number of Pareto solutions (default: Das and Dennis' method).

Main properties:
- nonconvex
- multifrontal
"""
function DTLZ3(m = 3, ref_dirs = gen_ref_dirs(m, 12))

    D = 10 + m - 1

    bounds = Array([zeros(D) ones(D)]')

    pf = generic_sphere(ref_dirs)
    pareto_set = [ generateChild(zeros(0), (fx, [0.0], [0.0])) for fx in pf ]

    return DTLZ3_f, bounds, pareto_set
end

function DTLZ4_f(x,m=3)
    g = DTLZ_g2(x, m)
    fx = fill(1.0 + g, m)
    DTLZ_hypersphere!(fx, x; α = 100)

    return fx, [0.0], [0.0]
end

"""
    DTLZ4(m = 3, ref_dirs = gen_ref_dirs(m, 12))

DTLZ4 returns `(f::function, bounds::Matrix{Float64}, pareto_set::Array{xFgh_indiv})`
where `f` is the objective function and `pareto_set` is an array with optimal Pareto solutions
with `n_solutions`.

### Parameters
- `m` number of objective functions
- `ref_dirs` number of Pareto solutions (default: Das and Dennis' method).

Main properties:
- nonconvex
- unifrontal
"""
function DTLZ4(m = 3, ref_dirs = gen_ref_dirs(m, 12))

    D = 10 + m - 1

    bounds = Array([zeros(D) ones(D)]')

    pf = generic_sphere(ref_dirs)
    pareto_set = [ generateChild(zeros(0), (fx, [0.0], [0.0])) for fx in pf ]

    return DTLZ4_f, bounds, pareto_set
end

function DTLZ5_f(x,m = 3)
    g = DTLZ_g2(x, m)
    fθ = fill(1.0 + g, m)
    θ = @. 1 / (2*(1 + g)) * (1 + 2g * x[1:m-1] )
    θ[1] = x[1]

    DTLZ_hypersphere!(fθ, θ)

    return fθ, [0.0], [0.0]
end

"""
    DTLZ5(m = 3)

DTLZ5 returns `(f::function, bounds::Matrix{Float64}, pareto_set::Array{xFgh_indiv})`
where `f` is the objective function and `pareto_set` is an array with optimal Pareto solutions
with `n_solutions`.

### Parameters
- `m` number of objective functions

"""
function DTLZ5(m = 3)

    D = 10 + m - 1

    bounds = Array([zeros(D) ones(D)]')

    if m == 3
        ref_dirs = gen_ref_dirs(2, 100)
        X = fill(0.5, length(ref_dirs), D)

        for i in eachindex(ref_dirs)
            X[i,1:2] = ref_dirs[i]
        end
        pareto_set = [ generateChild(X[i,:], DTLZ5_f(X[i,:])) for i in eachindex(ref_dirs) ]

    else
        pareto_set = []
    end

    return DTLZ5_f, bounds, pareto_set
end

function DTLZ6_f(x,m=3)
    g = sum(x[m:end] .^ 0.1)
    fθ = fill(1.0 + g, m)
    θ = @. 1 / (2*(1 + g)) * (1 + 2g * x[1:m-1] )
    θ[1] = x[1]

    DTLZ_hypersphere!(fθ, θ)

    return fθ, [0.0], [0.0]
end

"""
    DTLZ6(m = 3)

DTLZ6 returns `(f::function, bounds::Matrix{Float64}, pareto_set::Array{xFgh_indiv})`
where `f` is the objective function and `pareto_set` is an array with optimal Pareto solutions
with `n_solutions`.

### Parameters
- `m` number of objective functions

"""
function DTLZ6(m = 3)

    D = 10

    bounds = Array([zeros(D) ones(D)]')

    if m == 3
        ref_dirs = gen_ref_dirs(2, 100)
        X = fill(0.0, length(ref_dirs), D)

        for i in eachindex(ref_dirs)
            X[i,1:2] = ref_dirs[i]
        end
        pareto_set = [ generateChild(X[i,:], DTLZ6_f(X[i,:])) for i in eachindex(ref_dirs) ]

    else
        pareto_set = []
    end

    return DTLZ6_f, bounds, pareto_set
end

#=
####################################################################################
####################################################################################
####################################################################################
#                Constrained DTLZ
####################################################################################
####################################################################################
####################################################################################
=#

function C1_DTLZ1_f(x, m = 3)
    g = DTLZ_g1(x, m)
    D = length(x)

    fx = fill(0.5*(1 + g), m)
    DTLZ_hyperplane!(fx, x)

    c = fx[end]/0.6 + sum(fx[1:end-1]/0.5) - 1;
    return fx, [c], [0.0]
end

"""
    C1_DTLZ1(m = 3, ref_dirs = gen_ref_dirs(m, 12))

C1_DTLZ1 returns `(f::function, bounds::Matrix{Float64}, pareto_set::Array{xFgh_indiv})`
where `f` is the objective function and `pareto_set` is an array with optimal Pareto solutions
with `n_solutions`.

### Parameters
- `m` number of objective functions
- `ref_dirs` number of Pareto solutions (default: Das and Dennis' method).

Main properties:
- convex
- multifrontal
- constraints type 1
"""
function C1_DTLZ1(m=3, ref_dirs = gen_ref_dirs(m, 12))

    D = 10 + m - 1

    bounds = Array([zeros(D) ones(D)]')

    pf = 0.5ref_dirs
    pareto_set = [ generateChild(zeros(0), (fx, [0.0], [0.0])) for fx in pf ]

    return C1_DTLZ1_f, bounds, pareto_set
end

function C1_DTLZ3_f(x,m=3)
    g = DTLZ_g1(x, m)
    fx = fill(1 + g, m)
    DTLZ_hypersphere!(fx, x)

    if m == 2
        r = 6
    elseif m <= 3
        r = 9
    elseif m <= 8
        r = 12.5
    else
        r = 15
    end
    c = -(sum(fx .^ 2) - 16) * (sum(fx .^ 2) - r^2);

    return fx, [c], [0.0]
end

"""
    C1_DTLZ3(m = 3, ref_dirs = gen_ref_dirs(m, 12))

DTLZ3 returns `(f::function, bounds::Matrix{Float64}, pareto_set::Array{xFgh_indiv})`
where `f` is the objective function and `pareto_set` is an array with optimal Pareto solutions
with `n_solutions`.

### Parameters
- `m` number of objective functions
- `ref_dirs` number of Pareto solutions (default: Das and Dennis' method).

Main properties:
- nonconvex
- multifrontal
- constraints type 1
"""
function C1_DTLZ3(m = 3, ref_dirs = gen_ref_dirs(m, 12))

    D = 10 + m - 1

    bounds = Array([zeros(D) ones(D)]')

    pf = generic_sphere(ref_dirs)
    pareto_set = [ generateChild(zeros(0), (fx, [0.0], [0.0])) for fx in pf ]

    return C1_DTLZ3_f, bounds, pareto_set
end

function C2_DTLZ2_f(x, m = 3)
    g = DTLZ_g2(x, m)
    fx = fill(1.0 + g, m)
    DTLZ_hypersphere!(fx, x)

    if m == 3
        r = 0.4
    else
        r = 0.5
    end

    p1 = maximum( i -> (fx[i] .- 1).^2 + sum(fx[1:m].^2 .- fx[i]^2) - r^2,1:m)
    p2 = sum((fx .- 1 / sqrt(m)).^2) - r^2
    return fx, [-max(p1, p2)], [0.0]
end

"""
    C2_DTLZ2(m = 3, ref_dirs = gen_ref_dirs(m, 12))

DTLZ2 returns `(f::function, bounds::Matrix{Float64}, pareto_set::Array{xFgh_indiv})`
where `f` is the objective function and `pareto_set` is an array with optimal Pareto solutions
with `n_solutions`.

### Parameters
- `m` number of objective functions
- `ref_dirs` number of Pareto solutions (default: Das and Dennis' method).

Main properties:
- nonconvex
- unifrontal
- contraints type 2
"""
function C2_DTLZ2(m=3, ref_dirs = gen_ref_dirs(m, 12))

    D = 10 + m - 1

    bounds = Array([zeros(D) ones(D)]')

    pf = generic_sphere(ref_dirs)
    pareto_set = [ generateChild(zeros(0), (fx, [0.0], [0.0])) for fx in pf ]

    return C2_DTLZ2_f, bounds, pareto_set
end

function C3_DTLZ1_f(x, m = 3)
    g = DTLZ_g1(x, m)
    D = length(x)

    fx = fill(0.5*(1 + g), m)
    DTLZ_hyperplane!(fx, x)

    c = [sum(fx[j] .+ fx/0.5 .- fx[j]/0.5) - 1 for j in 1:m]
    return fx, -c, [0.0]
end

#=
"""
    C3_DTLZ1(m = 3, ref_dirs = gen_ref_dirs(m, 12))

C3_DTLZ1 returns `(f::function, bounds::Matrix{Float64}, pareto_set::Array{xFgh_indiv})`
where `f` is the objective function and `pareto_set` is an array with optimal Pareto solutions
with `n_solutions`.

### Parameters
- `m` number of objective functions
- `ref_dirs` number of Pareto solutions (default: Das and Dennis' method).

Main properties:
- convex
- multifrontal
- constraints type 3
"""
function C3_DTLZ1(m=3, ref_dirs = gen_ref_dirs(m, 12))

    D = m + 4

    bounds = Array([zeros(D) ones(D)]')

    @warn "C3_DTLZ1 is under development"

    pf = ref_dirs./([(0.5r -6*maximum(r)) for r in ref_dirs])
    pareto_set = [ generateChild(zeros(0), (fx, [0.0], [0.0])) for fx in pf ]

    return C3_DTLZ1_f, bounds, pareto_set
end
=#

function C3_DTLZ4_f(x,m=3)
    g = DTLZ_g2(x, m)
    fx = fill(1.0 + g, m)
    DTLZ_hypersphere!(fx, x; α = 100)

    c = [fx[j]^2/4  + sum(fx.^2 .- fx[j]^2) - 1 for j in 1:m]
    return fx, -c, [0.0]
end

"""
    C3_DTLZ4(m = 3, ref_dirs = gen_ref_dirs(m, 12))

C3_DTLZ4 returns `(f::function, bounds::Matrix{Float64}, pareto_set::Array{xFgh_indiv})`
where `f` is the objective function and `pareto_set` is an array with optimal Pareto solutions
with `n_solutions`.

### Parameters
- `m` number of objective functions
- `ref_dirs` number of Pareto solutions (default: Das and Dennis' method).

Main properties:
- nonconvex
- unifrontal
- constraints type 3
"""
function C3_DTLZ4(m = 3, ref_dirs = gen_ref_dirs(m, 12))

    D = m + 4

    bounds = Array([zeros(D) ones(D)]')

    pf = ref_dirs./([sqrt.(sum(r.^2)-3/4*maximum(r.^2)) for r in ref_dirs])

    pareto_set = [ generateChild(zeros(0), (fx, [0.0], [0.0])) for fx in pf ]

    return C3_DTLZ4_f, bounds, pareto_set
end


text:
"""
# Preconditions.jl

	function f(x)
		@require x > 0
	end

From [HTTP.jl/src/debug.jl](https://git.io/JWkNl)
"""
module Preconditions

using ReadmeDocs

const DEBUG_LEVEL = Ref(0)

export @require, @ensure

function method_name(bt)
    for f in bt
        for i in StackTraces.lookup(f)
            if !i.from_c &&
               i.linfo != nothing &&
               i.func != :precondition_error &&
               i.func != :postcondition_error &&
               i.func != :backtrace

               return i.func
            end
        end
    end
    return "unknown method"
end

@noinline function precondition_error(msg::String)
    msg = string(method_name(backtrace()), " requires ", msg)
    return ArgumentError(msg)
end

README"""
    @require precondition [message]

Throw `ArgumentError` if `precondition` is false.
"""
macro require(condition, msg = string(condition))
    esc(:($condition || throw(Preconditions.precondition_error($msg))))
end

@noinline function postcondition_error(msg::String, ls="", l="", rs="", r="")
    @nospecialize
    msg = string(method_name(backtrace()), " failed to ensure ", msg)
    if ls != ""
        msg = string(msg, "\n", ls, " = ", sprint(show, l),
                          "\n", rs, " = ", sprint(show, r))
    end
    return AssertionError(msg)
end

# Copied from stdlib/Test/src/Test.jl:get_test_result()
iscondition(ex) = isa(ex, Expr) &&
                  ex.head == :call &&
                  length(ex.args) == 3 &&
                  first(string(ex.args[1])) != '.' &&
                  (!isa(ex.args[2], Expr) || ex.args[2].head != :...) &&
                  (!isa(ex.args[3], Expr) || ex.args[3].head != :...) &&
                  (ex.args[1] === :(==) ||
                       Base.operator_precedence(ex.args[1]) ==
                           Base.operator_precedence(:(==)))

README"""
    @ensure postcondition [message]

Throw `ArgumentError` if `postcondition` is false.
"""
macro ensure(condition, msg = string(condition))

    if DEBUG_LEVEL[] < 0
        return :()
    end

    if iscondition(condition)
        l,r = condition.args[2], condition.args[3]
        ls, rs = string(l), string(r)
        return esc(quote
            if ! $condition
                # FIXME double-execution of condition l and r!
                throw(Preconditions.postcondition_error($msg,
                                                        $ls, $l, $rs, $r))
            end
        end)
    end

    esc(:(if ! $condition throw(Preconditions.postcondition_error($msg)) end))
end

end # module


text:
module TestRR

using Compat
using Compat.Test
using IntervalArithmetic
using Clp
using EAGODomainReduction

type t_solver
    LP_solver
end

type test_opt
    numVar
    numConstr
    f
    g
    solver::t_solver
    gL
    gU
    gL_loc
    gU_loc
end

f(x) = x[1]^2+x[2]^2
g(x) = [x[1]+x[2];
        x[1]-x[2]]

in_opt = t_solver(ClpSolver())
opt = test_opt(2,2,f,g,in_opt,[2 -1],[3 0],[1,2],[1,2])
X = [Interval(0.0,2.0) for i=1:2]
STD_Linear_RR!(X,[opt],7.0)

@test 0.5-1E-4 <= X[1].lo <= 0.5+1E-4
@test 1.25-1E-4 <= X[1].hi <= 1.25+1E-4
@test 1.0-1E-4 <= X[2].lo <= 1.0+1E-4
@test 1.75-1E-4 <= X[2].hi <= 1.75+1E-4

end


text:
using Symbolics

function symbolic_kinetic_energy(
    u,
    w,
    θ′,
    t;
    θ_t = θ_t,
    θ_k = θ_k,
    γ_OA = γ_OA,
)

    D = Differential(t)
    T = kinetic_energy(u, w, θ′, θ_t, D(θ_t), θ_k, D(θ_k), γ_OA, D(γ_OA))

    Du = Differential(u)
    Dw = Differential(w)
    Dθ′ = Differential(θ′)

    dT_du = Du(T)
    dT_dw = Dw(T)
    dT_dθ′ = -simplify(expand_derivatives(Dθ′(T)))

    # LHS of the equations
    L1 = D(dT_du) - dT_dw * θ′
    L2 = D(dT_dw) + dT_du * θ′
    L3 = D(dT_dθ′) + dT_du * w - dT_dw * u

    [L1, L2, L3]
end


text:
export Events, networkactivity, nodeactivity, cet_ncet, fixevents,
fixevents!, widenevents!, snapshot, snapshots, mintime, maxtime,
meannodeactivity, numevents, mergeevents, mergeevents!, flatten

# Functions for working with dynamic networks
import Base: zero, ==

"""
An dynamic edge in a dynamic network is defined as a set of events. An
event is an interaction between two nodes from start time t0 to stop
time t1, t0 <= t1, and is represent as a tuple (t0,t1). The set of
events is represented as a Vector, sorted wrt start time t0, and if
two start times are equal, then stop time t1.  Note: Unless specified,
it is assumed that two events in an event set will not overlap,
i.e. there exists no two tuples (t0,t1) and (s0,s1)
s.t. intersect([t0,t1],[s0,s1]) != null

A dynamic network can be represented using a
SparseMatrixCSC{Events} structure where an element in the
matrix (i.e. an edge) is represented using the following structure
"""
immutable Events
    timestamps:: Vector{Tuple{Float64,Float64}}
end
Events() = Events(Vector{Tuple{Float64,Float64}}())
zero(::Type{Events}) = Events()
(==)(x::Events, y::Events) = x.timestamps == y.timestamps
mergesorted(x::Events, y::Events) =
    Events(mergesorted(x.timestamps,y.timestamps))

numevents(G::SparseMatrixCSC{Events}) =
    sum(map(x->length(x.timestamps), G.nzval))

flatten(G::SparseMatrixCSC{<:Number,Int}) = G
flatten(G::SparseMatrixCSC{Events,Int}) =
    SparseMatrixCSC(G.m,G.n,G.colptr,G.rowval,ones(Int,length(G.nzval)))

"""
    Total time during which events in event set are active
    i.e. active duration of an edge
"""
nodeactivity(evs::Events) = sum(x -> x[2]-x[1], evs.timestamps)

meannodeactivity(G::SparseMatrixCSC{Events}) = mean(map(nodeactivity, G.nzval))
"""
Input: adj. matrix of a dynamic network
Output: adj. matrix of static network containing corresponding
 duration of each edge in the dynamic network
"""
nodeactivity(G::SparseMatrixCSC{Events}) =
    SparseMatrixCSC(G.m,G.n,copy(G.colptr),copy(G.rowval),
                    map(nodeactivity, G.nzval))

"""
Sum of active duration over all events in G
"""
networkactivity(G::SparseMatrixCSC{Events}) = sum(nodeactivity, G.nzval)

function mintime(G::SparseMatrixCSC{Events})
    t_min = Inf
    for events in G.nzval, x in events
        t_min = min(t_min,x[1])
    end
    t_min
end
function maxtime(G::SparseMatrixCSC{Events})
    t_max = -Inf
    for events in G.nzval, x in events
        t_max = max(t_max,x[2])
    end
    t_max
end

"""
    If event intersects with [t_s,t_e] then add edges to snapshot
"""
function snapshot(G::SparseMatrixCSC{Events},t_s::Real,t_e::Real)
    m,n = size(G)
    I = Int[]
    J = Int[]
    for col = 1:n
        for j = nzrange(G,col)
            row = G.rowval[j]
            events = G.nzval[j]
            Tc,Tn = cet_ncet(events,[(t_s,t_e)])
            if Tc > 0
                push!(I,row)
                push!(J,col)
            end
        end
    end
    sparse(I,J,1,m,n)
end

""" Convert to snapshots with window size t_w
"""
function snapshots(G::SparseMatrixCSC{Events},t_w::Real,
                   t_min::Real=-Inf,t_max::Real=Inf)
    if isinf(t_min)
        t_min = mintime(G)
    end
    if isinf(t_max)
        t_max = maxtime(G)
    end
    map(t -> snapshot(G,t,t+t_w), t_min:t_w:t_max)
end

"""
Make each event in network wider by adding making each event
occur `pre` time ahead and `post` time later
"""
function widenevents!(G::SparseMatrixCSC{Events},pre=0.5,post=0.5)
    m,n = size(G)
    for i = 1:n
        for j = nzrange(G,i)
            row = G.rowval[j]
            events = G.nzval[j]
            for k = 1:length(events)
                x = events[k]
                events[k] = (x[1]-pre,x[2]+post)
            end
            fixevents!(G.nzval[j])
        end
    end
    G
end

"""
Given a vector of timestamps (sorted wrt the start time),
with possibly overlapping durations,
return a vector of timestamps with no overlapping durations
"""
fixevents(x::Events) = Events(fixevents(x.timestamps))
function fixevents(tis::Vector)
    n = length(tis)
    tout = similar(tis,0)
    a,b = tis[1]
    i = 2
    while i <= n
        c,d = tis[i]
        i += 1
        if b < c
            push!(tout,(a,b))
            a,b = c,d
        else
            b = d
        end
    end
    push!(tout,(a,b))
    tout
end
mergeevents(x) = fixevents(x)

"""
Given a vector of timestamps (sorted wrt the start time),
with possibly overlapping durations,
modify vector of timestamps s.t. it has no overlapping durations
"""
function fixevents!(x::Events)
    fixevents!(x.timestamps)
    x
end
function fixevents!(tis::Vector)
    n = length(tis)
    a,b = tis[1]
    i = 2
    j = 1
    while i <= n
        c,d = tis[i]
        i += 1
        if b < c
            tis[j] = (a,b)
            j += 1
            a,b = c,d
        else
            b = d
        end
    end
    tis[j] = (a,b)
    resize!(tis,j)
end
mergeevents!(x) = fixevents!(x)

"""
Given two vectors of timestamps (sorted wrt the start time),
with possibly overlapping durations,
return a vector of timestamps with no overlapping durations
"""
fixevents(x::Events, y::Events) =
    Events(fixevents(x.timestamps, y.timestamps))
fixevents(tis1::Vector, tis2::Vector) =
    fixevents!(mergesorted(tis1,tis2))
mergeevents(x,y) = fixevents(x,y)

"""
Calculates CET and NCET of the first `n` elements in a sorted array of events
As in, the events from two node pairs are combined together
and `tis` is the combined events
Sortedness is not checked
"""
cet_ncet(x::Events, n::Int=length(x.timestamps)) =
    cet_ncet(x.timestamps, n)
function cet_ncet(tis::Array,n::Int=length(tis))
    T = eltype(eltype(tis))
    if n > length(tis)
        error("n is larger than array size")
    end
    n==0 && return (zero(T),zero(T))
    (T_c,T_n) = (zero(T),zero(T))
    a,b = tis[1]
    i = 2
    while i <= n
        c,d = tis[i]
        i += 1
        if b <= c
            T_n += b-a
            a,b = c,d
        elseif b > d
            T_n += c-a
            T_c += d-c
            a,b = d,b
        else
            T_n += c-a
            T_c += b-c
            a,b = b,d
        end
    end
    T_n += b-a
    (T_c,T_n)
end

"""
Calculates CET and NCET of the events in tsl and tsr
Assumes the events in each array are non-overlapping
Sortedness is not checked
"""
cet_ncet(x::Events,y::Events) =
    cet_ncet(x.timestamps, y.timestamps)
function cet_ncet(tsl::Array,tsr::Array)
    T = eltype(eltype(tsl))
    nl = length(tsl)
    nr = length(tsr)
    (T_c,T_n) = (zero(T),zero(T))

    il = 1
    ir = 1

    a,b = (zero(T),zero(T))

    if il<nl && ir<nr
        if tsl[il][1] <= tsr[ir][1]
            a,b = tsl[il]
            il += 1
        else
            a,b = tsr[ir]
            ir += 1
        end
    end

    while il <= nl && ir <= nr
        if tsl[il][1] <= tsr[ir][1]
            c,d = tsl[il]
            il += 1
        else
            c,d = tsr[ir]
            ir += 1
        end

        if b <= c
            T_n += b-a
            a,b = c,d
        elseif b > d
            T_n += c-a
            T_c += d-c
            a,b = d,b
        else
            T_n += c-a
            T_c += b-c
            a,b = b,d
        end
    end

    while il <= nl
        c,d = tsl[il]
        il += 1

        if b <= c
            T_n += b-a
            a,b = c,d
        elseif b > d
            T_n += c-a
            T_c += d-c
            a,b = d,b
        else
            T_n += c-a
            T_c += b-c
            a,b = b,d
        end
    end

    while ir <= nr
        c,d = tsr[ir]
        ir += 1

        if b <= c
            T_n += b-a
            a,b = c,d
        elseif b > d
            T_n += c-a
            T_c += d-c
            a,b = d,b
        else
            T_n += c-a
            T_c += b-c
            a,b = b,d
        end
    end

    T_n += b-a
    (T_c,T_n)
end


text:
mutable struct DSA3217 <: AbstractScanivalve
    ipaddr::IPv4
    port::Int
    params::Dict{Symbol,Any}
    buffer::CircMatBuffer{UInt8}
    task::DAQTask
    chans::Vector{Int}
    channames::Vector{String}
    conf::DAQConfig
end

"Returns the IP address of the device"
ipaddr(dev::DSA3217) = dev.ipaddr
portnum(dev::DSA3217) = dev.port

AbstractDAQs.isreading(dev::DSA3217) = isreading(dev.task)
AbstractDAQs.samplesread(dev::DSA3217) = samplesread(dev.task)

numstring(x::Integer, n=2) = string(10^n+x)[2:end]

function Base.show(io::IO, dev::DSA3217)
    println(io, "Scanivalve DSA3217")
    println(io, "    Dev Name: $(daqdevname(dev))")
    println(io, "    IP: $(string(dev.ipaddr))")
end

function openscani(ipaddr::IPv4, port=23,  timeout=5)

    sock = TCPSocket()
    t = Timer(_ -> close(sock), timeout)
    try
        connect(sock, ipaddr, port)
    catch e
        error("Could not connect to $ipaddr ! Turn on the device or set the right IP address!")
    finally
        close(t)
    end

    return sock
end

openscani(dev::DSA3217,  timeout=5) = openscani(ipaddr(dev), portnum(dev), timeout)

function openscani(fun::Function, ipaddr::IPv4, port=23, timeout=5)

    io = openscani(ipaddr, port, timeout)
    try
        fun(io)
    finally
        close(io)
    end

end

openscani(fun::Function, dev::DSA3217, timeout=5) =
    openscani(fun, ipaddr(dev), portnum(dev), timeout)

function DSA3217(devname="Scanivalve", ipaddr="191.30.80.131";
                 timeout=5, buflen=300_000, tag="", sn="")

    ip = IPv4(ipaddr)
    port = 23

    openscani(ip, port, timeout) do s
        buffer = zeros(Int8, 0, 0)
        println(s, "SET EU 1")
        println(s, "SET AVG 16")
        println(s, "SET PERIOD 500")
        println(s, "SET FPS 1")
        println(s, "SET BIN 1")
        println(s, "SET XSCANTRIG 0")
        println(s, "SET UNITSCAN PA")
        println(s, "SET TIME 1")
    end

    params = Dict{Symbol,Any}(:FPS=>1, :AVG=>16, :PERIOD=>500, :TIME=>1,:XSCANTRIG=>0, :EU=>1, :UNITSCAN=>"PA")

    ipars = Dict{String,Int}("FPS"=>1, "AVG"=>16, "PERIOD"=>500, "TIME"=>1,
                             "EU"=>1, "XSCANTRIG"=>0)
    spars = Dict{String,String}("UNITSCAN"=>"PA")
    fpars = Dict{String,Float64}()

    conf = DAQConfig(ipars, fpars, spars,
                   devname=devname, ip=ipaddr, model="DSA3217", sn=sn,tag=tag)
    task = DAQTask()
    buf = CircMatBuffer{UInt8}(112, buflen)
    chn = "P" .* numstring.(1:16)
    return DSA3217(ip, port, params, buf, task, collect(1:16), chn, conf)

end

function daqpacketsize(::Type{DSA3217}, eu=1, time=1)
    if time==0
        return (eu==0) ? 72 : 104
    else
        return (eu==0) ? 80 : 112
    end
end

daqparam(dev, param) = dev.params[param]

framesize(dev::DSA3217) = daqpacketsize(DSA3217, daqparam(dev,:EU), daqparam(dev,:TIME))

deltat(dev) = 16*daqparam(dev,:PERIOD)*1e-6 * daqparam(dev,:AVG)

function stopscan(dev::DSA3217)
    tsk = dev.task

    if isreading(tsk)
        tsk.stop = true
    end
end

function scan!(dev::DSA3217)

    tsk = dev.task
    isreading(tsk) && error("DSA is already reading!")

    cleartask!(tsk)

    buf = dev.buffer
    fps = daqparam(dev, :FPS)
    if fps > capacity(buf)
        resize!(buf, fps)
    end
    empty!(buf)

    fsize = framesize(dev)
    nfrs = length(buf)

    if fps == 0 # Read continously
        fps1 = typemax(Int32)
    else
        fps1 = fps
    end

    δt = deltat(dev) # Time per frame
    stopped = false

    openscani(dev) do sock
        println(sock, "SCAN")
        tsk.isreading = true
        t0 = time_ns()
        b = nextbuffer(buf)
        readbytes!(sock, b, fsize)
        t1 = time_ns()
        settiming!(tsk, t0, t1, 1)
        tsk.nread += 1
        for i in 2:fps1
            if tsk.stop
                stopped = true
                println(sock, "STOP")
                sleep(0.5)
                break
            end

            b = nextbuffer(buf)
            readbytes!(sock, b, fsize)
            tn = time_ns()
            tsk.nread += 1
            settiming!(tsk, t1, tn, i-1)

        end

        tsk.isreading = false

        # If we stopped, try to read another frame
        if stopped
            delay = max(2δt, 1.0)
            ev = Base.Event()
            Timer(_ -> begin
                      timeout=true
                      notify(ev)
                  end, delay)
            @async begin
                read(sock, fsize)
                read!(sock, buffer(tsk, idx))
                notify(ev)
            end
            wait(ev)
        end
    end
    return
end

const validunits = ["ATM", "BAR", "CMHG", "CMH2O", "DECIBAR", "FTH2O", "GCM2",
                    "INHG", "INH2O", "KPA", "KGCM2", "KGM2", "KIPIN2", "KNM2",
                    "MH2O", "MMHG", "MPA", "NCM2", "MBAR", "OZFT2", "OZIN2",
                    "PA", "PSF", "NM2", "PSI", "TORR"]

function AbstractDAQs.daqconfig(dev::DSA3217; kw...)
    cmds = String[]

    if haskey(kw, :avg)
        avg = round(Int, kw[:avg])
        daqconfigdev(dev, AVG=avg)
    else
        avg = dev.params[:AVG]
    end

    if haskey(kw, :rate) && haskey(kw, :dt)
        error("Parameters `rate` and `dt` can not be specified simultaneously!")
    elseif haskey(kw, :rate) || haskey(kw, :dt)
        if haskey(kw, :rate)
            rate = kw[:rate]
            period = round(Int, 1.0 / (rate*16e-6*avg))
        else
            dt = kw[:dt]
            period = round(Int, dt / (16e-6*avg))
        end
        daqconfigdev(dev, PERIOD=period)
    else
        period = dev.params[:PERIOD]
    end

    if haskey(kw, :nsamples) && haskey(kw, :time)
        error("Parameters `nsamples` and `time` can not be specified simultaneously!")
    elseif haskey(kw, :nsamples) || haskey(kw, :time)
        if haskey(kw, :nsamples)
            nsamples = kw[:nsamples]
        else
            tt = kw[:time]
            dt = period * 16e-6 * avg
            nsamples = round(Int, tt / dt)
        end
        daqconfigdev(dev, FPS=nsamples)
    else
        nsamples = dev.params[:FPS]
    end

    if haskey(kw, :trigger)
        trigger = kw[:trigger]
        daqconfigdev(dev, XSCANTRIG=trigger)
    end

end

function AbstractDAQs.daqstart(dev::DSA3217, usethread=false)
    if isreading(dev)
        error("Scanivalve already reading!")
    end

    if usethread
        tsk = Threads.@spawn scan!(dev)
    else
        tsk = @async scan!(dev)
    end
    dev.task.task = tsk
    return tsk
end

function AbstractDAQs.daqaddinput(dev::DSA3217, chans=1:16; channames="P")

    cmin, cmax = extrema(chans)
    if cmin < 1 || cmax > 16
        throw(ArgumentError("Only channels 1-16 are available to DSA3217"))
    end

    if isa(channames, AbstractString) || isa(channames, Symbol)
        chn = string(channames) .* numstring.(chans)
    elseif length(channames) == length(chans)
        chn = string.(channames) .* string.(chans)
    else
        throw(ArgumentError("Argument `channames` should have length 1 or the length of `chans`"))
    end

    dev.chans = collect(chans)
    dev.channames = chn

    return
end

function AbstractDAQs.daqstop(dev::DSA3217)

    tsk = dev.task.task
    if !istaskdone(tsk) && istaskstarted(tsk)
        stopscan(dev)
        wait(tsk)
    end
    dev.task.stop = false
    dev.task.isreading = false

end

function readpressure(dev::DSA3217)
    isreading(dev) && error("Scanivalve still acquiring data!")

    tsk = dev.task
    buf = dev.buffer
    nsamples = length(buf)
    δt = getdaqtime(dev, nsamples)

    if daqparam(dev, :EU) > 0
        press = read_eu_press(buf, dev.chans)
    else
        error("Reading data without engineering units (EU=1) not implemented yet!")
    end
    return press, 1.0/δt
end

function AbstractDAQs.daqread(dev::DSA3217)

    # Check if the reading is continous
    if daqparam(dev, :FPS) == 0
        # Stop reading!
        daqstop(dev)
        sleep(0.1)
    end
    # Wait for task to end
    if !istaskdone(dev.task.task) && istaskstarted(dev.task.task)
        wait(dev.task.task)
        sleep(0.1)
    end

    # Get the data:
    return readpressure(dev)
end

function AbstractDAQs.daqacquire(dev::DSA3217)
    scan!(dev)
    return readpressure(dev)
end

function read_eu_press(buf, chans)

    nt = length(buf)
    nch = length(chans)
    P = zeros(Float32, nch, nt)

    for i in 1:nt
        p1 = reinterpret(Float32, view(buf[i], 9:72))
        for k in 1:nch
            P[k,i] = p1[chans[k]]
        end
    end
    return P
end

function getdaqtime(dev, nfr)
    buf = dev.buffer
    b = buf[1,1] # Identify the packet

    # Sampling time from scan configuration and fallback value
    δt₀ = daqparam(dev, :AVG) * daqparam(dev, :PERIOD)*1e-6 * 16
    btype = daqparam(dev, :TIME)
    if b == 0x04 || b == 0x05 || nfr==1 # No time
        # Calculate from scan configuration
        return δt₀
    else
        if b == 0x07  # EU with time
            bt1 = reinterpret(Int32, buf[105:108,1])[1]
            bt2 = reinterpret(Int32, buf[105:108,nfr])[1]
        elseif b == 0x06
            bt1 = reinterpret(Int32, buf[73:76,1])[1]
            bt2 = reinterpret(Int32, buf[73:76,nfr])[1]
        else # This shouldn't happen! Fallback
            return δt₀
        end
        ft = (btype==1) ? 1e6 : 1e3
        return  (bt2 - bt1) / (ft*(nfr-1))
    end

end

function scan2press(dev::DSA3217, buf, info)

    nframes = info.FPS
    press = zeros(Float32, nframes, 16)

    if info.EU > 0
        for i in 1:nframes
            press[i,:] = reinterpret(Float32, view(buf, 9:72,i))
        end
    else
        for i in 1:nframes
            press[i,:] = Float32.(reinterpret(Int16, view(buf, 9:40, i)))
        end
    end
    return press

end

import Base.close

function close(scani::DSA3217)
    openscani(dev) do s
        println(s, "SET EU 1")
        println(s, "SET AVG 100")
        println(s, "SET PERIOD 500")
        println(s, "SET FPS 1")
        println(s, "SET BIN 0")
        println(s, "SET XSCANTRIG 0")
        println(s, "SET UNITSCAN PA")
        println(s, "SET TIME 0")
    end

end

numchans(scani::DSA3217) = 16
AbstractDAQs.numchannels(scani::DSA3217) = length(scani.chans)
AbstractDAQs.daqchannels(scani::DSA3217) = scani.channames

#socket(scani) = scani.socket

function scanpacksize(dev::DSA3217, TIME=0, EU=1)
    if TIME==0
        if EU==1
            return 104
        else
            return 72
        end
    else
        if EU==1
            return 112
        else
            return 80
        end
    end

end

checkfps(dev::DSA3217, fps) = clamp(fps, 1, 1_800_000)
checkperiod(dev::DSA3217, period) = clamp(period, 125, 65535)
checkavg(dev::DSA3217, avg) = clamp(avg, 1, 240)
checkxscantrig(dev::DSA3217, xtrig) = (xtrig != 1) ? 0 : 1
checktime(dev::DSA3217, t) = clamp(t, 0, 2)
checkeu(dev::DSA3217, eu) = (eu != 0) ? 1 : 0

const validparameters = [:FPS, :PERIOD, :AVG, :TIME, :EU, :UNITSCAN, :XSCANTRIG]

function AbstractDAQs.daqconfigdev(dev::DSA3217; kw...)

    k = keys(kw)
    cmds = String[]
    if :FPS ∈ k
        fps = kw[:FPS]
        if 0 ≤ fps < 1_000_000
            push!(cmds, "SET FPS $fps")
            dev.params[:FPS] = fps
        else
            throw(DomainError(fps, "FPS outside range (0 - 1000000)"))
        end
    end

    if :PERIOD ∈ k
        period = kw[:PERIOD]
        if 125 ≤ period ≤ 65535
            push!(cmds, "SET PERIOD $period")
            dev.params[:PERIOD] = period
        else
            throw(DomainError(period, "PERIOD outside range (126 - 65535)"))
        end
    end

    if :AVG ∈ k
        avg = kw[:AVG]
        if 1 ≤ avg ≤ 240
            push!(cmds, "SET AVG $avg")
            dev.params[:AVG] = avg
        else
            throw(DomainError(avg, "AVG outside range (1 - 240)"))
        end
    end

    if :TIME ∈ k
        tt = kw[:TIME]
        if 0 ≤ tt ≤ 2
            push!(cmds, "SET TIME $tt")
            dev.params[:TIME] = tt
        else
            throw(DomainError(tt, "TIME outside range (0-2)"))
        end
    end

    if :EU ∈ k
        eu = kw[:EU]
        if 0 ≤ eu ≤ 1
            push!(cmds, "SET EU $eu")
            dev.params[:EU] = eu
        else
            throw(DomainError(eu, "EU outside range (0 or 1)"))
        end
    end

    if :UNITSCAN ∈ k # User should check manually if the correct unit was used
        unitscan = kw[:UNITSCAN]
        if unitscan ∈ validunits
            push!(cmds, "SET UNITSCAN $unitscan")
        else
            throw(DomainError(unitscan, "Invalid unit!"))
        end
    end

    if :XSCANTRIG ∈ k
        xscantrig = kw[:EU]
        if 0 ≤ eu ≤ 1
            push!(cmds, "SET SCANTRIG $xscantrig")
            dev.params[:XSCANTRIG] = xscantrig
        else
            throw(DomainError(xscantrig, "XSCANTRIG should be either 0 or 1"))
        end
    end

    # Send commands to Scanivalve
    openscani(dev) do sock
        for c in cmds
            println(sock, c)
        end

        # Update conf field:
        updateconf!(dev::DSA3217)
    end
end

function updateconf!(dev::DSA3217)
    p = dev.params
    ipars = dev.conf.ipars
    spars = dev.conf.spars

    ipars["FPS"] = p[:FPS]
    ipars["AVG"] = p[:AVG]
    ipars["PERIOD"] = p[:PERIOD]
    ipars["TIME"] = p[:TIME]
    ipars["XSCANTRIG"] = p[:XSCANTRIG]
    ipars["EU"] = p[:EU]

    spars["UNITSCAN"] = p[:UNITSCAN]

    return

end

function AbstractDAQs.daqzero(dev::DSA3217; time=15)
    openscani(dev) do io
        println(io, "CALZ")
        sleep(time)
    end
end

function readmanylines(dev, cmd, delay=0.5)

    lst = String[]
    openscani(dev) do socket
        println(socket, cmd)
        timeout = false
        while !timeout
            ev = Base.Event()
            Timer(_ -> begin
                      timeout=true
                      notify(ev)
                  end, delay)
            @async begin
                line = readline(socket)
                push!(lst, line)
                notify(ev)
            end
            wait(ev)
        end
    end

    return lst

end

function readstatus(dev::DSA3217)
    openscani(dev) do s
        println(s, "STATUS")
        string(Char.(read(s, 180)[81:100])...)
    end
end

function readnlines(s, n=1)
    lst = String[]
    for i in 1:n
        push!(lst, readline(s))
    end
    return lst
end

function parse_dsa_set(lst)
    p = Dict{String,String}()
    for l in lst
        s = split(l, " ")
        p[s[2]] = s[3]
    end
    return p
end

function listany(dev::DSA3217, cmd, nparams)
    openscani(dev) do sock
        println(sock, "LIST $cmd\r")
        readnlines(sock, nparams)
    end
end

listanydict(scani::AbstractScanivalve, cmd, nparams) = parse_dsa_set(listany(scani, cmd, nparams))

function listanyval(scani::AbstractScanivalve,
                    cmd, nparams, ::Type{T}=Int) where {T<:Real}
    lst = listany(scani, cmd, nparams)
    nlst = length(lst)
    vals = zeros(T, nlst)

    for i in 1:nlst
        s = split(lst[i], " ")
        vals[i] = parse(T, s[3])
    end

    return vals
end

listscan(scani::DSA3217) = listanydict(scani, "S", 14)
listident(scani::DSA3217) = listanydict(scani, "I", 4)
listzero(scani::DSA3217) = listanyval(scani, "Z", numchans(scani), Int)
listoffset(scani::DSA3217) = listanyval(scani, "O", numchans(scani), Float64)
listdelta(scani::DSA3217) = listanyval(scani, "O", numchans(scani), Float64)
listgain(scani::DSA3217) = listanyval(scani, "G", numchans(scani), Float64)

setparam(scani, param, val) =
    openscani(scani) do sock
        println(sock, "SET $param $val\r")
    end


text:
module Grid
using DataStructures

struct GridEnsemble
    width::Int
    height::Int
    size::Int
    n_districts::Int
    district_pop::Int
    weights::AbstractArray{Int}
    plans::AbstractArray{Int, 2}
    n_plans::Int
end

"""Load enumerated plans from a @zschutzman-style CSV."""
function GridEnsemble(plans_file::AbstractString, width::Int, height::Int, n_districts::Int)
    n_plans = countlines(plans_file)
    plans = zeros(Int, n_plans, width * height)
    weights = ones(Int, n_plans)

    open(plans_file) do file
        for (plan_idx, plan) in enumerate(eachline(file))
            assignments = [parse(Int, a) for a in split(plan, ",")]
            @assert length(assignments) == width * height
            @simd for i in 1:length(assignments)
                plans[plan_idx, i] = assignments[i]
            end
        end
    end

    district_pop = Int((width * height) / n_districts)
    plan_size = width * height
    return GridEnsemble(width, height, plan_size, n_districts, district_pop,
                        weights, plans, size(plans, 1))
end

function unique_districts(ensemble::GridEnsemble)
    district_counts = DefaultDict{Tuple, Int}(0)
    for plan_idx in 1:ensemble.n_plans
        plan = ensemble.plans[plan_idx, :]
        districts = [[] for _ in 1:ensemble.n_districts]
        for (idx, node) in enumerate(plan)
            push!(districts[node], idx)
        end
        for dist in districts
            district_counts[tuple(sort(dist))] += 1
        end
    end
    return district_counts
end

const ensemble = GridEnsemble("enum.csv", 6, 6, 6)
const districts = unique_districts(ensemble)
const tie = Int(ensemble.district_pop / 2)
const win = tie + 1

function expected_seat_share(voters::Array{Int})
    binary_voters = zeros(Int, ensemble.size)
    for v in voters
        binary_voters[v] = 1
    end

    total_wins = 0
    for (dist, freq) in districts
        voters_in_dist = sum(binary_voters[dist[1]])
        if voters_in_dist == tie
            total_wins += 0.5 * freq
        elseif voters_in_dist >= win
            total_wins += freq
        end
    end
    total_wins / sum(values(districts))
end

function seat_hist(ensemble::GridEnsemble, voters::Array{Int})::AbstractDict
    binary_voters = zeros(Int, ensemble.size)
    for v in voters
        binary_voters[v] = 1
    end
    for r in eachrow(reshape(binary_voters, ensemble.width, ensemble.height)')
        println(r)
    end
    hist = DefaultDict{Float64, Int}(0)
    for i in 1:ensemble.n_plans
        binary_plan = zeros(Int, ensemble.n_districts, ensemble.size)
        for (j, assignment) in enumerate(ensemble.plans[i, :])
            binary_plan[assignment, j] = 1
        end
        counts = binary_plan * binary_voters
        wins = 0
        split = (ensemble.width * ensemble.height) / (ensemble.n_districts * 2)
        for c in counts
            if c == split
                wins += 0.5
            elseif c > split
                wins += 1
            end
        end
        if wins == 2
            println(ensemble.plans[i, :])
        end
        hist[wins] += 1
    end
    hist
end

export expected_seat_share, ensemble, districts, seat_hist

end


text:
"""
`LengthfulIter` is an iterator wrapper that allows length info known programmatically to be used.
It's used because iterator transformations, like filters, don't keep the length info,
which might however be known to the programmer.
"""
struct LengthfulIter{T,IterT}
  iter::IterT
  n::Int
end
LengthfulIter{T}(iter::IterT,n) where {T,IterT}= LengthfulIter{T,IterT}(iter,n)
@inline Base.length(li::LengthfulIter)= li.n
@inline Base.iterate(li::LengthfulIter)=    Base.iterate(li.iter)
@inline Base.iterate(li::LengthfulIter, s)= Base.iterate(li.iter,s)
@inline Base.eltype(::Type{LengthfulIter{T}}) where T= T
@inline Base.collect(li::LengthfulIter{T}) where T= _collect(li.iter,li.n,T)
# Stolen from @array.jl#600
function _collect(itr::Base.Generator,sz::Int, elT)
  _array_for(::Type{T}, itr) where {T} = Vector{T}(undef, sz)
  y= iterate(itr)
  y === nothing && return _array_for(elT, itr.iter)
  v1, st = y
  Base.collect_to_with_first!(_array_for(typeof(v1), itr.iter), v1, itr, st)
end
_collect(itr,n,elT)= Base.collect(itr)

"""
`TruesOf` is a non-allocating version of `findall`:
iterates over the trues of a `BitArray`

# Examples
```jldoctest; setup= :(import Random; Random.seed!(5))
julia> b= [true,true,true,false,false]|> BitVector
5-element BitVector:
 1
 1
 1
 0
 0

julia> foreach(i-> print(string(i)*" "), HTM.Truesof(b))
1 2 3
```
"""
struct Truesof
  b::BitArray
end
@inline Base.length(B::Truesof)= count(B.b)
@inline Base.eltype(::Type{Truesof})= Int
Base.iterate(B::Truesof, i::Int=1)= begin
  i= findnext(B.b, i)
  i === nothing ? nothing : (i, i+1)
end
Base.collect(B::Truesof)= collect(B.b)

"""
`sparse_foreach(f, s::SparseMatrixCSC,columnIdx)` iterates over the columns of the sparse matrix `s` specified by
`columnIdx` and applies `f(columnElt,columnRowIdx)` to each.
Column elements are a @view into the nonzero elements at the given column of `s`.

See also: [`sparse_map`](@ref)
"""
sparse_foreach(f, s::SparseMatrixCSC,columnIdx)=
  foreach(Truesof(columnIdx)) do c
    ci= nzrange(s,c)
    @views f(nonzeros(s)[ci], rowvals(s)[ci])
  end

"""
`sparse_map(f, s::SparseMatrixCSC,columnIdx)` maps each column of the sparse matrix `s` specified by
`columnIdx` to `f(columnElt,columnRowIdx)`.
Column elements are a @view into the nonzero elements at the given column of `s`.

See also: [`sparse_foreach`](@ref)
"""
sparse_map(f, s::SparseMatrixCSC,columnIdx)=
  map(Truesof(columnIdx)) do c
    ci= nzrange(s,c)
    f((@view nonzeros(s)[ci]),rowvals(s)[ci])
  end

"""
`@percolumn(f,a,b,k)`

Macro to fold a large vector `a` per columns `k` and apply binary operator `f` elementwise
- `a`: vector of size `Nc*k`, column-major (important because it will fold every `k` elements)
- `b`: vector of size `Nc`
"""
macro percolumn(f,a,b,k)
  esc(:( $f.(reshape($a,$k,:), $b') ))
end
"""
`@percolumn(reduce,a,k)`

Macro to `reduce(a)` per column `k`.
"""
macro percolumn(reduce,a,k)
  esc(:( $reduce(reshape($a,$k,:),dims=1)|> vec ))
end

"""
`bitarray(idx, dims)`

Create a bitarray with `true` only at `idx`.
"""
function bitarray(idx,dims)
  r= falses(dims)
  r[idx].= true
  return r
end
function bitarray(idx::Int,dims)
  r= falses(dims)
  r[idx]= true
  return r
end
padfalse(b::BitArray,dim)= [b;falses(dim-length(b))]
padfalse(b::Vector{T},dim) where T= [b;zeros(T,dim-length(b))]

"""
`Sequence(;δ,init)` is an easy way to define sequences with transition function `δ` and starting condition `init`
as an iterator.

# Examples

Fibonacci sequence:
```jldoctest exp; setup= :(import IterTools, Lazy)
fibonacci_δ(a,b)= b, a+b;
fibonacci_init= 0,1;
Lazy.@as _1 HTM.Sequence(fibonacci_δ, fibonacci_init) IterTools.take(_1, 13) foreach(println, _1)

# output
0
1
1
2
3
5
8
13
21
34
55
89
144
```

Factorial sequence:
```jldoctest exp
factorial_δ(a,b)= a*b, b+1;
factorial_init= 1,1;
Lazy.@as _1 HTM.Sequence(factorial_δ, factorial_init) IterTools.take(_1, 13) foreach(println, _1)

# output
1
1
2
6
24
120
720
5040
40320
362880
3628800
39916800
479001600
```
"""
struct Sequence{F <: Function, I}
  δ::F
  init::I
  Sequence(δ::F,init::I) where {F,I}= new{F,I}(δ,init)
  Sequence(;δ::F,init::I) where {F,I}= new{F,I}(δ,init)
end

Base.iterate(s::Sequence) = iterate(s, s.init)
Base.iterate(s::Sequence, state) = state[1], s.δ(state...)
Base.IteratorSize(::Sequence) = Base.SizeUnknown()
Base.IteratorEltype(::Sequence) = Base.EltypeUnknown()


text:
module UMAP

using Arpack
using Distances
using LinearAlgebra
using LsqFit: curve_fit
using NearestNeighborDescent
using SparseArrays

include("utils.jl")
include("embeddings.jl")
include("umap_.jl")

export umap, UMAP_, transform

end # module


text:
using DashBootstrapComponents, DashHtmlComponents

modal = html_div([
    dbc_button("Small modal", id = "open-sm", className = "me-1", n_clicks = 0),
    dbc_button("Large modal", id = "open-lg", className = "me-1", n_clicks = 0),
    dbc_button("Extra large modal", id = "open-xl", n_clicks = 0),
    dbc_modal(
        [dbc_modalheader(dbc_modaltitle("Header")), dbc_modalbody("A small modal.")],
        id = "modal-sm",
        size = "sm",
        is_open = false,
    ),
    dbc_modal(
        [dbc_modalheader(dbc_modaltitle("Header")), dbc_modalbody("A large modal.")],
        id = "modal-lg",
        size = "lg",
        is_open = false,
    ),
    dbc_modal(
        [dbc_modalheader(dbc_modaltitle("Header")), dbc_modalbody("An extra large modal.")],
        id = "modal-xl",
        size = "xl",
        is_open = false,
    ),
]);

function toggle_modal(n1, is_open)
    return n1 > 0 ? is_open == 0 : is_open
end;

callback!(
    app,
    Output("modal-sm", "is_open"),
    Input("open-sm", "n_clicks"),
    State("modal-sm", "is_open"),
) do n1, is_open
    return toggle_modal(n1, is_open)
end;

callback!(
    app,
    Output("modal-lg", "is_open"),
    Input("open-lg", "n_clicks"),
    State("modal-lg", "is_open"),
) do n1, is_open
    return toggle_modal(n1, is_open)
end;

callback!(
    app,
    Output("modal-xl", "is_open"),
    Input("open-xl", "n_clicks"),
    State("modal-xl", "is_open"),
) do n1, is_open
    return toggle_modal(n1, is_open)
end;

max_stars_repo_name: glsdown/dash-bootstrap-components

max_issues_repo_name: glsdown/dash-bootstrap-components

max_forks_repo_name: glsdown/dash-bootstrap-components


text:
export load_narrowpeak, filterPreds, combinePeaks, binarizePeaks, window_bed_file

function load_narrowpeak(stream, contigs, index; binSize=1000, loadp=true, mlogt=false, verbose=0)
    # Get number of bins
    numBins = ceil(Int64, sum(contigs.sizes) / binSize)
    # Create offsets difctionary
    chrOffsets = Dict{String,Int64}()
    for i in 1:contigs.count
        chrOffsets[contigs.names[i]] = contigs.offsets[i]
    end

    # mark all bins that are touched with 1
    binValues = falses(numBins)
    # also record the highest p-value
    pValues = zeros(numBins)
    # also keep track how many peaks were evaluated
    count = 0
    for line in eachline(stream)
        parts = split(rstrip(line), '\t')
        if haskey(chrOffsets, parts[1])
            count += 1
            if count < verbose
                println(parts)
            end
            startPos = ceil(Int64, (chrOffsets[parts[1]]+Int(parse(Float64, parts[2])))/binSize)
            endPos = ceil(Int64, (chrOffsets[parts[1]]+Int(parse(Float64, parts[3])))/binSize)
            for i in startPos:endPos
                assert(i!=0)
                # record bin that it was touched
                binValues[i] = true
                # loadp
                if loadp
                    # minus log 10 p-value transformation if necessary.
                    if mlogt
                        pval = -1*log(10, parse(Float64, parts[index]))
                    else
                        pval = parse(Float64, parts[index])
                    end
                    # record p-value if more significant.
                    if pValues[i] < pval
                        pValues[i] = pval
                    end
                end
            end
        end
    end
    close(stream)

    #assert(sum(binValues) == sum([i>0 for i in pValues]))
    binValues, pValues
end

function window_bed_file(stream, contigs; binSize=1000)
    numBins = ceil(Int64, sum(contigs.sizes) / binSize)
    chrOffsets = Dict{String,Int64}()
    for i in 1:contigs.count
        chrOffsets[contigs.names[i]] = contigs.offsets[i]
    end

    # mark all bins that are touched with 1
    binValues = falses(numBins)
    for line in eachline(stream)
        parts = split(line, '\t')
        if haskey(chrOffsets, parts[1])
            startPos = ceil(Int64, (chrOffsets[parts[1]]+parse(Int64, parts[2]))/binSize)
            endPos = ceil(Int64, (chrOffsets[parts[1]]+parse(Int64, parts[3]))/binSize)
            for i in startPos:endPos
                binValues[i] = true
            end
        end
    end

    binValues
end

# Get top x predictions with its truth values.
# Returns truth and preds.
function filterPreds(top, truth, pred)
    temp = Any[]
    for i in 1:length(truth)
        push!(temp, (truth[i], pred[i]))
    end
    temp = sort(temp, by= x -> x[2], rev=true)[1:top]
    [i[1] for i in temp], [i[2] for i in temp]
end

# Combine signal vector a and b by taking minimum element wise.
# Returns a vector.
function combinePeaks(a, b)
    ret = zeros(length(a))
    for i in 1:length(a)
        ret[i] = minimum((a[i], b[i]))
    end
    ret
end

# Extract top X peaks and output binarized vector of same size.
# Returns a vector.
function binarizePeaks(top, pred)
    temp = Any[]
    for i in 1:length(pred)
        push!(temp, (i, pred[i]))
    end
    sort!(temp, by=x->x[2], rev=true)

    ret = falses(length(pred))
    for i in 1:top
        ret[temp[i][1]] = true
    end
    ret
end


text:
@testset "RHS parse" begin
    ex = quote
        (a[1:end, :][1, 2][a, b'] + sin(cos(B))[b', (a')']) * (C * D + E)[(a')', a] * 3 * 2 + α * a[a, a] * K[1, 2, 3][a, a] -
        (((L[a, b] * P.P[b, c]) * (M[c, d] * N[d, e]) * Z[e, f]) * D[f, a]) * π
    end
    @test rhs_to_args(ex)[2] == [
        :(a[1:end, :][1, 2]),
        :(sin(cos(B))),
        :(C * D + E),
        :α,
        :a,
        :(K[1, 2, 3]),
        :L,
        :(P.P),
        :M,
        :N,
        :Z,
        :D,
        :π,
    ]

    ex = :(-a * (-b + (c + d) + (-e) - f - (g * h)) * i + j)
    @test make_only_product(ex, :a) == :(-a * (-b + (c + d) + (-e) - f - (g * h)) * i)
    @test make_only_product(ex, :b) == :(-a * -b * i)
    @test make_only_product(ex, :c) == :(-a * c * i)
    @test make_only_product(ex, :d) == :(-a * d * i)
    @test make_only_product(ex, :e) == :(-a * -e * i)
    @test make_only_product(ex, :f) == :(-a * -f * i)
    @test make_only_product(ex, :g) == :(-a * -(g * h) * i)
    @test make_only_product(ex, :h) == :(-a * -(g * h) * i)
    @test make_only_product(ex, :i) == :(-a * (-b + (c + d) + (-e) - f - (g * h)) * i)
    @test make_only_product(ex, :j) == :j

    ex = :(-a * (-conj(conj(b)) + conj(conj(c + d * e) + f) + -g))
    @test make_only_product(ex, :a) ==
          :(-a * (-conj(conj(b)) + conj(conj(c + d * e) + f) + -g))
    @test make_only_product(ex, :b) == :(-a * -conj(conj(b)))
    @test make_only_product(ex, :c) == :(-a * conj(conj(c)))
    @test make_only_product(ex, :d) == :(-a * conj(conj(d * e)))
    @test make_only_product(ex, :e) == :(-a * conj(conj(d * e)))
    @test make_only_product(ex, :f) == :(-a * conj(f))
    @test make_only_product(ex, :g) == :(-a * -(g))
end


text:
function fdc(order::Int64)

    # This function calculate FD coefficient for spatial finite difference
    # Currently we define the coefficient as the one for the 4th order
    # search "Finite difference coefficient" online for different order

    if order == 4
        FDC = [1/24, -9/8, 9/8, -1/24] # Finite-Difference Coefficient
        else error("please choose choose 4")
    end
    return FDC
end


text:
struct ThreadTask
    p::Ptr{UInt}
end
Base.pointer(tt::ThreadTask) = tt.p

@inline taskpointer(tid) = THREADPOOLPTR[] + tid*THREADBUFFERSIZE

function _call(p::Ptr{UInt})
    fptr = load(p + sizeof(UInt), Ptr{Cvoid})
    assume(fptr ≠ C_NULL)
    ccall(fptr, Cvoid, (Ptr{UInt},), p)
end
@inline function launch(f::F, tid::Integer, args::Vararg{Any,K}) where {F,K}
    p = taskpointer(tid)
    f(p, args...)
    state = _atomic_xchg!(p, TASK)
    state == WAIT && wake_thread!(tid % Int)
    return nothing
end

function (tt::ThreadTask)()
    p = pointer(tt)
    max_wait = one(UInt32) << 20
    wait_counter = max_wait
    GC.@preserve THREADPOOL begin
        while true
            if _atomic_state(p) == TASK
                _call(p)
                _atomic_store!(p, SPIN)
                wait_counter = zero(UInt32)
                continue
            end
            pause()
            if (wait_counter += one(UInt32)) > max_wait
                wait_counter = zero(UInt32)
                if _atomic_cas_cmp!(p, SPIN, WAIT)
                    Base.wait()
                    _call(p)
                    _atomic_cas_cmp!(p, TASK, SPIN)
                end
            end
        end
    end
end

# 1-based tid, pushes into task 2-nthreads()
# function wake_thread!(tid::T) where {T <: Unsigned}
function wake_thread!(tid::T) where {T <: Integer}
    tidp1 = tid + one(tid)
    assume(unsigned(length(Base.Workqueues)) > unsigned(tid))
    assume(unsigned(length(TASKS)) > unsigned(tidp1))
    @inbounds push!(Base.Workqueues[tidp1], TASKS[tid])
    ccall(:jl_wakeup_thread, Cvoid, (Int16,), tid % Int16)
end

# 1-based tid
@inline wait(tid::Integer) = wait(taskpointer(tid))
@inline function wait(p::Ptr{UInt})
    # note: based on relative values (SPIN = 0, WAIT = 1)
    # thus it should spin for as long as the task is doing anything else
    # while @show(stacktrace(), reinterpret(UInt, _atomic_load(p))) > reinterpret(UInt, WAIT)
    while _atomic_load(p) > reinterpret(UInt, WAIT)
    # while reinterpret(UInt, @show(_atomic_load(p))) > reinterpret(UInt, WAIT)
        pause()
    end
end

# function launch_thread(f::F, tid) where {F}
#     cfunc = @cfunction($mapper, Cvoid, ());

#     fptr = Base.unsafe_convert(Ptr{Cvoid}, cfunc)

#     ccall(fptr, Cvoid, ())

# end


text:
#!/usr/local/bin/julia

using ArgParse

function countNucleotides(DNAstring::String)
	A = count(c->c=='A', DNAstring)
	C = count(c->c=='C', DNAstring)
	G = count(c->c=='G', DNAstring)
	T = count(c->c=='T', DNAstring)

	return A, C, G, T
end

function main(args)
	# Argument parsing
	configs = ArgParseSettings()

	add_arg_group!(configs, "Options:")
	@add_arg_table! configs begin
		"--input", "-i"
			help	 = "Input filename"
			required = true
			action	 = :store_arg
			arg_type = String
		"--output", "-o"
			help	 = "Output filename"
			required = true
			action	 = :store_arg
			arg_type = String
		end
	# Store arguments to variables
	parsed_args = parse_args(args, configs)
	infile      = parsed_args["input"]
	outfile     = parsed_args["output"]

	# Read file and count nucleotides
	DNAstring  = open(f->read(f, String), infile)
	@show A, C, G, T = countNucleotides(DNAstring)
	open(f->write(f, "$A $C $G $T"), outfile, "w+")

end

main(ARGS)


text:
function find_mode(C::Chains;n=1000,k=10)
    Θ = sample(C,n)
    dk = [sort([norm(Θ[j] - Θ[i]) for j = 1:n])[k+1] for i = 1:n]
    return Θ[findmin(dk)[2]]
end

max_stars_repo_name: ap-browning/internalisation

max_issues_repo_name: ap-browning/internalisation

max_forks_repo_name: ap-browning/internalisation


text:
for filename in readdir("/foo/bar")
    if endswith(filename, ".mp3")
        print(filename)
    end
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# http://web.mit.edu/puzzle/www/2012/puzzles/okla_holmes_a/solution/

@testset "okla holmes-a meta" begin
    puzzle = wordlist(IOBuffer("""
        CARPAL
        THE SOUTH
        STERNO
        BYLINE
        SO CLOSE
        BUFFOON
        VESTIGE
        """))
    f = best_feature(puzzle)
    @test f.description == "can be completely broken down into chemical element symbols"
end


text:
include("modulos/bondad.jl")
include("modulos/graphics.jl")
include("modulos/datasets.jl")
include("modulos/rna.jl")
include("modulos/attributes_from_dataset.jl")

using Flux;

# Parametros principales de la RNA y del proceso de entrenamiento
topology = [4]; # Dos capas ocultas con 4 neuronas la primera y 3 la segunda
learningRate = 0.01; # Tasa de aprendizaje
numMaxEpochs = 1000; # Numero maximo de ciclos de entrenamiento
validationRatio = 0.2; # Porcentaje de patrones que se usaran para validacion
testRatio = 0.2; # Porcentaje de patrones que se usaran para test
maxEpochsVal = 15; # Numero de ciclos en los que si no se mejora el loss en el conjunto de validacion, se para el entrenamiento

#Si no está generado el dataset pues lo creamos
dataset_name="datasets/faces.data"
if (!isfile(dataset_name))
    (inputs, targets) = getInputs("datasets");
    println("Tamaños en la generación:")
    println(size(inputs))
    println(size(targets))
    write_dataset(dataset_name,inputs,targets)
end

# Cargamos el dataset
dataset = readdlm(dataset_name,',');

# Preparamos las entradas y las salidas deseadas
inputs = convert(Array{Float64,2}, dataset[:,1:42]);             #Array{Float64,2}
targets = oneHotEncoding(convert(Array{Any,1},dataset[:,43]));   #Array{Bool,2}

# Normalizamos las entradas, a pesar de que algunas se vayan a utilizar para test
#normalizeMinMax!(inputs);

# Creamos los indices de entrenamiento, validacion y test
(trainingIndices, validationIndices, testIndices) = holdOut(size(inputs,1), validationRatio, testRatio);

# Dividimos los datos
trainingInputs    = inputs[trainingIndices,:];
validationInputs  = inputs[validationIndices,:];
testInputs        = inputs[testIndices,:];
trainingTargets   = targets[trainingIndices,:];
validationTargets = targets[validationIndices,:];
testTargets       = targets[testIndices,:];

println("Empezando el entrenamiento.")
println()
# Y creamos y entrenamos la RNA con los parametros dados
(ann, trainingLosses, validationLosses, testLosses, trainingAccuracies,
validationAccuracies, testAccuracies) = trainClassANN(topology, trainingInputs, trainingTargets, validationInputs,
                                            validationTargets, testInputs, testTargets; maxEpochs=numMaxEpochs,
                                            learningRate=learningRate, maxEpochsVal=maxEpochsVal, showText=false);

print_train_results(trainingLosses, validationLosses, testLosses, trainingAccuracies, validationAccuracies, testAccuracies);

# esto podríamos meterlo dentro de la matriz de confusión
outputs=convert(Array{Float64,2},ann(inputs')');
outputs=classifyOutputs(outputs); #Array{Bool,2}
printConfusionMatrix(outputs, targets);


text:
# BBI
# ===
#
# The BigBed format is documented in
#   Kent, W. James, et al. "BigWig and BigBed: # enabling browsing of large
#   distributed datasets." Bioinformatics 26.17 (2010): 2204-2207.
#
# The low level details are documented in a series of tables in the supplement
# of that paper. The immutable types defined below are exactly the data layout
# described in those tables. They are labeled with the table they correspand to.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

include("types.jl")
include("reader.jl")
include("writer.jl")
include("intersection.jl")


text:
# Generate 3D data for the overthrust model
# Author: pwitte.slim@gmail.com
# Date: January 2017
#

using JUDI, HDF5

# Load overthrust model
if ~isfile("$(JUDI.JUDI_DATA)/overthrust_3D_true_model.h5")
    ftp_data("ftp://slim.gatech.edu/data/SoftwareRelease/WaveformInversion.jl/3DFWI/overthrust_3D_true_model.h5")
end
n, d, o, m = read(h5open("$(JUDI.JUDI_DATA)/overthrust_3D_true_model.h5","r"), "n", "d", "o", "m")

# Set up model structure
model = Model((n[1], n[2], n[3]), (d[1], d[2], d[3]), (o[1], o[2], o[3]), m)

# Set up source grid
nsrc = 97^2
xsrc = 400f0: 200f0: 19600f0	# 400 m space from boundary
ysrc = 400f0: 200f0: 19600f0
zsrc = 50f0
(xsrc, ysrc, zsrc) = setup_3D_grid(xsrc, ysrc, zsrc)
xsrc = convertToCell(xsrc)
ysrc = convertToCell(ysrc)
zsrc = convertToCell(zsrc)
timeS = 3000f0
dtS = 4f0
src_geometry = Geometry(xsrc, ysrc, zsrc; dt = dtS, t = timeS)

## Set up receivers
xmin = 100f0; xmax = 19900f0
ymin = 100f0; ymax = 19900f0
offset = 6000f0
nxrec = 240
nyrec = 240
xrec = Array{Any,1}(undef, nsrc)
yrec = Array{Any,1}(undef, nsrc)
zrec = Array{Any}(undef, nsrc)
for j = 1: nsrc
	# 6 km max. offset around current source
    xsrc[j] - offset < xmin ? xlocal1 = xmin : xlocal1 = xsrc[j] - offset
    xsrc[j] + offset > xmax ? xlocal2 = xmax : xlocal2 = xsrc[j] + offset
	ysrc[j] - offset < ymin ? ylocal1 = ymin : ylocal1 = ysrc[j] - offset
	ysrc[j] + offset > ymax ? ylocal2 = ymax : ylocal2 = ysrc[j] + offset
	xlocal = xlocal1: 50f0: xlocal2
	ylocal = ylocal1: 50f0: ylocal2
	zlocal = 500f0
	(xlocal,ylocal,zlocal) = setup_3D_grid(xlocal, ylocal, zlocal)
	xrec[j] = xlocal
	yrec[j] = ylocal
	zrec[j] = zlocal
end
timeR = 3000f0	# receiver recording time [ms]
dtR = 4f0	# receiver sampling interval
rec_geometry = Geometry(xrec, yrec, zrec; dt = dtR, t = timeR)

# Set up source
f0 = 0.008
wavelet = ricker_wavelet(timeS, dtS, f0)

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

# Set up modeling options
opt = Options(limit_m = true,
			  save_data_to_disk = true,
			  file_path = "path/to/shot/records",   # replace w/ path to directory in which data will be stored (~total of 1.7 TB)
 			  file_name = "overthrust_3D_shot_"   # adds x-y coordinates to name
              )

# Set up operators
Pr = judiProjection(rec_geometry)
Ps = judiProjection(src_geometry)
A_inv = judiModeling(model; options=opt)
q = judiVector(src_geometry, wavelet)

# Generate data and save as individual SEG-Y files to disk
d_obs = Pr*A_inv*adjoint(Ps)*q


text:
"""
    Res0dWithBN <: Block
Res0dWithBN is a residual whose neuron processing 0-d data. which has 2 BatchNorm0d blocks

# Structure
       Rⁱ → Rᵐ                  Rᵐ → Rⁱ
        ┌───┐   ┌────┐   ┌───┐   ┌───┐   ┌────┐   ┌───┐   ┌───┐
    ─┬─►│ W ├──►│ BN ├──►│f()├──►│ W ├──►│ BN ├──►│ + ├──►│f()├──►
     │  └───┘   └────┘   └───┘   └───┘   └────┘   └─┬─┘   └───┘
     └─────────────────────────►────────────────────┘
"""
mutable struct Res0dWithBN <: Block
    blocks::Vector
    function Res0dWithBN(i::Int, m::Int; type::Type=Array{Float32})
        a1 = Affine(i, m, type=type)
        b1 = BatchNorm0d(m)
        a2 = Affine(m, i, type=type)
        b2 = BatchNorm0d(i)
        new([a1, b1, a2, b2])
    end
    function Res0dWithBN()
        new(Vector(undef,4))
    end
end

@extend(Res0dWithBN, blocks)

function clone(this::Res0dWithBN; type::Type=Array{Float32})
    cloned = Res0dWithBN()
    cloned[1] = clone(this[1], type=type)
    cloned[2] = clone(this[2], type=type)
    cloned[3] = clone(this[3], type=type)
    cloned[4] = clone(this[4], type=type)
    return cloned
end

function forward(m::Res0dWithBN, x0)
    x1 = forward(m[1], x0) # Affine
    x2 = forward(m[2], x1) # batchnorm0d
    x3 =          relu(x2) # activation
    x4 = forward(m[3], x3) # Affine
    x5 = forward(m[4], x4) # batchnorm0d
    return relu(x5 + x0)   # short connectin + activation
end

function predict(m::Res0dWithBN, x0)
    x1 = predict(m[1], x0) # Affine
    x2 = predict(m[2], x1) # batchnorm0d
    x3 =          relu(x2) # activation
    x4 = predict(m[3], x3) # Affine
    x5 = predict(m[4], x4) # batchnorm0d
    return relu(x5 + x0)   # short connectin + activation
end

max_stars_repo_name: sonosole/Delta

max_issues_repo_name: sonosole/Delta

max_forks_repo_name: sonosole/Delta