text:
export @optic
export set, modify
export ∘, opcompose, var"⨟"
export Elements, Recursive, If, Properties
export setproperties
export constructorof
using ConstructionBase
using CompositionsBase
using Base: getproperty
using Base

const EXPERIMENTAL = """This function/method/type is experimental. It can be changed or deleted at any point without warning"""

"""
    modify(f, obj, optic)

Replace a part `x` of `obj` by `f(x)`. The `optic` argument selects
which part to replace.

```jldoctest
julia> using Accessors

julia> obj = (a=1, b=2); optic=@optic _.a; f = x -> "hello \$x";

julia> modify(f, obj, optic)
(a = "hello 1", b = 2)
```
See also [`set`](@ref).
"""
function modify end

"""
    set(obj, optic, val)

Replace a part according to `optic` of `obj` by `val`.

```jldoctest
julia> using Accessors

julia> obj = (a=1, b=2); lens=@optic _.a; val = 100;

julia> set(obj, lens, val)
(a = 100, b = 2)
```
See also [`modify`](@ref).
"""
function set end

"""
    optic₁ ⨟ optic₂

Compose optics `optic₁`, `optic₂`, ..., `opticₙ` to access nested objects.

# Example
```jldoctest
julia> using Accessors

julia> obj = (a = (b = (c = 1,),),);

julia> la = @optic _.a
       lb = @optic _.b
       lc = @optic _.c
       lens = la ⨟ lb ⨟ lc
(@optic _.c) ∘ (@optic _.b) ∘ (@optic _.a)

julia> lens(obj)
1
```
"""
opcompose

const BASE_COMPOSED_FUNCTION_HAS_SHOW = VERSION >= v"1.6.0-DEV.85"
const BASE_COMPOSED_FUNCTION_IS_PUBLIC = VERSION >= v"1.6.0-DEV.1037"
if !BASE_COMPOSED_FUNCTION_IS_PUBLIC
    using Compat: ComposedFunction
end
if !BASE_COMPOSED_FUNCTION_HAS_SHOW
    function show_composed_function(io::IO, c::ComposedFunction)
        show(io, c.outer)
        print(io, " ∘ ")
        show(io, c.inner)
    end
    function Base.show(io::IO, c::ComposedFunction)
        show_composed_function(io, c)
    end
    function Base.show(io::IO, ::MIME"text/plain", c::ComposedFunction)
        show_composed_function(io, c)
    end
end

const ComposedOptic{Outer, Inner} = ComposedFunction{Outer, Inner}
outertype(::Type{ComposedOptic{Outer, Inner}}) where {Outer, Inner} = Outer
innertype(::Type{ComposedOptic{Outer, Inner}}) where {Outer, Inner} = Inner

# TODO better name
# also better way to organize traits will
# probably only emerge over time
abstract type OpticStyle end
struct ModifyBased <: OpticStyle end
struct SetBased <: OpticStyle end
OpticStyle(obj) = OpticStyle(typeof(obj))
# defining lenses should be very lightweight
# e.g. only a single `set` implementation
# so we choose this as the default trait
OpticStyle(::Type{T}) where {T} = SetBased()

function OpticStyle(::Type{ComposedOptic{O,I}}) where {O,I}
    composed_optic_style(OpticStyle(O), OpticStyle(I))
end
composed_optic_style(::SetBased, ::SetBased) = SetBased()
composed_optic_style(::ModifyBased, ::SetBased) = ModifyBased()
composed_optic_style(::SetBased, ::ModifyBased) = ModifyBased()
composed_optic_style(::ModifyBased, ::ModifyBased) = ModifyBased()

@inline function set(obj, optic, val)
    _set(obj, optic, val, OpticStyle(optic))
end

function _set(obj, optic, val, ::SetBased)
    Optic = typeof(optic)
    error("""
    This should be unreachable. You probably need to overload
    `Accessors.set(obj, ::$Optic, val)
    """
   )
end

struct Constant{V}
    value::V
end
(o::Constant)(x) = o.value

@inline function _set(obj, optic, val, ::ModifyBased)
    modify(Constant(val), obj, optic)
end

@inline function _set(obj, optic::ComposedOptic, val, ::SetBased)
    inner_obj = optic.inner(obj)
    inner_val = set(inner_obj, optic.outer, val)
    set(obj, optic.inner, inner_val)
end

@inline function modify(f, obj, optic)
    _modify(f, obj, optic, OpticStyle(optic))
end

function _modify(f, obj, optic, ::ModifyBased)
    Optic = typeof(optic)
    error("""
          This should be unreachable. You probably need to overload:
          `Accessors.modify(f, obj, ::$Optic)`
          """)
end

function _modify(f, obj, optic::ComposedOptic, ::ModifyBased)
    otr = optic.outer
    inr = optic.inner
    modify(obj, inr) do o1
        modify(f, o1, otr)
    end
end

@inline function _modify(f, obj, optic, ::SetBased)
    set(obj, optic, f(optic(obj)))
end

"""
    Elements

Access all elements of a collection that implements `map`.

```jldoctest
julia> using Accessors

julia> obj = (1,2,3);

julia> set(obj, Elements(), 0)
(0, 0, 0)

julia> modify(x -> 2x, obj, Elements())
(2, 4, 6)
```
$EXPERIMENTAL
"""
struct Elements end
OpticStyle(::Type{<:Elements}) = ModifyBased()

function modify(f, obj, ::Elements)
    map(f, obj)
end

"""
    If(modify_condition)

Restric access to locations for which `modify_condition` holds.

```jldoctest
julia> using Accessors

julia> obj = (1,2,3,4,5,6);

julia> @set obj |> Elements() |> If(iseven) *= 10
(1, 20, 3, 40, 5, 60)
```

$EXPERIMENTAL
"""
struct If{C}
    modify_condition::C
end
OpticStyle(::Type{<:If}) = ModifyBased()

function modify(f, obj, w::If)
    if w.modify_condition(obj)
        f(obj)
    else
        obj
    end
end

"""
    mapproperties(f, obj)

Construct a copy of `obj`, with each property replaced by
the result of applying `f` to it.

```jldoctest
julia> using Accessors

julia> obj = (a=1, b=2);

julia> Accessors.mapproperties(x -> x+1, obj)
(a = 2, b = 3)
```
$EXPERIMENTAL
"""
function mapproperties(f, obj)
    # TODO move this helper elsewhere?
    # TODO should we use a generated function based on fieldnames?
    pnames = propertynames(obj)
    if isempty(pnames)
        return obj
    else
        ctor = constructorof(typeof(obj))
        new_props = map(pnames) do p
            f(getproperty(obj, p))
        end
        return ctor(new_props...)
    end
end

"""
    Properties()

Access all properties of an objects.

```jldoctest
julia> using Accessors

julia> obj = (a=1, b=2, c=3)
(a = 1, b = 2, c = 3)

julia> set(obj, Properties(), "hi")
(a = "hi", b = "hi", c = "hi")

julia> modify(x -> 2x, obj, Properties())
(a = 2, b = 4, c = 6)
```
Based on [`mapproperties`](@ref).

$EXPERIMENTAL
"""
struct Properties end
OpticStyle(::Type{<:Properties}) = ModifyBased()
modify(f, o, ::Properties) = mapproperties(f, o)

"""
    Recursive(descent_condition, optic)

Apply `optic` recursively as long as `descent_condition` holds.
```jldoctest
julia> using Accessors

julia> obj = (a=missing, b=1, c=(d=missing, e=(f=missing, g=2)))
(a = missing, b = 1, c = (d = missing, e = (f = missing, g = 2)))

julia> set(obj, Recursive(!ismissing, Properties()), 100)
(a = 100, b = 1, c = (d = 100, e = (f = 100, g = 2)))

julia> obj = (1,2,(3,(4,5),6))
(1, 2, (3, (4, 5), 6))

julia> modify(x -> 100x, obj, Recursive(x -> (x isa Tuple), Elements()))
(100, 200, (300, (400, 500), 600))
```
"""
struct Recursive{Descent, Optic}
    descent_condition::Descent
    optic::Optic
end
OpticStyle(::Type{Recursive{D,O}}) where {D,O} = ModifyBased() # Is this a good idea?

function _modify(f, obj, r::Recursive, ::ModifyBased)
    modify(obj, r.optic) do o
        if r.descent_condition(o)
            modify(f, o, r)
        else
            f(o)
        end
    end
end

################################################################################
##### Lenses
################################################################################
struct PropertyLens{fieldname} end

function (l::PropertyLens{field})(obj) where {field}
    getproperty(obj, field)
end

@inline function set(obj, l::PropertyLens{field}, val) where {field}
    patch = (;field => val)
    setproperties(obj, patch)
end

struct IndexLens{I <: Tuple}
    indices::I
end

Base.@propagate_inbounds function (lens::IndexLens)(obj)
    getindex(obj, lens.indices...)
end
Base.@propagate_inbounds function set(obj, lens::IndexLens, val)
    setindex(obj, val, lens.indices...)
end

struct DynamicIndexLens{F}
    f::F
end

Base.@propagate_inbounds function (lens::DynamicIndexLens)(obj)
    return obj[lens.f(obj)...]
end

Base.@propagate_inbounds function set(obj, lens::DynamicIndexLens, val)
    return setindex(obj, val, lens.f(obj)...)
end


text:
abstract type AbstractBackend; end

struct DummyBackend <: AbstractBackend; end

const backends = AbstractBackend[DummyBackend()]

get_backend() = last(backends)

set_backend!(t::AbstractBackend) = push!(backends, t)

reset_backend!() = pop!(backends)

ext: jl


text:
export AbstractHamiltonian, AbstractHamiltonianIterator

abstract type AbstractHamiltonian end
abstract type AbstractHamiltonianIterator end


text:
using SurrealNumbers
@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

@static if VERSION < v"0.7.0"
    const devnull = DevNull
end

clearcache()

# recursive construction and one and zero
x0 = SurrealFinite("0", ϕ,   ϕ)
x1 = SurrealFinite("1", [x0], ϕ)
x11 = SurrealFinite("-1", ϕ,  [x0])
x111  = SurrealDyadic("-1", ϕ,  [x0])
x4 = SurrealFinite( [x11],  [x1,x0])
x5 = SurrealFinite( [x11,x0,x4],  [x1])
x21  = SurrealFinite("2", [x1],  ϕ)
x22  = SurrealFinite( [x0, x1],  ϕ)
x23  = SurrealFinite( [x11,x1],  ϕ)
x24  = SurrealFinite( [x11,x0,x1],  ϕ)
x25 = [x11,x0,x1] ≀ ϕ
z = zero(x1)
z0 = zero(SurrealFinite)
o = one(x1)
o1 = one(SurrealFinite)
x00 = SurrealFinite([x11],[x1])

x3 = convert(SurrealFinite, 3)
x41 = convert(SurrealFinite, 4)
x42 = SurrealFinite( [x3],  ϕ)
x43 = convert(SurrealFinite, 2) * convert(SurrealFinite, 2)
x44 = SurrealFinite(4)
x6 = convert(SurrealFinite, 2) * convert(SurrealFinite, 3)

print("x1 = ")
pf(x1)
println()

println(" z = ", z)
println("x4 = ", x4)
println("x5 = ", x5)

s1 = convert(SurrealFinite, 1//2)
s2 = convert(SurrealFinite, 3//4)
s2a = convert(SurrealShort, 3//4)
s2b = convert(SurrealDyadic, 3//4)
print("s1 = ")
pf(s1)
println()
print("s2 = ")
pf(s2)
println()

println("x41 = ", x41, ", x42 = ", x42, ", x43 = ", x43)

pf(x43)
println()

spf(x43)
println()

@testset "constructor test" begin
    @test SurrealFinite("1/2", [0], [1]) ≅ 1/2
    @test_throws ErrorException SurrealFinite("1", [x1], [x0])
    @test_throws DomainError convert(SurrealFinite, NaN )
    @test_throws ErrorException SurrealFinite(1//3)
    @test_throws ErrorException SurrealFinite(100000 + 1//2)
    @test_throws ErrorException SurrealFinite(100000.000012)

    @test convert(SurrealFinite, 1.3 ) == 5854679515581645//4503599627370496
    @test all( convert(Array{SurrealFinite}, [1,2]) .== convert.(SurrealFinite, [1,2]) )
    @test z === zero(SurrealFinite)
    @test z0 === zero(SurrealFinite)
    @test o === one(SurrealFinite)
    @test o1 === one(SurrealFinite)
    @test s2a === s2
    @test s2b === s2
    @test x11 == x111
    @test dali(0) === SurrealZero
    @test dali(1) === SurrealOne
    @test dali(-1) == SurrealMinusOne
    @test dali(2) == SurrealTwo
end

@testset "comparisons" begin
    @test x0 <= x1
    @test x0 < x1
    @test x0 <= x0
    @test x1 >= x11
    @test x11 < x0
    @test !( x11 ≅ x1 )
    @test x21 ≅ x22 ≅ x23 ≅ x24 # should be an equivalence class
    @test x41 ≅ x42 ≅ x43
    @test x41 == x42
    @test x44 == x41
    @test x41 != x43

    @test x0 ≺ x5
    @test x4 ≺ x5
    @test x5 ≻ x4
    @test !(x5 ≺ x4)
    @test !(x4 ≺ x4)
    @test x0 ⪯ x5
    @test x5 ⪯ x5
    @test x41 ⪰ x3
    @test s1 ⪰ s1
    @test x1 ⪯ s2
    @test x21 ≺ x43
    @test isancestor(x21, x43)

    @test all( parents(x24) .== [x11,x0,x1] )
    @test all( parents(x5) .== sort([x11,x0,x4,x1]) )
end

@testset "basic operators" begin
    @test -x1 == x11
    @test - -x0 == x0
    @test - -x4 == x4
    # @test - -x4 === x4   # still flaky for version 0.7
    @test x1 + x11 ≅ x0
    @test x1 - x1 ≅ x0
    @test x1 + x0 == x1
    @test x1*x1 ≅ x1
    @test x0*x1 ≅ x0

    @test x4*x0 ≅ x0 # this case broke things before

    @test x11/one(x11) ≅ -x1
    @test x22 / x22 ≅ one(x22)
    @test x22/x1 ≅ x22
    @test x4/x1 ≅ x4
    @test float(x1/x4) == float(x1)/float(x4)
    @test_throws ErrorException x1/x0
    @test_throws ErrorException x1/x3

    @test convert(SurrealFinite, 6)/ convert(SurrealFinite, 3) ≅ convert(SurrealFinite, 2)
end

@testset "convert to surreal" begin
    @test convert(SurrealFinite,2)*convert(SurrealFinite,2) ≅ convert(SurrealFinite,4)
    @test convert(SurrealFinite,2)*convert(SurrealFinite,-3) ≅ convert(SurrealFinite,-6)
    @test convert(SurrealFinite,2)*x0 ≅ x0
    @test convert(SurrealFinite,2)*x1 ≅ convert(SurrealFinite,2)
    @test x4 ≅ -s1
    @test convert(SurrealFinite, 0.5) == convert(SurrealFinite, 1//2 )
    @test convert(SurrealFinite, 1.000000000001) == 1.000000000001
      # tests promotion
      #  float(562949953421875//562949953421312)

    # f =  1122342342.23422522
    # convert(SurrealFinite, 1122342342.23422522)
    @test convert(SurrealFinite, 0.5) == dali(0.5)
end

A = [x0, x1]
B = [x11, x4]
X = [ x0, x1, x1, x11, x1, x11, x23]
Y = copy(X)
unique2!(Y)
@testset "array operations" begin
    @test all( sort(X) .≅ convert.(SurrealFinite, [-1, -1, 0, 1, 1, 1, 2] ))
    @test all( unique(X) .≅ convert.(SurrealFinite, [0, 1, -1, 2] ))
    # @test all( Y .≅ convert.(SurrealFinite, [-1, 0, 1, 2] ))

    @test all( sort(float.(A + B)) .== [-1.0, -0.5, 0.0, 0.5])
    @test all( sort(float.(A - B)) .== [0.5, 1.0, 1.5, 2.0])
    @test A + ϕ == ϕ
    @test ϕ - B == ϕ

    @test ϕ <= A # can see from this that <= is not a proper ordering on 'sets' of surreals
    @test A <= ϕ # its just a convenience comparison operator
    @test ϕ < A
    @test A < ϕ
    @test B <= A
    @test B < A
    @test !(B > A)
    @test !(B >= A)

    @test all( x1*A .≅ A )
    @test all( x23*A .≅ A*x23 )

    @test hash(A) == 0x54b75ebf9d6abfa4
end

V = [-6.25, -4.0, -3.25, -2.5, -2.0, -1.625, -1.0, -0.25, 0.0, 1.0, 7.0, 0.5, 1.625, 2.25, 4.5, 6.75, 8.0, 12.125]
@testset "simple functions" begin
    @test sign(x0) == x0
    @test sign(x1) == x1
    @test sign(x4) == -one(x4)
    @test sign(-x4) == one(x4)

    @test round(x22) ≅ x22
    @test round(x23) ≅ x23
    @test round(x4) ≇ x4
    @test round(-x4) ≇ -x4
    @test round( convert(SurrealFinite, -0.75) ) ≅ round( -0.75 )
    @test round( convert(SurrealFinite, -1.25) ) ≅ round( -1.25 )

    @test floor(zero(SurrealFinite)) ≅ 0.0
    @test floor(one(SurrealFinite)) ≅ 1.0
    @test floor(x22) ≅ x22
    @test floor(x4) ≅ -one(x4)
    @test floor(x41) ≅ 4
    @test floor(x4) ≇ x4
    @test floor(-x4) ≇ -x4
    @test floor(-x4) == zero(x4)
    for v in V
        println("v = $v")
        @test floor( convert(SurrealFinite, v) ) ≅ floor( v )
        @test floor( Int, convert(SurrealFinite, v) ) == floor( v )
        @test ceil( convert(SurrealFinite, v) ) ≅ ceil( v )
        @test round( convert(SurrealFinite, v) ) ≅ round( v, RoundNearestTiesUp )
        @test trunc( convert(SurrealFinite, v) ) ≅ trunc( v )
        @test isinteger( convert(SurrealFinite, v) ) == isinteger( v )
    end
    @test_throws ErrorException floor( convert(SurrealFinite, 10000) )
    @test_throws ErrorException floor( convert(SurrealFinite, -10000) )

    # @test floor2(zero(SurrealFinite)) ≅ 0.0
    # @test floor2(one(SurrealFinite)) ≅ 1.0
    # @test floor2(x22) ≅ x22

    @test ceil(zero(SurrealFinite)) ≅ 0.0
    @test ceil(one(SurrealFinite)) ≅ 1.0
    @test ceil(x22) ≅ x22
    @test ceil(x4) ≅ zero(x4)
    @test ceil(x4) ≇ x4
    @test ceil(-x4) ≇ -x4
    @test ceil(-x4) == one(x4)

    @test mod(zero(SurrealFinite), x21) == 0
    @test mod(one(SurrealFinite), x21) == 1
    @test mod(x22, x21) ≅ 0.0
    @test mod(x41, x21) ≅ zero(x4)
    @test mod(x3, x43) ≅ 3
    @test mod(-x43, x3) ≅ 2
    @test mod(s1 + 2, x21) ≅ s1
    @test_throws ErrorException mod( x3, s2 )

    @test isinteger(x0) == true
    @test isinteger(x22) == true
    @test isinteger(s1) == false
    @test isinteger(x43) == true
    @test isinteger(x41) == true
    @test isinteger( convert(SurrealFinite,-4) ) == true
    @test isinteger( convert(SurrealFinite,-1//4) ) == false

    @test iseven(x0) == true
    @test iseven(x22) == true
    @test iseven(s1) == false
    @test iseven(x3) == false
    @test isodd(s1) == false
    @test isodd(x3) == true
    @test iseven( convert(SurrealFinite,-4) ) == true
    @test iseven( convert(SurrealFinite,-1) ) == false

    @test isinf(s1) == false
    @test isnan(s1) == false
    @test isfinite(s1)
end

@testset "simplifications" begin
    @test canonicalise(x22) ≅ x22
    @test canonicalise(x43) ≅ convert(SurrealFinite,4)
    @test canonicalise(x43) != convert(SurrealFinite,2)*convert(SurrealFinite,2)
    @test canonicalise(x43) ≅ convert(SurrealFinite,2)*convert(SurrealFinite,2)
    @test !iscanonical(x22)
    @test iscanonical(x21)
end

g3 = SurrealFinite( [ 7//16 ], [ 15//16 ] )
g4 = SurrealFinite( [convert(SurrealFinite,3)], [convert(SurrealFinite,17)] )
# g4 ≅ convert(SurrealFinite,4)

@testset "conversion from surreal back to real" begin
    @test convert(Int64, convert(SurrealFinite, -2)) == -2
    @test convert(Int32, convert(SurrealFinite, 6)) == 6
    @test convert(Int, convert(SurrealFinite, 0)) == 0
    # @test Int64( convert(SurrealFinite, 0) ) == 0
    @test convert(Int64, convert(SurrealFinite, -2)) == -2
    @test convert(UInt32, convert(SurrealFinite, 6)) == 6
    @test convert(Int8, convert(SurrealFinite, 0)) == 0
    @test_throws InexactError convert(UInt8, convert(SurrealFinite, -1))

    @test convert(Rational, convert(SurrealFinite, 1//8)) == 1//8
    @test convert(Rational, convert(SurrealFinite, 9//8)) == 9//8
    @test convert(Rational, convert(SurrealFinite, -9//8)) == -9//8
    @test convert(Rational, convert(SurrealFinite, -27//16) ) == -27//16
    @test convert(Rational, x4) == -1//2
    @test convert(Rational, convert(SurrealFinite,2)*convert(SurrealFinite,2) ) == 4

    @test g3 ≅ convert(SurrealFinite, 0.5)
    @test convert(Rational, SurrealFinite( [ 7//16 ], [ 15//16 ] ) ) == 1//2
    @test convert(Rational{Int32}, SurrealFinite( [ 7//16 ], [ 15//16 ] ) ) == 1//2

    @test convert(Float64, x43 ) == 4.0
    @test convert(Float64, x4) == -0.5
    @test convert(AbstractFloat, x43 ) == 4.0
    @test convert(AbstractFloat, x4) == -0.5
    @test float( x4 ) == -0.5
    @test float( s2 ) == 0.75
    @test float( g4 ) == 4.0

    @test all(convert.(SurrealFinite, [-1, 0, 1, 2] ) == [ convert(SurrealFinite, i) for i=-1:2 ])

    @test convert(String, x1) == "1"
    @test convert(String, x4) == "-1/2"
end

@testset "promotion" begin
    # lots of this has already been implicitly tested
    # note promoted type will be in canonical form, but other side might not
    @test x0 <= 1.0
    @test x1 == 1//1
    @test 1.0 + x1 == 2.0
    @test canonicalise(1//2 + x1) == 1.5
    @test 1//2 + x1 ≅ 1.5
    @test 1//2 + x1 ≇ 2.5
end

out_dir = "Data/"
a4 = convert(SurrealFinite, 13//16)
a5 = convert(SurrealFinite, -2//16)
s = "{ {{{{{|0}|{0|}}|{{0|}|}}|{{{0|}|}|}}|{{{{0|}|}|}|}} | }"    # from http://goodmath.scientopia.org/2006/08/21/arithmetic-with-surreal-numbers/
ss = convert(SurrealFinite, s)
@testset "I/O and string parsing" begin
    @test_throws Meta.ParseError convert(SurrealFinite, "{1//2| \\phi, 1}")
    @test convert(SurrealFinite, "{1//2 | 2,  1}") == SurrealFinite( [1//2], [2,1] )
    @test convert(SurrealFinite, raw"{-1//2 | {|\phi}}") ==  SurrealFinite( [-1//2], [zero(SurrealFinite)] )
    @test convert(SurrealFinite, "{-1//2 | 2.0, 3}") ==  SurrealFinite( [-1//2], [2,3] )
    R = read("$(out_dir)test_surreals.dat", SurrealFinite, 4)
    @test all( float(R[1:3]) .== [0.0, -1.0, 0.0] )
    @test_throws UndefRefError R[4]

    io_test_file = "$(out_dir)test_surreals_io.dat"
    io = open(io_test_file, "w")
    pf(io, a4)
    pf(io, a5)
    println(io, expand(a4))
    println(io, expand(a5))
    close(io)
    YY = read(io_test_file, SurrealFinite, 2)
    @test YY[1] == a4
    @test YY[2] == a5

    io_test_file = "$(out_dir)test_surreals_io.tex"
    io = open(io_test_file, "w")
    surreal2tex(io, a4; level = 0)
    surreal2tex(io, a4; level = 1)
    surreal2tex(io, a4; level = 2)
    surreal2tex(a4; level = 2)
    close(io)
    # should compare this to a calibration file
    #  but not sure if this might introduce potential for system dependencies that aren't real errors

    @test ss ≅ 4.0
    @test ss == convert(SurrealShort, 2.0) * convert(SurrealDyadic, 2.0)
end

@testset "structured output" begin
    @test surreal2dot(devnull, x0) == 1
    @test surreal2dot(x0) == 1
    @test surreal2dot(devnull, x23) == 5
    # @test surreal2dot(devnull, x43) == tree_nodes(x43)
    @test surreal2dag(devnull, x23) == 4
    # @test surreal2dag(devnull, x43) == nodes(x43)
    # @test surreal2dag(devnull, -x43) == nodes(x43)
    @test surreal2dag(devnull, -s2) == nodes(s2)

    # should compare these against a reference, but sort could get non-unique order, so potential for difference
end

a24 = dali(2) * dali(4);
(P,Q,U,V) = uniqueness(a24);
V2 = uniqueness_failure( U, V)
U2 = filter( (k,v) -> v>1 , U)
P2 = filter( (k,v) -> v==0xa09c3927cec6b030, P )

x32 = convert(SurrealFinite, 3/4) + convert(SurrealFinite, 3/4)
@testset "DAG statistics" begin
    @test dag_stats(x1) == SurrealDAGstats(2,2,1,1,ϕ,1,1,0,1)
    @test dag_stats(x41) == SurrealDAGstats(5,5,4,4,ϕ,1,4,0,4)
    @test dag_stats(x43) == SurrealDAGstats(11,21,14,6,ϕ,5,4,-1,4)
    @test dag_stats(x00) == SurrealDAGstats(4,5,4,2,ϕ,2,0,-1,1)
    @test dag_stats(x5) == SurrealDAGstats(5,12,9,3,ϕ,6,1//2,-1,1)
    @test dag_stats(s2) == SurrealDAGstats(4,7,5,3,ϕ,3,3//4,0,1)
    @test dag_stats(x6) == SurrealDAGstats(45,2574,82,12,ϕ,625,6,-3,7)
    @test dag_stats(x1, Dict{SurrealFinite,SurrealDAGstats}())[2] == Dict(x0=>SurrealDAGstats(1,1,0,0,ϕ,1,0,0,0),
                                                                          x1=>SurrealDAGstats(2,2,1,1,ϕ,1,1,0,1)    )
    @test dag_stats(x43).nodes == 11
    @test dag_stats(x43).nodes == nodes(x43)
    @test dag_stats(x43).tree_nodes == 21
    @test dag_stats(x43).tree_nodes == tree_nodes(x43)
    @test dag_stats(x43).edges == 14
    @test dag_stats(x43).edges == edges(x43)
    @test dag_stats(x43).generation == 6
    @test dag_stats(x43).paths == 5
    @test dag_stats(x43).value == 4
    @test dag_stats(x43).minval == -1
    @test dag_stats(x43).maxval == 4
    @test generation(x0) == 0

    @test generation(x1) == 1
    @test generation(x21) == 2
    @test generation(x22) == 2
    @test generation(x4) == 2
    @test generation(x5) == 3
    @test generation( x43 ) == 6
    @test generation( convert(SurrealFinite,4) ) == 4

    @test nodes(x0) == 1
    @test nodes(x1) == 2
    @test nodes(x41) == 5
    @test nodes(x43) == 11
    @test nodes(x5) == 5
    @test nodes(s2) == 4

    @test breadth(x1) == 1
    @test breadth(x43) == 5

    @test paths(x43) == 5

    @test dag_stats(x32).generation == 6
    @test dag_stats(x32).generation == generation(x32)

    # @test dag_stats(x41; LP=true) == SurrealDAGstats(5, 4, 4, [0, 1, 2, 3, 4], 1,4,0,4)
    @test dag_stats(x41; LP=true).longest_path == convert.(SurrealFinite, [0, 1, 2, 3, 4])

    @test uniqueness_max(U) == 1 # this test causes an error if we do say 2*4, without the caching in +
end

d33 = dali(33//2)
@testset "memory allocation and caching" begin
    # @test Base.summarysize(  ) == 943
    # 1775 was the old size before we made "dali" reuse cached numbers
    # 1743 comes up in v0.7, instead of the number above, not sure why versions are so different,
    # but if it is implementation dependent it shouldn't be a test

    # depends on previous results -- could be changed by earlier operations
    @test SurrealNumbers.ExistingSurreals[ SurrealNumbers.ExistingCanonicals[-1//2] ] == -1//2
    @test SurrealNumbers.ExistingSurreals[ SurrealNumbers.ExistingProducts[0x013b60c22ac142fa][0x84037d1bb0afff04] ] == -1
    clearcache()
    @test isempty(SurrealNumbers.ExistingSurreals)
    @test isempty(SurrealNumbers.ExistingCanonicals)
    @test isempty(SurrealNumbers.ExistingProducts)
    @test isempty(SurrealNumbers.ExistingSums)
    @test isempty(SurrealNumbers.ExistingNegations)
    Count['+'] == 0
    Count['*'] == 0

    m = dali(2) * dali(2)
    Count['+'] == 4
    Count['+'] == 13
    Count['*'] == 6
    Count['c'] == 20
end

@testset "new identity test" begin
    c1 = dali(1) + dali(3)
    c2 = dali(2) + dali(2)
    @test c1 == c2
    @test c1 === c2
end


text:
import Base.Threads: @spawn, @threads

function run_simulations(r::Rheology,p::Params,D₀::Vector,σ₂₂::Real,ϵ̇₁₁::Real,tspan ; multithreaded=true)
  # Vector variable
  vec_var = D₀

  # initial higher compressive stress multiplier (S = σ₁₁/σ₂₂)
  S = 1.0

  # Create a Rheology type instance containing elastic moduli and damage parameters. Change default values by  supplying keywords arguments.
  r_vec = [Rheology(r,(D₀=D0,)) for D0 in D₀]

  # time integration parameters :
  Δt_i = 0.0001 # initial timestep

  # initialize containers
  t_vec_vec = Vector{Vector{Float64}}(undef,length(vec_var))
  σᵢⱼ_vec_vec = Vector{Vector{SymmetricTensor{2,3}}}(undef,length(vec_var))
  ϵᵢⱼ_vec_vec = Vector{Vector{SymmetricTensor{2,3}}}(undef,length(vec_var))
  D_vec_vec =  Vector{Vector{Float64}}(undef,length(vec_var))

  # Build initial stress and strain tensors
  σᵢⱼ_i_vec = [build_principal_stress_tensor(r_vec[i],S,σ₂₂,D₀[i] ; abstol=1e-15) for i in eachindex(vec_var)]
  ϵᵢⱼ_i_vec = [compute_ϵij(r_vec[i],D₀[i],σᵢⱼ_i_vec[i]) for i in eachindex(vec_var)]

  if multithreaded
    @threads for i in eachindex(vec_var)
      # solve
      t_vec_vec[i], σᵢⱼ_vec_vec[i], ϵᵢⱼ_vec_vec[i], D_vec_vec[i] = adaptative_time_integration(r_vec[i],p,σᵢⱼ_i_vec[i],ϵᵢⱼ_i_vec[i],D₀[i],ϵ̇₁₁,Δt_i,tspan)

    end
  else
    for i in eachindex(vec_var)
      # solve
      t_vec_vec[i], σᵢⱼ_vec_vec[i], ϵᵢⱼ_vec_vec[i], D_vec_vec[i] = adaptative_time_integration(r_vec[i],p,σᵢⱼ_i_vec[i],ϵᵢⱼ_i_vec[i],D₀[i],ϵ̇₁₁,Δt_i,tspan)

    end
  end
  return t_vec_vec, σᵢⱼ_vec_vec, ϵᵢⱼ_vec_vec, D_vec_vec, r_vec
end

function run_simulations(r::Rheology,p::Params,D₀::Real,σ₂₂::Vector,ϵ̇₁₁::Real,tspan ; multithreaded=true)
  # Vector variable
  vec_var = σ₂₂

  # initial higher compressive stress multiplier (S = σ₁₁/σ₂₂)
  S = 1.0

  # Create a Rheology type instance containing elastic moduli and damage parameters. Change default values by  supplying keywords arguments.
  r_vec = [Rheology(r,(D₀=D₀,)) for i in eachindex(vec_var)]

  # time integration parameters :
  Δt_i = 0.0001 # initial timestep

  # initialize containers
  t_vec_vec = Vector{Vector{Float64}}(undef,length(vec_var))
  σᵢⱼ_vec_vec = Vector{Vector{SymmetricTensor{2,3}}}(undef,length(vec_var))
  ϵᵢⱼ_vec_vec = Vector{Vector{SymmetricTensor{2,3}}}(undef,length(vec_var))
  D_vec_vec =  Vector{Vector{Float64}}(undef,length(vec_var))

  # Build initial stress and strain tensors
  σᵢⱼ_i_vec = [build_principal_stress_tensor(r_vec[i],S,σ₂₂[i],D₀ ; abstol=1e-15) for i in eachindex(vec_var)]
  ϵᵢⱼ_i_vec = [compute_ϵij(r_vec[i],D₀,σᵢⱼ_i_vec[i]) for i in eachindex(vec_var)]

  if multithreaded
    @threads for i in eachindex(vec_var)
      # solve
      t_vec_vec[i], σᵢⱼ_vec_vec[i], ϵᵢⱼ_vec_vec[i], D_vec_vec[i] = adaptative_time_integration(r_vec[i],p,σᵢⱼ_i_vec[i],ϵᵢⱼ_i_vec[i],D₀,ϵ̇₁₁,Δt_i,tspan)

    end
  else
    for i in eachindex(vec_var)
      # solve
      t_vec_vec[i], σᵢⱼ_vec_vec[i], ϵᵢⱼ_vec_vec[i], D_vec_vec[i] = adaptative_time_integration(r_vec[i],p,σᵢⱼ_i_vec[i],ϵᵢⱼ_i_vec[i],D₀,ϵ̇₁₁,Δt_i,tspan)

    end
  end

  return t_vec_vec, σᵢⱼ_vec_vec, ϵᵢⱼ_vec_vec, D_vec_vec, r_vec
end

function run_simulations(r::Vector,p::Params,D₀::Real,σ₂₂::Real,ϵ̇₁₁::Real,tspan ; multithreaded=true)
  # Vector variable
  vec_var = r

  # initial higher compressive stress multiplier (S = σ₁₁/σ₂₂)
  S = 1.0

  # Create a Rheology type instance containing elastic moduli and damage parameters. Change default values by  supplying keywords arguments.
  r_vec = [Rheology(r[i],(D₀=D₀,)) for i in eachindex(vec_var)]

  # time integration parameters :
  Δt_i = 0.0001 # initial timestep

  # initialize containers
  t_vec_vec = Vector{Vector{Float64}}(undef,length(vec_var))
  σᵢⱼ_vec_vec = Vector{Vector{SymmetricTensor{2,3}}}(undef,length(vec_var))
  ϵᵢⱼ_vec_vec = Vector{Vector{SymmetricTensor{2,3}}}(undef,length(vec_var))
  D_vec_vec =  Vector{Vector{Float64}}(undef,length(vec_var))

  # Build initial stress and strain tensors
  σᵢⱼ_i_vec = [build_principal_stress_tensor(r_vec[i],S,σ₂₂,D₀ ; abstol=1e-15) for i in eachindex(vec_var)]
  ϵᵢⱼ_i_vec = [compute_ϵij(r_vec[i],D₀,σᵢⱼ_i_vec[i]) for i in eachindex(vec_var)]

  if multithreaded
    @threads for i in eachindex(vec_var)
      # solve
      t_vec_vec[i], σᵢⱼ_vec_vec[i], ϵᵢⱼ_vec_vec[i], D_vec_vec[i] = adaptative_time_integration(r_vec[i],p,σᵢⱼ_i_vec[i],ϵᵢⱼ_i_vec[i],D₀,ϵ̇₁₁,Δt_i,tspan)

    end
  else
    for i in eachindex(vec_var)
      # solve
      t_vec_vec[i], σᵢⱼ_vec_vec[i], ϵᵢⱼ_vec_vec[i], D_vec_vec[i] = adaptative_time_integration(r_vec[i],p,σᵢⱼ_i_vec[i],ϵᵢⱼ_i_vec[i],D₀,ϵ̇₁₁,Δt_i,tspan)

    end
  end
  return t_vec_vec, σᵢⱼ_vec_vec, ϵᵢⱼ_vec_vec, D_vec_vec, r_vec
end


text:
module util

function load_input(day::Int, strip::Bool = true)::String
    s = open(joinpath(dirname(@__DIR__), "inputs/day_$day.txt")) do f
        read(f, String)
    end
    if strip
        Base.strip(s)
    else
        s
    end
end

function input_lines(day::Int, strip::Bool = true)::Array{String}
    lines = split(load_input(day, true), '\n')
    if strip
        map(Base.strip, lines)
    else
        lines
    end
end

end


text:
"""
    RunnerInput(suite, graderfile, submissionfile, outfile)

A bundle of data that is used by a runner process to execute a test suite on a
code submission.
* `suite` is the `TestSuite` to be executed.
* `submissionfile` is the filename of the submission; it should be an absolute
    path.
* `outfile` is the filename to output the `serialize`'d `RunnerOutput`'.
"""
struct RunnerInput
    suite::TestSuite
    subfile::String
    funcname::Symbol
    outfile::String
end

"""
    RunnerOuput(result)

A bundle of data that is outputted by the runner processes.
This is currently just a vessel for `TestSuiteResult` but we could feasibly
add other information to this bundle that doesn't necessarily belong on the
`result` itself (e.g. execution metadata).
"""
struct RunnerOutput
    result::TestSuiteResult
end

"""
    runner_main(input; main=true)

Execute the runner.
If `main` is true, writes the `serialize`'d `RunnerOutput` to the file
specified in the `input` argument and the process will exit when this method
returns.
This function is potentially unsafe (it loads probably untrusted code) and so
should **never** be run outside of a sandbox.

This is mean to be the main entrypoint for runner subprocesses.
"""
function runner_main(input::RunnerInput; exit=true, write=true)
    try
        func = loadfunctionfromfile(
            input.subfile,
            input.funcname,
            nothing,
        )
        @debug "Loaded submission function." name=input.funcname file=input.subfile
        suite_result = runtestsuite(func, input.suite)
        output = RunnerOutput(suite_result)
        if write
            @debug "Writing runner output to file: $(input.outfile)"
            open(io -> serialize(io, output), input.outfile, "w")
        end
        exit && Base.exit(0)
        return output
    catch e
        @error "Unhandled error in runner." exception=(e, catch_backtrace())
        exit ? Base.exit(1) : rethrow(e)
    end
    # We should never get here™.
    error("Unanticipated control flow.")
end

function runner_main()
    input_path = get(ARGS, 1, nothing)
    if input_path === nothing
        @error "The serialized input file must be specified as the first command line argument."
        exit(1)
    end

    input = nothing
    try
        input = open(io -> deserialize(io), input_path, "r")
        typeassert(input, RunnerInput)
    catch exc
        @error "Unable to deserialize input file." exception=exc
        exit(1)
    end

    runner_main(input)

    # The `runner_main` function should always exit.
    error("Unanticipated control flow.")
end


text:
function part1(input)
    return count(n -> n in (2, 3, 4, 7), length.(reduce(vcat, last.(parse_input(input)))))
end

function part2(input)
    return sum(decode.(parse_input(input)))
end

parse_input(input) = parse_line.(readlines(input))

parse_line(line) = split.(split(line, " | "))

function decode(data)
    in_data, out_data = data
    display_digits = ["abcefg", "cf", "acdeg", "acdfg", "bcdf",
                      "abdfg", "abdefg", "acf", "abcdefg", "abcdfg"]
    encoded_digits = encode.(display_digits, 1, 2, 3, 4, 5, 6, 7)
    sorted_digits = sort(encoded_digits)

    for a in 1:7, b in 1:7, c in 1:7, d in 1:7, e in 1:7, f in 1:7, g in 1:7
        length(unique([a, b, c, d, e, f, g])) == 7 || continue
        digits = Int[]
        for D in in_data
            n = encode(D, a, b, c, d, e, f, g)
            push!(digits, n)
        end
        if sort(digits) == sorted_digits
            v = 0
            for D in out_data
                v *= 10
                n = encode(D, a, b, c, d, e, f, g)
                v += findfirst(encoded_digits .== n) - 1
            end
            return  v
        end
    end
    @assert false
end

function encode(D, a, b, c, d, e, f, g)
    n = 0
    for C in D
        C == 'a' && (n |= 1 << a)
        C == 'b' && (n |= 1 << b)
        C == 'c' && (n |= 1 << c)
        C == 'd' && (n |= 1 << d)
        C == 'e' && (n |= 1 << e)
        C == 'f' && (n |= 1 << f)
        C == 'g' && (n |= 1 << g)
    end
    return n
end


text:
# This file is a part of KCenters.jl
# License is Apache 2.0: https://www.apache.org/licenses/LICENSE-2.0.txt

using SimilaritySearch
export fftraversal, enet
using Distributed

function _ignore3(a, b, c)
end

"""
    fftraversal(callback::Function, dist::PreMetric, X::AbstractVector{T}, stop, callbackdist=_ignore3) where {T}

Selects a number of farthest points in `X`, using a farthest first traversal

- The callback function is called on each selected far point as `callback(centerID, dmax)` where `dmax` is the distance to the nearest previous reported center (the first is reported with typemax)
- The selected objects are far under the `dist` distance function with signature (T, T) -> Float64
- The number of points is determined by the stop criterion function with signature (Float64[], T[]) -> Bool
    - The first argument corresponds to the list of known distances (far objects)
    - The second argument corresponds to the database
- Check `criterions.jl` for basic implementations of stop criterions
- The callbackdist function is called on each distance evaluation between pivots and items in the dataset
    `callbackdist(index-pivot, index-item, distance)`
"""
function fftraversal(callback::Function, dist::PreMetric, X::AbstractVector{T}, stop, callbackdist=_ignore3) where {T}
    N = length(X)
    D = Vector{Float64}(undef, N)
    dmaxlist = Float64[]
    dset = [typemax(Float64) for i in 1:N]
    imax::Int = rand(1:N)
    dmax::Float64 = typemax(Float64)
    if N == 0
        return
    end

    k::Int = 0

    @inbounds while k <= N
        k += 1
        pivot = X[imax]
        push!(dmaxlist, dmax)
        callback(imax, dmax)
        dmax = 0.0
        ipivot = imax
        imax = 0

        D .= 0.0

        if nworkers() == 1
            Threads.@threads for i in 1:N
                D[i] = evaluate(dist, X[i], pivot)
            end
        else
            # only worths for very large datasets or very expensive distance functions
            D = pmap(obj -> evaluate(dist, obj, pivot), X)
        end

        for i in 1:N
            # d = evaluate(dist, X[i], pivot)
            d = D[i]
            callbackdist(i, ipivot, d)

            if d < dset[i]
                dset[i] = d
            end

            if dset[i] > dmax
                dmax = dset[i]
                imax = i
            end
        end

        if dmax == 0 || stop(dmaxlist, X)
            break
        end
    end
end

"""
    enet(dist::PreMetric, X::AbstractVector{T}, numcenters::Int, knr::Int=1; verbose=false) where T

Selects `numcenters` far from each other based on Farthest First Traversal.

- `dist` distance function
- `X` the objects to be computed
- `numcenters` number of centers to be computed
- `knr` number of nearest references per object (knr=1 defines a partition)

Returns a named tuple \$(nn, irefs, dmax)\$.

- `irefs` contains the list of centers (indexes to ``X``)
- `seq` contains the ``k`` nearest references for each object in ``X`` (in ``X`` order)
- `dmax` smallest distance among centers

"""
function enet(dist::PreMetric, X::AbstractVector{T}, numcenters::Integer, knr::Integer=1; verbose=false) where T
    # refs = Vector{Float64}[]
    irefs = Int32[]
    nn = [KnnResult(knr) for i in 1:length(X)]
    dmax = zero(Float32)

    function callback(c, _dmax)
        push!(irefs, c)
        dmax = convert(Float32, _dmax)
        verbose && println(stderr, "computing fartest point $(length(irefs)), dmax: $dmax, imax: $c")
    end

    function capturenn(i, refID, d)
        push!(nn[i], refID, d)
    end

    fftraversal(callback, dist, X, size_criterion(numcenters), capturenn)
    return (irefs=irefs, seq=nn, dmax=dmax)
end


text:
function _RobustBertsimas_latex()
    return """
        ```math
        \\left\\{ \\mu \\; \\middle| \\begin{array}{ll}
        s.t.  \\quad \\mu_i \\leq \\hat{r}_i + z_i \\Delta_i \\quad \\forall i = 1:\\mathcal{N} \\\\
        \\quad \\quad \\mu_i \\geq \\hat{r}_i - z_i \\Delta_i  \\quad \\forall i = 1:\\mathcal{N} \\\\
        \\quad \\quad z_i \\geq 0 \\quad \\forall i = 1:\\mathcal{N} \\\\
        \\quad \\quad z_i \\leq 1 \\quad \\forall i = 1:\\mathcal{N} \\\\
        \\quad \\quad \\sum_{i}^{\\mathcal{N}} z_i \\leq \\Gamma \\\\
        \\end{array}
        \\right\\} \\\\
        ```
        """
end

"""
    RobustBertsimas <: AbstractMeanVariance

Bertsimas's uncertainty set:

$(_RobustBertsimas_latex())

Atributes:
- `predicted_mean::Array{Float64,1}` (latex notation ``\\hat{r}``): Predicted mean of returns.
- `uncertainty_delta::Array{Float64,1}` (latex notation ``\\Delta``): Uncertainty around mean.
- `bertsimas_budjet::Array{Float64,1}` (latex notation ``\\Gamma``): Number of assets in worst case.
- `predicted_covariance::Array{Float64,2}`: Predicted covariance of returns.
"""
struct RobustBertsimas <: AbstractMeanVariance
    predicted_mean::Array{Float64,1}
    predicted_covariance::Array{Float64,2}
    uncertainty_delta::Array{Float64,1}
    bertsimas_budjet::Float64
    number_of_assets::Int
end

function RobustBertsimas(;
    predicted_mean::Array{Float64,1},
    predicted_covariance::Array{Float64,2},
    uncertainty_delta::Array{Float64,1},
    bertsimas_budjet::Float64,
)
    number_of_assets = size(predicted_mean, 1)
    if number_of_assets != size(predicted_covariance, 1)
        error("size of predicted mean ($number_of_assets) different to size of predicted covariance ($(size(predicted_covariance, 1)))")
    end
    if number_of_assets != size(uncertainty_delta, 1)
        error("size of predicted mean ($number_of_assets) different to size of uncertainty deltas ($(size(uncertainty_delta,1)))")
    end

    @assert prod(uncertainty_delta .>= 0)
    @assert issymmetric(predicted_covariance)
    return RobustBertsimas(
        predicted_mean, predicted_covariance, uncertainty_delta, bertsimas_budjet, number_of_assets
    )
end

function _portfolio_return_latex_RobustBertsimas_primal()
    return """
        ```math
        \\begin{aligned}
        \\min_{\\mu, z} \\quad & \\mu ' w \\\\
        s.t. \\quad & \\mu_i \\leq \\hat{r}_i + z_i \\Delta_i \\quad \\forall i = 1:\\mathcal{N} \\quad &: \\pi^-_i \\\\
        & \\mu_i \\geq \\hat{r}_i - z_i \\Delta_i  \\quad \\forall i = 1:\\mathcal{N} \\quad &: \\pi^+_i \\\\
        & z_i \\geq 0 \\quad \\forall i = 1:\\mathcal{N} \\\\
        & z_i \\leq 1 \\quad \\forall i = 1:\\mathcal{N} \\quad &: \\theta_i \\\\
        & \\sum_{i}^{\\mathcal{N}} z_i \\leq \\Gamma \\quad : \\lambda \\\\
        \\end{aligned}
        ```
        """
end

function _portfolio_return_latex_RobustBertsimas_dual()
    return """
        ```math
        \\begin{aligned}
        \\max_{\\lambda, \\pi^-, \\pi^+, \\theta} \\quad  \\sum_{i}^{\\mathcal{N}} (\\hat{r}_i (\\pi^+_i - \\pi^-_i) - \\theta_i ) - \\Gamma \\lambda\\\\
        s.t. \\quad &  w_i = \\pi^+_i - \\pi^-_i  \\quad \\forall i = 1:\\mathcal{N} \\\\
        &  \\Delta_i (\\pi^+_i + \\pi^-_i) - \\theta_i \\leq \\lambda \\quad \\forall i = 1:\\mathcal{N} \\\\
        & \\lambda \\geq 0 , \\; \\pi^- \\geq 0 , \\; \\pi^+ \\geq 0 , \\; \\theta \\geq 0 \\\\
        \\end{aligned}
        ```
        """
end

"""
    portfolio_return!(model::JuMP.Model, w, formulation::RobustBertsimas)

Returns worst case return (WCR) in Bertsimas's uncertainty set ([`RobustBertsimas`](@ref)).

WCR is defined by the following primal problem:

$(_portfolio_return_latex_RobustBertsimas_primal())

Which is equivalent to the following dual problem:

$(_portfolio_return_latex_RobustBertsimas_dual())

To avoid solving an optimization problem we enforece the dual constraints in
the upper level problem and return the objective expression (a lower bound of the optimum).

Arguments:
 - `model::JuMP.Model`: JuMP upper level portfolio optimization model.
 - `w`: portfolio optimization investment variable ("weights").
 - `formulation::RobustBertsimas`: Struct containing atributes of Bertsimas's uncertainty set.
"""
function portfolio_return!(model::JuMP.Model, w, formulation::RobustBertsimas)
    # parameters
    r̄ = formulation.predicted_mean
    numA = formulation.number_of_assets
    Δ = formulation.uncertainty_delta
    Λ = formulation.bertsimas_budjet
    # dual variables
    @variable(model, λ >= 0)
    @variable(model, π_neg[i=1:numA] >= 0)
    @variable(model, π_pos[i=1:numA] >= 0)
    @variable(model, θ[i=1:numA] >= 0)
    # constraints: from duality
    @constraints(
        model,
        begin
            constrain_dual1[i=1:numA], w[i] == π_pos[i] - π_neg[i]
        end
    )
    @constraints(
        model,
        begin
            constrain_dual2[i=1:numA], Δ[i] * (π_pos[i] + π_neg[i]) - θ[i] <= λ
        end
    )

    return sum(r̄[i] * (π_pos[i] - π_neg[i]) for i in 1:numA) - sum(θ[i] for i in 1:numA) - λ * Λ
end

##################################
function _BenTal_latex()
    return """
        ```math
        \\left\\{ \\mu \\; \\middle| \\begin{array}{ll}
        s.t.  \\quad \\sqrt{(\\hat{r} - \\mu) ' \\Sigma^{-1} (\\hat{r} - \\mu)} \\leq \\delta \\\\
        \\end{array}
        \\right\\} \\\\
        ```
        """
end
"""
    BenTal <: AbstractMeanVariance

BenTal's uncertainty set:

$(_BenTal_latex())

Atributes:
- `predicted_mean::Array{Float64,1}` (latex notation ``\\hat{r}``): Predicted mean of returns.
- `uncertainty_delta::Float64` (latex notation ``\\delta``): Uncertainty around mean.
- `predicted_covariance::Array{Float64,2}` (latex notation ``\\Sigma``): Predicted covariance of returns.
"""
struct RobustBenTal <: AbstractMeanVariance
    predicted_mean::Array{Float64,1}
    predicted_covariance::Array{Float64,2}
    uncertainty_delta::Float64
    number_of_assets::Int
end

function RobustBenTal(;
    predicted_mean::Array{Float64,1},
    predicted_covariance::Array{Float64,2},
    uncertainty_delta::Float64,
)
    number_of_assets = size(predicted_mean, 1)
    if number_of_assets != size(predicted_covariance, 1)
        error("size of predicted mean ($number_of_assets) different to size of predicted covariance ($(size(predicted_covariance, 1)))")
    end

    @assert uncertainty_delta >= 0
    @assert issymmetric(predicted_covariance)
    return RobustBenTal(
        predicted_mean, predicted_covariance, uncertainty_delta, number_of_assets
    )
end

function _portfolio_return_latex_RobustBenTal_primal()
    return """
        ```math
        \\begin{aligned}
        \\min_{\\mu, z} \\quad & \\mu ' w \\\\
        s.t. \\quad & ||Σ^{-\\frac{1}{2}} (\\mu - \\hat{r}) || \\leq \\delta \\quad &: \\theta \\\\
        \\end{aligned}
        ```
        """
end

function _portfolio_return_latex_RobustBenTal_dual()
    return """
        ```math
        \\begin{aligned}
        \\max_{\\theta} \\quad &  w ' \\hat{r} - \\theta \\delta \\\\
        s.t. \\quad & ||Σ^{\\frac{1}{2}}  w || \\leq \\theta \\\\
        \\end{aligned}
        ```
        """
end

"""
    portfolio_return!(model::JuMP.Model, w, formulation::RobustBenTal)

Returns worst case return (WCR) in BenTal's uncertainty set ([`RobustBenTal`](@ref)).

WCR is defined by the following primal problem:

$(_portfolio_return_latex_RobustBenTal_primal())

Which is equivalent to the following dual problem:

$(_portfolio_return_latex_RobustBenTal_dual())

To avoid solving an optimization problem we enforece the dual constraints in
the upper level problem and return the objective expression (a lower bound of the optimum).

Arguments:
    - `model::JuMP.Model`: JuMP upper level portfolio optimization model.
    - `w`: portfolio optimization investment variable ("weights").
    - `formulation::RobustBenTal`: Struct containing atributes of BenTal's uncertainty set.
"""
function portfolio_return!(model::JuMP.Model, w, formulation::RobustBenTal)
    # parameters
    r̄ = formulation.predicted_mean
    numA = formulation.number_of_assets
    δ = formulation.uncertainty_delta
    sqrt_Σ = sqrt(formulation.predicted_covariance)
    # dual variables
    @variable(model, θ)
    # constraints: from duality
    @constraint(model, [θ; sqrt_Σ * w] in JuMP.SecondOrderCone())

    return dot(w, r̄) - θ * δ
end


text:
using Test, DiffEqOperators
@testset "FD DO operators comparison" begin
    n = 100
    x=0.0:0.01:2π
    dx=diff(x)
    y = sin.(x);

    for dor in 1:4, aor in 2:2:8,
            LBC in (:None,:Dirichlet0, :Dirichlet, :Neumann0, :Neumann, :periodic),
            RBC in (:None,:Dirichlet0, :Dirichlet, :Neumann0, :periodic)

        D1 = DerivativeOperator{Float64}(dor,aor,dx[1],length(x),:None,:None)
        D2 = DiffEqOperators.FiniteDifference{Float64}(dor,aor,dx,length(x),:None,:None)

        @test convert(Array, D1) ≈ convert(Array, D2)
    end

end


text:
mutable struct P_bandpass{T<:Real}
	dfreq::Array{Complex{T}}
	dfftp::FFTW.rFFTWPlan
	difftp::FFTW.Plan
	filter::Array{T}
end

function P_bandpass(T=Float64;
		    fmin, fmax, order, length)

	l=length
	nrfft=div(l,2)+1
	dfftp=plan_rfft(zeros(T, l),[1])
	difftp=plan_irfft(complex.(zeros(T, nrfft)),l,[1])

	dfreq=complex.(zeros(T,nrfft))
	filter=zeros(T,nrfft)

	filter=abs.(rfft(impz(b)))
end

function filt!(y,x,pa::P_bandpass)
	mul!(pa.dfreq, pa.dfftp, x)
	mul!(y, pa.difftp, pa.dfreq)
end


text:
# Unit tests for permutations.jl
# Format of tests: TBD

module TestPermutations
    using MarkovTools
    using FactCheck # a different unit test framework
    using Base.Test

    # count_fixed_points
    p1 = []
    @test isequal(count_fixed_points(p1), 0)

    p1 = [1, 2, 3, 4, 5]
    @test isequal(count_fixed_points(p1), 5)

    p1 = [5, 3, 1, 2, 4]
    @test isequal(count_fixed_points(p1), 0)

    p1 = [5, 4, 3, 2, 1]
    @test isequal(count_fixed_points(p1), 1)

    N = 1e6
    p1 = [1:N]
    @test isequal(count_fixed_points(p1), N)

    facts("Testing count_fixed_points") do
        @fact count_fixed_points([1, 3, 2, 4, 6, 5]) => 2
        @fact count_fixed_points([]) => 0
    end

    # count_cycles
    p1 = []
    @test isequal(count_cycles(p1), 0)

    p1 = [1, 2, 3, 4, 5]
    @test isequal(count_cycles(p1), 5)

    p1 = [2, 1, 3, 5, 4]
    @test isequal(count_cycles(p1), 3)

    facts("Testing count_cycles") do
        @fact count_cycles([2, 3, 4, 5, 6, 7, 8, 1]) => 1
        @fact count_fixed_points([1, 2, 3, 4, 5, 6, 7, 8]) => 8
    end

end

max_stars_repo_name: PawelPamula/MarkovTools

max_issues_repo_name: PawelPamula/MarkovTools

max_forks_repo_name: PawelPamula/MarkovTools


text:
module NMF_VB

include("../Misc.jl")
using .Misc

using Distributions, SpecialFunctions

export standard_VB, online_VB

function standard_VB(X::Array{ℜ,2}, K::Ƶ ; μ::ℜ=nanmean(X), γ::ℜ=0.1, EPOCHS::Ƶ=1, ϵ::ℜ=1e-16) where {ℜ<:Real, Ƶ<:Int}
    I::Ƶ, J::Ƶ = size(X)

    a::ℜ, b::ℜ = I*J*γ, γ/μ
    α_K, α_JK, α_IK = fill(a/K,K), fill(a/(K*J),J,K), fill(a/(K*I),I,K)
    S₊ = sum(X)

    #log_λ = log(rand(Gamma(a,1.0/b)))
    log_λ = digamma(S₊ + a) - log(b + 1)
    log_θ_K = log.(rand(Dirichlet(α_K.+0.1)))
    log_θ_JK = reshape([log(θ_jk) for k=1:K for θ_jk=rand(Dirichlet(α_JK[:,k].+0.1))],J,K)
    log_θ_IK = reshape([log(θ_ik) for k=1:K for θ_ik=rand(Dirichlet(α_IK[:,k].+0.1))],I,K)

    S_K, S_JK, S_IK, S⁺::ℜ = zeros(ℜ,K), zeros(ℜ,J,K), zeros(ℜ,I,K), 0.0
    s, p, log_ρ, log_ρ₊ = Array{ℜ}(undef,K), Array{ℜ}(undef,K), Array{ℜ}(undef,K), 0.0

    ELBO::Array{ℜ} = zeros(ℜ,EPOCHS)
    X_full = similar(X)

    for eph=1:EPOCHS
        S_K .= 0.0
        S_JK .= 0.0
        S_IK .= 0.0
        S⁺ = 0.0

        for j=1:J, i=1:I #order of traversal is important
            log_ρ .= log_λ .+ log_θ_K .+ log_θ_JK[j,:] .+ log_θ_IK[i,:]
            log_ρ₊ = logsumexp(log_ρ)
            p .= exp.(log_ρ .- log_ρ₊)

            X_full[i,j] = isnan(X[i,j]) ? exp(log_ρ₊) : X[i,j]
            s .= X_full[i,j] .* p

            S_JK[j,:] .+= s
            S_IK[i,:] .+= s
            S_K[:] .+= s
            S⁺ += X_full[i,j]

            ELBO[eph] += isnan(X[i,j]) ? X_full[i,j] : X[i,j]*log_ρ₊ - lgamma(X[i,j] + 1.0)
        end
        ELBO[eph] += a*log(b) - (S⁺+a)*log(b+1.0) - sum(lgamma.(α_JK)) + sum(lgamma.(α_JK .+ S_JK))
        ELBO[eph] += sum(lgamma.(α_K)) - sum(lgamma.(α_IK)) - sum(lgamma.(α_K .+ S_K)) + sum(lgamma.(α_IK .+ S_IK))
        ELBO[eph] -= S⁺ * log_λ + sum(S_K .* log_θ_K) + sum(S_JK .* log_θ_JK) + sum(S_IK .* log_θ_IK)

        log_λ = digamma(S⁺+a) - log(b+1.0)
        log_θ_K .= digamma.(S_K.+α_K) .- digamma(S⁺+a)
        log_θ_JK .= digamma.(S_JK.+α_JK) .- digamma.(S_K.+α_K)'
        log_θ_IK .= digamma.(S_IK.+α_IK) .- digamma.(S_K.+α_K)'

    end

    return ELBO, X_full, log_λ, (log_θ_K, log_θ_IK, log_θ_JK)
end

function online_VB(X::Array{ℜ,2}, K::Ƶ, Ω::AbstractArray{ℜ,1}; μ::ℜ=nanmean(X), γ::ℜ=0.1, ϵ::ℜ=1e-16) where {ℜ<:Real, Ƶ<:Int}
    I::Ƶ, J::Ƶ = size(X)

    a::ℜ, b::ℜ = I*J*γ, γ/μ
    α_K, α_JK, α_IK = fill(a/K,K), fill(a/(K*J),J,K), fill(a/(K*I),I,K)

    log_λ = log(rand(Gamma(a,1.0/b)))
    log_θ_K = log.(rand(Dirichlet(α_K)))
    log_θ_JK = reshape([log(θ_jk) for k=1:K for θ_jk=rand(Dirichlet(α_JK[:,k]))],J,K)
    log_θ_IK = reshape([log(θ_ik) for k=1:K for θ_ik=rand(Dirichlet(α_IK[:,k]))],I,K)

    S_K, S_JK, S_IK, S⁺::ℜ = zeros(ℜ,K), zeros(ℜ,J,K), zeros(ℜ,I,K), 0.0
    s, p, log_ρ, log_ρ₊ = Array{ℜ}(undef,K), Array{ℜ}(undef,K), Array{ℜ}(undef,K), 0.0

    ELBO::Array{ℜ} = zeros(ℜ,length(Ω))
    X_full = similar(X)

    for (eph,ω) in enumerate(Ω)
        S_K .= 0.0
        S_JK .= 0.0
        S_IK .= 0.0
        S⁺ = 0.0

        for j=1:J, i=1:I #order of traversal is important
            log_ρ .= log_λ .+ log_θ_K .+ log_θ_JK[j,:] .+ log_θ_IK[i,:]
            log_ρ₊ = logsumexp(log_ρ)
            p .= exp.(log_ρ .- log_ρ₊)

            X_full[i,j] = isnan(X[i,j]) ? ω*exp(log_ρ₊) : rand(Binomial(Ƶ(X[i,j]),ω))
            s .= X_full[i,j] .* p

            S_JK[j,:] .+= s
            S_IK[i,:] .+= s
            S_K[:] .+= s
            S⁺ += X_full[i,j]

            ELBO[eph] += isnan(X[i,j]) ? X_full[i,j] : X[i,j]*log_ρ₊ - lgamma(X[i,j] + 1.0)
        end

        ELBO[eph] += a*log(b) - (S⁺+a)*log(b+1.0) - sum(lgamma.(α_JK)) + sum(lgamma.(α_JK .+ S_JK))
        ELBO[eph] += sum(lgamma.(α_K)) - sum(lgamma.(α_IK)) - sum(lgamma.(α_K .+ S_K)) + sum(lgamma.(α_IK .+ S_IK))
        ELBO[eph] -= S⁺ * log_λ + sum(S_K .* log_θ_K) + sum(S_JK .* log_θ_JK) + sum(S_IK .* log_θ_IK)

        log_λ = digamma(S⁺ + a) - log(b + ω) ## should be checked
        log_θ_K .= digamma.(S_K.+α_K) .- digamma(S⁺+a)
        log_θ_JK .= digamma.(S_JK.+α_JK) .- digamma.(S_K.+α_K)'
        log_θ_IK .= digamma.(S_IK.+α_IK) .- digamma.(S_K.+α_K)'

    end
    return ELBO, X_full, log_λ, (log_θ_K, log_θ_IK, log_θ_JK)
end

end

max_stars_repo_name: mehmetburakkurutmaz/bam

max_issues_repo_name: mehmetburakkurutmaz/bam

max_forks_repo_name: atcemgil/bam


text:
### Basic operators for any heat conduction system ###
# Note that the input vector `T` is temperature, and
# we expect the rhs of the equations to have units of dT/dt

export heatconduction_rhs!

# RHS of heat conduction equations
function heatconduction_rhs!(dT::Nodes{Primal,NX,NY},T::Nodes{Primal,NX,NY},sys::HeatConduction{NX,NY},t::Real) where {NX,NY}
  dT .= 0.0
  #_heatconduction_rhs_convectivederivative!(sys.Sc,T,sys,t)
  #_heatconduction_rhs_double_layer!(dT,sys,t)
  _heatconduction_rhs_double_layer_adiabatic!(dT,T,sys,sys.bctype)
  _heatconduction_rhs_forcing!(dT,sys,t)
  _heatconduction_rhs_qflux!(dT,sys)
  _heatconduction_rhs_qmodel!(dT,T,sys)
  _heatconduction_rhs_qline!(dT,sys)
  return dT
end

_heatconduction_rhs_qflux!(dT::Nodes{Primal,NX,NY},sys::HeatConduction{NX,NY}) where {NX,NY} =
      _heatconduction_rhs_qflux!(dT,sys,sys.qflux)

function _heatconduction_rhs_qflux!(dT::Nodes{Primal,NX,NY},sys::HeatConduction{NX,NY},qflux) where {NX,NY}
  @unpack params = sys
  @unpack ρ, c, d = params
  fact = 1.0/(ρ*c)
  fact = isnothing(d) ? fact : fact/d
  for q in qflux
    dT .+= fact.*q()
  end
  dT
end

_heatconduction_rhs_qflux!(dT::Nodes{Primal,NX,NY},sys::HeatConduction{NX,NY},::Nothing) where {NX,NY} = dT

_heatconduction_rhs_qmodel!(dT::Nodes{Primal,NX,NY},T::Nodes{Primal,NX,NY},sys::HeatConduction{NX,NY}) where {NX,NY} =
    _heatconduction_rhs_qmodel!(dT,T,sys,sys.qhdT)

function _heatconduction_rhs_qmodel!(dT::Nodes{Primal,NX,NY},T::Nodes{Primal,NX,NY},sys::HeatConduction{NX,NY},qhdT) where {NX,NY}
  @unpack params = sys
  @unpack ρ, c, d = params
  fact = 1.0/(ρ*c)
  fact = isnothing(d) ? fact : fact/d
  for q in qhdT
    dT .+= fact.*q(T)
  end
  dT
end

_heatconduction_rhs_qmodel!(dT::Nodes{Primal,NX,NY},T::Nodes{Primal,NX,NY},sys::HeatConduction{NX,NY},::Nothing) where {NX,NY} = dT

_heatconduction_rhs_qline!(dT::Nodes{Primal,NX,NY},sys::HeatConduction{NX,NY}) where {NX,NY} =
    _heatconduction_rhs_qline!(dT,sys,sys.qline)

function _heatconduction_rhs_qline!(dT::Nodes{Primal,NX,NY},sys::HeatConduction{NX,NY},qline) where {NX,NY}
  @unpack qline, params = sys
  @unpack ρ, c, d = params
  fact = 1.0/(ρ*c)
  fact = isnothing(d) ? fact : fact/d
  for ql in qline
    ql.cache1 .= ql.R*ql.q
    dT .-= fact.*ql.cache1
  end
  dT
end

_heatconduction_rhs_qline!(dT::Nodes{Primal,NX,NY},sys::HeatConduction{NX,NY},::Nothing) where {NX,NY} = dT

_heatconduction_rhs_forcing!(dT::Nodes{Primal,NX,NY},sys::HeatConduction{NX,NY},t) where {NX,NY} = _heatconduction_rhs_forcing!(dT,sys.qforce,cellsize(sys),t)

_heatconduction_rhs_forcing!(dT,::Nothing,Δx,t) = dT

function _heatconduction_rhs_forcing!(dT,qforce::Vector{<:ModulatedField},Δx,t)
  for q in qforce
    dT .+= q(t) #*Δx
  end
  dT
end

@inline _heatconduction_rhs_double_layer!(dT::Nodes{Primal,NX,NY},sys::HeatConduction{NX,NY,N,MT,
                                  ExternalInternalProblem},t::Real) where {NX,NY,N,MT} = dT

@inline _heatconduction_rhs_double_layer!(dT::Nodes{Primal,NX,NY},sys::HeatConduction{NX,NY,0},t::Real) where {NX,NY} = dT

function _heatconduction_rhs_double_layer!(dT::Nodes{Primal,NX,NY},sys::HeatConduction{NX,NY,N,MT,SD},t::Real) where {NX,NY,N,MT,SD}
    @unpack params = sys
    @unpack α = params
    Δx⁻¹ = 1/cellsize(sys)
    fact = α*Δx⁻¹
    surface_temperature_jump!(sys.ΔTs,sys,t)
    fill!(sys.Sc,0.0)
    sys.dlc(sys.Sc,sys.ΔTs)
    sys.Sc .*= fact
    dT .-= sys.Sc
end

@inline _heatconduction_rhs_double_layer_adiabatic!(dT::Nodes{Primal,NX,NY},T::Nodes{Primal,NX,NY},sys::HeatConduction{NX,NY,0},bctype) where {NX,NY} = dT

@inline _heatconduction_rhs_double_layer_adiabatic!(dT::Nodes{Primal,NX,NY},T::Nodes{Primal,NX,NY},sys::HeatConduction{NX,NY,N},bctype) where {NX,NY,N} = dT

function _heatconduction_rhs_double_layer_adiabatic!(dT::Nodes{Primal,NX,NY},T::Nodes{Primal,NX,NY},sys::HeatConduction{NX,NY,N},::Type{AdiabaticBC}) where {NX,NY,N}
    @unpack params, ΔTs, Sc, Cc = sys
    @unpack k = params
    fill!(ΔTs,0.0)
    heatflux_constraint_op!(ΔTs,T,sys)
    ΔTs .= -(Cc\ΔTs)
    ΔTs ./= k
    fill!(Sc,0.0)
    heatflux_constraint_force!(Sc,ΔTs,sys)
    dT .-= Sc
end

#=
function _heatconduction_rhs_convectivederivative!(u::Edges{Primal,NX,NY},w::Nodes{Dual,NX,NY},sys::HeatConduction{NX,NY},t) where {NX,NY}
    Δx⁻¹ = 1/cellsize(sys)
    fill!(sys.Vv,0.0)
    velocity!(sys.Vv,w,sys,t)
    _unscaled_convective_derivative!(sys.Vv,sys)
    sys.Vv .*= Δx⁻¹
    u .-= sys.Vv
end

_unscaled_convective_derivative!(u::Edges{Primal,NX,NY},sys::HeatConduction{NX,NY}) where {NX,NY} =
      _unscaled_convective_derivative!(u,sys.Vtf,sys.DVf,sys.VDVf)

# Operates in-place on `u`, which comes in with the velocity field and
# returns the unscaled convective derivative
function _unscaled_convective_derivative!(u::Edges{Primal,NX,NY},
                                          Vtf::EdgeGradient{Primal,Dual,NX,NY},
                                          DVf::EdgeGradient{Primal,Dual,NX,NY},
                                          VDVf::EdgeGradient{Primal,Dual,NX,NY}) where {NX,NY}
    transpose!(Vtf,grid_interpolate!(DVf,u))
    DVf .= 0.0
    grad!(DVf,u)
    product!(VDVf,Vtf,DVf)
    u .= 0.0
    grid_interpolate!(u,VDVf)
    u
end
=#


text:
function learn_rare!(model::Model, prior::GaussianPrior, opt::Opt, data; monitor = nothing)
	λm = [λmin(data.counts, t) for t in 1:model.T-1]
	for t in 1:model.T-1
		@assert model.λ[t] >= λm[t]
	end
	xinit = vcat(model.x, model.λ)
	logP = []
	properties = [[] for a in monitor]
	function f(x::Vector, grad::Vector, vals::Vector, prop::Vector)
		model.x .= x[1:end-model.T+1]
		model.λ .= x[end-model.T+2:end]
		if length(grad) > 0
			grad .= gradient_Likelihood_rare(model, data)
			@inbounds prior(grad[1:end-model.T+1], x[1:end-model.T+1])
		end
		l = Likelihood_rare(model, data) + prior(x[1:end-model.T+1])
		push!(vals, l)
		if monitor != nothing
			for i in 1:length(monitor)
				push!(prop[i], monitor[i](model, prior, data)[:])
			end
		end
		return l
	end

	f1(x::Vector, grad::Vector) = f(x, grad, logP, properties)
	opt.max_objective = f1
	(maxf, xmax, ret) = optimize(opt, xinit)
	@inbounds model.x .= xmax[1:end-model.T+1]
	@inbounds model.λ .= xmax[end-model.T+2:end]
	return (model, ret, logP, properties)
end


text:
# Range: 1 : 100,000,000
# Time:  36.494s

# Written by Alex Vear in 2019
# Public domain.  No rights reserved.

startval = 1
endval = 100000000

sqrs = [x^2 for x in startval:endval]
tris = [x*(x+1)/2 for x in startval:endval]
println(intersect(sqrs, tris))

max_stars_repo_name: axvr/codedump

max_issues_repo_name: axvr/codedump

max_forks_repo_name: axvr/codedump


text:
## Pipeline file to build an alignment from a full length sequence
struct SimpleAlign
    Header::Vector{String}
    Seqs::Vector{String}
    function SimpleAlign(H::Vector{String},S::Vector{String})
        l = length(S[1])
        length(H) == length(S) || throw(DimensionMismatch("header and seqs have different lenghts"))
        for i in eachindex(S)
            length(S[i]) == l || throw(DimensionMismatch("not an alignment"))
        end
        new(H,S)
    end
end

SimpleAlign(d::SummaryAlign) = SimpleAlign(d.Header, d.Seqs)

function idx_insert(d::SummaryAlign, thr_ins)
    gapfreq =  vec(1.0 .- sum(reshape(d.Pi,(d.q-1,d.N)),dims=1))
    return findall(gapfreq .> thr_ins)
end

function insertify(d::SummaryAlign, idx)
    headers, seqs = d.Header,d.Seqs
    newseqs = similar(seqs)
    for i in eachindex(seqs)
        charseq = collect(seqs[i])
        for j in idx
            _c = charseq[j]
            charseq[j] = isletter(_c) ? lowercase(_c) : '.'
        end
        newseqs[i] = String(charseq)
    end
    SimpleAlign(headers,newseqs)
end

function writefasta(fou::String, d::Dict{String,String})
    FastaWriter(fou) do fw
        for (k,v) in d
            writeentry(fw,k,v)
        end
    end
end

function writefasta(fou::String,d::SimpleAlign)
    FastaWriter(fou) do fw
        h,s = d.Header,d.Seqs
        for i in eachindex(h)
            writeentry(fw,h[i],s[i])
        end
    end
end

"""
    function align_seed_mafft(
            fin::String,
            fou::String;
            nthreads::Integer = 8,
            thr_ins::Real = 1.0
        )
Read `fin` full-length seed in FASTA format and produces `fou.fasta` and `fou.ins` alignments.

"""
function align_seed_mafft(
    fin::String,
    fou::String;
    nthreads::Integer = 1,
    thr_ins::Real = 0.5,
)
    mktempdir() do tmpdir
        @info "Working in $tmpdir"
        linsiout = joinpath(tmpdir, "all_linsi.fasta")
        @info "### Aligning $fin with mafft-linsi. This may take a while ... "
        linsi(fin, linsiout, nthreads = nthreads)
        @info "Alignment done!"
        rescor = summary_align(linsiout)
        tmp_align = if thr_ins < 1.0 # make insertion by hand
            idxins = idx_insert(rescor, thr_ins)
            insertify(rescor, idxins)
        else
            SimpleAlign(rescor)
        end
        @info "Converting in FASTA (+ insertions) format"
        file_tmp_align = joinpath(tmpdir, "tmp_linsi.fasta")
        writefasta(file_tmp_align, tmp_align)
        clean_align = extract_ins(file_tmp_align)
        writefasta(fou * ".ins", clean_align)
        clean_fasta, L = extract_align(clean_align)
        @info "Converting in FASTA format"
        writefasta(fou * ".fasta", clean_fasta)
        @info "L = $L"
        #@info "### Computing the HMM and re-aligning the seed"
        ####### why do we need .hmm ? why do we re-align the seed? #######
        #hmmbuild(file_tmp_align, fou * "_tmp.hmm", nthreads = 1)
        #file_tmp_stk = joinpath(tmpdir, "tmp_linsi.stk")
        #hmmalign(fou * "_tmp.hmm", fin, file_tmp_stk)
        # here file_tmp_align -> fin. we should align the original sequences
        # hmmalign(fou * ".hmm", file_tmp_align, file_tmp_stk)
        #d = read_stockholm(file_tmp_stk)
        #clean_seqs, g = extract_align(d)
        #writefasta(fou * "_tmp.fasta", clean_seqs)
        @info "Done!"
    end
end

"""
    function align_seed_pfam(
            fseqs::String,
            fout::String;
            nthreads::Integer = 8,
        )
Read `fseqs` full-length seed downloaded from Pfam, in Stockholm format. It produces `fou.fasta` and `fou.ins` alignments from the output of `hmmbuild`.

"""

function align_seed_pfam(
    fseqs::String,
    fou::String;
    nthreads::Integer = 1,
)

    mktempdir() do tmpdir
        @info "Working in $tmpdir"
        @info "### Computing the HMM and get aligned seed (with and without insertions)"
        hmmbuild(fseqs, joinpath(tmpdir, "tmp_stk.out"), joinpath(tmpdir, "tmp_stk.hmm"), nthreads = nthreads)
        file_tmp_stk = joinpath(tmpdir, "tmp_stk.out")
        d = read_stockholm(file_tmp_stk)
        #tmp, g = extract_align(d)
        writefasta(joinpath(tmpdir, "tmp.fasta"), d)
        clean_ins = extract_ins(joinpath(tmpdir, "tmp.fasta"))
        clean_fasta, L = extract_align(clean_ins)
        writefasta(fou * ".ins", clean_ins)
        writefasta(fou * ".fasta", clean_fasta)
        @info "L = $L"
        @info "Done!"
    end

end


text:
relaxed_tols = (default_tol_relax = 100,)
insts = OrderedDict()
insts["minimal"] = [
    ((2, 2), nothing, relaxed_tols),
    ]
insts["various"] = [
    ((2, 2), nothing, relaxed_tols),
    ((1, 4), nothing, relaxed_tols),
    ((4, 1), nothing, relaxed_tols),
    ]
return (EntanglementAssisted, insts)


text:
using DynamicPPL
using DynamicPPL: getargs_dottilde, getargs_tilde, get_type, get_symbol

@testset "getargs_dottilde" begin
    # Some things that are not expressions.
    @test getargs_dottilde(:x) === nothing
    @test getargs_dottilde(1.0) === nothing
    @test getargs_dottilde([1.0, 2.0, 4.0]) === nothing

    # Some expressions.
    @test getargs_dottilde(:(x ~ Normal(μ, σ))) === nothing
    @test getargs_dottilde(:((.~)(x, Normal(μ, σ)))) == (:x, :(Normal(μ, σ)))
    @test getargs_dottilde(:((~).(x, Normal(μ, σ)))) == (:x, :(Normal(μ, σ)))
    @test getargs_dottilde(:(@. x ~ Normal(μ, σ))) === nothing
    @test getargs_dottilde(:(@. x ~ Normal(μ, $(Expr(:$, :(sqrt(v))))))) === nothing
    @test getargs_dottilde(:(@~ Normal.(μ, σ))) === nothing
end

@testset "getargs_tilde" begin
    # Some things that are not expressions.
    @test getargs_tilde(:x) === nothing
    @test getargs_tilde(1.0) === nothing
    @test getargs_tilde([1.0, 2.0, 4.0]) === nothing

    # Some expressions.
    @test getargs_tilde(:(x ~ Normal(μ, σ))) == (:x, :(Normal(μ, σ)))
    @test getargs_tilde(:((.~)(x, Normal(μ, σ)))) === nothing
    @test getargs_tilde(:((~).(x, Normal(μ, σ)))) === nothing
    @test getargs_tilde(:(@. x ~ Normal(μ, σ))) === nothing
    @test getargs_tilde(:(@. x ~ Normal(μ, $(Expr(:$, :(sqrt(v))))))) === nothing
    @test getargs_tilde(:(@~ Normal.(μ, σ))) === nothing
end

@testset "get_type" begin
    @test get_type(:(::Type{T})) == :T
    @test get_type(:(a::Type{A})) == :A
    @test get_type(:(::Type{T < Float64})) === nothing
end

@testset "get_symbol" begin
    @test get_symbol(:(x::Int)) == :x
    @test get_symbol(:(a::Type{A})) == :a
    @test get_symbol(:(::Type{A})) === nothing
    @test get_symbol(:(y::Vector{Int})) == :y
end


text:
module MyLinearAlgebra

export poly

using LinearAlgebra

function poly(x)

	m,n = size(x);
	if m == n
	   #Characteristic polynomial (square x)
	   e = eigen(x).values;
	elseif (m==1) || (n==1)
	   e = x;
	end

	#Strip out infinities
	e = e[isfinite.(e)];

	#Expand recursion formula
	n = length(e);
	c = [1 zeros(1,n)];
	for j=1:n
	    c[2:(j+1)] = c[2:(j+1)] - e[j].*c[1:j];
	end

	return c

end

end

max_stars_repo_name: maalaria/Newbe

max_issues_repo_name: maalaria/Newbe

max_forks_repo_name: maalaria/Newbe


text:
# Autogenerated wrapper script for Gnuplot_jll for x86_64-w64-mingw32-cxx11
export gnuplot

using LibGD_jll
using Qt5Svg_jll
using Cairo_jll
using Libiconv_jll
JLLWrappers.@generate_wrapper_header("Gnuplot")
JLLWrappers.@declare_executable_product(gnuplot)
function __init__()
    JLLWrappers.@generate_init_header(LibGD_jll, Qt5Svg_jll, Cairo_jll, Libiconv_jll)
    JLLWrappers.@init_executable_product(
        gnuplot,
        "bin\\gnuplot.exe",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using DiffEqBase

struct CTKAlg <: DiffEqBase.AbstractSteadyStateAlgorithm end

"""
    solve(prob::DiffEqBase.AbstractSteadyStateProblem,
          alg::CTKAlg;
          nrm=norm,
          τstop=1e12*365*24*60*60,
          preprint="",
          maxItNewton=50)

Solves `prob` using the modified C.T.Kelley Shamanskii algorithm.
"""
function DiffEqBase.solve(prob::DiffEqBase.AbstractSteadyStateProblem,
                          alg::CTKAlg;
                          nrm=norm,
                          τstop=ustrip(u"s", 1e6u"Myr"),
                          preprint="",
                          maxItNewton=50)
    # Define the functions according to DiffEqBase.SteadyStateProblem type
    p = prob.p
    t = 0
    x0 = copy(prob.u0)
    dx = copy(x0)
    F(x) = prob.f.f(dx, x, p, t)
    ∇ₓF(x) = prob.f.jac(x, p, t)
    # Compute `u_steady` and `resid` as per DiffEqBase using my algorithm
    x_steady = NewtonChordShamanskii(F, ∇ₓF, nrm, x0, τstop; preprint=preprint, maxItNewton=maxItNewton)
    resid = F(x_steady)
    # Return the common DiffEqBase solution type
    DiffEqBase.build_solution(prob, alg, x_steady, resid; retcode=:Success)
end

"""
    SteadyStateProblem(F, ∇ₓF, x, p)

Returns the `SteadyStateProblem` defined by `F(x,p)=0`.
"""
function DiffEqBase.SteadyStateProblem(F, ∇ₓF, x, p::AbstractParameters)
    return SteadyStateProblem(ODEFunction(F, ∇ₓF, x, p), x, p)
end

function DiffEqBase.ODEProblem(F, ∇ₓF, x, p::AbstractParameters; tspan=(0.0, ustrip(u"s", 1000.0u"yr")))
    return ODEProblem(ODEFunction(F, ∇ₓF, x, p), x, tspan, p)
end

function DiffEqBase.ODEFunction(F, ∇ₓF, x, p::AbstractParameters)
    f = (du,u,p,t) -> du .= F(u,p)
    jac = (u,p,t) -> ∇ₓF(u,p)
    jp = ∇ₓF(x,p)
    return ODEFunction(f; jac=jac, jac_prototype=jp)
end

export solve, SteadyStateProblem, CTKAlg


text:
mutable struct GridWorld
    n
    rewards
    probabilities
end

function GridWorld(n)
    p = Dict(:rewards => ones(n * n, 4, 1), :states => zeros(n * n, 4, n * n))
    r = [-1]
    for i in 1:n
        for j in 1:n
            s_index = cartesian_to_linear(i, j, n)
            if index_to_state(s_index, n) == "x"
                # pass
            elseif index_to_state(s_index, n) == "top-right"
                p[:states][s_index, 1, cartesian_to_linear(i, j - 1, n)] = 1.0
                p[:states][s_index, 2, cartesian_to_linear(i, j, n)] = 1.0  # fixed
                p[:states][s_index, 3, cartesian_to_linear(i, j, n)] = 1.0  # fixed
                p[:states][s_index, 4, cartesian_to_linear(i + 1, j, n)] = 1.0
            elseif index_to_state(s_index, n) == "bottom-left"
                p[:states][s_index, 1, cartesian_to_linear(i, j, n)] = 1.0  # fixed
                p[:states][s_index, 2, cartesian_to_linear(i, j + 1, n)] = 1.0
                p[:states][s_index, 3, cartesian_to_linear(i - 1, j, n)] = 1.0
                p[:states][s_index, 4, cartesian_to_linear(i, j, n)] = 1.0  # fixed
            elseif index_to_state(s_index, n) == "left"
                p[:states][s_index, 1, cartesian_to_linear(i, j, n)] = 1.0  # fixed
                p[:states][s_index, 2, cartesian_to_linear(i, j + 1, n)] = 1.0
                p[:states][s_index, 3, cartesian_to_linear(i - 1, j, n)] = 1.0
                p[:states][s_index, 4, cartesian_to_linear(i + 1, j, n)] = 1.0
            elseif index_to_state(s_index, n) == "right"
                p[:states][s_index, 1, cartesian_to_linear(i, j - 1, n)] = 1.0
                p[:states][s_index, 2, cartesian_to_linear(i, j, n)] = 1.0  # fixed
                p[:states][s_index, 3, cartesian_to_linear(i - 1, j, n)] = 1.0
                p[:states][s_index, 4, cartesian_to_linear(i + 1, j, n)] = 1.0
            elseif index_to_state(s_index, n) == "top"
                p[:states][s_index, 1, cartesian_to_linear(i, j - 1, n)] = 1.0
                p[:states][s_index, 2, cartesian_to_linear(i, j + 1, n)] = 1.0
                p[:states][s_index, 3, cartesian_to_linear(i, j, n)] = 1.0  # fixed
                p[:states][s_index, 4, cartesian_to_linear(i + 1, j, n)] = 1.0
            elseif index_to_state(s_index, n) == "bottom"
                p[:states][s_index, 1, cartesian_to_linear(i, j - 1, n)] = 1.0
                p[:states][s_index, 2, cartesian_to_linear(i, j + 1, n)] = 1.0
                p[:states][s_index, 3, cartesian_to_linear(i - 1, j, n)] = 1.0
                p[:states][s_index, 4, cartesian_to_linear(i, j, n)] = 1.0  # fixed
            else  index_to_state(s_index, n) == "o"
                p[:states][s_index, 1, cartesian_to_linear(i, j - 1, n)] = 1.00
                p[:states][s_index, 2, cartesian_to_linear(i, j + 1, n)] = 1.00
                p[:states][s_index, 3, cartesian_to_linear(i - 1, j, n)] = 1.00
                p[:states][s_index, 4, cartesian_to_linear(i + 1, j, n)] = 1.00
            end
        end
    end
    grid = Array{String}(undef, n, n)
    for i in eachindex(grid)
        grid[i] = index_to_state(i, n)
    end
    return GridWorld(n, r, p), grid
end

function index_to_state(index, n)
    """Return `state` in n x n GridWorld given index."""

    if (index == 1) || (index == n * n)
        return "x"
    elseif div(index - 1, n) == n - 1 && rem(index, n) == 1
        return "top-right"
    elseif div(index - 1, n) == 0 && rem(index, n) == 0
        return "bottom-left"
    elseif div(index - 1, n) == 0
        return "left"
    elseif div(index - 1, n) == n - 1
        return "right"
    elseif rem(index, n) == 1
        return "top"
    elseif rem(index, n) == 0
        return "bottom"
    else
        return "o"
    end
end

function cartesian_to_linear(i, j, n)
    return (j - 1) * n + i
end

function random_policy(n)
    """Create random policy matrix for GridWorld(n)"""
    return 0.25 * ones(n * n, 4)
end

function policy_evaluation(pA, pR, pS, R, gamma; tol::Float64=1e-6)
    """Evaluate policy via dynamic programming.

        - `pA`: array of action probabilities, p(a | s).
        - `pR`: array of reward probabilities, p(r | s, a).
        - `pS`: array of action probabilities, p(s' | s, a).
        - `R`: array of reward values *corresponding* to p(r | s, a).

        # Returns
        - `V`: vector of policy values, V(s).
    """

    V = zeros(size(pA)[1])
    steps = 0
    delta = 0.0
    while true
        delta = 0.0
        for (v_idx, v) in enumerate(V)
            v_update = 0.0
            for (a_idx, pa) in enumerate(pA[v_idx, :])
                for (r_idx, pr) in enumerate(pR[v_idx, a_idx, :])
                    r = R[r_idx]
                    for (s_idx, ps) in enumerate(pS[v_idx, a_idx, :])
                        v_update += pa * pr * ps * (r + gamma * V[s_idx])
                    end
                end
            end
            V[v_idx] = v_update
            # println("V[$v_idx] = $(V[v_idx])")
            delta = max(delta, abs(v - v_update))  # update maximum change
        end
        # println(V)
        steps += 1
        # println("delta=$delta")
        if delta < tol
            break
        end
    end
    println("converged in $steps steps (delta=$delta)")
    return V
end

function policy_evaluation(V0, pA, pR, pS, R, gamma; tol::Float64=1e-6)
    """Evaluate policy via dynamic programming.

        - `V0`: initial guess of value.
        - `pA`: array of action probabilities, p(a | s).
        - `pR`: array of reward probabilities, p(r | s, a).
        - `pS`: array of action probabilities, p(s' | s, a).
        - `R`: array of reward values *corresponding* to p(r | s, a).

        # Returns
        - `V`: vector of policy values, V(s).
    """

    V = copy(V0)
    steps = 0
    delta = 0.0
    while true
        delta = 0.0
        for (v_idx, v) in enumerate(V)
            v_update = 0.0
            for (a_idx, pa) in enumerate(pA[v_idx, :])
                for (r_idx, pr) in enumerate(pR[v_idx, a_idx, :])
                    r = R[r_idx]
                    for (s_idx, ps) in enumerate(pS[v_idx, a_idx, :])
                        v_update += pa * pr * ps * (r + gamma * V[s_idx])
                    end
                end
            end
            V[v_idx] = v_update
            # println("V[$v_idx] = $(V[v_idx])")
            delta = max(delta, abs(v - v_update))  # update maximum change
        end
        # println(V)
        steps += 1
        # println("delta=$delta")
        if delta < tol
            break
        end
    end
    # println("converged in $steps steps (delta=$delta)")
    return V
end

function policy_sweep(V0, pA, pR, pS, R, gamma)
    V = copy(V0)
    for (v_idx, v) in enumerate(V)
        v_update = 0.0
        for (a_idx, pa) in enumerate(pA[v_idx, :])
            for (r_idx, pr) in enumerate(pR[v_idx, a_idx, :])
                r = R[r_idx]
                for (s_idx, ps) in enumerate(pS[v_idx, a_idx, :])
                    v_update += pa * pr * ps * (r + gamma * V[s_idx])
                end
            end
        end
        V[v_idx] = v_update
    end
    return V
end

# TODO: random assignment on ties
function policy_improvement(V, nactions, pR, pS, R, gamma)
    policy = zeros(length(V), nactions)  # policy[s, a] = π(a | s)
    for (v_idx, v) in enumerate(V)
        a_max = 0
        Q_max = -Inf
        for a_idx in 1:nactions
            Q = 0.0
            for (s_idx, ps) in enumerate(pS[v_idx, a_idx, :])
                v_next = V[s_idx]
                for (r_idx, pr) in enumerate(pR[v_idx, a_idx, :])
                    r = R[r_idx]
                    Q += ps * pr * (r + v_next)
                end
            end
            if Q > Q_max
                a_max = a_idx
                Q_max = Q
            end
        end
        # println("state={$v_idx}, new action={$a_max}")
        policy[v_idx, a_max] = 1.0
    end
    return policy
end

function policy_iteration(policy, pR, pS, R, gamma, tol=1e-6)
    nstates, nactions = size(policy)
    value = zeros(nstates)
    i = 0
    while true
        value_new = policy_evaluation(value, policy, pR, pS, R, gamma, tol=tol)
        policy = policy_improvement(value_new, nactions, pR, pS, R, gamma)
        i += 1
        println("iter=$i, max_diff=$(maximum(abs.(value_new - value)))")
        if isapprox(value, value_new)
            break
        else
            value = value_new
        end
    end
    return policy, value
end

function value_iteration(pR, pS, R, gamma; tol=1e-6)
    nstates, nactions, _ = size(pS)
    V = zeros(nstates)
    policy = zeros(nstates, nactions)  # policy[s, a] = π(a | s)
    Δ = Inf
    i = 0
    while Δ > tol
        Δ = 0.0
        for (v_idx, v) in enumerate(V)
            a_max = 0
            Q_max = -Inf
            for a_idx in 1:nactions
                Q = 0.0
                for (s_idx, ps) in enumerate(pS[v_idx, a_idx, :])
                    v_next = V[s_idx]
                    for (r_idx, pr) in enumerate(pR[v_idx, a_idx, :])
                        r = R[r_idx]
                        Q += ps * pr * (r + v_next)
                    end
                end
                if Q > Q_max
                    a_max = a_idx
                    Q_max = Q
                end
            end
            V[v_idx] = Q_max
            policy[v_idx, a_max] = 1.0
            Δ = max(Δ, abs(V[v_idx] - v))
        end
        i += 1
    end
    println("Converged in $i steps (Δ = $Δ).")
    return policy, value
end

world, grid = GridWorld(4);
pA = random_policy(4);
pR = world.probabilities[:rewards];
pS = world.probabilities[:states];
R = world.rewards;
gamma = 1.0;  # no discounting
policy, value = policy_iteration(pA, pR, pS, R, gamma);
policy, value = value_iteration(pA)

function value_grid(value)
    n = Int(sqrt(size(value)[1]))
    return reshape(value, n, n)
end

function policy_grid(policy)
    idx_to_action = Dict(
        1 => "left",
        2 => "right",
        3 => "up",
        4 => "down"
    )
    n = Int(sqrt(size(policy)[1]))
    actions = map(idx -> idx_to_action[idx[2]], argmax(policy, dims=2))
    actions[1] = "x"
    actions[end] = "x"
    return reshape(actions, n, n)
end


text:
# @test_throw, with additional testing for the exception message
macro test_throws_message(f, typ, ex...)
    quote
        msg = ""
        @test_throws $(esc(typ)) try
            $(esc(ex...))
        catch err
            msg = sprint(showerror, err)
            rethrow()
        end

        if !$(esc(f))(msg)
            # @test should return its result, but doesn't
            @error "Failed to validate error message\n$msg"
        end
        @test $(esc(f))(msg)
    end
end

# NOTE: based on test/pkg.jl::capture_stdout, but doesn't discard exceptions
macro grab_output(ex)
    quote
        mktemp() do fname, fout
            ret = nothing
            open(fname, "w") do fout
                redirect_stdout(fout) do
                    ret = $(esc(ex))
                end
            end
            ret, read(fname, String)
        end
    end
end

# Run some code on-device, returning captured standard output
macro on_device(ex)
    quote
        let
            function kernel()
                $(esc(ex))
                return
            end

            @cuda kernel()
            synchronize()
        end
    end
end

# helper function for sinking a value to prevent the callee from getting optimized away
@inline sink(i::Int32) =
    Base.llvmcall("""%slot = alloca i32
                     store volatile i32 %0, i32* %slot
                     %value = load volatile i32, i32* %slot
                     ret i32 %value""", Int32, Tuple{Int32}, i)
@inline sink(i::Int64) =
    Base.llvmcall("""%slot = alloca i64
                     store volatile i64 %0, i64* %slot
                     %value = load volatile i64, i64* %slot
                     ret i64 %value""", Int64, Tuple{Int64}, i)

function julia_script(code, args=``)
    # FIXME: this doesn't work when the compute mode is set to exclusive
    script = "using CUDAnative; import CUDAdrv; $code"
    out = Pipe()
    err = Pipe()
    cmd = `$(Base.julia_cmd()) -e $script`
    if Base.JLOptions().project != C_NULL
        # --project isn't preserved by julia_cmd()
        cmd = `$cmd --project=$(unsafe_string(Base.JLOptions().project))`
    end
    cmd = `$cmd $args`
    proc = run(pipeline(cmd, stdout=out, stderr=err), wait=false)
    close(out.in)
    close(err.in)
    wait(proc)
    proc.exitcode, read(out, String), read(err, String)
end

# a lightweight CUDA array type for testing purposes
## ctor & finalizer
mutable struct CuTestArray{T,N}
    buf::Mem.Buffer
    shape::NTuple{N,Int}
    function CuTestArray{T,N}(shape::NTuple{N,Int}) where {T,N}
        len = prod(shape)
        buf = Mem.alloc(len*sizeof(T))

        obj = new{T,N}(buf, shape)
        finalizer(unsafe_free!, obj)
        return obj
    end
end
function unsafe_free!(a::CuTestArray)
    CUDAdrv.isvalid(a.buf.ctx) && Mem.free(a.buf)
end
## memory copy operations
function CuTestArray(src::Array{T,N}) where {T,N}
    dst = CuTestArray{T,N}(size(src))
    Mem.upload!(dst.buf, pointer(src), length(src) * sizeof(T))
    return dst
end
function Base.Array(src::CuTestArray{T,N}) where {T,N}
    dst = Array{T,N}(undef, src.shape)
    Mem.download!(pointer(dst), src.buf, prod(src.shape) * sizeof(T))
    return dst
end
## conversions
using Adapt
function Adapt.adapt_storage(::CUDAnative.Adaptor, a::CuTestArray{T,N}) where {T,N}
    ptr = Base.unsafe_convert(CuPtr{T}, a.buf)
    devptr = CUDAnative.DevicePtr{T,AS.Global}(ptr)
    CuDeviceArray{T,N,AS.Global}(a.shape, devptr)
end


text:
using MOOSE

include("NonlinearForce.jl")

# Create the Mesh
mesh = buildSquare(0, 1, 0, 1, 2, 2)

# Create the System to hold the equations
diffusion_system = System{Float64}(mesh)

# Add a variable to solve for
u = addVariable!(diffusion_system, "u")

# Apply the Laplacian operator to the variable
addKernel!(diffusion_system, Diffusion(u))

# Apply the NonlinearForce operator
addKernel!(diffusion_system, NonlinearForce(u))

# u = 0 on the Left
addBC!(diffusion_system, DirichletBC(u, [4], 0.0))

# u = 1 on the Right
addBC!(diffusion_system, DirichletBC(u, [2], 1.0))

# Initialize the system of equations
initialize!(diffusion_system)

# Create a solver and solve
solver = PetscNonlinearImplicitSolver(diffusion_system)
solve!(solver, nl_max_its=5)

# Output
out = VTKOutput()
output(out, solver, "nonlinear_force_out")


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Many auxiliary functions for sun synchronous orbit computations.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

export compute_ss_orbit_by_ang_vel
export compute_ss_orbit_by_inclination
export compute_ss_orbit_by_semi_major_axis

"""
    function compute_ss_orbit_by_ang_vel(n::Number, e::Number)

Compute the Sun-synchronous orbit given the angular velocity `n` [rad/s] and the
eccentricity `e`.

# Returns

* The semi-major axis [m].
* The inclination [rad].
* The residues of the two functions.
* A boolean variable that indicates if the numerical algorithm converged.

"""
function compute_ss_orbit_by_ang_vel(n::Number, e::Number)
    # Check if the arguments are valid.
    (n <= 0)          && throw(ArgumentError("The angular velocity must be greater than 0."))
    !( 0. <= e < 1. ) && throw(ArgumentError("The eccentricity must be within the interval 0 <= e < 1."))

    # Tolerance for the Newton-Raphson method.
    tol = 1e-18

    # Maximum number of iterations for the Newton-Raphson method.
    maxIt = 30

    # Auxiliary variables.
    sqrt_m0 = sqrt(m0)
    aux_e   = 1-e^2
    sqrt_e2 = sqrt(aux_e)

    # Auxiliary constant to compute the functions.
    K1 = 3R0^2*J2*sqrt_m0/(4aux_e^2)

    # Solve for zeros of f1 and f2 using Newton-Raphson method
    # ==========================================================================
    #
    # The function `f1` is the residue related to the time derivative of RAAN,
    # and the function `f2` is the residue related to the angular velocity.

    # Initial guess based on the unperturbed model.
    a_k = (m0/n^2)^(1/3)
    i_k = acos( -ne*a_k^(3.5)/(2K1) )

    # By setting the initial values of `f1` and `f2` to `10tol`, we assure that
    # the loop will be executed at least one time.
    f1 = 10tol
    f2 = 10tol

    # Loop
    it = 0
    converged = true

    while (abs(f1) > tol) || (abs(f2) > tol)
        a_k_1 = a_k
        i_k_1 = i_k

        # Auxiliary variables.
        sin_i, cos_i = sincos(i_k_1)
        sin_i²       = sin_i^2
        cos_i²       = cos_i^2
        a_pm1d5      = a_k_1^(-1.5)
        a_pm2d5      = a_pm1d5*a_k_1^(-1) # --> a_k_1^(-2.5)
        a_pm3d5      = a_pm2d5*a_k_1^(-1) # --> a_k_1^(-3.5)
        a_pm4d5      = a_pm3d5*a_k_1^(-1) # --> a_k_1^(-4.5)
        aux_1        = K1*sqrt_e2*(3cos_i²-1)
        aux_2        = K1*(5cos_i²-1)

        # Compute the Jacobian.
        df1da = -7K1*a_pm4d5*cos_i
        df1di = -2K1*a_pm3d5*sin_i
        df2da = -1.5*sqrt_m0*a_pm2d5 - 3.5aux_1*a_pm4d5 - 3.5aux_2*a_pm4d5
        df2di = -6K1*sqrt_e2*a_pm3d5*cos_i*sin_i - 10K1*a_pm4d5*cos_i*sin_i

        J_k_1 = SMatrix{2,2}(df1da, df1di,
                             df2da, df2di)'

        # Compute the functions.
        f1 = ne + 2K1*a_pm3d5*cos_i
        f2 = sqrt_m0*a_pm1d5 + aux_1*a_pm3d5 + aux_2*a_pm3d5 - n

        # Compute the new estimate.
        vf    = SVector{2}(f1,f2)
        v_k_1 = SVector{2}(a_k_1, i_k_1)
        v_k   = v_k_1 - inv(J_k_1)*vf
        a_k   = v_k[1]
        i_k   = v_k[2]

        # Check if the maximum number of iterations has been reached.
        it += 1

        # If the maximum number of iterations allowed has been reached, then
        # indicate that the solution did not converged and exit loop.
        if (it >= maxIt)
            converged = false
            break
        end
    end

    # Return.
    a_k, i_k, f1, f2, converged
end

"""
    function compute_ss_orbit_by_inclination(i::Number, e::Number)

Compute the Sun-synchronous orbit given the inclination `i` [rad] and the
eccentricity `e`.

# Returns

The semi-major axis of the Sun-synchronous orbit [m].

"""
function compute_ss_orbit_by_inclination(i::Number, e::Number)
    !( 0. <= e < 1. ) && throw(ArgumentError("The eccentricity must be within the interval 0 <= e < 1."))

    # Auxiliary variables.
    sqrt_m0 = sqrt(m0)

    # Auxiliary constant.
    K1 = 3.0*R0^2*J2*sqrt_m0/(4.0*(1-e^2)^2)

    # Compute the semi-major axis.
    a = (-2*cos(i)*K1/ne)^(2.0/7.0)

    # Check if the orbit is valid.
    ( a*(1-e) <  R0 ) && throw(ErrorException("It was not possible to find a valid Sun-synchronous orbit with the inclination given."))

    # Return.
    a
end

"""
    function compute_ss_orbit_by_semi_major_axis(a::Number, e::Number)

Compute the Sun-synchronous orbit given the semi-major axis `a` [m] and the
eccentricity `e`.

# Returns

The inclination of the Sun-synchronous orbit [rad].

"""
function compute_ss_orbit_by_semi_major_axis(a::Number, e::Number)
    # Check if the arguments are valid.
    (a*(1-e) <= R0)   && throw(ArgumentError("The perigee must be larger than the Earth radius."))
    !( 0. <= e < 1. ) && throw(ArgumentError("The eccentricity must be within the interval 0 <= e < 1."))

    # Auxiliary variables.
    sqrt_m0 = sqrt(m0)

    # Auxiliary constant.
    K1 = 3R0^2*J2*sqrt_m0/(4(1-e^2)^2)

    # Compute the inclination.
    cos_i = -ne*a^3*sqrt(a)/(2K1)

    # Check if -1 <= cos_i <= 1.
    ( (cos_i < -1) || (cos_i > 1) ) && throw(ErrorException("It was not possible to find a Sun-synchronous orbit with the semi-major axis given."))

    # Return.
    acos(cos_i)
end


text:
@testset "simulate.jl" begin
    println("             > simulate.jl")
    N = 10
    p = 0.85
    q = 1.0
    σ = 0.2
    τ = 0.2
    sim = NestedFilterSimulation(
            N, p, q, σ, τ,
            1:20,
            LinRange(0.00, 1.00, 100),
            LinRange(0.00, 2.00, 100),
            LinRange(0.05, 2.00, 100),
            LinRange(0.05, 2.00, 100),
            2048, 512, 12,
            timestep = FixedTimestep(0.1)
          )

    for i in 1:10
        obs = propagate!(sim)
        @test obs.dt ≈ 0.1
    end

    times, epsps = run!(sim, T = 10)
    for i in eachindex(times)
        @test times[i] ≈ 0.1*i
    end
end


text:
@testset "VerifyOrder" begin

    using DataFrames
    using Tables

    @test Tables.rows(DataFrame(Chrom = ["chr1", "chr1", "chr1"],
                                Pos = [1,2,3])) |> verifyorder |> collect |> length == 3

    @test_throws ErrorException  Tables.rows(DataFrame(Chrom = ["chr1", "chr1", "chr1"],
                                                       Pos = [1,3,2])) |> verifyorder |> collect

    @test_throws ErrorException  Tables.rows(DataFrame(Chrom = ["chr1", "chr2", "chr1"],
                                                       Pos = [1,1,1])) |> verifyorder |> collect

end


text:
using BinaryBuilder, Printf

name = "PATHlib"
version = v"4.7.3"

# I guess they're hoping for double-digit numbers of patch releases.  /shrug
version_str = @sprintf("%d.%d.%02d", version.major, version.minor, version.patch)

# Collection of sources required to build
sources = [
    "https://github.com/ampl/pathlib/archive/$(version_str).tar.gz" =>
    "93244121cb03d1c726fcb4e33aa86e8cd59864c873eb03733b07aaef7f448ed8",
    "./bundled",
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/

mkdir -p ${prefix}/include ${prefix}/lib ${libdir}

function target_libdir()
{
    if [[ ${target} == x86_64-linux-gnu ]]; then
        echo -n linux64
    elif [[ ${target} == i686-linux-gnu ]]; then
        echo -n linux32
    elif [[ ${target} == x86_64-apple-* ]]; then
        echo -n osx
    elif [[ ${target} == x86_64-w64-* ]]; then
        echo -n win64
    elif [[ ${target} == i686-w64-* ]]; then
        echo -n win32
    else
        echo "ERROR: We don't have an upstream PATHlib release for ${target}" >&2
        exit 1
    fi
}

# pathlib comes as precompiled binaries, so we just copy the binaries and the include directory over
# We move the `dlext` stuff first as that sometimes goes into `bin`, and everything else into `lib`
mv -v pathlib-*/lib/$(target_libdir)/*.${dlext} ${libdir}
mv -v pathlib-*/lib/$(target_libdir)/* ${prefix}/lib
mv -v pathlib-*/include/*.h ${prefix}/include

# Add linkage args
LDFLAGS="${LDFLAGS} -lm -lgfortran"
LIBNAME="libpath47julia.${dlext}"

# By linking against the static libpath library, we are able to
LIBPATH_A="libpath47.a"
if [[ ${target} == *mingw* ]]; then
    LIBPATH_A="path47.lib"
fi

# Add compiler flags, some platform-dependent
CFLAGS="${CFLAGS} -I${prefix}/include -fPIC -shared"
if [[ ${target} == *linux* ]]; then
    CFLAGS="${CFLAGS} -Wl,-soname,${LIBNAME}"
elif [[ ${target} == *mingw* ]]; then
    CFLAGS="${CFLAGS} -DFNAME_LCASE_NODECOR -DUSE_OUTPUT_INTERFACE"
fi

cc -o ${libdir}/${LIBNAME} pathjulia.c ${prefix}/lib/${LIBPATH_A} ${CFLAGS} ${LDFLAGS}

# Install license file
install_license pathlib-*/LICENSE
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Platform("x86_64", "linux"; libc="glibc"),
    Platform("i686", "linux"; libc="glibc"),
    Platform("x86_64", "macos"),
    Platform("x86_64", "windows"),
    Platform("i686", "windows"),
]
platforms = expand_gfortran_versions(platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct("libpath47julia", :libpath47julia),
]

# Dependencies that must be installed before this package can be built
dependencies = [
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil