text:
"""
offonPDF(t::Vector,r::Vector,n::Int,nr::Int)

Active (ON) and Inactive (OFF) time distributions for GRSM model
Takes difference of ON and OFF time CDF to produce PDF
method = 1 uses DifferentialEquations.jl
"""
function offonPDF(t::Vector,r::Vector,n::Int,nr::Int,method)
    T,TA,TI = mat_GSR_DT(r,n,nr)
    pss = normalized_nullspace(T)
    SA=ontimeCDF(t,n,nr,TA,pss,method)
    SI=offtimeCDF(t,r,n,nr,TI,pss,method)
    return pdf_from_cdf(t,SI), pdf_from_cdf(t,SA)
end
function offonPDF_offeject(t::Vector,r::Vector,n::Int,nr::Int,method)
    T,TA,TI = mat_GSR_DT_offeject(r,n,nr)
    pss = normalized_nullspace(T)
    SA=ontimeCDF(t,n,nr,TA,pss,method)
    SI=offtimeCDF(t,r,n,nr,TI,pss,method)
    return pdf_from_cdf(t,SI), pdf_from_cdf(t,SA)
end

"""
pdf_from_cdf(t,S)
"""
function pdf_from_cdf(t,S)
    P = diff(S)
    P/(sum(P)*(t[2]-t[1]))
end

"""
ontimeCDF(tin::Vector,n::Int,nr::Int,rin::Vector,TA::Matrix,pss::Vector)
offtimeCDF(tin::Vector,n::Int,nr::Int,r::Vector,TI::Matrix,pss::Vector)

ON(OFF) dwell time distributions of GRS model
Found by computing accumulated probability into OFF(ON) states
where transitions out of OFF(ON) states are zeroed, starting from first instance of ON(OFF) state
weighted by steady state distribution (i.e. solving a first passage time problem)x
"""
function ontimeCDF(tin::Vector,n::Int,nr::Int,TA::Matrix,pss::Vector,method)
    t = [tin ; tin[end] + tin[2]-tin[1]] #add a time point so that diff() gives original length
    SAinit = init_prob(pss,n,nr)
    SA = time_evolve(t,TA,SAinit,method)
    accumulate(SA,n,nr)  # accumulated prob into OFF states
end
function offtimeCDF(tin::Vector,r::Vector,n::Int,nr::Int,TI::Matrix,pss::Vector,method)
    t = [tin ; tin[end] + tin[2]-tin[1]]
    nonzerosI = nonzero_rows(TI)  # only keep nonzero rows to reduce singularity of matrix
    TI = TI[nonzerosI,nonzerosI]
    SIinit = init_prob(pss,r,n,nr,nonzerosI)
    SI = time_evolve(t,TI,SIinit,method)
    accumulate(SI,n,nr,nonzerosI,length(pss)) # accumulated prob into ON states
end
function offtimeCDF(tin::Vector,r::Vector,n::Int,TI::Matrix,method)
    t = [tin ; tin[end] + tin[2]-tin[1]]
    SIinit = zeros(n+1)
    SIinit[n] = 1.
    SI = time_evolve(t,TI,SIinit,method)
    return SI[:,end]
end

"""
steady_state_offpath(rin::Vector,n::Int,nr::Int,nhist::Int,nalleles::Int)
GRS model where mRNA decay rate is accelerated to account for nonviability of off-pathway pre-mRNA
from RNA that is recursively spliced
"""
function steady_state_offpath(rin::Vector,n::Int,nr::Int,nhist::Int,nalleles::Int)
    r = copy(rin)
    nu = get_nu(r,n,nr)
    eta = get_eta(r,n,nr)
    # r[end] /= survival_fraction(nu,eta,nr)
    yieldfactor = survival_fraction(nu,eta,nr)
    mhist=steady_state(r,n,nr,nhist,nalleles)
    technical_loss(mhist,yieldfactor,nhist)
end
"""
steady_state(rin::Vector,n::Int,nr::Int,nhist::Int,nalleles::Int)
Steady state distribution of mRNA in GRM model (which is the same as GRSM model)
"""
function steady_state(rin::Vector,n::Int,nr::Int,nhist::Int,nalleles::Int)
    r = rin/rin[end]
    gammap,gamman = get_gamma(r,n)
    nu = get_nu(r,n,nr)
    T,B = transition_rate_mat(n,nr,gammap,gamman,nu)
    P = initial_pmf(T,nu[end],n,nr,nhist)
    mhist=steady_state(nhist,P,T,B)
    mhist = marginalize(mhist)
    allele_convolve(mhist[1:nhist],nalleles) # Convolve to obtain result for n alleles
end
"""
steady_state_offeject(rin::Vector,n::Int,nr::Int,nhist::Int,nalleles::Int)
Steady state distribution of mRNA in GRM model (which is the same as GRSM model)
"""
function steady_state_offeject(rin::Vector,n::Int,nr::Int,nhist::Int,nalleles::Int)
    r = rin/rin[end]
    gammap,gamman = get_gamma(r,n)
    nu = get_nu(r,n,nr)
    eta = get_eta(r,n,nr)
    T,B = transition_rate_mat_offeject(n,nr,gammap,gamman,nu,eta)
    P = initial_pmf(T,nu[end],n,nr,nhist)
    mhist=steady_state(nhist,P,T,B)
    mhist = marginalize(mhist)
    allele_convolve(mhist[1:nhist],nalleles) # Convolve to obtain result for n alleles
end
"""
steady_state(nhist::Int,nalleles::Int,ejectrate,P,T,B,tol = 1e-6)
Iterative algorithm for computing null space of truncated transition rate matrix
of Master equation of GR model to give steady state of mRNA in GRM model
for single allele
"""
function steady_state(nhist::Int,P::Matrix,T::Matrix,B::Matrix,tol = 1e-8,stepmax=10000)
    total = size(P,2)
    steps = 0
    err = 1.
    A = T - B
    while err > tol && steps < stepmax
        P0 = copy(P)
        P[:,1] = -A\P[:,2]
        for m = 2:total-1
            P[:,m] = @inbounds (-A + UniformScaling(m-1))\((B*P[:,m-1]) + m*P[:,m+1])
        end
        P[:,total] = (-A + UniformScaling(total-1))\((B*P[:,total-1]))
        P /=sum(P)
        err = norm((P-P0))
        steps += 1
    end
    return P  # marginalize over GR states with marginalize(P) = sum(P,dims=1)
end

"""
steady_state(r,n,nhist,nalleles)
Steady State of mRNA in G (telegraph) model
"""
function steady_state(r::Vector,yieldfactor::Float64,n::Int,nhist::Int,nalleles::Int)
    nh = nhist_loss(nhist,yieldfactor)
    mhist = steady_state(r,n,nh,nalleles)
    technical_loss(mhist,yieldfactor,nhist)
end

function steady_state(r::Vector,n::Int,nhist::Int,nalleles::Int)
    if n == 0
        return steady_state_n0(r,nhist,nalleles)
    else
        if nhist > 1800
            method = false
        else
            method = true
        end
        P = steady_state_full(r,n,nhist,method)
        return steady_state_rna(P,n,nhist,nalleles)
    end
end

"""
steady_state_n0(r::Vector,n::Int,nhist::Int)

steady state for n=0 (i.e. Poisson model)
"""
function steady_state_n0(r,nhist)
    mhist = Vector{Float64}(undef,nhist+2)
    k = r[1]/r[2]
    d = Poisson(k)
    for i in 1:nhist+2
        mhist[i] = pdf(d,i-1)
    end
    return mhist
end
steady_state_n0(r,nhist,nalleles) = allele_convolve(steady_state_n0(r,nhist),nalleles)[1:nhist]

"""
steady_state_full(r::Vector,n::Int,nhist::Int,method::Bool)
Steady State of full G (telegraph) model

method = true: use exact method
method = false: use fast method
"""

function steady_state_full(r::Vector,n::Int,nhist::Int,method::Bool)
    if method
        steady_state_full(r,n,nhist)
    else
        steady_state_fast(r,n,nhist)
    end
end

function steady_state_full(r::Vector,n::Int,nhist::Int)
    M = mat_GM(r,n,nhist)
    return normalized_nullspace(M)
end
"""
steady_state_rna(P,n,nhist,nalleles)
Steady State of rna distribution for G model accounting for number of alleles
"""
function steady_state_rna(P,n,nhist,nalleles)
    mhist = marginalize(P,n)
    allele_convolve(mhist,nalleles)[1:nhist]
end
"""
steady_state_fast(rin::Vector,n::Int,nhist::Int,nalleles::Int)
Steady State of G model accounting for number of alleles
using iterative block algorithm (may not converge to correct solution)
"""
# function steady_state_fast(rin::Vector,n::Int,nhist::Int,nalleles::Int)
#     mhist = marginalize(steady_state_fast(rin,n,nhist))
#     allele_convolve(mhist[1:nhist],nalleles) # Convolve to obtain result for n alleles
# end
function steady_state_fast(rin::Vector,n::Int,nhist::Int)
    r = rin/rin[2*n+2]
    gammap,gamman = get_gamma(r,n)
    nu = r[2*n+1]
    T,B = transition_rate_mat(n,gammap,gamman,nu)
    P = initial_pmf(T,nu,n,nhist)
    unfold(steady_state(nhist,P,T,B))
end

function checkP(A,B,P,total)
    err = abs.(A*P[:,1])
    for m = 2:total-1
        err += abs.(((A - UniformScaling(m-1))* P[:,m] + B*P[:,m-1] + m*P[:,m+1]))
    end
    err += abs.(((A - UniformScaling(total-1))* P[:,total] + B*P[:,total-1]))
    sum(err)/total
end
"""
nhist_loss(nhist,yieldfactor)

Compute length of pre-loss histogram
"""
nhist_loss(nhist,yieldfactor) = round(Int,nhist/yieldfactor)

"""
mat(r,n,nr,nhist)
Transition rate matrices for GSR model
return T,TA,TI
"""
function mat_GSR_DT(r,n,nr)
    gammap,gamman = get_gamma(r,n)
    nu = get_nu(r,n,nr)
    eta = get_eta(r,n,nr)
    T = transition_rate_mat_T(n,nr,gammap,gamman,nu,eta)
    transition_rate_mat(T,n,nr)
end
function mat_GSR_DT_offeject(r,n,nr)
    gammap,gamman = get_gamma(r,n)
    nu = get_nu(r,n,nr)
    eta = get_eta(r,n,nr)
    T,_ = transition_rate_mat_offeject(n,nr,gammap,gamman,nu,eta)
    transition_rate_mat(T,n,nr,2)
end
function mat_G_DT(r,n)
    gammap,gamman = get_gamma(r,n)
    T = transition_rate_mat(n,gammap,gamman)
    transition_rate_mat(T,n)
end
function mat_GSR_T(r,n,nr)
    gammap,gamman = get_gamma(r,n)
    nu = get_nu(r,n,nr)
    eta = get_eta(r,n,nr)
    transition_rate_mat_T(n,nr,gammap,gamman,nu,eta)
end
"""
mat_GM(r,n,nhist)
Transition rate matrix of GM model
"""
function mat_GM(r,n,nhist)
    if n == 0
        nu = r[1]/r[2]
        return transition_rate_mat(nu,nhist)
    else
        gammap,gamman = get_gamma(r,n)
        return transition_rate_mat(n,gammap,gamman, r[2*n+1],r[2*n+2],nhist)
    end
end
"""
get_rates_GSR(r,n,nr)
Convert rates in r to parameters used in the Transition rate matrices for GSR model
return gammap,gamman,nu,eta
"""
function get_rates_GSR(r,n,nr)
    gammap,gamman = get_gamma(r,n)
    nu = get_nu(r,n,nr)
    eta = get_eta(r,n,nr)
    return gammap,gamman,nu,eta
end

"""
transient(ts::Vector,r,n,nhist,nalleles,P0)

Compute mRNA pmf of GM model at times in vector ts starting
with initial condition P0
"""

"""
transient(t,r,yieldfactor,n,nhist,nalleles,P0::Vector,method)
transient(t,n,nhist,nalleles,P0,Mvals,Mvects)

Compute mRNA pmf of GM model at time t given initial condition P0
and eigenvalues and eigenvectors of model transition rate matrix
method = 1 uses JuliaDifferentialEquations.jl
method != 1 uses eigenvalue decomposition
"""
function transient(t,r,yieldfactor,n,nalleles,P0::Vector,method)
    mhist = transient(t,r,n,nalleles,P0,method)
    technical_loss!(mhist,yieldfactor)
    return mhist
end
function transient(t::Float64,r,n,nalleles,P0::Vector,method)
    P = transient(t,r,n,P0,method)
    mh = marginalize(P,n)
    allele_convolve(mh,nalleles)
end
function transient(t::Vector,r,n,nalleles,P0::Vector,method)
    P = transient(t,r,n,P0,method)
    # P += abs.(rand(size(P)))
    mhist = Array{Array,1}(undef,length(t))
    for i in 1:size(P,1)
        mh = marginalize(P[i,:],n)
        mhist[i]= allele_convolve(mh,nalleles)
    end
    return mhist
end
function transient(t,r,n,P0::Vector,method)
    nhist = Int(length(P0)/(n+1)) - 2
    M = mat_GM(r,n,nhist)
    time_evolve(t,M,P0,method)
end

function transient_delay(t::Vector,r0::Vector,r1::Vector,delay::Float64,n::Int,nalleles,P0::Vector)
    P = time_evolve_delay(t,r0,r1,delay,n,P0)
    mhist = Array{Array,1}(undef,length(t))
    for i in 1:size(P,1)
        mh = marginalize(P[i,1:end-1],n)
        mhist[i]= allele_convolve(mh,nalleles)
    end
    return mhist
end

"""
time_evolve(t,M::Matrix,Sinit::Vector)
Eigenvalue solution of Linear ODE with rate matrix T and initial vector Sinit
"""
function time_evolve(t,M::Matrix,S0::Vector,method)
    if method == 1
        return time_evolve_diff(t,M,S0)
    else
        return time_evolve(t,M,S0)
    end
end

function time_evolve(t,M::Matrix,Sinit::Vector)
    vals,vects = eig_decompose(M)
    weights = solve_vector(vects,Sinit)
    time_evolve(t,vals,vects,weights)
end
"""
time_evolve(t,vals::Vector,vects::Matrix,weights::Vector)

"""
function time_evolve(t::Float64,vals::Vector,vects::Matrix,weights::Vector)
    n = length(vals)
    S = zeros(n)
    for j = 1:n
        for i = 1:n
            S[j] += real(weights[i]*vects[j,i]*exp.(vals[i]*t))
        end
    end
    return S
end
function time_evolve(t::Vector,vals::Vector,vects::Matrix,weights::Vector)
    ntime = length(t)
    n = length(vals)
    S = Array{Float64,2}(undef,ntime,n)
    for j = 1:n
        Sj = zeros(ntime)
        for i = 1:n
            Sj += real(weights[i]*vects[j,i]*exp.(vals[i]*t))
        end
        S[:,j]=Sj
    end
    return S
end
"""
time_evolve_diff(t,M::Matrix,P0)

Solve transient problem using DifferentialEquations.jl
"""
function time_evolve_diff(t,M::Matrix,P0)
    global M_global = copy(M)
    tspan = (0.,t[end])
    prob = ODEProblem(fglobal,P0,tspan)
    # sol = solve(prob,saveat=t, lsoda(),abstol = 1e-4, reltol = 1e-4)
    sol = solve(prob,saveat=t, lsoda())
    return sol'
end

fglobal(u,p,t) = M_global*u

function time_evolve_delay(t,r0,r1,delay,n,P0)
    tspan = (0.,t[end])
    nhist = Int(length(P0)/(n+1)) - 2
    p  = [r0;r1;delay;n;nhist]
    P0 = [P0;1.]
    prob = ODEProblem(fdelay!,P0,tspan,p)
    sol = solve(prob,saveat=t)
    return sol'
end

function fdelay!(du,u,p,t)
    n = Int(p[end-1])
    nhist = Int(p[end])
    delay = p[end-3]
    r0 = p[1:2*(n+1)]
    r1 = p[2*(n+1)+1:4*(n+1)]
    r = r0*u[end] + (1-u[end])*r1
    M = mat_GM(r,n,nhist)
    du[1:end-1] = M*u[1:end-1]
    du[end] = -delay*u[end]
end

"""
technical_loss(mhist,yieldfactor)

Reduce counts due to technical loss
"""
function technical_loss!(mhist::Vector,yieldfactor)
    for i in eachindex(mhist)
        mhist[i] = technical_loss(mhist[i],yieldfactor,length(mhist[i]))
    end
end
"""
technical_loss(mhist,yieldfactor,nhist)

Reduce counts due to technical loss using Binomial sampling
"""
function technical_loss(mhist::Vector,yieldfactor,nhist)
    p = zeros(nhist)
    for m in eachindex(mhist)
        d = Binomial(m-1,clamp(yieldfactor,0.,1.))
        for c in 1:nhist
            p[c] += mhist[m]*pdf(d,c-1)
        end
    end
    normalize_histogram(p)
end
"""
technical_loss_poisson(mhist,yieldfactor,nhist)

Reduce counts due to technical loss using Poisson sampling
"""
function technical_loss_poisson(mhist,yieldfactor,nhist)
    p = zeros(nhist)
    for m in eachindex(mhist)
        d = Poisson(yieldfactor*(m-1))
        for c in 1:nhist
            p[c] += mhist[m]*pdf(d,c-1)
        end
    end
    normalize_histogram(p)
end

"""
additive_noise(mhist,noise,nhist)

Add Poisson noise to histogram
"""
function additive_noise(mhist,noise,nhist)
    p = zeros(nhist)
    d = Poisson(noise)
    for m in 1:nhist
        for n in 1:m
            p[m] += mhist[n]*pdf(d,m-n)
        end
    end
    normalize_histogram(p)
end
"""
threshold_noise(mhist,noise,yieldfactor,nhist)

Add Poisson noise to histogram then reduce counts due to technical loss

"""
function threshold_noise(mhist,noise,yieldfactor,nhist)
    h = additive_noise(mhist,noise,nhist)
    technical_loss(h,yieldfactor,nhist)
end

"""
solve_vector(A::Matrix,b::vector)
solve A x = b
If matrix divide has error higher than tol
use SVD and pseudoinverse with threshold
"""
function solve_vector(A::Matrix,b::Vector,th = 1e-16,tol=1e-1)
    x = A\b
    if norm(b-A*x,Inf) > tol
        M = svd(A)
        Sv = M.S
        Sv[abs.(Sv) .< th] .= 0.
        Sv[abs.(Sv) .>= th] = 1 ./ Sv[abs.(Sv) .>= th]
        x = M.V * diagm(Sv) * M.U' * b
    end
    return x[:,1] # return as vector
end
"""
initial_pmf(T,ejectrate,n,nr,nhist)

Tensor product of nullspace of T and Poisson density
with rate = mean ejection rate of mRNA
Used for initial condition of iteration algorithm to
compute steady state pmf of GRM Master Equation
"""
function initial_pmf(T,ejectrate,n,nr,nhist)
    t0 = normalized_nullspace(T)
    ejectrate *= sum(t0[(n+1)*2^(nr-1)+1:end])
    initial_pmf(t0,ejectrate,nhist)
end
function initial_pmf(T,ejectrate,n,nhist)
    t0 = normalized_nullspace(T)
    initial_pmf(t0,ejectrate,nhist)
end
function initial_pmf(t0,ejectrate,nhist)
    d = Poisson(ejectrate)
    nT = length(t0)
    total = nhist + 2
    P = Array{Float64,2}(undef,nT,total)
    for m = 1:total
        P[:,m] = t0 * pdf(d,m)
    end
    P
end
"""
get_gamma(r,n,nr)
G state forward and backward transition rates
for use in transition rate matrices of Master equation
(different from gamma used in Gillespie algorithms)
"""
function get_gamma(r,n)
    gammaf = zeros(n+2)
    gammab = zeros(n+2)
    for i = 1:n
        gammaf[i+1] = r[2*(i-1)+1]
        gammab[i+1] = r[2*i]
    end
    return gammaf, gammab
end
"""
get_nu(r,n,nr)
R step forward transition rates
"""
function get_nu(r,n,nr)
    r[2*n+1 : 2*n+nr+1]
end
"""
get_eta(r,n,nr)
Intron ejection rates at each R step
"""
function get_eta(r,n,nr)
    eta = zeros(nr)
    if length(r) > 2*n + 2*nr
        eta[1] = r[2*n + 1 + nr + 1]
        for i = 2:nr
            eta[i] = eta[i-1] + r[2*n + 1 + nr + i]
        end
    end
    return eta
end
"""
survival_fraction(nu,eta,nr)
Fraction of introns that are not spliced prior to ejection
"""
function survival_fraction(nu,eta,nr)
    pd = 1.
    for i = 1:nr
        pd *= nu[i+1]/(nu[i+1]+eta[i])
    end
    return pd
end

function reparametrize(r,n)
    gammaf,gammb = get_gamma(r,n)
    alpha = Vector{Float64}(undef,n)
    alpha[1] = 1
    for i in 2:n
        alpha[i] = alpha[i-1] * gammaf[i]/gammab[i]
    end
    return alpha
end

function invert_parameter(alpha,n)

end

"""
normalized_nullspace(M::AbstractMatrix)
Compute the normalized null space of a nxn matrix
of rank n-1 using QR decomposition with pivoting
"""
function normalized_nullspace(M::AbstractMatrix)
    F = qr(M,Val(true));  #QR decomposition with pivoting
    R = F.R
    m = size(M,1)
    # if rank(M) == m-1
    p = zeros(m)
    # Back substitution to solve R*p = 0
    p[end] = 1.
    for i in 1:m-1
        p[m-i] = -R[m-i,m-i+1:end]'*p[m-i+1:end]/R[m-i,m-i]
    end
    p = max.(p,eps(Float64))
    p /= sum(p);
    return F.P*p
    # else
    #     return zeros(m)
    # end
end
"""
    allele_convolve(mhist,nalleles)

    Convolve to compute distribution for contributions from multiple alleles
"""
function allele_convolve(mhist,nalleles)
    nhist = length(mhist)
    mhists = Array{Array{Float64,1}}(undef,nalleles)
    mhists[1] = float.(mhist)
    for i = 2:nalleles
        mhists[i] = zeros(nhist)
        for m = 0:nhist-1
            for m2 = 0:min(nhist-1,m)
                mhists[i][m+1] += mhists[i-1][m-m2+1]*mhist[m2+1]
            end
        end
    end
    return mhists[nalleles]
end
"""
    allele_deconvolve(mhist,nalleles)

    Deconvolve to compute distribution of one allele from contributions of multiple alleles
"""
allele_deconvolve(mhist,nalleles) = irfft((rfft(mhist)).^(1/nalleles),length(mhist))

"""
init_prob(pss,n,nr)
Initial condition for first passage time calculation of active (ON) state
(Sum over all states with first R step occupied weighted by equilibrium state
of all states with all R steps not occupied with reporters)
"""
function init_prob(pss,n,nr)
    l = length(pss)
    base = findbase(l,n,nr)
    SAinit = zeros(l)
    if base == 3
        for z=1:2^(nr-1)
            aSAinit = (n+1) + (n+1)*decimal(vcat(2,digits(z-1,base=2,pad=nr-1)))
            apss = (n+1) + (n+1)*decimal(vcat(0,digits(z-1,base=2,pad=nr-1)))
            SAinit[aSAinit] = pss[apss]
        end
    else
        for z=1:2^(nr-1)
            aSAinit = (n+1) + (n+1)*decimal(vcat(1,zeros(Int,nr-1)))
            apss = n+1
            SAinit[aSAinit] = pss[apss]
        end

    end
    SAinit / sum(SAinit)
end
"""
init_prob(pss,n,nr)
Initial condition for first passage time calculation of inactive (OFF) state
(Sum over all states with no R step introns occupied weighted by equilibrium state
of all states with a single R step intron not)
"""
function init_prob(pss,r,n,nr,nonzeros)
    l = length(pss)
    base = findbase(l,n,nr)
    SIinit = zeros(l)
    nu = get_nu(r,n,nr)
    eta = get_eta(r,n,nr)
    if base == 3
        # Start of OFF state by ejection
        for i in 1:n+1, z in 1:2^(nr-1)
            exons = digits(z-1,base=2,pad=nr-1)
            ainit = i + (n+1)*decimal(vcat(exons,0))
            apss = i + (n+1)*decimal(vcat(exons,2))
            SIinit[ainit] += pss[apss]*nu[nr+1]
        end
        # Start of OFF state by splicing
        for i in 1:n+1, z in 1:2^(nr)
            exons = digits(z-1,base=2,pad=nr)
            intronindex = findall(exons.==1)
            ainit = i + (n+1)*decimal(exons)
            for j in intronindex
                introns = copy(exons)
                introns[j] = 2
                apss  = i + (n+1)*decimal(introns)
                SIinit[ainit] += pss[apss]*eta[j]
            end
        end
    else
        # Start of OFF state by ejection
        for i in 1:n+1
            ainit = i
            apss = i + (n+1)*decimal(vcat(zeros(Int,nr-1),1),2)
            SIinit[ainit] += pss[apss]*nu[nr+1]
        end
        # Start of OFF state by splicing
        for i in 1:n+1
            ainit = i
            for j in 1:nr
                introns = zeros(Int,nr)
                introns[j] = 1
                apss  = i + (n+1)*decimal(introns,2)
                SIinit[ainit] += pss[apss]*eta[j]
            end
        end

    end
    SIinit = SIinit[nonzeros]
    SIinit/sum(SIinit)
end

"""
accumulate(SA::Matrix,n,nr)
Sum over all probability vectors accumulated into OFF states
"""
function accumulate(SA,n,nr)
    l,p = size(SA)
    base = findbase(p,n,nr)
    SAj = zeros(size(SA)[1])
    for i=1:n+1, z=1:base^nr
        zdigits = digits(z-1,base=base,pad=nr)
        if ~any(zdigits.> base-2)
            a = i + (n+1)*(z-1)
            SAj += SA[:,a]
        end
    end
    return SAj
end
"""
accumulate(SI::Matrix,n,nr,nonzeros)
Sum over all probability vectors accumulated into ON states
"""
function accumulate(SI,n,nr,nonzeros,p)
    # Sum over all probability vectors accumulated into ON states
    l = size(SI)[1]
    base = findbase(p,n,nr)
    SIj = zeros(size(SI)[1])
    for i=1:n+1, z=1:base^nr
        zdigits = digits(z-1,base=base,pad=nr)
        if any(zdigits.> base-2)
            a = i + (n+1)*(z-1)
            if a in nonzeros
                SIj += SI[:,findfirst(a .== nonzeros)]
            end
        end
    end
    return SIj
end
"""
marginalize(p::Vector,n,nhist)
marginalize(p::Vector,n)
marginalize(P::Matrix)

Marginalize over G states
"""
function marginalize(p::Vector,n,nhist)
    mhist = zeros(nhist)
    nT = n+1
    for m in 1:nhist
        i = (m-1)*nT
        mhist[m] = sum(p[i+1:i+nT])
    end
    return mhist
end

function marginalize(p::Vector,n)
    nhist = Int(length(p)/(n+1)) - 2
    marginalize(p,n,nhist)
end

marginalize(P::Matrix) = sum(P,dims=1)

unfold(P::Matrix) = reshape(P,length(P))

max_stars_repo_name: nih-niddk-mbs/Transcription

max_issues_repo_name: nih-niddk-mbs/Transcription

max_forks_repo_name: nih-niddk-mbs/Transcription


text:
export iauEra00
"""
Earth rotation angle (IAU 2000 model).

This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

Status:  canonical model.

Given:
   dj1,dj2   double    UT1 as a 2-part Julian Date (see note)

Returned (function value):
            double    Earth rotation angle (radians), range 0-2pi

Notes:

   1. The UT1 date dj1+dj2 is a Julian Date, apportioned in any
   convenient way between the arguments dj1 and dj2.  For example,
   JD(UT1)=2450123.7 could be expressed in any of these ways,
   among others:

            dj1            dj2

      2450123.7           0.0       (JD method)
      2451545.0       -1421.3       (J2000 method)
      2400000.5       50123.2       (MJD method)
      2450123.5           0.2       (date & time method)

   The JD method is the most natural and convenient to use in
   cases where the loss of several decimal digits of resolution
   is acceptable.  The J2000 and MJD methods are good compromises
   between resolution and convenience.  The date & time method is
   best matched to the algorithm used:  maximum precision is
   delivered when the dj1 argument is for 0hrs UT1 on the day in
   question and the dj2 argument lies in the range 0 to 1, or vice
   versa.

   2. The algorithm is adapted from Expression 22 of Capitaine et al.
   2000.  The time argument has been expressed in days directly,
   and, to retain precision, integer contributions have been
   eliminated.  The same formulation is given in IERS Conventions
   (2003), Chap. 5, Eq. 14.

Called:
   iauAnp       normalize angle into range 0 to 2pi

References:

   Capitaine N., Guinot B. and McCarthy D.D, 2000, Astron.
   Astrophys., 355, 398-405.

   McCarthy, D. D., Petit, G. (eds.), IERS Conventions (2003),
   IERS Technical Note No. 32, BKG (2004)

This revision:  2013 June 18

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# double iauEra00(double dj1, double dj2)
function iauEra00(dj1::Real, dj2::Real)
   return ccall((:iauEra00, libsofa_c), Cdouble,
                   (Cdouble, Cdouble),
                   convert(Float64, dj1), convert(Float64, dj2))
end

ext: jl


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

function ksstats(
        x::AbstractVector{T},
        y::AbstractVector{S},
        w_x::AbstractVector{D},
        w_y::AbstractVector{M}
    ) where {T<:Real, S<:Real, D<:Real, M<:Real}

    sort_idx = sortperm([x; y])
    pdf_diffs = [w_x/sum(w_x); -w_y/sum(w_y)][sort_idx]
    cdf_diffs = cumsum(pdf_diffs)
    δp = maximum(cdf_diffs)
    δn = -minimum(cdf_diffs)
    δ = max(δp, δn)
    (δ, δp, δn)
end

"""
    bat_compare(
        samples_1::DensitySampleVector,
        samples_2::DensitySampleVector;
        nsamples::Symbol=:effective
    )

Compares two `DensitySampleVector`s given by `samples_1` and `samples_2` applying the Kolmogorov-Smirnov test for all marginals.

`nsamples` specifies how to define a number of samples in the Kolmogorov-Smirnov
distribution. The default value is `nsamples=:effective`, which uses the
effective number of samples estimated by `bat_eff_sample_size`. The optimal keywords:

* `:length`  — length of the `DensitySamplesVector` is used

* `:weights` — the sum of the weights is used

Returns a NamedTuple of the shape

```julia
(result = X::TypedTables.Table, ...)
```
"""
function bat_compare(samples_1::DensitySampleVector, samples_2::DensitySampleVector;  nsamples::Symbol=:effective, nested_eff_ss=:0)

    n_params_1 = size(flatview(unshaped.(samples_1.v)))[1]
    n_params_2 = size(flatview(unshaped.(samples_2.v)))[1]

    if n_params_1 !== n_params_2
        throw(ArgumentError("Samples with different number of parameters cannot be compared"))
    end

    samples_1_flat = flatview(unshaped.(samples_1.v))
    samples_2_flat = flatview(unshaped.(samples_2.v))

    if nsamples == :effective
        @info "Using effective number of samples in the Kolmogorov–Smirnov distribution"
        smpl_size_1 = round.(Int64, bat_eff_sample_size(unshaped.(samples_1)).result)
        smpl_size_2 = round.(Int64, bat_eff_sample_size(unshaped.(samples_2)).result)
    elseif nsamples == :nested
        @info "Using effective number of samples of nested sampling algorithm"
        if (nested_eff_ss==0) throw(ArgumentError("You have to give nested_eff_ss for bat_compare")) end
        smpl_size_1 = round.(Int64, bat_eff_sample_size(unshaped.(samples_1)).result)
        smpl_size_2_pK = round.(Int64, nested_eff_ss/length(smpl_size_1))
        smpl_size_2 = []
        for i in 1:length(smpl_size_1)
            append!(smpl_size_2,smpl_size_2_pK)
        end
    elseif nsamples == :length
        @info "Using lengths of the vectors in the Kolmogorov–Smirnov distribution"
        length_1 = length(samples_1)
        length_2 = length(samples_1)
        smpl_size_1 = repeat([length_1], n_params_1)
        smpl_size_2 = repeat([length_2], n_params_2)
    elseif nsamples == :weights
        @info "Using sum of the weights in the Kolmogorov–Smirnov distribution"
        tot_weights_1 = sum(samples_1.weight)
        tot_weights_2 = sum(samples_2.weight)
        smpl_size_1 = repeat([tot_weights_1], n_params_1)
        smpl_size_2 = repeat([tot_weights_2], n_params_2)
    else
        throw(ArgumentError("Unknown keyword"))
    end

    p_values_array = map(i-> begin
        test_result = ksstats( samples_1_flat[i, :], samples_2_flat[i, :], samples_1.weight, samples_2.weight)
        ks_test = HypothesisTests.ApproximateTwoSampleKSTest(smpl_size_1[i], smpl_size_2[i], test_result...)
        HypothesisTests.pvalue(ks_test)
    end, Base.OneTo(n_params_1))

    table_result = TypedTables.Table(
            dim_ind = Base.OneTo(n_params_1), # indices of parameters
            nsamples_1 = smpl_size_1, # number of samples for `samples_1`
            nsamples_2 = smpl_size_2, # number of samples for `samples_2`
            ks_p_values = p_values_array, # Kolmogorov-Smirnov p-values
            # additional tests can be added to this table/function later
    )

    return (result=table_result, )
end

export bat_compare


text:
"""
`module  JAC.Continuum`
    ... a submodel of JAC that contains all functions and methods for computing, either in a given potential or with
        regard to a well-defined bound state.
"""
module Continuum

    using  GSL, Printf, SpecialFunctions
    using  ..Basics, ..Bsplines, ..Defaults, ..ManyElectron, ..Radial, ..Nuclear

    """
    `struct  Continuum.Settings`  ... defines a type for the parameters for computing continuum orbitals.

        + includeExchange         ::Bool    ... True, if the exchange is to be included and false otherwise.
        + mtp                     ::Int64   ... No of grid points for which the continuum orbital(s) are to be computed.
    """
    struct Settings
        includeExchange           ::Bool
        mtp                       ::Int64
    end

    """
    `Continuum.generateOrbitalForLevel(energy::Float64, sh::Subshell, level::Level, nm::Nuclear.Model, grid::Radial.Grid, basis::Basis,
                                       settings::Continuum.Settings)`
        ... to generate a continuum orbital for the (continuum) subshell sh, the energy and the effective charge within the given
            potential. The continuum orbital is generated orthogonal with regard to all subshells of the same symmetry in the basis.
            All further specifications about this generations are made by proper settings. A tupel of a (continuum)
            (orbital::Orbital, phase::Float64) is returned.
    """
    function generateOrbitalForLevel(energy::Float64, sh::Subshell, level::Level, nm::Nuclear.Model, grid::Radial.Grid, settings::Continuum.Settings)
        #
        # Generate a (local) potential for the given level
        nuclearPotential  = Nuclear.nuclearPotential(nm, grid)
        ## wp1 = compute("radial potential: core-Hartree", grid, wLevel)
        ## wp2 = compute("radial potential: Hartree-Slater", grid, wLevel)
        ## wp3 = compute("radial potential: Kohn-Sham", grid, wLevel)
        ## wp4 = compute("radial potential: Dirac-Fock-Slater", grid, wLevel)
        wp = compute("radial potential: Kohn-Sham", grid, level)
        pot = Basics.add(nuclearPotential, wp)
        Defaults.warn(AddWarning(), "All continuum orbitals are generated in a local (DFS) potential.")

        # Generate a continuum orbital due to the given solution method
        if      Defaults.GBL_CONT_SOLUTION  ==  ContBessel()
            cOrbital = Continuum.generateOrbitalBessel(energy, sh, grid::Radial.Grid, settings)
        elseif  Defaults.GBL_CONT_SOLUTION  ==  AsymptoticCoulomb()
            cOrbital = Continuum.generateOrbitalAsymptoticCoulomb(energy, sh, pot, settings)
        elseif  Defaults.GBL_CONT_SOLUTION  ==  NonrelativisticCoulomb()
            cOrbital = Continuum.generateOrbitalNonrelativisticCoulomb(energy, sh, pot.grid, settings)
        elseif  Defaults.GBL_CONT_SOLUTION  ==  BsplineGalerkin()
            cOrbital = Continuum.generateOrbitalGalerkin(energy, sh, pot, settings)
        else    error("stop a")
        end
        #
        #
        # Normalize the continuum orbital and determine its phase
        if      Defaults.GBL_CONT_NORMALIZATION  ==  PureSineNorm()
            cOrbital, phase = Continuum.normalizeOrbitalPureSine(cOrbital, pot.grid, settings)
        elseif  Defaults.GBL_CONT_NORMALIZATION  ==  CoulombSineNorm()
            cOrbital, phase = Continuum.normalizeOrbitalCoulombSine(cOrbital, pot, settings)
        elseif  Defaults.GBL_CONT_NORMALIZATION  ==  OngRussekNorm()
            cOrbital, phase = Continuum.normalizeOrbitalOngRussek(cOrbital, pot, settings)
        else    error("stop b")
        end

        return( cOrbital, phase )
    end

    """
    `Continuum.generateOrbitalLocalPotential(energy::Float64, sh::Subshell, pot::Radial.Potential, settings::Continuum.Settings)`
        ... to generate a continuum orbital for the (continuum) subshell sh, the energy and the given potential. The effective
            charge for the normalization of the continuum orbital is derived from the potential. All further specifications about
            this generations are made by proper settings; however, the function termintates if the settings.includeExchange = true.
            A tupel of a (continuum) (orbital::Orbital, phase::Float64, normFactor::Float64) is returned.
    """
    function generateOrbitalLocalPotential(energy::Float64, sh::Subshell, pot::Radial.Potential, settings::Continuum.Settings)
        settings.includeExchange   &&   error("Continuum orbital for local potential does not allow 'exchange'.")

        # Generate a continuum orbital due to the given solution method
        if      Defaults.GBL_CONT_SOLUTION  ==  ContBessel()
            cOrbital = Continuum.generateOrbitalBessel(energy, sh, pot.grid, settings)
        elseif  Defaults.GBL_CONT_SOLUTION  ==  ContSine()
            cOrbital = Continuum.generateOrbitalPureSine(energy, sh, pot.grid, settings)
        elseif  Defaults.GBL_CONT_SOLUTION  ==  AsymptoticCoulomb()
            cOrbital = Continuum.generateOrbitalAsymptoticCoulomb(energy, sh, pot, settings)
        elseif  Defaults.GBL_CONT_SOLUTION  ==  NonrelativisticCoulomb()
            cOrbital = Continuum.generateOrbitalNonrelativisticCoulomb(energy, sh, pot.grid, settings)
        elseif  Defaults.GBL_CONT_SOLUTION  ==  BsplineGalerkin()
            cOrbital = Continuum.generateOrbitalGalerkin(energy, sh, pot, settings)
        else    error("stop a")
        end
        #
        # Normalize the continuum orbital and determine its phase
        if      Defaults.GBL_CONT_NORMALIZATION  ==  PureSineNorm()
            cOrbital, phase, normFactor = Continuum.normalizeOrbitalPureSine(cOrbital, pot.grid, settings)
        elseif  Defaults.GBL_CONT_NORMALIZATION  ==  CoulombSineNorm()
            cOrbital, phase, normFactor = Continuum.normalizeOrbitalCoulombSine(cOrbital, pot, settings)
        elseif  Defaults.GBL_CONT_NORMALIZATION  ==  OngRussekNorm()
            cOrbital, phase, normFactor = Continuum.normalizeOrbitalOngRussek(cOrbital, pot, settings)
        else    error("stop b")
        end
        #
        return( cOrbital, phase, normFactor )
    end

    """
    `Continuum.generateOrbitalAsymptoticCoulomb(energy::Float64, sh::Subshell, pot::Radial.Potential, settings::Continuum.Settings)`
        ... to generate a simple relativistic and non-normalized continuum orbital but by assuming an asymptotic sin/cos behaviour
            for both components at all mesh points. For this behaviour, both components and their derivatives can be calculated
            analytically. The effective charge Zbar is determined for the given potential and the non-Coulombic phase shifts is set
            to zero. An non-normalized cOrbital::Orbital is returned.

            Warning: At present, only the real part eta.re is taken into account; check for consistency !!
    """
    function generateOrbitalAsymptoticCoulomb(energy::Float64, sh::Subshell, pot::Radial.Potential, settings::Continuum.Settings)
        P = zeros(settings.mtp);    Q = zeros(settings.mtp);   Pprime = zeros(settings.mtp);    Qprime = zeros(settings.mtp)

        phi0 = 0.3
        wc   = Defaults.getDefaults("speed of light: c");    q = sqrt( 2*energy + energy/(wc*wc));    kappa = sh.kappa
        Zbar = Radial.determineZbar(pot)
        y    = Zbar * (energy + wc^2) / (wc^2 * q);   gammaBar = sqrt(kappa^2 - Zbar^2 / wc^2)
        eta  = - (kappa - im*y) / (energy + wc^2) / (gammaBar + im*y) / (2im)
        NP   = sqrt( (energy + 2*wc^2)/(pi*wc^2*q) );   NQ = -sqrt( energy / (pi*wc^2*q) )

        ## println("Continuum-asymptotic Coulomb orbital for q = $q,  kappa = $kappa,  Zbar=$Zbar  y = $y,   gammaBar = $gammaBar,  eta = $eta")
        if  Defaults.GBL_PRINT_DEBUG  println("Continuum-asymptotic Coulomb orbital for q = $q,  kappa = $kappa,  Zbar=$Zbar,  eta = $eta")    end

        for  i = 2:settings.mtp
            theta = q * pot.grid.r[i]  +  y * log(2q * pot.grid.r[i])  - angle( SpecialFunctions.gamma(gammaBar + im*y) )  -  pi*gammaBar/2  +  eta.re
            P[i] = NP * cos(theta + phi0);      Pprime[i] = - NP * (q + y/pot.grid.r[i]) * sin(theta + phi0)
            Q[i] = - NQ * sin(theta + phi0);    Qprime[i] = - NP * (q + y/pot.grid.r[i]) * cos(theta + phi0)
        end

        cOrbital = Orbital( sh, false, true, energy, P, Q, Pprime, Qprime, Radial.Grid())

        return( cOrbital )
    end

    """
    `Continuum.generateOrbitalBessel(energy::Float64, sh::Subshell, grid::Radial.Grid, settings::Continuum.Settings)`
        ... to generate a simple, non-relativistic and non-normalized Bessel wave for the large component of the continuum
            orbital and a small component due to the kinetic-balance condition. A (non-normalized) orbital::Orbital is returned.
            While Pprime is obtained from the analytical expression of the large component, Qprime is set to zero here.
    """
    function generateOrbitalBessel(energy::Float64, sh::Subshell, grid::Radial.Grid, settings::Continuum.Settings)
        P = zeros(settings.mtp);    Q = zeros(settings.mtp);   Pprime = zeros(settings.mtp);    Qprime = zeros(settings.mtp)

        q = sqrt( 2*energy );    l = Basics.subshell_l(sh);   wc = Defaults.getDefaults("speed of light: c")
        for  i = 2:settings.mtp   P[i] = GSL.sf_bessel_jl(l, q * grid.r[i]) * grid.r[i]    end
        for  i = 2:settings.mtp   qz = q * grid.r[i];   Pprime[i] = (l/qz * GSL.sf_bessel_jl(l, qz) - GSL.sf_bessel_jl(l+1, qz)) + P[i]/grid.r[i]     end
        for  i = 2:settings.mtp   Q[i] = -1/(2 * wc) * Pprime[i]  +  sh.kappa/grid.r[i] * P[i]         end

        println("Continuum-Bessel orbital for q = $q,  l = $l")

        cOrbital = Orbital( sh, false, false, energy, P, Q, Pprime, Qprime, grid)

        return( cOrbital )
    end

    """
    `Continuum.generateOrbitalGalerkin(energy::Float64, sh::Subshell, pot::Radial.Potential, settings::Continuum.Settings)`
        ... to generate a non-normalized continuum orbital within the given local potential by using the Galerkin method and a
            given B-spline basis. A (non-normalized) orbital::Orbital is returned.
    """
    function generateOrbitalGalerkin(energy::Float64, sh::Subshell, pot::Radial.Potential, settings::Continuum.Settings)
        P = zeros(settings.mtp);   Q = zeros(settings.mtp);   Pprime = zeros(settings.mtp);    Qprime = zeros(settings.mtp)
        nsL = pot.grid.nsL - 1;    nsS = pot.grid.nsS - 1
        wa = Bsplines.generatePrimitives(pot.grid)
        wb = Bsplines.generateGalerkinMatrix(wa, nsL, nsS, energy, sh, pot)
        wc = adjoint(wb) * wb

        # Test for 'real-symmetric matrix' ... and symmetrize otherwise
        nx = 0
        for  i = 1:nsL+nsS
            for  j = i+1:nsL+nsS
                if  abs(  (wc[i,j] - wc[j,i])/(wc[i,j] + wc[j,i]) ) > 1.0e-7   nx = nx + 1
                    @show "generateOrbitalGalerkin", i, j, wc[i,j], wc[j,i]
                end
            end
        end
        wd = Basics.diagonalize("matrix: LinearAlgebra", wc) ## , range=1:1)
        ## println(">>> Galerkin-eigenvalues = $(wd.values[1]), $(wd.values[2]) for  $sh  with  energy = $energy")

        cOrbital = Bsplines.generateOrbitalFromPrimitives(energy, sh, settings.mtp, wd.vectors[1], wa)
        mtp = size(cOrbital.P,1)
        println(">> Continuum B-spline-Galerkin orbital for energy=" * @sprintf("%.4e",energy) * ",  kappa=$(sh.kappa) " *
                "[mpt=$mtp, r[mtp]=" * @sprintf("%.4e",pot.grid.r[mtp]) * ", smallest eigenvalue=" * @sprintf("%.4e",wd.values[1]) * "].")

        return( cOrbital )
    end

    """
    `Continuum.generateOrbitalNonrelativisticCoulomb(energy::Float64, sh::Subshell, Zeff::Float64, grid::Radial.Grid,
                                                     settings::Continuum.Settings)`
        ... to generate a simple, non-relativistic and non-normalized Coulomb wave for the (continuum) subshell sh,
            the small component is simply obtained by the kinetic-balance condition. A orbital::Orbital is returned.
            *****
            This function does not yet work because there is no Julia implementation of the hypergeometric function with
            complex arguments available.
    """
    function generateOrbitalNonrelativisticCoulomb(energy::Float64, sh::Subshell, Zeff::Float64, grid::Radial.Grid, settings::Continuum.Settings)
        println("\n error: presently, NO Julia implementation of the hypergeometric function with complex arguments available.")

        # Define the derivative of the confluent hypergeometric function
        function hgPrime(a,b,z)    wa = a/z * (GSL.hypergeom(a+1,b,z) - GSL.hypergeom(a,b,z));   return(wa)   end

        P = zeros(settings.mtp);    Q = zeros(settings.mtp);   Pprime = zeros(settings.mtp);    Qprime = zeros(settings.mtp)

        q  = sqrt( 2*energy );    l = Basics.subshell_l(sh);   wc = Defaults.getDefaults("speed of light: c");   np = Zeff/q;   ws = 1
        wa = 2 * sqrt(Zeff) / sqrt(1 - Base.MathConstants.e^(-2pi*np)) / factorial(2l + 1)
        for  s = 1:l    ws = ws * sqrt(s^2 + np^2)    end

        for  i = 2:settings.mtp   P[i] = wa * ws * (2q*grid.r[i])^l * Base.MathConstants.e^(-im*q*grid.r[i]) *
                                         GSL.hypergeom(l+1+im*np, 2l+2., 2im*q*grid.r[i])                     end
        for  i = 2:settings.mtp
            Pprime[i] = 2q*l * (2q*grid.r[i])^(l-1) * Base.MathConstants.e^(-im*q*grid.r[i]) * GSL.hypergeom(l+1+im*np, 2l+2., 2im*q*grid.r[i]) -
                        im*q * (2q*grid.r[i])^l * Base.MathConstants.e^(-im*q*grid.r[i]) * GSL.hypergeom(l+1+im*np, 2l+2., 2im*q*grid.r[i])     +
                        (2q*grid.r[i])^l * Base.MathConstants.e^(-im*q*grid.r[i]) * hgPrime(l+1+im*np, 2l+2., 2im*q*grid.r[i])
            Pprime[i] = wa * ws * Pprime[i]
        end
        for  i = 2:settings.mtp   Q[i] = -1/(2 * wc) * Pprime[i]  +  sh.kappa/grid.r[i] * P[i]                end

        cOrbital = Orbital( sh, false, true, energy, P, Q, Pprime, Qprime, Radial.Grid())

        return( cOrbital )
    end

    """
    + `(energy::Float64, sh::Subshell, pot::Radial.Potential, settings::Continuum.Settings)`
        ... to generate a simple, non-relativistic and non-normalized Coulomb wave for the (continuum) subshell sh,
            the small component is simply obtained by the kinetic-balance condition. A orbital::Orbital is returned.
    """
    function generateOrbitalNonrelativisticCoulomb(energy::Float64, sh::Subshell, pot::Radial.Potential, settings::Continuum.Settings)
        Zbar     = Radial.determineZbar(pot)
        cOrbital = Continuum.generateOrbitalNonrelativisticCoulomb(energy, sh, Zbar, pot.grid, settings)

        return( cOrbital )
    end

    """
    `Continuum.generateOrbitalPureSine(energy::Float64, sh::Subshell, grid::Radial.Grid, settings::Continuum.Settings)`
        ... to generate a simple, non-relativistic and non-normalized Bessel wave for the large component of the continuum
            orbital and a small component due to the kinetic-balance condition. A (non-normalized) orbital::Orbital is returned.
            While Pprime is obtained from the analytical expression of the large component, Qprime is set to zero here.
    """
    function generateOrbitalPureSine(energy::Float64, sh::Subshell, grid::Radial.Grid, settings::Continuum.Settings)
        P = zeros(settings.mtp);    Q = zeros(settings.mtp);   Pprime = zeros(settings.mtp);    Qprime = zeros(settings.mtp)

        phi0 = 0.0
        q = sqrt( 2*energy );    l = Basics.subshell_l(sh);   wc = Defaults.getDefaults("speed of light: c");    N = sqrt(2/(pi*q))
        for  i = 1:settings.mtp   P[i]      = N * sin(q * grid.r[i] - l*pi/2 + phi0)                  end
        for  i = 1:settings.mtp   Pprime[i] = N * q * cos(q * grid.r[i] - l*pi/2 + phi0)              end
        for  i = 2:settings.mtp   Q[i] = -1/(2 * wc) * Pprime[i]  +  sh.kappa/grid.r[i] * P[i]        end

        println(">> Continuum-pure sine orbital for q = $q,  l = $l")

        cOrbital = Orbital( sh, false, false, energy, P, Q, Pprime, Qprime, grid)

        return( cOrbital )
    end

    """
    `Continuum.gridConsistency(maxEnergy::Float64, grid::Radial.Grid)`
        ... to check the consistency of the given grid with the maximum energy of the required continuum electrons;
            an error message is issued if the grid.hp = 0.   or  15 * grid.hp < wavelength(maxEnergy)   or  if the grid has
            less than 600 grid points. The function also returns the recommended grid point where the normalization and phase
            is to be determined. This number is currently set to nrContinuum = grid.NoPoints - 200  ... to correct for the wrong
            'phase behaviour' at large r-values.
    """
    function gridConsistency(maxEnergy::Float64, grid::Radial.Grid)
        Defaults.warn(AddWarning(), "Continuum.gridConsistency(): Improve the grid point for normalization; currently 600.")

        wavenb      = sqrt( 2maxEnergy + maxEnergy * Defaults.getDefaults("alpha")^2 )
        wavelgth    = 2pi / wavenb
        nrContinuum = grid.NoPoints - 200

        if      grid.hp == 0.               error("Improper grid for continuum processes with grid.hp = 0.")
        elseif  15 * grid.hp > wavelgth     error("Improper grid for continuum processes with 15*grid.hp = $(15*grid.hp) > " *
                                                  "wavelength = $wavelgth .")
        elseif  grid.NoPoints < 600               error("Improper No of grid points for continuum processes; grid.NoPoints = $(grid.NoPoints).")
        elseif  grid.r[nrContinuum] < 2.0   error("Improper grid extent for continuum processes; grid.r[nrContinuum] = " *
                                                  "$(grid.r[nrContinuum]) < 2.")
        end

        return( nrContinuum )
    end

    """
    `Continuum.normalizeOrbitalPureSine(cOrbital::Orbital, grid::Radial.Grid, settings::Continuum.Settings)`
        ... to normalize the given continuum orbital with regard to a (asymptotic) pure sine-function.
            An (on-energy-scale-normalized) orbital::Orbital and its relative phase phi w.r.t.  sin(kr + phi) is returned.
    """
    function normalizeOrbitalPureSine(cOrbital::Orbital, grid::Radial.Grid, settings::Continuum.Settings)
        mtp = size( cOrbital.P, 1);    q = sqrt( 2*cOrbital.energy );        l = Basics.subshell_l(cOrbital.subshell)

        PPprime = cOrbital.P[mtp] / cOrbital.Pprime[mtp];    kr  = q * grid.r[mtp];   at = atan( PPprime * q )
        phi = at - kr + l*pi/2;
        phi = rem(phi+1000pi, pi)    ##x to bring phi in the interval 0. <= phi < pi
        # For normalization, determine the maximum of the large component within the last 300 point
        A = maximum( abs.(cOrbital.P[end-300:end]) );    N   = sqrt(2/(pi*q)) / A

        println(">> Pure-sine normalized continuum orbital with normalization constant N=" * @sprintf("%.4e",N) *
                " and phase phi=" * @sprintf("%.4e",phi) * " at r=" * @sprintf("%.4e",grid.r[mtp])   * " a.u." )

        if abs(N) > 30.0   error("Pure-sine normalized N=" * @sprintf("%.4e",N) * "==> enlarge box-size")   end

        P = N .* cOrbital.P;   Q = N .* cOrbital.Q;   Pprime = N .* cOrbital.Pprime;   Qprime = N .* cOrbital.Qprime;
        newOrbital = Orbital( cOrbital.subshell, cOrbital.isBound, cOrbital.useStandardGrid, cOrbital.energy,
                              P, Q, Pprime, Qprime, cOrbital.grid)

        normF = maximum( abs.(P[end-300:end]) )   ## determine the maximum value of the P component as asymptotic measure for the amplitude

        return( newOrbital, phi, normF )
    end

    """
    `Continuum.normalizeOrbitalCoulombSine(cOrbital::Orbital, pot::Radial.Potential, settings::Continuum.Settings)`
        ... to normalize the given continuum orbital with regard to a (asymptotic) pure sine-function.
            An (on-energy-scale-normalized) orbital::Orbital is returned.
    """
    function normalizeOrbitalCoulombSine(cOrbital::Orbital, pot::Radial.Potential, settings::Continuum.Settings)
        mtp = size( cOrbital.P, 1);   energy = cOrbital.energy

        if  false  ## use non-relativistic Coulomb-sine wave representation
            wc   = Defaults.getDefaults("speed of light: c");    q = sqrt( 2*cOrbital.energy );        l = Basics.subshell_l(cOrbital.subshell)
            Zbar = Radial.determineZbar(pot);
            PPprime = cOrbital.P[mtp] / cOrbital.Pprime[mtp];    kr  = q * pot.grid.r[mtp];   at = atan( PPprime * q )
            sigmac  = SpecialFunctions.gamma(l + 1 - im*Zbar/q)
            phi = at - kr + l*pi/2 - Zbar/q * log(2*q*pot.grid.r[mtp]) - atan(sigmac.im / sigmac.re);
            phi = rem(phi+1000pi, pi)    ##x to bring phi in the interval 0. <= phi < pi
            # For normalization, determine the maximum of the large component within the last 300 point
            A = maximum( abs.(cOrbital.P[end-300:end]) );    N   = sqrt(2/(pi*q)) / A

            println(">> (Non-relativistic) Coulomb-sine normalized continuum orbital with normalization constant N=" * @sprintf("%.4e",N) *
                    " and phase phi=" * @sprintf("%.4e",phi) * " at r=" * @sprintf("%.4e",pot.grid.r[mtp])   * " a.u." )
            #
        else      ## use relativistic Coulomb-sine wave representation
            wc   = Defaults.getDefaults("speed of light: c");    q = sqrt( 2*energy + energy/(wc*wc));      kappa = cOrbital.subshell.kappa
            Zbar = Radial.determineZbar(pot);    y = Zbar * (energy + wc^2) / (wc^2 * q);   gammaBar = sqrt(kappa^2 - Zbar^2 / wc^2)
            eta  = log( - (kappa - im*y) / (energy + wc^2) / (gammaBar + im*y)) / (2*im)
            NP   = sqrt( (energy + 2*wc^2)/(pi*wc^2*q) );   NQ = -sqrt( energy / (pi*wc^2*q) )

            # Select the 'ratio' on which the thetaprimephase is determined
            if  false   PprimeP = cOrbital.Pprime[mtp] / cOrbital.P[mtp];        at = atan( - PprimeP / thetaprime );   son = "P'/P"
            else        PprimeP = cOrbital.Q[mtp] * NP / cOrbital.P[mtp] / NQ;   at = atan( PprimeP );                  son = "Q/P"     end
            #
            argG  = SpecialFunctions.gamma(gammaBar + im*y)
            theta = q * pot.grid.r[mtp]  +  y * log(2q * pot.grid.r[mtp])  - atan(argG.im / argG.re)  -  pi*gammaBar/2 ##  +  eta.re
            phi = at - theta
            phi = rem(phi+1000pi, pi)    ##x to bring phi in the interval 0. <= phi < pi
            ## A   = cOrbital.P[mtp] / cos(theta + phi);           N   = NP / A
            A = maximum( abs.(cOrbital.P[end-300:end]) );    N   = NP / A

            println(">> (Relativistic) Coulomb-sine normalized continuum orbital with normalization constant N=" * @sprintf("%.4e",N) *
                    " and phase phi=" * @sprintf("%.4e",phi) * " at r=" * @sprintf("%.4e",pot.grid.r[mtp])   * " a.u." )
        end

        if abs(N) > 30.0   error("Coulomb-sine normalized N=" * @sprintf("%.4e",N) * "==> enlarge box-size")   end

        P = N .* cOrbital.P;   Q = N .* cOrbital.Q;   Pprime = N .* cOrbital.Pprime;   Qprime = N .* cOrbital.Qprime;
        newOrbital = Orbital( cOrbital.subshell, cOrbital.isBound, cOrbital.useStandardGrid, cOrbital.energy,
                              P, Q, Pprime, Qprime, cOrbital.grid)

        normF = maximum( abs.(P[end-300:end]) )   ## determine the maximum value of the P component as asymptotic measure for the amplitude

        return( newOrbital, phi, normF )
    end

    """
    `Continuum.normalizeOrbitalOngRussek(cOrbital::Orbital, pot::Radial.Potential, settings::Continuum.Settings)`
        ... to normalize the given continuum orbital with regard to a (asymptotic) pure sine-function.
            An (on-energy-scale-normalized) orbital::Orbital is returned.
    """
    function normalizeOrbitalOngRussek(cOrbital::Orbital, pot::Radial.Potential, settings::Continuum.Settings)
        mtp = size( cOrbital.P, 1);                         energy = cOrbital.energy;      kappa = cOrbital.subshell.kappa
        wc   = Defaults.getDefaults("speed of light: c");   E = wc^2 + cOrbital.energy;    q = sqrt( 2*energy + energy/(wc*wc))
        qq = sqrt(E^2/wc^2 - wc^2)

        dVdr = pot.Zr[mtp] / (pot.grid.r[mtp]^2);    V  = - pot.Zr[mtp] / pot.grid.r[mtp]
        Pe = cOrbital.P[mtp];   Pep = cOrbital.Pprime[mtp];       r0 = pot.grid.r[mtp]
        # Compute U-function
        U = wc /2. * ((E-V)*dVdr - wc^2 *kappa*(kappa+1) /r0^3) / (( (E-V)^2 -wc^4 - wc^2 *kappa*(kappa+1) / r0^2)^(3/2)) * Pe
        ## U = U - wc * (Pep + dVdr * Pe / (2*(E-V+ 2* wc^2)) ) / (( (E-V)^2 -wc^4 - wc^2 *kappa*(kappa+1)/r0^2)^(1/2))
        U = U - wc * (Pep + dVdr * Pe / (2*(E-V+ 2* wc^2)) ) / (( (E-V)^2 -wc^4 - wc^2 *kappa*(kappa+1)/r0^2)^(1/2))
        # Compute A-function
        A = (( (E-V)^2 -wc^4 - wc^2 *kappa*(kappa+1)/r0^2)^(1/2)) / (E-V+wc^2) * (Pe^2 + U^2)
        A = sqrt(A);              phir0 = atan(U, Pe)
        phi = phir0 - q*r0;       N  = 1 / (A * sqrt(pi*wc))
        phi = rem(phi+1000pi, pi)    ##x to bring phi in the interval 0. <= phi < pi

        println(">> WKB (Ong-Russek) normalized continuum orbital with normalization constant N=" * @sprintf("%.4e",N) *
                " and phase phi=" * @sprintf("%.4e",phi) *
                " at r=" * @sprintf("%.4e",pot.grid.r[mtp])   * " a.u.  U/Pe = $(U/Pe)" )

        if abs(N) > 100.0   error("WKB (Ong-Russek) normalized N=" * @sprintf("%.4e",N) * "==> enlarge box-size")   end

        P = N .* cOrbital.P;   Q = N .* cOrbital.Q;   Pprime = N .* cOrbital.Pprime;   Qprime = N .* cOrbital.Qprime;
        newOrbital = Orbital( cOrbital.subshell, cOrbital.isBound, cOrbital.useStandardGrid, cOrbital.energy,
                              P, Q, Pprime, Qprime, cOrbital.grid)

        normF = maximum( abs.(P[end-300:end]) )   ## determine the maximum value of the P component as asymptotic measure for the amplitude

        return( newOrbital, phi, normF )
    end

end # module


text:
# This file was generated, do not modify it. # hide
lrm = LinMod()
lr1 = machine(lrm, Xhp, y) # poly of degree 1 (line)
fit!(lr1, rows=train)

lrm.fs.features = [:hp1, :hp2] # poly of degree 2
lr2 = machine(lrm, Xhp, y)
fit!(lr2, rows=train)

lrm.fs.features = [:hp1, :hp2, :hp3] # poly of degree 3
lr3 = machine(lrm, Xhp, y)
fit!(lr3, rows=train)


text:
@testset "IfThenElse" begin
    mock = MOIU.MockOptimizer(MOIU.UniversalFallback(MOIU.Model{Int}()))
    config = MOIT.Config()

    MOIU.set_mock_optimize!(
        mock,
        (mock::MOIU.MockOptimizer) -> (MOIU.mock_optimize!(mock, [1, 1, 0])),
        (mock::MOIU.MockOptimizer) -> (MOIU.mock_optimize!(mock, [1, 1, 0])),
    )
    COIDT.ifthenelsetest(mock, config)
end


text:
# ---
# title: 528. Random Pick with Weight
# id: problem528
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Binary Search, Random
# link: <https://leetcode.com/problems/random-pick-with-weight/description/>
# hidden: true
# ---
#
# You are given an array of positive integers `w` where `w[i]` describes the
# weight of `i``th` index (0-indexed).
#
# We need to call the function `pickIndex()` which **randomly** returns an
# integer in the range `[0, w.length - 1]`. `pickIndex()` should return the
# integer proportional to its weight in the `w` array. For example, for `w = [1,
# 3]`, the probability of picking the index `0` is `1 / (1 + 3) = 0.25` (i.e
# 25%) while the probability of picking the index `1` is `3 / (1 + 3) = 0.75`
# (i.e 75%).
#
# More formally, the probability of picking index `i` is `w[i] / sum(w)`.
#
#
#
# **Example 1:**
#
#
#
#     **Input**
#     ["Solution","pickIndex"]
#     [[[1]],[]]
#     **Output**
#     [null,0]
#
#     **Explanation**
#     Solution solution = new Solution([1]);
#     solution.pickIndex(); // return 0. Since there is only one single element on the array the only option is to return the first element.
#
#
# **Example 2:**
#
#
#
#     **Input**
#     ["Solution","pickIndex","pickIndex","pickIndex","pickIndex","pickIndex"]
#     [[[1,3]],[],[],[],[],[]]
#     **Output**
#     [null,1,1,1,1,0]
#
#     **Explanation**
#     Solution solution = new Solution([1, 3]);
#     solution.pickIndex(); // return 1. It's returning the second element (index = 1) that has probability of 3/4.
#     solution.pickIndex(); // return 1
#     solution.pickIndex(); // return 1
#     solution.pickIndex(); // return 1
#     solution.pickIndex(); // return 0. It's returning the first element (index = 0) that has probability of 1/4.
#
#     Since this is a randomization problem, multiple answers are allowed so the following outputs can be considered correct :
#     [null,1,1,1,1,0]
#     [null,1,1,1,1,1]
#     [null,1,1,1,0,0]
#     [null,1,1,1,0,1]
#     [null,1,0,1,0,0]
#     ......
#     and so on.
#
#
#
#
# **Constraints:**
#
#   * `1 <= w.length <= 10000`
#   * `1 <= w[i] <= 10^5`
#   * `pickIndex` will be called at most `10000` times.
#
#
## @lc code=start
using LeetCode

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


text:
function MCMCChainSummary(
    chains::AbstractArray{T,3}, model::AbstractProbabilityModel{P}; threaded::Bool = Threads.nthreads() > 1
) where {T,P}
    P2, S, C = size(chains)
    @assert P == P2 "The model is of a $P-dimensional parameter space, but we have $P2 parameters in our chain."
    D = DistributionParameters.constrained_length(model)
    # May as well allocate it on our stack, because we're returning a summary of this chain rather than the chain
    # this way, its at least recoverable
    stride_ch1 = stride(chains, 2)
    # stride_ch2 = PaddedMatrices.calc_padding(D, T)
    # chainarray = PaddedMatrices.DynamicPtrArray{T,3}(pointer(STACK_POINTER_REF[],T), (D, samples, nchains), stride_ch2)
    stride_ch2 = D
    chainarray = Array{T}(undef, D, S, C)
    ptr_ch1 = pointer(chains)
    ptr_ch2 = pointer(chainarray)
    for _ ∈ 1:C*S
        s1 = PaddedMatrices.DynamicPtrVector{T}(ptr_ch1, P, stride_ch1)
        s2 = PaddedMatrices.DynamicPtrVector{T}(ptr_ch2, D, stride_ch2)
        DistributionParameters.constrain!(s2, model, s1)
        ptr_ch1 += sizeof(T) * stride_ch1
        ptr_ch2 += sizeof(T) * stride_ch2
    end
    MCMCChainSummary(
        chainarray,
        DistributionParameters.parameter_names(model),
        threaded = threaded
    )
end

@generated function store_named_tuple!(v::AbstractArray{T}, nt::NT) where {T, NT <: NamedTuple}
    P = first(NT.parameters)
    q = quote ptrv = pointer(v) end
    for p ∈ P
        θ = gensym(:θ)
        push!(q.args, :($θ = nt.$p; PaddedMatrices.unique_copyto!(ptrv, $θ); ptrv += $(sizeof(T)) * type_length($θ)))
    end
    push!(q.args, nothing)
    q
end

# The plan here is to take a function as an argument
# this function is to take keyword arguments with the same name as model parameters -- capturing the unused in _..., eg
# `f(; a, b, _...) = a * b
# and then return a named tuple.
# Calls MCMCChainSummary on the named tuple.
function MCMCChainSummary(
    f::F, chains::AbstractArray{T,3}, model::AbstractProbabilityModel{P}; threaded::Bool = Threads.nthreads() > 1
) where {T,P,F}

    P2, S, C = size(chains)
    @assert P == P2 "The model is of a $P-dimensional parameter space, but we have $P2 parameters in our chain."
    stride_ch1 = stride(chains, 2)
    ptr_ch1 = pointer(chains)
    s1 = PaddedMatrices.DynamicPtrVector{T}(ptr_ch1, P, stride_ch1)
    fnt = f(; constrain(data, s1)...)
    D = PaddedMatrices.type_length(fnt)
    chainarray = Array{T}(undef, D, S, C)

    param_names = DistributionParameters.parameter_names(fnt)

    stride_ch2 = D
    ptr_ch2 = pointer(chainarray)
    s2 = PaddedMatrices.DynamicPtrVector{T}(ptr_ch2, D, stride_ch2)
    store_named_tuple!(s2, fnt)
    for _ ∈ 2:C*S
        ptr_ch1 += sizeof(T) * stride_ch1
        ptr_ch2 += sizeof(T) * stride_ch2

        s1 = PaddedMatrices.DynamicPtrVector{T}(ptr_ch1, P, stride_ch1)
        s2 = PaddedMatrices.DynamicPtrVector{T}(ptr_ch2, D, stride_ch2)
        fnt = f(; constrain(data, s1)...)
        store_named_tuple!(s2, fnt)
    end
    MCMCChainSummary(
        chainarray,
        param_names,
        threaded = threaded
    )
end


text:
val = ccall((:plus_func, "./libfunc"), Int64, (Int64, Int64), 1, 2)
println(val)

max_stars_repo_name: kirisaki/julia-ffi

max_issues_repo_name: kirisaki/julia-ffi

max_forks_repo_name: kirisaki/julia-ffi


text:
let
	a = VecE3()
	@test typeof(a) <: VecE3
	@test typeof(a) <: AbstractVec
	@test a.x == 0.0
	@test a.y == 0.0
	@test a.z == 0.0

	b = VecE3(0.0,0.0,0.0)
	@test b.x == 0.0
	@test b.y == 0.0
	@test b.z == 0.0

	@test length(a) == 3
	@test a == b
	@test isequal(a, b)
	@test copy(a) == a
	@test convert(Vector{Float64}, a) == [0.0,0.0,0.0]
	@test convert(VecE3, [0.0,0.0,0.0]) == a

	@test isapprox(polar(1.0,0.0,0.0), VecE3(1.0,0.0,0.0))
	@test isapprox(polar(2.0,π/2,1.0), VecE3(0.0,2.0,1.0))
	@test isapprox(polar(3.0,-π/2,-0.5), VecE3(0.0,-3.0,-0.5))
	@test isapprox(polar(-0.5,1.0*π,1.0*π), VecE3(0.5,0.0,1.0*π))

	a = VecE3(0.0,1.0, 2.0)
	b = VecE3(0.5,2.0,-3.0)

	@test a != b

	@test a + b == VecE3(0.5, 3.0, -1.0)
	@test a + 1 == VecE3(1.0, 2.0,  3.0)
	@test a + 0.5 == VecE3(0.5, 1.5, 2.5)

	@test a - b == VecE3(-0.5, -1.0, 5.0)
	@test a - 1 == VecE3(-1.0,  0.0, 1.0)
	@test a - 0.5 == VecE3(-0.5, 0.5, 1.5)

	@test a * 2 == VecE3(0.0, 2.0, 4.0)
	@test a * 0.5 == VecE3(0.0, 0.5, 1.0)

	@test a / 2 == VecE3(0.0, 0.5, 1.0)
	@test a / 0.5 == VecE3(0.0, 2.0, 4.0)

	@test b.^2 == VecE3(0.25, 4.0, 9.0)
	@test VecE3(0.5, 2.0, 3.0).^0.5 == VecE3(0.5^0.5, 2.0^0.5, 3.0^0.5)

	@test isfinite(a)
	@test !isfinite(VecE3(Inf,0,0))
	@test !isfinite(VecE3(0,Inf,0))
	@test !isfinite(VecE3(0,-Inf,0))
	@test !isfinite(VecE3(Inf,Inf,Inf))

	@test !isinf(a)
	@test isinf(VecE3(Inf,0,0))
	@test isinf(VecE3(0,Inf,0))
	@test isinf(VecE3(0,-Inf,0))
	@test isinf(VecE3(Inf,Inf,Inf))

	@test !isnan(a)
	@test isnan(VecE3(NaN,0,NaN))
	@test isnan(VecE3(0,NaN,NaN))
	@test isnan(VecE3(NaN,NaN,NaN))

	@test round.(VecE3(0.25,1.75,0)) == VecE3(0.0,2.0,0.0)
	@test round.(VecE3(-0.25,-1.75,0)) == VecE3(-0.0,-2.0,0.0)

	@test floor.(VecE3(0.25,1.75,0.0)) == VecE3(0.0,1.0,0.0)
	@test floor.(VecE3(-0.25,-1.75,0.0)) == VecE3(-1.0,-2.0,0.0)

	@test ceil.(VecE3(0.25,1.75,0.0)) == VecE3(1.0,2.0,0.0)
	@test ceil.(VecE3(-0.25,-1.75,0.0)) == VecE3(-0.0,-1.0,0.0)

	@test trunc.(VecE3(0.25,1.75,0.0)) == VecE3(0.0,1.0,0.0)
	@test trunc.(VecE3(-0.25,-1.75,0.0)) == VecE3(-0.0,-1.0,0.0)

	@test clamp.(VecE3(1.0, 10.0, 0.5), 0.0, 5.0) == VecE3(1.0, 5.0, 0.5)
	@test clamp.(VecE3(-1.0, 4.0, 5.5), 0.0, 5.0) == VecE3(0.0, 4.0, 5.0)

	c = VecE3(3.0,4.0,5.0)

	@test isapprox(abs.(a), a)
	@test isapprox(abs.(c), c)
	@test isapprox(norm(c), sqrt(3^2 + 4^2 + 5^2))
	@test isapprox(normalize(a), VecE3(0.0,1.0/hypot(1.0,2.0), 2.0/hypot(1.0,2.0)))

	@test isapprox(dist(a,a), 0.0)
	@test isapprox(dist(b,b), 0.0)
	@test isapprox(dist(c,c), 0.0)
	@test isapprox(dist(a,b), norm([0.5, 1.0, 5.0]))
	@test isapprox(dist(a,c), norm([3.0, 3.0, 3.0]))

	@test isapprox(dist2(a,a), 0.0)
	@test isapprox(dist2(b,b), 0.0)
	@test isapprox(dist2(c,c), 0.0)
	@test isapprox(dist2(a,b), 0.5^2 + 1^2 + 5^2)
	@test isapprox(dist2(a,c), 3^2 + 3^2 + 3^2)

	@test isapprox(dot(a, b), 2.0 - 6.0)
	@test isapprox(dot(b, c), 1.5 + 8.0 - 15.0)
	@test isapprox(dot(c, c), norm(c)^2)

	@test isapprox(proj(b, a, Float64), -norm(VecE3(0.0, -0.8, -1.6)))
	@test isapprox(proj(c, a, Float64),  norm(VecE3(0.0,  2.8,  5.6)))

	@test isapprox(proj(b, a, VecE3), VecE3(0.0, -0.8, -1.6))
	@test isapprox(proj(c, a, VecE3), VecE3(0.0,  2.8,  5.6))

	@test isapprox(lerp(VecE3(0,0,0), VecE3(2,-2,2), 0.25), VecE3(0.5,-0.5,0.5))
	@test isapprox(lerp(VecE3(2,-2,2), VecE3(3,-3,3), 0.5), VecE3(2.5,-2.5,2.5))

	@test isapprox(rot(VecE3(1,2,3), VecE3(2,0,0), 0.0), VecE3(1,2,3))
	@test isapprox(rot(VecE3(1,2,3), VecE3(2,0,0), π/2), VecE3(1,-3,2))
	@test isapprox(rot(VecE3(1,2,3), VecE3(2,0,0), -π/2), VecE3(1,3,-2))
	@test isapprox(rot(VecE3(1,2,3), VecE3(2,0,0), 1.0π), VecE3(1,-2,-3))
	@test isapprox(rot(VecE3(1,2,3), VecE3(0,2,0), π/2), VecE3(3,2,-1))
	@test isapprox(rot(VecE3(1,2,3), VecE3(0,0,2), π/2), VecE3(-2,1,3))

	@test isapprox(rot_normalized(VecE3(1,2,3), VecE3(1,0,0), 0.0), VecE3(1,2,3))
	@test isapprox(rot_normalized(VecE3(1,2,3), VecE3(1,0,0), π/2), VecE3(1,-3,2))
	@test isapprox(rot_normalized(VecE3(1,2,3), VecE3(1,0,0), -π/2), VecE3(1,3,-2))
	@test isapprox(rot_normalized(VecE3(1,2,3), VecE3(1,0,0), 1.0π), VecE3(1,-2,-3))
	@test isapprox(rot_normalized(VecE3(1,2,3), VecE3(0,1,0), π/2), VecE3(3,2,-1))
	@test isapprox(rot_normalized(VecE3(1,2,3), VecE3(0,0,1), π/2), VecE3(-2,1,3))
end


text:
#
# This is a part of Manager module. It represents Manager's Remote API, which
# may be called from remote client or other manager. For this @rpc annotation
# is used for functions. The value, which is returned by these functions will
# be transferred to remote client (caller) as an answer. If @rpc function
# returns nothing, then response to the remove caller will not be transferred.
#
# @author DeadbraiN
#
import RpcApi
import Config
import Mutator
import World
import Event
import FastApi
import Server
import Helper
import Connection
import ManagerTypes
import Dots
#
# Sides positions of near instance(server), which is connected
# to the current one.
#
const _SIDE_LEFT  = "LEFT"
const _SIDE_RIGHT = "RIGHT"
const _SIDE_UP    = "UP"
const _SIDE_DOWN  = "DOWN"
#
# @rpc
# Grabs world's rectangle region and returns it. This function is optimized
# for fast transferring through network by obtaining small "Region" type
# instance. It works in a little complicated way:
# - first it normalizes coordinates. It means, that thay should be in a range 1:width
#   and 1:height of world (see Config.worldWidth/worldHright configs).
# - second, it splits obtained region into small 16x16 squares (see _calcBlocksAmount).
#   If the region isn't multiple to 16, then additional squares will be at the
#   right and at the bottom.
# - Every dot within one square may be described as 1 byte index (0..255, see
#   Region.Block.energy). This is how we store all energy dots (see RpcApi.Block.energy)
#   and organism dots (see RpcApi.Org types) within squares by index
# - empty (black) dots are skipped
# - In one RpcApi.Org type we may transfer any organism related info. It's possible
#   to add more fields to this type in future.
#
# @param man Manager data type
# @param x Start X coordinate of region
# @param y Start Y coordinate of region
# @param x1 End x. 0 means all width
# @param y1 End y. 0 means all height
#
function getRegion(man::ManagerTypes.ManagerData, x::Int = 1, y::Int = 1, x1::Int = 0, y1::Int = 0)
  x, y, x1, y1 = _normalizeRegion(man, x, y, x1, y1)

  local data::Array{UInt16, 2} = _cutRegion(man, x, y, x1, y1)
  local xBlocks::Int = _calcBlocksAmount(x, x1)
  local yBlocks::Int = _calcBlocksAmount(y, y1)
  local region::RpcApi.Region = RpcApi.Region(
    UInt8(xBlocks),
    UInt8(yBlocks),
    UInt16(x1-x+1),
    UInt16(y1-y+1),
    RpcApi.Block[],
    man.ips
  )
  local i::Int
  local xOffs::Int = 1
  local yOffs::Int = 1
  local xBlock::Int
  local yBlock::Int
  local curX::UInt16
  local curY::UInt16
  local block::RpcApi.Block
  local org::Creature.Organism
  local offs::UInt8
  local pos::Helper.Point = Helper.Point(0, 0)
  #
  # Creates all blocks
  #
  for i = 1:(xBlocks * yBlocks) push!(region.blocks, RpcApi.Block(UInt8[], RpcApi.Org[])) end
  #
  # Walks through all dots in a world and stores only energy and organisms
  #
  for curY = UInt16(y):UInt16(y1)
    pos.y = curY
    for curX = UInt16(x):UInt16(x1)
      pos.x = curX
      if _isOrganism(man, Int(curX), Int(curY)) _addOrg(man, region, curX, curY)
      elseif World.getEnergy(man.world, pos) > Dots.INDEX_EMPTY _addEnergy(man, region, curX, curY)
      end
    end
  end

  region
end
#
# @rpc
# Creates tasks and organisms according to Config. All tasks
# will be in man.tasks field.
# @param man Manager data type
#
function createOrganisms(man::ManagerTypes.ManagerData)
  local i::Int

  if man.cfg.modeQuiet < Config.MODE_QUIET_NO Helper.info("Creating organisms...") end
  man.totalOrganisms      = 0
  # TODO: for system with many managers id should be unique between them
  man.organismId          = UInt(1)
  man.cfg.orgErrors = 0
  for i = 1:man.cfg.orgStartAmount createOrganism(man) end

  nothing
end
#
# @rpc
# Creates one task and organism inside this task. Created
# task will be added to man.tasks array. Position
# may be set or random free position will be used.
# @param man Manager data type
# @param pos Position|nothing Position of the organism
# @return {Int} Organism id or false if organisms limit is riched
#
function createOrganism(man::ManagerTypes.ManagerData, pos::Helper.Point = Helper.Point(0, 0))
  if ManagerTypes.orgAmount(man) > man.cfg.worldMaxOrgs return nothing end
  Manager._createOrganism(man, nothing, pos)

  nothing
end
#
# @rpc
# Sets configuration value according to it's section and a key
# @param man Manager data type
# @param Unique config symbol
# @param value Custom value
#
function setConfig(man::ManagerTypes.ManagerData, name::Symbol, value::Any)
  setfield!(man.cfg, name, value)
  nothing
end
#
# @rpc
# Gets configuration value according to it's section and a key
# @param man Manager data type
# @param name config symbol
# @return {Any|nothing} value or nothing if invalid section or key
#
function getConfig(man::ManagerTypes.ManagerData, name::Symbol)
  getfield(man.cfg, name)
end
#
# @rpc
# Updates quiet mode of manager. In quiet mode we don't show
# info and warning messages to the terminal.
# @param man Manager data type
# @param mode New quiet mode
#
function setQuiet(man::ManagerTypes.ManagerData, mode::Int)
  man.cfg.modeQuiet = mode
  nothing
end
#
# @rpc
# Does one mutation for specified organism
# @param man Manager data type
# @param organismId Unique orgaism's ID
# @param percent Mutations percent from code size
#
function mutate(man::ManagerTypes.ManagerData, organismId::UInt, percent::Float64)
  if !haskey(man.organisms, organismId) return false end
  Mutator.mutate(man.cfg, man.organisms[organismId], percent)
end
#
# @rpc
# Returns current IPS (Iterations Per Second) value
# @param man Manager data type
# @return Float64
#
function getIps(man::ManagerTypes.ManagerData)
  man.ips
end
#
# @rpc
# Returns an organism by it's unique id
# @param man Manager data type
# @param orgId Unique organism's id
# @return Creature.Organism or false if no organism with this id
#
function getOrganism(man::ManagerTypes.ManagerData, orgId::UInt)
  if !haskey(man.organisms, orgId) return false end
  _createSimpleOrganism(orgId, man.organisms[orgId])
end
#
# @rpc
# Returns amount of organisms in a world. Died organisms are
# not calculated.
# @param man Manager data type
# @return {Int}
#
function getAmount(man::ManagerTypes.ManagerData)
  ManagerTypes.orgAmount(man)
end
#
# @rpc
# Returns organism's list in specified range or all ovailable. from
# and to indexes are not unique ids. This is just indexes in man.tasks
# list.
# @param man Manager data type
# @param from Start index in man.tasks list
# @param to   End index in man.tasks list. 0 - means last item
# @return {Array{SimpleOrganism}}
#
function getOrganisms(man::ManagerTypes.ManagerData, from::Int = 1, to::Int = 0)
  local orgs::Array{RpcApi.SimpleOrganism, 1} = Array{RpcApi.SimpleOrganism, 1}()
  local len::Int = length(man.tasks)
  local i::Int

  from = from < 1 || from > len ? 1 : from
  to   = to < 1 || to > len ? len : to

  for i = from:to push!(orgs, _createSimpleOrganism(man.tasks[i].id, man.tasks[i].organism)) end

  orgs
end
#
# @rpc
# Sets specified amount of energy into the point with x,y coordinates
# @param man Manager data type
# @param x X coordinate
# @param y Y coordinate
# @param energy Amount of energy
#
function setEnergy(man::ManagerTypes.ManagerData, x::Int, y::Int, energy::UInt16)
  local pos::Helper.Point = Helper.Point(x, y)

  if World.getEnergy(man.world, pos) === Dots.INDEX_EMPTY
    World.setEnergy(man.world, pos, energy)
  end

  nothing
end
#
# @rpc
# Sets custom amount of energy points randomly by world
# @param man Manager data type
# @param amount Amount of energy points
# @param energy Amount of energy within one point
#
function setRandomEnergy(man::ManagerTypes.ManagerData, amount::Int, energy::UInt16)
  local width::Int = div(man.world.width, Helper.fastRand(5) + 1)
  local height::Int = div(man.world.height, Helper.fastRand(9))
  local xOffset::Int = div(man.world.width, Helper.fastRand(9) + 1)
  local yOffset::Int = div(man.world.height, Helper.fastRand(9) + 1)

  if man.cfg.modeQuiet <= Config.MODE_QUIET_IMPORTANT Helper.info("Creating random energy...") end
  for i::Int = 1:amount
    setEnergy(man, Helper.fastRand(width) + div(xOffset, 10), Helper.fastRand(height) + div(yOffset, 10), energy)
  end

  nothing
end
#
# @rpc
# Makes a backup of currect app instance
# @param man Manager data type
#
function doBackup(man::ManagerTypes.ManagerData)
  # TODO: JBackup should be used
  Manager.backup(man)
  nothing
end
#
# @rpc
# Returns world statistics
# @param man Manager data type
#
function getStatistics(man::ManagerTypes.ManagerData)
  #
  # Only array should be copied, not organisms
  #
  local tasks::Array{ManagerTypes.OrganismTask, 1} = copy(man.tasks)
  local moreThenOne::Bool = length(tasks) > 1
  local minOrg::Creature.Organism
  local maxOrg::Creature.Organism

  if moreThenOne
    @inbounds sort!(tasks, alg = QuickSort, lt = (t1::ManagerTypes.OrganismTask, t2::ManagerTypes.OrganismTask) -> t1.organism.energy < t2.organism.energy)
    minOrg = tasks[1].organism
    maxOrg = tasks[end].organism
  end

  RpcApi.Statistics(
    length(man.organisms),
    man.ips,
    man.totalOrganisms,
    man.world.width,
    man.world.height,
    Config.format(man.cfg),
    _createSimpleOrganism(moreThenOne ? minOrg.id : UInt(0), moreThenOne ? minOrg : Creature.create(man.cfg)),
    _createSimpleOrganism(moreThenOne ? maxOrg.id : UInt(0), moreThenOne ? maxOrg : Creature.create(man.cfg))
  )
end
#
# @rpc
# Returns best organisms (with maximum energy)
# @param man Manager data type
# @param amount amount of best organisms
# @return {Array{RpcApi.SimpleOrganism, 1}}
#
function getBest(man::ManagerTypes.ManagerData, amount::Int)
  local best::Array{RpcApi.SimpleOrganism, 1} = []
  local i::Int
  local len::Int = getAmount(man)
  local org::Creature.Organism
  #
  # Only array should be copied, not organisms
  #
  local tasks::Array{ManagerTypes.OrganismTask, 1} = copy(man.tasks)

  amount = len > amount ? amount : len
  if amount > 0
    sort!(tasks, alg = QuickSort, lt = (t1, t2) -> t1.organism.energy < t2.organism.energy)
    i = 1
    while amount > 0 && i <= len
      org = tasks[len - i + 1].organism
      if org.alive
        amount -= 1
        push!(best, _createSimpleOrganism(org.id, org))
      end
      i += 1
    end
  end

  best
end
#
# @rpc
# TODO:
function setLeftWorld(man::ManagerTypes.ManagerData)
end
nothing
#
# @rpc
# TODO:
function setRightWorld(man::ManagerTypes.ManagerData)
  nothing
end
#
# @rpc
# TODO:
function setUpWorld(man::ManagerTypes.ManagerData)
  nothing
end
#
# @rpc
# TODO:
function setDownWorld(man::ManagerTypes.ManagerData)
  nothing
end
#
# @rpc
# Inits/Starts/Stops world streaming. It streams only differences to the
# clients. It adds/removes worlds change event handler and every time,
# when some world change occures, it makes a request to all streamed
# clients. First client should initialize streaming by calling this RPC
# function with state === STREAMING_INIT. In this case it should
# return Region() data with all dots in a world. After that server should
# wait next call of this RPC function with state === STREAMING_ON.
# It means that the client is ready for streaming. When all clients are
# ready we set man.cons.streamRun flag to true and all organisms
# will be run. If at least one client is waiting(state === STREAMING_INIT)
# then server should wait (skip running organisms).
# @param man Manager data type
# TODO: rewrite this description
function setStreaming(man::ManagerTypes.ManagerData)
  local socks::Array{Base.TCPSocket, 1} = man.cons.fastServer.socks
  local on::Int = 0

  for i::Int = 1:length(socks)
    if Helper.isopen(socks[i]) on += 1 end
  end
  #
  # Only one "fast" client is supported at the moment
  #
  if on > 1 return false end

  man.cons.streamInit = true
  getRegion(man)
end
#
# @rpc
# Changes organism color to specified. This is how we can mark
# organism visually to have a possibility to find it on a screen
# and do some observations. For returning original color back,
# you have to call this function again with previous color index.
# See return string value for previous color index
# @param man Manager data type
# @return {String} Message of color change operation. Is everything
# is okay, then message with original organism color will be returned
#
function markOrganism(man::ManagerTypes.ManagerData, orgId::UInt32, colorIndex::Int)
  if !haskey(man.organisms, orgId) return "Invalid organism id: 0x$(hex(orgId))" end
  local org::Creature.Organism = man.organisms[orgId]
  local prevColorIndex::UInt16 = org.color
  #
  # Organism should be alive and color index should be in valid range
  #
  if org.energy < 1 return "Organism 0x$(hex(orgId)) has died" end
  if colorIndex < Dots.INDEX_EMPTY || colorIndex > Dots.INDEX_MAX_ORG_COLOR
    return "Invalid color $colorIndex. Should be in range: $(Dots.INDEX_EMPTY):$(Dots.INDEX_MAX_ORG_COLOR)"
  end
  org.color = colorIndex

  "Organism 0x$(hex(orgId)) marked by $colorIndex color. Previous color was $prevColorIndex"
end

#
# Creates Org type instance according to specified coordinates and organism
# in these coordinates.
# @param man Manager data type
# @param region Region of a world we are working with
# @param x X coordinate
# @param y Y coordinate
#
function _addOrg(man::ManagerTypes.ManagerData, region::RpcApi.Region, x::UInt16, y::UInt16)
  local org::Creature.Organism = man.positions[Manager._getPosId(man, Helper.Point(x, y))]
  local block::RpcApi.Block
  local offs::UInt8

  block, offs = _findBlock(man, region, x, y)
  push!(block.orgs, RpcApi.Org(offs, org.id, org.color, org.age))
end
#
# Creates energy offset in found block for specified coordinates
# @param man Manager data type
# @param region Region of a world we are working with
# @param x X coordinate
# @param y Y coordinate
#
function _addEnergy(man::ManagerTypes.ManagerData, region::RpcApi.Region, x::UInt16, y::UInt16)
  local block::RpcApi.Block
  local offs::UInt8

  block, offs = _findBlock(man, region, x, y)
  push!(block.energy, offs)
end
#
# Creates Org type instance according to specified coordinates and organism
# in these coordinates.
# @param man Manager data type
# @param region Region of a world we are working with
# @param x X coordinate
# @param y Y coordinate
# @return {Tuple{RpcApi.Region, Int, Int}}
#
function _findBlock(man::ManagerTypes.ManagerData, region::RpcApi.Region, x::UInt16, y::UInt16)
  local xBlock::Int = div(x - 1, RpcApi.BLOCK_SIZE) + 1
  local yBlock::Int = div(y - 1, RpcApi.BLOCK_SIZE) + 1
  local xOffs::Int  = xBlock * RpcApi.BLOCK_SIZE - RpcApi.BLOCK_SIZE
  local yOffs::Int  = yBlock * RpcApi.BLOCK_SIZE - RpcApi.BLOCK_SIZE

  (region.blocks[(yBlock - 1) * region.xBlocks + xBlock], UInt8((y - yOffs - 1) * RpcApi.BLOCK_SIZE + x - xOffs - 1))
end
#
# Handler of EVENT_BEFORE_RESPONSE for "fast" pooling mode. Turns
# on "fast" streaming (pooling).
# @param man Manager data type
# @param sock Socket of one ready client
#
function _onFastStreaming(man::ManagerTypes.ManagerData, sock::Base.TCPSocket, data::Array{Any, 1}, ans::Connection.Answer)
  man.cons.streamInit = false
  man.dotCallback     = (pos::Helper.Point, color::UInt16)->_onDrawEnergy(man, pos, color)
  man.moveCallback    = (pos::Helper.Point, dir::Int, color::UInt16, orgId::UInt, outOfBorder::Bool)->_onDrawOrganism(man, pos, color, dir, orgId, outOfBorder)

  Event.off(man.world.obs, World.EVENT_DOT, man.dotCallback)
  Event.off(man.world.obs, World.EVENT_MOVE, man.moveCallback)
  Event.on(man.world.obs, World.EVENT_DOT, man.dotCallback)
  Event.on(man.world.obs, World.EVENT_MOVE, man.moveCallback)
end
#
# This is a handler on world dot change. It notify remote clients
# about these changes.
# @param man Manager data type
# @param pos Dot coordinates
# @param color Dot color
#
function _onDrawEnergy(man::ManagerTypes.ManagerData, pos::Helper.Point, color::UInt16)
  local socks::Array{Base.TCPSocket, 1} = man.cons.fastServer.socks
  local off::Bool = true

  for i::Int = 1:length(socks)
    if Helper.isopen(socks[i])
      off = false
      Server.request(socks[i], FastApi.API_DOT_COLOR, UInt16(pos.x), UInt16(pos.y), color < UInt16(1) ? Dots.INDEX_EMPTY : Dots.INDEX_ENERGY)
      Event.fire(man.obs, "dotrequest", man)
    end
  end
  #
  # This is how we push all active messages to the network
  # TODO: change to yieldto() do we need this?
  yield()
  Event.fire(man.obs, "yield", man)
  #
  # All "fast" clients were disconnected
  #
  if off Event.off(man.world.obs, World.EVENT_DOT, man.dotCallback) end
  if off Event.off(man.world.obs, World.EVENT_MOVE, man.moveCallback) end
end
#
# This is a handler on world dot change. It notify remote clients
# about these changes.
# @param man Manager data type
# @param pos Dot coordinates
# @param color Dot color
#
function _onDrawOrganism(man::ManagerTypes.ManagerData, pos::Helper.Point, color::UInt16, dir::Int, orgId::UInt, outOfBorder::Bool)
  local socks::Array{Base.TCPSocket, 1} = man.cons.fastServer.socks
  local dataIndex::UInt8
  local off::Bool = true
  local hasDirection::Bool = dir !== Dots.DIRECTION_NO
  local infoBits::UInt8 = UInt8(outOfBorder)
  #
  # Encodes direction of organism to the first(unused) half byte of color.
  # e.g.: 0x000F -> 0xF000
  #
  color |= (UInt16(dir) << 12)
  for i::Int = 1:length(socks)
    if Helper.isopen(socks[i])
      off = false
      Server.request(socks[i], FastApi.API_ORG_COLOR, UInt16(pos.x), UInt16(pos.y), color, orgId, infoBits)
      Event.fire(man.obs, "dotrequest", man)
    end
  end
  #
  # This is how we push all active messages to the network
  # TODO: change to yieldto() do we need this?
  yield()
  Event.fire(man.obs, "stepyield", man)
  #
  # All "fast" clients were disconnected
  #
  if off Event.off(man.world.obs, World.EVENT_DOT, man.dotCallback) end
  if off Event.off(man.world.obs, World.EVENT_MOVE, man.moveCallback) end
end
#
# Checks if in specified X and Y is an organism
# @param man Manager data type
# @param x X coordinate
# @param y Y coordinate
# @return {Bool}
#
function _isOrganism(man::ManagerTypes.ManagerData, x::Int, y::Int)
  haskey(man.positions, Manager._getPosId(man, Helper.Point(x, y)))
end
#
# Assembles RpcApi.SimpleOrganism type from wider Creature.Organism
# @return {RpcApi.SimpleOrganism}
#
function _createSimpleOrganism(id::UInt, org::Creature.Organism)
  RpcApi.SimpleOrganism(
    id,
    org.code,
    org.mutationProbabilities,
    org.mutationsOnClonePercent,
    org.mutationPeriod,
    org.mutationPercent,
    org.energy,
    org.mutationsFromStart,
    org.color,
    org.mem,
    [org.pos.x, org.pos.y],
    org.age,
    org.cloneEnergyPercent,
    org.codeSize
  )
end
#
# Creates client and returns it's ClientConnection type. It
# uses port number provided by "XXXport" command line
# argument or default one from Config module.
# @param man Manager data type
# @param side Server side ("LEFT", "RIGHT",...)
# @return ClientConnection object
#
function _createClient(man::ManagerTypes.ManagerData, side::String)
  local con::Client.ClientConnection
  local serverPort::Int = _getPort(man, getfield(Manager, Symbol("ARG_", side, "_SERVER_PORT")), Symbol("con", ucfirst(lowercase(side)), "ServerPort"))
  local serverIp::IPv4  = _getIp(man, getfield(Manager, Symbol("ARG_", side, "_SERVER_IP")), Symbol("con", ucfirst(lowercase(side)), "ServerIp"))
  local thisPort::Int   = _getPort(man, Manager.ARG_SERVER_PORT, :conServerPort)
  local thisIp::IPv4    = _getIp(man, Manager.ARG_SERVER_IP, :conServerIp)

  if serverPort === thisPort && serverIp === thisIp
    Helper.error(string("Error creating Client for ", side, " server. Remote port and IP are similar to current. Port: ", thisPort, ", IP: ", thisIp))
    # TODO: is zero based address + port correct?
    con = Client.create(IPv4(0), 0)
  elseif serverPort === 0
    # TODO: is zero based address + port correct?
    con = Client.create(IPv4(0), 0)
  else
    con = Client.create(serverIp, serverPort)
  end
  if Client.isOk(con)
    Event.on(con.observer, Connection.EVENT_AFTER_REQUEST, (ans::Connection.Answer)->_onAfterResponse(man, side, ans))
    Event.on(con.observer, Connection.EVENT_BEFORE_RESPONSE, (data::Connection.Command, ans::Connection.Answer)->_onBeforeResponse(man, side, data, ans))
  end
  con
end
#
# Creates all available server\instance connections data objects
# and returns them. They are: self server, left, right, up, and
# down clients.
# @param man Manager data type
# @return Manager.Connections object
#
function _createConnections(man::ManagerTypes.ManagerData)
  ManagerTypes.Connections(
    false,
    _createServer(man),
    _createServer(man, true),
    _createClient(man, _SIDE_LEFT),
    _createClient(man, _SIDE_RIGHT),
    _createClient(man, _SIDE_UP),
    _createClient(man, _SIDE_DOWN),
    Dict{UInt, Creature.Organism}()
  )
end
#
# Handler of answer from remote server(instance) as a result on
# our request to it. ans.id contains success code, ans.data contains
# organism unique id.
# @param man Manager data type
# @param side Side name of answering server
# @param ans Answer data object
#
function _onAfterResponse(man::ManagerTypes.ManagerData, side::String, ans::Connection.Answer)
  local org::Creature.Organism = man.cons.frozen[ans.data]
  local pos::Helper.Point

  delete!(man.cons.frozen, ans.data)
  if ans.id !== RpcApi.RPC_ORG_STEP_OK
    if World.getEnergy(man.world, org.pos) > Dots.INDEX_EMPTY
      pos = World.getNearFreePos(man.world, org.pos)
      #
      # If there is no free space, we have to kill this frozen organism
      #
      if pos === false return false end
      org.pos = pos
    end
    Manager._createOrganism(man, org, org.pos, true)
  end

  true
end
#
# Handler of input request from remote server to current client.
# Should fill answer data structure.
# @param man Manager data type
# @param side Remote server side
# @param data Input comand object
# @param ans Answer data object we have to fill
#
function _onBeforeResponse(man::ManagerTypes.ManagerData, side::String, data::Connection.Command, ans::Connection.Answer)
  local org::Creature.Organism
  local pos::Helper.Point
  #
  # This request means, that one organism wants to step outside
  # of it's current Manager/instance into this one. We have
  # to send OK response, which means, that we obtained an organism
  # and added him into the organisms pool. data.args[1] contains
  # organism. We also have to decrease organism's energy, because
  # moving between instances it's network performance issue.
  #
  org      = data.args[1]
  ans.data = org.id
  pos      = org.pos
  if     data.fn === RpcApi.RPC_ORG_STEP_RIGHT && side === _SIDE_LEFT  pos.x = 1
  elseif data.fn === RpcApi.RPC_ORG_STEP_LEFT  && side === _SIDE_RIGHT pos.x = man.world.width
  elseif data.fn === RpcApi.RPC_ORG_STEP_DOWN  && side === _SIDE_UP    pos.y = 1
  elseif data.fn === RpcApi.RPC_ORG_STEP_UP    && side === _SIDE_DOWN  pos.y = man.world.height
  end
  #
  # Something on a way of organism or maximum amount of organisms
  # were reached. So he can't move there at the moment.
  #
  if !Manager._isFree(man, pos) || length(man.tasks) > man.cfg.worldMaxOrgs
    ans.id = RpcApi.RPC_ORG_STEP_FAIL
    return false
  end
  #
  # Everything is okay, let's add an organism to the pool
  #
  Manager._createOrganism(man, org, org.pos, true)
  _decreaseOrganismEnergy(man, org, div(org.energy * man.cfg.conStepEnergySpendPercent, 100))
  ans.id = RpcApi.RPC_ORG_STEP_OK

  true
end
#
# Handler for commands obtained from all connected clients. All supported
# commands are in _rpcApi array. If current command is undefinedin _rpcApi
# then, false will be returned.
# TODO: describe arguments!
function _onClientCommand(man::ManagerTypes.ManagerData, sock::Base.TCPSocket, cmd::Connection.Command, ans::Connection.Answer)
  if cmd.fn > length(_rpcApi) || cmd.fn < 1 return nothing end
  #
  # Only streaming mode needs socket as first parameter. because we
  # have to know destination client we are working with.
  #
  ans.data = _rpcApi[cmd.fn](man, cmd.args...)

  nothing
end
#
# Creates server and returns it's ServerConnection type. It
# uses port number provided by "port" command line
# argument or default one from Config module. It also may
# create "fast" server, which uses special "fast" mode for
# pooling. In this case IP address is the same, but port
# is different.
# @param man Manager data type
# @param fast Switcher for fast server creation
# @return Server.ServerConnection object
#
function _createServer(man::ManagerTypes.ManagerData, fast::Bool = false)
  local con::Server.ServerConnection = Server.create(
    _getIp(man, Manager.ARG_SERVER_IP, :conServerIp),
    _getPort(man, fast ? Manager.ARG_FAST_SERVER_PORT : Manager.ARG_SERVER_PORT, fast ? :conFastServerPort : :conServerPort),
    fast
  )
  local fn::Function = fast ?
    (sock::Base.TCPSocket, data::Array{Any, 1}, ans::Connection.Answer)->_onFastStreaming(man, sock, data, ans) :
    (sock::Base.TCPSocket, cmd::Connection.Command, ans::Connection.Answer)->_onClientCommand(man, sock, cmd, ans)

  Event.on(con.observer, Connection.EVENT_BEFORE_RESPONSE, fn)

  con
end
#
# Cuts specified region and return 2D array of it
# @param man Manager data type
# @param x Start X coordinate of region
# @param y Start Y coordinate of region
# @param x1 End x. 0 means all width
# @param y1 End y. 0 means all height
#
function _cutRegion(man::ManagerTypes.ManagerData, x::Int, y::Int, x1::Int, y1::Int)
  copy(man.world.data[y:y1, x:x1])
end
#
# Returns amount of blocks depending on amount of pixels. For example, if
# we have 17x17 pixels block (for block size of 16x16) we have to have 2x2
# blocks for storing this.
# @param x Start coordinate (may be X or Y)
# @param x1 End coordinate (may be X or Y)
# @return {Int} Blocks amount
#
function _calcBlocksAmount(x::Int, x1::Int)
  local blocks::Int = x1 - x + 1
  div(blocks, RpcApi.BLOCK_SIZE) + (blocks % RpcApi.BLOCK_SIZE > 0 ? 1 : 0)
end
#
# Checks if specified coordinates are out of world. If so, then world
# coordinates will be returned. Otherwise the same coordinates will be
# returned
# @param man Manager data type
# @param x Start X coordinate of region
# @param y Start Y coordinate of region
# @param x1 End x. 0 means all width
# @param y1 End y. 0 means all height
# @return {Tuple}
#
function _normalizeRegion(man::ManagerTypes.ManagerData, x::Int = 1, y::Int = 1, x1::Int = 0, y1::Int = 0)
  local maxWidth::Int  = size(man.world.data)[2]
  local maxHeight::Int = size(man.world.data)[1]

  if (x1 < 1 || x1 > maxWidth)  x1 = maxWidth  end
  if (y1 < 1 || y1 > maxHeight) y1 = maxHeight end
  if (x  < 1 || x  > maxWidth)  x  = 1 end
  if (y  < 1 || y  > maxHeight) y  = 1 end

  (x, y, x1, y1)
end
#
# An API for remove clients. This manager will be a server for them.
# Only these functions may be called by clients. For calling them,
# you have to use "Client" module.
#
const _rpcApi = Function[
  getRegion,
  createOrganisms,
  createOrganism,
  setConfig,
  getConfig,
  setQuiet,
  mutate,
  getIps,
  getOrganism,
  getAmount,
  getOrganisms,
  setEnergy,
  setRandomEnergy,
  doBackup,
  getStatistics,
  getBest,
  setLeftWorld,
  setRightWorld,
  setUpWorld,
  setDownWorld,
  setStreaming,
  markOrganism
]

max_stars_repo_name: tmptrash/jevo

max_issues_repo_name: tmptrash/jevo

max_forks_repo_name: tmptrash/jevo