text:
using TunedSmoother,WaterLily,Plots
using GeometricMultigrid: Vcycle!

# shark parameters from test_waterlily
p = [0.0673015, -0.104415, -0.024194, -0.182964, -0.298915]

A,St,remeasure = 0.1,0.3,true
cycles,period,steps = 2,2A/St,120
t₀,Δt = 3,cycles*period/steps

# Set up and run simulation
sim = shark(6;A,St)
sim_step!(sim,t₀;remeasure)

# Get Poisson data over next period
data = [(TunedSmoother.sim_collect_poisson!(sim;Δt,remeasure),sim_time(sim)) for i ∈ 1:steps]

# Integrate pressure data at time t
get_forces(data;kw...) = map(data) do ((A,p,b),t)
	st = state(A,b;kw...)
	mg!(st;kw...)
	WaterLily.∮nds(st.x.data,sim.body,t*sim.L/sim.U)./(0.5*sim.L*sim.U^2)
end

begin
	reltol=1e-3
	thrust = plot(xlabel="scaled time", ylabel="scaled thrust force")
	side = plot(xlabel="scaled time", ylabel="scaled side force")
	forces = get_forces(data;reltol)
	plot!(thrust,last.(data),first.(forces),label="GS",c=palette(:default)[2])
	plot!(side,  last.(data), last.(forces),label="GS",c=palette(:default)[2])
	forces = get_forces(data;reltol,smooth! = pseudo!, p)
	plot!(thrust,last.(data),first.(forces),label="Ã⁻¹",line=:dash,c=palette(:default)[3])
	plot!(side,  last.(data), last.(forces),label="Ã⁻¹",line=:dash,c=palette(:default)[3])
end
plot(thrust,side,layout=(2,1))
savefig("forces.pdf")


text:
bsplines_exact(funs, x::Irrational, indices) = bsplines_exact(funs, big(x), indices)
bsplines_exact(funs, x::Irrational, indices, drv::AllDerivatives{N}) where N =
    bsplines_exact(funs, big(x), indices, drv)

function bsplines_exact(funs, x, indices)
    bigx = maybebig(x)
    OffsetArray([funs[i,1](bigx) for i=indices], indices)
end
function bsplines_exact(funs, x, indices, ::AllDerivatives{N}) where N
    bigx = maybebig(x)
    k = length(indices)
    k == size(funs, 2) || error("got $k indices, expected $(size(funs,2))")
    arr = [funs[i,d](bigx) for i=indices, d=1:min(k,N)]
    arr = hcat(arr, zeros(eltype(arr), k, max(0,N-k)))
    OffsetArray(arr, indices, 0:N-1)
end

# Enlarge x for calculating exact B-splines:
# * for Integer/Rational: use 64-bit to avoid overflow in rational arithmetic
# * for floating-point: use BigFloat to minimize error
maybebig(x::Integer) = Int64(x)
maybebig(x::Rational) = Rational{Int64}(x)
maybebig(x::AbstractFloat) = big(x)

@time @testset "bsplines" begin
    @testset "Return types" begin
        # b1 = BSplineBasis(1, 0:5)
        b1_types = [(Int,Float64), (Float64,Float64), (Float32,Float32),
                    (Rational{Int},Rational{Int}), (BigInt,BigFloat),
                    (BigFloat,BigFloat), (Rational{BigInt},Rational{BigInt}),
                    (Rational{Int16},Rational{Int})]
        # b2 = BSplineBasis(2, big(0):1:big(5))
        b2_types = [(Int,BigFloat), (Float32,BigFloat), (Rational{Int},Rational{BigInt})]
        # b3 = BSplineBasis(3, 0//1:5//1)
        b3_types = [(Int,Rational{Int}), (Float64,Float64), (BigInt,Rational{BigInt})]
        # b4 = BSplineBasis(4, Int128[0,1,2,3,4,5])
        b4_types = [(Int,Float64), (Float64,Float64), (Rational{Int},Rational{Int128})]
        # b5 = BSplineBasis(3, [-5.0, -3.5, -2.2, 1.0, 2.0, 3.6])
        b5_types = [(Int,Float64), (Float32,Float64), (Rational{Int},Float64)]
        # b6 = BSplineBasis(4, Rational{Int64}[-5//1, -7//2, -11//5, 1//1, 2//1, 18//5])
        b6_types = [(Int,Rational{Int64}), (Float32,Float32),
                    (Rational{Int},Rational{Int64})]
        # b7 = BSplineBasis(3, Int8[-6,-4,-2,0,0,1,2])
        b7_types = [(Int,Float64), (Float64,Float64), (Rational{Int16},Rational{Int16})]
        for (basis, types) = [(b1, b1_types), (b2, b2_types), (b3, b3_types),
                              (b4, b4_types), (b5, b5_types), (b6, b6_types),
                              (b7, b7_types)]
            leftknot = intervalindex(basis, 1)
            for (xtype, bspline_type) in types
                @test eltype(bsplines(basis, one(xtype), leftknot=leftknot)) == bspline_type
                @test eltype(bsplines(basis, one(xtype), AllDerivatives(2), leftknot=leftknot)) == bspline_type
                @test eltype(bsplines(basis, one(xtype), Derivative(0), leftknot=leftknot)) == bspline_type
            end
            for T = [Float32, Float64, BigFloat, Rational{Int128}]
                derivspace = Matrix{T}(undef, order(basis), order(basis))
                @test eltype(bsplines(basis, 1, AllDerivatives(2), leftknot=leftknot, derivspace=derivspace)) == T
                @test eltype(bsplines(basis, 1, Derivative(2), leftknot=leftknot, derivspace=derivspace)) == T
            end
        end

        function test_bsplines_nothing(basis, x)
            derivspace = Matrix{typeof(x)}(undef, order(basis), order(basis))
            @test bsplines(basis, x)                                                             === nothing
            @test bsplines(basis, x, leftknot=nothing)                                           === nothing
            @test bsplines(basis, x, AllDerivatives(2))                                          === nothing
            @test bsplines(basis, x, AllDerivatives(2), leftknot=nothing)                        === nothing
            @test bsplines(basis, x, AllDerivatives(2), derivspace=derivspace)                   === nothing
            @test bsplines(basis, x, AllDerivatives(2), leftknot=nothing, derivspace=derivspace) === nothing
            @test bsplines(basis, x, Derivative(2))                                              === nothing
            @test bsplines(basis, x, Derivative(2), leftknot=nothing)                            === nothing
            @test bsplines(basis, x, Derivative(2), derivspace=derivspace)                       === nothing
            @test bsplines(basis, x, Derivative(2), leftknot=nothing, derivspace=derivspace)     === nothing
        end
        test_bsplines_nothing(b1, NaN16)
        test_bsplines_nothing(b2, -Inf32)
        test_bsplines_nothing(b3, Inf64)
        test_bsplines_nothing(b4, 1//0)
        test_bsplines_nothing(b5, -1//0)
        test_bsplines_nothing(b6, -100)
        test_bsplines_nothing(b7, 10.0)
    end

    function test_bsplines(basis, basis_exact, x, isexact)
        NDERIV = 4
        leftknot = intervalindex(basis, x)
        range = leftknot-order(basis)+1:leftknot
        bsp0 = bsplines(basis, x, leftknot=leftknot)
        derivspace = Matrix{eltype(bsp0)}(undef, order(basis), order(basis))
        # AllDerivatives(NDERIVS)
        bspl_nderiv = bsplines(basis, x, AllDerivatives(NDERIV), leftknot=leftknot, derivspace=derivspace)
        bspl_nderiv_exact = bsplines_exact(basis_exact, x, range, AllDerivatives(NDERIV))
        if isexact
            @test bspl_nderiv == bspl_nderiv_exact
        else
            @test bspl_nderiv ≈ₑₗ bspl_nderiv_exact
        end
        @test bsp0 == view(bspl_nderiv,:,0)
        for N = 1:(NDERIV-1)
            # Derivative(N)
            @test bsplines(basis, x, Derivative(N), leftknot=leftknot, derivspace=derivspace) == view(bspl_nderiv,:,N)
            # AllDerivatives(N)
            if N == 1
                @test bsplines(basis, x, AllDerivatives(N), leftknot=leftknot) == OffsetArray(view(parent(bspl_nderiv),:,1:N), range, 0:N-1)
            else
                @test bsplines(basis, x, AllDerivatives(N), leftknot=leftknot, derivspace=derivspace) == OffsetArray(view(parent(bspl_nderiv),:,1:N), range, 0:N-1)
            end
        end
    end

    # b1 = BSplineBasis(1, 0:5)
    for x = Real[0, 1//2, 1.0f0, 4.1, Int8(5)]
        test_bsplines(b1, b1_exact, x, true)
    end
    # b2 = BSplineBasis(2, big(0):1:big(5))
    for x = Real[1.625, 2.6f0, big(5)]
        test_bsplines(b2, b2_exact, x, false)
    end
    for x = Real[22//7, big(9//2)]
        test_bsplines(b2, b2_exact, x, true)
    end
    # b3 = BSplineBasis(3, 0//1:5//1)
    for x = Real[0.123, 1.625f0, π, Float16(5.0)]
        test_bsplines(b3, b3_exact, x, false)
    end
    for x = Real[false, 1, 13//5]
        test_bsplines(b3, b3_exact, x, true)
    end
    # b4 = BSplineBasis(4, Int128[0,1,2,3,4,5])
    for x = Real[0, 2.6f0, 4.5, Int8(5)]
        test_bsplines(b4, b4_exact, x, false)
    end
    for x = Real[13//8]
        test_bsplines(b4, b4_exact, x, true)
    end
    # b5 = BSplineBasis(3, [-5.0, -3.5, -2.2, 1.0, 2.0, 3.6])
    for x = Real[-5, -4//1, -3.0f0, π, 3.6]
        test_bsplines(b5, b5_exact, x, false)
    end
    # b6 = BSplineBasis(4, Rational{Int64}[-5//1, -7//2, -11//5, 1//1, 2//1, 18//5])
    for x = Real[-3.0f0, 1.5]
        test_bsplines(b6, b6_exact, x, false)
    end
    for x = Real[-5, Int128(-2), Int8(-3), 18//5]
        test_bsplines(b6, b6_exact, x, true)
    end
    # b7 = BSplineBasis(3, Int8[-6,-4,-2,0,0,1,2])
    for x = Real[-6.0, -2.6f0, big(0.375), Int8(2)]
        test_bsplines(b7, b7_exact, x, false)
    end
    for x = Real[Rational{Int16}(0), 5//7]
        test_bsplines(b7, b7_exact, x, true)
    end

    leftknot = intervalindex(b5, 0)
    @test_throws ArgumentError bsplines(b5, 0, leftknot=leftknot-1)
    @test_throws ArgumentError bsplines(b5, 0, leftknot=leftknot+1)
    @test_throws ArgumentError bsplines(b5, 0, leftknot=nothing)
    @test_throws ArgumentError bsplines(b5, 0, Derivative(0), leftknot=leftknot-1)
    @test_throws ArgumentError bsplines(b5, 0, Derivative(1), leftknot=leftknot+1)
    @test_throws ArgumentError bsplines(b5, 0, Derivative(2), leftknot=nothing)
    @test_throws ArgumentError bsplines(b5, 0, AllDerivatives(1), leftknot=leftknot-1)
    @test_throws ArgumentError bsplines(b5, 0, AllDerivatives(2), leftknot=leftknot+1)
    @test_throws ArgumentError bsplines(b5, 0, AllDerivatives(3), leftknot=nothing)
    derivspace = Matrix{Float64}(undef, order(b5), order(b5))
    @test_throws ArgumentError bsplines(b5, 0, Derivative(0), derivspace=derivspace)
    @test_throws ArgumentError bsplines(b5, 0, AllDerivatives(1), derivspace=derivspace)
    @test_throws DimensionMismatch bsplines(b5, 0, Derivative(1), derivspace=zeros(order(b5)+1, order(b5)))
    @test_throws DimensionMismatch bsplines(b5, 0, AllDerivatives(2), derivspace=zeros(order(b5)))
    @test_throws DimensionMismatch bsplines(b5, 0, AllDerivatives(2), derivspace=zeros(order(b5), order(b5), 1))
    @test bsplines(b5, 0) == bsplines(b5, 0, leftknot=leftknot) == bsplines(b5, 0, Derivative(0), leftknot=leftknot)
    @test bsplines(b5, 0, Derivative(1)) == bsplines(b5, 0, Derivative(1), leftknot=leftknot) ==
                                            bsplines(b5, 0, Derivative(1), derivspace=derivspace) ==
                                            bsplines(b5, 0, Derivative(1), leftknot=leftknot, derivspace=derivspace)
    @test bsplines(b5, 0, AllDerivatives(3)) == bsplines(b5, 0, AllDerivatives(3), leftknot=leftknot) ==
                                                bsplines(b5, 0, AllDerivatives(3), derivspace=derivspace) ==
                                                bsplines(b5, 0, AllDerivatives(3), leftknot=leftknot, derivspace=derivspace)
end

@time @testset "bsplines!" begin
    function test_bsplines!_nothing(basis, x)
        dest     = Array{Float64}(undef, order(basis))
        dest_all = Array{Float32}(undef, order(basis), 2)
        derivspace     = Array{Float64}(undef, order(basis), order(basis))
        derivspace_all = Array{Float32}(undef, order(basis), order(basis))
        @test bsplines!(dest, basis, x) === nothing
        @test bsplines!(dest, basis, x, leftknot=nothing) === nothing
        @test bsplines!(dest_all, basis, x, AllDerivatives(2)) === nothing
        @test bsplines!(dest_all, basis, x, AllDerivatives(2), leftknot=nothing) === nothing
        @test bsplines!(dest_all, basis, x, AllDerivatives(2), derivspace=derivspace_all) === nothing
        @test bsplines!(dest_all, basis, x, AllDerivatives(2), leftknot=nothing, derivspace=derivspace_all) === nothing
        @test bsplines!(dest, basis, x, Derivative(2)) === nothing
        @test bsplines!(dest, basis, x, Derivative(2), leftknot=nothing) === nothing
        @test bsplines!(dest, basis, x, Derivative(2), derivspace=derivspace) === nothing
        @test bsplines!(dest, basis, x, Derivative(2), leftknot=nothing, derivspace=derivspace) === nothing
    end
    test_bsplines!_nothing(b1, NaN16)
    test_bsplines!_nothing(b2, -Inf32)
    test_bsplines!_nothing(b3, Inf64)
    test_bsplines!_nothing(b4, -100)
    test_bsplines!_nothing(b5, 10.0)

    function test_bsplines!(basis, basis_exact, T, x, isexact)
        NDERIV = 6
        leftknot = intervalindex(basis, x)
        offset = leftknot-order(basis)
        range = offset+1:leftknot
        derivspace = Matrix{T}(undef, order(basis), order(basis))

        # AllDerivatives(NDERIV)
        dest0 = fill(T(Inf), order(basis), NDERIV)
        bspl_nderiv_exact = bsplines_exact(basis_exact, x, range, AllDerivatives(NDERIV))
        bspl_nderiv = bsplines!(dest0, basis, x, AllDerivatives(NDERIV), leftknot=leftknot, derivspace=derivspace)
        if isexact
            @test bspl_nderiv == bspl_nderiv_exact
        else
            @test bspl_nderiv ≈ₑₗ bspl_nderiv_exact
        end
        @test parent(bspl_nderiv) === dest0

        # No derivative
        dest1 = fill(T(Inf), order(basis))
        bspl1 = bsplines!(dest1, basis, x, leftknot=leftknot)
        @test bspl1 == view(bspl_nderiv,:,0)
        @test parent(bspl1) === dest1

        # AllDerivatives(N) and Derivative(N) for N = 1:NDERIV-1
        for N = 1:NDERIV-1
            # Derivative(N)
            fill!(dest1, T(Inf))
            bspl1 = bsplines!(dest1, basis, x, Derivative(N), leftknot=leftknot, derivspace=derivspace)
            @test bspl1 == view(bspl_nderiv,:,N)
            @test parent(bspl1) === dest1
            # AllDerivatives(N)
            dest2 = fill(T(Inf), order(basis), N)
            if N == 1
                bspl2 = bsplines!(dest2, basis, x, AllDerivatives(N), leftknot=leftknot)
            else
                bspl2 = bsplines!(dest2, basis, x, AllDerivatives(N), leftknot=leftknot, derivspace=derivspace)
            end
            @test bspl2 == OffsetArray(view(parent(bspl_nderiv),:,1:N), offset, -1)
            @test parent(bspl2) === dest2
        end
    end

    # b1 = BSplineBasis(1, 0:5)
    for (T, x) = [(Float16, 0), (Float32, 0.5f0), (BigFloat, big(2)), (Rational{Int}, 5//2)]
        test_bsplines!(b1, b1_exact, T, x, true)
    end
    # b2 = BSplineBasis(2, big(0):1:big(5))
    for (T, x) = [(Rational{BigInt}, 8//9), (Rational{Int}, 22//7)]
        test_bsplines!(b2, b2_exact, T, x, true)
    end
    for (T, x) = [(Float32, 0.123), (Float64, 3//2), (BigFloat, 13//8)]
        test_bsplines!(b2, b2_exact, T, x, false)
    end
    # b3 = BSplineBasis(3, 0//1:5//1)
    for (T, x) = [(Rational{Int64}, 8//9), (Rational{Int32}, 1), (Rational{Int128}, 13//8)]
        test_bsplines!(b3, b3_exact, T, x, true)
    end
    for (T, x) = [(Float32, 0), (BigFloat, 2.6f0), (Float64, π)]
        test_bsplines!(b3, b3_exact, T, x, false)
    end
    # b4 = BSplineBasis(4, Int128[0,1,2,3,4,5])
    for (T, x) = [(Rational{Int}, 13//5), (Rational{BigInt}, 5)]
        test_bsplines!(b4, b4_exact, T, x, true)
    end
    for (T, x) = [(Float64, 1.625), (Float32, 9//2)]
        test_bsplines!(b4, b4_exact, T, x, false)
    end
    # b5 = BSplineBasis(3, [-5.0, -3.5, -2.2, 1.0, 2.0, 3.6])
    for (T, x) = [(Float32, -3.0f0), (Float64, Int8(-2)), (BigFloat, π)]
        test_bsplines!(b5, b5_exact, T, x, false)
    end
    # b6 = BSplineBasis(4, Rational{Int64}[-5//1, -7//2, -11//5, 1//1, 2//1, 18//5])
    for (T, x) = [(Rational{Int32}, -5//1), (Rational{Int64}, big(18//5))]
        test_bsplines!(b6, b6_exact, T, x, true)
    end
    for (T, x) = [(Float32, -4), (Float64, -1.0)]
        test_bsplines!(b6, b6_exact, T, x, false)
    end
    # b7 = BSplineBasis(3, Int8[-6,-4,-2,0,0,1,2])
    for (T,x) = [(Rational{Int16}, -6), (Rational{Int64}, big(1//8)),
                 (Rational{BigInt}, Int16(1))]
        test_bsplines!(b7, b7_exact, T, x, true)
    end
    for (T,x) = [(Float32, -9//2), (Float64, Int16(2))]
        test_bsplines!(b7, b7_exact, T, x, false)
    end

    leftknot = intervalindex(b5, 0)
    derivspace = zeros(order(b5), order(b5))
    @test_throws DimensionMismatch bsplines!(zeros(order(b5)+1), b5, 0)
    @test_throws DimensionMismatch bsplines!(zeros(order(b5)-1), b5, 0, leftknot=leftknot)
    @test_throws DimensionMismatch bsplines!(zeros(order(b5), 0), b5, 0, Derivative(1))
    @test_throws DimensionMismatch bsplines!(zeros(order(b5), 1), b5, 0, Derivative(1), leftknot=leftknot)
    @test_throws DimensionMismatch bsplines!(zeros(order(b5), 2), b5, 0, AllDerivatives(3))
    @test_throws DimensionMismatch bsplines!(zeros(order(b5), 3, 1), b5, 0, AllDerivatives(3), leftknot=leftknot)
    dest = zeros(order(b5))
    @test_throws ArgumentError bsplines!(dest, b5, 0, leftknot=leftknot-1)
    @test_throws ArgumentError bsplines!(dest, b5, 0, leftknot=leftknot+1)
    @test_throws ArgumentError bsplines!(dest, b5, 0, leftknot=nothing)
    @test_throws ArgumentError bsplines!(dest, b5, 0, Derivative(0), leftknot=leftknot-1)
    @test_throws ArgumentError bsplines!(dest, b5, 0, Derivative(1), leftknot=leftknot+1)
    @test_throws ArgumentError bsplines!(dest, b5, 0, Derivative(2), leftknot=nothing)
    @test_throws ArgumentError bsplines!(dest, b5, 0, Derivative(0), derivspace=derivspace)
    @test_throws DimensionMismatch bsplines!(dest, b5, 0, Derivative(1), derivspace=zeros(order(b5)+1, order(b5)))
    @test_throws DimensionMismatch bsplines!(dest, b5, 0, Derivative(2), derivspace=zeros(order(b5), order(b5), 1))
    dest = zeros(order(b5), 1)
    @test_throws ArgumentError bsplines!(dest, b5, 0, AllDerivatives(1), derivspace=derivspace)
    dest = zeros(order(b5), 3)
    @test_throws ArgumentError bsplines!(dest, b5, 0, AllDerivatives(3), leftknot=leftknot-1)
    @test_throws ArgumentError bsplines!(dest, b5, 0, AllDerivatives(3), leftknot=leftknot+1)
    @test_throws ArgumentError bsplines!(dest, b5, 0, AllDerivatives(3), leftknot=nothing)
    @test_throws DimensionMismatch bsplines!(dest, b5, 0, Derivative(1), derivspace=zeros(order(b5), order(b5)-1))
    @test_throws DimensionMismatch bsplines!(dest, b5, 0, Derivative(2), derivspace=zeros(order(b5)))
    dest1 = fill(NaN, order(b5))
    dest2 = fill(NaN, order(b5))
    bsp1 = bsplines!(dest1, b5, 0)
    bsp2 = bsplines!(dest2, b5, 0, leftknot=leftknot)
    @test bsp1 == bsp2
    @test parent(bsp1) === dest1
    @test parent(bsp2) === dest2
    fill!(dest1, NaN)
    fill!(dest2, NaN)
    dest3 = fill(NaN, order(b5))
    dest4 = fill(NaN, order(b5))
    bsp1 = bsplines!(dest1, b5, 0, Derivative(1))
    bsp2 = bsplines!(dest2, b5, 0, Derivative(1), leftknot=leftknot)
    bsp3 = bsplines!(dest3, b5, 0, Derivative(1), derivspace=derivspace)
    bsp4 = bsplines!(dest4, b5, 0, Derivative(1), leftknot=leftknot, derivspace=derivspace)
    @test bsp1 == bsp2 == bsp3 == bsp4
    @test parent(bsp1) === dest1
    @test parent(bsp2) === dest2
    @test parent(bsp3) === dest3
    @test parent(bsp4) === dest4
    dest1 = fill(NaN, order(b5), 3)
    dest2 = fill(NaN, order(b5), 3)
    dest3 = fill(NaN, order(b5), 3)
    dest4 = fill(NaN, order(b5), 3)
    bsp1 = bsplines!(dest1, b5, 0, AllDerivatives(3))
    bsp2 = bsplines!(dest2, b5, 0, AllDerivatives(3), leftknot=leftknot)
    bsp3 = bsplines!(dest3, b5, 0, AllDerivatives(3), derivspace=derivspace)
    bsp4 = bsplines!(dest4, b5, 0, AllDerivatives(3), leftknot=leftknot, derivspace=derivspace)
    @test bsp1 == bsp2 == bsp3 == bsp4
    @test parent(bsp1) === dest1
    @test parent(bsp2) === dest2
    @test parent(bsp3) === dest3
    @test parent(bsp4) === dest4
end

ext: jl


text:
abstract type CloudWatchLogsException <: Exception end

struct StreamNotFoundException <: CloudWatchLogsException
    stream::String
    group::String
    msg::Union{String, Nothing}
end

function Base.showerror(io::IO, exception::StreamNotFoundException)
    print(io, "Log stream ", exception.stream, " not found in group ", exception.group, ".")

    if exception.msg !== nothing
        print(io, " ", exception.msg)
    end
end

struct LogSubmissionException <: CloudWatchLogsException
    msg::String
end


text:
#f(x,y)=mean(KL(x,sigmoid(y))), where KL is the Kullback-Leibler divergence

FBinaryEntropyLossXsigmoidY(x,y)=([mean(x.*log(x)+(1.-x).*log(1.-x)+x.*y-log1pexp(y))],nothing)
FBinaryEntropyLossXsigmoidY_inplace(value,aux,x,y)=copy!(value,mean(x.*log(x)+(1.-x).*log(1.-x)+x.*y-log1pexp(y)))

DBinaryEntropyLossXsigmoidY[1]

function DBinaryEntropyLossXsigmoidY(derivativeIDX,f_c,faux_c,grad_c,grad_n,x,y)
    if derivativeIDX==1
        axpy!(grad_c[1]/length(x),y,grad_n)
    elseif derivativeIDX==2
        axpy!(grad_c[1]/length(x),x-sigmoid(y)),grad_n)
    end
end

Derivative[FBinaryEntropyLossXsigmoidY]=DBinaryEntropyLossXsigmoidY
Inplace[FBinaryEntropyLossXsigmoidY]=DBinaryEntropyLossXsigmoidY_inplace

BinaryKullbackLeiblerLossXsigmoidY(nx,ny)=ADnode(FbinaryEntropyLoss,[nx ny])
export BinaryKullbackLeiblerLossXsigmoidY

#TODO: GPU

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff


text:
#=
Copyright 2017 INSIGNEO Institute for in silico Medicine

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

   http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
=#

# A model in `openBF` consists of a collection of vessels,
# a heart, and a description of blood mechanical properties.
# Vessels should be defined to be of type `Vessel`, the inlet
# boundary condition is of type `Heart`, and the blood properties
# are contained in a `Blood` instance.
module BTypes

# Instances of `BTypes` data structures can be created as
#
#     BTypes.Vessel(args...)
#     BTypes.Heart(args...)
#     BTypes.Blood(args...)
#
# where the `args...` are specified for each data structure in the
# following. Functions in [initialise.jl](initialise.html) generate
# instances of `BTypes` data structures.
export Vessel, Heart, Blood

# Each vessel in the arterial system is represented by a single
# instance of type `Vessel`.
#
# --------------------------------------------------------------------
# Arguments:
# ---------------- ---------------------------------------------------
# `label`          Vessel clinical name. It is used in
#                  [`initialiseVessel`](initialise.html#initialiseVessel)
#                  to generate output file names.
#
# `ID`             Vessel ID within the `Graph` structure defined
#                  in [main.jl](main.html#grafo). Unique integer.
#
# `sn`             Source node ID within the `Graph` structure. Unique
#                  integer.
#
# `tn`             Terminal node ID within the `Graph` structure. Unique
#                  integer.
#
# `M`              Number of nodes along the vessel. This is user defined
#                  in the project `project_model.csv` file.
#
# `dx`             $\Delta x$ gives the spatial discretisation and it
#                  is computed once in [initialise.jl](initialise.html#dx).
#
# `Ccfl`           Courant-Friedrichs-Lewy (CFL) number. It is used to
#                  compute the $\Delta t$ by
#                  [`calculateDeltaT`](godunov.html#calculateDeltaT)
#                  at each time step. The CFL is always in the
#                  interval $[0,1]$.
#
# `beta`           Elastic constant used to compute the trans-mural
#                  [`pressure`](converter.html#pressure) `P`.
#
# `gamma`          Elastic constant used by
#                  [`waveSpeed`](converter.html#waveSpeed)
#                  to compute `c`.
#
# `R0`             Unstressed lumen radius.
#
# `A0`             Unstressed cross-sectional area.
#
# `A`              Cross-sectional area at each node along the vessel.
#
# `Q`              Flow-rate at each node along the vessel.
#
# `u`              Axial velocity at each node along the vessel.
#
# `c`              Wave speed at each node along the vessel.
#
# `P`              Trans-mural pressure at each node along the vessel.
#
# `W1M0`           Backward Riemann invariant at $x=L$.
#
# `W2M0`           Forward Riemann invariant.
#
# `U00A`           `A` in the ghost cell between 0 and -1. Ghost cells are
#                  used in [godunov.jl](godunov.html) to compute correctly
#                  the solution at boundaries cells.
#
# `U01A`           `A` in the ghost cell between -1 and -2.
#
# `UM1A`           `A` in the ghost cell between `M` and `M+1`.
#
# `UM2A`           `A` in the ghost cell between `M+1` and `M+2`. Ghost
#                  cells are defined for `Q` in the same manner.
#
# `temp_P_name`    Name of temporary file where to store `P` values
#                  while running the solution. Temporary files (`.temp`) are
#                  created for all the variables (`Q`, `A`, `u`, `c`).
#                  Temporary files are handled by
#                  functions defined in [IOutilities.jl](IOutilities.html).
#
# `out_P_name`     Name of final output file where to store `P` values.
#                  Every variable has its own `.out` file.
#
# `temp_P`         Julia handle to direct the output writing operation.
#                  Every variable has its own `IOStream`.
#
# `Rt`             Outlet reflection coefficient. See
#                  [boundary_conditions.jl](boundary_conditions.html).
#
# `R1`             Three element windkessel proximal resistance.
#
# `R2`             Peripheral resistance.
#
# `Cc`             Peripheral compliance.
#
# `Pc`             Pressure through peripheral compliance.
# --------------------------------------------------------------------
#
type Vessel
  label :: String

  #Topological notation
  ID :: Int64
  sn :: Int64
  tn :: Int64
  rn :: Int64

  #Numerical constants
  M       :: Int64
  dx      :: Float64
  invDx   :: Float64
  halfDx  :: Float64
  Ccfl    :: Float64

  #Physical constants
  beta  :: Array{Float64, 1}
  gamma :: Array{Float64, 1}
  gamma_ghost :: Array{Float64, 1}
  A0    :: Array{Float64, 1}
  dA0dx :: Array{Float64, 1}
  dTaudx :: Array{Float64, 1}
  Pext  :: Float64

  #Iterative solution
  A :: Array{Float64, 1}
  Q :: Array{Float64, 1}
  u :: Array{Float64, 1}
  c :: Array{Float64, 1}
  P :: Array{Float64, 1}

  #Riemann invariants
  W1M0 :: Float64
  W2M0 :: Float64

  #Ghost cells
  U00A :: Float64
  U00Q :: Float64
  U01A :: Float64
  U01Q :: Float64

  UM1A :: Float64
  UM1Q :: Float64
  UM2A :: Float64
  UM2Q :: Float64

  #Temporary files name
  temp_P_name :: String
  temp_Q_name :: String
  temp_A_name :: String
  temp_c_name :: String
  temp_u_name :: String

  last_P_name :: String
  last_Q_name :: String
  last_A_name :: String
  last_c_name :: String
  last_u_name :: String

  #Output files name
  out_P_name :: String
  out_Q_name :: String
  out_A_name :: String
  out_c_name :: String
  out_u_name :: String

  #Temporary files IOstreams
  temp_P :: IOStream
  temp_Q :: IOStream
  temp_A :: IOStream
  temp_c :: IOStream
  temp_u :: IOStream

  last_P :: IOStream
  last_Q :: IOStream
  last_A :: IOStream
  last_c :: IOStream
  last_u :: IOStream

  #Saving locations
  node2 :: Int64
  node3 :: Int64
  node4 :: Int64

  #Peripheral boundary condition
  Rt :: Float64

  R1 :: Float64
  R2 :: Float64
  Cc :: Float64
  Pc :: Float64

  #Slope
  slope :: Array{Float64, 1}

  #MUSCLArrays
  flux :: Array{Float64, 2}
  uStar :: Array{Float64, 2}

  vA :: Array{Float64, 1}
  vQ :: Array{Float64, 1}

  dU :: Array{Float64, 2}

  slopesA :: Array{Float64, 1}
  slopesQ :: Array{Float64, 1}

  Al :: Array{Float64, 1}
  Ar :: Array{Float64, 1}

  Ql :: Array{Float64, 1}
  Qr :: Array{Float64, 1}

  Fl :: Array{Float64, 2}
  Fr :: Array{Float64, 2}

end

# The inlet boundary condition is applied to the first vessel in the
# model. A `Heart` instance contains all inlet information.
#
# --------------------------------------------------------------------
# Arguments:
# ----------------- --------------------------------------------------
# `BC_switch`       Integer flag to select inlet boundary condition.
#                   It takes values from 1 to 4 and it is chosen by
#                   the user in the `project_constants.jl` file.
#
# `cardiac_T`       Cardiac period.
#
# `sys_T`           Length of systole.
#
# `initial_flow`    `Q` (ml/s) initial value to be set within all the vessels.
#
# `flow_amplitude`  `Q` maximum amplitude (ml/s) for the user-defined
#                   inlet boundary condition.
#
# `input_data`      Matrix defined as
#                   `[time::Array, flow::Array]`
#                   to be used as inlet boundary condition.
# --------------------------------------------------------------------
#
type Heart
  BC_switch  :: Int64
  inlet_type :: String

  cardiac_T :: Float64
  sys_T     :: Float64

  initial_flow   :: Float64
  flow_amplitude :: Float64

  input_data :: Array{Float64, 2}
end

# Blood mechanical properties.
#
# --------------------------------------------------------------------
# Arguments:
# -------  -----------------------------------------------------------
# `mu`     Dynamic viscosity (poise)
#
# `rho`    Density (kg/m$^3$)
#
# `Cf`     Viscous loss coefficient (see
#          [`loadGlobalConstants`](initialise.html#loadGlobalConstants)
#          and [`source`](godunov.html#source)).
# --------------------------------------------------------------------
#
type Blood
  mu  :: Float64
  rho :: Float64
  Cf  :: Float64
  gamma_profile :: Float64
end

#
end


text:
using mjWrap

################################################## Experiment Centric Functions

## accepts model,
## current and next states,
## current and next observations,
## applied controls that took us from current to next state

function rewardfunction(x::mjWrap.mjSet, tid::Int,
                                    s0::Vector{Float64}, # pass in states
                                    s1::Vector{Float64},
                                    o0::Vector{Float64}, # observations
                                    o1::Vector{Float64},
                                    ctrl::Vector{Float64}) # and controls
   pos0 = s0[1]
   pos1 = s1[1]
   height = s1[2]
   angle = s1[3]

   reward = (pos1 - pos0) / (x.dt*x.skip) # forward x direction

   upright_bonus = 3.0
   t_height = 0.8
   if height > t_height
      reward += upright_bonus
   end
   #reward -= upright_bonus * (height - t_height)^2
   reward -= 1e-2 * sum(ctrl.^2)

   return reward
end

function obsfunction(x::mjWrap.mjSet, tid::Int,
                                 s::Vector{Float64}, o::Vector{Float64})
   o[2] = mjWrap.wraptopi(o[2]) # global angle

   #o[:] = max(o[:], -10.0)
   #o[:] = min(o[:],  10.0) # clipped
   #scales = [0.5, 0.1, 0.5, 0.1, 0.5, 1.0, 1.0, 1.0, 3.0, 3.0, 3.0]

   #for i=1:length(o)
   #    o[i] = o[i] / scales[i]
   #end
   #o[:] = o./scales
   #return o
end

## accepts model,
## default state,
## preallocated matrix of states (nqnv x numT)
function initfunction!(x::mjWrap.mjSet,
                                    s::mjWrap.TrajSamples)
   numT = size(s.state, 3)

   nq, nv, nu, ns = mjWrap.modelparams(x)
   s0 = s.s0

   sTemp = zeros(nq+nv)
   c = zeros(nu)
   o = zeros(ns)
   tid = 1

   diverse = true #false

   init_state = view( s.state, :, 1, : )
   for t=1:numT
      init_state[:,t] .= s0
      if isodd(t) && diverse
         # print("D ")
         angle = pi
         init_state[2,t]            = s0[2]           +rand()  * 0.1 # - 0.05
         init_state[3,t]            = s0[3]           +rand()  * angle - angle/2.0

         mjWrap.reset(x, tid, init_state[:,t], c, o)
         for i=1:50
            mjWrap.step!(x, tid, c, sTemp, o)
         end

         if (init_state[2,t] < 1)
            println("Hopper height: ", init_state[2,t])
         end
         # s.state[:,t] = sTemp
      else
         # print("N ")
         s.state[2,t]            = s0[2]           +rand()  * 0.1 - 0.05
         init_state[3,t]          = s0[3]           +rand()  * 1.0 - 0.5
         init_state[(1+nq):end,t] = s0[(1+nq):end]+rand(nv) * 0.01 - 0.005
         if (init_state[2,t] < 1)
            println("Hopper height: ", init_state[2,t])
         end
      end
   end
end

function modelfunction!(x::mjWrap.mjSet, iter::Int=0)
end

## accepts model and ndarray of control noise to be applied at sample time
## Array is: nu x T x numT , same size as ctrl
function ctrlfunction!(x::mjWrap.mjSet, ctrl::Array{Float64})
   randn!(ctrl)
end

## evaluation function: reward function shapes exploration, evaluation
## function defines success or failure
function evalfunction(x::mjWrap.mjSet, d::mjWrap.TrajSamples)
   # swimmer success: max distance along x axis covered at end of T
   xaxis_index = 1
   return mean(d.state[xaxis_index,end,:]) # x position at time T, for all numT
end

using ExpFunctions
myfuncs = FunctionSet(modelfunction!,
                      initfunction!,
                      ctrlfunction!,
                      obsfunction,
                      rewardfunction,
                      evalfunction)


text:
module FrostHelperNoMovementTrigger

using ..Ahorn, Maple

@mapdef Trigger "FrostHelper/NoMovementTrigger" NoMovementTrigger(x::Integer, y::Integer, width::Integer=16, height::Integer=16)

const placements = Ahorn.PlacementDict(
	"No Movement (Frost Helper)" => Ahorn.EntityPlacement(
		NoMovementTrigger,
		"rectangle"
	)
)

end

max_stars_repo_name: JaThePlayer/FrostHelper

max_issues_repo_name: JaThePlayer/FrostHelper

max_forks_repo_name: JaThePlayer/FrostHelper


text:
"""
getΩS(As::Vector{CompoundFIDType{T,SST}}) where {T,SST}

returns ΩS::Vector{Vector{Vector{T}}}
"""
function getΩS(As::Vector{SHType{T}}) where T

    ΩS = Vector{Vector{Vector{T}}}(undef, length(As))

    for n = 1:length(As)

        ΩS[n] = Vector{Vector{T}}(undef, length(As[n].Ωs) + length(As[n].Ωs_singlets))
        for i = 1:length(As[n].Ωs)

            ΩS[n][i] = copy(As[n].Ωs[i])

        end

        for i = 1:length(As[n].Ωs_singlets)
            ΩS[n][length(As[n].Ωs)+i] = [ As[n].Ωs_singlets[i]; ]
        end
    end

    return ΩS
end

"""
getPs( ΩS::Vector{Vector{Vector{T}}}, hz2ppmfunc) where T <: Real

returns Ps::Vector{Vector{Vector{T}}}
"""
function getPs( ΩS::Vector{Vector{Vector{T}}}, hz2ppmfunc) where T <: Real

    N_compounds = length(ΩS)

    Ps = Vector{Vector{Vector{T}}}(undef, N_compounds)
    for n = 1:N_compounds

        Ps[n] = Vector{Vector{T}}(undef, length(ΩS[n]))
        for i = 1:length(ΩS[n])

            Ps[n][i] = Vector{T}(undef, length(ΩS[n][i]))
            for l = 1:length(ΩS[n][i])

                Ps[n][i][l] = hz2ppmfunc( ΩS[n][i][l]/(2*π) )
            end
        end
    end

    return Ps
end

#### for cost function.

# getPsnospininfo
function getPsnospininfo(As::Vector{NMRSpectraSimulator.SHType{T}}, hz2ppmfunc) where T

    ΩS_ppm = Vector{Vector{T}}(undef, length(As))

    for (n,A) in enumerate(As)

        ΩS_ppm[n] = hz2ppmfunc.( combinevectors(A.Ωs) ./ (2*π) )

        tmp = hz2ppmfunc.( A.Ωs_singlets ./ (2*π) )
        push!(ΩS_ppm[n], tmp...)
    end

    return ΩS_ppm
end

function combinevectors(x::Vector{Vector{T}})::Vector{T} where T

    if isempty(x)
        return Vector{T}(undef, 0)
    end

    N = sum(length(x[i]) for i = 1:length(x))

    y = Vector{T}(undef,N)

    st_ind = 0
    fin_ind = 0
    for i = 1:length(x)
        st_ind = fin_ind + 1
        fin_ind = st_ind + length(x[i]) - 1

        y[st_ind:fin_ind] = x[i]
    end

    return y
end

function isnumericallyclose(x::T, y::T, tol = eps(T)*2) where T
    if abs(x-y) < tol
        return true
    end

    return false
end

"""
    convertcompactdomain(x::T, a::T, b::T, c::T, d::T)::T

converts compact domain x ∈ [a,b] to compact domain out ∈ [c,d].
"""
function convertcompactdomain(x::T, a::T, b::T, c::T, d::T)::T where T <: Real

    return (x-a)*(d-c)/(b-a)+c
end

function convertcompactdomain(x::Vector{T}, a::T, b::T, c::T, d::T)::Vector{T} where T <: Real

    return collect( convertcompactdomain(x[i], a, b, c, d) for i = 1:length(x) )
end

### for plotting.
"""
prunelowsignalentries(x, threshold_factor::T, reduction_factor::Int)::Vector{Int} where T

keep the high signal portion and the reduced low signal portion.
keep_threshold = maximum(abs_x) * threshold_factor.
"""
function prunelowsignalentries(x, threshold_factor::T, reduction_factor::Int) where T

    abs_x = abs.(x)

    # find indices that have a low signal.
    keep_threshold = maximum(abs_x) * threshold_factor
    inds_for_reducing = findall(xx->xx<keep_threshold, abs_x)
    reduced_inds = inds_for_reducing[1:reduction_factor:length(inds_for_reducing)]

    # find indices that have a high signal.
    keep_inds = findall(xx->xx>keep_threshold, abs_x)

    # keep the high signal portion and the reduced low signal portion.
    inds = sort([keep_inds; reduced_inds])

    return inds, keep_inds, inds_for_reducing
end

### for reading config files.
function extractinfofromconfig( config_path::String,
    molecule_names::Vector{String}) where T <: Real

    N_compounds = length(molecule_names)

    # read.
    file_strings = readlines(config_path)

    cs_delta_group = Vector{Vector{Vector{Float64}}}(undef, N_compounds)
    #λ_group_labels = Vector{Vector{Vector{Int}}}(undef, N_compounds)

    for n = 1:N_compounds

        cs_delta_group[n] = extractmoleculeinfofromconfig(file_strings,
                                    molecule_names[n], config_path)
    end

    return cs_delta_group#, λ_group_labels
end

function extractmoleculeinfofromconfig( file_strings,
                                        name::String,
                                        file_label::String) where T <: Real

    # find the block of text for the input name.
    text_buffer = findblock(file_strings, name, "end compound", file_label)
    if isempty(text_buffer)
        println("Error processing config file.")
        return Vector{Vector{Float64}}(undef, 0), Vector{Vector{Int}}(undef, 0)
    end

    # get the number of groups.
    N_groups = length(filter(xx->occursin("Group",xx), text_buffer))

    # get default set up.
    cs_delta_group = Vector{Vector{Float64}}(undef, N_groups)
    #λ_group_labels = Vector{Vector{Int}}(undef, N_groups)

    for i = 1:N_groups

        # update position.
        cs_buffer = findblock(text_buffer, "Group $(i)", "end group", "Group $(i)")
        if isempty(cs_buffer)
            return Vector{Vector{Float64}}(undef, 0), Vector{Vector{Int}}(undef, 0)
        end
        N_cs = length(cs_buffer)

        cs_delta_group[i] = Vector{Float64}(undef, N_cs)
        #λ_group_labels[i] = Vector{Int}(undef, N_cs)

        for k = 1:length(cs_buffer)
            tokens = split(cs_buffer[k], (':',','))
            cs_delta_group[i][k] = tryparse(Float64, tokens[2])
            #λ_group_labels[i][k] = tryparse(Int, tokens[3])
        end
    end

    return cs_delta_group#, λ_group_labels
end

"""
uniqueinds(a_in::Vector{T}; atol = 1e-6) where T

Returns the unique values of `a_in`, with absolute tolerance `atol`, and the indices for each unique value.
"""
function uniqueinds(a_in::Vector{T}; atol = 1e-6) where T

    if length(a_in) < 2
        inds = Vector{Vector{Int}}(undef, length(a_in))

        if length(a_in) == 1
            inds[1] = Vector{Int}(undef, 1)
            inds[1][1] = 1
        end

        return copy(a_in), inds
    end

    a = copy(a_in)
    inds_a = collect(1:length(a))

    b = Vector{T}(undef, 0)
    inds_b = Vector{Vector{Int}}(undef, 0)

    while !isempty(a)

        target = pop!(a)
        inds = findall(xx->isapprox(target, xx; atol = atol), a)

        target_ind = pop!(inds_a)

        # save to output.
        push!(b, target)
        push!(inds_b, [target_ind; inds_a[inds]])

        # delete saved.
        deleteat!(a, inds)
        deleteat!(inds_a, inds)
    end

    return b, inds_b
end

"""
getpairs(inds::Vector{T})

get all exhaustive pairwise combos without symmetry of the 1D array `inds`.
"""
function getpairs(inds::Vector{T}) where T

    out = Vector{Tuple{T,T}}(undef, 0)
    for i = 1:length(inds)
        for j = i+1:length(inds)
            push!(out, (inds[i], inds[j]))
        end
    end

    return out
end

"""
isallsame(a::Vector{T}; atol::T = 1e-6) where T

returns true if the entries in `a` are all within an abolute tolerance of `atol`.
"""
function isallsame(a::Vector{T}; atol::T = 1e-6) where T
    if length(findall(xx->isapprox(a[1], xx; atol = atol), a)) == length(a)
        return true
    end

    return false
end

"""
Returns the edges of a connected path as specified by `vertices`.
"""
function getconnectpath(vertices::Vector{Int})
    N = length(vertices)

    out = Vector{Tuple{Int,Int}}(undef, N-1)

    for i = 1:N-1
        out[i] = (vertices[i], vertices[i+1])
    end

    return out
end

"""
removetargetintegers(C::Vector{Vector{Int}}, search_list::Vector{Int})

remove an entry from C if the any of its integer array values appear in any entry of `search_list`.
"""
function removetargetintegers(C::Vector{Vector{Int}}, search_list::Vector{Int})

    keep_flags = trues(length(C))
    for i = 1:length(C)

        if any(C[i][1] .== search_list)
            keep_flags[i] = false
        end
    end

    return C[keep_flags]
end

"""
keeptargetintegers(C::Vector{Vector{Int}}, search_list::Vector{Int})

keep an entry from C if the any of its integer array values appear in any entry of `search_list`.
"""
function keeptargetintegers(C::Vector{Vector{Int}}, search_list::Vector{Int})

    keep_flags = falses(length(C))
    for i = 1:length(C)

        if any(C[i][1] .== search_list)
            keep_flags[i] = true
        end
    end

    return C[keep_flags]
end


text:
using DataFrames, Gadfly

set_default_plot_size(3.3inch, 3.3inch)

srand(123)
D = convert(DataFrame, 99*rand(4, 4)+0.5)

xsc  = Scale.x_continuous(minvalue=0.0, maxvalue=100)
ysc  = Scale.y_continuous(minvalue=0.0, maxvalue=100)

plot(D, x=:x1, y=:x2, xend=:x3, yend=:x4, Geom.segment(arrow=true), xsc, ysc)
plot(D, x=:x1, y=:x2, xend=:x3, yend=:x4, Geom.vector, xsc, ysc)

plot(z=(x,y)->x*exp(-(x^2+y^2)),
        xmin=[-2], xmax=[2], ymin=[-2], ymax=[2],
        Geom.vectorfield(scale=0.4, samples=17),
        Scale.x_continuous(minvalue=-2.0, maxvalue=2.0),
        Scale.y_continuous(minvalue=-2.0, maxvalue=2.0),
        Guide.xlabel("x"), Guide.ylabel("y"), Guide.colorkey(title="z")
    )

Z = fill(1.0, 5, 5)
Z[2:4, 2:4] = 2.0
Z[3,3] = 3.0

plot(z=Z,
        Geom.vectorfield,
        Scale.x_continuous(minvalue=1.0, maxvalue=5.0),
        Scale.y_continuous(minvalue=1.0, maxvalue=5.0),
        Guide.xlabel("x"), Guide.ylabel("y"),
        Theme(key_position=:none)
    )


text:
function test{T,N}(x::Array{T,N})
  return N, T, size(x)
end

println(test([1,2,3,4]))
println(test([1 2 3 4]))
println(test([1 2; 3 4]))
println(test(zeros(Int, 3, 3, 3)))

function addone(x)
  y = zeros(Int, length(x))
  for i in 1:length(x)
    y[i] = x[i] + 1
  end
  return y
end

println(addone([1,2 ,3]))


text:
@fastmath @inline tripleToInt(i::UInt32, j::UInt32, k::UInt32, nB::UInt32, nB2::UInt32)::UInt32 = UInt32(i*nB2 + j*nB + k)
@fastmath @inline tripleToInt(i::UInt32, j::UInt32, k::UInt32, offset::UInt32, nB::UInt32, nB2::UInt32)::UInt32 = UInt32((i+offset)*nB2 + (j+offset)*nB + (k+offset))

tripleToInt(i,j,k, nB, nB2) = tripleToInt(UInt32(i), UInt32(j), UInt32(k), UInt32(nB), UInt32(nB2))
tripleToInt(i,j,k, offset,nB,nB2) = tripleToInt(UInt32(i), UInt32(j), UInt32(k), UInt32(offset), UInt32(nB), UInt32(nB2))
tripleToInt(i,j,k,offset,nB,nB2) = tripleToInt(i+offset,j+offset,k+offset,nB,nB2)

@fastmath @inline function intToTriple(::Type{T}, z::UInt32) where {T<:Integer}
    r,k = divrem(z,nB)
    i,j = divrem(r,nB)
    return (convert(T,i)-nBh,convert(T,j)-nBh,convert(T,k)-nBh)
end
intToTriple(z::UInt32) = intToTriple(Int64, z)


text:
# ------------------------ #
#  -- Robust Regression -- #
# ------------------------ #
# ->     r  = Xθ - y
# ->   ρ(r) is the robust penalty assoc with residual
# ->   ψ(r) = ρ'(r) (first deriv)
# ->   ϕ(r) = ψ'(r) (second deriv; may be discont)
# ->   Λ(r) = diag(ϕ(r))
# ->   f(θ) = ∑ρ.(r) + λ|θ|₂²
# ->  ∇f(θ) = X'ψ.(r) + λθ
# -> ∇²f(θ) = X'Λ(r)X + λI
# ---------------------------------------------------------

function fgh!(glr::GLR{RobustLoss{ρ},<:L2R}, X, y, scratch
              ) where ρ <: RobustRho1P{δ} where δ
    n, p = size(X)
    λ    = get_penalty_scale(glr, n)
    ψ_   = ψ(ρ)
    ϕ_   = ϕ(ρ)
    if glr.fit_intercept
        (f, g, H, θ) -> begin
            r  = scratch.n
            get_residuals!(r, X, θ, y)
            w  = scratch.n2
            w .= convert.(Float64, abs.(r) .<= δ)
            # gradient via ψ function
            g === nothing || begin
                ψr  = scratch.n3
                ψr .= ψ_.(r, w)
                apply_Xt!(g, X, ψr)
                g .+= λ .* θ
                glr.penalize_intercept || (g[end] -= λ * θ[end])
            end
            # Hessian via ϕ function
            H === nothing || begin
                # NOTE: Hessian allocates a ton anyway so use of scratch is a
                # bit pointless
                ϕr = ϕ_.(r, w)
                ΛX = ϕr .* X
                mul!(view(H, 1:p, 1:p), X', ΛX)
                ΛXt1 = sum(ΛX, dims=1)
                @inbounds for i in 1:p
                    H[i, end] = H[end, i] = ΛXt1[i]
                end
                H[end, end] = sum(ϕr)
                add_λI!(H, λ, glr.penalize_intercept)
            end
            # function value
            f === nothing || return glr.loss(r) + glr.penalty(view_θ(glr, θ))
        end
    else
        (f, g, H, θ) -> begin
            r = scratch.n
            get_residuals!(r, X, θ, y)
            w = scratch.n2
            w .= convert.(Float64, abs.(r) .<= δ)
            # gradient via ψ function
            g === nothing || begin
                ψr  = scratch.n3
                ψr .= ψ_.(r, w)
                apply_Xt!(g, X, ψr)
                g .+= λ .* θ
            end
            # Hessian via ϕ function
            H === nothing || (mul!(H, X', ϕ_.(r, w) .* X); add_λI!(H, λ))
            f === nothing || return glr.loss(r) + glr.penalty(θ)
        end
    end
end

function Hv!(glr::GLR{RobustLoss{ρ},<:L2R}, X, y, scratch
             ) where ρ <: RobustRho1P{δ} where δ
    n, p = size(X)
    λ    = get_penalty_scale(glr, n)
    ϕ_   = ϕ(ρ)
    # see d_logistic.jl for more comments on this (similar procedure)
    if glr.fit_intercept
        (Hv, θ, v) -> begin
            r  = scratch.n
            get_residuals!(r, X, θ, y)
            w  = scratch.n2
            w .= convert.(Float64, abs.(r) .<= δ)
            w .= ϕ_.(r, w)
            # views on first p rows (intercept row treated after)
            a    = 1:p
            Hvₐ  = view(Hv, a)
            vₐ   = view(v, a)
            XtΛ1 = view(scratch.p, a)     # we can recycle
            apply_Xt!(XtΛ1, X, w)
            vₑ   = v[end]
            # update for first p rows
            t    = scratch.n3
            apply_X!(t, X, vₐ)
            t  .*= w
            apply_Xt!(Hvₐ, X, t)
            Hvₐ .+= λ .* vₐ .+ XtΛ1 .* vₑ
            # update for the last row (intercept)
            Hv[end] = dot(XtΛ1, vₐ) + (sum(w) + λ_if_penalize_intercept(glr, λ)) * vₑ
        end
    else
        (Hv, θ, v) -> begin
            r  = scratch.n
            get_residuals!(r, X, θ, y)
            w  = scratch.n2
            w .= convert.(Float64, abs.(r) .<= δ)
            w .= ϕ_.(r, w)
            t  = scratch.n3
            apply_X!(t, X, v)
            t .*= w
            apply_Xt!(Hv, X, t)
            Hv .+= λ .* v
        end
    end
end

# For IWLS
function Mv!(glr::GLR{RobustLoss{ρ},<:L2R}, X, y, scratch;
             threshold=1e-6) where ρ <: RobustRho1P{δ} where δ
    n, p = size(X)
    λ    = get_penalty_scale(glr, n)
    ω_   = ω(ρ, threshold)
    # For one θ, we get one system of equation to solve
    # which we solve via an iterative method so, one θ
    # gives one way of applying the relevant matrix (X'ΛX+λI)
    (ωr, θ) -> begin
        r   = scratch.n
        get_residuals!(r, X, θ, y)
        w   = scratch.n2
        w  .= convert.(Float64, abs.(r) .<= δ)
        # ω = ψ(r)/r ; weighing factor for IWLS
        ωr .= ω_.(r, w)
        # function defining the application of (X'ΛX + λI)
        if glr.fit_intercept
            (Mv, v) -> begin
                a    = 1:p
                vₐ   = view(v, a)
                Mvₐ  = view(Mv, a)
                XtW1 = view(scratch.p, a)
                @inbounds for j in a
                    XtW1[j] = dot(ωr, view(X, :, j))
                end
                vₑ = v[end]
                t  = scratch.n
                apply_X!(t, X, vₐ)
                t .*= ωr
                mul!(Mvₐ, X', t)
                Mvₐ .+= λ .* vₐ .+ XtW1 .* vₑ
                Mv[end] = dot(XtW1, vₐ) +
                            (sum(ωr) + λ_if_penalize_intercept(glr, λ)) * vₑ
            end
        else
            (Mv, v) -> begin
                t  = scratch.n
                apply_X!(t, X, v)
                t .*= ωr
                mul!(Mv, X', t)
                Mv .+= λ .* v
            end
        end
    end
end

# this is a bit of an abuse in that in some cases the ρ is not everywhere
# differentiable
function smooth_fg!(glr::GLR{RobustLoss{ρ},<:ENR}, X, y, scratch
                    ) where ρ <: RobustRho1P{δ} where δ
    n, p = size(X)
    λ    = get_penalty_scale_l2(glr, n)
    ψ_   = ψ(ρ)
    (g, θ) -> begin
        r   = scratch.n
        get_residuals!(r, X, θ, y)
        w   = scratch.n2
        w  .= convert.(Float64, abs.(r) .<= δ)
        ψr  = scratch.n3
        ψr .= ψ_.(r, w)
        apply_Xt!(g, X, ψr)
        g .+= λ .* θ
        glr.fit_intercept && (glr.penalize_intercept || (g[end] -= λ * θ[end]))
        return glr.loss(r) + get_l2(glr.penalty)(view_θ(glr, θ))
    end
end


text:
# θ_t - kąt pomiędzy udem a osią x 0-45
# θ_k kąty pomiędzy tułowiem a osią x # 120 to 60
# kąty obrotu wiosla 30-120
# User interface
struct Boat_settings
    angle_ranges::Boat_angle_ranges
    timing::Boat_timing
    x1va::Float64
    x1air::Float64
end

Boat_settings(γ_minmax, θ_t_minmax, θ_k_minmax, timing, x1va, x1air) =
    Boat_settings(Boat_angle_ranges(γ_minmax, θ_t_minmax, θ_k_minmax), timing, x1va, x1air)

Boat_settings() = Boat_settings(Boat_angle_ranges(), Boat_timing(), 0.0, 0.0)


text:
@with_kw mutable struct AsynchronousData{T <: AbstractFloat}
    timestamp::Int = 1
    incumbent::Int = 1
    active::Int = 0
end

"""
    AsynchronousExecution

Functor object for using synchronous execution in an L-shaped algorithm (assuming multiple Julia cores are available). Create by supplying a [`Asynchronous`](@ref) object through `execution` in the `LShapedSolver` factory function and then pass to a `StochasticPrograms.jl` model.

"""
struct AsynchronousExecution{F <: AbstractFeasibility,
                             T <: AbstractFloat,
                             A <: AbstractVector,
                             S <: LQSolver} <: AbstractExecution
    data::AsynchronousData{T}
    subobjectives::Vector{A}
    model_objectives::Vector{A}
    finished::Vector{Int}
    subworkers::Vector{SubWorker{F,T,A,S}}
    work::Vector{Work}
    finalize::Vector{Work}
    metadata::Vector{MetaData}
    decisions::Decisions{A}
    cutqueue::CutQueue{T}
    active_workers::Vector{Future}
    triggered::Vector{Bool}
    added::Vector{Bool}
    max_active::Int
    κ::T

    function AsynchronousExecution(nscenarios::Integer, max_active::Int, κ::T, ::Type{F}, ::Type{T}, ::Type{A}, ::Type{S}) where {F <: AbstractFeasibility, T <: AbstractFloat, A <: AbstractVector, S <: LQSolver}
        return new{F,T,A,S}(AsynchronousData{T}(),
                            Vector{A}(),
                            Vector{A}(),
                            Vector{Int}(),
                            Vector{SubWorker{F,T,A,S}}(undef, nworkers()),
                            Vector{Work}(undef,nworkers()),
                            Vector{Work}(undef,nworkers()),
                            Vector{MetaData}(undef,nworkers()),
                            RemoteChannel(() -> DecisionChannel(Dict{Int,A}())),
                            RemoteChannel(() -> Channel{QCut{T}}(max_active*nworkers()*nscenarios)),
                            Vector{Future}(undef,nworkers()),
                            Vector{Bool}(),
                            Vector{Bool}(),
                            max_active,
                            κ)
    end
end

function initialize_subproblems!(execution::AsynchronousExecution,
                                 scenarioproblems::AbstractScenarioProblems,
                                 x::AbstractVector,
                                 subsolver::MPB.AbstractMathProgSolver)
    load_subproblems!(execution.subworkers, scenarioproblems, x, subsolver)
    return nothing
end

function finish_initilization!(lshaped::AbstractLShapedSolver, execution::AsynchronousExecution)
    # Load initial decision
    put!(execution.decisions, 1, copy(lshaped.x))
    # Prepare memory
    push!(execution.subobjectives, fill(1e10, nthetas(lshaped)))
    push!(execution.model_objectives, fill(-1e10, nthetas(lshaped)))
    push!(execution.finished, 0)
    push!(execution.triggered, false)
    push!(execution.added, false)
    push!(lshaped.Q_history, Inf)
    push!(lshaped.θ_history, -Inf)
    log_regularization!(lshaped)

    # Prepare work channels
    for w in workers()
        execution.work[w-1] = RemoteChannel(() -> Channel{Int}(execution.max_active+1), w)
        execution.finalize[w-1] = RemoteChannel(() -> Channel{Int}(1), w)
        execution.metadata[w-1] = RemoteChannel(() -> MetaChannel(), w)
        put!(execution.work[w-1], 1)
        put!(execution.metadata[w-1], 1, :gap, Inf)
    end
    return nothing
end

function start_workers!(lshaped::AbstractLShapedSolver, execution::AsynchronousExecution)
    for w in workers()
        worker_aggregator = remote_aggregator(lshaped.aggregation, scenarioproblems(lshaped.stochasticprogram), w)
        execution.active_workers[w-1] = remotecall(work_on_subproblems!,
                                                   w,
                                                   execution.subworkers[w-1],
                                                   execution.work[w-1],
                                                   execution.finalize[w-1],
                                                   execution.cutqueue,
                                                   execution.decisions,
                                                   execution.metadata[w-1],
                                                   worker_aggregator)
    end
    return nothing
end

function close_workers!(::AbstractLShapedSolver, execution::AsynchronousExecution)
    t = execution.data.timestamp
    map((w,aw)->!isready(aw) && put!(w,t), execution.finalize, execution.active_workers)
    map((w,aw)->!isready(aw) && put!(w,-1), execution.work, execution.active_workers)
    map(wait, execution.active_workers)
    return nothing
end

function resolve_subproblems!(::AbstractLShapedSolver, ::AsynchronousExecution)
    return nothing
end

function calculate_objective_value(lshaped::AbstractLShapedSolver, execution::AsynchronousExecution)
    return get_obj(lshaped)⋅decision(lshaped) + eval_second_stage(execution.subworkers, decision(lshaped))
end

function current_decision(lshaped::AbstractLShapedSolver, execution::AsynchronousExecution)
    t = timestamp(lshaped)
    return fetch(execution.decisions, t)
end

function timestamp(::AbstractLShapedSolver, execution::AsynchronousExecution)
    return execution.data.timestamp
end

function incumbent_decision(::AbstractLShapedSolver, t::Integer, regularization::AbstractRegularization, execution::AsynchronousExecution)
    return t > 1 ? fetch(execution.decisions, regularization.incumbents[t]) : regularization.ξ
end

function incumbent_objective(::AbstractLShapedSolver, t::Integer, regularization::AbstractRegularization, ::AsynchronousExecution)
    return t > 1 ? regularization.Q̃_history[regularization.incumbents[t]] : regularization.data.Q̃
end

function incumbent_trustregion(::AbstractLShapedSolver, t::Integer, rd::RegularizedDecomposition, ::AsynchronousExecution)
    return rd.σ_history[t]
end

function incumbent_trustregion(::AbstractLShapedSolver, t::Integer, tr::TrustRegion, ::AsynchronousExecution)
    return tr.Δ_history[t]
end

function readd_cuts!(lshaped::AbstractLShapedSolver, consolidation::Consolidation, execution::AsynchronousExecution)
    for i in eachindex(consolidation.cuts)
        for cut in consolidation.cuts[i]
            add_cut!(lshaped, cut, execution.θs[i], execution.subobjectives[i], sum(execution.subobjectives[i]), fetch(execution.decisions, i), check = false)
        end
        for cut in consolidation.feasibility_cuts[i]
            add_cut!(lshaped, cut, execution.θs[i], execution.subobjectives[i], Inf)
        end
    end
    return nothing
end

function subobjectives(lshaped::AbstractLShapedSolver, execution::AsynchronousExecution)
    t = timestamp(lshaped)
    return execution.subobjectives[t]
end

function set_subobjectives(lshaped::AbstractLShapedSolver, Qs::AbstractVector, execution::AsynchronousExecution)
    t = timestamp(lshaped)
    execution.subobjectives[t] .= Qs
    return nothing
end

function model_objectives(lshaped::AbstractLShapedSolver, execution::AsynchronousExecution)
    t = timestamp(lshaped)
    θs = t > 1 ? execution.model_objectives[t-1] : fill(-1e10, nthetas(lshaped))
    return θs
end

function set_model_objectives(lshaped::AbstractLShapedSolver, θs::AbstractVector, execution::AsynchronousExecution)
    t = timestamp(lshaped)
    ids = active_model_objectives(lshaped)
    execution.model_objectives[t][ids] .= θs[ids]
    return nothing
end

function fill_submodels!(lshaped::AbstractLShapedSolver, scenarioproblems, execution::AsynchronousExecution)
    return fill_submodels!(execution.subworkers, decision(lshaped), scenarioproblems)
end

function solve_master!(lshaped::AbstractLShapedSolver, execution::AsynchronousExecution)
    try
        solve_regularized_master!(lshaped, lshaped.regularization)
    catch
        # Master problem could not be solved for some reason.
        @unpack Q,θ = lshaped.data
        gap = abs(θ-Q)/(abs(Q)+1e-10)
        # Always print this warning
        @warn "Master problem could not be solved, solver returned status $(status(lshaped.mastersolver)). The following relative tolerance was reached: $(@sprintf("%.1e",gap)). Aborting procedure."
        return :StoppedPrematurely
    end
    if status(lshaped.mastersolver) == :Infeasible
        # Asynchronicity can sometimes yield an infeasible master. If so, try to continue
        if execution.max_active > 1 || execution.κ < 1.0
            # Ensure that there are still cuts to consider
            if isready(execution.cutqueue)
                return :Valid
            end
        end
        # Otherwise, abort
        @warn "Master is infeasible. Aborting procedure."
        return :Infeasible
    end
    if status(lshaped.mastersolver) == :Unbounded
        @warn "Master is unbounded. Aborting procedure."
        return :Unbounded
    end
    return :Optimal
end

function iterate!(lshaped::AbstractLShapedSolver, execution::AsynchronousExecution{F,T}) where {F <: AbstractFeasibility, T <: AbstractFloat}
    wait(execution.cutqueue)
    while isready(execution.cutqueue)
        newiterate = false
        # Add new cuts from subworkers
        t::Int, cut::SparseHyperPlane{T} = take!(execution.cutqueue)
        execution.data.timestamp = t
        # Break if any subproblem is infeasible or unbounded
        if infeasible(cut)
            @warn "Stochastic program is not second-stage feasible at the current decision. Rerun procedure with feasibility_cuts = true to use feasibility cuts."
            return :Infeasible
        end
        if !bounded(cut)
            return :Unbounded
        end
        # Otherwise, add new cut to master and update bookkeeping
        lshaped.execution.added[t] |= add_cut!(lshaped, cut)
        lshaped.execution.finished[t] += nsubproblems(cut)
        # Asynchronicity parameter should be 1 first iteration
        κ = t > 1 ? execution.κ : 1.0
        if execution.finished[t] == nthetas(lshaped)
            # All work from iteration t complete
            lshaped.Q_history[t] = current_objective_value(lshaped)
            lshaped.data.Q = lshaped.Q_history[t]
            lshaped.data.θ = t > 1 ? lshaped.θ_history[t-1] : -1e10
            # Update incumbent (if applicable)
            # lshaped.x .= fetch(execution.decisions, t)
            take_step!(lshaped)
            # Optimal if not using regularization and no cuts were added
            if lshaped.regularization isa NoRegularization && !execution.added[t]
                # Optimal, final log
                log!(lshaped, t; optimal = true)
                return :Optimal
            end
            # Consolidate (if applicable)
            consolidate!(lshaped, lshaped.consolidation)
            # Decrease number of active iterations
            execution.data.active = max(0, execution.data.active-1)
            newiterate |= !(lshaped.regularization isa NoRegularization && execution.triggered[t])
            newiterate |= t == 1
        elseif execution.finished[t] >= κ*nthetas(lshaped) && !execution.triggered[t] && execution.data.active < execution.max_active && execution.added[t]
            execution.triggered[t] = true
            newiterate |= true
        else
            newiterate |= false
        end
        # Generate new candidate decision
        if newiterate
            t = lshaped.data.iterations
            execution.data.timestamp = t
            # Resolve master
            status = solve_master!(lshaped)
            if status != :Optimal
                return status
            end
            # Update master solution
            update_solution!(lshaped, lshaped.mastersolver)
            θ = get_obj(lshaped)⋅lshaped.x + sum(execution.model_objectives[t])
            lshaped.data.θ = θ
            lshaped.θ_history[t] = θ
            # Check if optimal
            if check_optimality(lshaped)
                # Optimal, final log
                log!(lshaped, t)
                return :Optimal
            end
            # Project (if applicable)
            project!(lshaped)
            # Log progress at current timestamp
            log_regularization!(lshaped, t)
            # Update workers
            put!(execution.decisions, t+1, copy(lshaped.x))
            for w in workers()
                if !isready(execution.active_workers[w-1])
                    put!(execution.work[w-1], t+1)
                    put!(execution.metadata[w-1], t+1, :gap, gap(lshaped))
                end
            end
            # New active iteration
            execution.data.active += 1
            # Prepare memory for next iteration
            push!(execution.subobjectives, fill(1e10, nthetas(lshaped)))
            push!(execution.model_objectives, fill(-1e10, nthetas(lshaped)))
            push!(execution.finished, 0)
            push!(execution.triggered, false)
            push!(execution.added, false)
            # Log progress
            log!(lshaped)
            lshaped.θ_history[t+1] = -Inf
        end
    end
    # Just return a valid status for this iteration
    return :Valid
end

# API
# ------------------------------------------------------------
function (execution::Asynchronous)(nscenarios::Integer, ::Type{F}, ::Type{T}, ::Type{A}, ::Type{S}) where {F <: AbstractFeasibility, T <: AbstractFloat, A <: AbstractVector, S <: LQSolver}
    return AsynchronousExecution(nscenarios, execution.max_active, execution.κ, F, T, A, S)
end

function str(::Asynchronous)
    return "Asynchronous "
end