text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

isdefined(Main, :OffsetArrays) || @eval Main include("testhelpers/OffsetArrays.jl")
using .Main.OffsetArrays
using DelimitedFiles
using Random
using LinearAlgebra
using Statistics

const OAs_name = join(fullname(OffsetArrays), ".")

let
# Basics
v0 = rand(4)
v = OffsetArray(v0, (-3,))
h = OffsetArray([-1,1,-2,2,0], (-3,))
@test axes(v) == (-2:1,)
@test size(v) == (4,)
@test size(v, 1) == 4

A0 = [1 3; 2 4]
A = OffsetArray(A0, (-1,2))                   # IndexLinear
S = OffsetArray(view(A0, 1:2, 1:2), (-1,2))   # IndexCartesian
@test axes(A) == axes(S) == (0:1, 3:4)
@test size(A) == (2,2)
@test size(A, 1) == 2

# Scalar indexing
@test A[0,3] == A[1] == A[0,3,1] == S[0,3] == S[1] == S[0,3,1] == 1
@test A[1,3] == A[2] == A[1,3,1] == S[1,3] == S[2] == S[1,3,1] == 2
@test A[0,4] == A[3] == A[0,4,1] == S[0,4] == S[3] == S[0,4,1] == 3
@test A[1,4] == A[4] == A[1,4,1] == S[1,4] == S[4] == S[1,4,1] == 4
@test_throws BoundsError A[1,1]
@test_throws BoundsError S[1,1]
@test_throws BoundsError A[0,3,2]
@test_throws BoundsError S[0,3,2]
# partial indexing
S3 = OffsetArray(view(reshape(Vector(1:4*3*1), 4, 3, 1), 1:3, 1:2, :), (-1,-2,1))
@test S3[1,-1] == 2
@test S3[1,0] == 6
@test_throws BoundsError S3[1,1]
@test_throws BoundsError S3[1,-2]
S4 = OffsetArray(view(reshape(Vector(1:4*3*2), 4, 3, 2), 1:3, 1:2, :), (-1,-2,1))
@test S4[1,-1,2] == 2
@test S4[1,0,2] == 6
@test_throws BoundsError S4[1,1,2]
@test_throws BoundsError S4[1,-2,2]

# Vector indexing
@test A[:, 3] == S[:, 3] == OffsetArray([1,2], (A.offsets[1],))
@test A[:, 4] == S[:, 4] == OffsetArray([3,4], (A.offsets[1],))
@test_throws BoundsError A[:, 1]
@test_throws BoundsError S[:, 1]
@test A[0, :] == S[0, :] == OffsetArray([1,3], (A.offsets[2],))
@test A[1, :] == S[1, :] == OffsetArray([2,4], (A.offsets[2],))
@test_throws BoundsError A[2, :]
@test_throws BoundsError S[2, :]
@test A[0:1, 3] == S[0:1, 3] == [1,2]
@test A[[1,0], 3] == S[[1,0], 3] == [2,1]
@test A[0, 3:4] == S[0, 3:4] == [1,3]
@test A[1, [4,3]] == S[1, [4,3]] == [4,2]
@test A[:, :] == S[:, :] == A

A_3_3 = OffsetArray(Matrix{Int}(undef, 3,3), (-2,-1))
A_3_3[:, :] = reshape(1:9, 3, 3)
for i = 1:9 @test A_3_3[i] == i end
A_3_3[-1:1, 0:2] = reshape(1:9, 3, 3)
for i = 1:9 @test A_3_3[i] == i end
A_3_3[:, :] = 1:9
for i = 1:9 @test A_3_3[i] == i end
A_3_3[-1:1, 0:2] = 1:9
for i = 1:9 @test A_3_3[i] == i end
A_3_3[:] = 1:9
for i = 1:9 @test A_3_3[i] == i end
A_3_3[1:9] = 1:9
for i = 1:9 @test A_3_3[i] == i end

# CartesianIndexing
@test A[CartesianIndex((0,3))] == S[CartesianIndex((0,3))] == 1
@test_throws BoundsError A[CartesianIndex(1,1)]
@test_throws BoundsError S[CartesianIndex(1,1)]
@test eachindex(A) == 1:4
@test eachindex(S) == CartesianIndices(axes(S)) == CartesianIndices(map(Base.IdentityUnitRange, (0:1,3:4)))

# LinearIndices
# issue 27986
let a1 = [11,12,13], a2 = [1 2; 3 4]
    b1 = OffsetArray(a1, (-3,))
    i1 = LinearIndices(b1)
    @test i1[-2] == -2
    @test_throws BoundsError i1[-3]
    @test_throws BoundsError i1[1]
    @test i1[-2:end] === -2:0
    @test @inferred(i1[-2:0]) === -2:0
    @test_throws BoundsError i1[-3:end]
    @test_throws BoundsError i1[-2:1]
    b2 = OffsetArray(a2, (-3,5))
    i2 = LinearIndices(b2)
    @test i2[3] == 3
    @test_throws BoundsError i2[0]
    @test_throws BoundsError i2[5]
    @test @inferred(i2[2:3])   === 2:3
    @test @inferred(i2[1:2:4]) === 1:2:3
    @test_throws BoundsError i2[1:5]
    @test_throws BoundsError i2[1:2:5]
end

# logical indexing
@test A[A .> 2] == [3,4]
@test_throws BoundsError h[trues(2)]
@test_throws BoundsError h[trues(5)]
@test h[OffsetArray(trues(5), (-3,))] == parent(h)
@test h[OffsetArray([true,false,false,true,true], (-3,))] == parent(h)[[1,4,5]]
@test A[OffsetArray([true false; false true], A.offsets)] == [1,4]
@test A[OffsetArray([true true;  false true], A.offsets)] == [1,3,4]
@test_throws BoundsError A[[true true;  false true]]

# begin, end
a0 = rand(2,3,4,2)
a = OffsetArray(a0, (-2,-3,4,5))
@test a[begin,end,end,begin] == a0[begin,end,end,begin] ==
      a0[1,3,4,1] == a0[end-1,begin+2,begin+3,end-1]

# view
S = view(A, :, 3)
@test S == OffsetArray([1,2], (A.offsets[1],))
@test S[0] == 1
@test S[1] == 2
@test_throws BoundsError S[2]
@test axes(S) === (Base.IdentityUnitRange(0:1),)
S = view(A, 0, :)
@test S == OffsetArray([1,3], (A.offsets[2],))
@test S[3] == 1
@test S[4] == 3
@test_throws BoundsError S[1]
@test axes(S) === (Base.IdentityUnitRange(3:4),)
S = view(A, 0:0, 4)
@test S == [3]
@test S[1] == 3
@test_throws BoundsError S[0]
@test axes(S) === (Base.OneTo(1),)
S = view(A, 1, 3:4)
@test S == [2,4]
@test S[1] == 2
@test S[2] == 4
@test_throws BoundsError S[3]
@test axes(S) === (Base.OneTo(2),)
S = view(A, :, :)
@test S == A
@test S[0,3] == S[1] == 1
@test S[1,3] == S[2] == 2
@test S[0,4] == S[3] == 3
@test S[1,4] == S[4] == 4
@test_throws BoundsError S[1,1]
@test axes(S) === Base.IdentityUnitRange.((0:1, 3:4))
# https://github.com/JuliaArrays/OffsetArrays.jl/issues/27
g = OffsetArray(Vector(-2:3), (-3,))
gv = view(g, -1:2)
@test axes(gv, 1) === Base.OneTo(4)
@test collect(gv) == -1:2
gv = view(g, OffsetArray(-1:2, (-2,)))
@test axes(gv, 1) === Base.IdentityUnitRange(-1:2)
@test collect(gv) == -1:2
gv = view(g, OffsetArray(-1:2, (-1,)))
@test axes(gv, 1) === Base.IdentityUnitRange(0:3)
@test collect(gv) == -1:2

# iteration
for (a,d) in zip(A, A0)
    @test a == d
end

# show
io = IOBuffer()
show(io, v)
str = String(take!(io))
show(io, v0)
@test str == String(take!(io))
show(io, A)
str = String(take!(io))
@test str == "[1 3; 2 4]"
show(io, S)
str = String(take!(io))
@test str == "[1 3; 2 4]"
show(io, MIME("text/plain"), A)
strs = split(strip(String(take!(io))), '\n')
@test strs[2] == " 1  3"
@test strs[3] == " 2  4"
v = OffsetArray(rand(3), (-2,))
show(io, v)
str = String(take!(io))
show(io, parent(v))
@test str == String(take!(io))
smry = summary(v)
@test occursin("OffsetArray{Float64,1", smry)
@test occursin("with indices -1:1", smry)
function cmp_showf(printfunc, io, A; options = ())
    ioc = IOContext(io, :limit => true, :compact => true, options...)
    printfunc(ioc, A)
    str1 = String(take!(io))
    printfunc(ioc, parent(A))
    str2 = String(take!(io))
    @test str1 == str2
end
cmp_showf(Base.print_matrix, io, OffsetArray(rand(5,5), (10,-9)))       # rows&cols fit
cmp_showf(Base.print_matrix, io, OffsetArray(rand(10^3,5), (10,-9)))    # columns fit
cmp_showf(Base.print_matrix, io, OffsetArray(rand(5,10^3), (10,-9)))    # rows fit
cmp_showf(Base.print_matrix, io, OffsetArray(rand(10^3,10^3), (10,-9))) # neither fits
cmp_showf(Base.print_matrix, io, OffsetArray(reshape(range(-0.212121212121, stop=2/11, length=3*29), 3, 29), (-2, -15)); options=(:displaysize=>(53,210),))
cmp_showf(show, io, OffsetArray(collect(1:100), (100,)))   # issue #31641

targets1 = ["0-dimensional $OAs_name.OffsetArray{Float64,0,Array{Float64,0}}:\n1.0",
            "1-element $OAs_name.OffsetArray{Float64,1,Array{Float64,1}} with indices 2:2:\n 1.0",
            "1×1 $OAs_name.OffsetArray{Float64,2,Array{Float64,2}} with indices 2:2×3:3:\n 1.0",
            "1×1×1 $OAs_name.OffsetArray{Float64,3,Array{Float64,3}} with indices 2:2×3:3×4:4:\n[:, :, 4] =\n 1.0",
            "1×1×1×1 $OAs_name.OffsetArray{Float64,4,Array{Float64,4}} with indices 2:2×3:3×4:4×5:5:\n[:, :, 4, 5] =\n 1.0"]
targets2 = ["(fill(1.0), fill(1.0))",
            "([1.0], [1.0])",
            "([1.0], [1.0])",
            "([1.0], [1.0])",
            "([1.0], [1.0])"]
@testset "printing of OffsetArray with n=$n" for n = 0:4
    a = OffsetArray(fill(1.,ntuple(d->1,n)), ntuple(identity,n))
    show(IOContext(io, :limit => true), MIME("text/plain"), a)
    @test String(take!(io)) == targets1[n+1]
    show(IOContext(io, :limit => true), MIME("text/plain"), (a,a))
    @test String(take!(io)) == targets2[n+1]
end
P = OffsetArray(rand(8,8), (1,1))
PV = view(P, 2:3, :)
@test endswith(summary(PV), "with indices Base.OneTo(2)×2:9")

# Similar
B = similar(A, Float32)
@test isa(B, OffsetArray{Float32,2})
@test axes(B) === axes(A)
B = similar(A, (3,4))
@test isa(B, Array{Int,2})
@test size(B) == (3,4)
@test axes(B) === (Base.OneTo(3), Base.OneTo(4))
B = similar(A, (-3:3,1:4))
@test isa(B, OffsetArray{Int,2})
@test axes(B) === Base.IdentityUnitRange.((-3:3, 1:4))
B = similar(parent(A), (-3:3,1:4))
@test isa(B, OffsetArray{Int,2})
@test axes(B) === Base.IdentityUnitRange.((-3:3, 1:4))

# Indexing with OffsetArray indices
i1 = OffsetArray([2,1], (-5,))
i1 = OffsetArray([2,1], -5)
b = A0[i1, 1]
@test axes(b) === (Base.IdentityUnitRange(-4:-3),)
@test b[-4] == 2
@test b[-3] == 1
b = A0[1,i1]
@test axes(b) === (Base.IdentityUnitRange(-4:-3),)
@test b[-4] == 3
@test b[-3] == 1
v = view(A0, i1, 1)
@test axes(v) === (Base.IdentityUnitRange(-4:-3),)
v = view(A0, 1:1, i1)
@test axes(v) === (Base.OneTo(1), Base.IdentityUnitRange(-4:-3))

# copyto! and fill!
a = OffsetArray{Int}(undef, (-3:-1,))
fill!(a, -1)
copyto!(a, (1,2))   # non-array iterables
@test a[-3] == 1
@test a[-2] == 2
@test a[-1] == -1
fill!(a, -1)
copyto!(a, -2, (1,2))
@test a[-3] == -1
@test a[-2] == 1
@test a[-1] == 2
@test_throws BoundsError copyto!(a, 1, (1,2))
fill!(a, -1)
copyto!(a, -2, (1,2,3), 2)
@test a[-3] == -1
@test a[-2] == 2
@test a[-1] == 3
@test_throws BoundsError copyto!(a, -2, (1,2,3), 1)
fill!(a, -1)
copyto!(a, -2, (1,2,3), 1, 2)
@test a[-3] == -1
@test a[-2] == 1
@test a[-1] == 2

b = 1:2    # copy between AbstractArrays
bo = OffsetArray(1:2, (-3,))
@test_throws BoundsError copyto!(a, b)
fill!(a, -1)
copyto!(a, bo)
@test a[-3] == -1
@test a[-2] == 1
@test a[-1] == 2
fill!(a, -1)
copyto!(a, -2, bo)
@test a[-3] == -1
@test a[-2] == 1
@test a[-1] == 2
@test_throws BoundsError copyto!(a, -4, bo)
@test_throws BoundsError copyto!(a, -1, bo)
fill!(a, -1)
copyto!(a, -3, b, 2)
@test a[-3] == 2
@test a[-2] == a[-1] == -1
@test_throws BoundsError copyto!(a, -3, b, 1, 4)
am = OffsetArray{Int}(undef, (1:1, 7:9))  # for testing linear indexing
fill!(am, -1)
copyto!(am, b)
@test am[1] == 1
@test am[2] == 2
@test am[3] == -1
@test am[1,7] == 1
@test am[1,8] == 2
@test am[1,9] == -1

# map
dest = similar(am)
map!(+, dest, am, am)
@test dest[1,7] == 2
@test dest[1,8] == 4
@test dest[1,9] == -2

am = map(identity, a)
@test isa(am, OffsetArray)
@test am == a

# https://github.com/JuliaArrays/OffsetArrays.jl/issues/106
@test isequal(map(!, OffsetArray([true,missing],2)), OffsetArray([false, missing], 2))
@test isequal(map(!, OffsetArray([true missing; false true], 2, -1)), OffsetArray([false missing; true false], 2, -1))
P = view([true missing; false true; true false], 1:2:3, :)
@test IndexStyle(P) === IndexCartesian()
@test isequal(map(!, OffsetArray(P, 2, -1)), OffsetArray(map(!, P), 2, -1))

# dropdims
a0 = rand(1,1,8,8,1)
a = OffsetArray(a0, (-1,2,3,4,5))
@test @inferred(dropdims(a, dims=1)) == @inferred(dropdims(a, dims=(1,))) == OffsetArray(reshape(a, (1,8,8,1)), (2,3,4,5))
@test @inferred(dropdims(a, dims=5)) == @inferred(dropdims(a, dims=(5,))) == OffsetArray(reshape(a, (1,1,8,8)), (-1,2,3,4))
@test @inferred(dropdims(a, dims=(1,5))) == dropdims(a, dims=(5,1)) == OffsetArray(reshape(a, (1,8,8)), (2,3,4))
@test @inferred(dropdims(a, dims=(1,2,5))) == dropdims(a, dims=(5,2,1)) == OffsetArray(reshape(a, (8,8)), (3,4))
@test_throws ArgumentError dropdims(a, dims=0)
@test_throws ArgumentError dropdims(a, dims=(1,1))
@test_throws ArgumentError dropdims(a, dims=(1,2,1))
@test_throws ArgumentError dropdims(a, dims=(1,1,2))
@test_throws ArgumentError dropdims(a, dims=3)
@test_throws ArgumentError dropdims(a, dims=4)
@test_throws ArgumentError dropdims(a, dims=6)

# push!
v = OffsetArray(rand(4), (-3,))
v2 = copy(v)
@test push!(v2, 1) === v2
@test v2[axes(v, 1)] == v
@test v2[end] == 1
@test v2[begin] == v[begin] == v[-2]
v2 = copy(v)
@test push!(v2, 2, 1) === v2
@test v2[axes(v, 1)] == v
@test v2[end-1] == 2
@test v2[end] == 1

# append! from array
v2 = copy(v)
@test append!(v2, [2, 1]) === v2
@test v2[axes(v, 1)] == v
@test v2[lastindex(v)+1:end] == [2, 1]
# append! from HasLength iterator
v2 = copy(v)
@test append!(v2, (v for v in [2, 1])) === v2
@test v2[axes(v, 1)] == v
@test v2[lastindex(v)+1:end] == [2, 1]
# append! from SizeUnknown iterator
v2 = copy(v)
@test append!(v2, (v for v in [2, 1] if true)) === v2
@test v2[axes(v, 1)] == v
@test v2[lastindex(v)+1:end] == [2, 1]

# other functions
v = OffsetArray(v0, (-3,))
@test lastindex(v) == 1
@test v ≈ v
@test axes(v') === (Base.OneTo(1),Base.IdentityUnitRange(-2:1))
@test parent(v) == collect(v)
rv = reverse(v)
@test axes(rv) == axes(v)
@test rv[1] == v[-2]
@test rv[0] == v[-1]
@test rv[-1] == v[0]
@test rv[-2] == v[1]
cv = copy(v)
@test reverse!(cv) == rv

A = OffsetArray(rand(4,4), (-3,5))
@test lastindex(A) == 16
@test lastindex(A, 1) == 1
@test lastindex(A, 2) == 9
@test A ≈ A
@test axes(A') === Base.IdentityUnitRange.((6:9, -2:1))
@test parent(copy(A')) == copy(parent(A)')
@test collect(A) == parent(A)
@test maximum(A) == maximum(parent(A))
@test minimum(A) == minimum(parent(A))
@test extrema(A) == extrema(parent(A))
@test maximum(A, dims=1) == OffsetArray(maximum(parent(A), dims=1), A.offsets)
@test maximum(A, dims=2) == OffsetArray(maximum(parent(A), dims=2), A.offsets)
@test maximum(A, dims=1:2) == OffsetArray(maximum(parent(A), dims=1:2), A.offsets)
C = similar(A)
cumsum!(C, A, dims=1)
@test parent(C) == cumsum(parent(A), dims=1)
@test parent(cumsum(A, dims=1)) == cumsum(parent(A), dims=1)
cumsum!(C, A, dims=2)
@test parent(C) == cumsum(parent(A), dims=2)
R = similar(A, (1:1, 6:9))
maximum!(R, A)
@test parent(R) == maximum(parent(A), dims=1)
R = similar(A, (-2:1, 1:1))
maximum!(R, A)
@test parent(R) == maximum(parent(A), dims=2)
amin, iamin = findmin(A)
pmin, ipmin = findmin(parent(A))
@test amin == pmin
@test A[iamin] == amin
@test amin == parent(A)[ipmin]
amax, iamax = findmax(A)
pmax, ipmax = findmax(parent(A))
@test amax == pmax
@test A[iamax] == amax
@test amax == parent(A)[ipmax]
z = OffsetArray([0 0; 2 0; 0 0; 0 0], (-3,-1))
I = findall(!iszero, z)
@test I == [CartesianIndex(-1, 0)]
@test findall(!iszero,h) == [-2:1;]
@test findall(x->x>0, h) == [-1,1]
@test findall(x->x<0, h) == [-2,0]
@test findall(x->x==0, h) == [2]
@test mean(A_3_3) == median(A_3_3) == 5
@test mean(x->2x, A_3_3) == 10
@test mean(A_3_3, dims=1) == median(A_3_3, dims=1) == OffsetArray([2 5 8], A_3_3.offsets)
@test mean(A_3_3, dims=2) == median(A_3_3, dims=2) == OffsetArray(reshape([4,5,6],(3,1)), A_3_3.offsets)
@test var(A_3_3) == 7.5
@test std(A_3_3, dims=1) == OffsetArray([1 1 1], A_3_3.offsets)
@test std(A_3_3, dims=2) == OffsetArray(reshape([3,3,3], (3,1)), A_3_3.offsets)
@test sum(OffsetArray(fill(1,3000), -1000)) == 3000

# https://github.com/JuliaArrays/OffsetArrays.jl/issues/92
A92 = OffsetArray(reshape(1:27, 3, 3, 3), -2, -2, -2)
B92 = view(A92, :, :, -1:0)
@test axes(B92) == (-1:1, -1:1, 1:2)
@test sum(B92, dims=(2,3)) == OffsetArray(reshape([51,57,63], Val(3)), -2, -2, 0)
B92 = view(A92, :, :, Base.IdentityUnitRange(-1:0))
@test sum(B92, dims=(2,3)) == OffsetArray(reshape([51,57,63], Val(3)), -2, -2, -2)

@test norm(v) ≈ norm(parent(v))
@test norm(A) ≈ norm(parent(A))
@test dot(v, v) ≈ dot(v0, v0)

# Prior to its removal from Base, cumsum_kbn was used here. To achieve the same level of
# accuracy in the tests, we need to use BigFloats with enlarged precision.
@testset "high-precision array reduction" begin
    setprecision(BigFloat, 500) do
        v  = OffsetArray(BigFloat[1,1e100,1,-1e100], (-3,)) .* 1000
        v2 = OffsetArray(BigFloat[1,-1e100,1,1e100], ( 5,)) .* 1000
        @test isa(v, OffsetArray)
        cv  = OffsetArray(BigFloat[1, 1e100, 1e100,2], (-3,)) .* 1000
        cv2 = OffsetArray(BigFloat[1,-1e100,-1e100,2], ( 5,)) .* 1000
        @test cumsum(v) ≈ cv
        @test cumsum(v2) ≈ cv2
        @test sum(v) ≈ sum(parent(v))
    end
end

io = IOBuffer()
writedlm(io, A)
seek(io, 0)
@test readdlm(io, eltype(A)) == parent(A)

amin, amax = extrema(parent(A))
@test clamp.(A, (amax+amin)/2, amax).parent == clamp.(parent(A), (amax+amin)/2, amax)

@test unique(A, dims=1) == OffsetArray(parent(A), 0, first(axes(A, 2)) - 1)
@test unique(A, dims=2) == OffsetArray(parent(A), first(axes(A, 1)) - 1, 0)
v = OffsetArray(rand(8), (-2,))
@test sort(v) == OffsetArray(sort(parent(v)), v.offsets)
@test sortslices(A, dims=1) == OffsetArray(sortslices(parent(A), dims=1), A.offsets)
@test sortslices(A, dims=2) == OffsetArray(sortslices(parent(A), dims=2), A.offsets)
@test sort(A, dims=1) == OffsetArray(sort(parent(A), dims=1), A.offsets)
@test sort(A, dims=2) == OffsetArray(sort(parent(A), dims=2), A.offsets)
# Issue #33977
soa = OffsetArray([2,2,3], -2)
@test searchsorted(soa, 1) == -1:-2
@test searchsortedfirst(soa, 1) == -1
@test searchsortedlast(soa, 1) == -2
@test first(sort!(soa; alg=QuickSort)) == 2
@test first(sort!(soa; alg=MergeSort)) == 2
soa = OffsetArray([2,2,3], typemax(Int)-4)
@test searchsorted(soa, 1) == typemax(Int)-3:typemax(Int)-4
@test searchsortedfirst(soa, 2) == typemax(Int) - 3
@test searchsortedlast(soa, 2) == typemax(Int) - 2
@test first(sort!(soa; alg=QuickSort)) == 2
@test first(sort!(soa; alg=MergeSort)) == 2

@test mapslices(sort, A, dims=1) == OffsetArray(mapslices(sort, parent(A), dims=1), A.offsets)
@test mapslices(sort, A, dims=2) == OffsetArray(mapslices(sort, parent(A), dims=2), A.offsets)

@test rotl90(A) == OffsetArray(rotl90(parent(A)), A.offsets[[2,1]])
@test rotr90(A) == OffsetArray(rotr90(parent(A)), A.offsets[[2,1]])
@test reverse(A, dims=1) == OffsetArray(reverse(parent(A), dims=1), A.offsets)
@test reverse(A, dims=2) == OffsetArray(reverse(parent(A), dims=2), A.offsets)

@test A .+ 1 == OffsetArray(parent(A) .+ 1, A.offsets)
@test 2*A == OffsetArray(2*parent(A), A.offsets)
@test A+A == OffsetArray(parent(A)+parent(A), A.offsets)
@test A.*A == OffsetArray(parent(A).*parent(A), A.offsets)

@test circshift(A, (-1,2)) == OffsetArray(circshift(parent(A), (-1,2)), A.offsets)

src = reshape(Vector(1:16), (4,4))
dest = OffsetArray(Matrix{Int}(undef, 4,4), (-1,1))
circcopy!(dest, src)
@test parent(dest) == [8 12 16 4; 5 9 13 1; 6 10 14 2; 7 11 15 3]
@test dest[1:3,2:4] == src[1:3,2:4]

# reshape
A = OffsetArray(rand(4,4), (-3,5))
@test vec(A) == reshape(A, :) == reshape(A, 16) == reshape(A, Val(1)) == A[:] == vec(A.parent)
A = OffsetArray(view(rand(4,4), 1:4, 4:-1:1), (-3,5))
@test vec(A) == reshape(A, :) == reshape(A, 16) == reshape(A, Val(1)) == A[:] == vec(A.parent)
# issue #33614
A = OffsetArray(-1:0, (-2,))
@test reshape(A, :) === A
Arsc = reshape(A, :, 1)
Arss = reshape(A, 2, 1)
@test Arsc[1,1] == Arss[1,1] == -1
@test Arsc[2,1] == Arss[2,1] == 0
@test_throws BoundsError Arsc[0,1]
@test_throws BoundsError Arss[0,1]
A = OffsetArray([-1,0], (-2,))
Arsc = reshape(A, :, 1)
Arsc[1,1] = 5
@test first(A) == 5

# broadcast
a = [1]
b = OffsetArray(a, (0,))
@test @inferred(a .+ b) == [2]
a = OffsetArray([1, -2, 1], (-2,))
@test a .* a' == OffsetArray([ 1 -2  1;
                              -2  4 -2;
                               1 -2  1], (-2,-2))

end # let

# Check that similar throws a MethodError rather than a
# StackOverflowError if no appropriate method has been defined
# (#18107)
module SimilarUR
    using Test
    struct MyURange <: AbstractUnitRange{Int}
        start::Int
        stop::Int
    end
    ur = MyURange(1,3)
    a = Vector{Int}(undef, 2)
    @test_throws MethodError similar(a, ur)
    @test_throws MethodError similar(a, Float64, ur)
    @test_throws MethodError similar(a, Float64, (ur,))
    @test_throws MethodError similar(a, (2.0,3.0))
end

@testset "Issue 28101" begin
    A = OffsetArray(reshape(16:-1:1, (4, 4)), (-3,5))
    @test maximum(A, dims=1) == OffsetArray(maximum(parent(A), dims=1), A.offsets)
end

@testset "in-place reductions with mismatched dimensionalities" begin
    B = OffsetArray(reshape(1:24, 4, 3, 2), -5, 6, -7)
    for R in (fill(0, -4:-1), fill(0, -4:-1, 7:7), fill(0, -4:-1, 7:7, -6:-6))
        @test @inferred(maximum!(R, B)) == reshape(maximum(B, dims=(2,3)), axes(R)) == reshape(21:24, axes(R))
        @test @allocated(maximum!(R, B)) <= 800
        @test @inferred(minimum!(R, B)) == reshape(minimum(B, dims=(2,3)), axes(R)) == reshape(1:4, axes(R))
        @test @allocated(minimum!(R, B)) <= 800
    end
    for R in (fill(0, -4:-4, 7:9), fill(0, -4:-4, 7:9, -6:-6))
        @test @inferred(maximum!(R, B)) == reshape(maximum(B, dims=(1,3)), axes(R)) == reshape(16:4:24, axes(R))
        @test @allocated(maximum!(R, B)) <= 800
        @test @inferred(minimum!(R, B)) == reshape(minimum(B, dims=(1,3)), axes(R)) == reshape(1:4:9, axes(R))
        @test @allocated(minimum!(R, B)) <= 800
    end
    @test_throws DimensionMismatch maximum!(fill(0, -4:-1, 7:7, -6:-6, 1:1), B)
    @test_throws DimensionMismatch minimum!(fill(0, -4:-1, 7:7, -6:-6, 1:1), B)
    @test_throws DimensionMismatch maximum!(fill(0, -4:-4, 7:9, -6:-6, 1:1), B)
    @test_throws DimensionMismatch minimum!(fill(0, -4:-4, 7:9, -6:-6, 1:1), B)
    @test_throws DimensionMismatch maximum!(fill(0, -4:-4, 7:7, -6:-5, 1:1), B)
    @test_throws DimensionMismatch minimum!(fill(0, -4:-4, 7:7, -6:-5, 1:1), B)
end

max_stars_repo_name: Lun-Liu/SC-by-Default-Julia

max_issues_repo_name: Lun-Liu/SC-by-Default-Julia

max_forks_repo_name: Lun-Liu/SC-by-Default-Julia


text:
function CMAESOptimizer(arg0::jint, arg1::jdouble, arg2::jboolean, arg3::jint, arg4::jint, arg5::RandomGenerator, arg6::jboolean, arg7::ConvergenceChecker)
    return CMAESOptimizer((jint, jdouble, jboolean, jint, jint, RandomGenerator, jboolean, ConvergenceChecker), arg0, arg1, arg2, arg3, arg4, arg5, arg6, arg7)
end

function get_statistics_d_history(obj::CMAESOptimizer)
    return jcall(obj, "getStatisticsDHistory", List, ())
end

function get_statistics_fitness_history(obj::CMAESOptimizer)
    return jcall(obj, "getStatisticsFitnessHistory", List, ())
end

function get_statistics_mean_history(obj::CMAESOptimizer)
    return jcall(obj, "getStatisticsMeanHistory", List, ())
end

function get_statistics_sigma_history(obj::CMAESOptimizer)
    return jcall(obj, "getStatisticsSigmaHistory", List, ())
end

function optimize(obj::CMAESOptimizer, arg0::Vector{OptimizationData})
    return jcall(obj, "optimize", PointValuePair, (Vector{OptimizationData},), arg0)
end

ext: jl


text:
module SpikeNetOpt

using PyCall
using OrderedCollections
using LinearAlgebra
using SpikeSynchrony
using SpikingNeuralNetworks
SNN = SpikingNeuralNetworks
SNN.@load_units

using Evolutionary#, Test, Random
import DataStructures
using JLD

using Reexport

using UnicodePlots
#using Evolutionary, Test, Random
#using Pkg

#import Pkg
#using SpikingNeuralNetworks
#SNN = SpikingNeuralNetworks
using SpikeSynchrony
using Statistics
using JLD
using Distributed
using SharedArrays
using Plots
using UnicodePlots
#using Evolutionary
using Distributions
using LightGraphs
using Metaheuristics

# Export Function names and struct names found in file utils
export get_vm
function get_vm end
function checkmodel end
include("utils.jl")
include("current_search.jl")
#include("sciunit.jl")

end


text:
# Recurrent inference machines for seismic imaging
# Same notation as Putzky and Welling paper (2017)
# Author: Philipp Witte, pwitte3@gatech.edu
# Date: April 2020

using JUDI.TimeModeling, JUDI4Flux, LinearAlgebra, InvertibleNetworks
using PyPlot, Random, JLD, Flux
import Flux.Optimise.update!

Random.seed!(11)

# Training data
D = load("/data/pwitte3/models/overthrust_images_test.jld")
ntest = length(D["m"])

# Crop images and models to 400 x 120 (must be evenly dividable by 4)
for j=1:ntest
    D["m"][j] = D["m"][j][1:400, 1:120]
    D["m0"][j] = D["m0"][j][1:400, 1:120]
    D["dm"][j] = D["dm"][j][1:400, 1:120]
end

# Use one sample to set up operators
m = D["m"][1]
m0 = D["m0"][1]
dm = vec(D["dm"][1])

# Set up model structure
n = size(m0)
d = (25., 25.)
o = (0., 0.)

# Setup info and model structure
nsrc = 1
model0 = Model(n, d, o, m0)

# Set up receiver geometry
nxrec = 301
xrec = range(400f0, stop=9600f0, length=nxrec)
yrec = 0f0
zrec = range(250f0, stop=250f0, length=nxrec)

# receiver sampling and recording time
time = 2000f0   # receiver recording time [ms]
dt = 4f0    # receiver sampling interval [ms]

# Set up receiver structure
recGeometry = Geometry(xrec, yrec, zrec; dt=dt, t=time, nsrc=nsrc)

# Set up source geometry (cell array with source locations for each shot)
num_simsource = 21
xsrc = range(500f0, stop=9500f0, length=num_simsource)
ysrc = range(0f0, stop=0f0, length=num_simsource)
zsrc = range(20f0, stop=20f0, length=num_simsource)

# Set up source structure
srcGeometry = Geometry(xsrc, ysrc, zsrc; dt=dt, t=time)

# setup wavelet
f0 = 0.015f0     # MHz
wavelet = zeros(Float32, srcGeometry.nt[1], num_simsource)
for j=1:num_simsource
    wavelet[:, j] = ricker_wavelet(time, dt, f0) * randn(1)[1]/sqrt(num_simsource)
end
q = judiVector(srcGeometry, wavelet)
q_sim = zeros(Float32, size(wavelet))

# Set up info structure for linear operators
ntComp = get_computational_nt(srcGeometry, recGeometry, model0)
info = Info(prod(n), nsrc, ntComp)

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

# Return data as julia array
opt = Options(return_array=true)

# Setup operators
F0 = judiModeling(info, model0; options=opt)
Pr = judiProjection(info, recGeometry)
Ps = judiProjection(info, srcGeometry)

# Modeling operator and Jacobian
F0 = Pr*F0*Ps'  # initial model
J = judiJacobian(F0, q)

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

# Dimensions
n_in = 32
n_hidden = 64
batchsize = 1
maxiter = 4
Ψ(η) = identity(η)

# Unrolled loop
L = NetworkLoop(n[1], n[2], n_in, n_hidden, batchsize, maxiter, Ψ)

# Get network parameters and overwrite w/ saved values
iter = 200
P_curr = get_params(L)
P_save = load(join(["../data/network_params_iteration_", string(iter), ".jld"]))["P"]
for j=1:length(P_curr)
    P_curr[j].data = P_save[j].data
end

####################################################################################################
# Evaluate trained network for new data

# Step 1: Draw new m0, q_sim and d
i = 99 #randperm(ntest)[1]
η = D["dm"][i]
m0 = D["m0"][i]
for k=1:num_simsource
    q_sim[:,k] = ricker_wavelet(time, dt, f0) * randn(Float32, 1)[1]/sqrt(1f0*num_simsource)
end

# Set up Jacobian
J.model.m = m0
J.source[1] = q_sim
d = J*vec(η)    # Observed data
d = reshape(d, recGeometry.nt[1], nxrec, 1, batchsize)

# Step 2: Compute predicted image
η0 = zeros(Float32, n[1], n[2], 1, 1)   # zero initial guess
s0 = zeros(Float32, n[1], n[2], n_in-1, 1)
η_ = L.forward(η0, s0, d, J)[1]

# Comparison to RTM
rtm = J'*vec(d)

# Plot
η_ = η_ / norm(η_, 2)
η = η / norm(η, 2)
rtm = rtm / norm(rtm, 2)

figure(figsize=(5, 7))
subplot(3,1,1); imshow(reshape(η, model0.n)', cmap="gray", vmin=-2e-2, vmax=2e-2); title("True image")
subplot(3,1,2); imshow(reshape(rtm, model0.n)', cmap="gray", vmin=-2e-2, vmax=2e-2); title("RTM")
subplot(3,1,3); imshow(reshape(η_, model0.n)', cmap="gray", vmin=-2e-2, vmax=2e-2); title("i-RIM")

max_stars_repo_name: slimgroup/LoopUnrolledSeismicImaging

max_issues_repo_name: slimgroup/LoopUnrolledSeismicImaging

max_forks_repo_name: slimgroup/LoopUnrolledSeismicImaging


text:
function test_shiftableLoad()

    T = 4

    l = DR.DER.ShiftableLoad(
        T=T,
        cycle_begin_min=2,
        cycle_begin_max=3,
        cycle_length=2,
        cycle_load=[1.0, 2.0]
    )

    # GLPK does not support -in-Interval constraints, so we bridge
    model = MOI.Bridges.full_bridge_optimizer(GLPK.Optimizer(), Float64)

    lmin = zeros(T)
    lmax = 10.0 .* ones(T)
    price = collect(1:T)  # Cycle should start as early as possible

    # Instantiate initial model
    h = DR.DER.House(
        0, T,
        lmin, lmax, price,
        [l]
    )

    model, var2idx, con2idx = DR.DER.build_model!(h, model)

    MOI.optimize!(model)

    @test MOI.get(model, MOI.ObjectiveValue()) ≈ 8.0

    pnet = [
        var2idx[(:house, h.index, :pnet, t)]
        for t in 1:T
    ]

    @test MOI.get(model, MOI.VariablePrimal(), pnet[1]) ≈ 0.0
    @test MOI.get(model, MOI.VariablePrimal(), pnet[2]) ≈ 1.0
    @test MOI.get(model, MOI.VariablePrimal(), pnet[3]) ≈ 2.0
    @test MOI.get(model, MOI.VariablePrimal(), pnet[4]) ≈ 0.0

    return nothing
end

@testset "Shiftable" begin
    test_shiftableLoad()
end


text:
module Tutorial

# 1 Getting started

using Modia
@usingModiaPlot

# Define model
SimpleModel = Model(
    T = 0.4,
    x = Var(init=0.2),
    equation = :[T * der(x) + x = 1],
)

# Transform to ODE form
simpleModel = @instantiateModel(SimpleModel)

# Simulate with a default integrator of DifferentialEquations
simulate!(simpleModel, stopTime = 1.2)

# Simulate with a specific integrator (Tsit5) and use a unit for stopTime
simulate!(simpleModel, Tsit5(), stopTime = 1.2u"s")

# Produce a line plot with GLMakie
plot(simpleModel, ("x", "der(x)"), figure=1)

# 2.1 - Equation oriented models

using Modia

LowPassFilter = Model(
    T = 0.2,
    u = input,
    y = output | Var(:x),
    x = Var(init=0),
    equation = :[T * der(x) + x = u],
)

# 2.2 Merging models

HighPassFilter = LowPassFilter | Model( y = :(-x + u) )

setLogMerge(true)

LowAndHighPassFilter = LowPassFilter | Model(
        y = nothing,
        low = output | Var(:x),
        high = output | Var(:(-x + u)),
    )

setLogMerge(false)

@showModel LowAndHighPassFilter

# 2.3 Functions and tables

TestLowAndHighPassFilter = LowAndHighPassFilter | Model(
        u = :(sin( (time+1u"s")*u"1/s/s" * time)*u"V"),
        x = Var(init=0.2u"V")
    )

using Interpolations

table = CubicSplineInterpolation(0:0.5:2.0, [0.0, 0.7, 2.0, 1.8, 1.2])
TestLowAndHighPassFilter2 = TestLowAndHighPassFilter | Map(u = :(table(time*u"1/s")*u"V"))

# 2.4 Hierarchical modeling

TwoFilters = (
    high = HighPassFilter,
    low = LowPassFilter,
)

BandPassFilter = (
    u = input,
    y = output,
    high = HighPassFilter | Map(T=0.5, x=Var(init=0.1u"V")),
    low = LowPassFilter | Map(x=Var(init=0.2u"V")),
    equations = :[
        high.u = u,
        low.u = high.y,
        y = low.y]
)

TestBandPassFilter = BandPassFilter | Map(
        u = :(sin( (time+1u"s")*u"1/s/s" * time)*u"V")
    )

#bandPassFilter = @instantiateModel(TestBandPassFilter, logStateSelection=true)
#simulate!(bandPassFilter, Tsit5(), stopTime = 50u"s")
#plot(bandPassFilter, ["u", "y"], figure=2)

end


text:
### Common types

# type parameters:
#   T: element type
#   N: the number of dimensions
#   M: contiguous rank
#
# Note: M is the maximum dimension such that
# slices up to this dimension are contiguous.
#
# For example, given a 3D contiguous array a,
# the contiguous rank of a is 3, that of
# view(a, :, :, 1:3) is 2, and that of
# view(a, :, 1:2, 1:3) is 1, etc.
#
# This rank plays a key role in efficient
# linear indexing and type-stable subview
# calculation.
#
abstract StridedArrayView{T,N,M} <: DenseArray{T,N}
abstract ArrayView{T,N,M} <: StridedArrayView{T,N,M}
abstract UnsafeArrayView{T,N,M} <: StridedArrayView{T,N,M}

# a type for indicating contiguous rank (statically)
type ContRank{M} end

## auxiliary union types to simplify method definition
## (for internal use only)

@compat typealias Subs Union{Real,Colon,Range}
@compat typealias SubsNC Union{Real,Range}
@compat typealias SubsRange Union{Colon,Range}

### Common methods

iscontiguous(a::DenseArray) = false
iscontiguous(a::Array) = true

offset(a::Array) = 0

contiguousrank{T,N,M}(a::StridedArrayView{T,N,M}) = M

contrank{T,N}(a::Array{T,N}) = ContRank{N}
contrank{T,N,M}(a::StridedArrayView{T,N,M}) = ContRank{M}

length(a::StridedArrayView) = a.len
size(a::StridedArrayView) = a.shp

getdim{N}(s::NTuple{N,Int}, d::Integer) = (1 <= d <= N ? s[d] : 1)
size{T,N}(a::StridedArrayView{T,N}, d::Integer) = getdim(size(a), d)

## pointer conversion

if VERSION < v"0.4.0-dev+3768"
    convert{T}(::Type{Ptr{T}}, a::StridedArrayView{T}) = pointer(a)
else
    unsafe_convert{T}(::Type{Ptr{T}}, a::StridedArrayView{T}) = pointer(a)
end

## Create similar array

similar{T}(a::StridedArrayView, ::Type{T}, dims::Dims) = Array(T, dims)


text:
include("algoTypes.jl")
include("Comparer.jl")
include("ArrayComparer.jl")
include("CalculateComparer.jl")
include("FastComparer.jl")
include("Entity.jl")
include("BestEntity.jl")
include("Population.jl")
include("PopulationManager.jl")
include("Epoch.jl")
include("GA.jl")
include("CrossEntropy.jl")
include("PSO.jl")

max_stars_repo_name: unbornchikken/julia-ml-proto

max_issues_repo_name: unbornchikken/julia-ml-proto

max_forks_repo_name: unbornchikken/julia-ml-proto


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# EAGO
# A development environment for robust and global optimization
# See https://github.com/PSORLab/EAGO.jl
#############################################################################
# src/eeago_optimizer/optimize/optimize_lp_cone.jl
# Contains the optimize! routines for LP, SOCP, (and in the future MILP and
# MISOCP) type problems. This also includes functions to add variables,
# linear constraints, soc constraints, and unpack solutions.
#############################################################################

#=
LP          -> COPY TO RELAXED SOLVER AND SOLVE
MILP        -> COPY TO RELAXED SOLVER AND SOLVE
SOCP        -> COPY TO RELAXED SOLVER AND SOLVE
MISOCP      -> COPY TO RELAXED SOLVER AND SOLVE
DIFF_CVX    -> COPY TO NLP SOLVER AND SOLVE (POTENTIAL MULTISTART)
NS_CVX      -> COPY TO NLP SOLVER AND SOLVE (POTENTIAL MULTISTART)
DIFF_NCVX   -> APPLY GLOBAL SOLVER (UNLESS USER REQUEST LOCAL SOLVE THEN NLP)
NS_NCVX     -> APPLY GLOBAL SOLVER (UNLESS USER REQUEST LOCAL SOLVE THEN NLP)
MINCVX      -> APPLY GLOBAL SOLVER (LOCAL SOLVE OPTION FUTURE FEATURE)
=#

function add_variables(m::Optimizer, optimizer::T, variable_number::Int) where T

    variable_index = fill(VI(1), variable_number)
    for i = 1:variable_number
        @inbounds variable_index[i] = MOI.add_variable(optimizer)
        relaxed_variable = SV(@inbounds variable_index[i])
        v_info = @inbounds m._working_problem._variable_info[i]
        if v_info.is_integer && v_info.is_fixed
            MOI.add_constraint(optimizer, relaxed_variable, ET(v_info.lower_bound))
        elseif v_info.is_integer
            MOI.add_constraint(optimizer, relaxed_variable, ZO())
        elseif v_info.is_fixed
            MOI.add_constraint(optimizer, relaxed_variable, ET(v_info.lower_bound))
        elseif v_info.has_lower_bound && v_info.has_upper_bound
            MOI.add_constraint(optimizer, relaxed_variable, IT(v_info.lower_bound, v_info.upper_bound))
        elseif v_info.has_lower_bound
            MOI.add_constraint(optimizer, relaxed_variable, GT(v_info.lower_bound))
        elseif v_info.has_upper_bound
            MOI.add_constraint(optimizer, relaxed_variable, LT(v_info.upper_bound))
        end
    end

    return variable_index
end

### LP and MILP routines
function add_linear_constraints!(m::Optimizer, opt::T) where T

    # add linear constraints
    for (func, set) in m._input_problem._linear_leq_constraints
         MOI.add_constraint(opt, func, set)
    end
    for (func, set) in m._input_problem._linear_geq_constraints
        MOI.add_constraint(opt, func, set)
    end
    for (func, set) in m._input_problem._linear_eq_constraints
        MOI.add_constraint(opt, func, set)
    end

    return nothing
end

### LP and MILP routines
function add_soc_constraints!(m::Optimizer, opt::T) where T

    for (func, set) in m._input_problem._conic_second_order
         MOI.add_constraint(opt, func, set)
    end

    return nothing
end

function add_sv_or_aff_obj!(m::Optimizer, opt::T) where T

    if m._input_problem._objective_type === SINGLE_VARIABLE
        MOI.set(opt, MOI.ObjectiveFunction{SV}(), m._input_problem._objective_sv)
    elseif m._input_problem._objective_type === SCALAR_AFFINE
        MOI.set(opt, MOI.ObjectiveFunction{SAF}(), m._input_problem._objective_saf)
    end

    return nothing
end

function unpack_local_solve!(m::Optimizer, opt::T) where T

    m._maximum_node_id = 0

    m._termination_status_code = MOI.get(opt, MOI.TerminationStatus())
    m._result_status_code = MOI.get(opt, MOI.PrimalStatus())
    m._feasible_solution_found = m._result_status_code === MOI.FEASIBLE_POINT

    if MOI.get(opt, MOI.ResultCount()) > 0

        objective_value = MOI.get(opt, MOI.ObjectiveValue())

        # corrects for standard printing multiplier
        if m._input_problem._optimization_sense === MOI.MAX_SENSE
            objective_value *= -1.0
        end

        m._global_lower_bound = objective_value
        m._global_upper_bound = objective_value
        m._objective_value = objective_value
        m._best_upper_value = objective_value
        m._solution_value = objective_value
    end

    m._continuous_solution = zeros(m._input_problem._variable_count)
    for i = 1:m._input_problem._variable_count
         m._continuous_solution[i] = MOI.get(opt, MOI.VariablePrimal(),  m._relaxed_variable_index[i])
    end

    m._run_time = time() - m._start_time

    return nothing
end

function optimize!(::Val{LP}, m::Optimizer)

    relaxed_optimizer = m.relaxed_optimizer
    MOI.empty!(relaxed_optimizer)

    m._relaxed_variable_index = add_variables(m, relaxed_optimizer, m._input_problem._variable_count)
    add_linear_constraints!(m, relaxed_optimizer)
    add_sv_or_aff_obj!(m, relaxed_optimizer)
    MOI.set(relaxed_optimizer, MOI.ObjectiveSense(), m._input_problem._optimization_sense)

    if m._parameters.verbosity < 5
        MOI.set(relaxed_optimizer, MOI.Silent(), true)
    end
    m._parse_time = time() - m._start_time

    MOI.optimize!(relaxed_optimizer)

    unpack_local_solve!(m, relaxed_optimizer)

    return nothing
end

optimize!(::Val{MILP}, m::Optimizer) = optimize!(Val{LP}(), m)

function optimize!(::Val{SOCP}, m::Optimizer)

    relaxed_optimizer = m.relaxed_optimizer
    MOI.empty!(relaxed_optimizer)

    m._relaxed_variable_index = add_variables(m, relaxed_optimizer, m._input_problem._variable_count)
    add_linear_constraints!(m, relaxed_optimizer)
    add_soc_constraints!(m, relaxed_optimizer)
    add_sv_or_aff_obj!(m, relaxed_optimizer)
    MOI.set(relaxed_optimizer, MOI.ObjectiveSense(), m._input_problem._optimization_sense)

    if m._parameters.verbosity < 5
        MOI.set(relaxed_optimizer, MOI.Silent(), true)
    end
    m._parse_time = time() - m._start_time

    MOI.optimize!(relaxed_optimizer)

    unpack_local_solve!(m, relaxed_optimizer)

    return nothing
end

optimize!(::Val{MISOCP}, m::Optimizer) = optimize!(Val{SOCP}(), m)


text:
abstract type CIService end

struct GitHubActions <: CIService
    username::String
    email::String
    api_hostname::String
    clone_hostname::String

    function GitHubActions(;
        username="github-actions[bot]",
        email="41898282+github-actions[bot]@users.noreply.github.com",
        api_hostname="https://api.github.com",
        clone_hostname="github.com",
    )
        return new(username, email, api_hostname, clone_hostname)
    end
end

struct GitLabCI <: CIService
    username::String
    email::String
    api_hostname::String
    clone_hostname::String

    function GitLabCI(;
        username="gitlab-ci[bot]",
        email="gitlab-ci[bot]@gitlab.com",
        api_hostname="https://gitlab.com/api/v4",
        clone_hostname="gitlab.com",
    )
        return new(username, email, api_hostname, clone_hostname)
    end
end

ci_repository(::GitHubActions, env::AbstractDict=ENV) = env["GITHUB_REPOSITORY"]
ci_repository(::GitLabCI, env::AbstractDict=ENV) = env["CI_PROJECT_PATH"]

ci_token(::GitHubActions, env::AbstractDict=ENV) = env["GITHUB_TOKEN"]
ci_token(::GitLabCI, env::AbstractDict=ENV) = env["GITLAB_TOKEN"]

function get_api_and_repo(ci::GitHubActions; env=ENV)
    token = GitHub.Token(ci_token(ci, env))
    api = GitHub.GitHubAPI(; token=token, url=ci.api_hostname)
    repo, _ = @mock GitForge.get_repo(api, ci_repository(ci, env))

    return api, repo
end

function get_api_and_repo(ci::GitLabCI; env=ENV)
    token = GitLab.PersonalAccessToken(ci_token(ci, env))
    api = GitLab.GitLabAPI(; token=token, url=ci.api_hostname)
    repo, _ = @mock GitForge.get_repo(api, ci_repository(ci, env))

    return api, repo
end

function get_api_and_repo(ci::Any)
    err = "Unknown CI Config: $(typeof(ci))"
    @error(err)
    return throw(ErrorException(err))
end

function auto_detect_ci_service(; env::AbstractDict=ENV)
    if haskey(env, "GITHUB_REPOSITORY")
        return GitHubActions()
    elseif haskey(env, "GITLAB_CI")
        return GitLabCI()
    else
        throw(UnableToDetectCIService("Could not detect CI service"))
    end
end


text:
### Economic dispatch without renewable energy sources
function EconomicDispatch(time,
						  nodes::Nodes,
			   		      powerplants::PowerPlants,
			   			  )

	G = powerplants.unit
	T = time

	# Initialise JuMP model: Economic dispatch
	ED_mod = JuMP.Model(with_optimizer(Gurobi.Optimizer))

	# Variables
	@variable(ED_mod, P[T, G] >= 0)

	# Constraints
	# Market clearing/power balance
	@constraint(ED_mod, MarketClearing[t = T],
	    sum(P[t, g] for g in G) == nodes.systemload[t]);

	# Upper generation limit
	@constraint(ED_mod, GenerationLimitUp[g = G, t = T],
	    P[t, g] <= powerplants.pmax[g]);

	# Lower generation limit
	@constraint(ED_mod, GenerationLimitDown[g = G, t = T],
	    P[t, g] >= powerplants.pmin[g]);

	# Ramp-up limit
	@constraint(ED_mod, RampUp[g = G, t = T; t > T[1]],
		P[t, g] - P[t-1, g] <= powerplants.rup[g])

	# Ramp-down limit
	@constraint(ED_mod, RampDown[g = G, t = T; t > T[1]],
		-(P[t, g] - P[t-1, g]) <= powerplants.rdn[g])

	# Objective function
	@objective(ED_mod, Min,
	    sum(powerplants.mc[g][1] * P[t, g] for g in G, t in T));

	# Initiate optimisation process
	JuMP.optimize!(ED_mod)

	# Export results
	P_opt = JuMP.value.(P)
	price = JuMP.dual.(MarketClearing)

	return(ED_mod,
	       P_opt,
		   price
	       )
end

### Economic dispatch with renewable energy sources
function EconomicDispatch(time,
						  nodes::Nodes,
			   		      powerplants::PowerPlants,
						  renewables::Renewables,
			   			  )

	T = time
	G = powerplants.unit
	R = renewables.unit

	# Initialise JuMP model: Economic dispatch
	ED_mod = JuMP.Model(with_optimizer(Gurobi.Optimizer))

	# Variables
	@variable(ED_mod, P[T, G] >= 0)
	@variable(ED_mod, P_R[T, R] >= 0)

	# Constraints
	# Market clearing/power balance
	@constraint(ED_mod, MarketClearing[t = T],
	    sum(P[t, g] for g in G) +
			sum(P_R[t, r] for r in R) == nodes.systemload[t]);

	# Upper generation limit
	@constraint(ED_mod, GenerationLimitUp[g = G, t = T],
	    P[t, g] <= powerplants.pmax[g]);

	# Lower generation limit
	@constraint(ED_mod, GenerationLimitDown[g = G, t = T],
	    P[t, g] >= powerplants.pmin[g]);

	# Upper generation limit for renewables
	@constraint(ED_mod, ResGenerationLimitUp[r = R, t = T],
	    P_R[t, r] <= renewables.generation[r][t]);

	# Ramp-up limit
	@constraint(ED_mod, RampUp[g = G, t = T; t > T[1]],
		P[t, g] - P[t-1, g] <= powerplants.rup[g])

	# Ramp-down limit
	@constraint(ED_mod, RampDown[g = G, t = T; t > T[1]],
		-(P[t, g] - P[t-1, g]) <= powerplants.rdn[g])

	# Objective function
	@objective(ED_mod, Min,
	    sum(powerplants.mc[g][1] * P[t, g] for g in G, t in T) +
			sum(renewables.mc[r] * P_R[t, r] for r in R, t in T));

	# Initiate optimisation process
	JuMP.optimize!(ED_mod)

	# Export results
	P_opt = JuMP.value.(P)
	P_R_opt = JuMP.value.(P_R)
	price = JuMP.dual.(MarketClearing)

	# Correct RES dispatch
	# Total summed RES generation from dispatch for each t
	ResTotalGeneration = zeros(length(T), 1)

	for t in T
		ResTotalGeneration[t-(T[1]-1)] = sum(P_R_opt.data[t-(T[1]-1), :])
	end

	# Maximum output possible
	ResMaxGeneration = zeros(length(T), 1)
	for t in T
		ResMaxGeneration[t-(T[1]-1)] = sum(renewables.generation[r][t] for r in R)
	end

	# Redistribution among RES
	ResUnitGeneration = zeros(length(T), length(R))
	for t in T
		for r in R
			ResUnitGeneration[t-(T[1]-1), r] = (renewables.generation[r][t] / ResMaxGeneration[t-(T[1]-1)]) *
			ResTotalGeneration[t-(T[1]-1)]
		end
	end

	# Distributed, corrected P_R_opt
	P_R_opt_dist = JuMP.Containers.DenseAxisArray(ResUnitGeneration, T, R)

	return(ED_mod,
	       P_opt,
		   P_R_opt,
		   P_R_opt_dist,
		   price
	       )
end


text:
durations = [1.5, 2.0, 2.5, 3.0, 3.5, 4.5, 5.0, 6.0, 7.5, 9.0, 10.5, 12.0, 12.5, 15.0]
min_periods = [0.5, 0.52, 0.6517, 0.7824, 0.912, 1.178, 1.3056, 1.567, 1.952, 2.343, 2.75, 3.14, 3.257, 3.91]
max_periods = [50.045, 118.626, 231.69, 400.359, 635.76, 725, 725, 725, 725, 725, 725, 725, 725, 725]
num_dur = 14

function get_legal_durations(period::Float64,duration::Float64)
  min_duration = 0.0
  max_duration = 0.0
  i = 1
  #determine what maximum and minimum durations were searched for this period
  while min_duration == 0.0 || max_duration == 0.0
    if i > 14
      println("No durations match this period")
      return(0.0,0.0)
    end
    if period <= max_periods[i] && min_duration == 0.0
      min_duration = durations[i]
    end
    if period >= min_periods[num_dur+1-i] && max_duration == 0.0
      max_duration = durations[num_dur+1-i]
    end
    i+=1
  end
  if duration<=max_duration/24 && duration>=min_duration/24
    return duration
  elseif duration>=max_duration/24
    return max_duration/24
  elseif duration <=min_duration/24
    return min_duration/24
  end
end


text:
function newton_raphson(f_,h_!,x; ρ = 0.5, c = 1e-4, tol = 1e-4, max_iter = 1000)
    grad= zeros(length(x))
    hes= zeros(length(x),length(x))
    y = 0.
    for i = 1:max_iter
        y = h_!(x, grad, hes)
        search_dir = -(cholfact(Positive, hes)\grad)
        (norm(search_dir) > tol) || break
        step_size = 1.
        while f_(x+search_dir*step_size) > y+c*step_size*(grad' * search_dir)
            step_size *= ρ
            step_size > 1e-10 || error("Linesearch failed! Problematic Hessian or Gradient?")
        end
        x = x + search_dir*step_size
    end
    return x, y, grad, hes
end


text:
"""
    distaz!(Ev::SeisEvent)

Compute Δ, Θ by the Haversine formula. Compute Δ, Θ by the Haversine formula.
Updates `Ev.data` with distance, azimuth, and backazimuth for each channel,
written to Ev.data.dist, Ev.data.az, and Ev.data.baz, respectively.
"""
function distaz!(S::SeisEvent)
  TD = getfield(S, :data)
  ChanLoc = getfield(TD, :loc)
  SrcLoc = getfield(getfield(S, :hdr), :loc)
  N = getfield(TD, :n)
  rec = Array{Float64, 2}(undef, N, 2)
  for i = 1:N
    loc = getindex(ChanLoc, i)
    if typeof(loc) == GeoLoc
      rec[i,:] = [loc.lat loc.lon]
    else
      error(string(":loc for channel ", i, " is not a GeoLoc!"))
    end
  end
  D = gcdist(SrcLoc.lat, SrcLoc.lon, rec)
  @assert size(D,1) == N
  TD.dist = D[:,1]
  TD.az   = D[:,2]
  TD.baz  = D[:,3]
  note!(TD, "distaz!")
  return nothing
end


text:
cd(@__DIR__); runapprox = true; allNGWPs = false; include("setups/sunflower.jl")
gr(dpi = 200)

##
VM_NGWP_Lsym = vm_ngwp(𝚽sym, GP_dual_Lsym)
f = matread("../datasets/sunflower_barbara_voronoi.mat")["f_eye_voronoi"]
G_Sig.f = reshape(f, (N, 1))
dmatrix_VM_Lsym = ngwp_analysis(G_Sig, VM_NGWP_Lsym)
dvec_vm_ngwp_Lsym, BS_vm_ngwp_Lsym = ngwp_bestbasis(dmatrix_VM_Lsym, GP_dual_Lsym)
important_idx = sortperm(dvec_vm_ngwp_Lsym[:].^2; rev = true)
println("================feye-VM-NGWP-Lsym-top-basis-vectors=================")
for i in 2:17
    dr, dc = BS_vm_ngwp_Lsym.levlist[important_idx[i]]
    w = VM_NGWP_Lsym[dr, dc, :]
    j, k, l = NGWP_jkl(GP_dual_Lsym, dr, dc)
    print("(j, k, l) = ($(j), $(k), $(l))  ")
    if j == jmax_Lsym
        print("φ_{$(GP_dual_Lsym.ind[dr]-1)}")
    end
    println()
    scatter_gplot(X; marker = w, ms = LinRange(4.0, 14.0, N), c = :greys)
    plt = plot!(xlim = [-1.2, 1.2], ylim = [-1.2, 1.2], frame = :none,
                cbar = false, clims = (-0.15, 0.15))
    savefig(plt,
        "../figs/SunFlower_feye_DAG_VM_NGWP_Lsym_ibv$(lpad(i,2,"0"))_j$(j)_k$(k)_l$(l).png")
end


text:
# Problems taken from
# A. Lovison - "Singular Continuation: Generating Piecewise
# Linear Approximations to Pareto Sets via Global Analysis"
# (Thomann refers to these as 'lovison1-4')
# lovison1 is basically Two Parabolas.
# Note: signs are switched because we minimize
# Note: there are more problems in the paper
# TODO: optimal and critical points!

abstract type Lovison <: MOP end;

for i = 1 : 4
    typename = Symbol("Lovison", i);
    typename_constrained = Symbol( "Lovison",i, "Constrained" );
    @eval begin
        @doc $("""
                struct $typename <: Lovison <: MOP end;

            Problem taken from
            A. Lovison - "Singular Continuation: Generating Piecewise
            Linear Approximations to Pareto Sets via Global Analysis"
        """)
        struct $typename <: Lovison end;

        @doc $("""
                struct $typename_constrained <: Lovison <: MOP end;

            Problem taken from
            A. Lovison - "Singular Continuation: Generating Piecewise
            Linear Approximations to Pareto Sets via Global Analysis"
        """)
        struct $typename_constrained <: Lovison end;

        num_vars(::Union{$typename, $typename_constrained}) = 2;
        num_objectives(::Union{$typename, $typename_constrained}) = 2;

        export $typename, $typename_constrained;
    end
end

function get_objectives(::Lovison1)
    return [
        (x :: Vector{R} where R<:Real) -> sum( [1.05; 0.95].* x.^2 );
        (x :: Vector{R} where R<:Real) -> .99*(x[1]-3)^2 + 1.03*(x[2]-2.5)^2
    ];
end

function get_objectives(::Lovison2)
    return [
        (x :: Vector{R} where R<:Real) -> x[2];
        (x :: Vector{R} where R<:Real) -> (x[2]- x[1]^3)/(x[1]+1)
    ];
end

function get_objectives(::Lovison3)
    return [
        (x :: Vector{R} where R<:Real) -> sum( x.^2 );
        (x :: Vector{R} where R<:Real) -> (x[1]-6)^2 - (x[2]+.3)^2
    ];
end

function get_objectives(::Lovison4)
    return [
        (x :: Vector{R} where R<:Real) -> sum( x.^2 ) + 4 * exp( -(x[1]+2)^2 - x[2]^2 ) -  exp( -(x[1]-2)^2 - x[2]^2 )
        (x :: Vector{R} where R<:Real) -> (x[1]-6)^2 + (x[2]+.5)^2
    ];
end

# box constraints as given by Thomann (p.XII)
constraints( ::Lovison1Constrained ) = Box(
    zeros(2),
    fill(3.0, 2)
);

constraints( ::Lovison2Constrained ) = Box(
    fill(-.5, 2),
    [0.0, 0.5]
);

constaints( ::Lovison3Constrained ) = Box(
    [0.0, -4.0],
    [6.0, 4.0]
);

constaints( ::Lovison4Constrained ) = Box(
    [0.0, -1.0],
    [6.0, 1.0]
);


text:
function converging!(
        objective::Function,
        ip::Vector{Int64},
        population::Matrix{Float64},
        n_population::Int64,
        n_gene::Int64)::Matrix{Float64}
    n_children::Int8 = 10
    children::Matrix{Float64} = zeros(n_children, n_gene + 1)
    for i in 1:n_children
        ip[3:end] = sample(
            collect(1:n_population), n_gene, replace=false
        )
        children[i, :] = crossover(population[ip, :], n_gene)
    end

    family::Matrix{Float64} = zeros(n_children + 2, n_gene + 1)
    @inbounds for i in 1:n_gene + 1
        @simd for j in 1:n_children
            family[j, i] = children[j, i]
        end
        family[n_children + 1, i] = population[ip[1], i]
        family[n_children + 2, i] = population[ip[2], i]
    end

    family = sortslices(family, dims=1, by=x -> x[end])

    ic2::Int8 = rand(2:n_children + 2)
    @inbounds for i in 1:n_gene + 1
        population[ip[1], i] = family[1, i]  # Best
        population[ip[2], i] = family[ic2, i]  # Random
    end

    if 1e12 <= population[ip[2], end]
        population[ip[2], end] = objective(population[ip[2], 1:n_gene])
    end

    population = sortslices(population, dims=1, by=x -> x[end])

    return population
end

function crossover(parents::Matrix{Float64}, n_gene::Int64)::Vector{Float64}
    # Extended Normal Distribution Xover
    child::Vector{Float64} = zeros(n_gene + 1)

    ALPHA::Float64 = sqrt((1.0 - 2.0 * 0.35^2.0)) / 2.0
    BETA::Float64 = 0.35 / sqrt(n_gene - 1)

    t1::Vector{Float64} = (parents[2, 1:n_gene] - parents[1, 1:n_gene]) ./ 2.0
    t2::Vector{Float64} = randn() * ALPHA * (parents[2, 1:n_gene] - parents[1, 1:n_gene])
    t3::Vector{Float64} = zeros(n_gene)
    centroid::Vector{Float64} = reshape(
        mean(parents[3:end, 1:n_gene], dims=1), n_gene
    )

    for i in 1:n_gene
        t3 += randn() * BETA * (parents[i + 2, 1:n_gene] - centroid)
    end

    @inbounds for i in 1:n_gene
        child[i] = t1[i] + t2[i] + t3[i]
        child[i] = clamp(child[i], 0.0, 1.0)
    end

    child[end] = 1e12

    return child
end


text:
module ExampleTool

using ..Ahorn, Maple

# !NOT PART OF THE EXAMPLE
begin
    const output = true
    const logMouseMotion = false
    const logSelection = false
    function log(text)
        if !output
            return
        elseif startswith(text, "mouseMotion")
            if logMouseMotion
                println(text)
            end
        elseif startswith(text, "selection")
            if logSelection
                println(text)
            end
        else
            println(text)
        end
    end
end

"""
The name displayed in the tool list
"""
displayName = "Example"

"""
Determines how the tools are grouped in the tool list

Ahorn Groups: Brushes, Placements
"""
group = "Example"

# Ahorn.ListContainer: https://github.com/CelestialCartographers/Ahorn/blob/master/src/helpers/list_view_helper.jl

"""
Called when this tool is selected.

Use this function to populate the `subTools`, `layers`, and `materials` lists, and to set up anything your tool needs while active.
If your tool has a draw function, it should be assigned here to Ahorn.redrawingFuncs["tools"]
"""
toolSelected(subTools::Ahorn.ListContainer, layers::Ahorn.ListContainer, materials::Ahorn.ListContainer) = log("toolSelected: subTools, layers, materials")
toolSelected() = log("toolSelected")

"""
Called when this tool is deselected (when another tool is selected) or when a new room is selected.

Use this function to cleanup any temporary resources your tool uses.
"""
cleanup() = log("cleanup")

# Called when a subtool of this tool is selected
subToolSelected(list::Ahorn.ListContainer, selected::String) = log("subToolSelected: list, selected")

# Called when a layer of this tool is selected
# Appears to be called twice?
layerSelected(list::Ahorn.ListContainer,  materials::Ahorn.ListContainer, selected::String) = log("layerSelected: list, materials, selected=$selected")
layerSelected(list::Ahorn.ListContainer, selected::String) = log("layerSelected: list, selected=$selected")

# Called when a material of this tool is selected
materialSelected(list::Ahorn.ListContainer, selected::String) = log("materialSelected: list, selected=$selected")

# Called when a material of this tool is double clicked
materialDoubleClicked(list::Ahorn.ListContainer,  selected::String) = log("materialDoubleClicked: list, selected=$selected")
materialDoubleClicked(selected::String) = log("materialDoubleClicked: selected=$selected")

# Called when a new room is selected
roomChanged(room::Room) = log("roomChanged: room")
roomChanged(map::Map, room::Room) = log("roomChanged: room, map")

# Called when this tool is selected, or when a new room is selected
# `layers` is the array of layers that can be drawn to in the current room (see Ahorn.DrawableRoom)
# global toolsLayer = Ahorn.getLayerByName(layers, "tools")
layersChanged(layers::Array{Ahorn.Layer, 1}) = log("layersChanged: layers")

# Called when the search bar is used
materialFiltered(list::Ahorn.ListContainer, text::String) = log("toolSelected: list, text=$text")
materialFiltered(list::Ahorn.ListContainer) = log("toolSelected: list")

# Called to determine the materials to be pushed to the top of the materials column
getFavorites() = log("getFavorites")

# Grid Based coordinates
selectionMotion(rect::Ahorn.Rectangle) = log("selectionMotion: rect=$rect")
selectionMotion(x1::Number, y1::Number, x2::Number, y2::Number) = log("selectionMotion: x1=$x1, y1=$y1, x2=$x2, y2=$y2")

# Absolute coordinates
selectionMotionAbs(rect::Ahorn.Rectangle) = log("selectionMotionAbs: rect=$rect")
selectionMotionAbs(x1::Number, y1::Number, x2::Number, y2::Number) = log("selectionMotionAbs: x1=$x1, y1=$y1, x2=$x2, y2=$y2")

# Grid Based coordinates
selectionFinish(rect::Ahorn.Rectangle) = log("selectionFinish: rect=$rect")
selectionFinish(x1::Number, y1::Number, x2::Number, y2::Number) = log("selectionFinishAbs: x1=$x1, y1=$y1, x2=$x2, y2=$y2")

# Absolute coordinates
selectionFinishAbs(rect::Ahorn.Rectangle) = log("selectionFinishAbs: rect=$rect")
selectionFinishAbs(x1::Number, y1::Number, x2::Number, y2::Number) = log("selectionFinishAbs: x1=$x1, y1=$y1, x2=$x2, y2=$y2")

# From Ahorn/event_handler.jl
const mouseHandlers = [
   :leftClick,
   :middleClick,
   :rightClick,
   :backClick,
   :forwardClick
]
# Ahorn.eventMouse = Union{Gtk.GdkEventButton, Gtk.GdkEventMotion, Gtk.GdkEventCrossing}
for handler in mouseHandlers
    @eval $handler() = log($(string(handler)))
    @eval $handler(event::Ahorn.eventMouse, camera::Ahorn.Camera) = log($(string(handler)) * ": event=$(typeof(event)), camera")
    # Grid Based coordinates
    @eval $handler(x::Number, y::Number) = log($(string(handler)) * ": x=$x, y=$y")
    # Absolute coordinates
    @eval $(Symbol(string(handler) * "Abs"))(x::Number, y::Number) = log($(string(handler) * "Abs") * ": x=$x, y=$y")
end

mouseMotion(event::Ahorn.eventMouse, camera::Ahorn.Camera) = log("mouseMotion: event=$(typeof(event)), camera")
mouseMotion(x::Number, y::Number) = log("mouseMotion: x=$x, y=$y")
mouseMotionAbs(event::Ahorn.eventMouse, camera::Ahorn.Camera) = log("mouseMotionAbs: event=$(typeof(event)), camera")
mouseMotionAbs(x::Number, y::Number) = log("mouseMotionAbs: x=$x, y=$y")

# Ahorn.eventKey = Gtk.GdkEventKey
keyboard(event::Ahorn.eventKey) = log("keyboard: event=$(event.hardware_keycode)")

end


text:
data = reshape([4.939	2.432;
    7.976	7.237;
    10.976	5.826;
    13.976	3.874;
    16.976	6.697;
    19.976	7.778;
    22.976	6.937;
    26.012	5.736],2,8)

dist = 3.5/8


text:
# Currently there is an inefficiency in Julia v0.5 and possibly v0.6 with ntuple
# that will hopefully be fixed in the next version. To avoid
# slowdown, we implement a tuple constructor independent of Julia's
# metaprogramming.
#
# The method below is used for building the multi-dimensional
# derivative matrix:
@generated function make_cartesian_index{D}(i::Int, arr1::Int, arr2::CartesianIndex{D})
	levs = [:($q == i ? arr1 : arr2[$q]) for q=1:D]
	quote
		$(Expr(:meta, :inline))
		CartesianIndex{D}($(levs...))
	end
end

# Threshold method for full matrices, producing a sparse one
function threshold{T<:Real}(mat::Array{T, 2}; abs_tol=eps(T))
	for i in eachindex(mat)
		val = mat[i]
		if abs(val) < abs_tol
			mat[i] = 0
		end
	end
	return sparse(mat)
end

# Threshold method for sparse matrices
function threshold{T<:Real}(mat::SparseMatrixCSC{T, Int}; abs_tol=eps(T))
	(I, J, V) = findnz(mat)
	for (index, val) in enumerate(V)
		if abs(val) < abs_tol
			V[index] = zero(T)
		end
	end
	return dropzeros!(sparse(I, J, V))
end