text:
using SpecialFunctions

@testset "math" begin
    @testset "log10" begin
        @test testf(a->log10.(a), Float32[100])
    end

    @testset "pow" begin
        for T in (Float16, Float32, Float64, ComplexF32, ComplexF64)
            range = (T<:Integer) ? (T(5):T(10)) : T
            @test testf((x,y)->x.^y, rand(Float32, 1), rand(range, 1))
            @test testf((x,y)->x.^y, rand(Float32, 1), -rand(range, 1))
        end
    end

    @testset "isinf" begin
      for x in (Inf32, Inf, NaN32, NaN)
        @test testf(x->isinf.(x), [x])
      end
    end

    @testset "isnan" begin
      for x in (Inf32, Inf, NaN32, NaN)
        @test testf(x->isnan.(x), [x])
      end
    end

    for op in (exp, angle, exp2, exp10,)
        @testset "$op" begin
            for T in (Float16, Float32, Float64)
                @test testf(x->op.(x), rand(T, 1))
                @test testf(x->op.(x), -rand(T, 1))
            end
        end
    end

    for op in (expm1,)
        @testset "$op" begin
            # FIXME: add expm1(::Float16) to Base
            for T in (Float32, Float64)
                @test testf(x->op.(x), rand(T, 1))
                @test testf(x->op.(x), -rand(T, 1))
            end
        end
    end

    for op in (exp, abs, abs2, angle, log)
        @testset "Complex - $op" begin
            for T in (ComplexF16, ComplexF32, ComplexF64)
                @test testf(x->op.(x), rand(T, 1))
                @test testf(x->op.(x), -rand(T, 1))
            end

        end
    end
    @testset "mod and rem" begin
        for T in (Float16, Float32, Float64)
            @test testf(a->rem.(a, T(2)), T[0, 1, 1.5, 2, -1])
            @test testf(a->rem.(a, T(2), RoundNearest), T[0, 1, 1.5, 2, -1])
            @test testf(a->mod.(a, T(2)), T[0, 1, 1.5, 2, -1])
        end
    end

    @testset "rsqrt" begin
        # GPUCompiler.jl#173: a CUDA-only device function fails to validate
        function kernel(a)
            a[] = CUDA.rsqrt(a[])
            return
        end

        # make sure this test uses an actual device function
        @test_throws ErrorException kernel(ones(1))

        for T in (Float16, Float32)
            a = CuArray{T}([4])
            @cuda kernel(a)
            @test Array(a) == [0.5]
        end
    end

    @testset "fma" begin
        for T in (Float16, Float32, Float64)
            @test testf((x,y,z)->fma.(x,y,z), rand(T, 1), rand(T, 1), rand(T, 1))
            @test testf((x,y,z)->fma.(x,y,z), rand(T, 1), -rand(T, 1), -rand(T, 1))
        end
    end

    # something from SpecialFunctions.jl
    @testset "erf" begin
        @test testf(a->SpecialFunctions.erf.(a), Float32[1.0])
    end
    @testset "loggamma" begin
        @test testf(a->SpecialFunctions.loggamma.(a), Float32[1.0])
    end

    @testset "exp" begin
        # JuliaGPU/CUDA.jl#1085: exp uses Base.sincos performing a global CPU load
        @test testf(x->exp.(x), [1e7im])
    end

    @testset "fastmath" begin
        # libdevice provides some fast math functions
        a(x) = cos(x)
        b(x) = @fastmath cos(x)
        @test Array(map(a, cu([0.1,0.2]))) ≈ Array(map(b, cu([0.1,0.2])))

        # JuliaGPU/CUDA.jl#1352: some functions used to fall back to libm
        f(x) = log1p(x)
        g(x) = @fastmath log1p(x)
        @test Array(map(f, cu([0.1,0.2]))) ≈ Array(map(g, cu([0.1,0.2])))
    end

    @testset "byte_perm" begin
        bytes = UInt32[i for i in 1:8]
        x = bytes[4]<<24 | bytes[3]<<16 | bytes[2]<<8 | bytes[1]<<0
        y = bytes[8]<<24 | bytes[7]<<16 | bytes[6]<<8 | bytes[5]<<0
        sel = UInt32[4, 2, 4, 2]
        code = sel[4]<<12 | sel[3]<<8 | sel[2]<<4 | sel[1]<<0
        r = bytes[sel[4]+1]<<24 | bytes[sel[3]+1]<<16 | bytes[sel[2]+1]<<8 | bytes[sel[1]+1]<<0

        function kernel1(a)
            a[3] = CUDA.byte_perm(a[1], a[2], code % Int32)
            return
        end
        function kernel2(a)
            a[3] = CUDA.byte_perm(a[1], a[2], code % UInt16)
            return
        end

        for T in [UInt32, Int32]
            a = CuArray{T}([x, y, 0])
            @cuda kernel1(a)
            @test Array(a)[3] == r
            a = CuArray{T}([x, y, 0])
            @cuda kernel2(a)
            @test Array(a)[3] == r
        end
    end
end


text:
"""
Abstract type representing a function (mapping) that provides a cache and an in-place
evaluation for performance. This is the type to be used in the [`lazy_map`](@ref) function.

Derived types must implement the following method:

- [`evaluate!(cache,k,x...)`](@ref)

and optionally these ones:

- [`return_cache(k,x...)`](@ref)
- [`return_type(k,x...)`](@ref)

The mapping interface can be tested with the [`test_map`](@ref) function.

Note that most of the functionality implemented in terms of this interface
relies in duck typing. That is, it is not strictly needed to work with types
that inherit from `Map`. This is specially useful in order to accommodate
existing types into this framework without the need to implement a wrapper type
that inherits from `Map`. For instance, a default implementation is available
for `Function` objects.  However, we recommend that new types inherit from `Map`.

"""
abstract type Map <: GridapType end

"""
    return_cache(f,x...)

Returns the `cache` needed to lazy_map mapping `f` with arguments
of the same type as the objects in `x`.
This function returns `nothing` by default, i.e., no cache.
"""
return_cache(f,x...) = nothing

"""
    evaluate!(cache,f,x...)

Applies the mapping `f` at the arguments `x...` using
the scratch data provided in the given `cache` object. The `cache` object
is built with the [`return_cache`](@ref) function using arguments of the same type as in `x`.
In general, the returned value `y` can share some part of its state with the `cache` object.
If the result of two or more calls to this function need to be accessed simultaneously
(e.g., in multi-threading), create and use several `cache` objects (e.g., one cache
per thread).
"""
@noinline function evaluate!(cache,f,x...)
  @abstractmethod """
  Method

      Gridap.Arrays.evaluate!(cache,f,x...)

  has not been defined for the requested types (see stack trace).
  """
end

"""
    return_type(f,x...)

Returns the type of the result of calling mapping `f` with
arguments of the types of the objects `x`.
"""
return_type(f,x...) = typeof(return_value(f,x...))

return_value(f,x...) = evaluate(f,testargs(f,x...)...)

"""
    testargs(f,x...)

The default implementation of this function is `testargs(f,x...) = x`.
One can overload it in order to use `lazy_map` with 0-length array
and maps with non-trivial domains.
"""
testargs(f,x...) = x

"""
    evaluate(f,x...)

evaluates the mapping `f` at the arguments in `x` by creating a temporary cache
internally. This functions is equivalent to
```jl
cache = return_cache(f,x...)
evaluate!(cache,f,x...)
```
"""
function evaluate(f,x...)
  c = return_cache(f,x...)
  y = evaluate!(c,f,x...)
end

(m::Map)(x...) = evaluate(m,x...)

# Default implementation for Function
evaluate!(cache,f::Function,x...) = f(x...)

# Default implementation for constructors
evaluate!(cache,::Type{f},x...) where f = f(x...)

# Testing the interface
"""
    test_map(y,f,x...;cmp=(==))

Function used to test if the mapping `f` has been
implemented correctly. `f` is a `Map` sub-type, `x` is a tuple in the domain of the
mapping and `y` is the expected result. Function `cmp` is used to compare
the computed result with the expected one. The checks are done with the `@test`
macro.
"""
function test_map(y,f,x...;cmp=(==))
  z = evaluate(f,x...)
  @test cmp(z,y)
  @test typeof(z) == return_type(f,x...)
  cache = return_cache(f,x...)
  z = evaluate!(cache,f,x...)
  @test cmp(z,y)
  z = evaluate!(cache,f,x...)
  @test cmp(z,y)
  true
end

# Broadcast Functions

"""
    Broadcasting(f)

Returns a mapping that represents the "broadcasted" version of the
function `f`.

# Example

```jldoctest
using Gridap.Arrays

a = [3,2]
b = [2,1]

bm = Broadcasting(+)

c = evaluate(bm,a,b)

println(c)

# output
[5, 3]
```
"""
struct Broadcasting{F} <: Map
  f::F
end

return_cache(f::Broadcasting,x...) = nothing

evaluate!(cache,f::Broadcasting,x...) = broadcast(f.f,x...)

function return_value(f::Broadcasting,x...)
  broadcast( (y...) -> f.f(testargs(f.f,y...)...), x... )
end

function evaluate!(cache,f::Broadcasting,x::Union{Number,AbstractArray{<:Number}}...)
  r = _prepare_cache!(cache,x...)
  a = r.array
  broadcast!(f.f,a,x...)
  a
end

function evaluate!(cache,f::Broadcasting,x::AbstractArray{<:Number})
  setsize!(cache,size(x))
  a = cache.array
  @inbounds for i in eachindex(x)
    a[i] = f.f(x[i])
  end
  a
end

function evaluate!(cache,f::Broadcasting,args::Number...)
  f.f(args...)
end

function return_value(f::Broadcasting,x::Number...)
  return_value(f.f,x...)
end

function return_cache(f::Broadcasting,x::Number...)
  nothing
end

function return_value(f::Broadcasting,x::Union{Number,AbstractArray{<:Number}}...)
  s = map(_size_zero,x)
  bs = Base.Broadcast.broadcast_shape(s...)
  T = return_type(f.f,map(testitem,x)...)
  r = fill(testvalue(T),bs)
  r
end

function _size_zero(a)
  s = size(a)
  if length(a) == 0
    r = map(i-> (i==0 ? 1 : i) ,s)
  else
    r = s
  end
  r
end
_size_zero(a::Number) = (1,)

function return_cache(f::Broadcasting,x::Union{Number,AbstractArray{<:Number}}...)
  s = map(_size,x)
  bs = Base.Broadcast.broadcast_shape(s...)
  T = return_type(f.f,map(testitem,x)...)
  r = fill(testvalue(T),bs)
  cache = CachedArray(r)
  _prepare_cache!(cache,x...)
  cache
end

function _prepare_cache!(c,x...)
  s = map(_size,x)
  bs = Base.Broadcast.broadcast_shape(s...)
  if bs != size(c)
    setsize!(c,bs)
  end
  c
end

_size(a) = size(a)
_size(a::Number) = (1,)

"""
    OperationMap(f,args)

Returns a mapping that represents the result of applying the function `f`
to the arguments in the tuple `args`.
That is, `OperationMap(f,args)(x...)` is formally defined as
`f(map(a->a(x...),args)...)`
"""
struct OperationMap{K,L} <: Map
  k::K
  l::L
  function OperationMap(k,l)
    new{typeof(k),typeof(l)}(k,l)
  end
end

function return_cache(c::OperationMap,x...)
  cl = map(fi -> return_cache(fi,x...),c.l)
  lx = map(fi -> return_value(fi,x...),c.l)
  ck = return_cache(c.k,lx...)
  ck, cl
end

function evaluate!(cache,c::OperationMap,x...)
  ck, cf = cache
  lx = map((ci,fi) -> evaluate!(ci,fi,x...),cf,c.l)
  evaluate!(ck,c.k,lx...)
end

# Operations

"""
    Operation(op)

Retuns the map that results after applying an operation `f` over a set of map(s) `args`.
That is `Operation(f)(args)(x...)` is formally defined as
`f(map(a->a(x...),args)...)`.

# Example

```jldoctest
using Gridap.Arrays

fa(x) = x.*x
fb(x) = sqrt.(x)

x = collect(0:5)

fab = Operation(fa)(fb)
c = evaluate(fab,x)

println(c)

# output
[0.0, 1.0, 2.0, 3.0, 4.0, 5.0]
```
"""
struct Operation{T} <: Map
  op::T
end

evaluate!(cache,op::Operation,args::Map...) = OperationMap(op.op,args)
evaluate!(cache,op::Operation,args::Function...) = OperationMap(op.op,args)

"""
"""
function inverse_map(f)
  @unreachable """\n
  Function inverse_map is not implemented yet for objects of type $(typeof(f))
  """
end

ext: jl


text:
function rrule_chain_impl(ch::Type{Chain{Tl}}, x::Type{Tx}) where {Tl,Tx}
    # The number of layers
    N = length(Tl.parameters)

    # output of intermediate layers (x_n+1 = l_n(x_n))
    x_vars  = [Symbol("Y$(i-1)") for i=1:N+1]
    x_vars[1] = Symbol("x")
    # pullback of intermediate layers (pb_n = pb(l_n))
    pb_vars = [Symbol("∂Ȳ$i") for i=1:N]

    expr = Expr(:block)
    # Forward pass
    for i=1:N
        l = :(getfield(c.layers,$i))
        expr_i = :(($(x_vars[i+1]), $(pb_vars[i])) = rrule($l, $(x_vars[i])));
        push!(expr.args, expr_i)
    end

    # Work around closures not working in generated funcs
    push!(expr.args, :(chain_pb = PbClosure(tuple($(pb_vars...)))))

    push!(expr.args, :(return ($(x_vars[end]), chain_pb)))

    return expr
end

@generated function ChainRulesCore.rrule(c::Chain{C}, x::T) where {C,T}
    return rrule_chain_impl(c, x)
end


text:
@testset "SplitChannels" begin
    @test (SplitChannels <: Augmentor.AffineOperation) == false
    @test (SplitChannels <: Augmentor.ImageOperation) == false
    @test (SplitChannels <: Augmentor.Operation) == true

    @testset "constructor" begin
        @test typeof(@inferred(SplitChannels())) <: SplitChannels <: Augmentor.Operation
        @test str_show(SplitChannels()) == "Augmentor.SplitChannels()"
        @test str_showconst(SplitChannels()) == "SplitChannels()"
        @test str_showcompact(SplitChannels()) == "Split colorant into its color channels"
    end
    @testset "eager" begin
        @test Augmentor.supports_eager(SplitChannels) === false
        @test_throws MethodError Augmentor.applyeager(SplitChannels(), rand(2,2))
        res1 = channelview(reshape(rect, 1, 2, 3))
        res2 = channelview(rgb_rect)
        imgs = [
            (rect, res1),
            (Augmentor.prepareaffine(rect), res1),
            (OffsetArray(rect, -2, -1), res1),
            (view(rect, IdentityRange(1:2), IdentityRange(1:3)), res1),
            (rgb_rect, res2),
            (Augmentor.prepareaffine(rgb_rect), res2),
            (OffsetArray(rgb_rect, -2, -1), res2),
            (view(rgb_rect, IdentityRange(1:2), IdentityRange(1:3)), res2),
        ]
        @testset "single image" begin
            for (img_in, img_out) in imgs
                res = @inferred(Augmentor.applyeager(SplitChannels(), img_in))
                @test res == img_out
                @test typeof(res) == typeof(img_out)
            end
        end
        @testset "multiple images" begin
            for (img_in1, img_out1) in imgs, (img_in2, img_out2) in imgs
                img_in = (img_in1, img_in2)
                img_out = (img_out1, img_out2)
                res = @inferred(Augmentor.applyeager(SplitChannels(), img_in))
                @test res == img_out
                @test typeof(res) == typeof(img_out)
            end
        end
    end
    @testset "affine" begin
        @test Augmentor.supports_affine(SplitChannels) === false
    end
    @testset "affineview" begin
        @test Augmentor.supports_affineview(SplitChannels) === false
    end
    @testset "lazy" begin
        @test Augmentor.supports_lazy(SplitChannels) === true
        @test @inferred(Augmentor.supports_lazy(typeof(SplitChannels()))) === true
        @test_throws MethodError Augmentor.applylazy(SplitChannels(), rand(2,2))
        f1 = img->reshape(channelview(img), 1, 2, 3)
        f2 = img->Augmentor.plain_indices(channelview(img))
        imgs = [
            (rect, f1),
            (Augmentor.prepareaffine(rect), f1),
            (OffsetArray(rect, -2, -1), f1),
            (view(rect, IdentityRange(1:2), IdentityRange(1:3)), f1),
            (rgb_rect, f2),
            (Augmentor.prepareaffine(rgb_rect), f2),
            (OffsetArray(rgb_rect, -2, -1), f2),
            (view(rgb_rect, IdentityRange(1:2), IdentityRange(1:3)), f2),
        ]
        @testset "single image" begin
            for (img_in, f) in imgs
                img_out = f(img_in)
                res = @inferred(Augmentor.applylazy(SplitChannels(), img_in))
                @test res == img_out
                @test res == Augmentor.applyeager(SplitChannels(), img_in)
                @test typeof(res) == typeof(img_out)
            end
        end
        @testset "multiple images" begin
            for (img_in1, f1) in imgs, (img_in2, f2) in imgs
                img_out1 = f1(img_in1)
                img_out2 = f2(img_in2)
                img_in = (img_in1, img_in2)
                img_out = (img_out1, img_out2)
                res = @inferred(Augmentor.applylazy(SplitChannels(), img_in))
                @test res == img_out
                @test typeof(res) == typeof(img_out)
            end
        end
    end
    @testset "view" begin
        @test Augmentor.supports_view(SplitChannels) === false
    end
    @testset "stepview" begin
        @test Augmentor.supports_stepview(SplitChannels) === false
    end
    @testset "permute" begin
        @test Augmentor.supports_permute(SplitChannels) === false
    end
end

# --------------------------------------------------------------------

@testset "CombineChannels" begin
    @test (CombineChannels <: Augmentor.AffineOperation) == false
    @test (CombineChannels <: Augmentor.ImageOperation) == false
    @test (CombineChannels <: Augmentor.Operation) == true

    @testset "constructor" begin
        @test_throws MethodError CombineChannels()
        @test_throws MethodError CombineChannels(Float64)
        @test typeof(@inferred(CombineChannels(RGB))) <: CombineChannels <: Augmentor.Operation
        @test typeof(@inferred(CombineChannels(RGB{N0f8}))) <: CombineChannels <: Augmentor.Operation
        @test str_show(CombineChannels(RGB)) == "Augmentor.CombineChannels(RGB{Any})"
        @test str_show(CombineChannels(Gray{N0f8})) == "Augmentor.CombineChannels(Gray{N0f8})"
        @test str_showconst(CombineChannels(RGB{N0f8})) == "CombineChannels(RGB{N0f8})"
        @test str_showcompact(CombineChannels(Gray)) == "Combine color channels into colorant Gray{Any}"
    end
    @testset "eager" begin
        @test Augmentor.supports_eager(CombineChannels) === false
        @test_throws MethodError Augmentor.applyeager(CombineChannels(RGB), rand(RGB{N0f8},4,4))
        @test_throws MethodError Augmentor.applyeager(CombineChannels(Gray), rand(Gray{N0f8},4,4))
        @test_throws ArgumentError Augmentor.applyeager(CombineChannels(Gray), rand(4,4))
        @test_throws MethodError Augmentor.applyeager(CombineChannels(RGB), (rand(3,4,4), rand(RGB{N0f8},4,4)))
        @test_throws MethodError Augmentor.applyeager(CombineChannels(Gray), (rand(1,4,4), rand(Gray{N0f8},4,4)))
        @test_throws ArgumentError Augmentor.applyeager(CombineChannels(Gray), (rand(1,4,4), rand(4,4)))
        rect_split = reshape(channelview(rect), 1, 2, 3)
        rgb_rect_split = channelview(rgb_rect)
        imgs = [
            (rect_split, rect),
            (reshape(channelview(Augmentor.prepareaffine(rect)), 1, 2, 3), rect),
            (OffsetArray(rect_split, 0, -2, -1), rect),
            (view(rect_split, 1:1, IdentityRange(1:2), IdentityRange(1:3)), rect),
        ]
        rgb_imgs = [
            (rgb_rect_split, rgb_rect),
            (channelview(Augmentor.prepareaffine(rgb_rect)), rgb_rect),
            (OffsetArray(rgb_rect_split, 0, -2, -1), rgb_rect),
            (view(rgb_rect_split, IdentityRange(1:3), IdentityRange(1:2), IdentityRange(1:3)), rgb_rect),
        ]
        @testset "single image" begin
            for (img_in, img_out) in imgs
                res = @inferred(Augmentor.applyeager(CombineChannels(Gray), img_in))
                @test res == img_out
                @test typeof(res) == typeof(img_out)
                @test typeof(res) <: Array{Gray{N0f8}}
                res = @inferred(Augmentor.applyeager(CombineChannels(Gray{N0f8}), img_in))
                @test res == img_out
                @test typeof(res) == typeof(img_out)
                @test typeof(res) <: Array{Gray{N0f8}}
            end
            for (img_in, img_out) in rgb_imgs
                res = @inferred(Augmentor.applyeager(CombineChannels(RGB), img_in))
                @test res == img_out
                @test typeof(res) == typeof(img_out)
                @test typeof(res) <: Array{RGB{N0f8}}
                res = @inferred(Augmentor.applyeager(CombineChannels(RGB{N0f8}), img_in))
                @test res == img_out
                @test typeof(res) == typeof(img_out)
                @test typeof(res) <: Array{RGB{N0f8}}
            end
        end
        @testset "multiple images" begin
            for (img_in1, img_out1) in imgs, (img_in2, img_out2) in imgs
                img_in = (img_in1, img_in2)
                img_out = (img_out1, img_out2)
                res = @inferred(Augmentor.applyeager(CombineChannels(Gray), img_in))
                @test res == img_out
                @test typeof(res) == typeof(img_out)
                res = @inferred(Augmentor.applyeager(CombineChannels(Gray{N0f8}), img_in))
                @test res == img_out
                @test typeof(res) == typeof(img_out)
            end
            for (img_in1, img_out1) in rgb_imgs, (img_in2, img_out2) in rgb_imgs
                img_in = (img_in1, img_in2)
                img_out = (img_out1, img_out2)
                res = @inferred(Augmentor.applyeager(CombineChannels(RGB), img_in))
                @test res == img_out
                @test typeof(res) == typeof(img_out)
                res = @inferred(Augmentor.applyeager(CombineChannels(RGB{N0f8}), img_in))
                @test res == img_out
                @test typeof(res) == typeof(img_out)
            end
        end
    end
    @testset "affine" begin
        @test Augmentor.supports_affine(CombineChannels) === false
    end
    @testset "affineview" begin
        @test Augmentor.supports_affineview(CombineChannels) === false
    end
    @testset "lazy" begin
        @test Augmentor.supports_lazy(CombineChannels) === true
        @test @inferred(Augmentor.supports_lazy(typeof(CombineChannels(Gray)))) === true
        @test_throws MethodError Augmentor.applylazy(CombineChannels(RGB), rand(RGB{N0f8},4,4))
        @test_throws MethodError Augmentor.applylazy(CombineChannels(Gray), rand(Gray{N0f8},4,4))
        @test_throws ArgumentError Augmentor.applylazy(CombineChannels(Gray), rand(4,4))
        @test_throws MethodError Augmentor.applylazy(CombineChannels(RGB), (rand(3,4,4), rand(RGB{N0f8},4,4)))
        @test_throws MethodError Augmentor.applylazy(CombineChannels(Gray), (rand(1,4,4), rand(Gray{N0f8},4,4)))
        @test_throws ArgumentError Augmentor.applylazy(CombineChannels(Gray), (rand(1,4,4), rand(4,4)))
        rect_split = reshape(channelview(rect), 1, 2, 3)
        rgb_rect_split = channelview(rgb_rect)
        f1 = (T, img) -> colorview(T, reshape(img, 2, 3))
        f2 = (T, img) -> colorview(T, Augmentor.plain_indices(img))
        imgs = [
            (rect_split),
            (reshape(channelview(Augmentor.prepareaffine(rect)), 1, 2, 3)),
            (OffsetArray(rect_split, 0, -2, -1)),
            (view(rect_split, 1:1, IdentityRange(1:2), IdentityRange(1:3))),
        ]
        rgb_imgs = [
            (rgb_rect_split),
            (channelview(Augmentor.prepareaffine(rgb_rect))),
            (OffsetArray(rgb_rect_split, 0, -2, -1)),
            (view(rgb_rect_split, 1:3, IdentityRange(1:2), IdentityRange(1:3))),
        ]
        @testset "single image" begin
            for img_in in imgs, T in (Gray, Gray{N0f8})
                img_out = f1(T, img_in)
                res = @inferred(Augmentor.applylazy(CombineChannels(T), img_in))
                @test res == img_out
                @test res == Augmentor.applyeager(CombineChannels(T), img_in)
                @test typeof(res) == typeof(img_out)
                @test eltype(res) == Gray{N0f8}
            end
            for img_in in rgb_imgs, T in (RGB, RGB{N0f8})
                img_out = f2(T, img_in)
                res = @inferred(Augmentor.applylazy(CombineChannels(T), img_in))
                @test res == img_out
                @test res == Augmentor.applyeager(CombineChannels(T), img_in)
                @test typeof(res) == typeof(img_out)
                @test eltype(res) == RGB{N0f8}
            end
        end
        @testset "multiple images" begin
            for img_in1 in imgs, img_in2 in imgs, T in (Gray, Gray{N0f8})
                img_in = (img_in1, img_in2)
                img_out = f1.(T, img_in)
                res = @inferred(Augmentor.applylazy(CombineChannels(T), img_in))
                @test res == img_out
                @test typeof(res) == typeof(img_out)
                @test map(eltype,res) == (Gray{N0f8}, Gray{N0f8})
            end
            for img_in1 in rgb_imgs, img_in2 in rgb_imgs, T in (RGB, RGB{N0f8})
                img_in = (img_in1, img_in2)
                img_out = f2.(T, img_in)
                res = @inferred(Augmentor.applylazy(CombineChannels(T), img_in))
                @test res == img_out
                @test typeof(res) == typeof(img_out)
                @test map(eltype,res) == (RGB{N0f8}, RGB{N0f8})
            end
        end
    end
    @testset "view" begin
        @test Augmentor.supports_view(CombineChannels) === false
    end
    @testset "stepview" begin
        @test Augmentor.supports_stepview(CombineChannels) === false
    end
    @testset "permute" begin
        @test Augmentor.supports_permute(CombineChannels) === false
    end
end


text:
const Recalls = try
    Base.require(Base.PkgId(Base.UUID(0x30af7cf3eb4344c7afa733725b72a81e), "Recalls"))
catch
    nothing
end

istracing() = false

macro trace(args...)
    Recalls === nothing && return nothing
    quote
        if $istracing()
            $Recalls.@note($(args...),)
        end
        nothing
    end |> esc
end
# TODO: Use Preferences.jl to controll if Recalls.jl should be loaded? Since
# `args` may affect how much variables are captured in closures, the existence
# of Recalls.jl changes the program slightly.

function enable_tracing()
    if Recalls === nothing
        error(
            "Tracing requires Recalls.jl to be installed in the same environment.",
            " See: https://github.com/tkf/Recalls.jl",
        )
    end
    prev = istracing()
    @eval istracing() = true
    return prev
end

function disable_tracing()
    prev = istracing()
    @eval istracing() = false
    return prev
end

should_limit_retries() = istracing() || use_anchors() || use_retrylimit()

use_retrylimit() = false

function enable_retrylimit()
    prev = use_retrylimit()
    @eval use_retrylimit() = true
    return prev
end

function disable_retrylimit()
    prev = use_retrylimit()
    @eval use_retrylimit() = false
    return prev
end


text:
module Plotting

using DSP, FFTW
using Plots

export set_plot, plot_frec, plot_wave, plot_spectrogram, plot_zplane, plot_phase

# Sets backend renderer - PyPlot
function set_plot()
	pyplot()
end

# Plots wave in time
function plot_wave(s, fs;
		   plot_title="Modem Dialing",
		   image_name="wave")

	x_axis = (0:(length(s) - 1))./fs

	# Plot it
	p = plot(x_axis, s,
		 title="$plot_title - fs: $fs",
		 xlabel="Time [s]",
		 legend=false)

	# Save it in .png
	savefig(p, "$image_name.png")
end

# Plots wave in frecuency
function plot_frec(s, fs;
		   plot_title="Modem Dialing",
		   image_name="wave-frec")

	# FFT of the input wave
	X = abs.(fft(s))

	# Shift the returned spectrum
	Xshift = fftshift(X)

	X = X[1:div(end,2)]

	# In Hertz
	x_axis = (0:(1/length(X)):1) .* (fs/2)

	# In Omega (w)
	x_axis_omega = (-1:(1/length(Xshift)):1) .* fs

	x_neg = findfirst(x -> x >= -(fs/2), x_axis_omega)
	x_pos = findfirst(x -> x >= (fs/2), x_axis_omega)

	# Frecuency plot
	p = plot(x_axis, X,
		 title="$plot_title Frecuency Domain",
		 xlabel="Freq [Hz]",
		 legend=false);

	# Save it in .png
	savefig(p, "$image_name.png")

	# Omega plot
	p = plot(x_axis_omega[x_neg:x_pos], Xshift,
		 title="$plot_title Omega Frecuency Domain",
		 xlabel="Freq [w]",
		 legend=false);

	# Save it in .png
	savefig(p, "$image_name-omega.png")
end

# Plots the spectrogram
function plot_spectrogram(s, fs;
			  win=tukey(256, 0.5),
			  plot_title="Modem Dialing",
			  image_name="spec")

	# Transform it into a Vector
	s = vec(s)

	# Obtains the spectrogram
	spec = spectrogram(s, length(win), fs=fs, window=win)

	t, fr, pow = spec.time, spec.freq, spec.power

	# Rearrenge the negative frequencies
	neg_freq = findfirst(x -> x < 0, fr)

	if neg_freq !== nothing
		neg_pos = neg_freq - 1
		freq = [spec.freq[neg_freq:end]; spec.freq[1:neg_pos]]
		pow = [spec.power[neg_freq:end, :]; spec.power[1:neg_pos, :]]
	end

	# Plot it (HeatMap)
	hm = heatmap(t, fr, pow .+ eps() .|> log; seriescolor=:bluesreds,
		     title="$plot_title Spectrogram", xlabel="Time [s]", ylabel="Freq [Hz]")

	# Save it in .png
	savefig(hm, "$image_name.png")
end

function plot_zplane(h, fc)

	fil = PolynomialRatio(h, [1])

	zpg = ZeroPoleGain(fil)

	zs = zpg.z
	ps = zpg.p

	length(zs) > length(ps) && (push!(ps, 0))

	scatter(real.(zs), imag.(zs),
		color=[:black],
		marker=:circle,
		label="Zero",
		title="Zero-Pole diagram - Filter: $fc Hz")

	scatter!(real.(ps), imag.(ps),
		color=[:red],
		marker=:diamond,
		label="Pole")

	savefig("zero-pole-filter-$fc.png")
end

function plot_phase(h, fc, fs)

	X = fft(h)

	X = X[1:div(end, 2)]

	ph = unwrap(angle.(X))

	x_axis = (0:(1/length(ph)):1) .* (fs/2);

	p = plot(x_axis, ph,
		 xlabel="Frec [Hz]",
		 ylabel="Phase [rad]",
		 title="Phase diagram - Filter: $fc Hz",
		 legend=false,
		 xticks=(0:500:(fs/2)))

	savefig(p, "phase-filter-$fc.png")
end

end # module

max_stars_repo_name: PatricioIribarneCatella/TP-Telefono

max_issues_repo_name: PatricioIribarneCatella/TP-Telefono

max_forks_repo_name: PatricioIribarneCatella/TP-Telefono


text:
mutable struct SVMSurrogate{X,Y,M,L,U} <: AbstractSurrogate
    x::X
    y::Y
    model::M
    lb::L
    ub::U
end

function SVMSurrogate(x,y,lb::Number,ub::Number)
    xn = reshape(x,length(x),1)
    model = LIBSVM.fit!(SVC(),xn,y)
    SVMSurrogate(xn,y,model,lb,ub)
end

function (svmsurr::SVMSurrogate)(val::Number)
    return LIBSVM.predict(svmsurr.model,[val])
end

"""
SVMSurrogate(x,y,lb,ub)

Builds SVM surrogate.
"""
function SVMSurrogate(x,y,lb,ub)
    X = Array{Float64,2}(undef,length(x),length(x[1]))
    for j = 1:length(x)
        X[j,:] = vec(collect(x[j]))
    end
    model = LIBSVM.fit!(SVC(),X,y)
    SVMSurrogate(x,y,model,lb,ub)
end

function (svmsurr::SVMSurrogate)(val)
    n = length(val)
    return LIBSVM.predict(svmsurr.model,reshape(collect(val),1,n))[1]
end

function add_point!(svmsurr::SVMSurrogate,x_new,y_new)
     if length(svmsurr.lb) == 1
         #1D
         svmsurr.x = vcat(svmsurr.x,x_new)
         svmsurr.y = vcat(svmsurr.y,y_new)
         svmsurr.model = LIBSVM.fit!(SVC(),reshape(svmsurr.x,length(svmsurr.x),1),svmsurr.y)
     else
         n_previous = length(svmsurr.x)
         a = vcat(svmsurr.x,x_new)
         n_after = length(a)
         dim_new = n_after - n_previous
         n = length(svmsurr.x)
         d = length(svmsurr.x[1])
         tot_dim = n + dim_new
         X = Array{Float64,2}(undef,tot_dim,d)
         for j = 1:n
             X[j,:] = vec(collect(svmsurr.x[j]))
         end
         if dim_new == 1
             X[n+1,:] = vec(collect(x_new))
         else
             i = 1
             for j = n+1:tot_dim
                 X[j,:] = vec(collect(x_new[i]))
                 i = i + 1
             end
         end
         svmsurr.x = vcat(svmsurr.x,x_new)
         svmsurr.y = vcat(svmsurr.y,y_new)
         svmsurr.model = LIBSVM.fit!(SVC(),X,svmsurr.y)
     end
     nothing
 end


text:
Base.@kwdef struct PardisoJL <: SciMLLinearSolveAlgorithm
    nprocs::Union{Int, Nothing} = nothing
    solver_type::Union{Int, Pardiso.Solver, Nothing} = nothing
    matrix_type::Union{Int, Pardiso.MatrixType, Nothing} = nothing
    iparm::Union{Vector{Tuple{Int,Int}}, Nothing} = nothing
    dparm::Union{Vector{Tuple{Int,Int}}, Nothing} = nothing
end

MKLPardisoFactorize(;kwargs...) = PardisoJL(;kwargs...)
MKLPardisoIterate(;kwargs...) = PardisoJL(;kwargs...)

# TODO schur complement functionality

function init_cacheval(alg::PardisoJL, A, b, u, Pl, Pr, maxiters, abstol, reltol, verbose)
    @unpack nprocs, solver_type, matrix_type, iparm, dparm = alg

    if A isa DiffEqArrayOperator
        A = A.A
    end

    solver =
    if Pardiso.PARDISO_LOADED[]
        solver = Pardiso.PardisoSolver()
        solver_type !== nothing && Pardiso.set_solver!(solver, solver_type)

        solver
    else
        solver = Pardiso.MKLPardisoSolver()
        nprocs !== nothing && Pardiso.set_nprocs!(solver, nprocs)

        solver
    end

    Pardiso.pardisoinit(solver) # default initialization

    if matrix_type !== nothing
        Pardiso.set_matrixtype!(solver, matrix_type)
    else
        if eltype(A) <: Real
            Pardiso.set_matrixtype!(solver, Pardiso.REAL_NONSYM)
        elseif eltype(A) <: Complex
            Pardiso.set_matrixtype!(solver, Pardiso.COMPLEX_NONSYM)
        else
            error("Number type not supported by Pardiso")
        end
    end
    verbose && Pardiso.set_msglvl!(solver, Pardiso.MESSAGE_LEVEL_ON)

    # pass in vector of tuples like [(iparm::Int, key::Int) ...]
    if iparm !== nothing
        for i in length(iparm)
            Pardiso.set_iparm!(solver, iparm[i]...)
        end
    end

    if dparm !== nothing
        for i in length(dparm)
            Pardiso.set_dparm!(solver, dparm[i]...)
        end
    end

    Pardiso.set_phase!(solver, Pardiso.ANALYSIS)
    Pardiso.pardiso(solver, u, A, b)

    return solver
end

function SciMLBase.solve(cache::LinearCache, alg::PardisoJL; kwargs...)
    @unpack A, b, u = cache
    if A isa DiffEqArrayOperator
        A = A.A
    end

    if cache.isfresh
        Pardiso.set_phase!(cache.cacheval, Pardiso.NUM_FACT)
        Pardiso.pardiso(cache.cacheval, cache.u, cache.A, cache.b)
    end

    Pardiso.set_phase!(cache.cacheval, Pardiso.SOLVE_ITERATIVE_REFINE)
    Pardiso.pardiso(cache.cacheval, u, A, b)

    return SciMLBase.build_linear_solution(alg,cache.u,nothing,cache)
end

# Add finalizer to release memory
# Pardiso.set_phase!(cache.cacheval, Pardiso.RELEASE_ALL)

export PardisoJL, MKLPardisoFactorize, MKLPardisoIterate


text:
using HttpServer, DataFrames, Calendar, SQLite, Mustache, HttpCommon
require("templates.jl")

println("connecting to db")
conn = connect("blog.db")
posts = query("select * from posts",conn)
posts["url"] = map(encodeURI,posts["title"])

function myhandler(req,res)
  # Index Page
  if req.resource == "/" || req.resource == ""
    return Response(render(index_tmpl,{"title"=>"Contents","d"=>posts}))
  end
  # Bootstrap css/js files
  if ismatch(r"^/bootstrap/",req.resource)
    f = open(".$(req.resource)")
    h = HttpCommon.headers() #need to set headers correctly
    if ismatch(r".css$",req.resource)
      h["Content-Type"] = "text/css; charset=utf-8"
    elseif ismatch(r".css$",req.resource)
      h["Content-Type"] = "text/js; charset=utf-8"
    end
    return Response(readall(f),h)
  end
  # Posts
  i = 0
  for name in posts["url"]
    i += 1
    if ismatch(Regex("^/$name\$"),req.resource)
      return Response(render(post_tmpl,
        {"title"=>posts["title"][i],
         "date" =>posts["creationdate"][i],
         "summary"=>posts["summary"][i],
         "content"=>posts["contents"][i]}))
    end
  end
  return Response(404)
end

println("about to run")
run(Server(myhandler),8080)


text:
@doc raw"""
    general_form(xi::Point2D, xo::Point2D)

Returns the general form equation of a line that passes through points `xi` and `xo` as a
[`SVector`](@ref) which holds A, B and C such that:
```math
A \cdot x + B \cdot y +  C = 0.
```
"""
function general_form(xi::Point2D, xo::Point2D)
    A = xi[2] - xo[2]
    B = xo[1] - xi[1]
    C = xi[1] * xo[2] - xo[1] * xi[2]
    ABC = SVector(A, B, C)
    ABC /= norm(ABC) # for numerical reasons (?)
    return ABC
end
# the function above is faster than using linear algebra:
# [x1 y1; x2 y2] ⋅ [Â, B̂] = [-1, -1] yields to solution for Â and B̂ such that
# Â ⋅ x + B̂ ⋅ y + 1 = 0. Here we make it faster.
# F = transpose(hcat(xi, xo))
# b = SVector(-1, -1)
# x = F \ b
# ABC = vcat(b, 1) # and then we can normalize

# We could use LazySets.jl for some part of this particular function but it seems slow!
"""
    intersections(mesh::Mesh, cell_id::Int, track::Track)

Computes the entry and exit points of a given `track` that is known to cross element with id
`cell_id`.
"""
function intersections(
    mesh::Mesh, cell_id::Int32, track::Track{BIn,BOut,DFwd,DBwd,T}
) where {BIn,BOut,DFwd,DBwd,T}
    @unpack model, cell_nodes = mesh

    node_coordinates = get_node_coordinates(get_grid(model))
    node_ids = cell_nodes[cell_id]

    int_points = MVector{4,Point2D{T}}(
        Point2D(0, 0), Point2D(0, 0), Point2D(0, 0), Point2D(0, 0)
    )

    n_int = 0
    parallel_found = false
    for i in eachindex(node_ids)

        j = i == lastindex(node_ids) ? 1 : i + 1

        # get node coordinates and cast them to Point2D
        p1 = convert(Point2D{T}, node_coordinates[node_ids[i]])
        p2 = convert(Point2D{T}, node_coordinates[node_ids[j]])

        # compute general form equation for the selected element face
        ABC = general_form(p1, p2)

        # compute intersections between track and element face
        are_parallel, x_int = intersection(track.ABC, ABC)

        # if the track is parallel to the element face, it can be on top of it or do not
        # cross at all. Either way, we just ignore this case because we can compute the
        # intersections using the other faces.
        if are_parallel
            parallel_found = true
            continue

        # if the intersection is outside the face, avoid it
        elseif !point_in_segment(p1, p2, x_int)
            continue

        # this is a valid intersection, store it
        else
            n_int += 1
            int_points[n_int] = x_int
        end
    end

    # now let's check all the possible situations and handle extreme cases
    if n_int in (3, 4)
        ℓ = zero(T)
        # get the points that have the maximun distance between
        for i in 2:n_int, j in i:n_int
            x1 = int_points[i-1]
            x2 = int_points[j]
            ℓi = norm(x1 - x2)
            if ℓi > ℓ
                x_int1 = x1
                x_int2 = x2
                ℓ = ℓi
            end
        end
        return order_intersection_points(track, x_int1, x_int2)

    elseif n_int == 2 && parallel_found

        x_int1 = int_points[1]
        x_int2 = int_points[2]
        return order_intersection_points(track, x_int1, x_int2)

    elseif n_int == 2 && !parallel_found

        x_int1 = int_points[1]
        x_int2 = int_points[2]

        if isapprox(x_int1, x_int2)
            # do nothing and move a tiny step forward in the parent function, we are on a
            # vertex
            return x_int1, x_int2
        else
            return order_intersection_points(track, x_int1, x_int2)
        end
    elseif iszero(n_int) || isone(n_int)
        # the parent function will move a tiny step further
        return Point2D{T}(0, 0), Point2D{T}(0, 0)
        # error("This is an unexpected case. Please, submit an issue.")
    end
end

"""
    intersection(ABC1, ABC2)

Computes the intersection between two lines with general form equations given by `ABC1` and
`ABC2`. It also indicates if the lines are parallel.
"""
function intersection(ABC1::AbstractVector, ABC2::AbstractVector)
    a = ABC1[2] * ABC2[1]
    b = ABC2[2] * ABC1[1]
    x = y = zero(a)
    are_parallel = isapprox(a, b)
    if !are_parallel
        det = a - b
        x = (ABC1[3] * ABC2[2] - ABC2[3] * ABC1[2]) / det
        y = (ABC1[1] * ABC2[3] - ABC2[1] * ABC1[3]) / det
    end
    return are_parallel, Point2D(x, y)
end

# this one returns a matrix even if F is not invertible (have the same speed as the above)
# function intersection(ABC1, ABC2)
#     F = @SMatrix [ABC1[1] ABC1[2]
#                   ABC2[1] ABC2[2]]
#     b = @SVector [ABC1[3], ABC2[3]]

#     x = F \ (-b)

#     return Point2D(x[1], x[2])
# end

function order_intersection_points(track::Track, x1::Point2D, x2::Point2D)
    @unpack ϕ = track

    if isless(ϕ, π/2)
        if x1[1] < x2[1]
            xi = x1
            xo = x2
        else
            xi = x2
            xo = x1
        end
    else
        if x1[1] > x2[1]
            xi = x1
            xo = x2
        else
            xi = x2
            xo = x1
        end
    end
    return xi, xo
end


text:
using GeometryTypes: slice
@testset "Mesh Types" begin

@testset "Merging Mesh" begin
    baselen = 0.4f0
    dirlen = 2f0
    a = [
        (HyperRectangle{3,Float32}(Vec3f0(baselen), Vec3f0(dirlen, baselen, baselen)), RGBA(1f0,0f0,0f0,1f0)),
        (HyperRectangle{3,Float32}(Vec3f0(baselen), Vec3f0(baselen, dirlen, baselen)), RGBA(0f0,1f0,0f0,1f0)),
        (HyperRectangle{3,Float32}(Vec3f0(baselen), Vec3f0(baselen, baselen, dirlen)), RGBA(0f0,0f0,1f0,1f0))
    ]

    am = [GLNormalMesh(x) for x in a]
    axis = merge(am)
    axis2 = merge(am)
    @test axis == axis2 # test ==
    @test typeof(am[1]) == GLNormalColorMesh

    @test typeof(axis) == GLNormalAttributeMesh
    @test typeof(GLPlainMesh(axis)) == GLPlainMesh

    @test vertextype(axis) == Point{3, Float32}
    @test normaltype(axis) == Normal{3, Float32}
    @test facetype(axis) == Face{3, GLIndex}
    @test hasvertices(axis)
    @test hasfaces(axis)
    @test hasnormals(axis)
    @test !hascolors(axis)
end

@testset "Show" begin
    baselen = 0.4f0
    dirlen = 2f0
    m = GLNormalMesh(HyperRectangle{3,Float32}(Vec3f0(baselen),
                     Vec3f0(dirlen, baselen, baselen)))
    io = IOBuffer()
    show(io,m)
    seekstart(io)
    s = "HomogenousMesh(\n    faces: 12xGeometryTypes.Face{3,UInt32,-1},     normals: 24xGeometryTypes.Normal{3,Float32},     vertices: 24xFixedSizeArrays.Point{3,Float32}, )\n"
    #@fact readall(io) --> s #Win32 and Win64 have different ordering it seems.
end

@testset "construction" begin
    VT = vertextype(GLNormalMesh)
    FT = facetype(GLNormalMesh)
    vs = [VT(0., 0, 0), VT(1., 0, 0), VT(0., 1, 0)]
    fs = [FT(1, 2, 3)]

    # test for https://github.com/JuliaGeometry/GeometryTypes.jl/issues/92
    m = HomogenousMesh(vs, fs)
    @test HomogenousMesh(m) == m
end

@testset "Primitives" begin
    m = GLNormalMesh(Sphere(Point3f0(0), 1f0))
    @test length(vertices(m)) == 576
    @test length(faces(m)) == 1058
    @test length(normals(m)) == 576
end

@testset "Slice" begin
    test_mesh = HomogenousMesh(
        Point{3, Float64}[
            (0.0,0.0,10.0),
            (0.0,10.0,10.0),
            (0.0,0.0,0.0),
            (0.0,10.0,0.0),
            (10.0,0.0,10.0),
            (10.0,10.0,10.0),
            (10.0,0.0,0.0),
            (10.0,10.0,0.0)
        ],
        Face{3, Int}[
            (3,7,5),
            (1,3,5),
            (1,2,3),
            (2,4,3),
            (1,5,6),
            (2,1,6),
            (4,8,3),
            (3,8,7),
            (7,8,6),
            (7,6,5),
            (2,6,4),
            (4,6,8)
    ])
    s1 = slice(test_mesh, 1.0)
    s2 = slice(test_mesh, 2.0)
    s3 = slice(test_mesh, 0.0)

    @test length(s1) == 8
    @test length(s1) == length(s2)
    @test length(s3) == 0
    exp1 = [([0.0,1.0],[0.0,0.0]),([0.0,0.0],[1.0,0.0]),
                      ([1.0,0.0],[10.0,0.0]),([10.0,0.0],[10.0,1.0]),
                      ([10.0,1.0],[10.0,10.0]),([10.0,10.0],[1.0,10.0]),
                      ([1.0,10.0],[0.0,10.0]),([0.0,10.0],[0.0,1.0])]
    exp2 = [([0.0,2.0],[0.0,0.0]),([0.0,0.0],[2.0,0.0]),
                      ([2.0,0.0],[10.0,0.0]),([10.0,0.0],[10.0,2.0]),
                      ([10.0,2.0],[10.0,10.0]),([10.0,10.0],[2.0,10.0]),
                      ([2.0,10.0],[0.0,10.0]),([0.0,10.0],[0.0,2.0])]
    @test length(s1) == length(exp1)
    @test length(s2) == length(exp2)
end

@testset "checkbounds" begin
    m1 = HomogenousMesh([Point{3, Float64}(0.0,0.0,10.0),
                         Point{3, Float64}(0.0,10.0,10.0),
                         Point{3, Float64}(0.0,0.0,0.0)],
                        [Face{3, Int}(1,2,3)])
    @test checkbounds(m1)

    m2 = HomogenousMesh([Point{3, Float64}(0.0,0.0,10.0),
                         Point{3, Float64}(0.0,10.0,10.0),
                         Point{3, Float64}(0.0,0.0,0.0)],
                        [Face{3, GLIndex}(5,1,2)])
    @test !checkbounds(m2)

    # empty case
    m3 = HomogenousMesh([Point{3, Float64}(0.0,0.0,10.0),
                         Point{3, Float64}(0.0,10.0,10.0),
                         Point{3, Float64}(0.0,0.0,0.0)],
                        Face{3, GLIndex}[])
    @test checkbounds(m3)
end

@testset "vertex normals" begin
    test_mesh = HomogenousMesh(
    Point{3,Float64}[
        (0.0,0.0,10.0),
        (0.0,10.0,10.0),
        (0.0,0.0,0.0),
        (0.0,10.0,0.0),
        (10.0,0.0,10.0),
        (10.0,10.0,10.0),
        (10.0,0.0,0.0),
        (10.0,10.0,0.0),
    ], Face{3, Int}[
        (3,7,5),
        (1,3,5),
        (1,2,3),
        (3,2,4),
        (1,5,6),
        (2,1,6),
        (4,8,3),
        (3,8,7),
        (7,8,6),
        (5,7,6),
        (2,6,4),
        (4,6,8)
    ])
    ns = normals(test_mesh.vertices, test_mesh.faces, Normal{3, Float32})
    @test length(ns) == length(test_mesh.vertices)
    expect = Normal{3, Float32}[
        (-0.408248,-0.408248,0.816497),
        (-0.816497,0.408248,0.408248),
        (-0.57735,-0.57735,-0.57735),
        (-0.408248,0.816497,-0.408248),
        (0.408248,-0.816497,0.408248),
        (0.57735,0.57735,0.57735),
        (0.816497,-0.408248,-0.408248),
        (0.408248,0.408248,-0.816497)
    ]
    @test all(isapprox.(ns, expect))
end

@testset "conversion" begin
    VT = vertextype(GLNormalMesh)
    FT = facetype(GLNormalMesh)
    vs = [VT(0., 0, 0), VT(1., 0, 0), VT(0., 1, 0)]
    fs = [FT(1, 2, 3)]
    mesh = PlainMesh{eltype(VT), FT}(vertices=vs, faces=fs)
    @test convert(GLNormalMesh, mesh) == GLNormalMesh(vs, fs)
end
@testset "construction" begin
    VT = vertextype(GLNormalMesh)
    FT = facetype(GLNormalMesh)
    vs = [VT(0., 0, 0), VT(1., 0, 0), VT(0., 1, 0)]
    fs = [FT(1, 2, 3)]

    # test for https://github.com/JuliaGeometry/GeometryTypes.jl/issues/92
    m = HomogenousMesh(vs, fs)
    @test HomogenousMesh(m) == m
end

@testset "unique vertices" begin
    # geom <> with overlapping verts for middle points
    v = Point3f0[(0, 0, 0), (1, 1, 0), (1, -1, 0), (1, 1, 0), (2, 0, 0), (1, -1, 0)]
    f = GLTriangle[(1, 2, 3), (4, 5, 6)]
    m = GLPlainMesh(v, f)
    GeometryTypes.remove_overlap!(m)
    @test vertices(m) == Point3f0[(0, 0, 0), (1, 1, 0), (1, -1, 0), (2, 0, 0)]
    @test faces(m) == GLTriangle[(1, 2, 3), (2, 4, 3)]
end

end


text:
using Flux
using Flux: onehotbatch, onecold, flatten
using Flux.Losses: logitcrossentropy
using Flux.Data: DataLoader
using Parameters: @with_kw
using Statistics: mean
using CUDA
using MLDatasets: CIFAR10
using MLDataPattern: splitobs

if CUDA.has_cuda()
    @info "CUDA is on"
    CUDA.allowscalar(false)
end

function get_processed_data(args)
    x, y = CIFAR10.traindata()

    (train_x, train_y), (val_x, val_y) = splitobs((x, y), at=1-args.valsplit)

    train_x = float(train_x)
    train_y = onehotbatch(train_y, 0:9)
    val_x = float(val_x)
    val_y = onehotbatch(val_y, 0:9)

    return (train_x, train_y), (val_x, val_y)
end

function get_test_data()
    test_x, test_y = CIFAR10.testdata()

    test_x = float(test_x)
    test_y = onehotbatch(test_y, 0:9)

    return test_x, test_y
end

# VGG16 and VGG19 models
function vgg16()
    Chain(
        Conv((3, 3), 3 => 64, relu, pad=(1, 1), stride=(1, 1)),
        BatchNorm(64),
        Conv((3, 3), 64 => 64, relu, pad=(1, 1), stride=(1, 1)),
        BatchNorm(64),
        MaxPool((2,2)),
        Conv((3, 3), 64 => 128, relu, pad=(1, 1), stride=(1, 1)),
        BatchNorm(128),
        Conv((3, 3), 128 => 128, relu, pad=(1, 1), stride=(1, 1)),
        BatchNorm(128),
        MaxPool((2,2)),
        Conv((3, 3), 128 => 256, relu, pad=(1, 1), stride=(1, 1)),
        BatchNorm(256),
        Conv((3, 3), 256 => 256, relu, pad=(1, 1), stride=(1, 1)),
        BatchNorm(256),
        Conv((3, 3), 256 => 256, relu, pad=(1, 1), stride=(1, 1)),
        BatchNorm(256),
        MaxPool((2,2)),
        Conv((3, 3), 256 => 512, relu, pad=(1, 1), stride=(1, 1)),
        BatchNorm(512),
        Conv((3, 3), 512 => 512, relu, pad=(1, 1), stride=(1, 1)),
        BatchNorm(512),
        Conv((3, 3), 512 => 512, relu, pad=(1, 1), stride=(1, 1)),
        BatchNorm(512),
        MaxPool((2,2)),
        Conv((3, 3), 512 => 512, relu, pad=(1, 1), stride=(1, 1)),
        BatchNorm(512),
        Conv((3, 3), 512 => 512, relu, pad=(1, 1), stride=(1, 1)),
        BatchNorm(512),
        Conv((3, 3), 512 => 512, relu, pad=(1, 1), stride=(1, 1)),
        BatchNorm(512),
        MaxPool((2,2)),
        flatten,
        Dense(512, 4096, relu),
        Dropout(0.5),
        Dense(4096, 4096, relu),
        Dropout(0.5),
        Dense(4096, 10)
    )
end

@with_kw mutable struct Args
    batchsize::Int = 128
    lr::Float64 = 3e-4
    epochs::Int = 50
    valsplit::Float64 = 0.1
end

function train(; kws...)
    # Initialize the hyperparameters
    args = Args(; kws...)

    # Load the train, validation data
    train_data, val_data = get_processed_data(args)

    train_loader = DataLoader(train_data, batchsize=args.batchsize, shuffle=true)
    val_loader = DataLoader(val_data, batchsize=args.batchsize)

    @info("Constructing Model")
    m = vgg16() |> gpu

    loss(x, y) = logitcrossentropy(m(x), y)

    ## Training
    # Defining the optimizer
    opt = ADAM(args.lr)
    ps = Flux.params(m)

    @info("Training....")
    # Starting to train models
    for epoch in 1:args.epochs
        @info "Epoch $epoch"

        for (x, y) in train_loader
            x, y = x |> gpu, y |> gpu
            gs = Flux.gradient(() -> loss(x,y), ps)
            Flux.update!(opt, ps, gs)
        end

        validation_loss = 0f0
        for (x, y) in val_loader
            x, y = x |> gpu, y |> gpu
            validation_loss += loss(x, y)
        end
        validation_loss /= length(val_loader)
        @show validation_loss
    end

    return m
end

function test(m; kws...)
    args = Args(kws...)

    test_data = get_test_data()
    test_loader = DataLoader(test_data, batchsize=args.batchsize)

    correct, total = 0, 0
    for (x, y) in test_loader
        x, y = x |> gpu, y |> gpu
        correct += sum(onecold(cpu(m(x))) .== onecold(cpu(y)))
        total += size(y, 2)
    end
    test_accuracy = correct / total

    # Print the final accuracy
    @show test_accuracy
end

if abspath(PROGRAM_FILE) == @__FILE__
    m = train()
    test(m)
end

max_stars_repo_name: Ossifragus/model-zoo

max_issues_repo_name: Ossifragus/model-zoo

max_forks_repo_name: Ossifragus/model-zoo


text:
@inline function ArbFloat(x)
    prec = DEFAULT_PRECISION.x
    res  = ArbFloat{prec}(x)
    return res
end
@inline function ArbReal(x)
    prec = DEFAULT_PRECISION.x
    res  = ArbReal{prec}(x)
    return res
end

@inline function ArbComplex(x)
    prec = DEFAULT_PRECISION.x
    res  = ArbComplex{prec}(x)
    return res
end

@inline function ArbComplex(x, y)
    prec = DEFAULT_PRECISION.x
    res  = ArbComplex{prec}(x, y)
    return res
end

@inline function ArbComplex(x::I, y::T) where {I<:Integer, T<:AbstractFloat}
    prec = DEFAULT_PRECISION.x
    res  = ArbComplex{prec}(T(x), y)
    return res
end

@inline function ArbComplex(x::T, y::I) where {I<:Integer, T<:AbstractFloat}
    prec = DEFAULT_PRECISION.x
    res  = ArbComplex{prec}(x, T(y))
    return res
end

# IEEEFloat

# rounds up (widens)
function Float64(x::Mag)
    z = ccall(@libarb(mag_get_d), Float64, (Ref{Mag},), x)
    return z
end

function Float64(x::ArbFloat{P}) where {P}
    rounding = match_rounding_mode(RoundNearest)
    z = ccall(@libarb(arf_get_d), Float64, (Ref{ArbFloat}, Cint), x, rounding)
    return z
end
Float32(x::ArbFloat{P}) where {P} = Float32(Float64(x))
Float16(x::ArbFloat{P}) where {P} = Float16(Float64(x))

function Float64(x::ArbReal{P}) where {P}
    x = midpoint(x)
    y = ArbFloat{P}(x)
    Float64(y)
end
Float32(x::ArbReal{P}) where {P} = Float32(Float64(x))
Float16(x::ArbReal{P}) where {P} = Float16(Float64(x))

function Float64(x::ArbComplex{P}) where {P}
    x = midpoint(real(x))
    y = ArbFloat{P}(x)
    Float64(y)
end
Float32(x::ArbComplex{P}) where {P} = Float32(Float64(x))
Float16(x::ArbComplex{P}) where {P} = Float16(Float64(x))

function Float64(x::ArbFloat{P}, roundingmode::RoundingMode) where {P}
    rounding = match_rounding_mode(roundingmode)
    z = ccall(@libarb(arf_get_d), Float64, (Ref{ArbFloat}, Cint), x, rounding)
    return z
end
Float32(x::ArbFloat{P}, roundingmode::RoundingMode) where {P} = Float32(Float64(x, roundingmode))
Float16(x::ArbFloat{P}, roundingmode::RoundingMode) where {P} = Float16(Float64(x, roundingmode))

function Float64(x::ArbReal{P}, roundingmode::RoundingMode) where {P}
    x = midpoint(x)
    y = ArbFloat{P}(x)
    Float64(y, roundingmode)
end
Float32(x::ArbReal{P}, roundingmode::RoundingMode) where {P} = Float32(Float64(x, roundingmode))
Float16(x::ArbReal{P}, roundingmode::RoundingMode) where {P} = Float16(Float64(x, roundingmode))

function Float64(x::ArbComplex{P}, roundingmode::RoundingMode) where {P}
    x = midpoint(real(x))
    y = ArbFloat{P}(x)
    Float64(y, roundingmode)
end
Float32(x::ArbComplex{P}, roundingmode::RoundingMode) where {P} = Float32(Float64(x, roundingmode))
Float16(x::ArbComplex{P}, roundingmode::RoundingMode) where {P} = Float16(Float64(x, roundingmode))

function Complex{T}(x::ArbComplex{P}, roundingmode::RoundingMode) where {P, T<:IEEEFloat}
   re = T(real(x), roundingmode)
   im = T(imag(x), roundingmode)
   return Complex{T}(re, im)
end

Complex{T}(x::ArbComplex{P}) where {P, T<:IEEEFloat} = Complex{T}(x, RoundNearest)
Complex(x::ArbComplex{P}) where {P} = Complex{Float64}(x, RoundNearest)

# integers

for I in (:Int8, :Int16, :Int32, :Int64, :Int128)
  @eval begin

    function $I(x::ArbFloat{P}) where {P}
        !isinteger(x) && throw(InexactError("$(x) is not an integer"))
        bi = BigInt(BigFloat(x))
        !(typemin($I) <= bi <= typemax($I)) && throw(InexactError("$(x)"))
        return $I(bi)
    end

    function $I(x::ArbReal{P}) where {P}
        (!isexact(x) | !isinteger(x)) && throw(InexactError("$(x) is not an integer"))
        bi = BigInt(BigFloat(x))
        !(typemin($I) <= bi <= typemax($I)) && throw(InexactError("$(x)"))
        return $I(bi)
    end

    function $I(x::ArbComplex{P}) where {P}
        (!isexact(x) | !isinteger(x) | !iszero(imag(x))) && throw(InexactError("$(x) is not an integer"))
        bi = BigInt(BigFloat(x))
        !(typemin($I) <= bi <= typemax($I)) && throw(InexactError("$(x)"))
        return $I(bi)
    end

  end
end

# Rational

ArbFloat(x::T) where {S, T<:Rational{S}} = ArbFloat(x.num)/ArbFloat(x.den)
ArbReal(x::T) where {S, T<:Rational{S}} = ArbReal(x.num)/ArbReal(x.den)
ArbComplex(x::T) where {S, T<:Rational{S}} = ArbComplex(ArbReal(x), ArbReal(0))
ArbComplex(x::T, y::T) where {S, T<:Rational{S}} = ArbComplex(ArbReal(x), ArbReal(y))

# Irrational

function ArbReal{P}(x::Irrational{S}) where {P,S}
   mid = ArbFloat{P}(x)
   rad = ulp(mid)
   return setball(mid, rad)
end

function ArbReal(x::Irrational{S}) where {S}
   P = workingprecision(ArbReal)
   mid = ArbFloat{P}(x)
   rad = ulp(mid)
   return setball(mid, rad)
end

ArbComplex(x::Irrational{S}) where {S} = ArbComplex(ArbReal(x), ArbReal(0))
ArbComplex{P}(x::Irrational{S}) where {P,S} = ArbComplex{P}(ArbReal{P}(x), ArbReal{P}(0))
ArbComplex(x::Irrational{S}, y::Real) where {S} = ArbComplex(ArbReal(x), ArbReal(0))
ArbComplex{P}(x::Irrational{S}, y::Real) where {P,S} = ArbComplex{P}(ArbReal{P}(x), ArbReal{P}(0))

# fallback

ArbFloat{P}(x::T) where {P,T<:Real} = ArbFloat{P}(BigFloat(x))
ArbFloat(x::T) where {T<:Real} = ArbFloat{workingprecision(ArbFloat)}(BigFloat(x))
ArbFloat{P}(x::T) where {P,T<:Complex} = ArbFloat{P}(BigFloat(real(x)))
ArbFloat(x::T) where {T<:Complex} = ArbFloat{workingprecision(ArbFloat)}(BigFloat(real(x)))

ArbReal{P}(x::T) where {P,T<:Real} = ArbReal{P}(BigFloat(x))
ArbReal(x::T) where {T<:Real} = ArbReal{workingprecision(ArbReal)}(BigFloat(real(x)))
ArbReal{P}(x::T) where {P,T<:Complex} = ArbReal{P}(BigFloat(x))
ArbReal(x::T) where {T<:Complex} = ArbReal{workingprecision(ArbReal)}(BigFloat(real(x)))

ArbComplex{P}(x::BigInt) where {P} = ArbComplex{P}(ArbReal{P}(x))
ArbComplex{P}(x::BigFloat) where {P} = ArbComplex{P}(ArbReal{P}(x))
ArbComplex{P}(x::T) where {P,T<:Real} = ArbComplex{P}(BigFloat(x))
#ArbComplex(x::T) where {T<:Real} = ArbComplex{workingprecision(ArbComplex)}(BigFloat(x))
ArbComplex{P}(x::T, y::T) where {P,T<:Real} = ArbComplex{P}(BigFloat(x), BigFloat(y))
#ArbComplex(x::T, y::T) where {T<:Real} = ArbComplex{workingprecision(ArbComplex)}(BigFloat(x), BigFloat(y))
ArbComplex{P}(x::T) where {P,T<:Complex} = ArbComplex{P}(BigFloat(real(x)), BigFloat(imag(x)))
ArbComplex(x::T) where {T<:Complex} = ArbComplex{workingprecision(ArbComplex)}(BigFloat(real(x)), BigFloat(imag(x)))

# retype

ArbFloat(x::ArbReal{P}) where {P} = ArbFloat{P}(x)
ArbFloat(x::ArbComplex{P}) where {P} = ArbFloat{P}(real(x))
ArbReal(x::ArbFloat{P}) where {P} = ArbReal{P}(x)
ArbReal(x::ArbComplex{P}) where {P} = ArbReal{P}(real(x))
# ArbComplex(x::ArbFloat{P}) where {P} = ArbComplex{P}(ArbReal{P}(x))
# ArbComplex(x::ArbReal{P}) where {P} = ArbComplex{P}(x)

ArbFloat{Q}(x::ArbReal{P}) where {P,Q} = ArbFloat{Q}(ArbReal{Q}(x))
ArbFloat{Q}(x::ArbComplex{P}) where {P,Q} = ArbFloat{Q}(ArbReal{Q}(real(x)))
ArbReal{Q}(x::ArbFloat{P}) where {P,Q} = ArbReal{Q}(ArbFloat{Q}(x))
ArbReal{Q}(x::ArbComplex{P}) where {P,Q} = ArbReal{Q}(real(x))
ArbComplex{Q}(x::ArbFloat{P}) where {P,Q} = ArbComplex{Q}(ArbReal{Q}(ArbFloat{Q}(x)))
ArbComplex{Q}(x::ArbReal{P}) where {P,Q} = ArbComplex{Q}(ArbReal{Q}(x))

@inline convert(::Type{ArbFloat{Q}}, x::ArbReal{P}) where {P,Q} = ArbFloat{Q}(x)
@inline convert(::Type{ArbFloat{Q}}, x::ArbComplex{P}) where {P,Q} = ArbFloat{Q}(x)
@inline convert(::Type{ArbReal{Q}}, x::ArbFloat{P}) where {P,Q} = ArbReal{Q}(x)
@inline convert(::Type{ArbReal{Q}}, x::ArbComplex{P}) where {P,Q} = ArbReal{Q}(x)
@inline convert(::Type{ArbComplex{Q}}, x::ArbFloat{P}) where {P,Q} = ArbComplex{Q}(x)
@inline convert(::Type{ArbComplex{Q}}, x::ArbReal{P}) where {P,Q} = ArbComplex{Q}(x)

# change precision

function ArbFloat{P}(x::ArbFloat{Q}, roundingmode::RoundingMode) where {P,Q}
    rounding = match_rounding_mode(roundingmode)
    z = ArbFloat{P}()
    rnd = ccall(@libarb(arf_set_round), Cint,
                (Ref{ArbFloat}, Ref{ArbFloat}, Clong, Cint), z, x, P, rounding)
    return z
end
ArbFloat{P}(x::ArbFloat{Q}) where {P,Q} = ArbFloat{P}(x, RoundNearest)

function ArbReal{P}(x::ArbReal{Q}) where {P,Q}
    z = ArbReal{P}()
    ccall(@libarb(arb_set_round), Cvoid,
          (Ref{ArbReal}, Ref{ArbReal}, Clong), z, x, P)
    return z
end

function ArbComplex{P}(x::ArbComplex{Q}) where {P,Q}
    z = ArbComplex{P}()
    ccall(@libarb(acb_set_round), Cvoid,
          (Ref{ArbComplex}, Ref{ArbComplex}, Clong), z, x, P)
    return z
end

# widen

widen(::Type{ArbFloat{P}}) where {P} = ArbFloat{P+P+4}
widen(::Type{ArbReal{P}}) where {P} = ArbReal{P+P+4}
widen(::Type{ArbComplex{P}}) where {P} = ArbComplex{P+P+4}


text:
struct FileLogger <: AbstractLogger
    logger::SimpleLogger
    always_flush::Bool
end

function FileLogger(path; append=false, always_flush=true)
    filehandle = open(path, append ? "a" : "w")
    FileLogger(SimpleLogger(filehandle, BelowMinLevel), always_flush)
end

function handle_message(filelogger::FileLogger, args...; kwargs...)
    handle_message(filelogger.logger, args...; kwargs...)
    filelogger.always_flush && flush(filelogger.logger.stream)
end
shouldlog(filelogger::FileLogger, arg...) = true
min_enabled_level(filelogger::FileLogger) = BelowMinLevel
catch_exceptions(filelogger::FileLogger) = catch_exceptions(filelogger.logger)


text:
using AWS
using AWS.AWSServices: s3
using AWSS3
using Test
using Dates
using Retry
using HTTP
using FilePathsBase
using FilePathsBase: /, join
using FilePathsBase.TestPaths
using UUIDs: uuid4
using JSON3

aws = AWSConfig()

@testset "AWSS3.jl" begin
    include("s3path.jl")
    include("awss3.jl")
end


text:
using Term

prt(pan) = begin
    print(" " * hLine(pan.measure.w; style="red"))
    print(vLine(pan.measure.h; style="red") * pan)
    println(pan.measure, "   ", length(chars(pan.segments[1].plain)), "  ", length(pan.segments) )
end

# check that unreasonable widhts are ignored
prt(TextBox(
    "nofit"^25;
    width=1000
))

prt(TextBox(
    "truncate"^25;
    width=100,
    fit=:truncate
))

prt(TextBox(
    "truncate"^25;
    width=100,
))

prt(TextBox(
    "truncate"^8;
    fit=:fit
))

prt(TextBox(
    "[red]truncate[/red]"^8;
    fit=:fit
))

prt(TextBox(
    "[red]tru\nncate[/red]test"^1;
    fit=:fit
))

prt(
    TextBox(
    "[red]truncate[/red]test"^8;
    fit=:fit,
    justify=:left
)
)


text:
export SoftSphere

"""
    SoftSphere(; cutoff, nl_only, lorentz_mixing, force_units, energy_units, skip_shortcut)

The soft-sphere potential.
"""
struct SoftSphere{S, C, F, E} <: PairwiseInteraction
    cutoff::C
    nl_only::Bool
    lorentz_mixing::Bool
    force_units::F
    energy_units::E
end

function SoftSphere(;
                    cutoff=NoCutoff(),
                    nl_only=false,
                    lorentz_mixing=true,
                    force_units=u"kJ * mol^-1 * nm^-1",
                    energy_units=u"kJ * mol^-1",
                    skip_shortcut=false)
    return SoftSphere{skip_shortcut, typeof(cutoff), typeof(force_units), typeof(energy_units)}(
        cutoff, nl_only, lorentz_mixing, force_units, energy_units)
end

@inline @inbounds function force(inter::SoftSphere{S, C},
                                    dr,
                                    coord_i,
                                    coord_j,
                                    atom_i,
                                    atom_j,
                                    box_size) where {S, C}
    r2 = sum(abs2, dr)

    if !S && (iszero(atom_i.ϵ) || iszero(atom_j.ϵ) || iszero(atom_i.σ) || iszero(atom_j.σ))
        return ustrip.(zero(coord_i)) * inter.force_units
    end

    # Lorentz-Berthelot mixing rules use the arithmetic average for σ
    # Otherwise use the geometric average
    σ = inter.lorentz_mixing ? (atom_i.σ + atom_j.σ) / 2 : sqrt(atom_i.σ * atom_j.σ)
    ϵ = sqrt(atom_i.ϵ * atom_j.ϵ)

    cutoff = inter.cutoff
    σ2 = σ^2
    params = (σ2, ϵ)

    if cutoff_points(C) == 0
        f = force_divr_nocutoff(inter, r2, inv(r2), params)
    elseif cutoff_points(C) == 1
        r2 > cutoff.sqdist_cutoff && return ustrip.(zero(coord_i)) * inter.force_units

        f = force_divr_cutoff(cutoff, r2, inter, params)
    elseif cutoff_points(C) == 2
        r2 > cutoff.sqdist_cutoff && return ustrip.(zero(coord_i)) * inter.force_units

        if r2 < cutoff.sqdist_activation
            f = force_divr_nocutoff(inter, r2, inv(r2), params)
        else
            f = force_divr_cutoff(cutoff, r2, inter, params)
        end
    end

    return f * dr
end

@fastmath function force_divr_nocutoff(::SoftSphere, r2, invr2, (σ2, ϵ))
    six_term = (σ2 * invr2) ^ 3

    return (24ϵ * invr2) * 2 * six_term ^ 2
end

@inbounds function potential_energy(inter::SoftSphere{S, C},
                                    dr,
                                    coord_i,
                                    coord_j,
                                    atom_i,
                                    atom_j,
                                    box_size) where {S, C}
    r2 = sum(abs2, dr)

    if !S && (iszero(atom_i.ϵ) || iszero(atom_j.ϵ) || iszero(atom_i.σ) || iszero(atom_j.σ))
        return ustrip(zero(box_size[1])) * inter.energy_units
    end

    σ = inter.lorentz_mixing ? (atom_i.σ + atom_j.σ) / 2 : sqrt(atom_i.σ * atom_j.σ)
    ϵ = sqrt(atom_i.ϵ * atom_j.ϵ)

    cutoff = inter.cutoff
    σ2 = σ^2
    params = (σ2, ϵ)

    if cutoff_points(C) == 0
        potential(inter, r2, inv(r2), params)
    elseif cutoff_points(C) == 1
        r2 > cutoff.sqdist_cutoff && return ustrip(zero(box_size[1])) * inter.energy_units

        potential_cutoff(cutoff, r2, inter, params)
    elseif cutoff_points(C) == 2
        r2 > cutoff.sqdist_cutoff && return ustrip(zero(box_size[1])) * inter.energy_units

        if r2 < cutoff.sqdist_activation
            potential(inter, r2, inv(r2), params)
        else
            potential_cutoff(cutoff, r2, inter, params)
        end
    end
end

@fastmath function potential(::SoftSphere, r2, invr2, (σ2, ϵ))
    six_term = (σ2 * invr2) ^ 3

    return 4ϵ * (six_term ^ 2)
end


text:
# ---
# title: 1038. Binary Search Tree to Greater Sum Tree
# id: problem1038
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Binary Search Tree
# link: <https://leetcode.com/problems/binary-search-tree-to-greater-sum-tree/description/>
# hidden: true
# ---
#
# Given the `root` of a Binary Search Tree (BST), convert it to a Greater Tree
# such that every key of the original BST is changed to the original key plus
# sum of all keys greater than the original key in BST.
#
# As a reminder, a _binary search tree_ is a tree that satisfies these
# constraints:
#
#   * The left subtree of a node contains only nodes with keys  **less than**  the node's key.
#   * The right subtree of a node contains only nodes with keys  **greater than**  the node's key.
#   * Both the left and right subtrees must also be binary search trees.
#
# **Note:** This question is the same as 538:
# <https://leetcode.comhttps://leetcode.com/problems/convert-bst-to-greater-tree/>
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2019/05/02/tree.png)
#
#
#
#     Input: root = [4,1,6,0,2,5,7,null,null,null,3,null,null,null,8]
#     Output: [30,36,21,36,35,26,15,null,null,null,33,null,null,null,8]
#
#
# **Example 2:**
#
#
#
#     Input: root = [0,null,1]
#     Output: [1,null,1]
#
#
# **Example 3:**
#
#
#
#     Input: root = [1,0,2]
#     Output: [3,3,2]
#
#
# **Example 4:**
#
#
#
#     Input: root = [3,2,4,1]
#     Output: [7,9,4,10]
#
#
#
#
# **Constraints:**
#
#   * The number of nodes in the tree is in the range `[1, 100]`.
#   * `0 <= Node.val <= 100`
#   * All the values in the tree are **unique**.
#   * `root` is guaranteed to be a valid binary search tree.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
#!/usr/bin/env julia
using JSON
using Referencer

ref = Referencer.reference(Base)
open("test.json", "w") do io
    JSON.print(io, ref)
end


text:
const conf = Config(
	observation_shape = (3,3,3),
	action_space = collect(1:9),
	players = collect(1:2),
	stacked_observations = 1,
	num_workers = 2,
	max_moves = 9,
	num_unroll_steps = 5,
	td_steps = 5,
	PER = false,
	opponent= "human",
	training_steps = 10000,
	batch_size = 32,
	num_iters = 10
	)

const hyper = FeedForwardHP(
	width_hidden = 64,
	depth_representation = 3,
	depth_prediction = 3,
	depth_dynamics = 3,
	depth_policy = 1,
	depth_value = 1,
	depth_reward = 1,
	depth_state_head = 3,
	hidden_state_size = 27,
	reward_activation = tanh
)


text:
using Base.Broadcast: broadcast_shape

@timedtestset "zygote_rules" begin
    @timedtestset "Cholesky (getproperty)" begin
        rng, N = MersenneTwister(123456), 5
        A = randn(rng, N, N)
        S = A' * A + 1e-6I
        C = cholesky(S)

        # Check that non-differentiable ops run forwards and have `nothing` gradients.
        _, back = Zygote.pullback(C->C.info, C)
        @test back(1)[1] == (uplo=nothing, info=nothing, factors=nothing)
        _, back = Zygote.pullback(C->C.uplo, C)
        @test back(1)[1] == (uplo=nothing, info=nothing, factors=nothing)

        # Unit test retrieving the factors.
        @test_throws ErrorException Zygote.pullback(C->C.factors, C)

        # Test getproperty.
        adjoint_test(A->Cholesky(A, :U, 0).U, randn(rng, N, N), A)
        adjoint_test(A->Cholesky(A, :U, 0).L, randn(rng, N, N), A)
        adjoint_test(A->Cholesky(A, :L, 0).U, randn(rng, N, N), A)
        adjoint_test(A->Cholesky(A, :L, 0).L, randn(rng, N, N), A)
    end
    @timedtestset "colwise(::Euclidean, X, Y; dims=2)" begin
        rng, D, P = MersenneTwister(123456), 2, 3
        X, Y, D̄ = randn(rng, D, P), randn(rng, D, P), randn(rng, P)
        adjoint_test((X, Y)->colwise(Euclidean(), X, Y), D̄, X, Y)
    end
    @timedtestset "pairwise(::Euclidean, X, Y; dims=2)" begin
        rng, D, P, Q = MersenneTwister(123456), 2, 3, 5
        X, Y, D̄ = randn(rng, D, P), randn(rng, D, Q), randn(rng, P, Q)
        adjoint_test(
            (X, Y)->pairwise(Euclidean(), X, Y; dims=2), D̄, X, Y;
            rtol=1e-6, atol=1e-6,
        )
    end
    @timedtestset "pairwise(::Euclidean, X; dims=2)" begin
        rng, D, P = MersenneTwister(123456), 2, 3
        X, D̄ = randn(rng, D, P), randn(rng, P, P)
        adjoint_test(X->pairwise(Euclidean(), X; dims=2), D̄, X)
    end
    @timedtestset "Diagonal" begin
        rng, N = MersenneTwister(123456), 11
        adjoint_test(Diagonal, rand(rng, N, N), randn(rng, N))
        adjoint_test(x->Diagonal(x).diag, randn(rng, N), randn(rng, N))
    end
    @timedtestset "broadcast" begin
        @timedtestset "exp" begin
            rng, N = MersenneTwister(123456), 11
            adjoint_test(x->exp.(x), randn(rng, N), randn(rng, N))
        end
        @timedtestset "-" begin
            rng, N = MersenneTwister(123456), 11
            adjoint_test(x->.-x, randn(rng, N), randn(rng, N))
        end
    end
    @timedtestset "Pairwise when X ≈ Y" begin
        rng, D, P = MersenneTwister(13), 2, 3, 5
        X, D̄ = randn(rng, D, P), randn(rng, P, P)
        Y = X .+ 1e-3
        adjoint_test(
            (X, Y)->pairwise(Euclidean(Stheno.dtol), X, Y; dims=2), D̄, X, Y;
            rtol=1e-3, atol=1e-3, fdm=FiniteDifferences.Forward(2, 1)
        ) # relaxed test because of machine precision concerns with finite differences
    end
    @timedtestset "ldiv(::Diagonal, ::Matrix)" begin
        rng, P, Q = MersenneTwister(123456), 13, 15
        Ȳ = randn(rng, P, Q)
        d = randn(rng, P)
        X = randn(rng, P, Q)
        adjoint_test((d, X)->Diagonal(d) \ X, Ȳ, d, X)
    end
end


text:
mutable struct CSVExporter{T <: AbstractSchema} <: AbstractExporter
    isenabled::Bool
    schema::T
end

function isenabled(exporter::CSVExporter)
    return exporter.isenabled
end

function enable!(exporter::CSVExporter)
    exporter.isenabled = true
end

function disable!(exporter::CSVExporter)
    exporter.isenabled  = false
end

struct ExportContext{T₁ <: AbstractDialogControl,   T₂ <: AbstractDialogModel,  T₃ <: AbstractDisplayProperties, T₄ <: Union{AbstractExporter}} <: AbstractContext
    control::T₁
    model::T₂
    display_properties::T₃
    action::T₄
end

function (store::CSVExporter{<:AbstractSchema})(path::Path, li::LabelledIntervals)
    # Convert labelled intervals to a dataframe.
    labelled_intervals = get_labelled_intervals(li)
    df = DataFrame(label = String[], start = Int64[], stop = Int64[], first = Float64[], step = Float64[], last  = Float64[])
    for (key, labelled_interval) in pairs(labelled_intervals)
        nested_interval = labelled_interval.nested_interval
        start = get_start(nested_interval)
        stop = get_stop(nested_interval)
        interval = get_interval(nested_interval)
        push!(df, [key start stop first(interval) step(interval) last(interval)])
    end
    Base.display(df)
    store(path, df)
end

function (store::CSVExporter{<:AbstractSchema})(path::Path, data::DataFrame)
    disable!(store)
    path₁ = joinpath(get_directory(path), get_filename(path))
    # We are trying to create a new file.
    if !is_queryable_file(path₁)
        try
            CSV.write(path₁, data)
        catch e
            println(e)
            CImGui.OpenPopup("Do you have permission to create the file?")
        end
    # We are modifying an existing file.
    elseif is_writeable_file(path₁)
        try
            CSV.write(path₁, data)
        catch e
            println(e)
            CImGui.OpenPopup("Do you have permission to modify the file?")
        end
    else
        CImGui.OpenPopup("Do you have permission to modify the file?")
    end

    if CImGui.BeginPopupModal("Do you have permission to create the file?", C_NULL, CImGui.ImGuiWindowFlags_AlwaysAutoResize)
        CImGui.Text("Unable to write to the specified file.\nPlease verify that you have permission to create the file.\n\n")
        CImGui.Separator()
        CImGui.Button("OK", (120, 0)) && CImGui.CloseCurrentPopup()
        CImGui.SetItemDefaultFocus()
        CImGui.EndPopup()
    end

    if CImGui.BeginPopupModal("Do you have permission to modify the file?", C_NULL, CImGui.ImGuiWindowFlags_AlwaysAutoResize)
        CImGui.Text("Unable to write to the specified file.\nPlease verify that: \n   (1) the file exists; \n   (2) you have permission to modify the file.\n\n")
        CImGui.Separator()
        CImGui.Button("OK", (120, 0)) && CImGui.CloseCurrentPopup()
        CImGui.SetItemDefaultFocus()
        CImGui.EndPopup()
    end

    if CImGui.BeginPopupModal("Do you have permission to read the file?", C_NULL, CImGui.ImGuiWindowFlags_AlwaysAutoResize)
        CImGui.Text("Unable to access the specified file.\nPlease verify that: \n   (1) the file exists; \n   (2) you have permission to read the file.\n\n")
        CImGui.Separator()
        CImGui.Button("OK", (120, 0)) && CImGui.CloseCurrentPopup()
        CImGui.SetItemDefaultFocus()
        CImGui.EndPopup()
    end
end

function (context::ExportContext)(data::AbstractModel...)
        control = context.control
        model = context.model
        display_properties = context.display_properties
        action = context.action
        isenabled(control) ? control(model, display_properties, action) : nothing
        isenabled(action) ? action(get_path(model, ConfirmedStatus()), data...) : nothing
end

function isrunning(context::ExportContext)
        isenabled(context.control) || isenabled(context.action)
end

function enable!(context::ExportContext)
    enable!(context.control)
end


text:
using Gen
using JSON

# get map from character to integer in 1..26
alphabet = collect('a':'z')
push!(alphabet, ' ')
letter_to_int = Dict{Char,Int}()
for (i, letter) in enumerate(alphabet)
    letter_to_int[letter] = i
end

# load bigram data into a transition matrix
# bigrams.json from https://gist.github.com/lydell/c439049abac2c9226e53
data = JSON.parsefile("$(@__DIR__)/bigrams.json")
counts = zeros(Int64, (27, 27))
for (bigram, count) in data
    first = letter_to_int[bigram[1]]
    second = letter_to_int[bigram[2]]
    counts[second, first] = count
end
counts = counts .+ 1 # remove zeros
start_counts = [counts[i,27] for i=1:27]
const probs = start_counts / sum(start_counts)
const transition = counts ./ sum(counts, dims=1)

@gen function generate_original_text(alpha::Float64, len::Int)
    my_probs = probs .* (1 - alpha) .+ fill(1 / 27., 27) .* alpha
    my_transition = transition * (1 - alpha) .+ fill(1 / 27., 27, 27) .* alpha
    local cur::Int
    cur = @trace(categorical(my_probs), 1)
    for l=2:len
        cur = @trace(categorical(my_transition[:,cur]), l)
    end
end

@gen function model(alpha::Float64, len::Int)
    @trace(generate_original_text(alpha, len), :text)
    code = Int[@trace(uniform_discrete(1, 27), (:code, l)) for l=1:27]
    code
end

par_model = Map(model)

@gen function swap_proposal(prev_trace, m::Int)
    @trace(uniform_discrete(1, 27), :i)
    @trace(uniform_discrete(1, 27), :j)
end

function swap_involution(trace, fwd_choices::ChoiceMap, fwd_ret, proposal_args::Tuple)
    choices = get_choices(trace)
    model_args = get_args(trace)
    (replica,) = proposal_args
    (alphas, len_repeated) = model_args
    len::Int = len_repeated[1]
    bwd_choices = choicemap()
    i = fwd_choices[:i]
    j = fwd_choices[:j]
    bwd_choices[:j] = i
    bwd_choices[:i] = j
    constraints = choicemap()
    constraints[replica => (:code, i)] = choices[replica => (:code, j)]
    constraints[replica => (:code, j)] = choices[replica => (:code, i)]

    # update the latent letters
    for l=1:len
        local cur_char::Int
        local new_char::Int
        cur_char = choices[replica => :text => l]
        if cur_char == i
            new_char = j
        elseif cur_char == j
            new_char = i
        else
            new_char = cur_char
        end
        constraints[replica => :text => l] = new_char
    end

    (new_trace, weight, _, _) = update(trace, model_args, (NoChange(),), constraints)
    (new_trace, bwd_choices, weight)
end

@gen function exchange_proposal(prev_trace, m::Int)
    # exchange change m with (((m-1)+1)%n)+1
end

function exchange_involution(trace, fwd_choices::ChoiceMap, fwd_ret, proposal_args::Tuple)
    choices = get_choices(trace)
    model_args = get_args(trace)
    (alphas, _) = model_args
    (replica,) = proposal_args
    constraints = choicemap()
    num_replicas = length(alphas)
    replica_plus_one = (((replica-1)+1)%num_replicas)+1
    set_submap!(constraints, replica_plus_one, get_submap(choices, replica))
    set_submap!(constraints, replica, get_submap(choices, replica_plus_one))
    (new_trace, weight, _, _) = update(trace, model_args, (NoChange(),), constraints)
    (new_trace, EmptyChoiceMap(), weight)
end

to_string(text::Vector{Int}) = join(alphabet[text])

function get_original_text(choices::ChoiceMap, len::Int)
    Int[choices[:text => l] for l=1:len]
end

function get_output_text(original_text::Vector{Int}, code::Vector{Int})
    code[original_text]
end

function do_inference(encoded_text::AbstractString, num_iter::Int)
    len = length(encoded_text)
    choices = choicemap()

    alphas = collect(range(0, stop=1, length=10))
    num_replicas = length(alphas)

    for replica=1:num_replicas
        # set initial code to the identity
        for l=1:27
            choices[replica => (:code, l)] = l
        end

        # initialize original text
        for (l, char) in enumerate(encoded_text)
            choices[replica => :text => l] = letter_to_int[char]
        end
    end

    # initial trace
    (trace, _) = generate(par_model, (alphas, fill(len, num_replicas)), choices)

    # do MCMC
    for iter=1:num_iter

        # print state
        if (iter - 1) % 100 == 1
            println("iter: $iter")
            retval = get_retval(trace)
            choices = get_choices(trace)
            @assert length(retval) == length(alphas)
            println()
            println()
            for replica=1:num_replicas
                code = retval[replica]
                original_text = get_original_text(get_submap(choices, replica), len)
                output_text = get_output_text(original_text, code)
                println("$(to_string(original_text[1:120]))...")
            end
        end

        for replica=1:num_replicas
            (trace, _) = metropolis_hastings(trace, swap_proposal, (replica,), swap_involution)
            (trace, _) = metropolis_hastings(trace, exchange_proposal, (replica,), exchange_involution)
        end
    end
end

using Random: randperm, seed!

seed!(1)
seed!(2)
seed!(3)

original_text = join(split("""
to be or not to be that is the question
whether tis nobler in the mind to suffer
the slings and arrows of outrageous fortune
or to take arms against a sea of troubles
and by opposing end them to die to sleep
no more and by a sleep to say we end
the heartache and the thousand natural shocks
that flesh is heir to tis a consummation
devoutly to be wishd to die to sleep
to sleep perchance to dream ay theres the rub
for in that sleep of death what dreams may come
when we have shuffled off this mortal coil
must give us pause theres the respect
that makes calamity of so long life
for who would bear the whips and scorns of time
the oppressors wrong the proud mans contumely
the pangs of despised love the laws delay
the insolence of office and the spurns
that patient merit of the unworthy takes
when he himself might his quietus make
with a bare bodkin who would fardels bear
to grunt and sweat under a weary life
but that the dread of something after death
the undiscovered country from whose bourn
""", "\n"), " ")

println("original text:")
println(original_text)
code = randperm(27)
original_text_int = map((char) -> letter_to_int[char], collect(original_text))
encoded_text = join(map((letter_int) -> alphabet[letter_int], code[original_text_int]))
println("encoded text:")
println(encoded_text)
do_inference(encoded_text, 1000)

max_stars_repo_name: georgematheos/Gen

max_issues_repo_name: georgematheos/Gen

max_forks_repo_name: georgematheos/Gen


text:
const condadir = abspath(first(DEPOT_PATH), "conda")
const condadeps = joinpath(condadir, "deps.jl")

module DefaultDeps
    if isfile("deps.jl")
        include("deps.jl")
    elseif isfile(Main.condadeps)
        include(Main.condadeps)
    end
    if !isdefined(@__MODULE__, :MINICONDA_VERSION)
        const MINICONDA_VERSION = "3"
    end
    if !isdefined(@__MODULE__, :ROOTENV)
        const ROOTENV = joinpath(Main.condadir, MINICONDA_VERSION)
    end

    USE_MINIFORGE_DEFAULT = true
    if Sys.ARCH in [:x86, :i686]
        USE_MINIFORGE_DEFAULT = false
        @warn """The free/open-source Miniforge (i.e. the conda-forge channel) does not support this platform.
Using the Anaconda/defaults channel instead, which is free for non-commercial use but otherwise may require a license.
        """
    end
    if !isdefined(@__MODULE__, :USE_MINIFORGE)
        const USE_MINIFORGE = USE_MINIFORGE_DEFAULT
    end
    function default_conda_exe(ROOTENV)
        @static if Sys.iswindows()
            p = joinpath(ROOTENV, "Scripts")
            conda_bat = joinpath(p, "conda.bat")
            isfile(conda_bat) ? conda_bat : joinpath(p, "conda.exe")
        else
            joinpath(ROOTENV, "bin", "conda")
        end
    end

    if !isdefined(@__MODULE__, :CONDA_EXE)
        const CONDA_EXE = default_conda_exe(ROOTENV)
    end
end

MINICONDA_VERSION = get(ENV, "CONDA_JL_VERSION", DefaultDeps.MINICONDA_VERSION)
ROOTENV = get(ENV, "CONDA_JL_HOME") do
    root = DefaultDeps.ROOTENV

    # Ensure the ROOTENV uses the current MINICONDA_VERSION when not using a custom ROOTENV
    if normpath(dirname(root)) == normpath(condadir) && all(isdigit, basename(root))
        joinpath(condadir, MINICONDA_VERSION)
    else
        root
    end
end

USE_MINIFORGE = lowercase(get(ENV, "CONDA_JL_USE_MINIFORGE", DefaultDeps.USE_MINIFORGE ? "1" : "0")) in ("1","true","yes")

if isdir(ROOTENV) && MINICONDA_VERSION != DefaultDeps.MINICONDA_VERSION
    error("""Miniconda version changed, since last build.
However, a root enviroment already exists at $(ROOTENV).
Setting Miniconda version is not supported for existing root enviroments.
To leave Miniconda version as, it is unset the CONDA_JL_VERSION enviroment variable and rebuild.
To change Miniconda version, you must delete the root enviroment and rebuild.
WARNING: deleting the root enviroment will delete all the packages in it.
This will break many Julia packages that have used Conda to install their dependancies.
These will require rebuilding.
""")
end

CONDA_EXE = get(ENV, "CONDA_JL_CONDA_EXE") do
    if ROOTENV == DefaultDeps.ROOTENV
        DefaultDeps.CONDA_EXE
    else
        DefaultDeps.default_conda_exe(ROOTENV)
    end
end

if haskey(ENV, "CONDA_JL_CONDA_EXE")
    # Check to see if CONDA_EXE is an executable file
    if isfile(CONDA_EXE)
        if Sys.isexecutable(CONDA_EXE)
            @info "Executable conda located." CONDA_EXE
        else
            error("CONDA_JL_CONDA_EXE, $CONDA_EXE, cannot be executed by the current user.")
        end
    else
        error("CONDA_JL_CONDA_EXE, $CONDA_EXE, does not exist.")
    end
else
    if !isfile(CONDA_EXE)
        # An old CONDA_EXE has gone missing, revert to default in ROOTENV
        @info "CONDA_EXE not found. Reverting to default in ROOTENV" CONDA_EXE ROOTENV
        CONDA_EXE = DefaultDeps.default_conda_exe(ROOTENV)
    end
end

deps = """
const ROOTENV = "$(escape_string(ROOTENV))"
const MINICONDA_VERSION = "$(escape_string(MINICONDA_VERSION))"
const USE_MINIFORGE = $USE_MINIFORGE
const CONDA_EXE = "$(escape_string(CONDA_EXE))"
"""

mkpath(condadir)
mkpath(ROOTENV)

for depsfile in ("deps.jl", condadeps)
    if !isfile(depsfile) || read(depsfile, String) != deps
        write(depsfile, deps)
    end
end


text:
################################################################################
# Toy Example
################################################################################
T = Float64

# Define the dynamics of the system
p = 2 # Number of players
# model = DoubleIntegratorGame(p=p) # game with 3 players with unicycle dynamics
model = UnicycleGame(p=p) # game with 3 players with unicycle dynamics
model = BicycleGame(p=p) # game with 3 players with unicycle dynamics
n = model.n
m = model.m

# Define the horizon of the problem
N = 20 # N time steps
dt = 0.1 # each step lasts 0.1 second
probsize = ProblemSize(N,model) # Structure holding the relevant sizes of the problem

# Define the objective of each player
# We use a LQR cost
Q = [Diagonal(10*ones(SVector{model.ni[i],T})) for i=1:p] # Quadratic state cost
R = [Diagonal(0.1*ones(SVector{model.mi[i],T})) for i=1:p] # Quadratic control cost
# Desrired state
xf = [SVector{model.ni[1],T}([2,+0.4,0,0]),
      SVector{model.ni[2],T}([2, 0.0,0,0]),
      # SVector{model.ni[3],T}([3,-0.4,0,0]),
      # SVector{model.ni[4],T}([3,+0.8,0,0]),
      ]
# Desired control
uf = [zeros(SVector{model.mi[i],T}) for i=1:p]

# # Objectives of the game
game_obj = Algames.GameObjective(Q,R,xf,uf,N,model)
radius = 0.5*ones(p)
μ = 4.0*ones(p)
add_collision_cost!(game_obj, radius, μ)

# Define the constraints that each player must respect
game_con = Algames.GameConstraintValues(probsize)
# # Add collision avoidance
# radius = 0.05
# add_collision_avoidance!(game_con, radius)

# # Add control bounds
# u_max =  5*ones(SVector{m,T})
# u_min = -5*ones(SVector{m,T})
# add_control_bound!(game_con, u_max, u_min)

# # Add wall constraint
# walls = [Wall([-1.0,-1.75], [1.0,0.25], [1.,-1.]/sqrt(2))]
# add_wall_constraint!(game_con, walls)

# # Add circle constraint
# xc = [1., 2., 3.]
# yc = [1., 2., 3.]
# radius = [0.1, 0.2, 0.3]
# add_circle_constraint!(game_con, xc, yc, radius)

# Define the initial state of the system
x0 = SVector{model.n,T}([
    0.0, 0.0,# 0.5, #0.0,
   -0.4, 0.0,# 0.7, #0.6,
    0.0, 0.0,# 0.0, #0.0,
    0.0, 0.0,# 0.0, #0.0,
    ])

# Define the Options of the solver
opts = Options()
opts.ls_iter = 15
opts.outer_iter = 20
opts.inner_iter = 20
opts.ρ_0 = 1e0
opts.reg_0 = 1e-8
opts.α_dual = 1.0
opts.λ_max = 1.0*1e7
opts.ϵ_dyn = 1e-3
opts.ϵ_sta = 1e-3
opts.ϵ_con = 1e-3
opts.ϵ_opt = 1e-3
opts.regularize = true
# Define the game problem
prob = Algames.GameProblem(N,dt,x0,model,opts,game_obj,game_con)

# Solve the problem
@time newton_solve!(prob)
# @profiler newton_solve!(prob)

plot!(prob.model, prob.pdtraj.pr)
plot!(prob.stats)

prob.stats.outer_iter
# kick out of inner loop if not progress made
# add probsize to game_con struct
# add constructor for multiple radius of collision

function traj_dist(prob::Algames.GameProblem ;M::Int=10)
    traj = []
    for j = 1:M
        prob.opts.seed = j
        newton_solve!(prob)
        push!(traj, deepcopy(prob.pdtraj.pr))
    end
    plotset(prob.model, traj)
    return traj
end

function plotset(model::AbstractGameModel, traj)
    plt = plot(aspect_ratio=:equal)
    M = length(traj)
    N = length(traj[1])
    col = [:blue, :orange, :green]
    for j = 1:M
        for i = 1:model.p
            xi = [Algames.state(traj[j][k])[model.pz[i][1]] for k=1:N]
            yi = [Algames.state(traj[j][k])[model.pz[i][2]] for k=1:N]
            plot!(xi, yi, color=col[i], label=false)
            scatter!(xi, yi, color=col[i], label=false)
        end
    end
    display(plt)
    return nothing
end

function solver_progress(prob::Algames.GameProblem; M::Int=10, II_max::Int=2, IO::Int=1)
    # anim = @animate for ii = II_max:II_max
    anim = @animate for ii = 0:II_max
        prob.opts.outer_iter=IO
        prob.opts.inner_iter=ii
        traj_dist(prob, M=M)
    end
    gif(anim, "Documents/solver_progress_Uni_nncvx.gif", fps = 2)
    return nothing
end

prob.opts.amplitude_init = 3.0
prob.opts.f_init = x -> (rand(x) .- 0.5)
prob.opts.outer_iter = 0
prob.opts.inner_iter = 0
traj_dist(prob, M=100)

solver_progress(prob, M=100, II_max=20, IO=1)
solver_progress(prob, M=100, II_max=20, IO=20)


text:
# ---
# title: 1175. Prime Arrangements
# id: problem1175
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Math
# link: <https://leetcode.com/problems/prime-arrangements/description/>
# hidden: true
# ---
#
# Return the number of permutations of 1 to `n` so that prime numbers are at
# prime indices (1-indexed.)
#
# _(Recall that an integer  is prime if and only if it is greater than 1, and
# cannot be written as a product of two positive integers both smaller than
# it.)_
#
# Since the answer may be large, return the answer **modulo`10^9 + 7`**.
#
#
#
# **Example 1:**
#
#
#
#     Input: n = 5
#     Output: 12
#     Explanation: For example [1,2,5,4,3] is a valid permutation, but [5,2,3,4,1] is not because the prime number 5 is at index 1.
#
#
# **Example 2:**
#
#
#
#     Input: n = 100
#     Output: 682289015
#
#
#
#
# **Constraints:**
#
#   * `1 <= n <= 100`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
# tests AdhesionForce functions

# Functions for computing particle-wall contact forces of
# via LJ potentials and a quadrature

# loads in Point2D type
include("../../src/part/CellLists.jl")
# load wall types
include("../../src/part/ParticleTypes.jl")
# loads in LJ force function
include("../../src/part/CohesionForce.jl")
# load onSegment
include("../../src/part/isInside.jl")

# linewall
function NearestPoint!_linewall_TEST()
    point = Point2D(0.0,0.0)
    node = Point2D(0.0,1.0)
    linestart = Point2D(-1.0,0.0)
    lineend   = Point2D(1.0,0.0)
    wall = LineWall([linestart,lineend])

    #@code_warntype NearestPoint!(point,node,wall)
    @btime NearestPoint!($point,$node,$wall)

    nothing
end
# circlewall
function NearestPoint!_circlewall_TEST()
    point = Point2D(0.0,0.0)
    node = Point2D(0.0,1.0)
    center = Point2D(1.4,0.0)
    radius = 0.2
    wall = CircleWall(center,radius)

    #@code_warntype NearestPoint!(point,node,wall)
    @btime NearestPoint!($point,$node,$wall)

    println(point)

    nothing
end
# arcwall
function NearestPoint!_arcwall_TEST()
    point  = Point2D(0.0,0.0)
    node   = Point2D(0.0,1.0)

    r1     = Point2D(1.6,0.0)
    center = Point2D(1.4,0.0)
    r2     = Point2D(1.4,0.2)
    wall   = ArcWall(r1,center,r2)

    #@code_warntype NearestPoint!(point,node,wall)
    @btime NearestPoint!($point,$node,$wall)

    nothing
end

function isCloseEnough_TEST()
    node  = Point2D(0.0,1.0)
    point = Point2D(0.0,0.0)
    radius = 0.1
    k = 5.0
    rm = 1.5

    #@code_warntype isCloseEnough(point,node,k,rm,radius)
    @btime isCloseEnough($point,$node,$k,$rm,$radius)

    nothing
end

# linewall test
function generateQuadNodes!_linewall_TEST()
    node = Point2D(0.0,1.0)
    linestart = Point2D(-1.0,0.0)
    lineend   = Point2D(1.0,0.0)
    wall = LineWall([linestart,lineend])
    point = NearestPoint(node,wall)

    Nquad = 20
    k = 5.0
    rm = 1.5

    # initialize arrays
    xquad = zeros(Nquad); yquad = zeros(Nquad)

    #@code_warntype generateQuadNodes!(xquad,yquad,point,wall,node,k,rm)
    @btime generateQuadNodes!($xquad,$yquad,$point,$wall,$k,$rm)

    nothing
end
# circlewall
function generateQuadNodes!_circlewall_TEST()
    node   = Point2D(0.0,1.0)
    center = Point2D(1.4,0.0)
    radius = 0.2
    wall = CircleWall(center,radius)
    point = NearestPoint(node,wall)

    Nquad = 20
    k = 5.0
    rm = 1.5

    # initialize vectors to be assigned in-place
    xquad = zeros(Nquad)
    yquad = zeros(Nquad)

    #@code_warntype generateQuadNodes!(xquad,yquad,point,wall,k,rm)
    @btime generateQuadNodes!($xquad,$yquad,$point,$wall,$k,$rm)

    nothing
end
# arcwall
function generateQuadNodes!_arcwall_TEST()
    node   = Point2D(0.0,1.0)
    r1     = Point2D(1.6,0.0)
    center = Point2D(1.4,0.0)
    r2     = Point2D(1.4,0.2)
    wall   = ArcWall(r1,center,r2)
    point = NearestPoint(node,wall)

    Nquad = 20
    k = 5.0
    rm = 1.5

    # initialize vectors to be assigned in-place
    xquad = zeros(Nquad)
    yquad = zeros(Nquad)

    #@code_warntype generateQuadNodes!(xquad,yquad,point,wall,k,rm)
    @btime generateQuadNodes!($xquad,$yquad,$point,$wall,$k,$rm)

    nothing
end

# linewall
function isInLine_linewall_TEST()
    node = Point2D(0.0,1.0)
    linestart = Point2D(-1.0,0.0)
    lineend   = Point2D(1.0,0.0)
    wall = LineWall([linestart,lineend])
    point = NearestPoint(node,wall)

    Nquad = 20
    k = 5.0
    rm = 1.5

    # initialize vectors to be assigned in-place
    xquad = zeros(Nquad)
    yquad = zeros(Nquad)

    generateQuadNodes!(xquad,yquad,point,wall,k,rm)

    # initialize vectors so they don't have to keep being re-allocated
    s = Point2D(0.0,0.0)

    #@code_warntype isInLine(wall,xquad[2],yquad[2],s)
    @btime isInLine($wall,$xquad[2],$yquad[2],$s)

    nothing
end
# circlewall
function isInLine_circlewall_TEST()
    node   = Point2D(0.0,1.0)
    center = Point2D(1.4,0.0)
    radius = 0.2
    wall = CircleWall(center,radius)
    point = NearestPoint(node,wall)

    Nquad = 20
    k = 5.0
    rm = 1.5

    # initialize vectors to be assigned in-place
    xquad = zeros(Nquad)
    yquad = zeros(Nquad)

    generateQuadNodes!(xquad,yquad,point,wall,k,rm)

    # initialize vectors so they don't have to keep being re-allocated
    s = Point2D(0.0,0.0)

    #@code_warntype isInLine(wall,xquad[2],yquad[2],s)
    @btime isInLine($wall,$xquad[2],$yquad[2],s)

    nothing
end
# arcwall
function isInLine_arcwall_TEST()
    node   = Point2D(0.0,1.0)
    r1     = Point2D(1.6,0.0)
    center = Point2D(1.4,0.0)
    r2     = Point2D(1.4,0.2)
    wall   = ArcWall(r1,center,r2)
    point = NearestPoint(node,wall)

    Nquad = 20
    k = 5.0
    rm = 1.5

    # initialize vectors to be assigned in-place
    xquad = zeros(Nquad)
    yquad = zeros(Nquad)

    generateQuadNodes!(xquad,yquad,point,wall,k,rm)

    # initialize vectors so they don't have to keep being re-allocated
    s = Point2D(0.0,0.0)

    #@code_warntype isInLine(wall,xquad[2],yquad[2],s)
    @btime isInLine($wall,$xquad[2],$yquad[2],$s)

    nothing
end

# linewall
function wallTrapQuad_linewall_TEST()
    node = Point2D(0.0,0.2)
    linestart = Point2D(-1.0,0.0)
    lineend   = Point2D(1.0,0.0)
    wall = LineWall([linestart,lineend])
    point = NearestPoint(node,wall)

    particleradius = 0.1
    Nquad = 20
    k = 5.0
    rm = 1.5            # whats difference between rm and rc???
    ϵ = 1.0
    rc = 1.5

    # initialize vectors to be assigned in-place
    xquad = zeros(Nquad)
    yquad = zeros(Nquad)

    generateQuadNodes!(xquad,yquad,point,wall,k,rm)

    s = Point2D(0.0,0.0)

    fx,fy = wallTrapQuad(node,wall,xquad,yquad,ϵ,particleradius,rc,s)
    #@code_warntype wallTrapQuad(node,wall,xquad,yquad,ϵ,particleradius,rc,s)
    @btime wallTrapQuad($node,$wall,$xquad,$yquad,$ϵ,$particleradius,$rc,$s)

    return fx,fy
end

# circle wall
function wallTrapQuad_circlewall_TEST()
    node   = Point2D(1.1,0.0)
    center = Point2D(1.4,0.0)
    radius = 0.2
    wall = CircleWall(center,radius)
    point = NearestPoint(node,wall)

    particleradius = 0.1
    Nquad = 20
    k = 5.0
    rm = 1.5            # whats difference between rm and rc???
    ϵ = 1.0
    rc = 1.5

    # initialize vectors to be assigned in-place
    xquad = zeros(Nquad)
    yquad = zeros(Nquad)

    generateQuadNodes!(xquad,yquad,point,wall,k,rm)

    # initialize vectors so they don't have to keep being re-allocated
    s = Point2D(0.0,0.0)

    fx,fy = wallTrapQuad(node,wall,xquad,yquad,ϵ,particleradius,rc,s)
    #@code_warntype wallTrapQuad(node,wall,xquad,yquad,ϵ,particleradius,rc,s)
    @btime wallTrapQuad($node,$wall,$xquad,$yquad,$ϵ,$particleradius,$rc,$s)

    return fx,fy
end

# arc wall
function wallTrapQuad_arcwall_TEST()
    node   = Point2D(1.4,0.0)
    r1     = Point2D(1.6,0.0)
    center = Point2D(1.4,0.0)
    r2     = Point2D(1.4,0.2)
    wall   = ArcWall(r1,center,r2)
    point = NearestPoint(node,wall)

    particleradius = 0.1
    Nquad = 20
    k = 5.0
    rm = 1.5            # whats difference between rm and rc???
    ϵ = 1.0
    rc = 1.5

    # initialize vectors to be assigned in-place
    xquad = zeros(Nquad)
    yquad = zeros(Nquad)

    generateQuadNodes!(xquad,yquad,point,wall,k,rm)

    # initialize vectors so they don't have to keep being re-allocated
    s = Point2D(0.0,0.0)

    fx,fy = wallTrapQuad(node,wall,xquad,yquad,ϵ,particleradius,rc,s)
    #@code_warntype wallTrapQuad(node,wall,xquad,yquad,ϵ,particleradius,rc,s)
    @btime wallTrapQuad($node,$wall,$xquad,$yquad,$ϵ,$particleradius,$rc,$s)

    return fx,fy
end

max_stars_repo_name: pseastham/eFEMpart


text:
"""
---------------------------------------------------
Types and Constructors
---------------------------------------------------
"""

mutable struct NOMINAL_LANDING_SAMPLER <: SAMPLER
	z_init::Sampleable
	vx_init::Sampleable
	vz_init::Sampleable
	θ::Sampleable
	t_extra::Sampleable
end

function nominal_landing_sampler(;z_init = Uniform(300,600),
						vx_init = Uniform(30,45),
						vz_init = Uniform(-500fpm2fps,-300fpm2fps),
						θ = Truncated(Normal(6,3),2,15),
						t_extra = Uniform(0.1,0.3)) #Normal(10,3))
	return NOMINAL_LANDING_SAMPLER(z_init, vx_init, vz_init, θ, t_extra)
end

mutable struct NOMINAL_LANDING <: UAM_TRAJECTORY
	z_init::Float64
	vx_init::Float64
	vz_init::Float64
	θ::Float64
	t_extra::Float64
	sampler::SAMPLER
	dt::Float64
	vz_max::Float64
	az_max::Float64
	ax_max::Float64
	p::Array{Float64,2}
	v::Array{Float64,2}
end

function nominal_landing(;z_init = 0.0,
						vx_init = 0.0,
						vz_init = 0.0,
						θ = 0.0,
						t_extra = 0.0,
						sampler = nominal_landing_sampler(),
						dt = 1.0,
						vz_max = 550fpm2fps,
						az_max = 0.3g_ft,
						ax_max = 0.1g,
						p = zeros(1,3),
						v = zeros(1,3))
	return NOMINAL_LANDING(z_init, vx_init, vz_init, θ, t_extra, sampler, dt,
							vz_max, az_max, ax_max, p, v)
end

"""
---------------------------------------------------
Functions
---------------------------------------------------
"""

function solve_trajectory!(τ::NOMINAL_LANDING)
	# Get some needed info
	τ.vz_init = clamp(τ.vz_init, -τ.vz_max, Inf)
	ttot = abs(τ.z_init/τ.vz_init) + τ.t_extra*abs(τ.z_init/τ.vz_init)
	N = convert(Int64, floor(ttot/τ.dt))
	# Difference matrix
	D = zeros(N-1,N)
	for i = 1:(N-1)
	    D[i,i] = -1
	    D[i,i+1] = 1
	end
	D = D./τ.dt

	# Create and solve the optimization problem
	px = Variable(N)
	pz = Variable(N)
	vx = Variable(N)
	vz = Variable(N)

	objec = norm(D*vx) + λ*norm(vz*ft2m + vx*tand(τ.θ))

	constraints = [[px[i+1] == px[i] + vx[i]*τ.dt for i = 1:N-1];
				   [pz[i+1] == pz[i] + vz[i]*τ.dt for i = 1:N-1];
				   vx[1] == τ.vx_init; px[1] == 0;
				   vz[1] == τ.vz_init; pz[1] == τ.z_init;
				   vx[N] == 0; pz[N] == 0; vz[N] == 0;
				   vx ≥ 0; pz ≥ 0; vz ≤ 0; vz ≥ -τ.vz_max;
				   D*vx ≤ τ.ax_max; D*vz*ft2m ≤ τ.az_max]

	prob = minimize(objec, constraints)
	solve!(prob, ECOS.Optimizer(verbose=0))

	# Update the trajectory
	p = zeros(N,3)
	p[:,1] = px.value[:,1]
	p[:,3] = pz.value[:,1]
	v = zeros(N,3)
	v[:,1] = vx.value[:,1]
	v[:,3] = vz.value[:,1]

	τ.p = copy(p)
	τ.v = copy(v)

	return Int(prob.status) == 1
end

max_stars_repo_name: Airspace-Encounter-Models/em-model-uam

max_issues_repo_name: Airspace-Encounter-Models/em-model-uam

max_forks_repo_name: Airspace-Encounter-Models/em-model-uam


text:
# Inside make.jl
using Documenter, Proportions
makedocs(
         sitename = "Proportions.jl",
         modules  = [Proportions],
         doctest = true,
         pages=[
                "Home" => "index.md"
               ])
deploydocs(;
    repo="github.com/hillelawaskar/Proportions.jl",
)

max_stars_repo_name: hillelawaskar/Proportions

max_issues_repo_name: hillelawaskar/Proportions