text:
module SDMX

using JSON3, OrderedCollections
using Tables, PrettyTables

#--------------------------------------------------------------------
# Tables.jl implementation
#--------------------------------------------------------------------

struct Datatable <: Tables.AbstractRow #TODO: add attributes
    name::String
    headers::Vector{Symbol}
    data::Vector{Any}
    dimensions::OrderedDict{String,Vector{String}}
    attributes::OrderedDict{String,Vector{String}}
end

# Declare that SDMXdata is a table
Tables.istable(::Type{<:Datatable}) = true

# getter methods to avoud getproperty clash
name(dt::Datatable) = getfield(dt, :name)
headers(dt::Datatable) = getfield(dt, :headers)
data(dt::Datatable) = getfield(dt, :data)
dimensions(dt::Datatable) = getfield(dt, :dimensions)
attributes(dt::Datatable) = getfield(dt, :attributes)

# Tables.rows implementation. Fallback definitions are valid
Tables.rowaccess(::Type{<:Datatable}) = true
Tables.rows(dt::Datatable) = data(dt)

# Complete Abstractrow interface
Tables.getcolumn(dt::Datatable, i::Int) = [row[i] for row in data(dt)]
Tables.getcolumn(dt::Datatable, nm::Symbol) = [row[nm] for row in data(dt)]
Tables.columnnames(dt::Datatable) = headers(dt)

#--------------------------------------------------------------------
# Datatable prettyprinting
#--------------------------------------------------------------------

function Base.show(io::IO, dt::Datatable)
    println("\n ", length(headers(dt)), "x", length(data(dt)), " SDMX.Datatable")
    printstyled(" ", name(dt), "\n"; bold = true)
    pretty_table(dt; alignment = :l)
end

#--------------------------------------------------------------------
# Datatable generation
#--------------------------------------------------------------------

"""
    SDMX.read(js::Union{Vector{UInt8},String}; alldims::Bool = true)

Constructs SDMX.datatable compatible with Tables.jl

If 'alldims' is set to false only dimensions with more than one value or a specific role are included.

# Example
```julia
    HTTP.get("https://dataurl").body |> SDMX.read```

"""
# TODOs: select specific dimensions to include, select specific attributes to include, option for unpivoted series
function read(js::Union{Vector{UInt8},String};
            alldims::Union{Bool, Array{Symbol}} = true, unpivotseries::Bool = false)

    ds = JSON3.read(js)

    # Dimensions to include in Datatable
    dims = OrderedDict() # Collects 'Series' (if any) and 'Observations' dimensions
    mask = Dict()
    for (k, v) in ds.structure.dimensions
        if alldims
            dims[k] = v
            mask[k] = Colon()
        else
            m = [(length(i["values"]) > 1) | (haskey(i, "role")) for i in v] .== true # BitVector vector for filtering
            dims[k] = v[m]
            mask[k] = m
        end
    end

    # Headers
    headers = Vector{Symbol}()
    for v in values(dims)
        for i in v
            push!(headers, Symbol(i["name"]))
        end
    end
    push!(headers, :Value)

    # Rows of observations
    obs = []
    if haskey(dims, :series) # FIXME: this creates the tabla in "long" format. Unnecessary as can be selected in the query
        for (kₛ, vₛ) in ds.dataSets[1].series
            serieskey = parse.(Int, split(string(kₛ), ":")) .+ 1 #Transforms dimension into a 1-index Array
            serieskey = serieskey[mask[:series]] # Filters descriptive dimensions if alldims=false
            dimdesc = Vector{String}(undef, length(serieskey))
            if unpivotseries
                println("To be implemented")
            else
                for (i, n) in enumerate(serieskey)
                    dimdesc[i] = dims[:series][i][:values][n]["name"]
                end
                for (kₒ, vₒ) in vₛ["observations"]
                    obsdim = dims[:observation][1]["values"][parse(Int, string(kₒ))+1]["name"]
                    dimvalues = [dimdesc..., obsdim]
                    row = [dimvalues..., vₒ[1]]
                    push!(obs, (; zip(headers, row)...))
                end
            end
        end

    else
        for (k, v) in ds.dataSets[1].observations
            dimkey = parse.(Int, split(string(k), ":")) .+ 1 # Transforms dimension into a 1-index Array
            dimkey = dimkey[mask[:observation]] # Filters descriptive dimensions if alldims=false
            dimdesc = Vector{String}(undef, length(dimkey))
            for (i, n) in enumerate(dimkey)
                dimdesc[i] = dims[:observation][i]["values"][n]["name"]
            end
            value = v[1]
            push!(obs, (; zip(headers, [dimdesc..., value])...))  # Rows as Nametuples for Tables.jl
        end
    end

    # Collect info of all dimensions
    dims_all = OrderedDict()
    for v in values(ds.structure.dimensions)
        for i in v
            dims_all[i["name"]] = [j["name"] for j in i["values"]]
        end
    end

    # Collect info of all attributes

    attr_all = OrderedDict()
    for v in values(ds.structure.attributes)
        for i in v
            attr_all[i["name"]] = [j["name"] for j in i["values"]]
        end
    end

    # Create datatable
    Datatable(ds.structure.name,
        headers,
        obs,
        dims_all,
        attr_all) # Including headers, observation, all dimensions data, all attribute data
end

end #Module

max_stars_repo_name: lsaenzt/SDMX

max_issues_repo_name: lsaenzt/SDMX

max_forks_repo_name: lsaenzt/SDMX


text:
using Knet,AutoGrad,BenchmarkTools,Distributions
if !@isdefined(MODEL); MODEL=1; end

# Design choices:
# 1. time first vs layers fist
# 2. row major vs col major
# 3. merge inputs in t and single mmul Wx before lstm
# 4. do not merge inputs but cat(hidden,input) before mmul in lstm
# 5. do not merge and do Wx*input + Wh*hidden
# 6. combine gates in single mmul
# 7. merge outputs in t and single mmul Wy after each layer

#  mode=0  mode=1  mode=2   notes (min benchmark in ms with default args unless specified)
#  16.227  36.886  39.329   32b63d3 2017-03-25 32546 wps on aitest-gpu
# 249.359 505.947 552.944   32b63d3 2017-03-25 2315  wps on aitest-cpu
#  16.155  37.103  38.175   725b18b 2017-03-25 sum_outgrads uses axpy! 33529 best, 26743 sustained wps.
#   9.811  28.947  30.108   d3ea7d9 2017-03-26 41725 best, 35442 sustained wps. rnnlm with column major instances and split/merge inputs and outputs
#   3.269  11.644  12.728   vocab=100 MODEL=1 (column-major is worse with small vocab)
#   3.032  11.666  12.259   vocab=100 MODEL=2 (row-major is better with small vocab)
#  12.707  36.596  36.744   mean ms MODEL=1 vocab=10k seqlen=20; same as d3ea7d9 except mean instead of min benchmark
#  35.781  83.911 100.954   mean ms MODEL=1 vocab=10k seqlen=Gamma(4.32,4.88,mean=21.08)
#  34.245  91.380 111.394   mean ms MODEL=3 vocab=10k seqlen=Gamma(4.32,4.88,mean=21.08)
#  13.294  62.746  63.498   MODEL=5 seqlen=20,min
#  14.655  68.491  73.961   MODEL=5 seqlen=20,mean
#  15.109  71.534  73.894   MODEL=5 seqlen=Gamma(4.32,4.88),mean

# Usage:
#
# include("rnnlm.jl")
# m,s,o = main(iters=0)
# for i=1:2
# GC.gc(); @time main(model=m,state=s,optim=o,mode=0,iters=10)
# GC.gc(); @time main(model=m,state=s,optim=o,mode=1,iters=10)
# GC.gc(); @time main(model=m,state=s,optim=o,mode=2,iters=10)
# end
# GC.gc(); println(@benchmark main(model=$m,state=$s,optim=$o,mode=0))
# GC.gc(); println(@benchmark main(model=$m,state=$s,optim=$o,mode=1))
# GC.gc(); println(@benchmark main(model=$m,state=$s,optim=$o,mode=2))
# nothing

# Notes:
#
# sequence[t]: Vector{Int} token-minibatch input at time t
# TODO: consider learning the initial state instead of setting to 0
# pred3 = Array(pred3) #TODO: FIX BUGGY REDUCTION CODE FOR KNETARRAY IF PRED3 TOO BIG
# TODO: add gclip, gcheck

function main(;
              mode=0,
              iters=1,
              vocab=10000,
              batch=64,
              embed=128,
              hidden=[256],
              # seqlen=DiscreteUniform(20,20),
              seqlen=Gamma(4.32,4.88), # mean=4.32*4.88=21.08, var=4.32*4.88^2=102.88; mikolov ptb mean=21.10, var=102.92, excluding eos
              otype="Adam()",
              atype=(gpu()>=0 ? KnetArray{Float32} : Array{Float32}),
              model=initmodel(atype, hidden, vocab, embed),
              state=initstate(model,batch),
              optim=initoptim(model,otype),
              dropout=0,
              )
    if mode == 0
        for i in 1:iters
            sequence=randseq(vocab,batch,rand(seqlen))
            rnnlm(model, state, sequence; pdrop=dropout)
        end
    elseif mode == 1
        for i in 1:iters
            sequence=randseq(vocab,batch,rand(seqlen))
            rnnlmgrad(model, state, sequence; pdrop=dropout)
        end
    elseif mode == 2
        for i in 1:iters
            sequence=randseq(vocab,batch,rand(seqlen))
            grads = rnnlmgrad(model, state, sequence; pdrop=dropout) # 2:1830:3358  1:2585:3636
            update!(model, grads, optim)                             # 2:934 1:943
        end
    else
        error("mode=$mode")
    end
    return (model, state, optim)
end

# initoptim creates optimization parameters for each numeric weight
# array in the model.  This should work for a model consisting of any
# combination of tuple/array/dict.
initoptim(::KnetArray{<:Number},otype)=eval(Meta.parse(otype))
initoptim(::Array{<:Number},otype)=eval(Meta.parse(otype))
initoptim(a::AbstractDict,otype)=Dict(k=>initoptim(v,otype) for (k,v) in a)
initoptim(a,otype)=map(x->initoptim(x,otype), a)

# Create a random minibatch of sequences
function randseq(V,B,T)
    T = ceil(Int,T)
    s = Array{Vector{Int}}(undef,T)
    for t in 1:T
        s[t] = rand(1:V,B)
    end
    return s
end

if MODEL==1; @eval begin
@show MODEL

### Model 1: column-major, merge/split input/output
#  mode=0  mode=1  mode=2
#   9.811  28.947  30.108   vocab=10K
#   3.269  11.644  12.728   vocab=100 MODEL=1 (column-major is worse with small vocab)

function rnnlm(model, state, sequence, range=1:length(sequence)-1; pdrop=0) # 2:1830 1:2585
    index = vcat(sequence[range]...)
    input = Wm(model)[:,index]                          # 2:15
    for n = 1:nlayers(model)
        input = dropout(input, pdrop)
        input = Wx(model,n) * input                     # 2:26
        w,b,h,c = Wh(model,n),bh(model,n),hdd(state,n),cll(state,n)
        output = []
        j1 = j2 = 0
        for t in range
            j1 = j2 + 1
            j2 = j1 + length(sequence[t]) - 1
            input_t = input[:,j1:j2]                    # 2:35
            (h,c) = lstm(w,b,h,c,input_t)               # 2:991
            push!(output,h)
        end
        input = hcat(output...)                         # 2:39
    end
    pred1 = dropout(input,pdrop)
    pred2 = Wy(model) * pred1                           # 2:260  1:277:1132
    pred3 = pred2 .+ by(model)                          # 2:72  1:84:33
    nrows,ncols = size(pred3)
    golds = vcat(sequence[range .+ 1]...)
    golds += nrows*(0:(length(golds)-1))
    logp1 = logp(pred3,dims=1)                               # 2:354  1:1067:673
    logp2 = logp1[golds]
    logp3 = sum(logp2)
    return -logp3 / length(golds)
end

rnnlmgrad = grad(rnnlm)

function lstm(weight,bias,hidden,cell,input)            # 2:991  1:992:1617 (id:forw:back)
    gates   = weight * hidden .+ input .+ bias          # 2:312  1:434:499 (43+381+75) (cat+mmul+badd)
    h       = size(hidden,1)                            #
    forget  = sigm.(gates[1:h,:])                        # 2:134  1:98:99  (62+37) (index+sigm)
    ingate  = sigm.(gates[1+h:2h,:])                     # 2:99   1:73:123 (77+46)
    outgate = sigm.(gates[1+2h:3h,:])                    # 2:113  1:66:124 (87+37)
    change  = tanh.(gates[1+3h:4h,:])                    # 2:94   1:51:179 (130+49) replace end with 4h?
    cell    = cell .* forget + ingate .* change         # 2:137  1:106:202 (104+93+5) (bmul+bmul+add)
    hidden  = outgate .* tanh.(cell)                     # 2:100  1:69:194 (73+121) (tanh+bmul)
    return (hidden,cell)
end

nlayers(model)=div(length(model)-3,3)
Wm(model)=model[1]
Wx(model,n)=model[3n-1]
Wh(model,n)=model[3n]
bh(model,n)=model[3n+1]
Wy(model)=model[end-1]
by(model)=model[end]
hdd(state,n)=state[2n-1]
cll(state,n)=state[2n]

function initmodel(atype, hidden, vocab, embed)
    init(d...)=atype(xavier(Float32,d...))
    bias(d...)=atype(zeros(Float32,d...))
    N = length(hidden)
    model = Array{Any}(undef,3N+3)
    model[1] = init(embed,vocab) # Wm
    X = embed
    for n = 1:N
        H = hidden[n]
        model[3n-1] = init(4H,X) # Wx
        model[3n]   = init(4H,H) # Wh
        model[3n+1] = bias(4H,1) # bh
        model[3n+1][1:H] = 1     # forget gate bias = 1
        X = H
    end
    model[3N+2] = init(vocab,hidden[end]) # Wy
    model[3N+3] = bias(vocab,1)           # by
    return model
end

let blank = nothing; global initstate
function initstate(model, batch)
    N = nlayers(model)
    state = Array{Any}(undef,2N)
    for n = 1:N
        bias = bh(model,n)
        hidden = div(length(bias),4)
        if typeof(blank)!=typeof(bias) || size(blank)!=(hidden,batch)
            blank = fill!(similar(bias, hidden, batch),0)
        end
        state[2n-1] = state[2n] = blank
    end
    return state
end
end

end; elseif MODEL==2; @eval begin
@show MODEL

### Model 2: Row major, merge/split input/output
#  mode=0  mode=1  mode=2
#  16.501  36.352  37.112

function rnnlm(model, state, sequence, range=1:length(sequence)-1; pdrop=0) # 2:1830 1:2585
    index = vcat(sequence[range]...)
    input = Wm(model)[index,:]                          # 2:15
    for n = 1:nlayers(model)
        input = dropout(input, pdrop)
        input = input * Wx(model,n)                     # 2:26
        w,b,h,c = Wh(model,n),bh(model,n),hdd(state,n),cll(state,n)
        output = []
        j1 = j2 = 0
        for t in range
            j1 = j2 + 1
            j2 = j1 + length(sequence[t]) - 1
            input_t = input[j1:j2,:]                    # 2:35
            (h,c) = lstm(w,b,h,c,input_t)               # 2:991
            push!(output,h)
        end
        input = vcat(output...)                         # 2:39
    end
    pred1 = dropout(input,pdrop)
    pred2 = pred1 * Wy(model)                           # 2:260  1:277:1132
    pred3 = pred2 .+ by(model)                          # 2:72  1:84:33
    nrows,ncols = size(pred3)
    golds = vcat(sequence[range+1]...)
    golds = (golds-1)*nrows + (1:length(golds))
    logp1 = logp(pred3,dims=2)                               # 2:354  1:1067:673
    logp2 = logp1[golds]
    logp3 = sum(logp2)
    return -logp3 / length(golds)
end

rnnlmgrad = grad(rnnlm)

function lstm(weight,bias,hidden,cell,input)            # 2:991  1:992:1617 (id:forw:back)
    gates   = hidden * weight .+ input .+ bias          # 2:312  1:434:499 (43+381+75) (cat+mmul+badd)
    h       = size(hidden,2)                            #
    forget  = sigm.(gates[:,1:h])                        # 2:134  1:98:99  (62+37) (index+sigm)
    ingate  = sigm.(gates[:,1+h:2h])                     # 2:99   1:73:123 (77+46)
    outgate = sigm.(gates[:,1+2h:3h])                    # 2:113  1:66:124 (87+37)
    change  = tanh.(gates[:,1+3h:4h])                    # 2:94   1:51:179 (130+49) replace end with 4h?
    cell    = cell .* forget + ingate .* change         # 2:137  1:106:202 (104+93+5) (bmul+bmul+add)
    hidden  = outgate .* tanh.(cell)                     # 2:100  1:69:194 (73+121) (tanh+bmul)
    return (hidden,cell)
end

nlayers(model)=div(length(model)-3,3)
Wm(model)=model[1]
Wx(model,n)=model[3n-1]
Wh(model,n)=model[3n]
bh(model,n)=model[3n+1]
Wy(model)=model[end-1]
by(model)=model[end]
hdd(state,n)=state[2n-1]
cll(state,n)=state[2n]

function initmodel(atype, hidden, vocab, embed)
    init(d...)=atype(xavier(Float32,d...))
    bias(d...)=atype(zeros(Float32,d...))
    N = length(hidden)
    model = Array{Any}(undef,3N+3)
    model[1] = init(vocab,embed) # Wm
    X = embed
    for n = 1:N
        H = hidden[n]
        model[3n-1] = init(X,4H) # Wx
        model[3n]   = init(H,4H) # Wh
        model[3n+1] = bias(1,4H) # bh
        model[3n+1][1:H] = 1     # forget gate bias = 1
        X = H
    end
    model[3N+2] = init(hidden[end],vocab) # Wy
    model[3N+3] = bias(1,vocab)           # by
    return model
end

let blank = nothing; global initstate
function initstate(model, batch)
    N = nlayers(model)
    state = Array{Any}(undef,2N)
    for n = 1:N
        bias = bh(model,n)
        hidden = div(length(bias),4)
        if typeof(blank)!=typeof(bias) || size(blank)!=(batch,hidden)
            blank = fill!(similar(bias, batch,hidden),0)
        end
        state[2n-1] = state[2n] = blank
    end
    return state
end
end

end; elseif MODEL==3; @eval begin
@show MODEL

### Model 3: column-major, concat hidden/input, merge output
#  mode=0  mode=1  mode=2
#  10.458  32.010  33.764   vocab=10K
#   3.837  12.536  13.083   vocba=100

function rnnlm(model, state, sequence, range=1:length(sequence)-1; pdrop=0) # 2:1830 1:2585
    input = []
    for t in range
        push!(input, Wm(model)[:,sequence[t]])
    end
    for n = 1:nlayers(model)
        w,b,h,c = Wh(model,n),bh(model,n),hdd(state,n),cll(state,n)
        for i in 1:length(range)
            t = range[i]
            input_t = input[i]
            input_t = dropout(input_t, pdrop)
            (h,c) = lstm(w,b,h,c,input_t)               # 2:991
            input[i] = h
        end
    end
    pred0 = hcat(input...)
    pred1 = dropout(pred0,pdrop)
    pred2 = Wy(model) * pred1                           # 2:260  1:277:1132
    pred3 = pred2 .+ by(model)                          # 2:72  1:84:33
    nrows,ncols = size(pred3)
    golds = vcat(sequence[range+1]...)
    golds += nrows*(0:(length(golds)-1))
    logp1 = logp(pred3,dims=1)                               # 2:354  1:1067:673
    logp2 = logp1[golds]
    logp3 = sum(logp2)
    return -logp3 / length(golds)
end

rnnlmgrad = grad(rnnlm)

function lstm(weight,bias,hidden,cell,input)            # 2:991  1:992:1617 (id:forw:back)
    gates   = weight * vcat(hidden, input) .+ bias      # 2:312  1:434:499 (43+381+75) (cat+mmul+badd)
    h       = size(hidden,1)                            #
    forget  = sigm.(gates[1:h,:])                        # 2:134  1:98:99  (62+37) (index+sigm)
    ingate  = sigm.(gates[1+h:2h,:])                     # 2:99   1:73:123 (77+46)
    outgate = sigm.(gates[1+2h:3h,:])                    # 2:113  1:66:124 (87+37)
    change  = tanh.(gates[1+3h:4h,:])                    # 2:94   1:51:179 (130+49) replace end with 4h?
    cell    = cell .* forget + ingate .* change         # 2:137  1:106:202 (104+93+5) (bmul+bmul+add)
    hidden  = outgate .* tanh.(cell)                     # 2:100  1:69:194 (73+121) (tanh+bmul)
    return (hidden,cell)
end

nlayers(model)=div(length(model)-3,2)
Wm(model)=model[1]
Wx(model,n)=nothing
Wh(model,n)=model[2n]
bh(model,n)=model[2n+1]
Wy(model)=model[end-1]
by(model)=model[end]
hdd(state,n)=state[2n-1]
cll(state,n)=state[2n]

function initmodel(atype, hidden, vocab, embed)
    init(d...)=atype(xavier(Float32,d...))
    bias(d...)=atype(zeros(Float32,d...))
    N = length(hidden)
    model = Array{Any}(undef,2N+3)
    model[1] = init(embed,vocab) # Wm
    X = embed
    for n = 1:N
        H = hidden[n]
        model[2n]   = init(4H,H+X) # Wh
        model[2n+1] = bias(4H,1) # bh
        model[2n+1][1:H] = 1     # forget gate bias = 1
        X = H
    end
    model[2N+2] = init(vocab,hidden[end]) # Wy
    model[2N+3] = bias(vocab,1)           # by
    return model
end

let blank = nothing; global initstate
function initstate(model, batch)
    N = nlayers(model)
    state = Array{Any}(undef,2N)
    for n = 1:N
        bias = bh(model,n)
        hidden = div(length(bias),4)
        if typeof(blank)!=typeof(bias) || size(blank)!=(hidden,batch)
            blank = fill!(similar(bias, hidden, batch),0)
        end
        state[2n-1] = state[2n] = blank
    end
    return state
end
end

end; elseif MODEL==4; @eval begin
@show MODEL

### Model 4: row-major, concat hidden/input, merge output
#  mode=0  mode=1  mode=2
#  16.298  39.139  39.879  vocab=10K
#   3.173  11.715  12.223  vocab=100

function rnnlm(model, state, sequence, range=1:length(sequence)-1; pdrop=0) # 2:1830 1:2585
    input = []
    for t in range
        push!(input, Wm(model)[sequence[t],:])
    end
    for n = 1:nlayers(model)
        w,b,h,c = Wh(model,n),bh(model,n),hdd(state,n),cll(state,n)
        for i in 1:length(range)
            t = range[i]
            input_t = input[i]
            input_t = dropout(input_t, pdrop)
            (h,c) = lstm(w,b,h,c,input_t)               # 2:991
            input[i] = h
        end
    end
    pred0 = vcat(input...)
    pred1 = dropout(pred0,pdrop)
    pred2 = pred1 * Wy(model)                           # 2:260  1:277:1132
    pred3 = pred2 .+ by(model)                          # 2:72  1:84:33
    nrows,ncols = size(pred3)
    golds = vcat(sequence[range+1]...)
    golds = (golds-1)*nrows + (1:length(golds))
    logp1 = logp(pred3,dims=2)                               # 2:354  1:1067:673
    logp2 = logp1[golds]
    logp3 = sum(logp2)
    return -logp3 / length(golds)
end

rnnlmgrad = grad(rnnlm)

function lstm(weight,bias,hidden,cell,input)            # 2:991  1:992:1617 (id:forw:back)
    gates   = hcat(hidden, input) * weight .+ bias      # 2:312  1:434:499 (43+381+75) (cat+mmul+badd)
    h       = size(hidden,2)                            #
    forget  = sigm.(gates[:,1:h])                        # 2:134  1:98:99  (62+37) (index+sigm)
    ingate  = sigm.(gates[:,1+h:2h])                     # 2:99   1:73:123 (77+46)
    outgate = sigm.(gates[:,1+2h:3h])                    # 2:113  1:66:124 (87+37)
    change  = tanh.(gates[:,1+3h:4h])                    # 2:94   1:51:179 (130+49) replace end with 4h?
    cell    = cell .* forget + ingate .* change         # 2:137  1:106:202 (104+93+5) (bmul+bmul+add)
    hidden  = outgate .* tanh.(cell)                     # 2:100  1:69:194 (73+121) (tanh+bmul)
    return (hidden,cell)
end

nlayers(model)=div(length(model)-3,2)
Wm(model)=model[1]
Wx(model,n)=nothing
Wh(model,n)=model[2n]
bh(model,n)=model[2n+1]
Wy(model)=model[end-1]
by(model)=model[end]
hdd(state,n)=state[2n-1]
cll(state,n)=state[2n]

function initmodel(atype, hidden, vocab, embed)
    init(d...)=atype(xavier(Float32,d...))
    bias(d...)=atype(zeros(Float32,d...))
    N = length(hidden)
    model = Array{Any}(undef,2N+3)
    model[1] = init(vocab,embed) # Wm
    X = embed
    for n = 1:N
        H = hidden[n]
        model[2n]   = init(H+X,4H) # Wh
        model[2n+1] = bias(1,4H) # bh
        model[2n+1][1:H] = 1     # forget gate bias = 1
        X = H
    end
    model[2N+2] = init(hidden[end],vocab) # Wy
    model[2N+3] = bias(1,vocab)           # by
    return model
end

let blank = nothing; global initstate
function initstate(model, batch)
    N = nlayers(model)
    state = Array{Any}(undef,2N)
    for n = 1:N
        bias = bh(model,n)
        hidden = div(length(bias),4)
        if typeof(blank)!=typeof(bias) || size(blank)!=(batch,hidden)
            blank = fill!(similar(bias, batch, hidden),0)
        end
        state[2n-1] = state[2n] = blank
    end
    return state
end
end

end; elseif MODEL==5; @eval begin
@show MODEL

### Model 5: Column-major, no split/merge, concat hidden/input
#  mode=0  mode=1  mode=2
#  13.294  62.746  63.498  seqlen=20,min
#  14.655  68.491  73.961  seqlen=20,mean
#  15.109  71.534  76.200  seqlen=Gamma(4.32,4.88),mean

function rnnlm(model, state, sequence; pdrop=0)
    T = length(sequence)
    N = nlayers(model)
    state = copy(state)
    total = count = 0
    for t in 1:(T-1)
        input = Wm(model)[:,sequence[t]]
        for n in 1:N
            w,b,h,c = Wh(model,n),bh(model,n),hdd(state,n),cll(state,n)
            input = dropout(input,pdrop)
            (h,c) = lstm(w,b,h,c,input)
            input = h
            state[2n-1] = h
            state[2n] = c
        end
        input = dropout(input,pdrop)
        logp0 = Wy(model) * input .+ by(model)
        logp1 = logp(logp0,dims=1)
        golds = sequence[t+1]
        index = golds + size(logp1,1)*(0:(length(golds)-1))
        logp2 = logp1[index]
        total += sum(logp2)
        count += length(golds)
    end
    return -total/count
end

#=  ### alternative implementation
function rnnlm(model, state, sequence; pdrop=0)
    T = length(sequence)
    input = Array{Any}(undef,T-1)
    for t in 1:(T-1)
        input[t] = Wm(model)[:,sequence[t]]
    end
    for n = 1:nlayers(model)
        w,b,h,c = Wh(model,n),bh(model,n),hdd(state,n),cll(state,n)
        for t in 1:(T-1)
            input_t = input[t]
            input_t = dropout(input_t, pdrop)
            (h,c) = lstm(w,b,h,c,input_t)               # 2:991
            input[t] = h
        end
    end
    total = count = 0
    for t in 1:(T-1)
        pred0 = input[t]
        pred1 = dropout(pred0,pdrop)
        pred2 = Wy(model) * pred1                           # 2:260  1:277:1132
        pred3 = pred2 .+ by(model)                          # 2:72  1:84:33
        nrows,ncols = size(pred3)
        golds = sequence[t+1]
        index = golds + nrows*(0:(length(golds)-1))
        logp1 = logp(pred3,dims=1)                               # 2:354  1:1067:673
        logp2 = logp1[index]
        logp3 = sum(logp2)
        total += logp3
        count += length(golds)
    end
    return -total/count
end
=#

rnnlmgrad = grad(rnnlm)

function lstm(weight,bias,hidden,cell,input)            # 2:991  1:992:1617 (id:forw:back)
    gates   = weight * vcat(hidden, input) .+ bias      # 2:312  1:434:499 (43+381+75) (cat+mmul+badd)
    h       = size(hidden,1)                            #
    forget  = sigm.(gates[1:h,:])                        # 2:134  1:98:99  (62+37) (index+sigm)
    ingate  = sigm.(gates[1+h:2h,:])                     # 2:99   1:73:123 (77+46)
    outgate = sigm.(gates[1+2h:3h,:])                    # 2:113  1:66:124 (87+37)
    change  = tanh.(gates[1+3h:4h,:])                    # 2:94   1:51:179 (130+49) replace end with 4h?
    cell    = cell .* forget + ingate .* change         # 2:137  1:106:202 (104+93+5) (bmul+bmul+add)
    hidden  = outgate .* tanh.(cell)                     # 2:100  1:69:194 (73+121) (tanh+bmul)
    return (hidden,cell)
end

nlayers(model)=div(length(model)-3,2)
Wm(model)=model[1]
Wx(model,n)=nothing
Wh(model,n)=model[2n]
bh(model,n)=model[2n+1]
Wy(model)=model[end-1]
by(model)=model[end]
hdd(state,n)=state[2n-1]
cll(state,n)=state[2n]

function initmodel(atype, hidden, vocab, embed)
    init(d...)=atype(xavier(Float32,d...))
    bias(d...)=atype(zeros(Float32,d...))
    N = length(hidden)
    model = Array{Any}(undef,2N+3)
    model[1] = init(embed,vocab) # Wm
    X = embed
    for n = 1:N
        H = hidden[n]
        model[2n]   = init(4H,H+X) # Wh
        model[2n+1] = bias(4H,1) # bh
        model[2n+1][1:H] = 1     # forget gate bias = 1
        X = H
    end
    model[2N+2] = init(vocab,hidden[end]) # Wy
    model[2N+3] = bias(vocab,1)           # by
    return model
end

let blank = nothing; global initstate
function initstate(model, batch)
    N = nlayers(model)
    state = Array{Any}(undef,2N)
    for n = 1:N
        bias = bh(model,n)
        hidden = div(length(bias),4)
        if typeof(blank)!=typeof(bias) || size(blank)!=(hidden,batch)
            blank = fill!(similar(bias, hidden, batch),0)
        end
        state[2n-1] = state[2n] = blank
    end
    return state
end
end

### Model X: Column major, no split/merge, no concat hidden/input

### Model X: Column major, no split/merge, no concat hidden/input, no merging gates

end; end # if MODEL==1

nothing

### MODEL=2: row-major, merge output, concat hidden/input

# Forward pass profile (2585/7166 of total)
# 2585 ./<missing>:0; (::#kw##rnnlm)(::Array{Any,1}, ::#rnnlm, ::AutoGrad.R...
#  86   .../dyuret/knet/master/prof/rnnlm.jl:158; input = embed[sequence[t],:]
#  999  .../dyuret/knet/master/prof/rnnlm.jl:159; pred,state = predict(model,state,input; pdrop=pdrop)
#   992 .../dyuret/knet/master/prof/rnnlm.jl:147; (newstate[2k-1],newstate[2k]) = lstm(model[2k-1],model[2k],state[2k-1],state[2k],input)
#    505 .../dyuret/knet/master/prof/rnnlm.jl:129; gates   = hcat(input,hidden) * weight .+ bias
#    98  .../dyuret/knet/master/prof/rnnlm.jl:131; forget  = sigm(gates[:,1:hsize])
#    73  .../dyuret/knet/master/prof/rnnlm.jl:132; ingate  = sigm(gates[:,1+hsize:2hsize])
#    66  .../dyuret/knet/master/prof/rnnlm.jl:133; outgate = sigm(gates[:,1+2hsize:3hsize])
#    51  .../dyuret/knet/master/prof/rnnlm.jl:134; change  = tanh(gates[:,1+3hsize:end])
#    106 .../dyuret/knet/master/prof/rnnlm.jl:135; cell    = cell .* forget + ingate .* change
#    69  .../dyuret/knet/master/prof/rnnlm.jl:136; hidden  = outgate .* tanh(cell)
#  51   .../dyuret/knet/master/prof/rnnlm.jl:165; pred0 = vcat(preds...)
#  277  .../dyuret/knet/master/prof/rnnlm.jl:167; pred2 = pred1 * model[end-1]
#  84   .../dyuret/knet/master/prof/rnnlm.jl:168; pred3 = pred2 .+ model[end]
#  1067 .../dyuret/knet/master/prof/rnnlm.jl:169; logp1 = logp(pred3,2)
#   424 ...ret/.julia/v0.5/Knet/src/unary.jl:176; x1 = maximum(x,d...)
#   53  ...ret/.julia/v0.5/Knet/src/unary.jl:177; x2 = x .- x1
#   88  ...ret/.julia/v0.5/Knet/src/unary.jl:178; x3 = exp(x2)
#   414 ...ret/.julia/v0.5/Knet/src/unary.jl:179; x4 = sum(x3,d...)
#   1   ...ret/.julia/v0.5/Knet/src/unary.jl:180; x5 = log(x4)
#   85  ...ret/.julia/v0.5/Knet/src/unary.jl:181; x6 = x2 .- x5
#  17   .../dyuret/knet/master/prof/rnnlm.jl:174; index[i] = i + (golds[i]-1)*nrows

# Backward pass profile (3636/7166 of total)
# 3636 ...et/.julia/v0.5/AutoGrad/src/core.jl:40; (::AutoGrad.##gradfun#1#3{#rnnlm,Int64})(::Array{Any,1...
#  2948 core.jl:231 og = r.func(Grad{i},n.outgrad,r.value,r.args...;r.kwargs...)
#   1512 ./<missing>:0; *(::Type{AutoGrad.Grad{2}}, ::Knet.KnetArray{Float32...
#   679  ./<missing>:0; logp(::Type{AutoGrad.Grad{1}}, ::Knet.KnetArray{Floa...
#   286  ./<missing>:0; .*(::Type{AutoGrad.Grad{1}}, ::Knet.KnetArray{Float3...
#   109  ./<missing>:0; sigm(::Type{AutoGrad.Grad{1}}, ::Knet.KnetArray{Floa...
#   95   ./<missing>:0; tanh(::Type{AutoGrad.Grad{1}}, ::Knet.KnetArray{Floa...
#   92   ./<missing>:0; .+(::Type{AutoGrad.Grad{2}}, ::Knet.KnetArray{Float3...
#   60   ...AutoGrad/src/base/abstractarray.jl:85; cat(::Type{AutoGrad.Grad{2}}, ::Knet.KnetArray{Float...
#   13   ./<missing>:0; getindex(::Type{AutoGrad.Grad{1}}, ::Knet.KnetArray{...
#   13   ./essentials.jl:216; vector_any()
#  684  ...et/.julia/v0.5/AutoGrad/src/core.jl:233; backward_pass(::AutoGrad.Rec{Array{Any,1}}, ::AutoGra...
#   271 ...uret/.julia/v0.5/Knet/src/karray.jl:995; sum_outgrads(::Knet.KnetArray{Float32,2}, ::AutoGrad....
#   168 ...lia/v0.5/AutoGrad/src/interfaces.jl:71; sum_outgrads(::Void, ::AutoGrad.UngetIndex)
#   167 ...lia/v0.5/AutoGrad/src/interfaces.jl:92; sum_outgrads(::Array{Any,1}, ::AutoGrad.UngetIndex)
#   72  ...uret/.julia/v0.5/Knet/src/karray.jl:992; sum_outgrads(::Knet.KnetArray{Float32,2}, ::Knet.Knet...

# Update profile (943/7166 of total)
# 940 ...yuret/.julia/v0.5/Knet/src/update.jl:404; update!(wi,gi,pi)
#  31  ...uret/.julia/v0.5/Knet/src/update.jl:346; scale!(p.beta1, p.fstm)
#  53  ...uret/.julia/v0.5/Knet/src/update.jl:347; axpy!(1-p.beta1, g, p.fstm)
#  35  ...uret/.julia/v0.5/Knet/src/update.jl:348; scale!(p.beta2, p.scndm)
#  82  ...uret/.julia/v0.5/Knet/src/update.jl:349; axpy!(1-p.beta2, g .* g, p.scndm)
#  29  ...uret/.julia/v0.5/Knet/src/update.jl:350; fstm_corrected = p.fstm / (1 - p.beta1 ^ p.t)
#  27  ...uret/.julia/v0.5/Knet/src/update.jl:351; scndm_corrected = p.scndm / (1 - p.beta2 ^ p.t)
#  683 ...uret/.julia/v0.5/Knet/src/update.jl:352; axpy!(-p.lr, (fstm_corrected ./ (sqrt(scndm_corrected) + p.eps)), w)

# using Knet
# include(Knet.dir("prof/rnnlm.jl"))
# (m,s,o) = main(iters=0)
# for i in 1:10
#     @time main(model=m,optim=o,state=s,mode=0,iters=10)
#     @time main(model=m,optim=o,state=s,mode=1,iters=10)
#     @time main(model=m,optim=o,state=s,mode=2,iters=10)
#     println()
# end

# commit 359d3646 2018-09-22, julia 1.0.0   vs commit 4aa5f92f 2018-08-14, julia 0.6.4
# GPU:V100, CPU:Intel(R) Xeon(R) CPU E5-2680 v4 @ 2.40GHz
# MODEL = 1
#
# 0.043899 seconds (79.10 k allocations: 2.382 MiB)			0.042887 seconds (50.72 k allocations: 1.868 MiB)
# 0.256060 seconds (280.02 k allocations: 8.785 MiB, 35.70% gc time)    0.219261 seconds (160.46 k allocations: 7.654 MiB)
# 0.294070 seconds (335.07 k allocations: 11.161 MiB, 28.79% gc time)   0.308199 seconds (197.60 k allocations: 7.673 MiB, 55.45% gc time)

ext: jl


text:
# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.

using Test

@testset "UnitCommitment" begin
    include("instance_test.jl")
    include("model_test.jl")
    include("sensitivity_test.jl")
    include("screening_test.jl")
    include("convert_test.jl")
    include("validate_test.jl")
    include("initcond_test.jl")
end


text:
#=
predict or center stepper
=#

mutable struct PredOrCentStepper{T <: Real} <: Stepper{T}
    use_adjustment::Bool
    use_curve_search::Bool
    max_cent_steps::Int
    pred_prox_bound::T
    searcher_options

    prev_alpha::T
    cent_count::Int
    rhs::Point{T}
    dir::Point{T}
    temp::Point{T}
    dir_noadj::Point{T}
    dir_adj::Point{T}
    dir_temp::Vector{T}

    searcher::StepSearcher{T}
    unadj_only::Bool
    unadj_alpha::T

    function PredOrCentStepper{T}(;
        use_adjustment::Bool = true,
        use_curve_search::Bool = use_adjustment,
        max_cent_steps::Int = 4,
        pred_prox_bound::T = T(0.0332), # from Alfonso solver
        searcher_options...
        ) where {T <: Real}
        stepper = new{T}()
        if use_curve_search
            # can only use curve search if using adjustment
            @assert use_adjustment
        end
        stepper.use_adjustment = use_adjustment
        stepper.use_curve_search = use_curve_search
        stepper.max_cent_steps = max_cent_steps
        stepper.pred_prox_bound = pred_prox_bound
        stepper.searcher_options = searcher_options
        return stepper
    end
end

function load(stepper::PredOrCentStepper{T}, solver::Solver{T}) where {T <: Real}
    model = solver.model
    if stepper.use_adjustment && !any(Cones.use_dder3, model.cones)
        # model has no cones that use third order deriv
        stepper.use_adjustment = false
    end

    stepper.prev_alpha = one(T)
    stepper.cent_count = 0
    stepper.rhs = Point(model)
    stepper.dir = Point(model)
    stepper.temp = Point(model)
    stepper.dir_noadj = Point(model, ztsk_only = true)
    if stepper.use_adjustment
        stepper.dir_adj = Point(model, ztsk_only = true)
    end
    stepper.dir_temp = zeros(T, length(stepper.rhs.vec))

    stepper.searcher = StepSearcher{T}(model; stepper.searcher_options...)
    stepper.unadj_only = false
    stepper.unadj_alpha = 0

    return stepper
end

function step(stepper::PredOrCentStepper{T}, solver::Solver{T}) where {T <: Real}
    point = solver.point
    model = solver.model
    rhs = stepper.rhs
    dir = stepper.dir
    dir_noadj = stepper.dir_noadj
    cones = model.cones

    # update linear system solver factorization
    solver.time_upsys += @elapsed update_lhs(solver.syssolver, solver)

    # decide whether to predict or center
    is_pred = (stepper.cent_count >= stepper.max_cent_steps) ||
        (stepper.searcher.prox < stepper.pred_prox_bound)
    stepper.cent_count = (is_pred ? 0 : stepper.cent_count + 1)
    rhs_fun_noadj = (is_pred ? update_rhs_pred : update_rhs_cent)

    # get unadjusted direction
    solver.time_uprhs += @elapsed rhs_fun_noadj(solver, rhs)
    solver.time_getdir += @elapsed get_directions(stepper, solver)
    copyto!(dir_noadj.vec, dir.vec)
    try_noadj = true

    if stepper.use_adjustment
        # get adjustment direction
        rhs_fun_adj = (is_pred ? update_rhs_predadj : update_rhs_centadj)
        dir_adj = stepper.dir_adj
        solver.time_uprhs += @elapsed rhs_fun_adj(solver, rhs, dir)
        solver.time_getdir += @elapsed get_directions(stepper, solver)
        copyto!(dir_adj.vec, dir.vec)

        if stepper.use_curve_search
            # do single curve search with adjustment
            stepper.unadj_only = false
            solver.time_search += @elapsed alpha = search_alpha(point, model,
                stepper)

            if iszero(alpha)
                # try not using adjustment
                @warn("very small alpha in curve search; trying without adjustment")
            else
                # step
                update_stepper_points(alpha, point, stepper, false)
                stepper.prev_alpha = alpha
                return true
            end
        else
            # do two line searches, first for unadjusted alpha, then for corrected alpha
            try_noadj = false
            stepper.unadj_only = true
            solver.time_search += @elapsed alpha = search_alpha(point, model,
                stepper)
            stepper.unadj_alpha = alpha
            unadj_sched = stepper.searcher.prev_sched

            if !iszero(alpha)
                stepper.unadj_only = false
                solver.time_search += @elapsed alpha = search_alpha(point,
                    model, stepper)

                if iszero(alpha)
                    # use unadjusted direction: start at alpha found
                    # during unadjusted direction search
                    stepper.unadj_only = true
                    solver.time_search += @elapsed alpha = search_alpha(point,
                        model, stepper, sched = unadj_sched)
                    # check alpha didn't decrease more
                    @assert stepper.searcher.prev_sched == unadj_sched
                end

                update_stepper_points(alpha, point, stepper, false)
                stepper.prev_alpha = alpha
                return true
            end
        end
    end

    if try_noadj
        # do line search in unadjusted direction
        stepper.unadj_only = true
        solver.time_search += @elapsed alpha = search_alpha(point, model, stepper)
    end

    if iszero(alpha)
        @warn("very small alpha in line search; terminating")
        solver.status = NumericalFailure
        stepper.prev_alpha = alpha
        return false
    end

    # step
    update_stepper_points(alpha, point, stepper, false)
    stepper.prev_alpha = alpha

    return true
end

expect_improvement(stepper::PredOrCentStepper) = iszero(stepper.cent_count)

function update_stepper_points(
    alpha::T,
    point::Point{T},
    stepper::PredOrCentStepper{T},
    ztsk_only::Bool,
    ) where {T <: Real}
    if ztsk_only
        cand = stepper.temp.ztsk
        copyto!(cand, point.ztsk)
        dir_noadj = stepper.dir_noadj.ztsk
    else
        cand = point.vec
        dir_noadj = stepper.dir_noadj.vec
    end

    @. cand += alpha * dir_noadj
    if !stepper.unadj_only
        dir_adj = (ztsk_only ? stepper.dir_adj.ztsk : stepper.dir_adj.vec)
        adj_factor = (stepper.use_curve_search ? abs2(alpha) :
            alpha * stepper.unadj_alpha)
        @. cand += adj_factor * dir_adj
    end

    return
end

print_header_more(stepper::PredOrCentStepper, solver::Solver) =
    @printf("%5s %9s", "step", "alpha")

function print_iteration_more(stepper::PredOrCentStepper, solver::Solver)
    step = (iszero(stepper.cent_count) ? "pred" : "cent")
    @printf("%5s %9.2e", step, stepper.prev_alpha)
    return
end


text:
splitstring(str::AbstractString, at::Union{Char,AbstractString}) = String.(split(str, at))

function Base.Unicode.isuppercase(s::AbstractString)
    result = true
    for ch in s
      if !isuppercase(ch)
          result = false
          break
      end
    end
    return result
end

function uppercasefirst(x::AbstractString)
    s = strip(x)
    return if isempty(s) || isuppercase(x)
               s
           else
               string(uppercase(x[1]), x[2:end])
           end
end


text:
#   Copyright (c) 2020-2022 Velexi Corporation
#
#   Licensed under the Apache License, Version 2.0 (the "License");
#   you may not use this file except in compliance with the License.
#   You may obtain a copy of the License at
#
#       http://www.apache.org/licenses/LICENSE-2.0
#
#   Unless required by applicable law or agreed to in writing, software
#   distributed under the License is distributed on an "AS IS" BASIS,
#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#   See the License for the specific language governing permissions and
#   limitations under the License.

"""
Blade_operators_mgs.jl defines versions of operations on the Blade type based
on modified Gram-Schmidt orthogonalization (as opposed to Householder
triangularization).

Notes
-----
* Compared with algorithms based on Householder triangularization, algorithms
  based on modified Gram-Schmidt orthogonalization typically have higher
  computational performance at the cost of lower accuracy and higher memory
  usage.
"""

# --- Exports

export wedge_mgs
export dual_mgs

# --- Method definitions

using LinearAlgebra: norm

"""
    wedge_mgs(B::Blade, C::Blade)::AbstractBlade

Return the outer product of `B` and `C`.

Implementation
==============
* The extension of the basis of `B` to a basis for span(B, C) is computed
  using modified Gram-Schmidt orthogonalization.
"""
function wedge_mgs(B::Blade, C::Blade)
    # Compute the rejections of the basis of `C` from the subspace represented
    # by `B`
    rejections = rejection(basis(C), B)

    # Orthogonalize the rejections using the modified Gram-Schmidt algorithm
    for j in 1:size(rejections, 2)
        norm = LinearAlgebra.norm(rejections[:, j])
        if norm > 0
            rejections[:, j] /= norm
        else
            return zero(B)
        end

        column = rejections[:, j]
        for j_inner in (j + 1):size(rejections, 2)
            rejections[:, j_inner] -=
                (rejections[:, j_inner] ⋅ column) * column
        end
    end

    # Construct the Blade representing the outer product
    B_wedge_C_basis = hcat(basis(B), rejections)
    Blade{typeof(B.volume)}(dim(B), grade(B) + grade(C),
                            B_wedge_C_basis, volume(B) * volume(C),
                            copy_basis=false)

#    rejections = rejection(basis(C), B)
#    Blade(hcat(basis(B), rejections), volume=volume(B) * volume(C))
end

"""
    dual_mgs(B::AbstractBlade, C::AbstractBlade)

Return the dual `B` relative to `C`.

Notes
=====

* `dual(B, C)` is only defined if (1) `B` and `C` are extended from real
  vector spaces of the same dimension and (2) the subspace represented by `B`
  is contained in subspace represented by `C`.

* The volume of `C` is ignored.

Implementation
==============

* The extension of `basis(B)` to `basis(C)` is computed using modified
  Gram-Schmidt orthogonalization.
"""
function dual_mgs(B::Blade, C::Blade)
    # --- Handle edge cases

    # Check that B and C are extended from the real vector spaces of the same
    # dimension
    if dim(B) != dim(C)
        throw(DimensionMismatch("`dim(B)` not equal to `dim(C)`"))
    end

    # Check that B is contained in C
    projection_coefficients = transpose(basis(C)) * basis(B)
    if LinearAlgebra.norm(projection_coefficients)^2 ≉ grade(B)
        throw(DimensionMismatch("`B` not contained in `C`"))
    end

    # Subspaces represented by B and C are the same
    if grade(B) == grade(C)
        return Scalar(volume(B))
    end

    # --- Compute dual using the basis vectors of `C` with the smallest
    #     projections (i.e., largest rejections) onto the basis of `B`.

    permutation = sortperm(sum(abs.(projection_coefficients), dims=2)[:, 1])
    dual_basis = rejection(basis(C)[:, permutation[1:grade(C) - grade(B)]], B)
    Blade(dual_basis, volume=volume(B))
end


text:
"""
```
Module      : FinanceLib
Description : Implement Base modules for the FinanceLib library
Copyright   : (c) 2022 Kishaloy Neogi
License     : MIT
Maintainer  : Kishaloy Neogi
Email       : nkishaloy@yahoo.com
```
The module describes the base modules of Finance like npv,xnpv,irr,xirr,time value of money etc.

PV is mentioned as PV,Future value as FV and Terminal value as TV

You may see the github repository at <https://github.com/n-kishaloy/FinanceLib.jl>
"""
module FinanceLib

import Dates

"""
`yearFrac(d0,d1) = Day difference (d1 - d0) in fraction of a year`

* d0 = reference date
* d1 = target date

This function is similar to its counterpart in MS Excel except it divides the day
by 365.25 instead of more complicated rules followed in Excel.
"""
yearFrac(d0,d1) = Dates.value(d1 - d0)/365.25

"""
`invYearFrac(t0, yF) = Reverse of yearFrac with anchor Date and YearFrac to calculate final Date.`

* t0 = Date from which to calculate.
* yF = Year Fraction as Float64.
"""
invYearFrac(t0, yF) = t0 + Dates.Day(round(yF * 365.25))

"""
PeriodSeries = Series of Tuples with 1st value gives period value in Float64
"""
const PeriodSeries   = Vector{Tuple{Float64, Float64}}

"""
DateSeries = Series of Tuples with 1st value gives period value in Date
"""
const DateSeries  = Vector{Tuple{Dates.Date, Float64}}

dateToPeriodSeries(d0, pd) = (((x,y),) -> (yearFrac(d0, x), y)).(pd)

"""
`disFactAnnual(r) = Discount factor for 1 period = 1/(1+r)`

* r = rate for 1 period
"""
disFactAnnual(r) = 1/(1+r)

"""
`disFact(r, n) = Discount factor = 1/(1+r)^n`

* r = rate for 1 period
* n = period given as Float64
"""
disFact(r, n) = 1/(1+r)^n

"""
`xdisFact(r, d0, d1) = Discount Factor between period = 1/(1+r)^yearFrac(d0, d1)`

* r  = rate for 1 year during period (d0,d1)
* d0 = begin date
* d1 = end date
"""
xdisFact(r, d0, d1) = 1/(1+r)^yearFrac(d0, d1)

"""
`fwdDisFact((r0,t0), (r1,t1))
    = Forward rate between t0 and t1 given as Float64
    = disFact(r1,t1) / disFact(r0,t0)`

* (r0, t0) = Tuple of Rate and Time in (0,t0)
* (r1, t1) = Tuple of Rate and Time in (0,t1)
"""
fwdDisFact((r0,t0), (r1,t1)) = disFact(r1,t1) / disFact(r0,t0)

"""
`tMul(r,n) = Time Multiplier`

* r = rate of increase / time period
* n = nos of time periods
"""
tMul(r,n) = (1+r)^n

"""
`tMul(r,n,m) = Time Multiplier`

*r = rate of increase / time period
*n = nos of time periods
*m = nos of period / time period
"""
tMul(r,n,m) = (1+r/m)^(n*m)

"""
`rateGwth(f,p,n) = Rate which gives a growth over a specified period`

* f = FV
* p = PV
* n = nos of time periods
"""
rateGwth(f,p,n) = (f/p)^(1/n) - 1

"""
`periodGwth(f,p,r) = Number of time periods to increase from p to f`

* f = FV
* p = PV
* r = rate of return
"""
periodGwth(f,p,r) = log(f/p)/log(1+r)

"""
`pv(fv,r,n) = PV of a Future cash flow`

*fv = Future cash flow
*r  = rate of return
*n  = number of periods
"""
pv(fv,r,n) = fv/(1+r)^n

"""
`pv(fv,r,n,m) = PV of a Future cash flow with multiple compounding per period`

* fv = Future cash flow
* r  = rate of return
* n  = number of periods
* m  = number of compounding per period
"""
pv(fv,r,n,m) = pv(fv,r/m,n*m)

"""
`pvr(fv,r,n) = PV of continuous growth`

* fv = FV
* r  = rate of return in year on year term
* n  = number of periods
"""
pvr(fv,r,n) = fv/r^n

"""
`pvc(fv,r,n) = PV of continuous expontial growth`

* fv = FV
* r  = rate of return in exponential term
* n  = number of periods
"""
pvc(fv,r,n) = fv/exp(r*n)

"""
`pvAnnuity(pmt,r,n,m) = PV of an annuity with multiple payments per period`

* pmt = payment made in each transaction
* r   = rate of return
* n   = number of periods (say, years)
* m   = number of payments per period (say, monthly where `m = 12`)
"""
pvAnnuity(pmt,r,n,m) = pmt/(r/m)*(1 - 1/(1 + r/m)^(n*m))

"""
`pvAnnuity(pmt,r,n) = PV of an annuity with single payments per period`

* pmt = payment made in each transaction
* r   = rate of return
* n   = number of periods (say, years)
"""
pvAnnuity(pmt,r,n) = pmt/r*(1 - 1/(1 + r)^n)

"""
`pvAnnuity(pmt,r) = PV of infinite payments`

* pmt = payments made at a time
* r   = rate of return between /payments/
"""
pvAnnuity(pmt,r) = pmt/r

"""
`fvAnnuity(pmt,r,n,m) = FV of an annuity with multiple payments per period`

* pmt = payment made in each transaction
* r   = rate of return
* n   = number of periods (say, years)
* m   = number of payments per period (say, monthly where *m = 12*)
"""
fvAnnuity(pmt,r,n,m) = pmt/(r/m) * ((1 + r/m)^(n*m) - 1)

"""
`fvAnnuity(pmt,r,n) = FV of an annuity with single payments per period`

* pmt = payment made in each transaction
* r   = rate of return
* n   = number of periods (say, years)
"""
fvAnnuity(pmt,r,n) = pmt/r * ((1 + r)^n - 1)

"""
`pmt(pv,r,n,m) = Payment to cover the PV of an Annuity`

* pv = PV of Annuity
* r  = rate of return
* n  = number of periods (say, years)
* m  = number of payments per period (say, monthly where m = 12)
"""
pmt(pv,r,n,m) = pv * (r/m) / (1 - 1/(1 + r/m)^(n*m))

"""
`pmt(pv,r,n,m) = Payment to cover the PV of an Annuity`

* pv = PV of Annuity
* r  = rate of return
* n  = number of periods (say, years)
"""
pmt(pv,r,n) = pv * r / (1 - 1/(1 + r)^n)

"""
`fmt(fv,r,n,m) = Payment to cover the FV of an Annuity`

* fv = FV of Annuity
* n  = number of periods (say, years)
* m  = number of payments per period (say, monthly where m = 12)
* r  = rate of return
"""
fmt(fv,r,n,m) = fv * (r/m) / ((1 + r/m)^(n*m) - 1)

"""
`fmt(fv,r,n) = Payment to cover the FV of an Annuity`

* fv = FV of Annuity
* n  = number of periods (say, years)
* r  = rate of return
"""
fmt(fv,r,n) = fv * r / ((1 + r)^n - 1)

"""
`fv(pv,r,n,m) = FV of a Future cash flow with multiple compounding per period`

* pv = Present cash flow
* r  = rate of return
* n  = number of periods
* m  = number of compounding per period
"""
fv(pv,r,n,m) = pv*(1+r/m)^(n*m)

"""
`fv(pv,r,n) = FV of a Future cash flow with multiple compounding per period`

* pv = Present cash flow
* r  = rate of return
* n  = number of periods
"""
fv(pv,r,n) = pv*(1+r)^n

"""
`fvc(pv,r,n) = FV of continuous expontial growth`

* pv = PV
* r = rate of return
* n = number of periods
"""
fvc(pv,r,n) = pv*exp(r*n)

"""
`effR(r,m) = effective rate of return for multiple compounding per period`

  * r = nominal rate of return in a period
  * m = number of compounding per period
"""
effR(r,m) = (1.0 + r/m)^m - 1.0

"""
`expToEffR(r) = Exp rate to effective rate`

* r = exponential rate
"""
expToEffR(r) = exp(r) - 1.0

"""
`expToNomR(r,m) = Exp rate to nominal rate`

* r = exponential rate
* m = number of compounding per period
"""
expToNomR(r,m) = (exp(r/m) - 1.0)*m

"""
`nominalRate(r,m) = nominal rate of return for multiple compounding per period`

  * r = effective rate of return in a period
  * m = number of compounding per period
"""
nomR(r,m) = ((1.0 + r)^(1.0/m) - 1.0)*m

"""
`expR(r) = real rate of return for continuous exponential compounding`

  * r = effective rate of return in a period
"""
expR(r) = log(1.0 + r)

"""
`expR(r,m) = real rate of return for continuous exponential compounding`

  * m = number of compounding per period
  * r = nominal rate of return in a period
"""
expR(r,m) = m*log(1.0 + r/m)

"""
`npv(r,tim,cf,t0) = NPV of cash flows against time given in periods`

* r   = rate of return across the periods
* tim = vector of time of cash flows given as Float64
* cf  = vector of corresponding cash flows
* t0  = time period at which the NPV is sought. Essentially, NPV(ti - t0)
"""
function npv(r,tim,cf,t0) r += 1.0; sum(cf ./ r .^ tim) * r^t0 end

"""
`npv(r,cf,t0) = NPV of cash flows against time given in periods`

* r   = rate of return across the periods
* cf  = vector of tuple of (time in Float64, cash flows)
* t0  = time period at which the NPV is sought. Essentially, NPV(ti - t0)
"""
function npv(r,cf,t0) r += 1.0; sum( (((t,x),) -> x/r^t).(cf) ) * r^t0 end

"""
`npv(tr::PeriodSeries,cf,t0) = NPV of cash flows with rate given by period`

* tr  = vector of Tuple (time in Float64, rate between periods @ rate / period)
* cf  = vector of corresponding cash flows
* t0  = time period at which the NPV is sought

Note: this function assumes that the time is given in ascending order, thus between 2 consecutive period given in the vector tr, tr[i-1] and tr[i], the rate / period is assumed to be r[i] between period t[i-1] to t[i], essentially f(i-1,i) = tr[i][2] ==> Forward rate per period.

Also the time t0 is assumed to be **less than** 1st time period => t0 < tr[1][1]
"""
npv(tr::PeriodSeries,cf,t0) = (foldl(((y, i)-> (y+cf[i])/(1 + tr[i][2])^(tr[i][1] - tr[i-1][1])), length(tr):-1:2; init=0.0) + cf[1])/(1 + tr[1][2])^(tr[1][1] - t0)

import Roots

"""
`irr(tim,cf) = IRR of cash flow against time given in periods`

* tim = vector of time of cash flows given as Float64
* cf  = vector of corresponding cash flows
"""
irr(tim,cf) = Roots.find_zero(r -> npv(r,tim,cf,0.0),(0.0,0.5))

"""
`irr(cf) = IRR of cash flow against time given in periods`

* cf  = vector of tuple of (time in Float64, cash flows)
"""
irr(cf) = Roots.find_zero(r -> npv(r,cf,0.0),(0.0,0.5))

"""
`xnpv(r,tim,cf,t0) = NPV of cash flows against time given by Date`

* r   = rate of return across the years
* tim = vector of time of cash flows given as Date
* cf  = vector of corresponding cash flows
* t0  = Date at which the NPV is sought.
"""
function xnpv(r,tim,cf,t0) r += 1.0; sum(cf./(t -> r^yearFrac(t0,t)).(tim)) end

"""
`xnpv(r,cf,t0) = NPV of cash flows against time given by Date`

* r   = rate of return across the years
* cf  = vector of tuple of (time in Date, cash flows)
* t0  = Date at which the NPV is sought.
"""
function xnpv(r,cf,t0) r += 1.0; sum((((t,x),) -> x/r^yearFrac(t0,t)).(cf)) end

"""
`xnpv(tr::DateSeries,cf,t0) = NPV of cash flows with rate given by Date`

* tr  = vector of Tuple (time in Date, rate between Datess @ rate / year)
* cf  = vector of corresponding cash flows
* t0  = Date at which the NPV is sought

Note: this function assumes that the time is given in ascending order, thus between 2 consecutive period given in the vector tr, tr[i-1] and tr[i], the rate / period is assumed to be r[i] between period t[i-1] to t[i], essentially f(i-1,i) = tr[i][2] ==> Forward rate per period.

Also the time t0 is assumed to be **less than** 1st time period => t0 < tr[1][1]
"""
xnpv(tr::DateSeries,cf,t0) = (foldl(((y, i)-> (y+cf[i])/(1 + tr[i][2])^yearFrac(tr[i-1][1], tr[i][1])), length(tr):-1:2; init=0.0) + cf[1])/(1 + tr[1][2])^yearFrac(t0, tr[1][1])

"""
`xirr(tim,cf) = IRR of cash flow against time given in Dates`

* tim = vector of time of cash flows given as Dates
* cf  = vector of corresponding cash flows
"""
xirr(tim,cf) = Roots.find_zero(r -> xnpv(r,tim,cf,Dates.Date(2020,1,1)),(0.0,0.5))

"""
`xirr(cf) = IRR of cash flow against time given in Dates`

* cf  = vector of tuple of (time in Dates, cash flows)
"""
xirr(cf) = Roots.find_zero(r -> xnpv(r,cf,Dates.Date(2020,1,1)),(0.0,0.5))

"""
`sharpe(rf,ra,sa) = (ra - rf)/sa`

Sharpe ratio where:

* rf = risk free rate
* ra = rate of return of portfolio 'a'
* sa = std dev of portfolio 'a'
"""
sharpe(rf,ra,sa) = (ra - rf)/sa

"""
`rBeta(rf,rm,r) = (r - rf)/(rm - rf)`

Map from rate to beta given a rf and rm
"""
rBeta(rf,rm,r) = (r - rf)/(rm - rf)

"""
`betaR(rf,rm,b) = rf + b*(rm - rf)`

Map from beta to rate given a rf and rm
"""
betaR(rf,rm,b) = rf + b*(rm - rf)

abstract type RateType end

abstract type NomRate <: RateType end
abstract type EffRate <: RateType end
abstract type ExpRate <: RateType end

"""
`RateCurve{T<:RateType} : struct having rates at different periods`
* T    = Can be NomRate, EffRate or ExpRate
* rate = Vector of rates at periods
* freq = Frequency at which rates are represented

Thus if freq = 2, then rate[3] = rate @ 1.5 years
"""
struct RateCurve{T<:RateType}
  rate :: Vector{Float64}
  freq :: Int64
end

"""
`DiscountFactor : struct having discount factor of all period`
* factor  = Vector of factor at periods
* freq    = Frequency at which factors are represented
"""
struct DiscountFactor
  factor  :: Vector{Float64}
  freq    :: Int64
end

"""
`estimR(rt, y) = estimate rate from a RateCurve and period through Interpolation`

* rt  = RateCurve
* y   = period for which rate is to be estimated
"""
estimR(rt::RateCurve{NomRate}, y::Float64) = estimR(rt.rate, rt.freq, y)
estimR(rt::RateCurve{EffRate}, y::Float64) = estimR(rt.rate, rt.freq, y)
estimR(rt::RateCurve{ExpRate}, y::Float64) = estimR(rt.rate, rt.freq, y)

function estimR(rx, fq, y)
  pt = y*fq; fl = unsafe_trunc(Int64, pt); f0 = Int64(fl); r0 = rx[f0]
  fl == pt ? r0 : r0 + (rx[f0+1] - r0)*(pt-fl)
end

"""
`discFactorToNominalRate(dF) = Convert Discount rate dF to Spot rate`
"""
discFactorToNominalRate(dF) = RateCurve{NomRate}( (i -> ((1/dF.factor[i])^(1/i) - 1)*dF.freq).(axes(dF.factor,1)), dF.freq)

effR(rC::RateCurve{NomRate}) = RateCurve{EffRate}((x -> effR(x, rC.freq)).(rC.rate), rC.freq)

effR(rC::RateCurve{ExpRate}) = RateCurve{EffRate}((x -> expToEffR(x)).(rC.rate), rC.freq)

nomR(rC::RateCurve{EffRate}) = RateCurve{NomRate}((x->nomR(x,rC.freq)).(rC.rate),rC.freq)

nomR(rC::RateCurve{ExpRate}) = RateCurve{NomRate}((x->expToNomR(x,rC.freq)).(rC.rate),rC.freq)

expR(rC::RateCurve{NomRate}) = RateCurve{ExpRate}((x -> expR(x, rC.freq)).(rC.rate), rC.freq)

expR(rC::RateCurve{EffRate}) = RateCurve{ExpRate}((x -> expR(x)).(rC.rate), rC.freq)

include("FixedIncomes/mod.jl")
include("Derivatives/mod.jl")
include("Statements/mod.jl")

end # module


text:
# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.

"""
    repair!(instance)

Verifies that the given unit commitment instance is valid and automatically
fixes some validation errors if possible, issuing a warning for each error
found. If a validation error cannot be automatically fixed, issues an
exception.

Returns the number of validation errors found.
"""
function repair!(instance::UnitCommitmentInstance)::Int
    n_errors = 0

    for g in instance.units

        # Startup costs and delays must be increasing
        for s in 2:length(g.startup_categories)
            if g.startup_categories[s].delay <= g.startup_categories[s-1].delay
                prev_value = g.startup_categories[s].delay
                new_value = g.startup_categories[s-1].delay + 1
                @warn "Generator $(g.name) has non-increasing startup delays (category $s). " *
                      "Changing delay: $prev_value → $new_value"
                g.startup_categories[s].delay = new_value
                n_errors += 1
            end

            if g.startup_categories[s].cost < g.startup_categories[s-1].cost
                prev_value = g.startup_categories[s].cost
                new_value = g.startup_categories[s-1].cost
                @warn "Generator $(g.name) has decreasing startup cost (category $s). " *
                      "Changing cost: $prev_value → $new_value"
                g.startup_categories[s].cost = new_value
                n_errors += 1
            end
        end

        for t in 1:instance.time
            # Production cost curve should be convex
            for k in 2:length(g.cost_segments)
                cost = g.cost_segments[k].cost[t]
                min_cost = g.cost_segments[k-1].cost[t]
                if cost < min_cost - 1e-5
                    @warn "Generator $(g.name) has non-convex production cost curve " *
                          "(segment $k, time $t). Changing cost: $cost → $min_cost"
                    g.cost_segments[k].cost[t] = min_cost
                    n_errors += 1
                end
            end

            # Startup limit must be greater than min_power
            if g.startup_limit < g.min_power[t]
                new_limit = g.min_power[t]
                prev_limit = g.startup_limit
                @warn "Generator $(g.name) has startup limit lower than minimum power. " *
                      "Changing startup limit: $prev_limit → $new_limit"
                g.startup_limit = new_limit
                n_errors += 1
            end
        end
    end

    return n_errors
end

export repair!


text:
"""
	inbin!(traj::Trajectory, fb::FrequencyBin)

Determine whether `traj` goes through the bin `fb` at any time point.
Set the `timetofreq` field of `traj` accordingly.
"""
function inbin!(traj::Trajectory, fb::FrequencyBin)
	for (t,f) in traj
		if in(f, fb)
			traj.timetofreq[fb] = t
			return true
		end
	end
	return false
end
"""
	inbin(traj::Trajectory, fb::FrequencyBin)

Determine whether `traj` goes through the bin `fb` at any time point.
"""
function inbin(traj, fb)
	for (t,f) in traj
		if in(f, fb)
			return true
		end
	end
	return false
end

function filter!(trajectories::Vector{Trajectory}, fb::FrequencyBin)
	idx = findall(traj -> inbin!(traj, fb), trajectories)
	return trajectories[idx]
end
filter!(trajectories, f, df) = filter!(trajectories, FrequencyBin(f, df))

"""
	filter(f::Function, trajectories::Vector{Trajectory})
"""
function filter(f::Function, trajectories::Vector{Trajectory})
	idx = findall(f, trajectories)
	return trajectories[idx]
end
"""
	filter(trajectories::Vector{Trajectory}, fb::FrequencyBin)
"""
function filter(trajectories::Vector{Trajectory}, fb::FrequencyBin)
	return filter(traj -> inbin(traj, fb), trajectories)
end

"""
	has_fitness_above(traj, ϕmin)

Does `traj` have a fitness *strictly* above `ϕmin` for all its duration?
"""
function has_fitness_above(traj::Trajectory, ϕmin)
	for ϕ in traj.ϕpos
		if ismissing(ϕ) || ϕ <= ϕmin
			return false
		end
	end
	return true
end
"""
	has_fitness_above(traj, ϕmin, fb::FrequencyBin)

Does `traj` have a fitness *strictly* above `ϕmin` when it enters `fb`?
Return `false` if `traj` never enters `fb`.
"""
function has_fitness_above(traj::Trajectory, ϕmin::Float64, fb::FrequencyBin)
	for (i, (t,f)) in enumerate(traj)
		if in(f, fb)
			return traj.ϕpos[i] > ϕmin
		end
	end
	return false
end


text:
module test_julia_egal

using Test

@test 1 == 1.0

@test "" === ""
@test 1 !== 1.0
@test [] !== []

struct A
    val
end

@test A("") === A("")
@test A([]) !== A([])

end # module test_julia_egal


text:
export
abs,
align,
any,
argsort,
clip,
clip_lower,
clip_upper,
corr,
corrwith,
count,
cov,
cummax,
cummin,
cumprod,
cumsum,
drop_duplicates,
duplicated,
filter,
first,
head,
hist,
idxmax,
idxmin,
iloc,
index,
join,
last,
loc,
plot,
reindex,
reindex_axis,
reindex_like,
rename,
resample,
reset_index,
rolling_count,
select,
set_index,
sort_index,
sum,
tail,
to_csv,
to_latex,
truncate,
xs,
GroupBy,
Rolling,
Iloc,
Index,
Loc,
MultiIndex,
Series,
agg,
aggregate,
get_group,
groups,
indices,
mean,
median,
ohlc,
pivot_table,
read_csv,
std,
transform,
var,
melt,
@>,
@query,
app,
columns,
concat,
crosstab,
cut,
deletecolumn!,
ewma,
expanding_apply,
expanding_corr,
expanding_corr_pairwise,
expanding_count,
expanding_cov,
expanding_kurt,
expanding_max,
expanding_mean,
expanding_median,
expanding_min,
expanding_quantile,
expanding_skew,
expanding_std,
expanding_sum,
expanding_var,
from_arrays,
from_tuples,
get_dummies,
groupby,
rolling,
index!,
isin,
name,
name!,
qcut,
read_excel,
read_html,
read_json,
read_pickle,
read_table,
rolling_apply,
rolling_corr,
rolling_corr_pairwise,
rolling_count,
rolling_cov,
rolling_kurt,
rolling_max,
rolling_mean,
rolling_median,
rolling_min,
rolling_quantile,
rolling_quantile,
rolling_skew,
rolling_std,
rolling_sum,
rolling_sum,
rolling_var,
rolling_window,
sample,
save,
setcolumns!,
setname!,
siz,
stats,
to_clipboard,
to_numeric,
values,
DataFrame,
T,
asfreq,
asof,
bdate_range,
date_range,
delevel,
describe,
diff,
dropna,
ewma,
ewmcorr,
ewmcov,
ewmstd,
ewmvar,
fillna,
first_valid_index,
last_valid_index,
mean,
median,
min,
mode,
order,
pct_change,
period_range,
pivot,
quantile,
query,
rank,
reodrer_levels,
replace,
resample,
shift,
skew,
sort,
sortlevel,
stack,
std,
swaplevel,
to_datetime,
to_dense,
to_dict,
to_excel,
to_gbq,
to_hdf,
to_html,
to_json,
to_msgpack,
to_panel,
to_pickle,
to_records,
to_sparse,
to_sql,
to_string,
to_timedelta,
tz_conert,
tz_localize,
unstack,
var,
weekday,
read_sql,
read_hdf,
HDFStore,
info,
put,
walk

if !isdefined(Base, :drop)
    export drop
end

if !isdefined(Base, :take)
    export take
end


text:
export save_gen, load_gen!, string

## NEATIndiv to JSON
function get_name(func)
    for (name, f) in pairs(activ_dict)
        if f == func
            return name
        end
    end
end

function get_name(dict::Dict)
    new_dict = Dict()
    for (k, v) in pairs(dict)
        new_dict[k]=get_name(v)
    end
    new_dict
end

function fitness_value(f::Float64)
	if f == -Inf
		"-Inf"
	end
	string(f)
end

function Base.string(indiv::NEATIndividual)
    d = Dict(
	"Type"=>"NEAT",
    "Genes"=>indiv.genes,
    "Fitness"=>fitness_value.(indiv.fitness),
    "Activ_func"=>get_name(indiv.activ_functions)
    )
    JSON.json(d)
end

function Base.string(indiv::HyperNEATIndividual)
    d = Dict(
	"Type"=>"HyperNEAT",
    "Genes"=>indiv.genes,
    "Fitness"=>fitness_value.(indiv.fitness),
    "Activ_func"=>get_name(indiv.activ_functions)
    )
    JSON.json(d)
end

## Load
"""
Function to load an entire population from a gen folder.
The evolution need to be initialized first.
"""
function load_gen!(e::Evolution, path::String)
	individualNameList = readdir("gens/$path")
	individualList::Array{NEATIndiv} = []
	for i in eachindex(individualNameList)
		indString = read("gens/$path/$(individualNameList[i])", String)
		ind = NEATIndiv(e.cfg, indString)
		build!(ind)
		push!(individualList,ind)
	end
	e.population = individualList
end


text:
module TestRandom
    using Test
    using BioEnergeticFoodWebs

    @test_throws AssertionError nichemodel(10, 150)

    A = nichemodel(10, 20)
    @test size(A, 1) == 10
    @test size(A, 2) == 10

    A = nichemodel(10, 0.12)
    @test size(A, 1) == 10
    @test size(A, 2) == 10

    A = nichemodel(10, 0.12, toltype=:rel)
    A = nichemodel(10, 0.12, toltype=:abs)

    A = [0 1; 1 0]
    @test BioEnergeticFoodWebs.connectance(A) == 0.5

end

max_stars_repo_name: SimmonsBI/BioEnergeticFoodWebs

max_issues_repo_name: SimmonsBI/BioEnergeticFoodWebs

max_forks_repo_name: SimmonsBI/BioEnergeticFoodWebs


text:
export MHT, MHT1, MHT2, MHT3,MHT4,MHT5,MHT6,MHT7;
export CheapFunction, ExpensiveFunction;

# Wrappers for Heterogenous Problems
struct CheapFunction <: Function
    f :: Function
end

struct ExpensiveFunction <: Function
    f :: Function
end

(F::CheapFunction)( args...; kwargs... ) = F.f( args...; kwargs... );
(F::ExpensiveFunction)( args...; kwargs... ) = F.f( args...; kwargs...);

####### General Setup ###########
mht_ref_string = "Thomann, 2019,
“Doctoral Thesis: A Trust Region Approach for
Multi-Objective Heterogenous Optimization”";
@doc """
    abstract type MHT <: MOP end

Abstract super type for the set of test problems MHT1-MHT7
(called T1-T7 in the source) from
$(mht_ref_string)
"""
abstract type MHT <: MOP end;

########## MHT1 ##############
@doc """
    MHT1()

Bi-Objective, bi-variate problem MHT1 (T1) from
$(mht_ref_string).

Unconstrained, convex objectives, convex Pareto Front.
``f₁`` expensive, ``f₂`` cheap.

The objectives are given as
```math
\\begin{aligned}
    f_1( \\mathbf{x} ) &= \\frac{1}{2}x_1^2 + x_2^2 - 10 x_1 - 100,\\\\
    f_2( \\mathbf{x} ) &= x_1^2 + \\frac{1}{2}x_2^2 - 10 x_2 - 100.
\\end{aligned}
```
"""
struct MHT1 <: MHT end;

constraints( mop :: MHT1 ) = nothing;
num_vars( mop :: MHT1 ) = 2;
num_objectives( mop :: MHT1 ) = 2;

function get_objectives( mop :: MHT1 )
    f1 = ExpensiveFunction(
        function( x :: Vector{R} where R<:Real)
            0.5 * x[1]^2 + x[2]^2 - 10 * x[1] - 100
        end
    );
    f2 = CheapFunction(
        function( x :: Vector{R} where R<:Real)
            x[1]^2 + 0.5 * x[2]^2 - 10 * x[2] - 100
        end
    );
    return [f1, f2]
end

get_ideal_point( mop :: MHT1 ) = fill(-150.0, 2);

mht1_points = [
    [-34.24; 47.06],
    [45.72; -1.46],
    [30.03; -35.81],
    [-7.82; 41.57],
    [29.22; 45.95],
    [15.75; -46.43],
    [34.91; 43.4],
    [17.87; 25.77],
    [24.31; -10.78],
    [15.55; -32.88],
]
function get_starting_point( mop :: MHT1, N :: Int64 )
    n = 1 + (N-1) % 10
    return mht1_points[n]
end

########## MHT2 ##############
@doc """
    MHT2()

Bi-Objective, bi-variate problem MHT2 (T2) from
$(mht_ref_string).

Unconstrained, non-convex objectives, non-convex Pareto Front.
``f₁`` cheap, ``f₂`` expensive.

The objectives are given as
```math
\\begin{aligned}
    f_1( \\mathbf{x} ) &= \\sin( x₂ ),\\\\
    f_2( \\mathbf{x} ) &= 1 -
        \\exp \\left(
            - \\left( x₁ - \\frac{1}{\\sqrt{2}} \\right)^2
            - \\left( x₂ - \\frac{1}{\\sqrt{2}} \\right)^2
        \\right)
\\end{aligned}
```
"""
struct MHT2 <: MHT end;

constraints( mop :: MHT2 ) = nothing;
num_vars( mop :: MHT2 ) = 2;
num_objectives( mop :: MHT2 ) = 2;

function get_objectives( mop :: MHT2 )
    f1 = CheapFunction(
        function( x :: Vector{R} where R<:Real)
            sin( x[2] )
        end
    );
    f2 = ExpensiveFunction(
        function( x :: Vector{R} where R<:Real)
            1 - exp( -( x[1] - 1/sqrt(2) )^2 - ( x[2] - 1/sqrt(2) )^2 )
        end
    );
    return [f1, f2]
end

get_ideal_point( mop :: MHT2 ) = [-1.0; 0.0];

mht2_points = [
    [20.6; -46.82],
    [-22.31; -45.38],
    [-40.29; 32.35],
    [-.4; 1.0],
    [45.02; -46.56],
    [-6.13; -11.84],
    [-0.01; -0.51],
    [-31.31; -1.02],
    [0.51; 1.00],
    [.9; 1.51],
]
get_starting_point( ::MHT2, N :: Int64) = mht2_points[ 1 + (N-1) % 10 ];

########## MHT3 ##############
@doc """
    MHT3()

Bi-Objective, bi-variate problem MHT3 (T3) from
$(mht_ref_string).

Box constraints [-2,2]², convex objectives, convex Pareto Front.
``f₁`` cheap, ``f₂`` expensive.

The objectives are given as
```math
\\begin{aligned}
    f_1( \\mathbf{x} ) &= x₁ + 2,\\\\
    f_2( \\mathbf{x} ) &= x₁ - 2 + x₂.
\\end{aligned}
```
"""
struct MHT3 <: MHT end;

constraints( mop :: MHT3 ) = Box(
    fill(-2.0, 2),
    fill(2.0, 2)
);
num_vars( mop :: MHT3 ) = 2;
num_objectives( mop :: MHT3 ) = 2;

function get_objectives( mop :: MHT3 )
    f1 = CheapFunction(
        function( x :: Vector{R} where R<:Real)
            x[1] + 2
        end
    );
    f2 = ExpensiveFunction(
        function( x :: Vector{R} where R<:Real)
            x[1] - 2 + x[2]
        end
    );
    return [f1, f2]
end

get_ideal_point( mop :: MHT3 ) = [ 0.0; -6.0];

mht3_points = [
    2.5767 -0.0793 -0.3193 0.0511 1.9058 0.8659 1.8695 -0.8956 2.2557 0.7349;
    1.6543 -0.3848 -1.1619 0.0646 1.7690 -0.7283 0.197 2.634 0.3009 0.5223
]
get_starting_point( ::MHT3, N :: Int64) = mht3_points[:,1 + (N-1) % 10 ];

########## MHT4 ##############
@doc """
    MHT4(n)

Bi-Objective, `n`-variate problem MHT4 (T4) from
$(mht_ref_string).

Box constraints [-10,10]ⁿ, convex objectives, convex Pareto Front.
``f₁`` expensive, ``f₂`` cheap.

The objectives are given as
```math
\\begin{aligned}
    f_1( \\mathbf{x} ) &= Σ_{i=1}^{n-1} xᵢ² + 2, \\\\
    f_2( \\mathbf{x} ) &= Σ_{i=1}^{n} xᵢ - 2
\\end{aligned}
```
"""
@with_kw struct MHT4 <: MHT
    n_vars :: Int64 = 2;
end;

constraints( mop :: MHT4 ) = Box(
    fill(-10.0, mop.n_vars),
    fill(10.0, mop.n_vars)
);
num_vars( mop :: MHT4 ) = mop.n_vars;
num_objectives( mop :: MHT4 ) = 2;

function get_objectives( mop :: MHT4 )
    f1 = ExpensiveFunction(
        function( x :: Vector{R} where R<:Real )
            sum( x[1:end-1].^2 .+ 2.0 )
        end
    );
    f2 = CheapFunction(
        function( x :: Vector{R} where R<:Real)
            sum( x .- 2.0 )
        end
    );
    return [f1, f2]
end

get_ideal_point( mop :: MHT4 ) = [
    (mop.n_vars - 1)*2.0;
    -12.0 * mop.n_vars
];

mht4_points = [
    8.2401 -5.8961 0.0861 3.7051 -0.4188 2.4755 -9.0298 9.8070 2.5461 -8.1104;
    2.7963 9.4970 -4.333 -8.8077 0.7261 6.2651 3.0539 -0.4131 5.1145 -2.6839
]
function get_starting_point( mop::MHT4, N :: Int64)
    if mop.n_vars == 2
        return mht4_points[:,1 + (N-1) % 10 ];
    else
        return get_random_point( mop )
    end
end

########## MHT5 ##############
@doc """
    MHT5()

Bi-Objective, bi-variate problem MHT5 (T5) from
$(mht_ref_string).

Box constraints [1e-12,30] × [0,30],
convex objectives, convex Pareto Front.

``f₁`` expensive, ``f₂`` cheap.

The objectives are given as
```math
\\begin{aligned}
    f_1( \\mathbf{x} ) &= x₁ \\ln( x₁ ) + x₂²,\\\\
    f_2( \\mathbf{x} ) &= x₁² + x₂⁴.
\\end{aligned}
```
"""
struct MHT5 <: MHT end;

constraints( mop :: MHT5 ) = Box(
    [1e-12; 0.0],
    [30.0; 30.0]
);
num_vars( mop :: MHT5 ) = 2;
num_objectives( mop :: MHT5 ) = 2;

function get_objectives( mop :: MHT5 )
    f1 = ExpensiveFunction(
        function( x :: Vector{R} where R<:Real)
            x[1] * log( x[1] ) + x[2]^2
        end
    );
    f2 = CheapFunction(
        function( x :: Vector{R} where R<:Real)
            x[1]^2 + x[2]^4
        end
    );
    return [f1, f2]
end

get_ideal_point( mop :: MHT5 ) = [ -1/ℯ; 0.0];

mht5_points = [
    25.9816 4.2130 7.7573 2.1431 10.4985 22.1153 18.7311 29.6945 25.6892 14.5115;
    2.6608 19.5013 24.4573 8.8787 6.5516 9.8276 9.1069 23.0979 25.2217 3.8922
]
get_starting_point( ::MHT5, N :: Int64) = mht5_points[:,1 + (N-1) % 10 ];

########## MHT6 ##############
@doc """
    MHT6()

Bi-Objective, bi-variate problem MHT6 (T6) from
$(mht_ref_string).

Box constraints [1e-12,100]²,
convex objectives, convex Pareto Front.

``f₁`` expensive, ``f₂`` cheap.

The objectives are given as
```math
\\begin{aligned}
    f_1( \\mathbf{x} ) &= -\\ln( x₁ ) - \\ln( x₂ ),\\\\
    f_2( \\mathbf{x} ) &= x₁² + x₂.
\\end{aligned}
```
"""
struct MHT6 <: MHT end;

constraints( mop :: MHT6 ) = Box(
   fill(1e-12,2),
   fill(100.0,2)
);

num_vars( mop :: MHT6 ) = 2;
num_objectives( mop :: MHT6 ) = 2;

function get_objectives( mop :: MHT6 )
    f1 = ExpensiveFunction(
        function( x :: Vector{R} where R<:Real)
            - log( x[1] ) - log( x[2] )
        end
    );
    f2 = CheapFunction(
        function( x :: Vector{R} where R<:Real)
            x[1]^2 + x[2]
        end
    );
    return [f1, f2]
end

get_ideal_point( mop :: MHT6 ) = [-2*log(100);0.0];

mht6_points = [
   9.5896 95.6433 61.6860 23.3495 44.4905 5.8716 72.6812 25.441 42.2296 44.4367;
   9.4335 20.8035 50.6963 6.1688 2.7671 52.5858 58.64 7.8104 35.9188 17.2381
]
get_starting_point( ::MHT6, N :: Int64) = mht6_points[:,1 + (N-1) % 10 ];

########## MHT7 ##############
@doc """
    MHT7()

Bi-Objective, tri-variate problem MHT7 (T7) from
$(mht_ref_string).

Box constraints [0, 30]^3
convex objectives, convex Pareto Front.
9
``f₁`` expensive, ``f₂`` cheap.

The objectives are given as
```math
\\begin{aligned}
    f_1( \\mathbf{x} ) &= Σᵢ xᵢ⁴ + Σᵢ xᵢ³,\\\\
    f_2( \\mathbf{x} ) &= Σᵢ xᵢ
\\end{aligned}
```
"""
struct MHT7 <: MHT end;

constraints( mop :: MHT7 ) = Box(
   zeros(3),
   fill(30.0,3)
);

num_vars( mop :: MHT7 ) = 3;
num_objectives( mop :: MHT7 ) = 2;

function get_objectives( mop :: MHT7 )
    f1 = ExpensiveFunction(
        function( x :: Vector{R} where R<:Real)
            sum( x.^4 ) + sum( x.^3 )
        end
    );
    f2 = CheapFunction(
        function( x :: Vector{R} where R<:Real)
            sum( x )
        end
    );
    return [f1, f2]
end

get_ideal_point( mop :: MHT7 ) = zeros(2);

mht7_points = [
    24.4417 27.4013 8.3549 28.9467 28.7150 4.2566 23.7662 1.0714 20.3621 11.7668;
    27.1738 18.9708 16.4064 4.7284 14.5613 12.6528 28.7848 25.4739 22.7322 19.6643;
    3.8096 2.9262 28.7252 29.1178 24.0084 27.4721 19.6722 28.0198 22.2940 5.1356
];
get_starting_point( ::MHT7, N :: Int64) = mht7_points[:,1 + (N-1) % 10 ];

##########################################################################
# Pareto Sets when individual minima / other points are connected by a line

function get_pareto_set( mop :: M where M<:MHT )
    return SamplingFunction( mop, :ParetoSet, [:regular, :random] );
end
function get_pareto_front( mop :: M where M<:MHT )
    return SamplingFunction( mop, :ParetoFront, [:regular, :random] );
end

function individual_minima(::MHT1)
    x_min1 = [ 10.0 ;  0 ];
    x_min2 = [ 0; 10.0 ];
    x_diff = x_min2 .- x_min1
    return x_min1, x_min2, x_diff
end

function individual_minima(::MHT2)
    x_min1 = 1/√2 .* ones(2);
    x_min2 = [1/√2; -π/2];
    x_diff = x_min2 .- x_min1
    return x_min1, x_min2, x_diff
end

function individual_minima(::MHT3)
    x_min1 = [ -2.0 ;  2.0 ];
    x_min2 = [ -2.0 ; -2.0 ];
    x_diff = x_min2 .- x_min1
    return x_min1, x_min2, x_diff
end

function individual_minima(mop::MHT4)
    x_min1 = [ zeros(mop.n_vars-1); -10.0 ];
    x_min2 = fill( -10.0, mop.n_vars );
    x_diff = x_min2 .- x_min1
    return x_min1, x_min2, x_diff
end

function individual_minima(mop::MHT5)
    x_min1 = [ 1/ℯ, 0.0];
    x_min2 = [1e-12, 0.0];
    x_diff = x_min2 .- x_min1
    return x_min1, x_min2, x_diff
end

function get_points( sample_func :: SamplingFunction, mop :: Union{MHT1,MHT2,MHT3,MHT4,MHT5},
    ::Val{:ParetoSet}, num_points :: Int; method = :regular )

    if method == :regular
        t_range = num_points == 1 ? [0.5,] : range(0, 1.0; length = num_points );
    elseif method == :random
        t_range = rand(num_points);
    end

    x_min1, x_min2, x_diff = individual_minima( mop );

    return [ x_min1 .+ τ .* x_diff for τ in t_range ]
end

function get_points( sample_func :: SamplingFunction, mop :: MHT6,
    ::Val{:ParetoSet}, num_points :: Int; method = :regular )
    α_min = 1/101;
    α_max = 1/(1e-12 + 1);

    num_p1 = floor(Int, num_points/2);
    num_p2 = num_points - num_p1;

    if method == :regular
        α_range = num_p1 == 1 ? [0.5,] : range(α_min, α_max; length = num_p1 );
    elseif method == :random
        α_range = max.( α_min .+ rand(num_p1), α_max ) ;
    end

    ps₁ = [ [sqrt( (1-α)/(2*α) ); 1/α-1] for α ∈ α_range ];

    α_min₂ = 1/(2*100^2 + 1);
    α_range₂ = num_p2 == 1 ? [0.5,] : range(α_min₂, α_min; length = num_p2 );

    ps₂ = [ [sqrt( (1-α)/(2*α) ); 100.0] for α ∈ α_range₂ ];

    return [ ps₁; ps₂ ]
end

function get_points( sample_func :: SamplingFunction, mop :: MHT7,
    ::Val{:ParetoSet}, num_points :: Int; method = :regular )
    @warn "Problem has a single Pareto optimum.
    Returning it $(num_points) times for compatibility."
    return [ zeros(2) for i = 1 : num_points ];
end

function get_points( sample_func :: SamplingFunction, mop :: Union{MHT1,MHT2,MHT3,MHT4,MHT5,MHT6,MHT7},
    ::Val{:ParetoFront}, num_points :: Int; method = :regular )

    pset = get_points( sample_func, mop, Val(:ParetoSet), num_points; method );
    F = get_vector_objective( mop );
    return F.( pset );
end


text:
using RoME
using Statistics
# , Distributions
using Test

import Base: convert

##

@testset "test sampling from BearingRange factor..." begin

##

p2br = Pose2Point2BearingRange(Normal(0,0.1),Normal(20.0,1.0))

fg = initfg()
addVariable!(fg, :x0, Pose2)
addVariable!(fg, :x1, Point2)
addFactor!(fg, [:x0;:x1], p2br, graphinit=false)

meas = sampleFactor(IIF._getCCW(fg, :x0x1f1), 100)

##

# meas = getSample(p2br, 100)
@test abs(Statistics.mean(meas[1][1,:])) < 0.1
@test 0.05 < abs(Statistics.std(meas[1][1,:])) < 0.2

@test abs(Statistics.mean(meas[1][2,:]) - 20.0) < 1.0
@test 0.5 < abs(Statistics.std(meas[1][2,:])) < 1.5

##

end

@testset "test BearingRange factor residual function..." begin

##

# dummy variables
fg = initfg()
X0 = addVariable!(fg, :x0, Pose2)
X1 = addVariable!(fg, :x1, Point2)

##

p2br = Pose2Point2BearingRange(Normal(0,0.1),Normal(20.0,1.0))

xi = zeros(3)
li = zeros(2); li[1] = 20.0;
zi = (zeros(2),); zi[1][2] = 20.0

res = testFactorResidualBinary( p2br,
                                Pose2,
                                Point2,
                                xi,
                                li,
                                zi )
#

@show res
@test norm(res) < 1e-14

##

xi = zeros(3)
li = zeros(2); li[2] = 20.0;
zi = (zeros(2),); zi[1][:] = [pi/2;20.0]

# idx = 1
# res = zeros(2)
# p2br(res, fmd, idx, zi, xi, li)

res = testFactorResidualBinary( p2br,
                                Pose2,
                                Point2,
                                xi,
                                li,
                                zi )

@show res
@test norm(res) < 1e-14

xi = zeros(3); xi[3] = pi/2
li = zeros(2); li[2] = 20.0;
zi = (zeros(2),); zi[1][:] = [0.0;20.0]

res = testFactorResidualBinary( p2br,
                                Pose2,
                                Point2,
                                xi,
                                li,
                                zi )
#
@show res
@test norm(res) < 1e-14

##

xi = zeros(3); xi[3] = -pi/2
li = zeros(2); li[1] = 20.0;
# zi = ([0.0;pi/2],[0.0;20.0],)
zi = (zeros(2),); zi[1][:] = [pi/2;20.0]

# idx = 2
# res = zeros(2)
# p2br(res, fmd, idx, zi, xi, li)

res = testFactorResidualBinary( p2br,
                                Pose2,
                                Point2,
                                xi,
                                li,
                                zi )

@show res
@test norm(res) < 1e-14

##
# #TODO Update to new CalcFactor
# # test parametric Pose2Point2BearingRange
# f = Pose2Point2BearingRange(Normal(0.0,1), Normal(10.0,1))
# @test isapprox(f([0.,0,0], [10.,0]), 0, atol = 1e-9)
# @test isapprox(f([0,0,pi/2], [0.,10]), 0, atol = 1e-9)

# f = Pose2Point2BearingRange(Normal(pi/2,1), Normal(10.0,1))
# @test isapprox(f([0.,0,0], [0.,10]), 0, atol = 1e-9)
# @test isapprox(f([0,0,pi/2], [-10.,0]), 0, atol = 1e-9)

# f = Pose2Point2BearingRange(Normal(pi,1), Normal(10.0,1))
# @test isapprox(f([0.,0,0], [-10.,0]), 0, atol = 1e-9)
# @test isapprox(f([0,0,pi/2], [0.,-10]), 0, atol = 1e-9)

# f = Pose2Point2BearingRange(Normal(-pi/2,1), Normal(10.0,1))
# @test isapprox(f([0.,0,0], [0.,-10]), 0, atol = 1e-9)
# @test isapprox(f([0,0,pi/2], [10.,0]), 0, atol = 1e-9)

end

@testset "test unimodal bearing range factor, solve for landmark..." begin

##

# Start with an empty graph
N = 1
fg = initfg()

#add pose with partial constraint
addVariable!(fg, :x0, Pose2)
addFactor!(fg, [:x0], PriorPose2(MvNormal(zeros(3), 0.01*Matrix{Float64}(LinearAlgebra.I, 3,3))), graphinit=false)
# force particular initialization
setVal!(fg, :x0, zeros(3,1))

##----------- sanity check that predictbelief plumbing is doing the right thing
pts, = predictbelief(fg, :x0, ls(fg, :x0), N=75)
@test sum(abs.(Statistics.mean(pts,dims=2)) .< [0.1; 0.1; 0.1]) == 3
@test sum([0.05; 0.05; 0.05] .< Statistics.std(pts,dims=2) .< [0.15; 0.15; 0.15]) == 3
#------------

# Add landmark
addVariable!(fg, :l1, Point2, tags=[:LANDMARK;])
li = zeros(2,1); li[1,1] = 20.0;
setVal!(fg, :l1, li)

# Add bearing range measurement between pose and landmark
p2br = Pose2Point2BearingRange(Normal(0,0.1),Normal(20.0,1.0))
addFactor!(fg, [:x0; :l1], p2br, graphinit=false)

# there should be just one (the bearingrange) factor connected to :l1
@test length(ls(fg, :l1)) == 1
# drawGraph(fg, show=true)

# check the forward convolution is working properly
pts, = predictbelief(fg, :l1, ls(fg, :l1), N=75)
@test sum(abs.(Statistics.mean(pts,dims=2) - [20.0; 0.0]) .< [2.0; 2.0]) == 2
@test sum([0.1; 0.1] .< Statistics.std(pts,dims=2) .< [3.0; 3.0]) == 2

# using Gadfly, KernelDensityEstimate, KernelDensityEstimatePlotting
#
# pl = plotKDE(kde!(pts))
# pl.coord = Coord.Cartesian(xmin=-5,xmax=25, ymin=-10.0,ymax=10)
# pl

##

end

@testset "test unimodal bearing range factor, solve for pose..." begin

##

# Start with an empty graph
N = 1
fg = initfg()

# Add landmark
addVariable!(fg, :l1, Point2, tags=[:LANDMARK;])
addFactor!(fg, [:l1], PriorPoint2(MvNormal([20.0;0.0], Matrix(Diagonal([1.0;1.0].^2)))),  graphinit=false ) # could be IIF.Prior
li = zeros(2,1); li[1,1] = 20.0;
setVal!(fg, :l1, li)

#add pose with partial constraint
addVariable!(fg, :x0, Pose2)
# force particular initialization
setVal!(fg, :x0, zeros(3,1))

# Add bearing range measurement between pose and landmark
p2br = Pose2Point2BearingRange(Normal(0,0.1),Normal(20.0,1.0))
addFactor!(fg, [:x0; :l1], p2br, graphinit=false)

# there should be just one (the bearingrange) factor connected to :l1
@test length(ls(fg, :x0)) == 1
# writeGraphPdf(fg)

# check the forward convolution is working properly
pts, = predictbelief(fg, :x0, ls(fg, :x0), N=75)

pts[3,:] .= TU.wrapRad.(pts[3,:])
@show abs.(Statistics.mean(pts,dims=2))
@test sum(abs.(Statistics.mean(pts,dims=2)) .< [2.0; 2.0; 2.0]) == 3
@show Statistics.std(pts,dims=2)
@test sum([0.1; 2.0; 0.01] .< Statistics.std(pts,dims=2) .< [5.0; 10.0; 2.0]) == 3

##

end

@testset "Testing Pose2Point2Bearing Initialization and Packing" begin

##

p2p2b = Pose2Point2Bearing( MvNormal([0.2,0.2,0.2], [1.0 0 0;0 1 0;0 0 1]) )
packed = convert(PackedPose2Point2Bearing, p2p2b)
p2p2bTest = convert(Pose2Point2Bearing, packed)
@test p2p2b.bearing.μ == p2p2bTest.bearing.μ
@test p2p2b.bearing.Σ.mat == p2p2bTest.bearing.Σ.mat

##

end


text:
using TypeDBClient
using DataFrames
using CSV
using Profile
using ProfileSVG
using ProfileView
using Dates

using TypeDBClient: AbstractCoreSession,
                    CoreClient,
                    CoreSession,
                    AbstractCoreTransaction,
                    Proto,
                    commit,
                    insert

include("define_schema.jl")
include("utils.jl")

database_name = "biograkn"

client = CoreClient("127.0.0.1",1729)

delete_db_create_db(client, database_name)
@info "Deleted all databases, create database $database_name and the schema"

function uniprot(client::TypeDBClient.AbstractCoreClient, database_name::String, make_inserts::Bool=true)
    # Uniprot
        @info "1. Uniprot"
        df = CSV.File(joinpath(@__DIR__, "../dataset/uniprot/uniprot.csv"), delim="\t") |> DataFrame

        start_preparing = time()
        query_strings = String[]
        data_pair = Pair[]
        for row in eachrow(df)
            tmp_string = " \$_ isa protein, has uniprot-id \"$(row[1])\", has uniprot-name \"$(row[4])\";"
            push!(query_strings, tmp_string)
            push!(data_pair, (row[1]=>row[4]))
        end
        end_preparing = time()

        @info "Time for preparing: $(end_preparing - start_preparing)"

        sess = CoreSession(client, database_name , Proto.Session_Type.DATA, request_timeout=Inf)

        time1 = time()
        if make_inserts
            res_insert = insert_data(sess, query_strings, 100, 50)
        end
        time2 = time()

        @info "Uniprot took $(round(time2 - time1, digits = 2)) seconds"

        return res_insert
end

function transaction_insert(trans::AbstractCoreTransaction,
            data,
            count_channel::Channel,
            nr_channel,
            error_channel::Channel)

    count_transactions = length(trans.session.transactions)
    @info "Channel Nr: $nr_channel is open. Transactions open: $count_transactions"
    for strings in data
        ins_str = "insert \n" * join(strings, "\n")
        res = insert(trans, ins_str)
        if isempty(res)
            @async put!(error_channel, ins_str)
            @info ins_str
        end
    end
    commit(trans)
    put!(count_channel, nr_channel)
end

function insert_data(session::AbstractCoreSession,
            data::Vector{String},
            intervall::Int64 = 100,
            sub_intervall::Int64 = 50)

    ranges = collect(Iterators.partition(Iterators.partition(data, sub_intervall), intervall))

    length_ranges = length(ranges)
    @info "Length ranges $(length_ranges)"

    count_channel = Channel(10)
    error_channel = Channel(10)
    err_queries = String[]

    for i in 1:length_ranges
        tmp_client = CoreClient("127.0.0.1", 1729)
        tmp_session = CoreSession(tmp_client, database_name , Proto.Session_Type.DATA, request_timeout=Inf, error_time = 20)
        trans_in = transaction(tmp_session, Proto.Transaction_Type.WRITE)
        @async transaction_insert(trans_in, ranges[i], count_channel, i, error_channel)
    end

    @async begin
        while isopen(error_channel)
            yield()
            if isready(error_channel)
                push!(err_queries, take!(error_channel))
            end
        end
    end

    for _ in 1:length_ranges
       tmp = take!(count_channel)
       @info "Task $tmp is ready"
    end

    return err_queries
end

rest = uniprot(client, database_name)


text:
k_col = fill(k₁, ne(s))
k_col[cyl] .= k₂
k = diagm(k_col)

funcs[:mcopy] = Dict(:operator => I(ne(sd)), :type => MatrixFunc())
funcs[:k] = Dict(:operator => k, :type => MatrixFunc())
funcs[:e2p] = Dict(:operator => e2p * I(nv(sd)), :type => MatrixFunc())
funcs[:e2t] = Dict(:operator => e2t * I(nv(sd)), :type => MatrixFunc())
funcs[:t2e] = Dict(:operator => (1/e2t) * I(nv(sd)), :type => MatrixFunc())
funcs[:half] = Dict(:operator => 0.5 * I(nv(sd)), :type => MatrixFunc())
funcs[:third] = Dict(:operator => (1/3) * I(ne(sd)), :type => MatrixFunc())
funcs[:R₀] = Dict(:operator => R₀ * I(nv(sd)), :type => MatrixFunc())
funcs[:kᵥ] = Dict(:operator => kᵥ * I(ne(sd)), :type => MatrixFunc())
funcs[:kᵨ] = Dict(:operator => kᵨ * I(ne(sd)), :type => MatrixFunc())
funcs[:dneg₁] = Dict(:operator => -1 * I(ne(sd)), :type => MatrixFunc())
funcs[:mult₀] = Dict(:operator => (x′,x,y) -> (x′ .= x .* y), :type => InPlaceFunc())
funcs[:plus] = Dict(:operator => (x′,x,y) -> (x′ .= x .+ y), :type => InPlaceFunc())
funcs[:sum₀] = Dict(:operator => (x′,x,y) -> (x′ .= x .+ y), :type => InPlaceFunc())
funcs[:sum₁] = Dict(:operator => (x′,x,y) -> (x′ .= x .+ y), :type => InPlaceFunc())
funcs[:sum₂] = Dict(:operator => (x′,x,y) -> (x′ .= x .+ y), :type => InPlaceFunc())
funcs[:sum₀̃] = Dict(:operator => (x′,x,y) -> (x′ .= x .+ y), :type => InPlaceFunc())
funcs[:sum₁̃] = Dict(:operator => (x′,x,y) -> (x′ .= x .+ y), :type => InPlaceFunc())
funcs[:sum₂̃] = Dict(:operator => (x′,x,y) -> (x′ .= x .+ y), :type => InPlaceFunc())
funcs[:sub₀] = Dict(:operator => (x′,x,y) -> (x′ .= x .- y), :type => InPlaceFunc())
funcs[:div₀] = Dict(:operator => (x′,x,y) -> (x′ .= x ./ y), :type => InPlaceFunc())
funcs[:div₁] = Dict(:operator => (x′,x,y) -> (x′ .= x ./ y), :type => InPlaceFunc())
funcs[:avg₀₁] = Dict(:operator => avg_mat(Val{(0,1)}, s), :type => MatrixFunc())
funcs[:neg₁] = Dict(:operator => -1 * I(ne(sd)), :type => MatrixFunc())
funcs[:neg₀] = Dict(:operator => -1 * I(nv(sd)), :type => MatrixFunc())
funcs[:avg] = Dict(:operator => -1 * I(nv(sd)), :type => MatrixFunc())

funcs[:∂ₗ₀₊] = Dict(:operator => (x′,x) -> (x′ .= x; x′[∂ₗ₀₊] .= 0), :type => InPlaceFunc())
funcs[:∂ₜₐ] = Dict(:operator => (x′,x) -> (x′ .= x; x′[cyl_inner] .= 0), :type => InPlaceFunc())
funcs[:∂ₜ] = Dict(:operator => (x′,x) -> (x′ .= x; x′[∂ₑ₀₊] .= 0; x′[∂ₒ₀] .= 0), :type => InPlaceFunc())
funcs[:∂ₚ] = Dict(:operator => (x′,x) -> (x′ .= x; x′[∂ₑ₀₊] .= 0; x′[cyl_inner] .= 0), :type => InPlaceFunc())
funcs[:∂ᵥ] = Dict(:operator => (x′,x) -> (x′ .= x; x′[cyl_edge] .= 0; x′[∂ₑ₁₊] .= 0), :type => InPlaceFunc())

const ramp_up = 0.1 # seconds
# Final temperatures should be
temperature = t2e * 172.89

# Pre-allocate boundary condition values
const ∂ₗ₀_z = zeros(length(∂ₗ₀))
const ∂ᵣ₀_z = zeros(length(∂ᵣ₀))
const ∂ₗ₀_t = fill(temperature / ramp_up, length(∂ₗ₀))
const ∂ᵣ₀_t = fill(-(temperature / ramp_up), length(∂ᵣ₀))
ramp_up_temp(x′, x, t) = begin
  x′ .= x
  if t < ramp_up
    setindex!(x′, ∂ₗ₀_t, ∂ₗ₀)
    setindex!(x′, ∂ᵣ₀_t, ∂ᵣ₀)
  else
    setindex!(x′, ∂ₗ₀_z, ∂ₗ₀)
    setindex!(x′, ∂ᵣ₀_z, ∂ᵣ₀)
  end
end

const ∂₁₊_z = zeros(length(∂₁₊))
const ∂ₗ₁₊_z = zeros(length(∂ₗ₁₊))
const ∂ᵣ₁₊_z = zeros(length(∂ᵣ₁₊))
const ∂ₑ₀₊_z = zeros(length(∂ₑ₀₊))
funcs[:∂₀ₛ] = Dict(:operator => ramp_up_temp, :type => TDInPlaceFunc())
funcs[:∂₁ₛ] = Dict(:operator => (x′,x) -> (x′ .= x; setindex!(x′, ∂₁₊_z ,  ∂₁₊)), :type => InPlaceFunc())
funcs[:∂₁ₑ] = Dict(:operator => (x′,x) -> (x′ .= x; setindex!(x′, ∂ₗ₁₊_z, ∂ₗ₁₊); setindex!(x′, ∂ᵣ₁₊_z, ∂ᵣ₁₊)), :type => InPlaceFunc())
funcs[:mask₀ₗ] = Dict(:operator => (x′,x) -> (x′ .= x; setindex!(x′, ∂ₑ₀₊_z, ∂ₑ₀₊)), :type => InPlaceFunc())

wedge_cache = init_wedge_ops(sd)
v2comp = comp_support(sd);
cache_mat = Dict(:t2c => tri2comp(s, v2comp), :e2c => edge2comp(s, v2comp), :cross => changes(sd, v2comp),
                 :α_cache => zeros(ntriangles(sd)*3), :β_cache => zeros(ntriangles(sd)*3))

funcs[:∧₁₀′] = Dict(:operator => (x′, α, β) -> (cp_2_1!(x′, β, α, cache_mat)), :type => InPlaceFunc())
funcs[:∧₁₁′] = Dict(:operator => (x′, α, β) -> (pd_wedge!(x′, Val{(1,1)},s, α, β; wedge_cache...)), :type => InPlaceFunc())

const hodge_LU = lu(⋆(Val{1}, sd))
funcs[:⋆₁⁻¹] = Dict(:operator => (x′,x) -> (x′ .= -1 * (hodge_LU \ x)), :type => InPlaceFunc())

max_stars_repo_name: AlgebraicJulia/DECAPODES-Benchmarks

max_issues_repo_name: AlgebraicJulia/DECAPODES-Benchmarks

max_forks_repo_name: AlgebraicJulia/DECAPODES-Benchmarks


text:
# This file is a part of SimilaritySearch.jl

using SearchModels, Random
using StatsBase
import SearchModels: combine, mutate
export OptimizeParameters, optimize!, MinRecall, ParetoRecall, ParetoRadius

abstract type ErrorFunction end
@with_kw struct MinRecall <: ErrorFunction
    minrecall = 0.9
end

struct ParetoRecall <: ErrorFunction end
struct ParetoRadius <: ErrorFunction end

function runconfig end
function setconfig! end

function create_error_function(index::AbstractSearchContext, gold, knnlist::Vector{KnnResult}, queries, ksearch, verbose)
    n = length(index)
    nt = Threads.nthreads()
    vmin = Vector{Float64}(undef, nt)
    vmax = Vector{Float64}(undef, nt)
    vacc = Vector{Float64}(undef, nt)
    covradius = Vector{Float64}(undef, length(knnlist))
    pools = getpools(index)
    R = [Set{Int32}() for _ in knnlist]

    function lossfun(conf)
        vmin .= typemax(eltype(vmin))
        vmax .= typemin(eltype(vmax))
        vacc .= 0.0

        searchtime = @elapsed begin
            Threads.@threads for i in 1:length(queries)
                knnlist[i] = reuse!(knnlist[i], ksearch)
                _, v = runconfig(conf, index, queries[i], knnlist[i], pools)
                ti = Threads.threadid()
                vmin[ti] = min(v, vmin[ti])
                vmax[ti] = max(v, vmax[ti])
                vacc[ti] += v
            end
        end

        v = minimum(vmin), sum(vacc)/length(knnlist), maximum(vmax)

        for i in eachindex(knnlist)
            res = knnlist[i]
            covradius[i] = length(res) == 0 ? typemax(Float32) : maximum(res)
        end

        rmin, rmax = extrema(covradius)
        ravg = mean(covradius)

        recall = if gold !== nothing
            for (i, res) in enumerate(knnlist)
                empty!(R[i])
                union!(R[i], res.id)
            end

            macrorecall(gold, R)
        else
            nothing
        end

        verbose && println(stderr, "error_function> config: $conf, searchtime: $searchtime, recall: $recall, length: $(length(index))")
        (visited=v, radius=(rmin, ravg, rmax), recall=recall, searchtime=searchtime/length(knnlist))
    end
end

_kfun(x) = 1.0 - 1.0 / (1.0 + x)

"""
    optimize!(
        index::AbstractSearchContext,
        kind::ErrorFunction=ParetoRecall(),
        space::AbstractSolutionSpace=optimization_space(index);
        queries=nothing,
        ksearch=10,
        numqueries=64,
        initialpopulation=8,
        verbose=false,
        params=SearchParams(; maxpopulation=8, bsize=4, mutbsize=8, crossbsize=2, tol=-1.0, maxiters=8, verbose)
    )

Tries to configure the `index` to achieve the specified performance (`kind`). The optimization procedure is an stochastic search over the configuration space yielded by `kind` and `queries`.

# Arguments
- `index`: the index to be optimized
- `kind`: The kind of optimization to apply, it can be `ParetoRecall()`, `ParetoRadius()` or `MinRecall(r)` where `r` is the expected recall (0-1, 1 being the best quality but at cost of the search time)

# Keyword arguments

- `queries`: the set of queries to be used to measure performances, a validation set. It can be an `AbstractDatabase` or nothing.
- `queries_ksearch`: the number of neighbors to retrieve for `queries`
- `queries_size`: if `queries===nothing` then a sample of the already indexed database is used, `queries_size` is the size of the sample.
- `initialpopulation`: the initial sample for the optimization procedure
- `space`: defines the search space
- `params`: the parameters of the solver, see [`search_models` function from `SearchModels.jl`](https://github.com/sadit/SearchModels.jl) package for more information.
- `verbose`: controls if the procedure is verbose or not
"""
function optimize!(
            index::AbstractSearchContext,
            kind::ErrorFunction=ParetoRecall(),
            space::AbstractSolutionSpace=optimization_space(index);
            queries=nothing,
            ksearch=10,
            numqueries=64,
            initialpopulation=16,
            verbose=false,
            params=SearchParams(; maxpopulation=16, bsize=4, mutbsize=16, crossbsize=8, tol=-1.0, maxiters=16, verbose)
    )

    if queries === nothing
        sample = rand(1:length(index), numqueries) |> unique
        queries = SubDatabase(index.db, sample)
    end

    knnlist = [KnnResult(ksearch) for i in eachindex(queries)]
    gold = nothing
    if kind isa ParetoRecall || kind isa MinRecall
        db = @view index.db[1:length(index)]
        seq = ExhaustiveSearch(index.dist, db)
        searchbatch(seq, queries, knnlist; parallel=true)
        gold = [Set(res.id) for res in knnlist]
    end

    M = Ref(0.0)
    R = Ref(0.0)
    function inspect_pop(space, params, population)
        if M[] == 0.0
            for (c, p) in population
                M[] = max(p.visited[end], M[])
                R[] = max(p.radius[end], R[])
            end
        end
    end

    errorfun = create_error_function(index, gold, knnlist, queries, ksearch, verbose)

    function geterr(p)
        cost = p.visited[2] / M[]
        if kind isa ParetoRecall
            cost^2 + (1.0 - p.recall)^2
        elseif kind isa MinRecall
            p.recall < kind.minrecall ? 3.0 - 2 * p.recall : cost
        else
            _kfun(cost) + _kfun(p.radius[2] / R[])
        end
    end

    bestlist = search_models(
        errorfun,
        space,
        initialpopulation,
        params;
        inspect_population=inspect_pop,
        geterr=geterr)

    config, perf = bestlist[1]
    verbose && println(stderr, "== finished opt. $(typeof(index)): search-params: $(params), opt-config: $config, perf: $perf, kind=$(kind), length=$(length(index))")
    setconfig!(config, index, perf)
    bestlist
end

max_stars_repo_name: sadit/NNS

max_issues_repo_name: sadit/NNS

max_forks_repo_name: sadit/NNS


text:
module TestDoctest

import ShowGraphviz
using Documenter: doctest
using Test

function test_doctest()
    doctest(ShowGraphviz; manual = false)
end

end  # module