text:
export race

import Base: show

"""
    race(sources...)
    race(sources::S) where { S <: Tuple }

Combines multiple Observables to create an Observable that mirrors the output of the first Observable to emit an item.
Essentially it subscribes to the observable that was the first to start emitting.

# Arguments
- `sources`: input sources

# Examples
```jldoctest
using Rocket

source1 = of(1)
source2 = of(2)

subscribe!(race(source1, source2), logger())
;

# output
[LogActor] Data: 1
[LogActor] Completed
```

```
using Rocket

source1 = of(1) |> async()
source2 = of(2)

subscribe!(race(source1, source2), logger())
;

# output
[LogActor] Data: 2
[LogActor] Completed
```
See also: [`Subscribable`](@ref), [`subscribe!`](@ref)
"""
race()                                 = error("race operator expects at least one inner observable on input")
race(args...)                          = race(args)
race(sources::S) where { S <: Tuple }  = RaceObservable{union_type(sources), S}(sources)
race(sources::V) where { V <: Vector } = RaceObservable{union_type(sources), V}(sources)

##

struct RaceInnerActor{L, W, I} <: Actor{L}
    wrapper :: W
end

Base.show(io::IO, inner::RaceInnerActor{L, W, I}) where { L, W, I } = print(io, "RaceInnerActor($L, $I)")

on_next!(actor::RaceInnerActor{L, W, I}, data::L) where { L, W, I } = next_received!(actor.wrapper, data, Val{I}())
on_error!(actor::RaceInnerActor{L, W, I}, err)    where { L, W, I } = error_received!(actor.wrapper, err, Val{I}())
on_complete!(actor::RaceInnerActor{L, W, I})      where { L, W, I } = complete_received!(actor.wrapper, Val{I}())

##

mutable struct RaceActorWrapper{A}
    actor               :: A
    subscriptions       :: Vector{Teardown}
    first_emmited_index :: Union{Nothing, Int}

    RaceActorWrapper{A}(actor::A) where A = new(actor, Vector{Teardown}(), nothing)
end

has_emmited(wrapper::RaceActorWrapper)                     = get_first_emmited_index(wrapper) !== nothing
get_first_emmited_index(wrapper::RaceActorWrapper)         = wrapper.first_emmited_index
set_first_emmited_index!(wrapper::RaceActorWrapper, index) = wrapper.first_emmited_index = index

function next_received!(wrapper::RaceActorWrapper, data, index::Val{I}) where I
    first = get_first_emmited_index(wrapper)
    if first === nothing
        set_first_emmited_index!(wrapper, I)
        dispose_except(wrapper, index)
        next_received!(wrapper, data, index)
    elseif first === I
        next!(wrapper.actor, data)
    end
end

function error_received!(wrapper::RaceActorWrapper, err, index::Val{I}) where I
    first = get_first_emmited_index(wrapper)
    if first === nothing
        set_first_emmited_index!(wrapper, I)
        error_received!(wrapper, err, index)
    elseif first === I
        dispose(wrapper)
        error!(wrapper.actor, err)
    end
end

function complete_received!(wrapper::RaceActorWrapper, index::Val{I}) where I
    first = get_first_emmited_index(wrapper)
    if first === nothing
        set_first_emmited_index!(wrapper, I)
        complete_received!(wrapper, index)
    elseif first === I
        dispose(wrapper)
        complete!(wrapper.actor)
    end
end

dispose(wrapper::RaceActorWrapper) = foreach(s -> unsubscribe!(s), wrapper.subscriptions)

function dispose_except(wrapper::RaceActorWrapper, index::Val{I}) where I
    for (i, subscription) in enumerate(wrapper.subscriptions)
        if i !== I
            unsubscribe!(subscription)
        end
    end
end

##

@subscribable struct RaceObservable{T, S} <: Subscribable{T}
    sources :: S
end

function on_subscribe!(observable::RaceObservable, actor::A) where { A }
    wrapper = RaceActorWrapper{A}(actor)
    for (index, source) in enumerate(observable.sources)
        push!(wrapper.subscriptions, subscribe!(source, RaceInnerActor{eltype(source), typeof(wrapper), index}(wrapper)))
        if has_emmited(wrapper)
            break
        end
    end
    return RaceSubscription(wrapper)
end

##

struct RaceSubscription{W} <: Teardown
    wrapper :: W
end

as_teardown(::Type{ <: RaceSubscription }) = UnsubscribableTeardownLogic()

function on_unsubscribe!(subscription::RaceSubscription)
    dispose(subscription.wrapper)
    return nothing
end

Base.show(io::IO, ::RaceObservable{D}) where D  = print(io, "RaceObservable($D)")
Base.show(io::IO, ::RaceSubscription)           = print(io, "RaceSubscription()")


text:
using HDF5,JLD,KUparser,KUnet,Base.Test
include("pxyequal.jl")
infoln(x)=(println();info(x))

infoln("Loading data")
@date @load "conll07.tst.jld4"
@show nbatch = 10
@show ndeps = length(deprel)
@show ft = Flist.acl11eager
@show pt = ArcEager13
p0 = pt(1,ndeps)
s1 = corpus[1]

infoln("Let's train a reasonable model")
h0 = 20000
net = [Mmul(h0), Bias(), Relu(),
       Mmul(p0.nmove), Bias(), Logp(),
       LogpLoss()]
@date (p,x,y)=oparse(pt, corpus, ndeps, ft)
@date train(net, hcat(x...), hcat(y...))

infoln("0. Use gparser as reference for beam=1")
@date pxy0=(p0,x0,y0)=gparse(pt, corpus, ndeps, ft, net; returnxy=true)
@show map(size, pxy0)
@show evalparse(p0, corpus)

println()
@show nbeam = 1

infoln("1. nbeam=1 Single sentence")
@date p1=bparse(pt, [s1], ndeps, ft, net, nbeam)
@test @show isequal(p1[1],p0[1])

infoln("2. nbeam=1 Single sentence, xy")
@date pxy2=(p2,x2,y2)=bparse(pt, [s1], ndeps, ft, net, nbeam; returnxy=true)
@show map(size, pxy2)
@show I2 = 1:size(x2[1],2)
infoln("these may not be equal because of early stop")
@show isequal(p2[1],p0[1])
@show isequal(x2[1],x0[1][:,I2])
@show isequal(y2[1],y0[1][:,I2])

infoln("3. nbeam=1 Multiple sentences")
@date p3=bparse(pt, corpus, ndeps, ft, net, nbeam)
@test @show isequal(p3,p0)

infoln("4. nbeam=1 Multiple sentences, xy")
@date pxy4=(p4,x4,y4)=bparse(pt, corpus, ndeps, ft, net, nbeam; returnxy=true)
@show map(size, pxy4)
@show isequal(p4,p0)
@test @show isequal(p4[1],p2[1])
@test @show isequal(x4[1][:,I2],x2[1])
@test @show isequal(y4[1][:,I2],y2[1])

infoln("5. nbeam=1 Batch processing")
@date p5=bparse(pt, corpus, ndeps, ft, net, nbeam; nbatch=nbatch)
@test @show isequal(p5,p0)

infoln("6. nbeam=1 Batch processing, xy")
@date pxy6=(p6,x6,y6)=bparse(pt, corpus, ndeps, ft, net, nbeam; nbatch=nbatch, returnxy=true)
@show map(size,pxy6)
@test @show pxyequal(pxy6, pxy4)

if nworkers() > 1
    infoln("7. nbeam=1 Multi-cpu processing")
    @date p7=bparse(pt, corpus, ndeps, ft, net, nbeam; nbatch=nbatch, usepmap=true)
    @test @show isequal(p7,p0)
    infoln("8. nbeam=1 Multi-cpu processing, xy")
    @date pxy8=bparse(pt, corpus, ndeps, ft, net, nbeam; nbatch=nbatch, usepmap=true, returnxy=true)
    @test @show pxyequal(pxy8, pxy4)
else
    infoln("Use 'julia -p n' to test multi-cpu.")
end

println()
@show nbeam = 10

infoln("1. nbeam=10 Single sentence")
@date q1=bparse(pt, [s1], ndeps, ft, net, nbeam)
@show evalparse(q1, [s1])

infoln("2. nbeam=10 Single sentence, xy")
@date qxy2=(q2,x2,y2)=bparse(pt, [s1], ndeps, ft, net, nbeam; returnxy=true)
@show map(size, (q2,x2,y2))
@show I2 = 1:size(x2[1],2)
@show isequal(q2,q1) # may not be equal due to early stop

infoln("3. nbeam=10 Multiple sentences")
@date q3=bparse(pt, corpus, ndeps, ft, net, nbeam)
@show evalparse(q3, corpus)
@test @show isequal(q3[1],q1[1])

infoln("4. nbeam=10 Multiple sentences, xy")
@date qxy4=(q4,x4,y4)=bparse(pt, corpus, ndeps, ft, net, nbeam; returnxy=true)
@show map(size, qxy4)
@show isequal(q4, q3)           # may not equal, earlystop
@test @show isequal(q4[1], q2[1])
@test @show isequal(x4[1][:,I2], x2[1])
@test @show isequal(y4[1][:,I2], y2[1])

infoln("5. nbeam=10 Batch processing")
@date q5=bparse(pt, corpus, ndeps, ft, net, nbeam; nbatch=nbatch)
@test @show isequal(q5,q3)

infoln("6. nbeam=10 Batch processing, xy")
@date qxy6=(q6,x6,y6)=bparse(pt, corpus, ndeps, ft, net, nbeam; nbatch=nbatch, returnxy=true)
@show map(size, qxy6)
@test @show pxyequal(qxy6, qxy4)

if nworkers() > 1
    infoln("7. nbeam=10 Multi-cpu processing")
    @date q7=bparse(pt, corpus, ndeps, ft, net, nbeam; nbatch=nbatch, usepmap=true)
    @test @show isequal(q7, q3)
    infoln("8. nbeam=10 Multi-cpu processing, xy")
    @date qxy8=bparse(pt, corpus, ndeps, ft, net, nbeam; nbatch=nbatch, usepmap=true, returnxy=true)
    @test @show pxyequal(qxy8, qxy6)
else
    infoln("Use 'julia -p n' to test multi-cpu.")
end

ext: jl


text:
struct HospitalInfectionParams
    ishealthcare::BitVector
    hospital_staff_ids::Vector{PersonIdx}
    kernel_constant::Float64

    healthcare_detection_prob::Float64
    healthcare_detection_delay::Float64
  end

ishealthcare(params::HospitalInfectionParams, person_id::Integer) = params.ishealthcare[person_id]

function saveparams(dict, p::HospitalInfectionParams, prefix::AbstractString="")
  dict[prefix*"kernel_constant"] = p.kernel_constant
  dict[prefix*"detection_prob"] = p.healthcare_detection_prob
  dict[prefix*"healthcare_detection_delay"] = p.healthcare_detection_delay
end


text:
module testInvestments

using Base.Test
using DataFrames
using Datetime
using TimeData
include(string(Pkg.dir("AssetMgmt"), "/src/AssetMgmt.jl"))

println("\n Running investments tests\n")

## using Markowitz

vals = rand(8, 4)
vals = AssetMgmt.makeWeights(vals)
valsDf = DataFrame(vals)

invs = AssetMgmt.Investments(valsDf, [1:8])

######################
## evolving weights ##
######################

wgts = [0.4 0.3 0.3; 0.3 0.4 0.3]
rets = [0.08 0.02 0.06; 0.08 0.02 0.06]
pRets = AssetMgmt.invRetCore(wgts, rets)
expWgts = [0.4090909090909091 0.28977272727272724 0.30113636363636365;
           0.30857142857142855 0.38857142857142857 0.3028571428571428]
@test_approx_eq(expWgts, AssetMgmt.evolWgtsCore(wgts, rets))

end


text:
@testset "Utilities" begin

    @testset "Group actions" begin
        f1 = s -> (s * s,);
        f2 = s-> (2s, -s, 5s);
        f3 = s -> (s + 1,);
        action1 = GroupActions(f1)
        action2 = GroupActions(f1, f2)
        action3 = GroupActions(f1, f2, f3)
        r1 = Int[]
        HC.apply_actions(action1, 3) do s
            push!(r1, s)
            false
        end
        @test r1 == [9]

        r2 = Int[]
        HC.apply_actions(action2, 3) do s
            push!(r2, s)
            false
        end
        @test r2 == [9, 18, -9, 45, 6, -3, 15]

        r3 = Int[]
        HC.apply_actions(action3, 3) do s
            push!(r3, s)
            false
        end
        @test r3 == [9, 18, 19, -9, -8, 45, 46, 10, 6, 7, -3, -2, 15, 16, 4]
        @test action3(3) == [3, 9, 18, 19, -9, -8, 45, 46, 10, 6, 7, -3, -2, 15, 16, 4]

        # also test with arrays
        g1 = s -> [s * s];
        g2 = s-> [2s, -s, 5s];
        g3 = s -> [s + 1];
        action1 = GroupActions(g1)
        action2 = GroupActions(g1, g2)
        action3 = GroupActions(g1, g2, g3)
        r1 = Int[]
        HC.apply_actions(action1, 3) do s
            push!(r1, s)
            false
        end
        @test r1 == [9]

        r2 = Int[]
        HC.apply_actions(action2, 3) do s
            push!(r2, s)
            false
        end
        @test r2 == [9, 18, -9, 45, 6, -3, 15]

        r3 = Int[]
        HC.apply_actions(action3, 3) do s
            push!(r3, s)
            false
        end
        r3
        @test r3 == [9, 18, 19, -9, -8, 45, 46, 10, 6, 7, -3, -2, 15, 16, 4]
    end

    @testset "UniquePoints" begin
        Random.seed!(1234)
        X = [randn(ComplexF64, 10) for _ = 1:2_000]
        indices = HC.unique!(rand(1:2000, 20))
        data = HC.UniquePoints(X)

        test_show_juno(data)
        @test string(data) == "UniquePoints{Array{Complex{Float64},1},Float64,typeof(euclidean_distance),Nothing} with 2000 points"

        @test length(data) == 2_000

        for i ∈ indices
            @test HC.is_contained(data, X[i])
            @test HC.is_contained(data, X[i], Val(true)) == i
            @test HC.is_contained(data, X[i] .+ 1e-4) == false
            @test HC.is_contained(data, X[i] .+ 1e-9, Val(true)) == i
            @test HC.is_contained(data, X[i] .+ 1e-9) == true
            @test HC.add!(data, X[i]) == false
            @test HC.add!(data, X[i], Val(true)) == i
            @test data[i] == X[i]
            @test HC.add!(data, X[i] .+ 1e-4) == true
            @test HC.add!(data, X[i] .- 1e-4, Val(true)) == -1
        end

        empty!(data)
        @test length(points(data)) == 0

        # Test many points with nearly indentical distance to the inserted point
        p = shuffle!([[cis(k/100*2π)] for k=0:99])
        data = HC.UniquePoints(p)
        @test HC.is_contained(data, [0.0im]) == false

        # Test with group action
        x = randn(ComplexF64, 4)
        permutation1(x) = ([x[2]; x[1]; x[3]; x[4]],)
        permutation2(x) = ([x[1]; x[2]; x[4]; x[3]],)
        X = [v for v in GroupActions(permutation1, permutation2)(x)]

        # One group action
        data = UniquePoints(X, group_action = permutation1)
        @test length(data) == 2

        # Two group actions
        data = HC.UniquePoints(X, group_actions = GroupActions(permutation1, permutation2))
        @test length(data) == 1

        # Group action and reality check
        x = randn(4)
        X = [im.*x, x, randn(ComplexF64, 4)]

        data = HC.UniquePoints(X, group_action = x -> (im.*x, (-1).*x, (-im).*x), check_real=false)
        @test HC.is_real_vector(points(data)[1]) == false
        @test length(points(data)) == 2

        data = HC.UniquePoints(X, group_action = x -> (im.*x, (-1).*x, (-im).*x), check_real = true)
        @test HC.is_real_vector(points(data)[1]) == true
        @test length(points(data)) == 2

        data = HC.UniquePoints(randn(ComplexF64, 4), check_real = true)
        k = add!(data, randn(4), Val{true}())
        @test k == -2
        k = add!(data, randn(ComplexF64, 4), Val{true}())
        @test k == -1

        pts = [[1.0,0.5], [1.0,0.5], [0.5,1.0]]
        @test points(UniquePoints(pts)) == [[1.0, 0.5], [0.5, 1.0]]
        @test unique_points(pts) == [[1.0, 0.5], [0.5, 1.0]]
        @test unique_points(pts; group_action = x -> ([x[2],x[1]],)) == [[1.0, 0.5]]
        @test unique_points(pts; group_action = x -> [x[2],x[1]]) == [[1.0, 0.5]]
    end

    @testset "Multiplicities" begin
        # test empty input
        @test isempty(multiplicities(Vector{Vector{ComplexF64}}()))

        V = [randn(3) + im.* randn(3) for i in 1:10]
        W = [map(v -> v + [1e-7 * v[1] ;1.0;1.0], V); V; map(v -> v + [1e-6 * v[1]; 0; 0], V)]
        U = [[cis(rand()) for _=1:3] for _=1:20]
        X = [[0.499064-0.707804im, -0.164223+0.0933025im, -0.00695332-0.139516im, 0.000716617-0.135055im, 0.009846-0.12984im, 0.020204-0.123609im, 0.34326+0.161402im], [-0.795428-0.0975931im, -0.20942-0.382162im, -0.119855-0.0414809im, -0.123105-0.0293171im, -0.128161-0.0148951im, -0.134641+0.00272884im, 0.258239+0.185188im], [0.605937-0.289401im, 0.42308-0.188036im, -0.057496+0.200812im,0.284444+0.0477843im, 0.191189-0.225824im, -0.233732-0.181827im, -0.0939606+0.154403im], [-0.493372+0.512561im, -0.385166+0.204525im, -0.120853-0.142237im, 0.0742939+0.163197im, -0.298456+0.202635im, 0.100108-0.219369im, 0.193075+0.105514im], [0.483174-0.581023im, 0.23432-0.0528133im, -0.111264-0.067633im, -0.292867-0.014453im, 0.24805-0.267034im, 0.141232-0.13631im, 0.168603+0.261794im], [-0.455116+0.4482im, -0.366004+0.264082im, -0.0578486+0.110708im, -0.296085-0.0619238im, 0.216907-0.327015im, 0.207934+0.242906im, -0.121251+0.101268im], [0.526416-0.465825im, 0.39214-0.0908531im, -0.260222+0.00391701im, -0.110348-0.1151im, 0.101057-0.190265im, 0.308379-0.137148im, 0.0488701+0.29665im], [0.147886-0.750864im, -0.235168-0.257487im, 0.0394847+0.194512im, -0.0325228-0.227443im, -0.111526-0.214717im, 0.0245065+0.162303im, 0.304189-0.150105im], [-0.372175-0.538886im, -0.117563-0.341994im, 0.0752294-0.0141672im, 0.143671-0.249302im, -0.374165-0.228281im, -0.0971581+0.346633im, 0.168499+0.0394346im], [-0.503759+0.48im, -0.447123+0.262057im, -0.214529+0.109917im, 0.148493-0.0330798im, -0.0510215-0.262778im, -0.110165+0.198317im, 0.200704+0.0508037im], [0.623685-0.513102im, 0.0522679+0.0223332im, -0.188478-0.10101im, 0.229926-0.157889im, -0.146424-0.0562049im, 0.209112-0.185147im, 0.205229+0.276164im], [0.127017-0.655061im, -0.147053-0.569588im, -0.134863+0.0901295im, 0.152888+0.049255im, -0.124224+0.0624774im, 0.152924+0.0778815im, 0.105575-0.310906im], [0.575711+0.384346im, 0.336129+0.0531313im,0.155698-0.152359im, -0.104555+0.197167im, -0.254208-0.215823im, 0.327345+0.0728999im, -0.0274605+0.288691im], [0.479139-0.593553im, -0.0960412-0.315975im, 0.0438803-0.245038im, 0.000516985+0.193173im, -0.0828034+0.172297im, -0.00479366-0.21881im, 0.353615-0.0168729im], [-0.37103+0.609109im, -0.201984+0.259807im, -0.158043+0.252162im, -0.211405-0.106454im, 0.0274292-0.262351im, 0.333102+0.0455649im, -0.0847849+0.220536im], [-0.200697+0.682536im, -0.247073+0.359358im, -0.172902-0.0702972im, 0.141297+0.110399im, -0.171049+0.317616im, -0.0125869-0.240803im, 0.219982+0.00435721im], [-0.0896611-0.805671im, -0.0393637-0.183565im, 0.224508-0.167715im, 0.0918438-0.180778im, -0.0879493-0.157071im, -0.250607-0.0700599im, 0.0521971+0.291749im], [-0.226054+0.605127im, -0.210133+0.55141im,0.218162-0.0576634im, 0.08772-0.0936im, -0.0820928-0.101495im, -0.22419-0.0507194im, -0.0146449+0.311527im], [0.498516+0.360158im, 0.219394+0.0657585im, -0.0836376+0.0106629im, 0.00702255-0.322935im, 0.481634+0.13367im, -0.133322+0.397338im, -0.171875-0.0532251im], [0.440673-0.62568im, -0.110019-0.330907im, -0.043096+0.193744im, 0.0628576-0.220991im, -0.0144861-0.24152im, -0.0436694+0.158228im, 0.338948-0.0133119im], [0.622749+0.219955im, 0.628838+0.216896im, 0.0331851+0.145177im, -0.123675-0.0636893im, 0.0348129+0.117408im, -0.133898-0.0420959im, 0.155277-0.144857im], [0.652186-0.0569317im, 0.552174-0.169926im, -0.0977967-0.0521929im, -0.161552-0.169677im, 0.00437006+0.237406im, -0.0164646+0.116656im, 0.310812-0.0599935im], [0.448293+0.574542im, 0.247779+0.106134im, -0.0955677-0.215494im, 0.0935008+0.252581im, 0.336093-0.0655203im, -0.242972-0.0256941im, -0.0177608+0.301502im], [0.445422+0.671948im, -0.0726023+0.313409im, -0.044775+0.21726im,0.154346+0.0261763im, -0.0305679+0.188435im, 0.20006+0.0391233im, -0.0810809-0.297159im], [-0.00979382-0.81029im, 0.106391-0.195197im, 0.108622-0.139201im, 0.234362-0.1806im, -0.20136-0.131455im, -0.11823-0.118552im, -0.0308846+0.297292im], [0.318611+0.584183im, 0.105547+0.221636im, 0.248174+0.0583734im, -0.00256072-0.290164im, -0.302039-0.0824955im, -0.114841+0.409322im, 0.25213+0.0736154im], [0.696901+0.0684572im, 0.505426+0.266373im, -0.152493+0.0711982im, 0.197532-0.0412531im, 0.214441+0.0165324im, -0.096866+0.0847148im, -0.103987-0.201388im], [0.115188-0.657244im, -0.157297-0.566844im, -0.133216+0.092546im, 0.153751+0.0464909im, -0.123077+0.0647066im, 0.154303+0.0751121im, 0.0999536-0.312759im], [0.653199-0.277909im, 0.556407-0.00960311im, 0.204776-0.112518im, -0.112199+0.0771487im, -0.0695128+0.163961im, 0.167741-0.0723464im, -0.16148-0.148729im], [0.411538+0.550657im, 0.413916+0.420472im, -0.0555052+0.089203im, -0.0609812+0.227821im, -0.0913546-0.187807im, -0.0607857-0.0922333im, 0.237526+0.0198403im], [0.602049+0.184366im, 0.695878+0.270759im, -0.0420156+0.0359236im, -0.0431941+0.0331572im, -0.0446403+0.0302374im, -0.046326+0.027106im, 0.0979775-0.156785im], [-0.270799+0.567571im, -0.113714+0.377218im, 0.0592799-0.219174im, 0.226587+0.091218im, 0.0269123+0.326853im, -0.402766+0.100383im, -0.018186-0.240642im]]

        M = multiplicities(V, distance=infinity_distance)
        @test length(M) == 0

        N = multiplicities(W, distance=infinity_distance)
        sort!(N, by=first)
        @test length(N) == 10
        @test unique([length(m) for m in N]) == [2]
        @test N == [[i, i+10] for i=11:20]

        O = multiplicities([U;U], distance=infinity_distance)
        @test length(O) == 20

        P = multiplicities(X, distance=(x,y) -> 1-abs(LinearAlgebra.dot(x,y)))
        @test length(P) == 3

        # Test with group action
        x = randn(ComplexF64, 4)
        permutation1(x) = ([x[2]; x[1]; x[3]; x[4]],)
        permutation2(x) = ([x[1]; x[2]; x[4]; x[3]],)
        X = GroupActions(permutation1, permutation2)(x)

        # One group action
        m = multiplicities(X, group_action = permutation1) |> sort
        @test m == [[1, 2], [3, 4]]

        # Two group actions
        m = multiplicities(X, group_actions = GroupActions(permutation1, permutation2))
        @test m == [[1, 2, 3, 4]]

        # docs
        X = [[1, 2, 3, 4], [2,1,3,4], [1,2,4,3], [2,1,4,3]]
        permutation(x) = [x[2], x[1], x[3], x[4]]
        m = sort(multiplicities(X, group_action = permutation))
        @test m == [[1,2], [3,4]]
    end

    @testset "Polynomials" begin
        @polyvar p q a b c x y z
        e = [p + 1]
        f = [a * b * c ]
        g = [x+y, y + z, x + z]
        @test expand(f ∘ g) == [x^2*y + x^2*z + x*y^2 + 2*x*y*z + x*z^2 + y^2*z + y*z^2]
        @test expand(e ∘ f ∘ g) == [x^2*y + x^2*z + x*y^2 + 2*x*y*z + x*z^2 + y^2*z + y*z^2 + 1]
        @test expand(e ∘ (f ∘ g)) == [x^2*y + x^2*z + x*y^2 + 2*x*y*z + x*z^2 + y^2*z + y*z^2 + 1]
        @test expand((e ∘ f) ∘ g) == [x^2*y + x^2*z + x*y^2 + 2*x*y*z + x*z^2 + y^2*z + y*z^2 + 1]
        @test e ∘ f == e ∘ f

        @test validate(f ∘ g) == true
        @test validate(f ∘ g, parameters=[z]) == true
        @test validate(f ∘ g, parameters=[c]) == false
        @test validate(g ∘ f) == false

        @test HC.nvariables(f ∘ g) == 3
        @test HC.nvariables(f ∘ g, parameters=[z]) == 2

        @polyvar x y z
        @test is_homogeneous([x^2+y^2+x*y, x^5])
        @test is_homogeneous([x^2+y^2+x*y, x^4+1]) == false

        #test weighted degree
        @test is_homogeneous(x^3+x*y, [(x, 2), (y, 4)])
        @test homogenize(x+x*y, [(x, 2), (y, 4)], z) == x*z^4+x*y

        @test is_homogeneous(homogenize([x^2+y^2+x*y, x^4+1]))
        @test is_homogeneous([x^2+z^2 + y, x^4+z^4], [x,z]) == false
        @test is_homogeneous(homogenize([x^2+z^2*y, x^4+z^4*y], [x,z]), [x,z]) == true

        @test is_homogeneous(f ∘ g) == true
        h = [a * b * c  + p * c^3] ∘ [x+y, y + z, x + z]
        @test is_homogeneous(h, parameters=[p]) == true
        h2 = [a * b * c  + p] ∘ [x+y, y + z, x + z]
        @test is_homogeneous(h2, parameters=[p]) == false

        #homogenize
        @test homogenize(x^2+y+1, z) == x^2+y*z + z^2
        # This needs to be an array due to a compiler bug
        @test homogenize([x^2+y+p], z, parameters=[p]) == [x^2+y*z + z^2*p]

        @test homogenize([x^3+p, 1+y], z, parameters=[p]) == [x^3+p*z^3, z+y]

        h2 = [a * b * c  + p] ∘ [x+y, y + z, x + q]
        @test validate(homogenize(h2, parameters=[p, q]), parameters=[p, q])
        h3 = [a * b * c  + 1] ∘ [x+y, y + z, x + 1]
        @test validate(homogenize(h3))

        # Weylnorm
        @polyvar x y z
        f = 3.0x^2 + 2x*y - y^2
        g = (-2.5+2im) * x^2 - 3.0*x*y + 4y^2
        @test HC.weyldot(f, g) == 3.0 * conj(-2.5 + 2im) + 2.0 * (-3.0) / 2 + (-1.0) * 4.0
        @test HC.weyldot(f, f) == 9.0 + 4.0  / 2 + 1.0
        @test HC.weylnorm(f)^2 ≈ HC.weyldot(f,f)
        @test HC.weyldot([f, f], [g, g]) == 2 * HC.weyldot(f, g)
        @test HC.weylnorm([f, f]) == √HC.weyldot([f, f], [f, f])

        # linear system
        @polyvar x y z
        f = [2x+3y+z+5, -1x+2y+4z-2]
        A, b = linear_system(f)
        @test A == [2 3 1; -1 2 4]
        @test b == [-5, 2]
        @test A * [x,y,z] - b == f
    end

    @testset "Multihomogeneous" begin
        @polyvar x y
        f = [x*y-2, x^2-4]

        affine_hominfo = HC.HomogenizationInformation(variable_groups = ((x,), (y,)))
        @test HC.is_homogeneous(f, affine_hominfo) == false

        @polyvar v w
        g = [x * y - 2v * w, x^2-4v^2]
        hominfo = HC.HomogenizationInformation(variable_groups = ((x,v), (y,w)))
        @test HC.is_homogeneous(g, hominfo) == true

        hominfo = HC.HomogenizationInformation(homvars=(v,w), variable_groups = ((x,v), (y,w)))
        @test HC.is_homogeneous(g, hominfo) == true

        @test HC.homogenize(f, hominfo) == g
        HC.multidegrees(f, ([x], [y])) == [1 2; 1 0]
    end

    @testset "Linear algebra" begin
        # square LU
        A = rand(Complex{Float64}, 12, 12)
        b = rand(Complex{Float64}, 12)
        x = zeros(ComplexF64, 12)
        jac = HC.Jacobian(A)
        HC.update!(jac, A)
        C, d = copy(A), copy(b)
        HC.solve!(x, jac, b)
        @test norm(x - C \ d) < 1e-10

        # square QR / update infos (full rank)
        HC.update!(jac, A; update_infos=true)
        @test jac.corank == 0
        @test jac.cond > 1
        HC.solve!(x, jac, b; update_digits_lost=false)
        @test norm(x - C \ d) < 1e-10

        x .= 0
        HC.solve!(x, jac, b; update_digits_lost=true)
        @test norm(x - C \ d) < 1e-10

        # square QR / not full rank
        A = rand(Complex{Float64}, 12, 12)
        b = rand(Complex{Float64}, 12)
        x = zeros(ComplexF64, 12)
        jac = HC.Jacobian(A)

        # make a corank 1 matrix
        U, s, V = svd(A)
        s[end] = 1e-16
        Â = U * diagm(0 => s) * V'
        HC.update!(jac, Â; update_infos=true)
        @test jac.corank_proposal == 1
        HC.update_rank!(jac)
        @test jac.corank == 1
        @test jac.cond > 1e14
        HC.solve!(x, jac, copy(b); update_digits_lost=false)
        @test norm(x - qr(Â ./ jac.D, Val(true)) \ (b ./ jac.D)) < 1e-10

        x .= 0
        HC.solve!(x, jac, copy(b); update_digits_lost=true)
        @test norm(x - qr(Â ./ jac.D, Val(true)) \ (b ./ jac.D)) < 1e-10
        @test jac.digits_lost == nothing

        # make a corank 1 matrix with not super clear rank
        U, s, V = svd(A)
        s[end] = 1e-15
        Â = U * diagm(0 => s) * V'
        HC.update!(jac, Â; update_infos=true)
        @test jac.corank_proposal == 1
        HC.update_rank!(jac)
        @test jac.corank == 1
        @test jac.cond > 1e13
        HC.solve!(x, jac, copy(b); update_digits_lost=false)
        @test norm(x - qr(Â ./ jac.D, Val(true)) \ (b ./ jac.D)) < 1e-10

        # make a corank 2 matrix
        U, s, V = svd(A)
        s[end-1:end] .= 1e-16
        Â = U * diagm(0 => s) * V'
        HC.update!(jac, Â; update_infos=true)
        @test jac.corank_proposal == 2
        HC.update_rank!(jac)
        @test jac.corank == 2
        @test jac.cond > 1e14
        HC.solve!(x, jac, copy(b); update_digits_lost=false)
        @test norm(x - qr(Â ./ jac.D, Val(true)) \ (b ./ jac.D)) < 1e-10
        x .= 0
        HC.solve!(x, jac, copy(b); update_digits_lost=true)
        @test norm(x - qr(Â ./ jac.D, Val(true)) \ (b ./ jac.D)) < 1e-10
        @test jac.digits_lost == nothing

        # overdetermined
        A = rand(Complex{Float64}, 15, 12)
        b = rand(Complex{Float64}, 15)
        C, d = copy(A), copy(b)
        jac = HC.Jacobian(A)
        HC.update!(jac, A)
        jac.D .= 1
        x = zeros(ComplexF64, 12)
        HC.solve!(x, jac, b)
        @test norm(x - qr(C, Val(true)) \ d) < 1e-10

        x .= 0
        HC.solve!(x, jac, b; update_digits_lost=true)
        @test norm(x - qr(C, Val(true)) \ d) < 1e-10
        @test jac.digits_lost == nothing

        # make a corank 1 matrix
        U, s, V = svd(A)
        s[end] = 1e-15
        Â = U * diagm(0 => s) * V'
        HC.update!(jac, Â; update_infos=true)
        @test jac.corank_proposal == 1
        HC.update_rank!(jac)
        @test jac.corank == 1
        @test jac.cond > 1e14
        HC.solve!(x, jac, copy(b); update_digits_lost=false)
        @test norm(x - qr(Â ./ jac.D, Val(true)) \ (b ./ jac.D)) < 1e-10

        x .= 0
        HC.solve!(x, jac, copy(b); update_digits_lost=true)
        @test norm(x - qr(Â ./ jac.D, Val(true)) \ (b ./ jac.D)) < 1e-10
        @test jac.digits_lost == nothing
    end

    @testset "Misc" begin
        @test HC.unpack(5, 2) == 5
        @test HC.unpack(nothing, 2) == 2

        x = rand()
        @test HC.nthroot(x, 1) == x
        @test HC.nthroot(x, 0) == one(x)

        x = rand(6)
        @test HC.is_real_vector(x)

        segment = HC.ComplexSegment(2, 4)
        @test length(segment) ≈ 2
        test_show_juno(segment)
        @test string(segment) == "ComplexSegment(2.0 + 0.0im, 4.0 + 0.0im)"
    end

    @testset "Hermite Normal Form" begin
        A = [0 3 1 0; -2 2 -1 2; 1 -1 2 3; -3 3 3 2]
        H, U = HC.hnf(A)
        @test A * U == H

        D = zeros(Int, 1,1)
        D[1,1] = -2
        H, U = HC.hnf(D)
        @test H[1,1] == 2
        @test U[1,1] == -1
    end
end


text:
function square(x::ArbFloat{P}; roundingmode::RoundingMode=RoundNearest) where {P}
    z = ArbFloat{P}()
    rounding = match_rounding_mode(roundingmode)
    ccall(@libarb(arf_mul_rnd_any), Cvoid, (Ref{ArbFloat}, Ref{ArbFloat}, Ref{ArbFloat}, Clong, Cint), z, x, x, P, rounding)
    return z
end

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

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

function cube(x::ArbFloat{P}) where {P}
    x1 = ArbReal{P}(x)
    w = square(x1)
    w *= x1
    z = midpoint_byref(w)
    return z
end

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

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

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

function sqrt(x::ArbReal{P}, roundingmode::RoundingMode=RoundNearest) where {P}
    rounding = match_rounding_mode(roundingmode)
    z = ArbReal{P}()
    ccall(@libarb(arb_sqrt), Cvoid, (Ref{ArbReal}, Ref{ArbReal}, Clong, Cint), z, x, P, rounding)
    return z
end

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

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

function rsqrt(x::ArbReal{P}, roundingmode::RoundingMode=RoundNearest) where {P}
    rounding = match_rounding_mode(roundingmode)
    z = ArbReal{P}()
    ccall(@libarb(arb_rsqrt), Cvoid, (Ref{ArbReal}, Ref{ArbReal}, Clong, Cint), z, x, P, rounding)
    return z
end

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

cbrt(x::ArbFloat{P}) where {P} = root(x, 3)
cbrt(x::ArbReal{P}) where {P} = root(x, 3)
cbrt(x::ArbComplex{P}) where {P} = root(x, 3)

rcbrt(x::ArbFloat{P}) where {P} = root(1/x, 3)
rcbrt(x::ArbReal{P}) where {P} = root(1/x, 3)
rcbrt(x::ArbComplex{P}) where {P} = root(1/x, 3)

function hypot(x::ArbFloat{P}, y::ArbFloat{P}) where {P}
    x1 = ArbReal{P}(x)
    y1 = ArbReal{P}(y)
    w = hypot(x1, y1)
    z = midpoint_byref(w)
    return z
end

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

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

function pow(x::ArbFloat{P}, y::ArbFloat{P}) where {P}
    x1 = ArbReal{P}(x)
    y1 = ArbReal{P}(y)
    w = pow(x1, y1)
    z = midpoint_byref(w)
    return z
end

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

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

function pow(x::ArbFloat{P}, y::T) where {P, T<:Integer}
    x1 = ArbReal{P}(x)
    w = pow(x1, y)
    z = midpoint_byref(w)
    return z
end

function pow(x::ArbReal{P}, y::T) where {P, T<:Integer}
    z = ArbReal{P}()
    u = Culong(y)
    ccall(@libarb(arb_pow_ui), Cvoid, (Ref{ArbReal}, Ref{ArbReal}, Culong, Clong), z, x, u, P)
    return z
end

function pow(x::ArbComplex{P}, y::T) where {P, T<:Integer}
    z = ArbComplex{P}()
    u = Culong(y)
    ccall(@libarb(acb_pow_ui), Cvoid, (Ref{ArbComplex}, Ref{ArbComplex}, Culong, Clong), z, x, u, P)
    return z
end

pow(x::ArbFloat{P}, y::T) where {P, T<:Union{Integer,AbstractFloat}} = pow(x, ArbFloat{P}(y))
pow(x::ArbReal{P}, y::T) where {P, T<:Union{Integer,AbstractFloat}} = pow(x, ArbReal{P}(y))
pow(x::ArbComplex{P}, y::T) where {P, T<:Union{Integer,AbstractFloat}} = pow(x, ArbComplex{P}(y))

function root(x::ArbFloat{P}, y::T) where {P, T<:Integer}
    y < 0 && return pow(x, abs(y))
    x1 = ArbReal{P}(x)
    w = root(x1, y)
    z = midpoint_byref(w)
    return z
end
root(x::ArbFloat{P}, y::T) where {P, T<:AbstractFloat} = pow(x, -y)

function root(x::ArbReal{P}, y::T) where {P, T<:Integer}
    y < 0 && return pow(x, abs(y))
    z = ArbReal{P}()
    u = Culong(y)
    ccall(@libarb(arb_root_ui), Cvoid, (Ref{ArbReal}, Ref{ArbReal}, Culong, Clong), z, x, u, P)
    return z
end
root(x::ArbReal{P}, y::T) where {P, T<:AbstractFloat} = pow(x, -y)

function root(x::ArbComplex{P}, y::T) where {P, T<:Integer}
    y < 0 && return pow(x, abs(y))
    z = ArbComplex{P}()
    u = Culong(y)
    ccall(@libarb(acb_root_ui), Cvoid, (Ref{ArbReal}, Ref{ArbReal}, Culong, Clong), z, x, u, P)
    return z
end
root(x::ArbComplex{P}, y::T) where {P, T<:AbstractFloat} = pow(x, -y)

function factorial(x::Signed) where {P}
    z = ArbReal{P}()
    u = Culong(x)
    ccall(@libarb(arb_fac_ui), Cvoid, (Ref{ArbReal}, Culong, Clong), z, u, P)
    return z
end

function factorial(x::ArbFloat{P}) where {P}
    z = ArbReal{P}()
    u = Culong(Clong(x))
    ccall(@libarb(arb_fac_ui), Cvoid, (Ref{ArbReal}, Culong, Clong), z, u, P)
    return midpoint_byref(z)
end

function factorial(x::ArbReal{P}) where {P}
    z = ArbReal{P}()
    u = Culong(Clong(x))
    ccall(@libarb(arb_fac_ui), Cvoid, (Ref{ArbReal}, Culong, Clong), z, u, P)
    return z
end

function doublefactorial(x::Signed) where {P}
    z = ArbReal{P}()
    u = Culong(x)
    ccall(@libarb(arb_doublefac_ui), Cvoid, (Ref{ArbReal}, Culong, Clong), z, u, P)
    return z
end

function doublefactorial(x::ArbFloat{P}) where {P}
    z = ArbReal{P}()
    u = Culong(Clong(x))
    ccall(@libarb(arb_doublefac_ui), Cvoid, (Ref{ArbReal}, Culong, Clong), z, u, P)
    return midpoint_byref(z)
end

function doublefactorial(x::ArbReal{P}) where {P}
    z = ArbReal{P}()
    u = Culong(Clong(x))
    ccall(@libarb(arb_doublefac_ui), Cvoid, (Ref{ArbReal}, Culong, Clong), z, u, P)
    return z
end

function risingfactorial(x::Signed, n::Signed) where {P}
    z = ArbReal{P}()
    ux = Culong(x)
    un = Culong(n)
    ccall(@libarb(arb_rising_ui), Cvoid, (Ref{ArbReal}, Culong, Culong, Clong), z, ux, un, P)
    return z
end

function risingfactorial(x::ArbFloat{P}, n::ArbFloat{P}) where {P}
    z = ArbReal{P}()
    ux = Culong(Clong(x))
    un = Culong(Clong(n))
    ccall(@libarb(arb_rising_ui), Cvoid, (Ref{ArbReal}, Culong, Culong, Clong), z, ux, un, P)
    return midpoint_byref(z)
end

function risingfactorial(x::ArbReal{P}, n::ArbFloat{P}) where {P}
    z = ArbReal{P}()
    u = Culong(Clong(x))
    ccall(@libarb(arb_fac_ui), Cvoid, (Ref{ArbReal}, Ref{ArbReal}, Ref{ArbReal}, Clong), z, x, n, P)
    return z
end

function binomial(n::ArbFloat{P}, k::ArbFloat{P}) where {P}
    z = ArbReal{P}()
    un = Culong(Clong(n))
    uk = Culong(Clong(k))
    ccall(@libarb(arb_bin_uiui), Cvoid, (Ref{ArbReal}, Culong, Culong, Clong), z, ux, un, P)
    return midpoint_byref(z)
end

function binomial(n::ArbReal{P}, k::ArbReal{P}) where {P}
    z = ArbReal{P}()
    un = Culong(Clong(n))
    uk = Culong(Clong(k))
    ccall(@libarb(arb_bin_uiui), Cvoid, (Ref{ArbReal}, Culong, Culong, Clong), z, ux, un, P)
    return z
end


text:
module Morbit

export morbit_formatter;
export MixedMOP, optimize, AlgorithmConfig, add_objective!, add_vector_objective!;
export ExactConfig, TaylorConfig, RbfConfig, LagrangeConfig;
export save_config, save_database, save_iter_data;
export load_config, load_database, load_iter_data;

using Printf: @sprintf

# steepest descent and directed search
using LinearAlgebra: norm, pinv
import JuMP;
import OSQP;

# LagrangeModels and PS
import NLopt;

using Parameters: @with_kw, @unpack, @pack!
using MathOptInterface;
const MOI = MathOptInterface;

import ThreadSafeDicts: ThreadSafeDict;
using Memoize;

import UUIDs;

import FiniteDiff#erences
const FD = FiniteDiff#erences

import ForwardDiff
const AD = ForwardDiff

import Logging: LogLevel, @logmsg
import Logging

include("custom_logging.jl")

include("shorthands.jl");
include("Interfaces.jl");
include("diff_wrappers.jl");

# implementations (order should not matter)
include("VectorObjectiveFunction.jl");
include("MixedMOP.jl");

include("ResultImplementation.jl")
include("DataBaseImplementation.jl")
include("IterDataImplementation.jl")
include("SurrogatesImplementation.jl");

include("ConfigImplementations.jl")

# utilities
include("adding_objectives.jl");
include("descent.jl")
include("saving.jl")
include("utilities.jl")

function shrink_radius( ac :: AbstractConfig, Δ :: NumOrVec, steplength :: NumOrVec)
    if radius_update_method(ac) == :standard
        return Δ .* _gamma_shrink(ac);
    elseif radius_update_method(ac) == :steplength
        return steplength .* _gamma_shrink(ac);
    end
end

function shrink_radius_much( ac :: AbstractConfig, Δ :: NumOrVec, steplength :: NumOrVec)
    if radius_update_method(ac) == :standard
        return Δ .* _gamma_shrink_much(ac);
    elseif radius_update_method(ac) == :steplength
        return steplength .* _gamma_shrink_much(ac);
    end
end

function grow_radius( ac :: AbstractConfig, Δ :: NumOrVec, steplength :: NumOrVec)
    if radius_update_method(ac) == :standard
        return min( Δᵘ(ac), _gamma_grow(ac) .* Δ )
    elseif radius_update_method(ac) == :steplength
        return min( Δᵘ(ac), ( _gamma_grow .+ steplength ./ Δ ) .* Δ );
    end
end

"""
Perform initialization of the data passed to `optimize` function.
"""
function initialize_data( mop :: AbstractMOP, x0 :: RVec, fx0 :: Vec = Float32[];
    algo_config :: Union{AbstractConfig, Nothing} = nothing,
    populated_db :: Union{AbstractDB, Nothing} = nothing )

    if num_objectives(mop) == 0
        error("`mop` has no objectives!")
    end

    @warn "The evaluation counter of `mop` is reset."
    reset_evals!( mop );
    # for backwards-compatibility with unconstrained problems:
    if num_vars(mop) == 0
        MOI.add_variables(mop, length(x0))
    end

    # initialize first iteration site
    @assert !isempty( x0 ) "Please provide a non-empty feasible starting point `x0`."

    x_scaled = scale( x0, mop );
    tfn = TransformerFn(mop)

    # initalize first objective vector
    if isempty( fx0 )
        # if no starting function value was provided, eval objectives
        fx_sorted = eval_all_objectives( mop, tfn, x_scaled );
    else
        fx_sorted = apply_internal_sorting( fx0, mop );
    end

    # ensure at least half-precision
    F = Base.promote_eltype( x_scaled, fx_sorted, Float32 )
    x = F.(x_scaled)
    fx = F.(fx_sorted)

    if isnothing( algo_config )
        ac = DefaultConfig{F}()
    else
        ac = algo_config
    end

    # initialize database
    if !isnothing(populated_db)
        # has a database been provided? if yes, prepare
        data_base = populated_db;
        transform!( data_base, mop );
    else
        data_base = init_db( use_db(ac) );
        set_transformed!(data_base, true)
    end

    # make sure, x & fx are in database
    start_res = init_res(Result{F}, x, fx )
    ensure_contains_result!(data_base, start_res);

    iter_data = init_iter_data(IterData, x, fx, Δ⁰(ac));
    xᵗ_index!(iter_data, get_id(start_res));

    # initialize surrogate models
    sc = init_surrogates( mop, iter_data, ac );

    return (mop, tfn, ac, iter_data, sc)
end

function iterate!( iter_data :: AbstractIterData, mop :: AbstractMOP, tfn :: TransformerFn,
    sc :: SurrogateContainer, algo_config :: AbstractConfig )

    # config params
    ν_succ = _nu_success( algo_config );
    ν_acc = _nu_accept( algo_config );
    μ = _mu( algo_config );
    _beta = _beta( algo_config );
    _eps_crit = _eps_crit( algo_config );
    _gamma_crit = _gamma_crit( algo_config );
    count_nonlin = count_nonlinear_iterations( algo_config );

    # read iter data to handy variables
    x = xᵗ(iter_data);
    fx = fxᵗ(iter_data);
    Δ = Δᵗ(iter_data);

    # check some stopping conditions
    if num_iterations(iter_data) >= max_iter(algo_config)
        @logmsg loglevel1 "Stopping. Maximum number of iterations reached."
        return true;
    end

    if !_budget_okay(mop, algo_config)
        @logmsg loglevel1 "Stopping. Computational budget is exhausted."
        return true; # TODO stop codes
    end

    if Δ_abs_test( Δ, algo_config )
        return true;
    end

    # set iteration counter
    if it_stat(iter_data) != MODELIMPROVING || count_nonlin
        inc_iterations!(iter_data)
        set_model_improvements!(iter_data, 0);
    else
        inc_model_improvements!(iter_data)
    end

    @logmsg loglevel1 """\n
        |--------------------------------------------
        |Iteration $(num_iterations(iter_data)).$(num_model_improvements(iter_data))
        |--------------------------------------------
        |  Current trust region radius is $(Δ).
        |  Current number of function evals is $(num_evals(mop)).
        |  Iterate is $(_prettify(unscale(x, mop)))
        |  Values are $(_prettify(reverse_internal_sorting(fx,mop)))
        |--------------------------------------------
    """

    # update surrogate models
    if num_iterations(iter_data) > 1
        if it_stat == MODELIMPROVING
            improve_surrogates!( sc, mop, iter_data, algo_config; ensure_fully_linear = false );
        else
            update_surrogates!( sc, mop, iter_data, algo_config; ensure_fully_linear = false );
        end
    end

    # calculate descent step and criticality value
    ω, x₊, mx₊, steplength = compute_descent_step( algo_config, mop, iter_data, sc )
    @logmsg loglevel1 "Criticality is ω = $(ω)."

    # stop at the end of the this loop run?
    if ω_Δ_rel_test(ω, Δ, algo_config) || ω_abs_test( ω, algo_config )
        return true;
    end

    # Criticallity test
    _fully_linear = fully_linear(sc)
    num_critical_loops = 0;

    if ω <= _eps_crit && (!_fully_linear || all(Δ .> mu * ω))
        @logmsg loglevel1 "Entered Criticallity Test."
        if !_fully_linear
            @logmsg loglevel1 "Ensuring all models to be fully linear."
            update_surrogates!( sc, mop, iter_data, algo_config; ensure_fully_linear = true );

            ω, x₊, mx₊, steplength = compute_descent_step(algo_config,mop,iter_data,sc);
            if !fully_linear(sc)
                @logmsg loglevel1 "Could not make all models fully linear. Trying one last descent step."
                @goto MAIN;
            end
        end

        while all(Δᵗ(iter_data) .> mu * ω)
            @logmsg loglevel1 "Criticality loop $(num_critical_loops + 1)."
            if num_critical_loops >= max_critical_loops(algo_config)
                @logmsg loglevel1 "Maximum number ($(max_critical_loops(algo_config))) of critical loops reached. Exiting..."
                @goto EXIT_MAIN
            end
            if !_budget_okay(mop, algo_config)
                @logmsg loglevel1 "Computational budget exhausted. Exiting…"
                @goto EXIT_MAIN
            end

            # shrink radius
            Δᵗ!( iter_data, Δᵗ(iter_data) .* _gamma_crit );
            # make model linear
            update_surrogates!( sc, mop, iter_data, algo_config; ensure_fully_linear = true );
            # (re)calculate criticality
            # TODO make backtracking optional and don't do here
            ω, data = get_criticality

            if Δ_abs_test( Δᵗ(iter_data) , algo_config ) ||
                ω_Δ_rel_test(ω, Δᵗ(iter_data), algo_config) || ω_abs_test( ω, algo_config )

                @goto EXIT_MAIN
            end

            if !fully_linear(sc)
                @logmsg loglevel1 "Could not make all models fully linear. Trying one last descent step."
                @goto MAIN;
            end

            num_critical_loops += 1;
        end
        "Exiting after $(num_critical_loops) loops with ω = $(ω) and Δ = $(Δᵗ(iter_data))."
        @goto MAIN
        @label EXIT_MAIN
        return true;
    end# Crit test if

    @label MAIN # re-entry point after criticality test

    mx = eval_models(sc, x);
    fx₊ = eval_all_objectives(mop, tfn, x₊);

    if strict_acceptance_test( algo_config )
        ρ = minimum( (fx .- fx₊) ./ (mx .- mx₊) )
    else
        ρ = (maximum(fx) - maximum( fx₊ ))/(maximum(mx)-maximum(mx₊))
    end

    @logmsg loglevel2 """\n
    Attempting descent of length $steplength.
    | f(x)  | $(_prettify(fx))
    | f(x₊) | $(_prettify(fx₊))
    | m(x)  | $(_prettify(mx))
    | m(x₊) | $(_prettify(mx₊))
    The error betwenn f(x) and m(x) is $(sum(abs.(fx .- mx))).
    $(strict_acceptance_test(algo_config) ? "All" : "One") of the components must decrease.
    Thus, ρ is $ρ.
    """

    # update trust region radius
    ρ = isnan(ρ) ? -Inf : ρ;
    old_Δ = copy(Δᵗ(iter_data));  # if it was changed in criticality test
    if ρ >= ν_succ
        if Δ < _beta * ω
            Δᵗ!(iter_data, grow_radius(algo_config, Δ, steplength) );
        end
        it_stat!(iter_data, SUCCESSFULL)
    else
        if fully_linear(sc)
            if ρ < ν_acc
                Δᵗ!(iter_data, shrink_radius_much(algo_config, Δ, steplength) );
                it_stat!(iter_data, INACCEPTABLE)
            else
                Δᵗ!(iter_data, shrink_radius(algo_config, Δ, steplength) );
                it_stat!(iter_data, ACCEPTABLE)
            end
        else
            it_stat!(iter_data, MODELIMPROVING)
        end
    end

    # accept x?
    old_x_index = xᵗ_index(iter_data); # for stamp!ing
    if it_stat(iter_data) == SUCCESSFULL || it_stat(iter_data) == ACCEPTABLE
        trial_point_index = set_next_iterate!(iter_data, x₊, fx₊);
    else
        trial_point_index = keep_current_iterate!(iter_data, x₊, fx₊);
    end
    @logmsg loglevel1 """\n
        The trial point was $( (it_stat(iter_data) == SUCCESSFULL) || (it_stat(iter_data) == ACCEPTABLE ) ? "" : "not ")accepted.
        The iteration is $(it_stat(iter_data)).
        Moreover, the radius was updated as below:
        old radius : $(old_Δ)
        new radius : $(Δᵗ(iter_data)) ($(round(Δᵗ(iter_data)/old_Δ * 100; digits=1)) %)
    """

    stamp!(db(iter_data), Dict(
            "iter_status" => it_stat(iter_data),
            "xᵗ_index" => old_x_index,
            "xᵗ₊_index" => trial_point_index,
            "ρ" => ρ,
            "Δ" => Δᵗ(iter_data),
            "ω" => ω,
            "num_critical_loops" => num_critical_loops,
            "stepsize" => steplength,
            "model_meta" => [ deepcopy(sw.meta) for sw ∈ sc.surrogates ]
        )
    );

    if x_tol_rel_test( x, x₊, algo_config  ) ||
        x_tol_abs_test( x, x₊, algo_config ) ||
        f_tol_rel_test( fx, fx₊, algo_config  ) ||
        f_tol_abs_test( fx, fx₊, algo_config )
        return true;
    end
    return false
end

function get_return_values( iter_data :: AbstractIterData, mop :: AbstractMOP)
    ret_x = unscale(xᵗ(iter_data),mop);
    ret_fx = reverse_internal_sorting(fxᵗ(iter_data),mop);
    # TODO make tests for this: are returen values the same before and after `finalize_iter_data!` ?
    return (ret_x, ret_fx)
end

function finalize_iter_data!( iter_data :: AbstractIterData, mop :: AbstractMOP )
    if !isnothing( db(iter_data) )
        untransform!(db(iter_data), mop)
    end
    nothing
end

############################################
function optimize( mop :: AbstractMOP, x0 :: Vec, fx0 :: Vec = Float16[];
    algo_config :: Union{Nothing, AbstractConfig} = nothing,
    populated_db :: Union{AbstractDB, Nothing} = nothing # TODO make passing of AbstractIterData possible
    )

    mop, tfn, ac, iter_data, sc = initialize_data( mop, x0, fx0; algo_config, populated_db )

    @logmsg loglevel1 _stop_info_str( ac, mop )

    @logmsg loglevel1 "Entering main optimization loop."
    while true
        abbort = iterate!(iter_data, mop, tfn, sc, ac)
        abbort && break;
    end# while

    ret_x, ret_fx = get_return_values( iter_data, mop )
    # unscale sites and re-sort values to return to user
    finalize_iter_data!(iter_data, mop)

    # @assert all((ret_x, ret_fx) .== get_return_values( iter_data, mop ))

    @logmsg loglevel1 _fin_info_str(iter_data, mop)

    return ret_x, ret_fx, iter_data
end# function optimize

end#module


text:
function test_antidomain_singlevariable(
    model::MOI.ModelLike,
    config::MOIT.Config{T},
) where {T <: Real}
    @MOIT.requires MOI.supports_constraint(model, MOI.VariableIndex, MOI.Integer) # x1, x2
    @MOIT.requires MOI.supports_constraint(model, MOI.VariableIndex, MOI.EqualTo{T}) # c1
    @MOIT.requires MOI.supports_constraint(model, MOI.VariableIndex, MOI.GreaterThan{T}) # c2
    @MOIT.requires MOI.supports_constraint(model, MOI.VariableIndex, MOI.LessThan{T}) # c3
    @MOIT.requires MOI.supports_constraint(model, MOI.VectorOfVariables, CP.AllDifferent) # c4
    @MOIT.requires MOI.supports_constraint(model, MOI.VariableIndex, CP.AntiDomain{T}) # c5

    x1, _ = MOI.add_constrained_variable(model, MOI.Integer())
    x2, _ = MOI.add_constrained_variable(model, MOI.Integer())

    c1 = MOI.add_constraint(model, x1, MOI.EqualTo(T(1)))
    c2 = MOI.add_constraint(model, x2, MOI.GreaterThan(T(1)))
    c3 = MOI.add_constraint(model, x2, MOI.LessThan(T(3)))
    c4 = MOI.add_constraint(model, MOI.VectorOfVariables([x1, x2]), CP.AllDifferent(2))

    c5 = MOI.add_constraint(model, x2, CP.AntiDomain(Set(T[3])))

    @test MOI.is_valid(model, x1)
    @test MOI.is_valid(model, x2)
    @test MOI.is_valid(model, c1)
    @test MOI.is_valid(model, c2)
    @test MOI.is_valid(model, c3)
    @test MOI.is_valid(model, c4)
    @test MOI.is_valid(model, c5)

    if MOIT._supports(config, MOI.optimize!)
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMIZE_NOT_CALLED
        MOI.optimize!(model)
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMAL
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

        @test MOI.get(model, MOI.ResultCount()) >= 1
        @test MOI.get(model, MOI.VariablePrimal(), x1) == 1
        @test MOI.get(model, MOI.VariablePrimal(), x2) == 2
    end
end

function MOIT.setup_test(
    ::typeof(test_antidomain_singlevariable),
    mock::MOIU.MockOptimizer,
    ::Config{T},
) where {T <: Real}
    MOIU.set_mock_optimize!(
        mock,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, T[1, 2]),
    )
    return
end

function test_antidomain_scalaraffinefunction(
    model::MOI.ModelLike,
    config::MOIT.Config{T},
) where {T <: Real}
    @MOIT.requires MOI.supports_constraint(model, MOI.VariableIndex, MOI.Integer) # x1, x2
    @MOIT.requires MOI.supports_constraint(model, MOI.VariableIndex, MOI.EqualTo{Int}) # c1
    @MOIT.requires MOI.supports_constraint(model, MOI.VariableIndex, MOI.GreaterThan{Int}) # c2
    @MOIT.requires MOI.supports_constraint(model, MOI.VariableIndex, MOI.LessThan{Int}) # c3
    @MOIT.requires MOI.supports_constraint(model, MOI.ScalarAffineFunction{Int}, CP.AllDifferent) # c4
    @MOIT.requires MOI.supports_constraint(model, MOI.VariableIndex, CP.AntiDomain{Int}) # c5

    x1, _ = MOI.add_constrained_variable(model, MOI.Integer())
    x2, _ = MOI.add_constrained_variable(model, MOI.Integer())

    c1 = MOI.add_constraint(model, x1, MOI.EqualTo(T(1)))
    c2 = MOI.add_constraint(model, x2, MOI.GreaterThan(T(1)))
    c3 = MOI.add_constraint(model, x2, MOI.LessThan(T(3)))
    c4 = MOI.add_constraint(model, MOIU.vectorize(one(T) .* [x1, x2]), CP.AllDifferent(2))

    c5 = MOI.add_constraint(
        model,
        one(T) * x2,
        CP.AntiDomain(Set([3])),
    )

    @test MOI.is_valid(model, x1)
    @test MOI.is_valid(model, x2)
    @test MOI.is_valid(model, c1)
    @test MOI.is_valid(model, c2)
    @test MOI.is_valid(model, c3)
    @test MOI.is_valid(model, c4)
    @test MOI.is_valid(model, c5)

    if MOIT._supports(config, MOI.optimize!)
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMIZE_NOT_CALLED
        MOI.optimize!(model)
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMAL
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

        @test MOI.get(model, MOI.ResultCount()) >= 1
        @test MOI.get(model, MOI.VariablePrimal(), x1) == 1
        @test MOI.get(model, MOI.VariablePrimal(), x2) == 2
    end
end

function MOIT.setup_test(
    ::typeof(test_antidomain_scalaraffinefunction),
    mock::MOIU.MockOptimizer,
    ::Config{T},
) where {T <: Real}
    MOIU.set_mock_optimize!(
        mock,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, T[1, 2]),
    )
    return
end


text:
#    Problem : GROUP A
#    *********
#	 A nondiagonal quartic test problem.
#
#	 This problem has an arrow-head type Hessian with a tridiagonal
#    central part and a border of width 1.
#    The Hessian is singular at the solution.
#
#    Origonal SIF Source:
#    A.R. Conn, N.I.M. Gould, M. Lescrenier and Ph.L. Toint,
#    "Performance of a multi-frontal scheme for partially separable
#    optimization"
#    Report 88/4, Dept of Mathematics, FUNDP (Namur, B), 1988.
#
#    NONDIA.jl.SIF classification OUR2-AN-V-0
#
#    Number of variables is variable
#
# Daniel Henderson, 08/2021

f = (x) -> begin
	fx = (x[1]-1)^2
	for i in 2:lastindex(x)
		fx += 100*(x[1]-x[i-1]^2)^2
	end
    return fx
end

g! = (g, x) -> begin
	g[1] = 2*(x[1]-1)
	for i in 2:lastindex(x)
		g[1]   += 200(x[1]-x[i-1]^2)
		g[i-1] += -400x[i-1]*(x[1] - x[i-1]^2)
	end
    return g
end

fg! = (g, x) -> begin
	fx = (x[1]-1)^2
	g[1] = 2*(x[1]-1)
	for i in 2:lastindex(x)
		fx += 100*(x[1]-x[i-1]^2)^2
		g[1]   += 200(x[1]-x[i-1]^2)
		g[i-1] += -400x[i-1]*(x[1] - x[i-1]^2)
	end
    return g
end

init = (n::Int=5000) -> begin
	n < 2 && @warn("ARWHEAD: number of variables must be ≥ 2")
	n = max(n, 2)

    return n, -1.0ones(n)
end

TestSet["NONDIA"] = UncProgram("NONDIA", f, g!, fg!, init)


text:
module ContingencyTest

using Test
using ForneyLab
using ForneyLab: vague, dims

@testset "Contingency ProbabilityDistribution and Message construction" begin
    @test ProbabilityDistribution(Multivariate, Contingency, p=[0.1 0.4; 0.3 0.2]) == ProbabilityDistribution{Multivariate, Contingency}(Dict(:p=>[0.1 0.4; 0.3 0.2]))
    @test_throws Exception ProbabilityDistribution(Univariate, Contingency)
    @test ProbabilityDistribution(Contingency, p=[0.1 0.4; 0.3 0.2]) == ProbabilityDistribution{Multivariate, Contingency}(Dict(:p=>[0.1 0.4; 0.3 0.2]))
    @test ProbabilityDistribution(Contingency) == ProbabilityDistribution{Multivariate, Contingency}(Dict(:p=>(1/9)*ones(3,3)))
end

@testset "dims" begin
    @test dims(ProbabilityDistribution(Contingency, p=[0.1 0.4; 0.3 0.2])) == 2
end

@testset "vague" begin
    @test vague(Contingency) == ProbabilityDistribution(Contingency, p=(1/9)*ones(3,3))
    @test vague(Contingency, (2,3)) == ProbabilityDistribution(Contingency, p=(1/6)*ones(2,3))
end

@testset "differentialEntropy" begin
    @test differentialEntropy(ProbabilityDistribution(Multivariate, Contingency, p=[0.2 0.1 0.7; 0.4 0.3 0.3; 0.1 0.6 0.3])) == 2.78866425274534
    @test differentialEntropy(ProbabilityDistribution(Multivariate, Contingency, p=[0.3*[0.2 0.1 0.7; 0.4 0.3 0.3; 0.1 0.6 0.3], 0.7*[0.2 0.1 0.7; 0.4 0.3 0.3; 0.1 0.6 0.3]])) == 4.621257158910021
end

end # module


text:
"""
    f∂f(f,x)

returns f and df/fx evaluated in `x`, using `ForwardDiff.jl`, `DiffResults.jl` and `StaticArrays.jl` to calculate f and dfdx in one pass
"""
function f∂f(f,x::T) where T
    _f(z) = f(only(z))
    x_vec =   SVector(x)
    ∂result = DiffResults.GradientResult(x_vec)
    _∂f =  ForwardDiff.gradient!(∂result, _f,x_vec)
    fx =  DiffResults.value(_∂f)
    ∂f∂x = only(DiffResults.gradient(_∂f))
    return fx,fx/∂f∂x
end

"""
    ad_newton(f,x0;kwargs...)

    Performs newton root finding (via Roots.jl), but the derivatives are calculated with ForwardDiff.jl
    kwargs are passed to `Roots.find_zero`
"""
function ad_newton(f,x0;kwargs...)
    function f0(x)
        fx,∂f∂x = f∂f(f,x)
        return fx,fx/∂f∂x
    end
    return Roots.find_zero(f0,x0,Roots.Newton();kwargs...)
end


text:
"""
$(TYPEDEF)
$(TYPEDFIELDS)

"""
struct QuantNormal <: AbstractDistDelayed
    mu::NestedFloat
    sigma::NestedFloat
    q::NestedFloat
end

function quantnormal(mu::Float64, sigma::Float64, q::Float64)::Float64
    return round(normal(mu, sigma) / q) * q
end


text:
module FakeBlock

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Fake Block" => Ahorn.EntityPlacement(
        Maple.FakeBlock,
        "rectangle",
        Dict{String, Any}(),
        Ahorn.tileEntityFinalizer
    ),
)

Ahorn.editingOptions(entity::Maple.FakeBlock) = Dict{String, Any}(
    "tiletype" => Ahorn.tiletypeEditingOptions()
)

Ahorn.minimumSize(entity::Maple.FakeBlock) = 8, 8
Ahorn.resizable(entity::Maple.FakeBlock) = true, true

Ahorn.selection(entity::Maple.FakeBlock) = Ahorn.getEntityRectangle(entity)

Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::Maple.FakeBlock, room::Maple.Room) = Ahorn.drawTileEntity(ctx, room, entity, alpha=0.7)

end


text:
@deprecate param(x) x
@deprecate data(x) x


text:
using FractionalDiffEq
using Test
using MittagLeffler

@testset "FractionalDiffEq.jl" begin
    include("test.jl")
end


text:
# ------------ Call sites ------------ #

@inline function (ctx::SimulateContext)(c::typeof(markov),
                                        addr::Address,
                                        call::Function,
                                        len::Int,
                                        args...)
    visit!(ctx, addr => 1)
    ps = get_subparameters(ctx, addr)
    ret, cl = simulate(ps, call, args...)
    v_ret = Vector{typeof(ret)}(undef, len)
    v_cl = Vector{typeof(cl)}(undef, len)
    v_ret[1] = ret
    v_cl[1] = cl
    for i in 2:len
        visit!(ctx, addr => i)
        ret, cl = simulate(ps, call, v_ret[i-1]...)
        v_ret[i] = ret
        v_cl[i] = cl
    end
    sc = sum(map(v_cl) do cl
                 get_score(cl)
             end)
    add_call!(ctx, addr, VectorizedCallSite{typeof(markov)}(VectorizedTrace(v_cl), sc, call, len, args, v_ret))
    return v_ret
end


text:
"""
    fit!(learner, phases)

Fit `learner` with [`Phase`](#)s `phases`. See [`./docs/training/basics.md`] for more
info on the training loop.
"""
function fit!(learner::Learner, phases::AbstractVector{<:Phase})
    if !(learner.callbacks.initialized)
        handle(Init(), learner, phases[1])
        learner.callbacks.initialized = true
    end
    try
        for phase in phases
            fitepoch!(learner, phase)
        end
    catch e
        if e isa CancelFittingException
            @debug "Fitting was cancelled" error=e
        else
            rethrow()
        end
    end
    return learner
end

fit!(learner::Learner, phase::Phase)::Learner = fit!(learner, [phase])

"""
    fit!(learner, n)

Shorthand for `fit!(learner, repeat([TrainingPhase(), ValidationPhase()], n))`, i.e.
a very basic training loop of `n` epochs of training followed by validation.
"""
fit!(learner, n::Int)::Learner = fit!(learner, repeat([TrainingPhase(), ValidationPhase()], n))

function fitepoch!(learner, phase)
    try
        fitepochphase!(learner, phase)
    catch e
        if e isa CancelEpochException
            @debug "Epoch was cancelled" error=e
        else
            rethrow()
        end
    end

    return learner
end

function fitbatch!(learner, batch, phase)
    learner.batch = BatchState()
    try
        fitbatchphase!(learner, batch, phase)
    catch e
        if e isa CancelBatchException
            @debug "Batch was cancelled" error=e
        else
            rethrow()
        end
    end
end

"""
    fitepoch!(learner, phase)

Fit `learner` for one epoch.

Customize by deriving custom phase from `Phase`.
"""
function fitepochphase!(
        learner::Learner,
        phase::Union{TrainingPhase, ValidationPhase, TestPhase},
    )
    dataiter = getdataiter(phase, learner)
    if dataiter === nothing
        throw(CancelEpochException("No data found for phase $(typeof(phase))"))
    end

    handle(EpochBegin(), learner, phase)

    for batch in dataiter
        fitbatch!(learner, batch, phase)
    end

    handle(EpochEnd(), learner, phase)
end

function fitbatchphase!(
        learner::Learner,
        batch,
        phase::AbstractTrainingPhase,
    )

    b = learner.batch
    b.xs, b.ys = batch

    handle(BatchBegin(), learner, phase)

    b.grads = gradient(learner.params) do
        b.ŷs = learner.model(b.xs)

        handle(LossBegin(), learner, phase)
        b.loss = learner.lossfn(b.ŷs, b.ys)

        handle(BackwardBegin(), learner, phase)
        return b.loss
    end
    handle(BackwardEnd(), learner, phase)

    update!(learner.optimizer, learner.params, b.grads)

    handle(BatchEnd(), learner, phase)
    return learner
end

function fitbatchphase!(
        learner::Learner,
        batch,
        phase::ValidationPhase,
    )
    b = learner.batch
    b.xs, b.ys = batch
    handle(BatchBegin(), learner, phase)

    b.ŷs = learner.model(b.xs)

    handle(LossBegin(), learner, phase)
    b.loss = learner.lossfn(b.ŷs, b.ys)

    handle(BatchEnd(), learner, phase)
end


text:
# Based on https://arxiv.org/abs/1409.0473
include("0-data.jl")
using Flux: flip, logitcrossentropy, reset!, throttle
using Parameters: @with_kw
using StatsBase: wsample

@with_kw mutable struct Args
    lr::Float64 = 1e-3      # learning rate
    Nin::Int = 0            # size of input layer, will be assigned as length(alphabet)
    Nh::Int = 30            # size of hidden layer
    phones_len::Int = 0     # length of phonemes
    throttle::Int = 30      # throttle timeout
end

function build_model(args)
    # A recurrent model which takes a token and returns a context-dependent
    # annotation.
    forward  = LSTM(args.Nin, args.Nh÷2)
    backward = LSTM(args.Nin, args.Nh÷2)
    encode(tokens) = vcat.(forward.(tokens), flip(backward, tokens))

    alignnet = Dense(2*args.Nh, 1)

    # A recurrent model which takes a sequence of annotations, attends, and returns
    # a predicted output token.
    recur   = LSTM(args.Nh+args.phones_len, args.Nh)
    toalpha = Dense(args.Nh, args.phones_len)
    return (forward, backward, alignnet, recur, toalpha), encode
end

align(s, t, alignnet) = alignnet(vcat(t, s .* Int.(ones(1, size(t, 2)))))

function asoftmax(xs)
  xs = [exp.(x) for x in xs]
  s = sum(xs)
  return [x ./ s for x in xs]
end

function decode1(tokens, phone, state)
    # Unpack models
    forward, backward, alignnet, recur, toalpha = state
    weights = asoftmax([align(recur.state[2], t, alignnet) for t in tokens])
    context = sum(map((a, b) -> a .* b, weights, tokens))
    y = recur(vcat(Float32.(phone), context))
    return toalpha(y)
end

decode(tokens, phones, state) = [decode1(tokens, phone, state) for phone in phones]

function model(x, y, state, encode)
    # Unpack models
    forward, backward, alignnet, recur, toalpha = state
    ŷ = decode(encode(x), y, state)
    reset!(state)
    return ŷ
end

function predict(s, state, encode, alphabet, phones)
    ts = encode(tokenise(s, alphabet))
    ps = Any[:start]
    for i = 1:50
      dist = softmax(decode1(ts, onehot(ps[end], phones), state))
      next = wsample(phones, vec(dist))
      next == :end && break
      push!(ps, next)
    end
    reset!(state)
    return ps[2:end]
end

function train(; kws...)
    # Initialize Hyperparameters
    args = Args(; kws...)
    @info("Loading Data...")
    data,alphabet,phones = getdata(args)

    # The full model
    # state = (forward, backward, alignnet, recur, toalpha)
    @info("Constructing Model...")
    state, encode = build_model(args)

    loss(x, yo, y) = sum(logitcrossentropy.(model(x, yo, state, encode), y))
    evalcb = () -> @show loss(data[500]...)
    opt = ADAM(args.lr)
    @info("Training...")
    Flux.train!(loss, params(state), data, opt, cb = throttle(evalcb, args.throttle))
    return state, encode, alphabet, phones
end

cd(@__DIR__)
state, encode, alphabet, phones = train()
@info("Testing...")
predict("PHYLOGENY", state, encode, alphabet, phones)

max_stars_repo_name: shawn-kern/model-zoo

max_issues_repo_name: shawn-kern/model-zoo

max_forks_repo_name: shawn-kern/model-zoo


text:
# 'cstrt_reinhardt_ngon' provides the vertices coordinates (a,b) of a
# Reinhardt (m*k)-gon
function cstrt_reinhardt_ngon(m,k)
if mod(m,2) == 1 && m>=3 && k>=1
    a = zeros(m*k-1); b = zeros(m*k-1);
    for i = 1:m-1
        a[i*k] = sin(2*pi*i/m)/(2*cos(pi/(2*m)));
        b[i*k] = (1 - cos(2*pi*i/m))/(2*cos(pi/(2*m)));
    end
    if k>1
        for j = 1:k-1
            a[Int(k*(m-1)/2+j)] = cos((1+(m-3)/2+j/k)*(pi/m));
            b[Int(k*(m-1)/2+j)] = sin((1+(m-3)/2+j/k)*(pi/m));
            for i = 1:Int((m-1)/2)
                a[Int((i+(m-1)/2)*k+j)] = a[i*k] + cos((2*i+1+(m-3)/2+j/k)*(pi/m));
                b[Int((i+(m-1)/2)*k+j)] = b[i*k] + sin((2*i+1+(m-3)/2+j/k)*(pi/m));
                a[Int(m*k-((i+(m-1)/2)*k+j))] = -a[Int((i+(m-1)/2)*k+j)];
                b[Int(m*k-((i+(m-1)/2)*k+j))] = b[Int((i+(m-1)/2)*k+j)];
            end
        end
    end
    return(a,b)
end
end

max_stars_repo_name: cbingane/optigon

max_issues_repo_name: cbingane/optigon

max_forks_repo_name: cbingane/optigon


text:
using MPI
using Oceananigans.Fields: location
import Oceananigans.Fields: Field
import Oceananigans.Grids: AbstractGrid, size

function Field(X, Y, Z, arch::AbstractMultiArchitecture, grid::AbstractGrid,
                bcs = FieldBoundaryConditions(grid, (X, Y, Z)),
               data = new_data(eltype(grid), arch, grid, (X, Y, Z)))

    communicative_bcs = inject_halo_communication_boundary_conditions(bcs, arch.local_rank, arch.connectivity)

    return Field(X, Y, Z, child_architecture(arch), grid, communicative_bcs, data)
end


text:
using Documenter, QuantumInformation.ShortNames

makedocs(
    format=:html,
    sitename="QuantumInformation.jl",
    modules=[QuantumInformation],
    pages=["index.md", "example.md", "latex.md", "api.md"]
)

deploydocs(
    repo="github.com/goropikari/QuantumInformation.jl.git",
    julia="0.6",
    target="build",
    osname = "linux",
    deps=nothing,
    make=nothing
)


text:
struct Quaternion{T<:Number}
    v::Vector3{T}
    w::T

    Quaternion{T}(v, w) where T = new{T}(v, w)

    Quaternion{T}() where T = new{T}(Vector3{T}(), one(T))
end

Quaternion(v::Vector3{U}, w::T) where {T,U} = Quaternion{promote_type(T,U)}(v, w)