text:
module TrueSkill

global const BETA = 1.0::Float64
global const MU = 0.0::Float64
global const SIGMA = (BETA * 6)::Float64
global const GAMMA = (BETA * 0.03)::Float64
global const P_DRAW = 0.0::Float64
global const EPSILON = 1e-6::Float64
global const ITERATIONS = 30::Int64
global const sqrt2 = sqrt(2)
global const sqrt2pi = sqrt(2*pi)
global const PI = 1/(SIGMA^2)
global const TAU = MU*PI
global const minInt64 = (-9223372036854775808)::Int64
global const maxInt64 = ( 9223372036854775807)::Int64

function erfc(x::Float64)
    """Complementary error function (thanks to http://bit.ly/zOLqbc)"""
    z = abs(x)
    t = 1.0 / (1.0 + z / 2.0)
    r = begin
        a = -0.82215223 + t * 0.17087277
        b =  1.48851587 + t * a
        c = -1.13520398 + t * b
        d =  0.27886807 + t * c
        e = -0.18628806 + t * d
        f =  0.09678418 + t * e
        g =  0.37409196 + t * f
        h =  1.00002368 + t * g
        t * exp(-z * z - 1.26551223 + t * h)
        end
    if x < 0
        r = 2.0 - r
    end
    return r
end
function erfcinv(y::Float64)
    """The inverse function of erfc."""
    if y >= 2
        return -Inf
    elseif y < 0
        throw(DomainError(y, "argument must be nonnegative"))
    elseif y == 0
        return Inf
    end
    zero_point = y < 1
    if ! zero_point
        y = 2 - y
    end
    t = sqrt(-2 * log(y / 2.0))
    x = -0.70711 * ((2.30753 + t * 0.27061) / (1.0 + t * (0.99229 + t * 0.04481)) - t)
    for _ in 0:2
        err = erfc(x) - y
        x += err / (1.12837916709551257 * exp(-(x^2)) - x * err)
    end
    if zero_point
        r = x
    else
        r = -x
    end
    return r
end
function tau_pi(mu::Float64, sigma::Float64)
    if sigma > 0.
        _pi = sigma^-2
        _tau = _pi * mu
    elseif (sigma + 1e-5) < 0.0
        error("sigma should be greater than 0")
    else
        _pi = Inf
        _tau = Inf
    end
    return _tau, _pi
end
function mu_sigma(_tau::Float64, _pi::Float64)
    if _pi > 0.0
        sigma = sqrt(1/_pi)
        mu = _tau / _pi
    elseif (_pi + 1e-5) < 0.
        error("Precision should be greater than 0")
    else
        sigma = Inf
        mu = 0.0
    end
    return mu, sigma
end

struct Gaussian
    mu::Float64
    sigma::Float64
    function Gaussian(mu::Float64=MU, sigma::Float64=SIGMA)
        if sigma>=0.0
            return new(mu, sigma)
        else
            throw(error("Require: (sigma >= 0.0)"))
        end
    end
    function Gaussian(;mu::Float64=MU, sigma::Float64=SIGMA)
        return new(mu, sigma)
    end
end

global const N01 = Gaussian(0.0, 1.0)
global const Ninf = Gaussian(0.0, Inf)
global const N00 = Gaussian(0.0, 0.0)

Base.show(io::IO, g::Gaussian) = print("Gaussian(mu=", round(g.mu,digits=3), ", sigma=", round(g.sigma,digits=3), ")")
function _pi_(N::Gaussian)
    if N.sigma > 0.
        return N.sigma^-2
    else
        return Inf
    end
end
function _tau_(N::Gaussian)
    if N.sigma > 0.
        return N.mu * (N.sigma^-2)
    else
        return Inf
    end
end
function cdf(N::Gaussian, x::Float64)
    z = -(x - N.mu) / (N.sigma * sqrt2)
    return (0.5 * erfc(z))::Float64
end
function pdf(N::Gaussian, x::Float64)
    normalizer = (sqrt(2*pi) * N.sigma)^-1
    functional = exp( -((x - N.mu)^2) / (2*N.sigma ^2) )
    return (normalizer * functional)::Float64
end
function ppf(N::Gaussian, p::Float64)
    return N.mu - N.sigma * sqrt2  * erfcinv(2 * p)
end
function approx(;N::Gaussian, margin::Float64, tie::Bool)
    approx(N, margin, tie)
end
function approx(N::Gaussian, margin::Float64, tie::Bool)
    #The range is [alpha, beta]
    if !tie
        _alpha = (margin-N.mu)/N.sigma
        v = pdf(N01,-_alpha) / cdf(N01,-_alpha)
        w = v * (v + (-_alpha))
    else
        _alpha = (-margin-N.mu)/N.sigma
        _beta  = (margin-N.mu)/N.sigma
        v = (pdf(N01,_alpha)-pdf(N01,_beta))/(cdf(N01,_beta)-cdf(N01,_alpha))
        u = (_alpha*pdf(N01,_alpha)-_beta*pdf(N01,_beta))/(cdf(N01,_beta)-cdf(N01,_alpha))
        w =  - ( u - v^2 )
    end
    mu = N.mu + N.sigma * v
    sigma = N.sigma*sqrt(1-w)
    return Gaussian(mu, sigma)
end
function delta(N::Gaussian, M::Gaussian)
    return abs(N.mu - M.mu) , abs(N.sigma - M.sigma)
end
function exclude(N::Gaussian,M::Gaussian)
    return Gaussian(N.mu - M.mu, sqrt(N.sigma^2 - M.sigma^2) )
end
function Base.:+(N::Gaussian, M::Gaussian)
    mu = N.mu + M.mu
    sigma = sqrt(N.sigma^2 + M.sigma^2)
    return Gaussian(mu, sigma)
end
function Base.:-(N::Gaussian, M::Gaussian)
    mu = N.mu - M.mu
    sigma = sqrt(N.sigma^2 + M.sigma^2)
    return Gaussian(mu, sigma)
end
function Base.:*(N::Gaussian, M::Gaussian)
    _pi = _pi_(N) + _pi_(M)
    _tau = _tau_(N) + _tau_(M)
    mu, sigma = mu_sigma(_tau, _pi)
    return Gaussian(mu, sigma)
end
function Base.:/(N::Gaussian, M::Gaussian)
    _pi = _pi_(N) - _pi_(M)
    _tau = _tau_(N) - _tau_(M)
    mu, sigma = mu_sigma(_tau, _pi)
    return Gaussian(mu, sigma)
end
function Base.isapprox(N::Gaussian, M::Gaussian, atol::Real=0)
    return (abs(N.mu - M.mu) < atol) & (abs(N.sigma - M.sigma) < atol)
end
function forget(N::Gaussian, gamma::Float64, t::Int64)
    return Gaussian(N.mu, sqrt(N.sigma^2 + t*gamma^2))
end
function compute_margin(p_draw::Float64, sd::Float64)
    _N = Gaussian(0.0, sd )
    res = abs(ppf(_N, 0.5-p_draw/2))
    return res
end
struct Player
    prior::Gaussian
    beta::Float64
    gamma::Float64
    draw::Gaussian
    function Player(prior::Gaussian=Gaussian(MU,SIGMA), beta::Float64=BETA, gamma::Float64=GAMMA, draw::Gaussian=Ninf)
        new(prior, beta, gamma, draw)
    end
    function Player(;prior::Gaussian=Gaussian(MU,SIGMA), beta::Float64=BETA, gamma::Float64=GAMMA, draw::Gaussian=Ninf)
        Player(prior, beta, gamma, draw)
    end
end
Base.show(io::IO, r::Player) = print("Player(Gaussian(mu=", round(r.prior.mu,digits=3),", sigma=", round(r.prior.sigma,digits=3), "), beta=", round(r.beta, digits=3), ", gamma=" , round(r.gamma, digits=3),")")
function performance(R::Player)
    return Gaussian(R.prior.mu, sqrt(R.prior.sigma^2 + R.beta^2))
end

function likelihood(R::Player, likelihood_performance::Gaussian)
    return Gaussian(likelihood_performance.mu, sqrt(likelihood_performance.sigma^2 + R.beta^2))
end

mutable struct team_messages
    prior::Gaussian
    likelihood_lose::Gaussian
    likelihood_win::Gaussian
    likelihood_draw::Gaussian
end
function p(tm::team_messages)
    return tm.prior*tm.likelihood_lose*tm.likelihood_win*tm.likelihood_draw
end
function posterior_win(tm::team_messages)
    return tm.prior*tm.likelihood_lose*tm.likelihood_draw
end
function posterior_lose(tm::team_messages)
    return tm.prior*tm.likelihood_win*tm.likelihood_draw
end
function posterior_draw(tm::team_messages)
    return tm.prior*tm.likelihood_win*tm.likelihood_lose
end
function likelihood(tm::team_messages)
    return tm.likelihood_win*tm.likelihood_lose*tm.likelihood_draw
end
mutable struct draw_messages
    prior::Gaussian
    prior_team::Gaussian
    likelihood_lose::Gaussian
    likelihood_win::Gaussian
end
function p(um::draw_messages)
    return um.prior_team*um.likelihood_lose*um.likelihood_win
end
function posterior_win(um::draw_messages)
    return um.prior_team*um.likelihood_lose
end
function posterior_lose(um::draw_messages)
    return um.prior_team*um.likelihood_win
end
function likelihood(um::draw_messages)
    return um.likelihood_win*um.likelihood_lose
end
mutable struct diff_messages
    prior::Gaussian
    likelihood::Gaussian
end
function p(dm::diff_messages)
    return dm.prior*dm.likelihood
end

function Base.max(tuple1::Tuple{Float64,Float64}, tuple2::Tuple{Float64,Float64})
    return max(tuple1[1],tuple2[1]), max(tuple1[2],tuple2[2])
end
function Base.:>(tuple::Tuple{Float64,Float64}, threshold::Float64)
    return (tuple[1] > threshold) | (tuple[2] > threshold)
end

mutable struct Game
    teams::Vector{Vector{Player}}
    result::Vector{Float64}
    p_draw::Float64
    likelihoods::Vector{Vector{Gaussian}}
    evidence::Float64
    function Game(teams::Vector{Vector{Player}}, result::Vector{Float64}=Float64[],p_draw::Float64=0.0)
        ((0.0 > p_draw) | (1.0 <= p_draw)) &&  throw(error("0.0 <= Draw probability < 1.0"))
        (length(result)>0) && (length(teams)!= length(result)) && throw(error("(length(result)>0) & (length(teams)!= length(result))"))
        (p_draw == 0.0) && (length(result)>0) && (length(unique(result)) != length(result)) && throw(error("(p_draw == 0.0) && (length(result)>0) && (length(unique(result)) != length(result))"))

        _g = new(teams,result,p_draw,[],0.0)
        likelihoods(_g)
        return _g
    end
    function Game(teams::Vector{Vector{Player}}, result::Vector{Float64}=Float64[];p_draw::Float64=0.0)
        Game(teams, result, p_draw)
    end
    function Game(teams::Vector{Vector{Player}}; result::Vector{Float64}=Float64[],p_draw::Float64=0.0)
        Game(teams, result, p_draw)
    end
end
Base.length(G::Game) = length(G.teams)
function size(G::Game)
    return [length(team) for team in g.teams]
end
function performance(team::Vector{Player})
    p = N00
    for r in team
        p += performance(r)
    end
    return p
end
#TODO: reducir perf(game)
function performance(G::Game,i::Int64)
    res = N00
    for r in G.teams[i]
        res += performance(r)
    end
    return res
end
function draw_performance(G::Game,i::Int64)
    res = N00
    for r in G.teams[i]
        res += r.draw.sigma < Inf ? approx(r.draw,0.,false) : Ninf
    end
    return res
end
function likelihood_teams(g::Game)
    r = g.result == Float64[] ? [Float64(i) for i in length(g.teams):-1:1] : g.result
    o = sortperm(r, rev=true)
    t = [team_messages(performance(g,o[e]), Ninf, Ninf, Ninf) for e in 1:length(g)]
    d = [diff_messages(t[e].prior - t[e+1].prior, Ninf) for e in 1:length(g)-1]
    tie = [r[o[e]]==r[o[e+1]] for e in 1:length(d)]
    margin = [ g.p_draw==0.0 ?  0.0 :
               compute_margin(g.p_draw, sqrt( sum([a.beta^2 for a in g.teams[o[e]]]) + sum([a.beta^2 for a in g.teams[o[e+1]]]) ))
               for e in 1:length(d)]
    g.evidence = 1
    for e in 1:length(d)
        g.evidence *= !tie[e] ? 1-cdf(d[e].prior, margin[e]) : cdf(d[e].prior, margin[e])-cdf(d[e].prior, -margin[e])
    end
    step = (Inf, Inf)::Tuple{Float64,Float64}; iter = 0::Int64
    while (step > 1e-6) & (iter < 10)
        step = (0., 0.)
        for e in 1:length(d)-1#e=1
            d[e].prior = posterior_win(t[e]) - posterior_lose(t[e+1])
            d[e].likelihood = approx(d[e].prior,margin[e],tie[e])/d[e].prior
            likelihood_lose = posterior_win(t[e]) - d[e].likelihood
            step = max(step,delta(t[e+1].likelihood_lose,likelihood_lose))
            t[e+1].likelihood_lose = likelihood_lose
        end
        for e in length(d):-1:2
            d[e].prior = posterior_win(t[e]) - posterior_lose(t[e+1])
            d[e].likelihood = approx(d[e].prior,margin[e],tie[e])/d[e].prior
            likelihood_win = (posterior_lose(t[e+1]) + d[e].likelihood)
            step = max(step,delta(t[e].likelihood_win,likelihood_win))
            t[e].likelihood_win = likelihood_win
        end
        iter += 1
    end
    if length(d)==1
        d[1].prior = posterior_win(t[1]) - posterior_lose(t[2])
        d[1].likelihood = approx(d[1].prior,margin[1],tie[1])/d[1].prior
    end
    t[1].likelihood_win = (posterior_lose(t[2]) + d[1].likelihood)
    t[end].likelihood_lose = (posterior_win(t[end-1]) - d[end].likelihood)

    return [ likelihood(t[o[e]]) for e in 1:length(t)]
end
function likelihoods(g::Game)
    m_t_ft = likelihood_teams(g)
    g.likelihoods = [[
        likelihood(g.teams[e][i],
                   m_t_ft[e] - exclude(performance(g,e),performance(g.teams[e][i])))
        for i in 1:length(g.teams[e])] for e in 1:length(g)]
    return g.likelihoods
end
function posteriors(g::Game)
    return [[ g.likelihoods[e][i] * g.teams[e][i].prior for i in 1:length(g.teams[e])] for e in 1:length(g)]
end
mutable struct Skill
    forward::Gaussian
    backward::Gaussian
    likelihood::Gaussian
    elapsed::Int64
    online::Gaussian
    function Skill(forward::Gaussian=Ninf, backward::Gaussian=Ninf, likelihood::Gaussian=Ninf, elapsed::Int64=0)
        return new(forward, backward, likelihood, elapsed)
    end
    function Skill(;forward::Gaussian=Ninf, backward::Gaussian=Ninf, likelihood::Gaussian=Ninf, elapsed::Int64=0)
        return new(forward, backward, likelihood, elapsed)
    end
end
mutable struct Agent
    player::Player
    message::Gaussian
    last_time::Int64
end
function receive(agent::Agent, elapsed::Int64)
    if agent.message != Ninf
        res = forget(agent.message, agent.player.gamma, elapsed)
    else
       res = agent.player.prior
    end
    return res
end
function clean(agents::Dict{String,Agent}, last_time::Bool=false)
    for a in keys(agents)
        agents[a].message = Ninf
        if last_time agents[a].last_time = minInt64 end
    end
end
mutable struct Item
    agent::String
    likelihood::Gaussian
end
mutable struct Team
    items::Vector{Item}
    output::Float64
end
mutable struct Event
    teams::Vector{Team}
    evidence::Float64
end
function outputs(event::Event)
    return [ team.output for team in event.teams]
end
function get_composition(events::Vector{Event})
    return [[[i.agent for i in t.items] for t in e.teams] for e in events]
end
function get_results(events::Vector{Event})
    return [ [t.output for t in e.teams ] for e in events]
end
function compute_elapsed(last_time::Int64, actual_time::Int64)
    return last_time == minInt64 ? 0 : ( last_time == maxInt64 ? 1 : (actual_time - last_time))
end
mutable struct Batch
    time::Int64
    events::Vector{Event}
    skills::Dict{String,Skill}
    agents::Dict{String,Agent}
    p_draw::Float64
    function Batch(composition::Vector{Vector{Vector{String}}}, results::Vector{Vector{Float64}}=Vector{Vector{Float64}}() ,time::Int64=0, agents::Dict{String,Agent}=Dict{String,Agent}(), p_draw::Float64 = 0.0)
        (length(results)>0) & (length(composition)!= length(results)) && throw(error("(length(results)>0) & (length(composition)!= length(results))"))

        this_agents = Set(vcat((composition...)...))
        elapsed = Dict([ (a, compute_elapsed(agents[a].last_time, time) ) for a in this_agents  ])
        skills = Dict([ (a, Skill(receive(agents[a],elapsed[a]) ,Ninf ,Ninf , elapsed[a])) for a in this_agents  ])
        events = [Event([Team([Item(composition[e][t][a], Ninf) for a in 1:length(composition[e][t]) ]
                              , length(results) > 0 ? results[e][t] : Float64(length(composition[e])-t)  ) for t in 1:length(composition[e]) ]
                        ,0.0) for e in 1:length(composition) ]

        b = new(time, events , skills, agents, p_draw)

        iteration(b)
        return b
    end
    function Batch(;composition::Vector{Vector{Vector{String}}}, results::Vector{Vector{Float64}}
                 ,time::Int64=0 , agents::Dict{String,Agent}=Dict{String,Agent}(), p_draw::Float64 = 0.0)
        Batch(composition, results, time, agents, p_draw)
    end
end

Base.show(io::IO, b::Batch) = print("Batch(time=", b.time, ", events=", b.events, ")")
Base.length(b::Batch) = length(b.events)

function add_events(b::Batch, composition::Vector{Vector{Vector{String}}}, results::Vector{Vector{Float64}}=Vector{Vector{Float64}}())
    this_agents = Set(vcat((composition...)...))
    for a in this_agents#a="c"
        elapsed = compute_elapsed(b.agents[a].last_time , b.time )
        if !haskey(b.skills,a)
            b.skills[a] = Skill(receive(b.agents[a],elapsed) ,Ninf ,Ninf , elapsed)
        else
            b.skills[a].elapsed = elapsed
            b.skills[a].forward = receive(b.agents[a],elapsed)
        end
    end
    from = length(b)+1
    for e in 1:length(composition)
        event = Event([Team([Item(composition[e][t][a], Ninf) for a in 1:length(composition[e][t]) ]
                              , length(results) > 0 ? results[e][t] : Float64(length(composition[e])-t) ) for t in 1:length(composition[e]) ] , 0.0)
        push!(b.events, event)
    end
    iteration(b, from)
end
function posterior(b::Batch, agent::String)#agent="a_b"
    return b.skills[agent].likelihood*b.skills[agent].backward*b.skills[agent].forward
end
function posteriors(b::Batch)
    res = Dict{String,Gaussian}()
    for (a, s) in b.skills
        res[a] = posterior(b,a)
    end
    return res
end
function within_prior(b::Batch, item::Item, online = false)
    r = b.agents[item.agent].player
    if online
        return Player(b.skills[item.agent].forward,r.beta,r.gamma)
    else
        wp = posterior(b,item.agent)/item.likelihood
        return Player(wp,r.beta,r.gamma)
    end
end
function within_priors(b::Batch, event::Int64; online = false)#event=1
    return [ [within_prior(b,item,online) for item in team.items ] for team in b.events[event].teams ]
end
function iteration(b::Batch, from::Int64 = 1)
    for e in from:length(b)#e=1

        g = Game(within_priors(b, e), outputs(b.events[e]), b.p_draw)

        for (t, team) in enumerate(b.events[e].teams)#(t,team) = (2, b.events[e].teams[2])
            for (i, item) in enumerate(team.items)#(i, item) = (2, team.items[2])
                b.skills[item.agent].likelihood = (b.skills[item.agent].likelihood / item.likelihood) * g.likelihoods[t][i]
                item.likelihood = g.likelihoods[t][i]
            end
        end

        b.events[e].evidence = g.evidence
    end
end
function log_evidence2(b::Batch; online::Bool = false, agents::Vector{String} = Vector{String}())
    if isempty(agents)
        if online
            return sum([log(Game(within_priors(b, e, online=true), outputs(b.events[e]), b.p_draw).evidence) for e in 1:length(b)])
        else
            return sum([log(event.evidence) for event in b.events])
        end
    else
        filter = [!isdisjoint(vcat((comp...)...),agents) for comp in get_composition(b.events)]
        if online
            return sum([log(Game(within_priors(b, e, online=true), outputs(b.events[e]), b.p_draw).evidence) for e in 1:length(b) if filter[e] ])
        else
            return sum([log(b.events[e].evidence) for e in 1:length(b.events) if filter[e]])
        end
    end
end
function convergence(b::Batch, epsilon::Float64=1e-6, iterations::Int64 = 20)
    iter = 0::Int64
    step = (Inf, Inf)
    while (step > epsilon) & (iter < iterations)
        old = copy(posteriors(b))
        iteration(b)
        step = diff(old, posteriors(b))
        iter += 1
    end
    return iter
end
function forward_prior_out(b::Batch, agent::String)
    return b.skills[agent].forward * b.skills[agent].likelihood
end
function backward_prior_out(b::Batch, agent::String)
    N = b.skills[agent].likelihood*b.skills[agent].backward
    return forget(N, b.agents[agent].player.gamma, b.skills[agent].elapsed)
end
function new_backward_info(b::Batch)
    for a in keys(b.skills)
        b.skills[a].backward = b.agents[a].message
    end
    return iteration(b)
end
function new_forward_info(b::Batch)
    for a in keys(b.skills)
        b.skills[a].forward = receive(b.agents[a], b.skills[a].elapsed)
    end
    return iteration(b)
end
mutable struct History
    size::Int64
    batches::Vector{Batch}
    agents::Dict{String,Agent}
    time::Bool
    mu::Float64
    sigma::Float64
    beta::Float64
    gamma::Float64
    p_draw::Float64
    online::Bool
    function History(composition::Vector{Vector{Vector{String}}}, results::Vector{Vector{Float64}}=Vector{Vector{Float64}}(), times::Vector{Int64}=Int64[], priors::Dict{String,Player}=Dict{String,Player}(); mu::Float64=MU, sigma::Float64=SIGMA, beta::Float64=BETA, gamma::Float64=GAMMA, p_draw::Float64=P_DRAW, online::Bool=false)
        (length(results) > 0) & (length(composition) != length(results)) && throw(error("(length(times) > 0) & (length(composition) != length(results))"))
        (length(times) > 0) & (length(composition) != length(times)) && throw(error("length(times) > 0) & (length(composition) != length(times))"))

        agents = Dict([ (a, Agent(haskey(priors, a) ? priors[a] : Player(Gaussian(mu, sigma), beta, gamma), Ninf, minInt64)) for a in Set(vcat((composition...)...)) ])
        h = new(length(composition), Vector{Batch}(), agents, length(times)>0, mu, sigma, beta, gamma, p_draw, online)
        trueskill(h, composition, results, times, online)
        return h
    end
    function History(;composition::Vector{Vector{Vector{String}}},results::Vector{Vector{Float64}}=Vector{Vector{Float64}}(),times::Vector{Int64}=Int64[],priors::Dict{String,Player}=Dict{String,Player}(), mu::Float64=MU, sigma::Float64=SIGMA, beta::Float64=BETA, gamma::Float64=GAMMA, p_draw::Float64=P_DRAW, online::Bool=false)
        History(composition, results, times, priors, mu=mu, sigma=sigma, beta=beta, gamma=gamma, p_draw=p_draw, online=online)
    end
end

# History(composition::Vector{Vector{Vector{String}}}, results::Vector{Vector{Int64}}=Vector{Vector{Int64}}(), times::Vector{Int64}=Int64[], priors::Dict{String,Player}=Dict{String,Player}(); mu::Float64=MU, sigma::Float64=SIGMA, beta::Float64=BETA, gamma::Float64=GAMMA, p_draw::Float64=P_DRAW, online::Bool=false) = History(composition, convert(Vector{Vector{Float64}},results), times, priors, mu=mu, sigma=sigma, beta=beta, gamma=gamma, p_draw=p_draw, online=online)
#
# History(;composition::Vector{Vector{Vector{String}}}, results::Vector{Vector{Int64}}=Vector{Vector{Int64}}(), times::Vector{Int64}=Int64[], priors::Dict{String,Player}=Dict{String,Player}(), mu::Float64=MU, sigma::Float64=SIGMA, beta::Float64=BETA, gamma::Float64=GAMMA, p_draw::Float64=P_DRAW, online::Bool=false) = History(composition, convert(Vector{Vector{Float64}},results), times, priors, mu=mu, sigma=sigma, beta=beta, gamma=gamma, p_draw=p_draw, online=online)

Base.length(h::History) = h.size
Base.show(io::IO, h::History) = print("History(Events=", h.size
                                     ,", Batches=", length(h.batches)
                                    ,", Agents=", length(h.agents), ")")
function trueskill(h::History, composition::Vector{Vector{Vector{String}}},results::Vector{Vector{Float64}}, times::Vector{Int64}, online::Bool)
    o = length(times)>0 ? sortperm(times) : [i for i in 1:length(composition)]
    i = 1::Int64
    last = 0.0
    while i <= length(h)
        j, t = i, length(times) == 0 ? i : times[o[i]]
        while ((length(times)>0) & (j < length(h)) && (times[o[j+1]] == t)) j += 1 end
        if length(results)>0
            b = Batch(composition[o[i:j]],results[o[i:j]], t, h.agents, h.p_draw)
        else
            b = Batch(composition[o[i:j]], results , t, h.agents, h.p_draw)
        end
        push!(h.batches,b)
        if online
            new = round(100*(i/length(h)))
            if new != last
                print("\r",new,"%")
                last = new
            end
            for a in keys(b.skills)
                b.skills[a].online = b.skills[a].forward
            end
            convergence(h)
        end
        for a in keys(b.skills)
            h.agents[a].last_time = length(times) == 0 ? maxInt64 : t
            h.agents[a].message = forward_prior_out(b,a)
        end
        i = j + 1
    end
    if online println("\r100.0%") end
end
function diff(old::Dict{String,Gaussian}, new::Dict{String,Gaussian})
    step = (0., 0.)
    for a in keys(old)
        step = max(step, delta(old[a],new[a]))
    end
    return step
end
function iteration(h::History)
    step = (0., 0.)

    clean(h.agents)
    for j in length(h.batches)-1:-1:1# j = 2
        for a in keys(h.batches[j+1].skills)
            h.agents[a].message = backward_prior_out(h.batches[j+1],a)
        end
        old = copy(posteriors(h.batches[j]))
        new_backward_info(h.batches[j])
        step = max(step, diff(old, posteriors(h.batches[j])))
    end

    clean(h.agents)
    for j in 2:length(h.batches)# j = 3
        for a in keys(h.batches[j-1].skills)
            h.agents[a].message = forward_prior_out(h.batches[j-1],a)
        end
        old = copy(posteriors(h.batches[j]))
        new_forward_info(h.batches[j])
        step = max(step, diff(old, posteriors(h.batches[j])))
    end

    if (length(h.batches) == 1)
        old = copy(posteriors(h.batches[1]))
        iteration(h.batches[1])
        step = max(step, diff(old, posteriors(h.batches[1])))
    end

    return step
end
function convergence(h::History; epsilon::Float64=EPSILON, iterations::Int64=ITERATIONS, verbose = true)
    step = (Inf, Inf)::Tuple{Float64,Float64}
    iter = 1::Int64
    while (step > epsilon) & (iter <= iterations)
        verbose && print("Iteration = ", iter)
        step = iteration(h)
        iter += 1
        verbose && println(", step = ", step)
    end
    verbose && println("End")
    return step, iter
end
# function posteriors(h::History)
#     res = Vector{Vector{Vector{Gaussian}}}()
#     for b in h.batches
#         for e in b.events
#             ps = Vector{Vector{Gaussian}}()
#             for a in get_composition(e)
#
#                 t_p = (b.time, posterior(b,a))
#                 res[a] = haskey(res, a) ? push!(res[a],t_p) : [t_p]
#         end
#     end
# end
function learning_curves(h::History)
    res = Dict{String,Array{Tuple{Int64,Gaussian}}}()
    for b in h.batches
        for a in keys(b.skills)
            t_p = (b.time, posterior(b,a))
            res[a] = haskey(res, a) ? push!(res[a],t_p) : [t_p]
        end
    end
    return res
end
function log_evidence(h::History; online::Bool = false, agents::Vector{String} = Vector{String}() )
    return sum([log_evidence2(b, online=online, agents = agents) for b in h.batches])
end
function add_events(h::History,composition::Vector{Vector{Vector{String}}};results::Vector{Vector{Float64}}=Vector{Vector{Float64}}(),times::Vector{Int64}=Int64[],priors::Dict{String,Player}=Dict{String,Player}())
    add_events(h, composition, results, times, priors)
end
function add_events(h::History,composition::Vector{Vector{Vector{String}}},results::Vector{Vector{Float64}}=Vector{Vector{Float64}}(),times::Vector{Int64}=Int64[],priors::Dict{String,Player}=Dict{String,Player}())

    (length(times)>0) & !h.time && throw(error("length(times)>0 but !h.time"))
    (length(times)==0) & h.time && throw(error("length(times)==0 but h.time"))
    (length(results) > 0) & (length(composition) != length(results)) && throw(error("(length(results) > 0) & (length(composition) != length(results))"))
    (length(times) > 0) & (length(composition) != length(times)) && throw(error("length(times) > 0) & (length(composition) != length(times))"))

    this_agents = Set(vcat((composition...)...))
    for a in this_agents
        if !haskey(h.agents,a)
            h.agents[a] = Agent(haskey(priors, a) ? priors[a] : Player(Gaussian(h.mu, h.sigma), h.beta, h.gamma), Ninf, minInt64)
        end
    end

    clean(h.agents,true)
    n = length(composition)
    o = length(times)>0 ? sortperm(times) : [i for i in 1:length(composition)]
    i = 1::Int64; k = 1::Int64
    while i <= n
        j, t = i, length(times) == 0 ? i : times[o[i]]
        while ((length(times)>0) & (j < n) && (times[o[j+1]] == t)) j += 1 end
        while (!h.time & (h.size > k) ) || (h.time && (length(h.batches) >= k) && (h.batches[k].time < t))
            b = h.batches[k]
            if (k>1) new_forward_info(b) end
            for a in intersect(keys(b.skills), this_agents)#a ="a"
                b.skills[a].elapsed = compute_elapsed(b.agents[a].last_time , b.time )
                h.agents[a].last_time = length(times) == 0 ? maxInt64 : b.time
                h.agents[a].message = forward_prior_out(b,a)
            end
            k += 1
        end
        if (h.time && (length(h.batches) >= k) && (h.batches[k].time == t))
            b = h.batches[k]
            if length(results)>0
                add_events(b, composition[o[i:j]],results[o[i:j]])
            else
                add_events(b, composition[o[i:j]])
            end
        else
            if !h.time k = k + 1 end
            if length(results)>0
                b =  Batch(composition[o[i:j]],results[o[i:j]], t, h.agents, h.p_draw)
            else
                b =  Batch(composition[o[i:j]], results , t, h.agents, h.p_draw)
            end
            insert!(h.batches,  k , b)
            if h.time k = k + 1 end
        end
        for a in keys(b.skills)#a="a"
            h.agents[a].last_time = length(times) == 0 ? maxInt64 : t
            h.agents[a].message = forward_prior_out(b,a)
        end
        i = j + 1
    end
    while h.time && (length(h.batches) >= k)
        b = h.batches[k]
        new_forward_info(b)
        for a in intersect(keys(b.skills), this_agents)#a ="a"
            b.skills[a].elapsed = compute_elapsed(b.agents[a].last_time , b.time )
            h.agents[a].last_time = length(times) == 0 ? maxInt64 : b.time
            h.agents[a].message = forward_prior_out(b,a)
        end
        k += 1
    end
    h.size = h.size + n
    iteration(h)
end

end # module


text:
function init(randomGenerator::Any, seed::Any, μ = 1, σ = 0)
    """
      the First (Functional) Function

      Inputs:

      Mean μ [1st-Order  moment-generating function]  (fixed, #TODO: can it be a Variable?)
      STD σ  [2nd-Order moment-generating function] (fixed , #TODO: can it be Variable?)

      RNG seed
      Function randomGenerator(seed)

      outputs:

      # TODO: Does the seed `remain` the same? can it be(come) changeable afterwards?
      """
    μ = μ
    σ = σ
    seed = seed

    randomGenerator = randomGenerator(seed)

    return randomGenerator
end

function init(randomGenerator, seed, μ = 1, σ = 0)
    """
    inputs:

    Mean μ [1st-Order  moment-generating function]  (fixed, #TODO: can it be a Variable?)
    STD σ  [2nd-Order moment-generating function] (fixed , #TODO: can it be Variable?)

    seed: RNG seed
    Function randomGenerator(seed)

    # TODO: Does the seed `remain` the same? can it be(come) changeable afterwards?
    """

    model.μ = μ
    model.σ = σ
    model.seed = seed

    model.randomGenerator = randomGenerator(seed)

    return randomGenerator()
end

function init(randomGenerator, seed, μ = 1, σ = 0)
    """
    inputs:

    Mean μ [1st moment] (fixed, #TODO: can it be a Variable?)
    STD σ   [2nd moment] (fixed , #TODO: can it be a Variable?)
    RNG seed
    Function randomGenerator(seed)

    # TODO:
    Q1. Does the seed `remain` the same? OR
    Q2.can it be(come) changeable afterwards?
    """

    model.μ = μ
    model.σ = σ
    model.seed = seed

    model.randomGenerator = randomGenerator(seed)

    return randomGenerator

end

#---  arrays = convertions -

function array2rational(T::Any)
    return rationalize(vec(Float64.(T)))

end

function array2irrational(T::Any)
    return irr(vec(Float64.(T)))
end

function array2vector(T::Any)
    return
    Float64.(vcat(T[1], vec(T[2])))
end

#working
function randomArray(min=0,max=1,n=100,y = exp(1 / -x))
    a=[];b=[];X=[]
    for i in enumerate(n)
        x = [i in rand(min,max)];X=[X x]
        c=[i in rand(1, y)]; a = [a c]
        d=[j in rand(1, y)]; b = [b d]
    end
    return X,a,b
end

ext: jl

max_stars_repo_name: adamwillisMastery/DeepLearner

max_issues_repo_name: adamwillisMastery/DeepLearner

max_forks_repo_name: adamwillisXanax/DeepLearner


text:
@testset "Util" begin

    str = "This is a 寬度不一的 Char array."
    str = DS.padding(str, 50)
    @test textwidth(str) == 50

    str = ""
    str = DS.padding(str, 50)
    @test textwidth(str) == 50

    str = "apple\norange"
    str = DS.padding_vertical(str, 3)
    @test str == "apple\norange\n"

end


text:
"""
    sme(RR::Vector{<:Real}, SBP::Vector{<:Real}; thresholdRR::Float64 = 4.0, thresholdSBP::Float64 = 1.0, seqLen::Int = 3, delay::Int = 1, minCor::Float64 = 0.8)

Calculate the Baroreflex Sensitivity (BRS) for a serious of RR intervals and respective systolic blood pressure using the Sequence Method (SME).

# Keywords:

* thresholdRR: Threshold for change in RR interval to qualify for a valid sequence (literature: 4ms)
* thresholdSBP: Threshold for change in SBP to qualify for a valid sequence (literature: 1mmHg)
* seqLen: Minimum length of a valid sequence (literature: 3)
* delay: Delay between RR and SBP (literature: 1)
* minCor: The minimal correlation between RR and SBP in a sequence to qualify as a valid sequence

# Return:

Returns a SME struct for which a plotting recipe is provided, so that found sequences can be inspected by calling `plot()` when using Plots.jl

"""
function sme(RR::Vector{<:Real}, SBP::Vector{<:Real}; thresholdRR::Float64 = 5.0, thresholdSBP::Float64 = 1.0, seqLen::Int = 3, delay::Int = 1, minCor::Float64 = 0.85)
    # apply delay
    if delay > 0
        RR = RR[1+delay:end]
        SBP = SBP[1:end-delay]
    else
        SBP = SBP[1+delay:end]
        RR = RR[1:end-delay]
    end

    # find and store sequences
    difRR = diff(RR)
    difSBP = diff(SBP)
    start = Vector{Int}() # stores all sequence starts
    n = Vector{Int}() # stores respective sequence length
    direction = Vector{Bool}() # stores the direction of the sequence, true = Up, false = Down

    #temporary variables
    tempStart = -1
    tempN = -1
    up = false # indicator if current sequence is ascending
    down = false # indicator if current sequence is descending
    for i in 1:length(difRR)
        if (difRR[i] > thresholdRR) & (difSBP[i] > thresholdSBP) & up
            tempN += 1
        elseif (difRR[i] < -thresholdRR) & (difSBP[i] < -thresholdSBP) & down
            tempN += 1
        else # means that a sequence gets interrupted
            if tempN >= seqLen
                push!(start, tempStart)
                push!(n, tempN)
                push!(direction, up)
            end
            # initialize new starting sequence
            if ((difRR[i] > thresholdRR) & (difSBP[i] > thresholdSBP)) || ((difRR[i] < -thresholdRR) & (difSBP[i] < -thresholdSBP))
            tempStart = i
            tempN = 2
            up = difRR[i] > 0
            down = !up
            else
                up = false
                down = false
                tempStart = -1
                tempN = -1
            end
        end
    end
    # push last sequence if present
    if tempN >= seqLen
        push!(start, tempStart)
        push!(n, tempN)
        push!(direction, up)
    end

    # is sequence valid (cor > threshold)
    valid = Vector{Bool}(undef, length(start))
    for i in 1:length(start)
        valid[i] = cor(RR[start[i]:start[i]+n[i]-1], SBP[start[i]:start[i]+n[i]-1]) > minCor
    end
    # select valid ones
    start = start[valid]
    n = n[valid]
    direction = direction[valid]

    # calculate BRS as slope of a sequences, with robust regression
    slopes = Vector{Float64}(undef, length(start))
    #slopes2 = Vector{Float64}(undef, length(start))
    for i in 1:length(start)
            # slopes[i] = mean(diag(SBP[start[i]:start[i]+n[i]-1] \ RR[start[i]:start[i]+n[i]-1]))
            slopes[i] = theilSenEstimator(SBP[start[i]:start[i]+n[i]-1], RR[start[i]:start[i]+n[i]-1])[1]
        # df = DataFrame(x = SBP[start[i]:start[i]+n[i]-1], y = RR[start[i]:start[i]+n[i]-1])
        # slopes2[i] = coef(lm(@formula(y ~ x), df))[2]
    end
    # some more return values
    nUp = sum(direction)
    nDown = length(direction) - nUp
    slopeUp = slopes[direction]
    slopeDown = slopes[direction .== false]
    # return (mean(slopes), mean(slopes2))
    result = SME(value = mean(slopes), sBRSup= mean(slopeUp), sBRSdown= mean(slopeDown), nUp = nUp, nDown = nDown, start = start, n = n, direction = direction, SBP = SBP, RR = RR)
    return result
    # return (mean(slopes), mean(slopeUp), mean(slopeDown), nUp, nDown, start, n, direction)
    # return (median(slopes), median(slopeUp), median(slopeDown), nUp, nDown)
end

"""
    rmssdr(RR::Vector{<:Real}, SBP::Vector{<:Real})

Returns the RMSSD ratio := RMSSD(RR) / RMSSD(SBP)
"""
function rmssdr(RR::Vector{<:Real}, SBP::Vector{<:Real})
    rmssdr = rms(RR) / rms(SBP)
end

"""
Struct that stores all information related to the sequence method. The main result is stored in 'value'. It can be plotted for visual inspection.
"""
@with_kw mutable struct SME
    value::Real = 0.0
    sBRSup::Real = 0.0
    sBRSdown::Real = 0.0
    delay::Int = 0
    nUp::Real = 0
    nDown::Real = 0
    start::Vector{<:Int} = Int[]
    n::Vector{<:Int} = Int[]
    direction::Vector{<:Bool} = Bool[]
    RR = Real[]
    SBP = Real[]
end

@recipe function f(res::SME)

    title := "sBRS = $(round(res.value, digits = 2))[ms/mmHg], #sequences = $(length(res.n))"
    layout := (2,1)

    @series begin
        subplot := 1
        yguide := "SBP [mmHg]"
        label := ""
        color := :grey

        res.SBP
    end

    @series begin
        subplot := 2
        yguide := "RR [ms]"
        label := ""
        color := :grey

        res.RR
    end
    starts = res.start
    for (i, start) in enumerate(starts)
        @series begin
            subplot := 1
            yguide := "SBP [mmHg]"
            label := ""
            seriescolor := res.direction[i] ? :red : :green
            linewidth := 2

            start:(start + res.n[i]-1), res.SBP[start:(start + res.n[i]-1)]
            # start+1:(start+1 + res.n[i]-1), res.SBP[start+1:(start+1 + res.n[i]-1)]
        end

        @series begin
            subplot := 2
            yguide := "RR [ms]"
            label := ""
            seriescolor := res.direction[i] ? :red : :green
            linewidth := 2

            start:(start + res.n[i]-1), res.RR[start:(start + res.n[i]-1)]
        end
    end
end


text:
module WIAGwebData

using MySQL
using Infiltrator
using DataFrames

dbwiag = nothing

function setDBWIAG(pwd = missing, host = "127.0.0.1", user = "wiag", db = "wiag")
    global dbwiag
    if ismissing(pwd)
        println("Passwort für User ", user)
        pwd = readline()
    end
    dbwiag = DBInterface.connect(MySQL.Connection, host, user, pwd, db = db)
end

"""
    updatenamevariant()

obsolete: see namelookup
"""
function updatenamevariant(fieldsrc::AbstractString, tablename::AbstractString)::Int
    dbwiag = DBInterface.connect(MySQL.Connection, "localhost", "wiag", "Wogen&Wellen", db="wiag");

    DBInterface.execute(dbwiag, "DELETE FROM " * tablename);

    # Do not use the same DB connection with an open cursor to insert data
    dfsrc = DBInterface.execute(dbwiag,
                                "SELECT wiagid, " * fieldsrc
                                * " FROM person") |> DataFrame;

    tblid = 1
    for row in eachrow(dfsrc)
        id, fns = row
        # println("id: ", id)
        # println(fns)

        if ismissing(fns) || fns == "" continue end
        for nv in split(fns, r",|;")
            insertstmt = ("INSERT INTO " * tablename * " VALUES ("
                          * string(tblid) * ","
                          * string(id) * ","
                          * "'" * strip(nv) * "')")
            # println(insertstmt)
            DBInterface.execute(dbwiag, insertstmt)
            tblid += 1;
        end
    end
    return tblid;
end

"""
    fillera(tblera::AbstractString, tblperson::AbstractString, tbloffice::AbstractString, colnameid::AbstractString, datereference=false)::Int

Compute earliest and latest date for each person, identified by `colnameid` and `colnameidinoffice`.

Take fields `date_hist_first` and `date_hist_last` into account if `datereference` is set to `true`.
"""
function fillera(tblera::AbstractString,
                 tblperson::AbstractString,
                 tbloffice::AbstractString,
                 colnameid = "id",
                 colnameidinoffice = "id_person";
                 datereference = false)::Int
    global dbwiag
    if isnothing(dbwiag)
        error("There is no valid database connection. Use `setDBWIAG'.")
    end

    DBInterface.execute(dbwiag, "DELETE FROM " * tblera);

    if datereference
        sqlselect = "SELECT " * colnameid * " as idperson, " *
            " date_birth, date_death, date_hist_first, date_hist_last " *
            " FROM " * tblperson;
    else
        sqlselect = "SELECT " * colnameid * " as idperson, date_birth, date_death " * " FROM " * tblperson;
    end

    dfperson = DBInterface.execute(dbwiag, sqlselect) |> DataFrame;

    rgxyear = r"[1-9][0-9][0-9]+";
    rgxcentury = r"([1-9][0-9])?\. [Jahrh|Jhd]"
    tblid = 0;

    # officestmt = DBInterface.prepare(dbwiag, "SELECT date_start, date_end FROM office"
    #                                  * " WHERE wiagid_person = ?")

    sqlselect = "SELECT " * colnameidinoffice * ", date_start, date_end " * " FROM " * tbloffice
    dfoffice = DBInterface.execute(dbwiag, sqlselect) |> DataFrame

    function parsemaybe(s)::Union{Missing, Int}
        r = missing
        if !ismissing(s)
            rgm = match(rgxyear, s)
            if !isnothing(rgm)
                r = parse(Int, rgm.match)
            else
                rgm = match(rgxcentury, s)
                if !isnothing(rgm) && !isnothing(rgm[1])
                    r = parse(Int, rgm[1]) * 100 - 50
                end
            end
        end
        return r
    end

    csqlvalues = String[]
    for row in eachrow(dfperson)
        erastart = Inf
        eraend = -Inf

        idperson, datebirth, datedeath = row[[:idperson, :date_birth, :date_death]]

        vcand = parsemaybe(datebirth)
        if !ismissing(vcand) erastart = vcand end

        vcand = parsemaybe(datedeath)
        if !ismissing(vcand) eraend = vcand end

        if datereference
            datehistfirst, datehistlast = row[[:date_hist_first, :date_hist_last]]
            vcand = parsemaybe(datehistfirst)
            if !ismissing(vcand) && vcand < erastart
                erastart = vcand
            end
            vcand = parsemaybe(datehistlast)
            if !ismissing(vcand) && vcand > eraend
                eraend = vcand
            end
        end

        # println(wiagid, " ", typeof(dfoffice[:wiagid_person]))
        ixperson = dfoffice[:, colnameidinoffice] .== string(idperson)

        dfofficeperson = dfoffice[ixperson, :];
        for oc in eachrow(dfofficeperson)
            datestart = oc[:date_start]
            dateend = oc[:date_end]

            vcand = parsemaybe(datestart)
            if !ismissing(vcand) && vcand < erastart
                erastart = vcand
            end

            vcand = parsemaybe(dateend)
            if !ismissing(vcand) && vcand > eraend
                eraend = vcand
            end

        end

        if erastart == Inf && eraend != -Inf
            erastart = eraend
        elseif erastart != Inf && eraend == -Inf
            eraend = erastart
        end

        sqlerastart = erastart == Inf ? "NULL" : "'" * string(erastart) * "'"
        sqleraend = eraend == -Inf ? "NULL" : "'" * string(eraend) * "'"

        # if !ismissing(erastartdb) && !(typeof(erastartdb) == Int)
        #     println("start: ", erastartdb)
        # end
        # if !ismissing(eraenddb) && !(typeof(eraenddb) == Int)
        #     println("end: ", eraenddb)
        # end

        push!(csqlvalues, "('" * string(idperson) * "', " * sqlerastart * ", " * sqleraend * ")")
        tblid += 1
    end

    sqlvalues = join(csqlvalues, ", ")
    DBInterface.execute(dbwiag, "INSERT INTO " * tblera * " VALUES " * sqlvalues)

    # println(sqlvalues);

    return tblid
end

"""
    fillofficedate(tblofficedate::AbstractString, tbloffice::AbstractString, colnameid::AbstractString)::Int

Extract dates as an integer values.
"""
function fillofficedate(tblofficedate::AbstractString,
                        tbloffice::AbstractString,
                        colnameid::AbstractString = "id")::Int

    global dbwiag
    if isnothing(dbwiag)
        setDBWIAG()
    end

    rgxyear = r"[1-9][0-9][0-9]+";
    rgxcentury = r"([1-9][0-9])?\. [Jahrh|Jhd]"

    function parsemaybe(s)::String
        r = "NULL"
        if !ismissing(s)
            rgm = match(rgxyear, s)
            if !isnothing(rgm)
                r = parse(Int, rgm.match) |> string
            else
                rgm = match(rgxcentury, s)
                if !isnothing(rgm) && !isnothing(rgm[1])
                    r = parse(Int, rgm[1]) * 100 - 50 |> string
                end
            end
        end
        return r
    end

    sql = "SELECT " * colnameid * ", date_start, date_end FROM " * tbloffice
    dfoffice = DBInterface.execute(dbwiag, sql) |> DataFrame;

    csqlvalues = String[]
    tblid = 0;
    for row in eachrow(dfoffice)
        id, date_start, date_end = row

        numdate_start = parsemaybe(date_start)
        numdate_end = parsemaybe(date_end)

        push!(csqlvalues, "(" * id * ", " * numdate_start * ", " * numdate_end * ")")

        tblid += 1
        # if tblid > 25 break end
    end

    DBInterface.execute(dbwiag, "DELETE FROM " * tblofficedate);

    sqlvalues = join(csqlvalues, ", ")
    DBInterface.execute(dbwiag, "INSERT INTO " * tblofficedate * " VALUES " * sqlvalues)

    return tblid
end

"""
    fillnamelookup(tablename::AbstractString)::Int

Fill `tablename` with combinations of givenname and familyname and their variants.
"""
function fillnamelookup(tbllookup::AbstractString,
                        tblperson::AbstractString,
                        colnameid::AbstractString = "id")::Int
    msg = 200
    if isnothing(dbwiag)
        error("There is no valid database connection. Use `setDBWIAG'.")
    end

    DBInterface.execute(dbwiag, "DELETE FROM " * tbllookup);

    dfperson = DBInterface.execute(dbwiag,
                                   "SELECT " * colnameid * " as id_person, " *
                                   "givenname, prefix_name, familyname, givenname_variant, familyname_variant " *
                                   "FROM " * tblperson * " person") |> DataFrame;

    # SQL
    # INSERT INTO dsttable VALUES (NULL, 'id_person1', 'givenname1', 'prefix_name1', 'familyname1'),
    # ('NULL', 'id_person2', 'givenname2', 'prefix_name2', 'familyname2');
    #
    # structure
    # gn[:] prefix fn|fnv
    # gn[1] prefix fn|fnv
    # gnv[:] prefix fn|fnv
    # gnv[1] prefix fn|fnv

    # In the web application choose a version with or without prefix.

    imsg = 0
    csqlvalues = String[]
    appendtosqlrow(row) = append!(csqlvalues, row)
    for row in eachrow(dfperson)
        idperson = row[:id_person]
        gn = row[:givenname]
        prefix = row[:prefix_name]
        fn = row[:familyname]
        gnv = row[:givenname_variant]
        fnv = row[:familyname_variant]

        fillnamelookupgn(idperson, gn, prefix, fn, fnv) |> appendtosqlrow

        if !ismissing(gnv) && gnv != ""
            # sets of givennames
            cgnv = split(gnv, r", *")
            for gnve in cgnv
                fillnamelookupgn(idperson, gnve, prefix, fn, fnv) |> appendtosqlrow
            end
        end
        imsg += 1

        if imsg % msg == 0
            println("write row: ", imsg)
        end
    end
    sqlvalues = join(csqlvalues, ", ")

    irowout = length(csqlvalues)

    DBInterface.execute(dbwiag, "INSERT INTO " * tbllookup * " VALUES " * sqlvalues)

    return irowout

end

"""
    sqlstring(s::AbstractString)::AbstractString

remove labels in data fields ("Taufname: Karl") and escape apostrophes
"""
function sqlstring(s::AbstractString)::AbstractString
    poslabel = findfirst(':', s)
    if !isnothing(poslabel)
        s = strip(s[poslabel + 1:end])
    end
    s = replace(s, "'" => "''")
    s = "'" * s * "'"
    return s
end

function fillnamelookupgn(id_person, gn, prefix, fn, fnv)
    csql = String[]

    isnull(s) = ismissing(s) || s == "NULL"
    function pushcsql(gni, fni)
        sgni = isnull(gni) ? "NULL" : sqlstring(gni)
        sfni = isnull(fni) ? "NULL" : sqlstring(fni)
        prefixi = isnull(prefix) ? "NULL" : sqlstring(prefix)
        if !ismissing(id_person)
            id_person_sql = sqlstring(id_person)
            values = "(" * "NULL, " * id_person_sql * ", " * sgni * ", " * prefixi * ", " * sfni * ")"
            push!(csql, values)
        else
            @warn "Missing ID for ", gni
        end
    end

    pushcsql(gn, fn)
    cgn = split(gn);
    # more than one givenname -> write a version with the first givenname only
    if length(cgn) > 1
        pushcsql(cgn[1], fn)
    end

    # familyname variants
    if !ismissing(fnv) && strip(fnv) != ""
        cfnv = split(fnv, r", *")
        for fnve in cfnv
            pushcsql(gn, fnve)
            # more than one givenname
            if length(cgn) > 1
                pushcsql(cgn[1], fnve)
            end
        end
    end

    return csql
end

end

max_stars_repo_name: WIAG-ADW-GOE/WIAGweb

max_issues_repo_name: WIAG-ADW-GOE/WIAGweb

max_forks_repo_name: WIAG-ADW-GOE/WIAGweb


text:
using IcanProj
using JLD2

function check_waterbalance(res_vic::VicRes)

    # Precipitation

    data_prec = get_variable(res_vic, "prec")

    # Evapotranspiration

    data_evap = get_variable(res_vic, "evap")

    # data_trans_veg = get_variable(res_vic, "transp_veg")

    # Sublimation snow

    data_sub_snow = get_variable(res_vic, "sub_snow")

    data_sub_canop = get_variable(res_vic, "sub_canop")

    # Runoff

    data_baseflow = get_variable(res_vic, "baseflow")

    data_runoff = get_variable(res_vic, "runoff")

    # Check waterbalance

    inputs = sum(data_prec)

    outputs = sum(data_evap) + sum(data_sub_snow) + sum(data_sub_canop) + sum(data_baseflow) + sum(data_runoff)

    # Print results

    println("Inputs: $(round(inputs)) mm")
    println("Outputs: $(round(outputs)) mm")
    println("Difference: $(round(inputs)-round(outputs)) mm")

end

# Loop over all simulations

opt = get_options()

watersheds = opt["stat_sel"]

resolutions = opt["resolutions"]

for watershed in watersheds

    for resolution in resolutions

        file = joinpath(opt["eval_folder"], watershed, "res_$(resolution).jld2")

        @load file res_vic

        println()

        println("Watershed: $(watershed)")

        check_waterbalance(res_vic)

        println()

    end

end


text:
function get_free_state_parameters(obj::AbstractGradientConverter)
    return jcall(obj, "getFreeStateParameters", jint, ())
end

function get_parameters(obj::NumericalGradientConverter, arg0::FieldSpacecraftState, arg1::ForceModel)
    return jcall(obj, "getParameters", Vector{Gradient}, (FieldSpacecraftState, ForceModel), arg0, arg1)
end

function get_state(obj::NumericalGradientConverter, arg0::ForceModel)
    return jcall(obj, "getState", FieldSpacecraftState, (ForceModel,), arg0)
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

abstract AbstractElement

type Element{E<:AbstractElement}
    id :: Int
    connectivity :: Vector{Int}
    integration_points :: Vector{IP}
    fields :: Dict{AbstractString, Field}
    properties :: E
end

function Element{E<:AbstractElement}(::Type{E}, id::Int64, connectivity::Vector{Int64})
    return Element{E}(id, connectivity, [], Dict(), E())
end

function Element{E<:AbstractElement}(::Type{E}, connectivity::Vector{Int64})
    return Element{E}(-1, connectivity, [], Dict(), E())
end

function getindex(element::Element, field_name::AbstractString)
    return element.fields[field_name]
end

function setindex!(element::Element, data::Field, field_name)
    element.fields[field_name] = data
end

function setindex!(element::Element, data::Function, field_name)
    if method_exists(data, Tuple{Element, Vector, Float64})
        # create enclosure to pass element as argument
        function wrapper_(ip, time)
            return data(element, ip, time)
        end
        field = Field(wrapper_)
    else
        field = Field(data)
    end
    element.fields[field_name] = field
end

function setindex!(element::Element, data, field_name)
    element.fields[field_name] = Field(data)
end

function call(element::Element, field_name)
    return element[field_name]
end

function call(element::Element, field_name, time)
    return element[field_name](time)
end

function last(element::Element, field_name::AbstractString)
    return last(element[field_name])
end

function call(element::Element, ip, time::Float64=0.0)
    return get_basis(element, ip, time)
end

function call(element::Element, ip, time::Float64, ::Type{Val{:Jacobian}})
    X = element("geometry", time)
    dN = get_dbasis(element, ip, time)
    nbasis = length(element)
    if isa(X.data, Vector)
        J = sum([kron(dN[:,i], X[i]') for i=1:nbasis])
    else
        c = get_connectivity(element)
        J = sum([kron(dN[:,i], X[c[i]]') for i=1:nbasis])
    end
    return J
end

function call(element::Element, ip, time::Float64, ::Type{Val{:detJ}})
    J = element(ip, time, Val{:Jacobian})
    n, m = size(J)
    if n == m  # volume element
        return det(J)
    end
    JT = transpose(J)
    if size(JT, 2) == 1  # boundary of 2d problem, || ∂X/∂ξ ||
        return norm(JT)
    else # manifold on 3d problem, || ∂X/∂ξ₁ × ∂X/∂ξ₂ ||
        return norm(cross(JT[:,1], JT[:,2]))
    end
end

function call(element::Element, ip, time::Float64, ::Type{Val{:Grad}})
    J = element(ip, time, Val{:Jacobian})
    return inv(J)*get_dbasis(element, ip, time)
end

function call(element::Element, field_name::AbstractString, ip, time::Float64, ::Type{Val{:Grad}})
    return element(ip, time, Val{:Grad})*element[field_name](time)
end

function call(element::Element, field::Field, time)
    return field(time)
end

function call(element::Element, field::DCTI, time)
    return field.data
end

function call(element::Element, field_name::AbstractString, time)
    field = element[field_name]
    return element(field, time)
end

function call(element::Element, field_name::AbstractString, ip, time::Float64)
    field = element[field_name]
    return element(field, ip, time)
end

function call(element::Element, field::DCTI, ip, time::Float64)
    return field.data
end

function call(element::Element, field::DCTV, ip, time::Float64)
    return field(time).data
end

function call(element::Element, field::CVTV, ip, time::Float64)
    return field(ip, time)
end

function call(element::Element, field::Field, ip, time::Float64)
    field_ = field(time)
    basis = element(ip, time)
    n = length(element)
    if isa(field_.data, Vector)
        m = length(field_)
        if n != m
            error("Error when trying to interpolate field $field at coords $ip and time $time: element length is $n and field length is $m, f = Nᵢfᵢ makes no sense!")
        end
        return sum([field_[i]*basis[i] for i=1:n])
    else
        c = get_connectivity(element)
        return sum([field_[c[i]]*basis[i] for i=1:n])
    end
end

function size(element::Element, dim)
    return size(element)[dim]
end

""" Update element field based on a dictionary of nodal data and connectivity information.

Examples
--------
julia> data = Dict(1 => [0.0, 0.0], 2 => [1.0, 2.0])
julia> element = Seg2([1, 2])
julia> update!(element, "geometry", data)

As a result element now have time invariant (variable) vector field "geometry" with data ([0.0, 0.0], [1.0, 2.0]).

"""
function update!(element::Element, field_name, data::Dict)
    #element[field_name] = Field(data)
    element[field_name] = [data[i] for i in get_connectivity(element)]
end

function update!{K,V}(element::Element, field_name, data::Pair{Float64, Dict{K, V}})
    time, field_data = data
    element_data = V[field_data[i] for i in get_connectivity(element)]
    update!(element, field_name, time => element_data)
    #if haskey(element, field_name)
    #    update!(element[field_name], data)
    #else
    #    element[field_name] = Field(data)
    #end
end

function update!(element::Element, field_name::AbstractString, datas::Union{Real, Vector, Pair{Float64, Union{Float64, Real, Vector{Any}}}}...)
    for data in datas
        if haskey(element, field_name)
            update!(element[field_name], data)
        else
            if length(data) != length(element)
                update!(element, field_name, DCTI(data))
            else
                element[field_name] = data
            end
        end
    end
end

#=
function update!(element::Element, field_name, data::Pair...)
    for data in datas
        update!(element, field_name, data)
    end
end
=#

function update!(element::Element, field_name, data::Pair{Float64, Vector{Any}})
    if haskey(element, field_name)
        update!(element[field_name], data)
    else
        element[field_name] = data
    end
end

function update!(element::Element, field_name, data::Pair{Float64, Vector{Int64}})
    if haskey(element, field_name)
        update!(element[field_name], data)
    else
        element[field_name] = data
    end
end

function update!(element::Element, field_name, data::Pair{Float64, Vector{Float64}})
    if haskey(element, field_name)
        update!(element[field_name], data)
    else
        element[field_name] = data
    end
end

function update!(element::Element, field_name, data::Pair{Float64, Vector{Vector{Float64}}})
    if haskey(element, field_name)
        update!(element[field_name], data)
    else
        element[field_name] = data
    end
end

function update!(element::Element, field_name, data::Pair{Float64, Float64})
    if haskey(element, field_name)
        update!(element[field_name], data)
    else
        element[field_name] = data
    end
end

function update!(element::Element, field_name::AbstractString, data::Union{Float64, Vector})
    if haskey(element, field_name)
        update!(element[field_name], data)
    else
        if length(data) != length(element)
            update!(element, field_name, DCTI(data))
        else
            element[field_name] = data
        end
    end
end

function update!(element::Element, datas::Pair...)
    for (field_name, data) in datas
        if haskey(element, field_name)
            update!(element[field_name], data)
        else
            element[field_name] = data
        end
    end
end

function update!(element::Element, field_name, data::Function)
    element[field_name] = data
end

function update!(element::Element, field_name, field::Field)
    element[field_name] = field
end

function update!(elements::Vector, field_name, data)
    for element in elements
        update!(element, field_name, data)
    end
end

""" Check existence of field. """
function haskey(element::Element, field_name)
    haskey(element.fields, field_name)
end

function get_connectivity(element::Element)
    return element.connectivity
end

function get_integration_points(element::Element)
    # first time initialize default integration points
    if length(element.integration_points) == 0
        ips = get_integration_points(element.properties)
        element.integration_points = [IP(i, w, xi) for (i, (w, xi)) in enumerate(ips)]
    end
    return element.integration_points
end

""" This is a special case, temporarily change order
of integration scheme mainly for mass matrix.
"""
function get_integration_points(element::Element, change_order::Int)
    order = get_integration_order(element.properties)
    order += change_order
    ips = get_integration_points(element.properties, order)
    return [IP(i, w, xi) for (i, (w, xi)) in enumerate(ips)]
end

function get_gdofs(element::Element)
    return get_gdofs(element, 1)
end

""" Return dual basis transformation matrix Ae. """
function get_dualbasis(element::Element, time::Float64, order=1)
    nnodes = length(element)
    De = zeros(nnodes, nnodes)
    Me = zeros(nnodes, nnodes)
    for ip in get_integration_points(element, order)
        detJ = element(ip, time, Val{:detJ})
        w = ip.weight*detJ
        N = element(ip, time)
        De += w*diagm(vec(N))
        Me += w*N'*N
    end
    return De, Me, De*inv(Me)
end

""" Find inverse isoparametric mapping of element. """
function get_local_coordinates(element::Element, X::Vector, time::Float64; max_iterations=10, tolerance=1.0e-6)
    haskey(element, "geometry") || error("element geometry not defined, cannot calculate inverse isoparametric mapping")
    dim = size(element, 1)
    dim == length(X) || error("manifolds not supported.")
    xi = zeros(dim)
    dX = element("geometry", xi, time) - X
    for i=1:max_iterations
        J = element(xi, time, Val{:Jacobian})'
        xi -= J \ dX
        dX = element("geometry", xi, time) - X
        norm(dX) < tolerance && return xi
    end
    info("X = $X, dX = $dX, xi = $xi")
    error("Unable to find inverse isoparametric mapping for element $element for X = $X")
end

""" Test is X inside element. """
function inside{E}(element::Element{E}, X, time)
    xi = get_local_coordinates(element, X, time)
    return inside(E, xi)
end


text:
#! format: off
get_variable_multiplier(_, ::Type{<:PSY.Area}, ::AbstractAGCFormulation) = NaN
########################## ActivePowerVariable, Area ###########################

########################## SteadyStateFrequencyDeviation ##################################
get_variable_binary(::SteadyStateFrequencyDeviation, ::Type{<:PSY.Area}, ::AbstractAGCFormulation) = false

get_variable_binary(::ActivePowerVariable, ::Type{<:PSY.Area}, ::AbstractAGCFormulation) = false

########################## SmoothACE, AggregationTopology ###########################

get_variable_binary(::SmoothACE, ::Type{<:PSY.AggregationTopology}, ::AbstractAGCFormulation) = false

########################## DeltaActivePowerUpVariable, Area ###########################

get_variable_binary(::DeltaActivePowerUpVariable, ::Type{<:PSY.Area}, ::AbstractAGCFormulation) = false
get_variable_lower_bound(::DeltaActivePowerUpVariable, ::PSY.Area, ::AbstractAGCFormulation) = 0.0

########################## DeltaActivePowerDownVariable, Area ###########################

get_variable_binary(::DeltaActivePowerDownVariable, ::Type{<:PSY.Area}, ::AbstractAGCFormulation) = false
get_variable_lower_bound(::DeltaActivePowerDownVariable, ::PSY.Area, ::AbstractAGCFormulation) = 0.0

########################## AdditionalDeltaPowerUpVariable, Area ###########################

get_variable_binary(::AdditionalDeltaActivePowerUpVariable, ::Type{<:PSY.Area}, ::AbstractAGCFormulation) = false
get_variable_lower_bound(::AdditionalDeltaActivePowerUpVariable, ::PSY.Area, ::AbstractAGCFormulation) = 0.0

########################## AdditionalDeltaPowerDownVariable, Area ###########################

get_variable_binary(::AdditionalDeltaActivePowerDownVariable, ::Type{<:PSY.Area}, ::AbstractAGCFormulation) = false
get_variable_lower_bound(::AdditionalDeltaActivePowerDownVariable, ::PSY.Area, ::AbstractAGCFormulation) = 0.0

########################## AreaMismatchVariable, Area ###########################
get_variable_binary(::AreaMismatchVariable, ::Type{<:PSY.Area}, ::AbstractAGCFormulation) = false

########################## LiftVariable, Area ###########################
get_variable_binary(::LiftVariable, ::Type{<:PSY.Area}, ::AbstractAGCFormulation) = false
get_variable_lower_bound(::LiftVariable, ::PSY.Area, ::AbstractAGCFormulation) = 0.0

initial_condition_default(::AreaControlError, d::PSY.AGC, ::AbstractAGCFormulation) = PSY.get_initial_ace(d)
initial_condition_variable(::AreaControlError, d::PSY.AGC, ::AbstractAGCFormulation) = AreaMismatchVariable()

get_variable_multiplier(::SteadyStateFrequencyDeviation, d::PSY.AGC, ::AbstractAGCFormulation) = -10 * PSY.get_bias(d)

#! format: on

"""
Steady State deviation of the frequency
"""
function add_variables!(
    container::OptimizationContainer,
    ::Type{T},
) where {T <: SteadyStateFrequencyDeviation}
    time_steps = get_time_steps(container)
    variable = add_variable_container!(container, T(), PSY.Area, time_steps)
    for t in time_steps
        variable[t] = JuMP.@variable(container.JuMPmodel, base_name = "ΔF_{$(t)}")
    end
end

########################## Initial Condition ###########################

function _get_variable_initial_value(
    d::PSY.Component,
    key::ICKey,
    ::AbstractAGCFormulation,
    ::Nothing,
)
    return _get_ace_error(d, key)
end

function add_constraints!(
    container::OptimizationContainer,
    ::Type{T},
    ::Type{LiftVariable},
    areas::IS.FlattenIteratorWrapper{U},
    ::ServiceModel{PSY.AGC, V},
) where {T <: AbsoluteValueConstraint, U <: PSY.Area, V <: PIDSmoothACE}
    time_steps = get_time_steps(container)
    area_names = PSY.get_name.(areas)
    container_lb =
        add_constraints_container!(container, T(), U, area_names, time_steps, meta="lb")
    container_ub =
        add_constraints_container!(container, T(), U, area_names, time_steps, meta="ub")
    mismatch = get_variable(container, AreaMismatchVariable(), U)
    z = get_variable(container, LiftVariable(), U)
    jump_model = get_jump_model(container)

    for t in time_steps, a in area_names
        container_lb[a, t] = JuMP.@constraint(jump_model, mismatch[a, t] <= z[a, t])
        container_ub[a, t] = JuMP.@constraint(jump_model, -1 * mismatch[a, t] <= z[a, t])
    end
    return
end

"""
Expression for the power deviation given deviation in the frequency. This expression allows
updating the response of the frequency depending on commitment decisions
"""
function add_constraints!(
    container::OptimizationContainer,
    ::Type{T},
    ::Type{SteadyStateFrequencyDeviation},
    areas::IS.FlattenIteratorWrapper{U},
    ::ServiceModel{PSY.AGC, V},
    sys::PSY.System,
) where {T <: FrequencyResponseConstraint, U <: PSY.Area, V <: PIDSmoothACE}
    time_steps = get_time_steps(container)
    area_names = PSY.get_name.(areas)

    frequency_response = 0.0
    for area in PSY.get_components(PSY.Area, sys)
        frequency_response += PSY.get_load_response(area)
    end

    for g in PSY.get_components(PSY.RegulationDevice, sys, x -> PSY.get_available(x))
        d = PSY.get_droop(g)
        response = 1 / d
        frequency_response += response
    end

    IS.@assert_op frequency_response >= 0.0

    # This value is the one updated later in simulation based on the UC result
    inv_frequency_reponse = 1 / frequency_response

    area_balance = get_variable(container, ActivePowerVariable(), U)
    frequency = get_variable(container, SteadyStateFrequencyDeviation(), U)
    R_up = get_variable(container, DeltaActivePowerUpVariable(), U)
    R_dn = get_variable(container, DeltaActivePowerDownVariable(), U)
    R_up_emergency = get_variable(container, AdditionalDeltaActivePowerUpVariable(), U)
    R_dn_emergency = get_variable(container, AdditionalDeltaActivePowerUpVariable(), U)

    const_container = add_constraints_container!(container, T(), PSY.System, time_steps)

    for t in time_steps
        system_balance = sum(area_balance.data[:, t])
        for a in area_names
            JuMP.add_to_expression!(system_balance, R_up[a, t])
            JuMP.add_to_expression!(system_balance, -1 * R_dn[a, t])
            JuMP.add_to_expression!(system_balance, R_up_emergency[a, t])
            JuMP.add_to_expression!(system_balance, -1 * R_dn_emergency[a, t])
        end
        const_container[t] = JuMP.@constraint(
            container.JuMPmodel,
            frequency[t] == -inv_frequency_reponse * system_balance
        )
    end
    return
end

function add_constraints!(
    container::OptimizationContainer,
    ::Type{T},
    ::Type{SteadyStateFrequencyDeviation},
    areas::IS.FlattenIteratorWrapper{U},
    ::ServiceModel{PSY.AGC, V},
    sys::PSY.System,
) where {T <: SACEPIDAreaConstraint, U <: PSY.Area, V <: PIDSmoothACE}
    services = get_available_components(PSY.AGC, sys)
    time_steps = get_time_steps(container)
    area_names = [PSY.get_name(PSY.get_area(s)) for s in services]
    RAW_ACE = get_expression(container, RawACE(), U)
    SACE = get_variable(container, SmoothACE(), U)
    SACE_pid = add_constraints_container!(
        container,
        SACEPIDAreaConstraint(),
        U,
        area_names,
        time_steps,
    )

    jump_model = get_jump_model(container)
    for (ix, service) in enumerate(services)
        kp = PSY.get_K_p(service)
        ki = PSY.get_K_i(service)
        kd = PSY.get_K_d(service)
        Δt = convert(Dates.Second, container.resolution).value
        a = PSY.get_name(PSY.get_area(service))
        for t in time_steps
            if t == 1
                ACE_ini = get_initial_condition(container, AreaControlError(), PSY.AGC)[ix]
                ace_exp = get_value(ACE_ini) + kp * ((1 + Δt / (kp / ki)) * (RAW_ACE[a, t]))
                SACE_pid[a, t] = JuMP.@constraint(jump_model, SACE[a, t] == ace_exp)
                continue
            end
            SACE_pid[a, t] = JuMP.@constraint(
                jump_model,
                SACE[a, t] ==
                SACE[a, t - 1] +
                kp * (
                    (1 + Δt / (kp / ki) + (kd / kp) / Δt) * (RAW_ACE[a, t]) +
                    (-1 - 2 * (kd / kp) / Δt) * (RAW_ACE[a, t - 1])
                )
            )
        end
    end
    return
end

function add_constraints!(
    container::OptimizationContainer,
    ::Type{T},
    ::Type{SmoothACE},
    areas::IS.FlattenIteratorWrapper{U},
    ::ServiceModel{PSY.AGC, V},
    sys::PSY.System,
) where {T <: BalanceAuxConstraint, U <: PSY.Area, V <: PIDSmoothACE}
    time_steps = get_time_steps(container)
    area_names = PSY.get_name.(areas)
    aux_equation = add_constraints_container!(
        container,
        BalanceAuxConstraint(),
        PSY.System,
        area_names,
        time_steps,
    )
    area_mismatch = get_variable(container, AreaMismatchVariable(), PSY.Area)
    SACE = get_variable(container, SmoothACE(), PSY.Area)
    R_up = get_variable(container, DeltaActivePowerUpVariable(), PSY.Area)
    R_dn = get_variable(container, DeltaActivePowerDownVariable(), PSY.Area)
    R_up_emergency =
        get_variable(container, AdditionalDeltaActivePowerUpVariable(), PSY.Area)
    R_dn_emergency =
        get_variable(container, AdditionalDeltaActivePowerUpVariable(), PSY.Area)

    for t in time_steps, a in area_names
        aux_equation[a, t] = JuMP.@constraint(
            container.JuMPmodel,
            -1 * SACE[a, t] ==
            (R_up[a, t] - R_dn[a, t]) +
            (R_up_emergency[a, t] - R_dn_emergency[a, t]) +
            area_mismatch[a, t]
        )
    end
    return
end

function objective_function!(
    container::OptimizationContainer,
    areas::IS.FlattenIteratorWrapper{T},
    ::ServiceModel{<:PSY.AGC, U},
) where {T <: PSY.Area, U <: PIDSmoothACE}
    add_proportional_cost!(container, LiftVariable(), areas, U())
    return
end

# Defined here so we can dispatch on PIDSmoothACE
function add_feedforward_arguments!(
    container::OptimizationContainer,
    model::ServiceModel{PSY.AGC, PIDSmoothACE},
    areas::IS.FlattenIteratorWrapper{PSY.AGC},
)
    for ff in get_feedforwards(model)
        @debug "arguments" ff V _group = LOG_GROUP_FEEDFORWARDS_CONSTRUCTION
        add_feedforward_arguments!(container, model, areas, ff)
    end
    return
end

function add_feedforward_constraints!(
    container::OptimizationContainer,
    model::ServiceModel{PSY.AGC, PIDSmoothACE},
    areas::IS.FlattenIteratorWrapper{PSY.AGC},
)
    for ff in get_feedforwards(model)
        @debug "arguments" ff V _group = LOG_GROUP_FEEDFORWARDS_CONSTRUCTION
        add_feedforward_constraints!(container, model, areas, ff)
    end
    return
end


text:
module Williamson2ShallowWaterIMEXTests

using FillArrays
using Test
using WriteVTK
using Gridap
using GridapGeosciences

include("Williamson2InitialConditions.jl")

# Solves the steady state Williamson2 test case for the shallow water equations on a sphere
# of physical radius 6371220m. Involves a modified coriolis term that exactly balances
# the potential gradient term to achieve a steady state
# reference:
# D. L. Williamson, J. B. Drake, J. J.HackRüdiger Jakob, P. N.Swarztrauber, (1992)
# J Comp. Phys. 102 211-224

l2_err_u = [0.009874787793428196, 0.0028512787026009232, 0.0007426759094111235]
l2_err_h = [0.005610917625840937, 0.0014555948928746177, 0.0003683350733739   ]

order  = 1
degree = 4

for i in 1:3
  n      = 2*2^i
  nstep  = 5*n
  Uc     = sqrt(g*H₀)
  dx     = 2.0*π*rₑ/(4*n)
  dt     = 0.05*dx/Uc
  println("timestep: ", dt)   # gravity wave time step

  model = CubedSphereDiscreteModel(n; radius=rₑ)
  hf, uf = shallow_water_imex_time_stepper(model, order, degree,
                                      h₀, u₀, f₀, g,
                                      dt, 0.0, nstep;
                                      write_solution=false,
                                      write_solution_freq=5,
                                      write_diagnostics=true,
                                      write_diagnostics_freq=1,
                                      dump_diagnostics_on_screen=true)

  Ω     = Triangulation(model)
  dΩ    = Measure(Ω, degree)
  hc    = CellField(h₀, Ω)
  e     = h₀-hf
  err_h = sqrt(sum(∫(e⋅e)*dΩ))/sqrt(sum(∫(hc⋅hc)*dΩ))
  uc    = CellField(u₀, Ω)
  e     = u₀-uf
  err_u = sqrt(sum(∫(e⋅e)*dΩ))/sqrt(sum(∫(uc⋅uc)*dΩ))
  println("n=", n, ",\terr_u: ", err_u, ",\terr_h: ", err_h)

  @test abs(err_u - l2_err_u[i]) < 10.0^-12
  @test abs(err_h - l2_err_h[i]) < 10.0^-12
end

end

max_stars_repo_name: santiagobadia/GridapGeosciences

max_issues_repo_name: santiagobadia/GridapGeosciences

max_forks_repo_name: santiagobadia/GridapGeosciences


text:
# function used to perform the chi square test of independence upon
# the coappearance matrix Q
function HYOPTHESISCHI(Q::Array{Int64,2},
                       I::Int64,
                       J::Int64,
                       sig::Float64,
                       sum_cond_sku::Array{Int64,1})
    # Initialise chi square test
    M = convert(Int64,((I^2-I)/2))

    ## Determine the acceptance level of the test
    accept = cquantile(Chisq(1), sig/M)

    ## Create the arrays for the export of the chi-square test
    chi_yy = convert.(Float64,Q)
    norm   = Matrix{Float64}(undef,I,I) .= 0
    dep    = Matrix{Float64}(undef,I,I) .= 0

    @inbounds for i = 2:I
        for j = 1:i-1
            dep[i,j], norm[i,j] = chi_values(chi_yy[i,j],
                                            sum_cond_sku[i],
                                            sum_cond_sku[j],
                                            J, accept)
        end
    end
    dep = dep' + dep
    norm = norm' + norm
    return dep::Array{Float64,2},
           norm::Array{Float64,2}
end

function chi_part(a::Float64,b::Float64)
    (a-b)^2/b
end

function chi_part(a::Int64,b::Float64)
    (a-b)^2/b
end

function chi_values(chi_yy::Float64,
                    sum_cond_sku_i::Int64,
                    sum_cond_sku_j::Int64,
                    J::Int64,
                    accept::Float64)
    chi = 0.0
    dep = 0.0
    norm = (sum_cond_sku_i * sum_cond_sku_j)/J
    if chi_yy > norm
        chi_nr = J - sum_cond_sku_i
        chi_nd = J - sum_cond_sku_j
        chi_yn = sum_cond_sku_i - chi_yy
        chi_ny = sum_cond_sku_j - chi_yy
        chi_nn = chi_nd - chi_yn
        ind_yn = (sum_cond_sku_i * chi_nd)/J
        ind_ny = (sum_cond_sku_j * chi_nr)/J
        ind_nn = (chi_nr * chi_nd)/J
        chi += chi_part(chi_yy,norm)
        chi += chi_part(chi_yn,ind_yn)
        chi += chi_part(chi_ny,ind_ny)
        chi += chi_part(chi_nn,ind_nn)
        if chi > accept
            dep = chi_yy - norm
        else
            norm = chi_yy
        end
    end
    return dep, norm
end

# function to calculate the weight of each warehouse. It shows us the density of
# the independent coappearances in each warehouse if we were to allocate
# all SKUs according to the highest independent coappearances.
function WHWEIGHT(capacity::Array{Int64,1},
                  sum_nor::Array{Float64,1})
    sort_sum_nor = sort!(copy(sum_nor), rev = true)
    sort_sum_nor = vcat(sort_sum_nor,Array{Float64,1}(undef,sum(capacity)-I).=0)
    weight = Array{Float64,1}(undef,size(capacity)) .= 0
    weight[1] = sum(sort_sum_nor[1:capacity[1]])/sum(sort_sum_nor)
    for k = 2:size(capacity,1)
        weight[k] = sum(sort_sum_nor[sum(capacity[1:k-1]):sum(capacity[1:k])])/
                    sum(sort_sum_nor)
    end
    return weight::Array{Float64,1}
end

# function to find the largest warehouse that still has
# space left
function WHSPACE(cap_left::Array{Int64,1},
                 space::Int64)
    k_ind = 0
    for k in 1:size(cap_left,1)
        if cap_left[k] >= space
            k_ind = k
            break
        end
    end
    return k_ind::Int64
end

# function to select the best SKU for an allocation
function SELECTIK(sum_dep::Array{Float64,1},
                  sum_nor::Array{Float64,1},
                  weight::Array{Float64,1},
                  cap_left::Array{Int64,1},
                  X::Array{Bool,2},
                  dep::Array{Float64,2},)
    i       = findmax(sum_nor)[2]
    if sum(X[i,:]) > 0
        alloc = sum(X, dims = 2)
        for j = 1:size(X,1)
            if alloc[j] == 0
                i = j
                break
            end
        end
    end
    k_ind   = WHSPACE(cap_left,1)
    k_max   = findmax(cap_left)[2]
    pot_dep = WHPOTDEP(cap_left,i,X,dep)
    k_dep   = findmax(pot_dep)[2]
    if pot_dep[k_dep] > 0 && k_ind != k_ind
        sum_dep[i] + sum_nor[i] * weight[k_dep] >
        sum_nor[i] * weight[k_ind]
            k = k_dep
        elseif k_max != k_ind && sum_dep[i] + sum_nor[i] * weight[k_max] >
            sum_nor[i] * weight[k_ind]
                k = k_max
        else
            k = k_ind
    end
    return i::Int64,
           k::Int64
end

# function to check for each warehouse with free space whether
# SKU i has significant dependencies to other already allocated SKUs
function WHPOTDEP(cap_left::Array{Int64,1},
                  i::Int64,
                  X::Array{Bool,2},
                  dep::Array{Float64,2})
    pot_dep = Array{Float64,1}(undef,size(cap_left,1)) .= 0
    for k in 1:size(cap_left,1)
        if cap_left[k] > 0
            pot_dep[k] = CALCVAL(X,dep,i,k)
        end
    end
    return pot_dep::Array{Float64,1}
end

# function to remove every so far assigned dependent SKU-pair from
# the coappearance matrix dep to prevent the allocation bias described
# in our article.
function REMOVEALLOC(X::Array{Bool,2},
                     cap_left::Array{Int64,1},
                     nor::Array{Float64,2},
                     dep::Array{Float64,2})
    Qs = copy(dep) .= 0
    for a = 1:size(Qs,1)
        for b = 1:size(Qs,1)
            for k = 1:size(cap_left,1)
                if Qs[a,b] > 0
                    if X[a,k] * X[b,k] == 1
                        Qs[a,b] = 0
                    end
                end
            end
        end
    end
    Qs = nor .+ Qs
    return Qs::Array{Float64,2}
end

# function to check for all unallocated SKUs whether they have positive
# dependencies to the SKUs in the warehouse k the last SKU was allocated
# to. If so, check whether the dependencies are expected to dominate the
# independent coapperances. If yes, allocate the corresponding SKUs to
# the warehouse k.
function ADDDEPENDENT!(X::Array{Bool,2},
                       cap_left::Array{Int64,1},
                       k::Int64,
                       dep::Array{Float64,2},
                       nor::Array{Float64,2},
                       sum_dep::Array{Float64,1},
                       sum_nor::Array{Float64,1})
    add = 1
    while add == 1 && cap_left[k] > 0 && sum(X) < size(X,1)
        add = 0
        pot_dep = Array{Float64,1}(undef,size(X,1)) .= 0
        pot_nor = Array{Float64,1}(undef,size(X,1)) .= 0
        FINDDEP!(X::Array{Bool,2},
                dep::Array{Float64,2},
                nor::Array{Float64,2},
                k::Int64,
                pot_dep::Array{Float64,1},
                pot_nor::Array{Float64,1})
        i = findmax(pot_dep)[2]
        if findmax(pot_dep)[1] > 0
            if pot_dep[i] >= findmax(pot_nor)[1]
                ALLOCATEONE!(X::Array{Bool,2},
                             sum_dep::Array{Float64,1},
                             sum_nor::Array{Float64,1},
                             cap_left::Array{Int64,1},
                             i::Int64,
                             k::Int64)
                add = 1
            end
        end
    end
end

# function to allocate a selcted product to a selected warehouse
function ALLOCATEONE!(X::Array{Bool,2},
                      sum_dep::Array{Float64,1},
                      sum_nor::Array{Float64,1},
                      cap_left::Array{Int64,1},
                      i::Int64,
                      k::Int64)
    if sum(X[i,:]) == 0
        X[i,k] = 1
        sum_dep[i] = 0
        sum_nor[i] = 0
        cap_left[k] -= 1
    else
        error("This product is already allocated!")
    end
end

## check whether all warehouse except the last one are already full. If
## that is the case just allocate the remaining SKUs yet not allocated there.
function FILLLAST!(X::Array{Bool,2},cap_left::Array{Int64,1})
    go = 1
    for i = 1:size(cap_left,1)-1
        if cap_left[i] == 0 && go == 1
            go = 1
        else
            go = 0
        end
    end
    if go == 1
        allocated = sum(X,dims=2)
        for i = 1:length(allocated)
            if allocated[i] == 0
                X[i,size(cap_left,1)] = 1
                cap_left[size(cap_left,1)] -= 1
            end
        end
    end
end

# function to check the dependencies to already allocated SKUs
function FINDDEP!(X::Array{Bool,2},
                  dep::Array{Float64,2},
                  nor::Array{Float64,2},
                  k::Int64,
                  pot_dep::Array{Float64,1},
                  pot_nor::Array{Float64,1})
    allocated = sum(X,dims=2)
    @inbounds @simd for i in 1:size(X,1)
        if allocated[i] == 0
            pot_dep[i] = CALCVAL(X,dep,i,k)
            pot_nor[i] = CALCVAL(X,nor,i,k)
            if pot_dep[i] > 0
                pot_dep[i] += pot_nor[i]
            end
        end
    end
    return pot_dep::Array{Float64,1},
           pot_nor::Array{Float64,1}
end

function CALCVAL(X::Matrix{Bool},T::Matrix{Float64},i::Int64,k::Int64)
    out = 0
    @avxt for y = 1:size(X,1)
        out += X[y,k] * T[y,i]
    end
    return out
end

# function to allocate the SKUs with the highest potential allocation
# value to each warehouse with leftover storage space until it is full
function FILLUP!(X::Array{Bool,2},
                Q::Array{Float64,2},
                capacity_left::Array{Int64,1})
    for d = 1:size(capacity_left,1)
        while capacity_left[d] > 0
            best_allocation = Array{Float64,1}(undef,size(Q,1)) .= 0
            for i = 1:size(Q,1)
                if X[i,d] == 0
                    best_allocation[i] = CALCVAL(X,Q,i,d)
                end
            end
            if findmax(best_allocation)[1] > 0
                X[findmax(best_allocation)[2],d] = 1
                capacity_left[d] -= 1
            else
                capacity_left[d] = 0
            end
        end
    end
end

# function to apply a pair-wise exchange local search on the allocation
# of the CHI heuristic
function LOCALSEARCHCHI(W::Array{Int64,2},
                        Q::Array{Int64,2})
    coapp_sort = vec(sum(Q,dims = 2))
    coapp_sort = sortperm(coapp_sort,rev=true)
    for trial = 1:3
        improvement = 0
        for k = 1:size(W,2)-1
            @inbounds for i in coapp_sort
                if W[i,k] == 1 && W[i,k+1] == 0
                    for j in coapp_sort
                        if W[j,k+1] == 1 && W[j,k] == 0
                            iteration = POTENTIALCHI!(i,j,k,k+1,W,Q)
                            if iteration > 0
                                improvement += 1
                                break
                            end
                        end
                    end
                end
            end
        end
    end
   return W::Array{Int64,2}
end

function POTENTIALCHI!(i::Int64,
                       j::Int64,
                       k::Int64,
                       g::Int64,
                       X::Array{Int64,2},
                       Q::Array{Int64,2})
    potential = 0
    @avxt for y = 1:size(Q,1)
        potential += (X[y,g] - X[y,k]) * Q[y,i] + (X[y,k] - X[y,g]) * Q[y,j]
    end
    if potential > 0
        impr = 1
        X[i,k] = 0
        X[j,g] = 0
        X[i,g] = 1
        X[j,k] = 1
    else
        impr = 0
    end
    return impr::Int64
end

max_stars_repo_name: beyondatlas/splitdeliveries

max_issues_repo_name: beyondatlas/splitdeliveries

max_forks_repo_name: beyondatlas/splitdeliveries


text:
"""
# Usage

    solve(prob::MultiTermsFODEProblem, ::ClosedFormHankelM)

Use Closed-Form Hankel matrix algorithm to obtain numerical solution at zero initial condition.
"""
struct ClosedFormHankelM <: FractionalDiffEqAlgorithm end

function solve(prob::MultiTermsFODEProblem, h, ::ClosedFormHankelM)
    @unpack parameters, orders, rightfun, rparameters, rorders, u0, tspan = prob
    t0 = tspan[1]; T = tspan[2]
    t = collect(t0:h:T)
    u = rightfun.(t)
    u = u[:]
    A, B = 0, 0

    g = genfun(1)
    nt = length(t)
    n = length(parameters)
    m = length(rparameters)
    for i=1:n
        A = A .+ getvec(orders[i], nt, g)*parameters[i]/(h^orders[i])
    end

    for i=1:m
        B = B .+ getvec(rorders[i], nt, g)*rparameters[i]/(h^rorders[i])
    end

    A = rotl90(newhankel(A[end:-1:1]))
    B = rotl90(newhankel(B[end:-1:1]))

    y = B*inv(A)*u
    return FODESolution(t, y)
end

function newhankel(v)
    n = length(v)
    v = v[:]

    hankelm = zeros(n, n)
    for i=1:length(v)
        hankelm[i, 1:end-i+1] = v[i:end]
    end

    return hankelm
end


text:
# we write a function that solves Ly = b for y where L is lower triangular and b is a vector
function Lsolve(L, b)

    # we get the sizes of L and b
    LRows, LCols = size(L)
    bRows = size(b)[1]

    # we check that L and b ahve appropriate sizes
    if LRows != LCols || bRows != LRows
        error("The size of L or b is not appropriate")
    end

    # we initialize y to b
    y = copy(b);

    # we loop over the rows of L
    for k = 1:LRows

        # we check that we dont divide by zero
        if L[k,k] == 0
            error("There are zeros on the diagonal of L")
        end

        # we modify y accordingly
        for j = 1:k-1
            y[k] -= L[k,j] * y[j];
        end
        y[k] /= L[k,k];
    end

    # we return y
    return y
end

max_stars_repo_path: 1a/LU Decomposition/Lsolve.jl

max_issues_repo_path: 1a/LU Decomposition/Lsolve.jl

max_forks_repo_path: 1a/LU Decomposition/Lsolve.jl


text:
@inline ∇_dot_qᶜ(i, j, k, grid::AbstractGrid{FT}, ::Nothing, args...) where FT = zero(FT)
@inline ∂ⱼ_τ₁ⱼ(i, j, k, grid::AbstractGrid{FT}, ::Nothing, args...) where FT = zero(FT)
@inline ∂ⱼ_τ₂ⱼ(i, j, k, grid::AbstractGrid{FT}, ::Nothing, args...) where FT = zero(FT)
@inline ∂ⱼ_τ₃ⱼ(i, j, k, grid::AbstractGrid{FT}, ::Nothing, args...) where FT = zero(FT)

calculate_diffusivities!(K, arch, grid, closure::Nothing, args...) = nothing


text:
## Transmission Loss Demonstration

## Preamble
using Plots

## Lloyds Mirror
include("../LloydsMirror.jl")

f = 150.
r_src = 0.
z_src = 25.
z = range(0, 4e2, length = 500)
r = range(0, 5e3, length = 500)
c = 1500.

TL(r, z) = LloydsMirror.lloydsmirror_singlereflection.(c, f, r_src, r, z_src, z)

ptField = heatmap(r, z, TL,
	seriescolor = cgrad(:jet, rev = true),
	title = "TL (dB)",
	yaxis = (:flip, "Depth (m)"),
	xaxis = ("Range (m)"))
display(ptField)
savefig(pt, "img/SoundField_LloydsMirror_Simple.png")

ptSlice = plot(r, r -> TL(r, 200.),
	yaxis = :flip)
plot!(r, r -> 20log10(r))
display(ptSlice)
savefig(ptSlice, "img/SoundField_LloydsMirror_Simple_Slice.png")


text:
incx = 1
incy = 1


text:
tutorials_dir = joinpath(pkgdir(TuringGLM), "docs", "src", "tutorials")

"""Run all Pluto notebooks (".jl" files) in `tutorials_dir` and write outputs to HTML files."""
function build()
    println("Building tutorials")
    # Evaluate notebooks in the same process to avoid having to recompile from scratch each time.
    # This is similar to how Documenter and Franklin evaluate code.
    # Note that things like method overrides and other global changes may leak between notebooks!
    use_distributed = false
    output_format = documenter_output
    bopts = BuildOptions(tutorials_dir; use_distributed, output_format)
    parallel_build(bopts)
    return nothing
end

"Return Markdown file links which can be passed to Documenter.jl."
function markdown_files()
    md_files = map(tutorials) do tutorial
        file = lowercase(replace(tutorial, " " => '_'))
        return joinpath("tutorials", "$file.md")
    end
    return md_files
end


text:
using LinearMaps
"""
A(σ₀)δψ=-A(δσ)ψ₀
Therefore,
perturbed data: δψ=-A⁻¹(σ₀)A(δσ)ψ₀
"""
function bornmod!(pa::ParamExpt, δσ, σ0)
	updateA!(pa.paσ, δσ)
	updateA!(pa.paσ0, σ0)
	# pa.paσ.A[end,:] .= 0.0 # is it necessary? A=A(σ₀)+A(δσ), except for the last row
	qrA0=factorize(pa.paσ0.A)
	@showprogress 5 "time loop ψ\t" for it in 1:length(pa.tgrid)
		snap_in=view(pa.LP,:,:,it)
		# generate ψ0
		forwψ!(pa.ψ0,snap_in,pa.paσ0,qrA0)

		# compute A(δσ)*ψ0, store this in adjψ
		applyA!(pa.adjψ, pa.ψ0, pa.paσ; A=pa.paσ.A)

		# compute A(σ0)⁻¹ * adjψ
		# do we need to mute boundary source?! --- otherwise it seems to bug!
		applyinvA!(pa.ψ, pa.adjψ, pa.paσ0; A=qrA0, mute_boundary_source=false)

		# record
		dat_slice=view(pa.data,it,:)
		mul!(dat_slice,pa.ACQ,pa.ψ) # record
		rmul!(dat_slice, -1.0)
	end
	return nothing
end

# just want to extract the data out into δd
function Fborn_map!(δd, δσ, pa, σ0)
	bornmod!(pa, δσ, σ0)
	copyto!(δd, pa.data)
end

function Fadj_map!(δσ, δd, pa, σ0)
	updateA!(pa.paσ0, σ0)
	qrAt0=factorize(pa.paσ0.At)
	qrA0=factorize(pa.paσ0.A)
	fill!(pa.g, 0.0)
	@showprogress 10 "time loop ψ\t" for it in 1:length(pa.tgrid)
		snap_in=view(pa.LP,:,:,it)
		# generate ψ0
		forwψ!(pa.ψ0,snap_in,pa.paσ0,qrA0)

		# prepare the adjoint source
		for ir in 1:size(pa.ACQ,1)
			pa.data_misfit[ir]=δd[it+(ir-1)*length(pa.tgrid)]
		end
		# migrate the data!; using this routine will also test its consistency
		adjψ_core!(pa.adjψ, pa.adjψ2, pa.ψ0,
		      pa.adjsrc, pa.g, pa.gtemp,
		      pa.ACQ, pa.data_misfit, pa.paσ0, qrAt0)

	end
	copyto!(δσ,pa.g)
	return nothing
end

"""
Return a `LinearMap` object to perform linearized modeling and its transpose.

# Arguments
* `pa::PoissonExpt`
* `σ0::Array`
"""
function operator_Born(pa::ParamExpt, σ0)
	@assert length(σ0)==prod(length.(pa.mgrid))
	fw=(y,x)->Fborn_map!(y, x, pa, σ0)
	bk=(y,x)->Fadj_map!(y, x, pa, σ0)

	return LinearMap(fw, bk,
		  length(pa.data),  # length of output
		  length(pa.σ), # length of input
		  ismutating=true)
end


text:
export G

"""
Instantiate a Geometric Algebra with 'p' positive, 'q' negative and 'r' zero dimensions.

'q' and 'r' defaults to zero if omitted.

The basis will be named with prefix `e` and indexed starting from 1.

For example:

```julia
# Basic
Hyper   = G(1)       # Hyperbolic numbers.
ℂ       = G(0,1)     # Complex numbers.
Dual    = G(0,0,1)   # Dual numbers.
ℍ       = G(0,2)     # Quaternions.

# Clifford
Cl2 =       G(2)     # Clifford algebra for 2D vector space.
Cl3 =       G(3)     # Clifford algebra for 3D vector space.
Spacetime = G(1,3)   # Clifford algebra for timespace vectors.

# Geometric
PGA2D = G(2,0,1)     # Projective Euclidean 2D plane. (dual)
PGA3D = G(3,0,1)     # Projective Euclidean 3D space. (dual)
CGA2D = G(3,1)       # conformal 2D space.
CGA3D = G(4,1)       # Conformal 3D space.

# High-Dimensional GA
DCGA3D = G(6,2)      # Double Conformal 3D Space.
TCGA3D = G(9,3)      # Tripple Conformal 3D Space.
DCGSTA = G(4,8)      # Double Conformal Geometric Space Time Algebra.
QCGA   = G(9,6)      # Quadric Conformal Geometric Algebra.
```
"""
function G(p::Integer, q::Integer, r::Integer)
    total = p + q + r
    basis = "e*" * join([string(i) for i in 1:total], "|")
    metric = [fill(1, p) ; fill(-1, q) ; fill(0, r)]
    galgebra.ga.Ga(basis, g=metric)
end

G(p::Integer, q::Integer) = G(p, q, 0)

G(p::Integer) = G(p, 0, 0)


text:
function read_grammar()
      f = open(joinpath(@__DIR__, "grammar.txt"), "r")
      buffer = read(f, String)
      close(f)
      split(buffer)
end

function isLiteral(word::AbstractString)
      startswith(word, "\"") && endswith(word, "\"")
end

function isTerminal(word::AbstractString)
      startswith(word, "<") && endswith(word, ">")
end

struct Grammar
      𝐕ₜ::Set{Terminal}
      𝐕ₙ::Set{Nonterminal}
      𝐏::Dict{Nonterminal,Set{Statement}}
end

function show(io::IO,x::Grammar)
      print(io,"𝐆( 𝐕ₜ, 𝐕ₙ, 𝐏 ), ")
      print(io,"|𝐕ₜ| = $(length(x.𝐕ₜ)) ,")
      print(io,"|𝐕ₙ| = $(length(x.𝐕ₙ)) ,")
      ps=sum(x.𝐏) do d
            length(d.second)
      end
      println(io,"|𝐏| = $ps")
end

function prettyshow(io::IO,x::Grammar)
      show(io,x)
      print(io,"Terminal\t")
      foreach(x->print(io,x," "),x.𝐕ₜ)
      print(io,"\nNonterminal\t")
      foreach(x->print(io,x," "),x.𝐕ₙ)
      println(io)
      for (k,v) in x.𝐏
            print(io,k)
            for (i,p) in enumerate(v)
                  if i==1 print(io,"\n\t:  ")
                  else  print(io,"\n\t|  ") end
                  for w in p
                        print(io,string(w)," ")
                  end
            end
            println(io,"\n;")
      end
end

function chooseToken!(𝐕ₜ, 𝐕ₙ, word)
      sym=:ϵ
      if isTerminal(word)
            sym=Symbol(word)
            push!(𝐕ₜ, sym)
      elseif isLiteral(word)
            sym=Symbol(word[2:end - 1])
            push!(𝐕ₜ, sym)
      else
            sym=Symbol(word)
            push!(𝐕ₙ, sym)
      end
      sym
end

function Grammar(buffer)
      𝐕ₜ = Set{Terminal}()
      𝐕ₙ = Set{Nonterminal}()
      𝐏 = Dict{Nonterminal,Set{Statement}}()
      stack = 0
      current = :ϵ
      production = Symbol[]
      for word in buffer
            if word == "{"
                  stack+=1
            elseif word == "}"
                  stack-=1
            elseif stack!=0 continue
            elseif current==:ϵ
                  sym=chooseToken!(𝐕ₜ,𝐕ₙ,word)
                  current=Symbol(sym)
                  𝐏[current]=Set{Statement}()
            elseif word ==":" continue
            elseif word =="|"
                  push!(𝐏[current],production)
                  production=Symbol[]
            elseif word ==";"
                  push!(𝐏[current],production)
                  production=Symbol[]
                  current=:ϵ
            else
                  sym=chooseToken!(𝐕ₜ,𝐕ₙ,word)
                  push!(production,sym)
            end
      end
      Grammar(𝐕ₜ,𝐕ₙ,𝐏)
end

max_stars_repo_name: yeruoforever/SimpleGrammar

max_issues_repo_name: yeruoforever/SimpleGrammar

max_forks_repo_name: yeruoforever/SimpleGrammar


text:
module CSFML

include("LibCSFML.jl")
using .LibCSFML

end # module


text:
load("debug.jl")

module TestGlobalScope
import Base.*
import Debug.*

function debug_hook(line::Int, file, scope::Scope)
    print(line, ":")
    if line < 30
        if (line >  20) debug_eval(scope, :( print("\tx = ", x) )) end
        if (line >  21) debug_eval(scope, :( print("\ty = ", y) )) end
        if (line >  20) debug_eval(scope, :( print("\tz = ", z) )) end
        if (line == 21) debug_eval(scope, :( (x,y,z) = (7,8,9)  )) end
    end
    println()
end

@debug begin
    x, z = 1, 2  # line 20
    y = 2        #      21
    x            #      22
end

println("(x, y, z) = ", (x, y, z))

try
    @debug begin
        x=1
    end
    error("Shouldn't be able to use @debug in local scope!")
catch e
    print("@debug was disallowed in local scope, as it should.")
end

end  # module