text:
export TCS, seqtcs
export TCSGreedySearch
export TCSGreedySearchWithTimestamp

"""
    seqtcs(seq, background::Int=1, foreground::Int=2) -> newseq
expand `seq` with `background` and `foreground`'s indexes. For example, if `seq` is [i,j,k], then
`newseq` is [B,F,i,B,F,j,B,F,k,B], of which B is `background` index and F is `foreground` index.

# Example
    julia> seqtcs([7,3,5], 2, 4)'
    1×10 LinearAlgebra.Adjoint{Int64,Array{Int64,1}}:
     2  4  7  2  4  3  2  4  5  2
"""
function seqtcs(seq, background::Int=1, foreground::Int=2)
    L = length(seq)       # sequence length
    N = 3 * L + 1         # topology length
    label = zeros(Int, N)
    label[1:3:N] .= background
    label[2:3:N] .= foreground
    label[3:3:N] .= seq
    return label
end

"""
    TCS(p::Array{T,2}, seqlabel; background::Int=1, foreground::Int=2) -> target, lossvalue
# Inputs
+ `p`        : probability of softmax output
+ `seqlabel` : like [i,j,k], i/j/k is neither background state nor foreground state. If `p` has no label (e.g. pure noise or oov) then `seq` is [].
# Outputs
+ `target`    : target of softmax's output
+ `lossvalue` : negative log-likelyhood
"""
function TCS(p::Array{TYPE,2}, seqlabel; background::Int=1, foreground::Int=2) where TYPE
    seq  = seqtcs(seqlabel, background, foreground)
    Log0 = LogZero(TYPE)   # approximate -Inf of TYPE
    ZERO = TYPE(0)         # typed zero,e.g. Float32(0)
    S, T = size(p)         # assert p is a 2-D tensor
    L = length(seq)        # topology length
    r = fill!(Array{TYPE,2}(undef,S,T), ZERO)  # 𝜸 = p(s[k,t] | x[1:T]), k in softmax's indexing

    if L == 1
        r[seq[1],:] .= TYPE(1)
        return r, - sum(log.(p[seq[1],:]))
    end

    a = fill!(Array{TYPE,2}(undef,L,T), Log0)  # 𝜶 = p(s[k,t], x[1:t]), k in TCS topology's indexing
    b = fill!(Array{TYPE,2}(undef,L,T), Log0)  # 𝛃 = p(x[t+1:T] | s[k,t]), k in TCS topology's indexing
    a[1,1] = log(p[seq[1],1])  # background entrance
    a[2,1] = log(p[seq[2],1])  # foreground entrance
    b[L-1,T] = ZERO
    b[L-0,T] = ZERO

    # --- forward in log scale ---
	for t = 2:T
	    for s = 1:L
	        if s!=1
				R = mod(s,3)
	            if R==1 || s==2 || R==0
	                a[s,t] = LogSum2Exp(a[s,t-1], a[s-1,t-1])
	            elseif R==2
	                a[s,t] = LogSum3Exp(a[s,t-1], a[s-1,t-1], a[s-2,t-1])
	            end
	        else
	            a[s,t] = a[s,t-1]
	        end
	        a[s,t] += log(p[seq[s],t])
	    end
	end

    # --- backward in log scale ---
	for t = T-1:-1:1
		for s = L:-1:1
			Q = b[s,t+1] + log(p[seq[s],t+1])
			if s!=L
				R = mod(s,3)
				V = b[s+1,t+1] + log(p[seq[s+1],t+1])
				if R==1 || R==2 || s==L-1
					b[s,t] = LogSum2Exp(Q, V)
				elseif R==0
					b[s,t] = LogSum3Exp(Q, V, b[s+2,t+1] + log(p[seq[s+2],t+1]))
				end
			else
				b[s,t] = Q
			end
		end
	end

    # loglikely of TCS
    logsum = LogSum3Exp(Log0, a[1,1] + b[1,1], a[2,1] + b[2,1])

    # log weight --> normal probs
	g = exp.((a + b) .- logsum)

    # Reduce First line
    r[seq[1],:] .+= g[1,:]
    # Reduce other lines
    for n = 1:div(L-1,3)
        s = 3*n
        r[seq[s-1],:] .+= g[s-1,:]  # reduce forground states
        r[seq[s  ],:] .+= g[s,  :]  # reduce labels' states
        r[seq[s+1],:] .+= g[s+1,:]  # reduce background state
    end

    return r, -logsum
end

"""
    TCSGreedySearch(x::Array; background::Int=1, foreground::Int=2, dims=1) -> hypothesis
remove repeats and background/foreground of argmax(x, dims=dims)
"""
function TCSGreedySearch(x::Array; background::Int=1, foreground::Int=2, dims=1)
    hyp = Vector{Int}(undef, 0)
    idx = argmax(x,dims=dims)
    for t = 1:length(idx)
        previous = idx[t≠1 ? t-1 : t][1]
        current  = idx[t][1]
        if !((current==previous && t≠1) ||
             (current==background) ||
             (current==foreground))
            push!(hyp, current)
        end
    end
    return hyp
end

"""
    TCSGreedySearch(x::Array; background::Int=1, foreground::Int=2, dims=1) -> hypothesis, timestamp
"""
function TCSGreedySearchWithTimestamp(x::Array; background::Int=1, foreground::Int=2, dims=1)
    hyp = Vector{Int}(undef, 0)
    stp = Vector{Float32}(undef, 0)
    idx = argmax(x,dims=dims)
    T   = length(idx)
    for t = 1:T
        previous = idx[t≠1 ? t-1 : t][1]
        current  = idx[t][1]
        if !((current==previous && t≠1) ||
             (current==background) ||
             (current==foreground))
            push!(hyp, current)
            push!(stp, t / T)
        end
    end
    return hyp, stp
end

max_stars_repo_name: sonosole/Delta

max_issues_repo_name: sonosole/Delta

max_forks_repo_name: sonosole/Delta


text:
@inbounds for i = ix:ix-1+n
    x[i] = y[i] + z[i]*a
end

ext: jl


text:
using Distributions
function stan_string(p::Union{Type{Bernoulli},Bernoulli})
  try
    parameters = (params(p)[1])
    return string("bernoulli($parameters)")
  catch
    return string("bernoulli")
  end
end
function stan_string(p::Union{Type{Binomial},Binomial})
  try
    parameters = (params(p)[1],params(p)[2])
    return string("binomial$parameters")
  catch
    return string("binomial")
  end
end
function stan_string(p::Union{Type{BetaBinomial},BetaBinomial})
  try
    parameters = (params(p)[1],params(p)[2],params(p)[3])
    return string("beta_binomial$parameters")
  catch
    return string("beta_binomial")
  end
end
function stan_string(p::Union{Type{Hypergeometric},Hypergeometric})
  try
    parameters = (params(p)[1],params(p)[2],params(p)[3])
    return string("hypergeometric$parameters")
  catch
    return string("hypergeometric")
  end
end
function stan_string(p::Union{Type{Categorical},Categorical})
  try
    parameters = (params(p)[1])
    return string("categorical($parameters)")
  catch
    return string("categorical")
  end
end
function stan_string(p::Union{Type{NegativeBinomial},NegativeBinomial})
  try
    parameters = (params(p)[1],params(p)[2])
    return string("neg_binomial$parameters")
  catch
    return string("neg_binomial")
  end
end
function stan_string(p::Union{Type{Poisson},Poisson})
  try
    parameters = (params(p)[1])
    return string("poisson($parameters)")
  catch
    return string("poisson")
  end
end
function stan_string(p::Union{Type{Normal},Normal})
  try
    parameters = (params(p)[1],params(p)[2])
    return string("normal$parameters")
  catch
    return string("normal")
  end
end
function stan_string(p::Union{Type{TDist},TDist})
  try
    parameters = (params(p)[1])
    return string("student_t($parameters,0,1)")
  catch
    return string("student_t")
  end
end
function stan_string(p::Union{Type{Cauchy},Cauchy})
  try
    parameters = (params(p)[1],params(p)[2])
    return string("cauchy$parameters")
  catch
    return string("cauchy")
  end
end
function stan_string(p::Union{Type{Laplace},Laplace})
  try
    parameters = (params(p)[1],params(p)[2])
    return string("double_exponential$parameters")
  catch
    return string("double_exponential")
  end
end
function stan_string(p::Union{Type{Distributions.Logistic},Distributions.Logistic})
  try
    parameters = (params(p)[1],params(p)[2])
    return string("logistic$parameters")
  catch
    return string("logistic")
  end
end
function stan_string(p::Union{Type{Gumbel},Gumbel})
  try
    parameters = (params(p)[1],params(p)[2])
    return string("gumbel$parameters")
  catch
    return string("gumbel")
  end
end
function stan_string(p::Union{Type{LogNormal},LogNormal})
  try
    parameters = (params(p)[1],params(p)[2])
    return string("lognormal$parameters")
  catch
    return string("lognormal")
  end
end
function stan_string(p::Union{Type{Chisq},Chisq})
  try
    parameters = (params(p)[1])
    return string("chi_square($parameters)")
  catch
    return string("chi_square")
  end
end
function stan_string(p::Union{Type{Exponential},Exponential})
  try
    parameters = (params(p)[1])
    return string("exponential($parameters)")
  catch
    return string("exponential")
  end
end
function stan_string(p::Union{Type{Gamma},Gamma})
  try
    parameters = (params(p)[1],params(p)[2])
    return string("gamma$parameters")
  catch
    return string("gamma")
  end
end
function stan_string(p::Union{Type{InverseGamma},InverseGamma})
  try
    parameters = (params(p)[1],params(p)[2])
    return string("inv_gamma$parameters")
  catch
    return string("inv_gamma")
  end
end
function stan_string(p::Union{Type{Weibull},Weibull})
  try
    parameters = (params(p)[1],params(p)[2])
    return string("weibull$parameters")
  catch
    return string("weibull")
  end
end
function stan_string(p::Union{Type{Frechet},Frechet})
  try
    parameters = (params(p)[1],params(p)[2])
    return string("frechet$parameters")
  catch
    return string("frechet")
  end
end
function stan_string(p::Union{Type{Rayleigh},Rayleigh})
  try
    parameters = (params(p)[1])
    return string("rayleigh($parameters)")
  catch
    return string("rayleigh")
  end
end
function stan_string(p::Union{Type{Pareto},Pareto})
  try
    parameters = (params(p)[1],params(p)[2])
    return string("pareto$parameters")
  catch
    return string("pareto")
  end
end
function stan_string(p::Union{Type{GeneralizedPareto},GeneralizedPareto})
  try
    parameters = (params(p)[1],params(p)[2],params(p)[3])
    return string("pareto_type_2$parameters")
  catch
    return string("pareto_type_2")
  end
end
function stan_string(p::Union{Type{Beta},Beta})
  try
    parameters = (params(p)[1],params(p)[2])
    return string("beta$parameters")
  catch
    return string("beta")
  end
end
function stan_string(p::Union{Type{Uniform},Uniform})
  try
    parameters = (params(p)[1],params(p)[2])
    return string("uniform$parameters")
  catch
    return string("uniform")
  end
end
function stan_string(p::Union{Type{VonMises},VonMises})
  try
    parameters = (params(p)[1],params(p)[2])
    return string("von_mises$parameters")
  catch
    return string("von_mises")
  end
end
function stan_string(p::T) where {T<:Truncated}
  lower = p.lower
  upper = p.upper
  raw_string = stan_string(p.untruncated)
  return string(raw_string," T[$lower,$upper]")
end


text:
using LinearAlgebra
using Statistics
using Printf
using StatsFuns
#using SortingAlgorithms
#using Distributions
using StatsBase

################################################################################
# Gibbs algorithms
################################################################################

"""
Gibbs sampler using Kalman filter for the OU SDEMEM model.
"""
function gibbs_exact(R::Int,
                     y::Array,
                     dt::Real,
                     Σ_i_σ_ϵ::Real,
                     Σ_i_ϕ::Array,
                     γ_ϕ_0::Real,
                     γ_σ_ϵ_0::Real,
                     μ_i_ϕ::Array,
                     μ_i_σ_ϵ::Real,
                     α_star_ϕ::Real,
                     α_star_σ_ϵ::Real,
                     log_λ_i_ϕ::Vector,
                     log_λ_i_σ_ϵ::Real,
                     update_interval::Real,
                     start_update::Real,
                     α_power::Real,
                     startval_ϕ::Array,
                     startval_σ_ϵ::Real,
                     prior_parameters_η::Array,
                     prior_parameters_σ_ϵ::Array,
					 start_from_true_param::Bool=false)

    println("Starting Gibbs sampler using Kalman.")

    # dim for problem
    M,N = size(y)

    # initilization
    chain_ϕ = zeros(M, 3, R)
    chain_σ_ϵ = zeros(1, R)
    chain_η = zeros(6, R)
    ϕ_star = zeros(size(startval_ϕ))
    accept_vec = zeros(3,R)
    accept = false

    # first iteration
    chain_ϕ[:,:,1] = startval_ϕ
    chain_σ_ϵ[1] = startval_σ_ϵ
    print_interval = 1000

	if start_from_true_param
		chain_η[:,1] = η
	else
		# sample star values for η
		for j = 1:3

		    # get paramaters
		    μ_0_j, M_0_j, α_j, β_j = prior_parameters_η[j,:]

		    chain_η[j,1] = μ_0_j # set start value to mean
		    chain_η[j+3,1] = (α_j - 1)/β_j

		end

	end
    accept_vec[:,1] = [M;1;3] # we accapt the start values

    α_prob_accept_σ_ϵ = 0.
    α_prob_accept_ϕ = zeros(M)

    # compute likelihood value
    loglik_old_y = kalman(y, startval_σ_ϵ, startval_ϕ, dt)

    # main loop
    for r = 2:R

        # first gibbs stage, update \phi
        for m = 1:M
            ϕ_star[m,:] = chain_ϕ[m,:,r-1] + cholesky(exp(log_λ_i_ϕ[m])*Σ_i_ϕ[m]).U*randn(3)
        end

        loglik_star_y = kalman(y, chain_σ_ϵ[r-1], ϕ_star, dt)

        loglik_old_ϕ = calc_loglik_ϕ(chain_ϕ[:,:,r-1], chain_η[:,r-1])
        loglik_star_ϕ = calc_loglik_ϕ(ϕ_star, chain_η[:,r-1])

        for m = 1:M

            # TODO remove jacobian!
            #jacobian_old =  sum(chain_ϕ[m,:,r-1])
            #jacobian_star = sum(ϕ_star[m,:])

            #log_α = (loglik_star_y[m] + loglik_star_ϕ[m] + jacobian_star) - (loglik_old_y[m] + loglik_old_ϕ[m] - jacobian_old)
            log_α = (loglik_star_y[m] + loglik_star_ϕ[m]) - (loglik_old_y[m] + loglik_old_ϕ[m])

            α_prob_accept_ϕ[m] = min(1, exp(log_α))

            if any(isnan,ϕ_star[m,:]) == true || loglik_star_y[m] == NaN
                log_α = log(0)
                α_prob_accept_ϕ[m] = min(1, exp(log_α))
            end

            accept = log(rand()) < log_α

            if accept
                chain_ϕ[m,:,r] = ϕ_star[m,:]
                accept_vec[1,r] = accept_vec[1,r] +1
                loglik_old_y[m] = loglik_star_y[m]
            else
                chain_ϕ[m,:,r] = chain_ϕ[m,:,r-1]
            end

        end

        # secound gibbs stage, update σ_ϵ
        σ_ϵ_star = chain_σ_ϵ[r-1] + sqrt(exp(log_λ_i_σ_ϵ)*Σ_i_σ_ϵ)*randn()

        loglik_star_y = kalman(y, σ_ϵ_star, chain_ϕ[:,:,r], dt)

        prior_old_σ_ϵ = gammalogpdf(prior_parameters_σ_ϵ[1],1/prior_parameters_σ_ϵ[2],chain_σ_ϵ[r-1])
        prior_star_σ_ϵ = gammalogpdf(prior_parameters_σ_ϵ[1],1/prior_parameters_σ_ϵ[2],σ_ϵ_star)

        log_α = (sum(loglik_star_y) + prior_star_σ_ϵ) - (sum(loglik_old_y) + prior_old_σ_ϵ)
        α_prob_accept_σ_ϵ = min(1, exp(log_α))

        # correct if we have  NaNs
        if any(isnan,loglik_star_y) == true || σ_ϵ_star == NaN
            log_α = log(0)
            α_prob_accept_σ_ϵ = min(1, exp(log_α))
        end

        accept = log(rand()) < log_α

        if accept
            chain_σ_ϵ[r] = σ_ϵ_star
            accept_vec[2,r] = accept_vec[2,r] + 1
            loglik_old_y = loglik_star_y
        else
            chain_σ_ϵ[r] = chain_σ_ϵ[r-1]
        end

        # third stage, update η
        for j = 1:3

            # get paramaters
            μ_0_j, M_0_j, α_j, β_j = prior_parameters_η[j,:]

            ϕ_j = chain_ϕ[:,j,r-1]
            ϕ_j_mean = 1/M*sum(ϕ_j)

            α = α_j + div(M,2)
            β = β_j + 0.5*sum((ϕ_j .- ϕ_j_mean).^2) + (M*M_0_j)/(2*(M+M_0_j))*(ϕ_j_mean-μ_0_j)^2

            τ_j = rand_gamma(α)
            τ_j = τ_j/β

            μ = (M*τ_j*ϕ_j_mean)/(M*τ_j + M_0_j*τ_j) + (M_0_j*τ_j*μ_0_j)/(M*τ_j + M_0_j*τ_j)
            τ = M*τ_j + M_0_j*τ_j

            μ_j = μ + sqrt(1/τ)*randn()

            chain_η[j,r] = μ_j
            chain_η[j+3,r] = τ_j

            accept_vec[3,r] = accept_vec[3,r] +1

        end

        print_progress(r,print_interval,accept_vec,M)

        log_λ_i_σ_ϵ,μ_i_σ_ϵ,Σ_i_σ_ϵ = adaptive_update_σ_ϵ(r,update_interval,start_update,log_λ_i_σ_ϵ,μ_i_σ_ϵ,Σ_i_σ_ϵ,γ_σ_ϵ_0,α_power,chain_σ_ϵ,α_prob_accept_σ_ϵ,α_star_σ_ϵ)
        adaptive_update_ϕ!(r,update_interval,start_update,log_λ_i_ϕ,μ_i_ϕ,Σ_i_ϕ,γ_ϕ_0,α_power,M,chain_ϕ,α_prob_accept_ϕ,α_star_ϕ)

    end

    println("Ending Gibbs sampler using Kalman.")

    return chain_ϕ, chain_σ_ϵ, chain_η, accept_vec

end

"""
Gibbs sampler using CPMMH for the OU SDEMEM model.

The naive version of the Gibbs scheme (Gibbs scheme #1) where we do not integrate out the random numbers.
"""
function naive_gibbs_cpmmh(R::Int,
                           y::Array,
                           dt::Real,
                           Σ_i_σ_ϵ::Real,
                           Σ_i_ϕ::Array,
                           γ_ϕ_0::Real,
                           γ_σ_ϵ_0::Real,
                           μ_i_ϕ::Array,
                           μ_i_σ_ϵ::Real,
                           α_star_ϕ::Real,
                           α_star_σ_ϵ::Real,
                           log_λ_i_ϕ::Vector,
                           log_λ_i_σ_ϵ::Real,
                           update_interval::Real,
                           start_update::Real,
                           α_power::Real,
                           startval_ϕ::Array,
                           startval_σ_ϵ::Real,
                           prior_parameters_η::Array,
                           prior_parameters_σ_ϵ::Array,
                           nbr_particles::Int,
                           ρ::Real)

    println("Starting Gibbs sampler using CPMMH.")

    # dim for problem
    M,N = size(y)

    # initilization
    chain_ϕ = zeros(M, 3, R)
    chain_σ_ϵ = zeros(1, R)
    chain_η = zeros(6, R)
    ϕ_star = zeros(size(startval_ϕ))
    accept_vec = zeros(3,R)
    accept = false

    # first iteration
    chain_ϕ[:,:,1] = startval_ϕ
    chain_σ_ϵ[1] = startval_σ_ϵ

    if ρ > 0
        run_sort = true
    else
        run_sort = false
    end

    # sample star values for η
    for j = 1:3

        # get paramaters
        μ_0_j, M_0_j, α_j, β_j = prior_parameters_η[j,:]

        chain_η[j,1] = μ_0_j # set start value to mean
        chain_η[j+3,1] = (α_j - 1)/β_j # set start value to mode

    end

    accept_vec[:,1] = [M;1;3] # we accapt the start values

    α_prob_accept_σ_ϵ = 0.
    α_prob_accept_ϕ = zeros(M)

    print_interval = 100
    # TODO check if the random numbers are correct, now, we have ONE set of random numbers
    # but they are updated for each block

    # set random numbers
    u_prop_old = randn(nbr_particles,N+1,M)
    u_resample_old = randn(N,2,M)

    # estimate likelihood valie
    loglik_old_y = cpf(y, startval_σ_ϵ, startval_ϕ, dt, u_prop_old, u_resample_old,  nbr_particles, run_sort)

    # main loop
    for r = 2:R

        # first gibbs stage, update \phi
        for m = 1:M
            ϕ_star[m,:] = chain_ϕ[m,:,r-1] + cholesky(exp(log_λ_i_ϕ[m])*Σ_i_ϕ[m]).U*randn(3)
        end

        # update random numbers for block 1
        u_prop_new_block_1 = randn(nbr_particles,N+1,M)
        u_resample_new_block_1 = randn(N,2,M_subjects)

        u_prop_star_block_1 = ρ*u_prop_old + sqrt(1-ρ^2)*u_prop_new_block_1 # TODO check this!
        u_resample_star_block_1 = ρ*u_resample_old + sqrt(1-ρ^2)*u_resample_new_block_1

        # calc log-lik
        loglik_star_y = cpf(y, chain_σ_ϵ[r-1], ϕ_star, dt, u_prop_star_block_1, u_resample_star_block_1,  nbr_particles, run_sort)
        loglik_old_ϕ = calc_loglik_ϕ(chain_ϕ[:,:,r-1], chain_η[:,r-1])
        loglik_star_ϕ = calc_loglik_ϕ(ϕ_star, chain_η[:,r-1])

        for m = 1:M

            # TODO remove jacobian!
            #jacobian_old =  sum(chain_ϕ[m,:,r-1])
            #jacobian_star = sum(ϕ_star[m,:])

            #log_α = (loglik_star_y[m] + loglik_star_ϕ[m] + jacobian_star) - (loglik_old_y[m] + loglik_old_ϕ[m] - jacobian_old)
            log_α = (loglik_star_y[m] + loglik_star_ϕ[m]) - (loglik_old_y[m] + loglik_old_ϕ[m])

            α_prob_accept_ϕ[m] = min(1, exp(log_α))

            if any(isnan,ϕ_star[m,:]) == true || loglik_star_y[m] == NaN
                log_α = log(0)
                α_prob_accept_ϕ[m] = min(1, exp(log_α))
            end

            accept = log(rand()) < log_α

            if accept
                chain_ϕ[m,:,r] = ϕ_star[m,:]
                accept_vec[1,r] = accept_vec[1,r] +1
                loglik_old_y[m] = loglik_star_y[m]
                u_prop_old[:,:,m] = u_prop_star_block_1[:,:,m] # update random numbers for block 1
                u_resample_old[:,:,m] = u_resample_star_block_1[:,:,m]
            else
                chain_ϕ[m,:,r] = chain_ϕ[m,:,r-1]
            end

        end

        # secound gibbs stage, update σ_ϵ
        σ_ϵ_star = chain_σ_ϵ[r-1] + sqrt(exp(log_λ_i_σ_ϵ)*Σ_i_σ_ϵ)*randn()

        # update random numbers for block 2
        u_prop_new_block_2 = randn(nbr_particles,N+1,M)
        u_resample_new_block_2 = randn(N,2,M_subjects)

        u_prop_star_block_2 = ρ*u_prop_old + sqrt(1-ρ^2)*u_prop_new_block_2
        u_resample_star_block_2 = ρ*u_resample_old + sqrt(1-ρ^2)*u_resample_new_block_2

        # calc loglik
        loglik_star_y = cpf(y, σ_ϵ_star, chain_ϕ[:,:,r], dt, u_prop_star_block_2, u_resample_star_block_2,  nbr_particles, run_sort)

        prior_old_σ_ϵ = gammalogpdf(prior_parameters_σ_ϵ[1],1/prior_parameters_σ_ϵ[2],chain_σ_ϵ[r-1])
        prior_star_σ_ϵ = gammalogpdf(prior_parameters_σ_ϵ[1],1/prior_parameters_σ_ϵ[2],σ_ϵ_star)

        log_α = (sum(loglik_star_y) + prior_star_σ_ϵ) - (sum(loglik_old_y) + prior_old_σ_ϵ)
        α_prob_accept_σ_ϵ = min(1, exp(log_α))

        # correct if we have  NaNs
        if any(isnan,loglik_star_y) == true || σ_ϵ_star == NaN
            log_α = log(0)
            α_prob_accept_σ_ϵ = min(1, exp(log_α))
        end

        accept = log(rand()) < log_α

        if accept
            chain_σ_ϵ[r] = σ_ϵ_star
            accept_vec[2,r] = accept_vec[2,r] + 1
            loglik_old_y = loglik_star_y
            u_prop_old[:,:,:] = u_prop_star_block_2[:,:,:] # update random numbers for block 2
            u_resample_old[:,:,:] = u_resample_star_block_2[:,:,:]
        else
            chain_σ_ϵ[r] = chain_σ_ϵ[r-1]
        end

        # third stage, update η
        for j = 1:3

            # get paramaters
            μ_0_j, M_0_j, α_j, β_j = prior_parameters_η[j,:]

            ϕ_j = chain_ϕ[:,j,r-1]
            ϕ_j_mean = 1/M*sum(ϕ_j)

            α = α_j + div(M,2)
            β = β_j + 0.5*sum((ϕ_j .- ϕ_j_mean).^2) + (M*M_0_j)/(2*(M+M_0_j))*(ϕ_j_mean-μ_0_j)^2

            τ_j = rand_gamma(α) # TODO check scaling for Gamma
            τ_j = τ_j/β

            μ = (M*τ_j*ϕ_j_mean)/(M*τ_j + M_0_j*τ_j) + (M_0_j*τ_j*μ_0_j)/(M*τ_j + M_0_j*τ_j)
            τ = M*τ_j + M_0_j*τ_j

            μ_j = μ + sqrt(1/τ)*randn()

            chain_η[j,r] = μ_j
            chain_η[j+3,r] = τ_j

            accept_vec[3,r] = accept_vec[3,r] +1

        end

        print_progress(r,print_interval,accept_vec,M)

        log_λ_i_σ_ϵ,μ_i_σ_ϵ,Σ_i_σ_ϵ = adaptive_update_σ_ϵ(r,update_interval,start_update,log_λ_i_σ_ϵ,μ_i_σ_ϵ,Σ_i_σ_ϵ,γ_σ_ϵ_0,α_power,chain_σ_ϵ,α_prob_accept_σ_ϵ,α_star_σ_ϵ)
        adaptive_update_ϕ!(r,update_interval,start_update,log_λ_i_ϕ,μ_i_ϕ,Σ_i_ϕ,γ_ϕ_0,α_power,M,chain_ϕ,α_prob_accept_ϕ, α_star_ϕ)

    end

    println("Ending Gibbs sampler using CPMMH.")

    return chain_ϕ, chain_σ_ϵ, chain_η, accept_vec

end

"""
Gibbs sampler using CPMMH for the OU SDEMEM model.

Gibbs scheme #2 where we integrate out the random numbers.
"""
function gibbs_cpmmh(R::Int,
                     y::Array,
                     dt::Real,
                     Σ_i_σ_ϵ::Real,
                     Σ_i_ϕ::Array,
                     γ_ϕ_0::Real,
                     γ_σ_ϵ_0::Real,
                     μ_i_ϕ::Array,
                     μ_i_σ_ϵ::Real,
                     α_star_ϕ::Real,
                     α_star_σ_ϵ::Real,
                     log_λ_i_ϕ::Vector,
                     log_λ_i_σ_ϵ::Real,
                     update_interval::Real,
                     start_update::Real,
                     α_power::Real,
                     startval_ϕ::Array,
                     startval_σ_ϵ::Real,
                     startval_η::Array,
                     prior_parameters_η::Array,
                     prior_parameters_σ_ϵ::Array,
                     nbr_particles::Int,
                     ρ::Real)

    println("Starting Gibbs sampler using CPMMH.")

    # dim for problem
    M,N = size(y)

    # initilization
    chain_ϕ = zeros(M, 3, R)
    chain_σ_ϵ = zeros(1, R)
    chain_η = zeros(6, R)
    ϕ_star = zeros(size(startval_ϕ))
    accept_vec = zeros(3,R)
    accept = false

    # first iteration
    chain_ϕ[:,:,1] = startval_ϕ
    chain_σ_ϵ[1] = startval_σ_ϵ

    if ρ > 0
        run_sort = true
    else
        run_sort = false
    end

    chain_η[:,1] = startval_η

    accept_vec[:,1] = [M;1;3] # we accapt the start values

    α_prob_accept_σ_ϵ = 0.
    α_prob_accept_ϕ = zeros(M)

    print_interval = 100

    # set random numbers
    u_prop_old = randn(nbr_particles,N+1,M)
    u_resample_old = randn(N,2,M)

    # estimate likelihood valie
    loglik_old_y = cpf(y, startval_σ_ϵ, startval_ϕ, dt, u_prop_old, u_resample_old,  nbr_particles, run_sort)

    # main loop
    for r = 2:R

        for m = 1:M
            ϕ_star[m,:] = chain_ϕ[m,:,r-1] + cholesky(exp(log_λ_i_ϕ[m])*Σ_i_ϕ[m]).U*randn(3)
        end

        u_prop_new = randn(nbr_particles,N+1,M)
        u_resample_new = randn(N,2,M_subjects)

        u_prop_star = ρ*u_prop_old + sqrt(1-ρ^2)*u_prop_new # TODO check this!
        u_resample_star = ρ*u_resample_old + sqrt(1-ρ^2)*u_resample_new

        loglik_star_y = cpf(y, chain_σ_ϵ[r-1], ϕ_star, dt, u_prop_star, u_resample_star,  nbr_particles, run_sort)

        loglik_old_ϕ = calc_loglik_ϕ(chain_ϕ[:,:,r-1], chain_η[:,r-1])
        loglik_star_ϕ = calc_loglik_ϕ(ϕ_star, chain_η[:,r-1])

        for m = 1:M

            # TODO remove jacobian!
            # jacobian correct
            #jacobian_old =  sum(chain_ϕ[m,:,r-1])
            #jacobian_star = sum(ϕ_star[m,:])

            # no prior here since we have a normal prior on the log scale
            #log_α = (loglik_star_y[m] + loglik_star_ϕ[m] + jacobian_star) - (loglik_old_y[m] + loglik_old_ϕ[m] - jacobian_old)
            log_α = (loglik_star_y[m] + loglik_star_ϕ[m]) - (loglik_old_y[m] + loglik_old_ϕ[m])

            α_prob_accept_ϕ[m] = min(1, exp(log_α))

            if any(isnan,ϕ_star[m,:]) == true || loglik_star_y[m] == NaN
                log_α = log(0)
                α_prob_accept_ϕ[m] = min(1, exp(log_α))
            end

            accept = log(rand()) < log_α

            if accept
                chain_ϕ[m,:,r] = ϕ_star[m,:]
                accept_vec[1,r] = accept_vec[1,r] +1
                loglik_old_y[m] = loglik_star_y[m]
                u_prop_old[:,:,m] = u_prop_star[:,:,m]
                u_resample_old[:,:,m] = u_resample_star[:,:,m]
            else
                chain_ϕ[m,:,r] = chain_ϕ[m,:,r-1]
            end

        end

        # secound gibbs stage, update σ_ϵ
        σ_ϵ_star = chain_σ_ϵ[r-1] + sqrt(exp(log_λ_i_σ_ϵ)*Σ_i_σ_ϵ)*randn()

        loglik_star_y = cpf(y, σ_ϵ_star, chain_ϕ[:,:,r], dt, u_prop_old, u_resample_old,  nbr_particles, run_sort)

        prior_old_σ_ϵ = gammalogpdf(prior_parameters_σ_ϵ[1],1/prior_parameters_σ_ϵ[2],chain_σ_ϵ[r-1])
        prior_star_σ_ϵ = gammalogpdf(prior_parameters_σ_ϵ[1],1/prior_parameters_σ_ϵ[2],σ_ϵ_star)

        # TODO fix???
        log_α = (sum(loglik_star_y) + prior_star_σ_ϵ) - (sum(loglik_old_y) + prior_old_σ_ϵ)
        α_prob_accept_σ_ϵ = min(1, exp(log_α))

        # correct if we have  NaNs
        if any(isnan,loglik_star_y) == true || σ_ϵ_star == NaN
            log_α = log(0)
            α_prob_accept_σ_ϵ = min(1, exp(log_α))
        end

        accept = log(rand()) < log_α

        if accept
            chain_σ_ϵ[r] = σ_ϵ_star
            accept_vec[2,r] = accept_vec[2,r] + 1
            loglik_old_y = loglik_star_y
        else
            chain_σ_ϵ[r] = chain_σ_ϵ[r-1]
        end

        # third stage, update η
        for j = 1:3

            # get paramaters
            μ_0_j, M_0_j, α_j, β_j = prior_parameters_η[j,:]

            ϕ_j = chain_ϕ[:,j,r-1]
            ϕ_j_mean = 1/M*sum(ϕ_j)

            α = α_j + div(M,2)
            β = β_j + 0.5*sum((ϕ_j .- ϕ_j_mean).^2) + (M*M_0_j)/(2*(M+M_0_j))*(ϕ_j_mean-μ_0_j)^2

            τ_j = rand_gamma(α) # TODO check scaling for Gamma
            τ_j = τ_j/β

            μ = (M*τ_j*ϕ_j_mean)/(M*τ_j + M_0_j*τ_j) + (M_0_j*τ_j*μ_0_j)/(M*τ_j + M_0_j*τ_j)
            τ = M*τ_j + M_0_j*τ_j

            μ_j = μ + sqrt(1/τ)*randn()

            chain_η[j,r] = μ_j
            chain_η[j+3,r] = τ_j

            accept_vec[3,r] = accept_vec[3,r] +1

        end

        print_progress(r,print_interval,accept_vec,M)

        log_λ_i_σ_ϵ,μ_i_σ_ϵ,Σ_i_σ_ϵ = adaptive_update_σ_ϵ(r,update_interval,start_update,log_λ_i_σ_ϵ,μ_i_σ_ϵ,Σ_i_σ_ϵ,γ_σ_ϵ_0,α_power,chain_σ_ϵ,α_prob_accept_σ_ϵ,α_star_σ_ϵ)
        adaptive_update_ϕ!(r,update_interval,start_update,log_λ_i_ϕ,μ_i_ϕ,Σ_i_ϕ,γ_ϕ_0,α_power,M,chain_ϕ,α_prob_accept_ϕ,α_star_ϕ)

    end

    println("Ending Gibbs sampler using CPMMH.")

    return chain_ϕ, chain_σ_ϵ, chain_η, accept_vec

end

################################################################################
# help functions for Gibbs algorithms
################################################################################

# print progress to consol
function print_progress(r::Int,print_interval::Int,accept_vec::Array,M::Int)

    if mod(r-1,print_interval) == 0
        println("-------------")
        @printf "Percentage done:  %.2f %%\n" r/R*100
        @printf "Accaptance rate on %.0f to %.0f for ϕ: %.2f %%:\n" r-print_interval r sum(accept_vec[1,r-print_interval:r])/(M*print_interval)*100
        @printf "Accaptance rate on %.0f to %.0f for σ_ϵ: %.2f %%:\n" r-print_interval r sum(accept_vec[2,r-print_interval:r])/print_interval*100
    end

end

# Generlized AM adaptaton for the sigma_epsilon parameter,
# algorithm 4 in https://people.eecs.berkeley.edu/~jordan/sail/readings/andrieu-thoms.pdf)
function adaptive_update_σ_ϵ(r::Int,update_interval::Int,start_update::Int,
    log_λ_i_σ_ϵ::Real,μ_i_σ_ϵ::Real,Σ_i_σ_ϵ::Real,γ_σ_ϵ_0::Real,α_power::Real,
    chain_σ_ϵ::Array,α_prob_accept_σ_ϵ::Real,α_star_σ_ϵ::Real)

    if mod(r-1,update_interval) == 0 && r-1 >= start_update && r < 60000

        #γ = γ_σ_ϵ_0/(r-start_update)^α_power
        γ = γ_σ_ϵ_0/(r)^α_power

        log_λ_i_σ_ϵ = log_λ_i_σ_ϵ + γ*(α_prob_accept_σ_ϵ - α_star_σ_ϵ)
        μ_i_σ_ϵ = μ_i_σ_ϵ + γ*(chain_σ_ϵ[r] - μ_i_σ_ϵ)
        Σ_i_σ_ϵ = Σ_i_σ_ϵ + γ*((chain_σ_ϵ[r]-μ_i_σ_ϵ)^2 - Σ_i_σ_ϵ)

    end

    return log_λ_i_σ_ϵ,μ_i_σ_ϵ,Σ_i_σ_ϵ
end

# Generlized AM adaptaton for the random effects ϕ,
# algorithm 4 in https://people.eecs.berkeley.edu/~jordan/sail/readings/andrieu-thoms.pdf)
function adaptive_update_ϕ!(r::Int,update_interval::Int,start_update::Int,
    log_λ_i_ϕ::Vector,μ_i_ϕ::Array,Σ_i_ϕ::Array,γ_ϕ_0::Real,α_power::Real,
    M::Int,chain_ϕ::Array,α_prob_accept_ϕ::Vector,α_star_ϕ::Real)

    if mod(r-1,update_interval) == 0 && r-1 >= start_update && r < 60000

        #γ = γ_ϕ_0/(r-start_update)^α_power
        γ = γ_ϕ_0/(r)^α_power

        for m in 1:M
            Σ_i_ϕ_tmp = Σ_i_ϕ[m]
            log_λ_i_ϕ[m] =  log_λ_i_ϕ[m]  + γ*(α_prob_accept_ϕ[m] - α_star_ϕ)
            μ_i_ϕ[m,:] = μ_i_ϕ[m,:] + γ*(chain_ϕ[m,:,r] - μ_i_ϕ[m,:])
            Σ_i_ϕ_tmp = Σ_i_ϕ_tmp + γ*((chain_ϕ[m,:,r]-μ_i_ϕ[m,:])*(chain_ϕ[m,:,r]-μ_i_ϕ[m,:])' - Σ_i_ϕ_tmp)
            Σ_i_ϕ[m] = (Σ_i_ϕ_tmp + Σ_i_ϕ_tmp')/2 # add the transpose and diviade by 2 for numerical stability (to make sure that Σ_i_ϕ[m] is symmertic)
        end

    end

end

# calc loglik for ϕ
function calc_loglik_ϕ(ϕ, η)

    loglik = zeros(size(ϕ,1))

    for m = 1:size(ϕ,1)
        for i in 1:size(ϕ,2)
            loglik[m] = loglik[m] + normlogpdf(η[i], sqrt(1/η[i+3]), ϕ[m,i])
        end
    end

    return loglik

end

# Sample from Gamma(α,1), code adapted from https://www.cs.toronto.edu/~radford/csc2541.F04/gamma.html
function rand_gamma(α::Real)

    x = 0

    # TODO: We should have some check here such that alpha is a natrual number
    for i = 1:Int64(α)
        x = x + log(rand())
    end

    return -x

end

################################################################################
# Kalman/pf
################################################################################

# TODO rewrite the kalman and cpf in the same "style" as for the neuronal data model, to see if parallelization works

"""
Parallel Kalman filter
"""
function kalman(y_full::Array, σ_ϵ::Real, ϕ::Array, τ_s::Real, diagnostics::Bool=false)

    M = size(y_full,1)

    loglik_est = zeros(M)

    if diagnostics; path = Vector{Float64}[]; end #[zeros(length(100))]; end

    for m in 1:M

        if diagnostics
            loglik_est[m],x = kalman_filter(y_full[m,:], σ_ϵ, ϕ[m,:], τ_s,diagnostics)
            append!(path,[x])
        else
            loglik_est[m] = kalman_filter(y_full[m,:], σ_ϵ, ϕ[m,:], τ_s,diagnostics)
        end

    end

    diagnostics == true ? (return loglik_est, path) : (return loglik_est)

end

# Kalman filter. Estimate the likelihood for the OU process (LTI) system. The code follows the notation in "Grey-box pharmacokinetic/pharmacodynamic modelling of euglycaemic clamp study".
function kalman_filter(y::Vector, σ_ϵ::Real, ϕ::Vector, τ_s::Real, diagnostics::Bool=false)

    T = length(y)

    if diagnostics; x_hat_kalman = zeros(T); end

    #println(loglik_est[m])
    θ_1 = exp(ϕ[1])
    θ_2 = exp(ϕ[2])
    θ_3 = exp(ϕ[3])

    # set model
    B = θ_1*θ_2
    A = -θ_1
    σ = θ_3
    C = 1
    S = σ_ϵ^2

    # start values for the kalman filter
    P_start = var(y)
    x_hat_start = 0 # TODO: This is not correct we can fix this by sampling the start values form N(0, σ_ϵ)

    P_k = P_start
    x_k = x_hat_start

    loglik_est = 0.

    # main loop
    for k = 1:T

        x_k = exp(A*τ_s)*x_k + (1/A)*(exp(τ_s*A)-1)*B
        P_k = exp(A*τ_s)*P_k*exp(A*τ_s) + σ^2*(1/(2*A))*(exp(2*A*τ_s)-1)

        ϵ_k = y[k]-C*x_k
        R_k = C*P_k*C + S
        K = P_k*C*inv(R_k)
        x_k = x_k + K*ϵ_k
        P_k = P_k - K*R_k*K

        loglik_est = loglik_est -0.5*(log(det(R_k)) + ϵ_k*inv(R_k)*ϵ_k)

        if diagnostics; x_hat_kalman[k] = x_k; end

    end

    diagnostics == true ? (return loglik_est, x_hat_kalman) : (return loglik_est)

end

"""
Parallel cpf
"""
function cpf(y_full::Array, σ_ϵ::Real, ϕ::Array, dt::Real, u_prop::Array, u_resample::Array, N::Int, run_sort::Bool,diagnostics::Bool=false)

    M = size(y_full,1)

    loglik_est = zeros(M)

    if diagnostics; path = Matrix{Float64}[]; end

    for m in 1:M

        if diagnostics
            loglik_est[m],x = particlefilter(y_full[m,:], σ_ϵ, ϕ[m,:], dt, u_prop[:, :, m], u_resample[:, :, m], N, run_sort, diagnostics)
            push!(path, x)
        else
            loglik_est[m] = particlefilter(y_full[m,:], σ_ϵ, ϕ[m,:], dt, u_prop[:, :, m], u_resample[:, :, m], N, run_sort, diagnostics)
        end

    end

    diagnostics == true ? (return loglik_est, path) : (return loglik_est)

end

# Estimates the loglikelihood using the correlated pf.
function particlefilter(y::Vector, σ_ϵ::Real, ϕ::Array, dt::Real, u_prop::Array, u_resample_temp::Array, N::Real, run_sort::Bool, diagnostics::Bool=false)

    # pre-allocation

    loglik = 0.

    T = length(y)
    u_resample_calc = zeros(T)

    θ_1 = exp(ϕ[1]) # set model parameters
    θ_2 = exp(ϕ[2])
    θ_3 = exp(ϕ[3])

    # convert standard normal to standard uniforms
    #u_resample_temp = u_resample
    u_prop_temp = u_prop[:,2:end]
    #u_prop_temp_init = u_prop[:,1]

    for i in 1:T
      u_resample_calc[i] = exp(-(u_resample_temp[i,1]^2+u_resample_temp[i,2]^2)/2)
    end

    # pre-allocation
    x = zeros(N,T) # particels
    w = zeros(N,T) # weigts

    # set start values
    xinit = zeros(N) + std(y)*u_prop[:,1] #(std(y[m,:]) + σ_ϵ)*randn(N)

    #==
    if N > 20
        xinit = zeros(N) + std(y)*u_prop[:,1] #(std(y[m,:]) + σ_ϵ)*randn(N)
    elseif N <= 20 && N > 5
        xinit = zeros(N) + std(y)/10*u_prop[:,1] #(std(y[m,:]) + σ_ϵ)*randn(N)
    else
        xinit = zeros(N) + std(y)/20*u_prop[:,1] #(std(y[m,:]) + σ_ϵ)*randn(N)
    end
    ==#

    for t in 1:T

        if t == 1 # first iteration

          #x[:,1] = xinit

          # propagate particelsend
          #x[:,1] = stateprop(xinit, θ_1, θ_2, θ_3, dt)
          #for i = 1:N; w[i,1] = 1/N; end;
          #x[:,1] = stateprop(xinit, θ_1, θ_2, θ_3, dt)

          # propagate particels
          x[:,t] = stateprop(xinit, θ_1, θ_2, θ_3, dt, u_prop_temp[:,t]) #r*x_resample.*exp.(-x_resample .+ e[:,t])

        else

            # sort particles and weigths
            # alg=RadixSort
            if run_sort == true
                sorted_idx = sortperm(x[:,t-1]; alg = QuickSort) # ensure that we use QuickSort
                w[:,t-1] = w[sorted_idx,t-1]
                x[:,t-1] = x[sorted_idx,t-1]
            end

            # resample particels using systematic resampling
            ind = sysresample2(w[:,t-1],N,u_resample_calc[t])
            x_resample = x[ind,t-1]

            # propagate particels
            x[:,t] = stateprop(x_resample, θ_1, θ_2, θ_3, dt, u_prop_temp[:,t]) #r*x_resample.*exp.(-x_resample .+ e[:,t])

        end

        # calc weigths and update loglik
        loglik = loglik + calc_weigths(w,x[:,t],y[t],σ_ϵ,t)

    end

    diagnostics == true ? (return loglik, x) : (return loglik)

end

################################################################################
# help functions for kalman/pf
################################################################################

# sample exactly from the OU process
function stateprop(x_old::Vector, θ_1::Real, θ_2::Real, θ_3::Real, dt::Real, u::Array)

    x_prop = zeros(length(x_old))

    σ_cond = sqrt((θ_3^2/(2*θ_1))*(1-exp(-2*θ_1*dt)))

    for i in 1:length(x_old)
        μ_cond = θ_2 + (x_old[i] - θ_2)*exp(-θ_1*dt)
        x_prop[i] = μ_cond + σ_cond*u[i]
    end

    return x_prop

end

# calc weigths and update loglik
function calc_weigths(w::Array, x::Array, y::Real, σ_ϵ::Real, t::Int)

    N = length(x) # nbr particels
    logw = zeros(N)
    w_temp = zeros(N)

    # calc w
    for i in 1:N; logw[i] = log_normalpdf(x[i], σ_ϵ, y); end

    # find largets wegith
    constant = maximum(logw)

    # subtract largets weigth
    for i in 1:N; w_temp[i] = exp(logw[i] - constant); end

    # calc sum of weigths
    w_sum = sum(w_temp)

    # normalize weigths
    for i in 1:N; w_temp[i] = w_temp[i]/w_sum; end

    w[:,t] = w_temp

    # return loglik
    return constant + log(w_sum) - log(N)

end

# calc weigths and update loglik
function calc_weigths!(w::Array, loglik::Array, x::Array, y::Real, σ_ϵ::Real, m::Int, t::Int)

    N = length(x) # nbr particels
    logw = zeros(N)
    w_temp = zeros(N)

    # calc w
    @inbounds for i in 1:N; logw[i] = log_normalpdf(x[i], σ_ϵ, y); end

    # find largets wegith
    constant = maximum(logw)

    # subtract largets weigth
    @inbounds for i in 1:N; w_temp[i] = exp(logw[i] - constant); end

    # calc sum of weigths
    w_sum = sum(w_temp)

    # update loglik
    loglik[m] =  loglik[m] + constant + log(w_sum) - log(N)

    # normalize weigths
    @inbounds for i in 1:N; w_temp[i] = w_temp[i]/w_sum; end

    w[:,t] = w_temp

end

# Log-pdf for the normal distribution
function log_normalpdf(μ::Real, σ::Real, x::Real)

  R = σ^2
  ϵ = x-μ
  return -0.5*(log(det(R)) + ϵ*inv(R)*ϵ)
end

# Systematic resampling. Code adapted from Andrew Golightly.
function sysresample2(wts::Array,N::Int64,uni::Real)

  vec = zeros(Int64,N)
  wsum = sum(wts)
  k = 1
  u = uni/N
  wsumtarg = u
  wsumcurr=wts[k]/wsum
  delta = 1/N

  for i = 1:N
    while wsumcurr < wsumtarg
      k = k+1
      wsumcurr=wsumcurr+wts[k]/wsum
    end
    vec[i]=k
    wsumtarg=wsumtarg+delta
  end

  return vec

end

max_stars_repo_path: src/SDEMEM OU process/mcmc.jl

max_issues_repo_path: src/SDEMEM OU process/mcmc.jl

max_forks_repo_path: src/SDEMEM OU process/mcmc.jl


text:
using BenchmarkTools
using TrajectoryOptimization
using SNOPT7, Ipopt
using Suppressor
using Logging
const TO = TrajectoryOptimization

const suite = BenchmarkGroup()

suite["block"] = BenchmarkGroup()
suite["pendulum"] = BenchmarkGroup()
suite["acrobot"] = BenchmarkGroup()
suite["cartpole"] = BenchmarkGroup()
suite["car"] = BenchmarkGroup()

function run_benchmarks!(suite::BenchmarkGroup, prob::Problem, opts::Vector{<:TO.AbstractSolverOptions})
    for opt in opts
        println("Solver: $(solver_name(opt))")
        @time @suppress solve(prob, opt)
        suite[solver_name(opt)] = @benchmarkable @suppress solve($prob, $opt)
    end
    nothing
end

benchmarks = filter!(readdir(@__DIR__)) do file
    endswith(file, "_benchmarks.jl")
end
for file in benchmarks
    include(file)
end

block_benchmarks!(suite["block"])
pendulum_benchmarks!(suite["pendulum"])
cartpole_benchmarks!(suite["cartpole"])
car_benchmarks!(suite["car"])
# quadrotor_benchmarks!(suite["quadrotor"])
disable_logging(Logging.Info)

paramspath = joinpath(dirname(@__FILE__), "tune.json")

if isfile(paramspath)
    loadparams!(suite, BenchmarkTools.load(paramspath)[1], :evals);
else
    tune!(suite)
    BenchmarkTools.save(paramspath, params(suite));
end

SUITE = suite

max_issues_repo_name: MehrdadZareian/GMARA

max_forks_repo_name: MehrdadZareian/GMARA


text:
using BlackBoxOptim

# Lets start and then restart an optimization of 2D rosenbrock
# as in the README:
rosenbrock2d(x) = abs2(1.0 - x[1]) + 100.0 * abs2(x[2] - x[1]^2)

# If you want to restart optimization you need to get a handle
# to the optimizer, problem and params. Use the function
# setup_bboptimize instead of bboptimize, but with the same
# parameters. We just run it for 10 steps though:
opt = bbsetup(rosenbrock2d; SearchRange = (-5.0, 5.0), NumDimensions = 2,
    MaxSteps = 10, TraceMode = :silent);

# Now we can run it (this will do 10 steps):
res10 = bboptimize(opt)

# And lets run it again (another 10 steps):
res20 = bboptimize(opt)

# Lets run 980 steps this time to try to get more progress:
res1000 = bboptimize(opt; MaxSteps = 980)

# And finish it off with another 9000 steps in batches of 1000:
results = [bboptimize(opt; MaxSteps = 1000) for _ in 1:9]
res10k = results[end]

# Now compare fitness progress:
println("Fitness progress: ", map(best_fitness, [res10, res20, res1000, res10k]))

# And note that the elapsed time was much larger for first call since things
# were compiled and that the elapsed time for last two entries should be about the same
# (since they were with 980 and 1000 steps, respectively):
println("Elapsed time: ", map(BlackBoxOptim.elapsed_time, [res10, res20, res1000, res10k]))

# And results should approach the Rosenbrock optimum at [1.0, 1.0]:
println("Best solution progress: ", map(best_candidate, [res10, res20, res1000, res10k]))


text:
using Roots

# Given
i = 0.09
# compounded annually
# double the amount
log(2) / log(1 + i)

f(x) = 2 - (1 + i)^x
find_zero(f, 1)


text:
using OMEinsum, LightGraphs, BenchmarkTools

function uniformsize(@nospecialize(code::EinCode{ixs,iy}), size::Int) where {ixs, iy}
    Dict([c=>size for c in [OMEinsum.flatten(ixs)..., iy...]])
end
uniformsize(ne::OMEinsum.NestedEinsum, size::Int) = uniformsize(OMEinsum.flatten(ne), size)

function random_regular_eincode(n, k)
	g = LightGraphs.random_regular_graph(n, k)
	ixs = [minmax(e.src,e.dst) for e in LightGraphs.edges(g)]
	code = EinCode((ixs..., [(i,) for i in LightGraphs.vertices(g)]...), ())
end

number_of_nodes = 200
code = random_regular_eincode(number_of_nodes, 3)
optcode = @benchmark optimize_greedy($code, $(uniformsize(code, 2)); nrepeat=10, method=OMEinsum.MinSpaceOut())


text:
include("../../src/ToyModelsMD.jl")
using .ToyModelsMD
using DelimitedFiles
using Plots

function Afunc(system::AbstractSystem,mm::MixedModel,ind::Int64,dt::Float64)
  lam::Float64 = 2.0e3/dt
  pot::Float64 = potential(system,mm.potentials[1])
  ind>=1500 && system.x[2]>0.7 && pot<-145.0 ? lam : 0.0
end

function cbfunc(xstore::Array,system::AbstractSystem,args::Tuple{MixedModel, Tuple{Int64, Float64}})
  xstore[args[2][1],:] .= copy(system.x)
end

### define system ###
dt = 0.0001
steps = 1500
tf = steps*dt

# mueller-brown potential
mb = MullerBrown()

# gaussian force model
Mx = 21
My = 21
Mt = 21
shape = (Mx,My,Mt,2)
# Define gaussian basis centers and variance
mux = zeros(Float64,Mx)
muy = zeros(Float64,My)
mut = zeros(Float64,Mt)
for j=1:Mx
  mux[j] = -1.5+(j-1)*(1.5-(-1.5))/(Mx-1)
end
for j=1:My
  muy[j] = -0.5+(j-1)*(2-(-0.5))/(My-1)
end
for j in 1:Mt
  mut[j] = 1+(j*steps/(Mt-1)) # 1,1+2sigamt,1+4sigmat,...,1+2(Mt-1)sigmat which is steps+1
  mut[j] = (j*steps/(Mt-1)) # 1,1+2sigamt,1+4sigmat,...,1+2(Mt-1)sigmat which is steps+1
  mut[j] *= dt
end
sigx = zeros(Float64,Mx)
sigy = zeros(Float64,My)
sigt = zeros(Float64,Mt)
sigx[:] .= ((1.5-(-1.5))/(Mx-1)/2.0)^2 #variance
sigy[:] .= ((2-(-0.5))/(My-1)/2.0)^2 #variance
sigt[:] .= ((steps*1.0/(Mt-1))/2.0)^2 #variance
sigt[:] .*= dt^2
# read in optimal coefficients
x = readdlm("coeff0x.txt")
y = readdlm("coeff0y.txt")
cs = zeros(Float64,shape)
for i in 1:21
  cs[:,:,i,1] .= reshape(x[:,i],(21,21))'
  cs[:,:,i,2] .= reshape(y[:,i],(21,21))'
end
cs = reshape(cs, (21*21*21,2))
# define model - yay!
gm = GaussianModel2D(Mx,My,Mt,mux,muy,mut,sigx,sigy,sigt,cs)

# mixed model
model = MixedModel(2,[mb,gm])

# thermostat and system
langevin = Langevin(1.0,1.0)
system = ThermostattedSystem(2,0.0,[0.63,0.03],zeros(Float64,2),langevin)

# define value baseline function
shape = (Mx,My,Mt)
vmux = zeros(Float64,Mx)
vmuy = zeros(Float64,My)
vmut = zeros(Float64,Mt)
for j=1:Mx
  vmux[j] = -1.5+(j-1)*(1.5-(-1.5))/(Mx-1)
end
for j=1:My
  vmuy[j] = -0.5+(j-1)*(2-(-0.5))/(My-1)
end
for j in 1:Mt
  vmut[j] = 1+(j-1)*steps*1.0/(Mt-1) # 1,1+2sigamt,1+4sigmat,...,1+2(Mt-1)sigmat which is steps+1
  vmut[j] *= dt
end
vsigx = zeros(Float64,Mx)
vsigy = zeros(Float64,My)
vsigt = zeros(Float64,Mt)
vsigx[:] .= ((1.5-(-1.5))/(Mx-1)/2.0)^2 #variance
vsigy[:] .= ((2-(-0.5))/(My-1)/2.0)^2 #variance
vsigt[:] .= ((steps*1.0/(Mt-1))/2.0)^2 #variance
vsigt[:] .*= dt^2
# read in optimal coefficients
v = readdlm("coeff0V.txt")
vcs = zeros(Float64,shape)
for i in 1:21
  vcs[:,:,i] .= reshape(v[:,i],(21,21))'
end
vcs = reshape(vcs, (21*21*21))
vbl = GaussianValueBaseline2D(Mx,My,Mt,mux,muy,mut,sigx,sigy,sigt,vcs)

# define callback function
mcvb = MCVBCallback(Afunc,gm,vbl)
# set up integrator
integrator = StochasticEuler(system,model)

# set up optimization loop
nepochs = 500
ntraj = 100
lrf = dt*0.2
lrv = dt*5.0
dkls = zeros(Float64,nepochs)
avgAs = zeros(Float64,nepochs)
for i in 1:nepochs
  println("Epoch: $i")
  dkl = 0.0
  avgA = 0.0
  gradf = zeros(Float64,2*Mx*My*Mt)
  gradv = zeros(Float64,Mx*My*Mt)
  for j in 1:ntraj
    println("Traj: $j")
    system.x = [0.63,0.03]
    system.t = 0.0
    runtraj!(steps,dt,integrator,mcvb)
    gradf += mcvb.gradf
    gradv += mcvb.gradv
    dkl += mcvb.dkl
    avgA += Afunc(system,model,1500,dt)
  end
  # average quantities
  gradf /= ntraj
  gradv /= ntraj
  dkl /= ntraj
  avgA /= ntraj
  # update coefficients
  model.potentials[2].cs += (lrf .* reshape(gradf,(2,Mx*My*Mt))')
  mcvb.vbl.cs += (lrv .* gradv)
  dkls[i] = dkl
  avgAs[i] = avgA
  #println(dkl*dt)
end
writedlm("dkls.txt",dkls)
writedlm("avgAs.txt",avgAs)

# define callback function
xstore = zeros(steps,2)
cbf(system,args...) = cbfunc(xstore,system,args)
cb = Callback(1,cbf)
for i in 1:100
  println("Trajectory $i")
  xstore .= 0.0
  system.x = [0.63,0.03]
  system.t = 0.0
  runtraj!(steps,dt,integrator,cb)
  println(Afunc(system,model,1500,dt))
  writedlm("traj_$i.txt",xstore)
end


text:
# Method - one implementation of a Function
# Functions can have many methods

f(x, y) = 2x + y
f(x) = x + 1

println(methods(f))

println(f(1))
println(f(1, 2))

# Method ambiguity

g(x::Float64, y) = x + y
g(x, y::Float64) = x + y

println(g(1.0, 2))
println(g(1, 2.0))
# g(1, 1) # throws MethodError because no method is found with types
# g(1.0, 1.0) # throws MethodError because function call is ambiguous

max_stars_repo_name: cadamsmith/julia-programming

max_issues_repo_name: cadamsmith/julia-programming

max_forks_repo_name: cadamsmith/julia-programming


text:
using KitePodModels
using Documenter

DocMeta.setdocmeta!(KitePodModels, :DocTestSetup, :(using KitePodModels); recursive=true)

makedocs(;
    modules=[KitePodModels],
    authors="Uwe Fechner <fechner@aenarete.eu> and contributors",
    repo="https://github.com/aenarete/KitePodModels.jl/blob/{commit}{path}#{line}",
    sitename="KitePodModels.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://aenarete.github.io/KitePodModels.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Types" => "types.md",
        "Functions" => "functions.md",
        "Parameters" => "parameters.md",
    ],
)

deploydocs(;
    repo="github.com/aenarete/KitePodModels.jl",
    devbranch="main",
)


text:
module Clapeyron
using StaticArrays
using LinearAlgebra

#for the assoc solver
import PackedVectorsOfVectors
const PackedVofV = PackedVectorsOfVectors.PackedVectorOfVectors

#for non allocating vectors of zeros and ones
using FillArrays: FillArrays

using Roots: Roots
using NLSolvers
using DiffResults, ForwardDiff
using Scratch
include("solvers/Solvers.jl")
using .Solvers
using .Solvers: log
using Unitful
import LogExpFunctions
include("constants.jl")
include("models/basetools.jl") #type hierarchy
include("utils/ParamOptions.jl")
include("utils/vectors.jl")
include("utils/ClapeyronParam.jl")

include("utils/macros.jl")
using CSV, Tables
include("utils/database.jl")
include("utils/misc.jl")

include("models/combiningrules.jl")

include("models/eos.jl")
include("utils/visualisation.jl")
include("utils/split_model.jl")
include("utils/UserReader.jl")

include("models/eos/ideal/BasicIdeal.jl") #before macros, because its used there
include("models/eos/ideal/MonomerIdeal.jl")
include("models/eos/ideal/ReidIdeal.jl")
include("models/eos/ideal/WalkerIdeal.jl")
include("models/eos/ideal/JobackIdeal.jl")

include("models/eos/SAFT/PCSAFT/PCSAFT.jl")
include("models/eos/SAFT/PCSAFT/variants/sPCSAFT.jl")

include("models/eos/SAFT/ogSAFT/ogSAFT.jl")
include("models/eos/SAFT/CPA/CPA.jl")
include("models/eos/SAFT/SAFTVRSW/SAFTVRSW.jl")
include("models/eos/SAFT/LJSAFT/LJSAFT.jl")
include("models/eos/SAFT/softSAFT/softSAFT.jl")
include("models/eos/SAFT/SAFTVRMie/SAFTVRMie.jl")
include("models/eos/SAFT/SAFTVRMie/variants/SAFTVRQMie.jl")
include("models/eos/SAFT/SAFTgammaMie/SAFTgammaMie.jl")

include("models/eos/SAFT/CKSAFT/CKSAFT.jl")
include("models/eos/SAFT/CKSAFT/variants/sCKSAFT.jl")
include("models/eos/SAFT/BACKSAFT/BACKSAFT.jl")

include("models/eos/SAFT/equations.jl")

include("models/eos/cubic/equations.jl")
include("models/eos/cubic/vdW.jl")
include("models/eos/cubic/RK/RK.jl")
include("models/eos/cubic/PR/PR.jl")

include("models/eos/Activity/Wilson/Wilson.jl")
include("models/eos/Activity/NRTL/NRTL.jl")
include("models/eos/Activity/UNIQUAC/UNIQUAC.jl")
include("models/eos/Activity/UNIFAC/UNIFAC.jl")
include("models/eos/Activity/UNIFAC/variants/ogUNIFAC.jl")

include("models/eos/Activity/COSMOSAC/utils.jl")
include("models/eos/Activity/COSMOSAC/COSMOSAC02.jl")
include("models/eos/Activity/COSMOSAC/COSMOSAC10.jl")
include("models/eos/Activity/COSMOSAC/COSMOSACdsp.jl")

include("models/eos/cubic/alphas/alphas.jl")
include("models/eos/cubic/mixing/mixing.jl")
include("models/eos/cubic/translation/translation.jl")

include("models/eos/cubic/RK/variants/SRK.jl")
include("models/eos/cubic/RK/variants/PSRK.jl")
include("models/eos/cubic/PR/variants/PR78.jl")
include("models/eos/cubic/PR/variants/VTPR.jl")
include("models/eos/cubic/PR/variants/UMRPR.jl")

include("models/eos/Activity/equations.jl")

include("models/eos/EmpiricHelmholtz/IAPWS95.jl")
include("models/eos/EmpiricHelmholtz/PropaneRef.jl")
include("models/eos/EmpiricHelmholtz/GERG2008/GERG2008.jl")

include("models/eos/LatticeFluid/SanchezLacombe/SanchezLacombe.jl")

include("models/eos/SPUNG/SPUNG.jl")

include("models/eos/cached/CachedEoS.jl")

include("methods/methods.jl")

end # module


text:
"""
    rand_graph(n, m; bidirected=true, seed=-1, kws...)

Generate a random (Erdós-Renyi) `GNNGraph` with `n` nodes
and `m` edges.

If `bidirected=true` the reverse edge of each edge will be present.
If `bidirected=false` instead, `m` unrelated edges are generated.
In any case, the output graph will contain no self-loops or multi-edges.

Use a `seed > 0` for reproducibility.

Additional keyword arguments will be passed to the [`GNNGraph`](@ref) constructor.

# Usage

```juliarepl
julia> g = rand_graph(5, 4, bidirected=false)
GNNGraph:
    num_nodes = 5
    num_edges = 4

julia> edge_index(g)
([1, 3, 3, 4], [5, 4, 5, 2])

# In the bidirected case, edge data will be duplicated on the reverse edges if needed.
julia> g = rand_graph(5, 4, edata=rand(16, 2))
GNNGraph:
    num_nodes = 5
    num_edges = 4
    edata:
        e => (16, 4)

# Each edge has a reverse
julia> edge_index(g)
([1, 3, 3, 4], [3, 4, 1, 3])

```
"""
function rand_graph(n::Integer, m::Integer; bidirected=true, seed=-1, kws...)
    if bidirected
        @assert iseven(m) "Need even number of edges for bidirected graphs, given m=$m."
    end
    m2 = bidirected ? m÷2 : m
    return GNNGraph(Graphs.erdos_renyi(n, m2; is_directed=!bidirected, seed); kws...)
end


text:
@testset "affine" begin

    let rng = MersenneTwister(123456)

        @test params(DiagAffine(1.0, 1.0)) == [[1.0], [1.0]]
        @test params(DiagAffine(1.0, 1.0)) == params(DiagAffine([1.0], [1.0]))

        # Test that the application yields sane results.
        x1, x2 = randn(rng, 2, 5), randn(rng, 1, 5)
        @test apply(DiagAffine(0.0, 0.0), x2) == x2
        @test apply(DiagAffine(1.0, 0.0), x2) == x2 + 1
        @test apply(DiagAffine(ones(2), zeros(2)), x1) == x1 + ones(x1)
        @test apply(DiagAffine(1.0, 0.0), [1.0]) ≈ [2.0]

        # Test that the logdetJ yields some sane results.
        @test logdetJ(DiagAffine(0.0, 0.0), [5.0]) == zero(5.0)
        @test logdetJ(DiagAffine(0.0, 1.0), [5.0]) == one(5.0)
        @test logdetJ(DiagAffine(0.0, 2.0), [5.0]) ≈ 2 * one(5.0)
        @test logdetJ(DiagAffine(0, 0), randn(1, 5)) == zeros(5)
        @test logdetJ(DiagAffine(0, 1), randn(1, 7)) == ones(7)
    end
end


text:
type ModifySource <: AbstractSource
  original::AbstractSource
  cmd::Base.AbstractCmd
end

function download(package::AbstractPackage, source::ModifySource, workdir::String; kwargs...)
  const name = package_name(package)
  const build = builddir(package, workdir)
  download(package, source.original, workdir; kwargs...)

  cd(joinpath(build, "source", name)) do
    run(source.cmd)
  end
end

max_stars_repo_name: UCL/DebbyPacker

max_issues_repo_name: UCL/DebbyPacker

max_forks_repo_name: UCL/DebbyPacker


text:
using StatGeochem
using Documenter

DocMeta.setdocmeta!(StatGeochem, :DocTestSetup, :(using StatGeochem); recursive=true)

makedocs(;
    modules=[StatGeochem],
    authors="C. Brenhin Keller",
    repo="https://github.com/brenhinkeller/StatGeochem.jl/blob/{commit}{path}#{line}",
    sitename="StatGeochem.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://brenhinkeller.github.io/StatGeochem.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/brenhinkeller/StatGeochem.jl",
    devbranch = "main",
)


text:
# Autogenerated wrapper script for HAMdetector_model_binaries_jll for x86_64-linux-gnu-cxx11
export libtbb, libtbbmalloc, libtbbmalloc_proxy, model_1, model_2, model_3, model_4

JLLWrappers.@generate_wrapper_header("HAMdetector_model_binaries")
JLLWrappers.@declare_library_product(libtbb, "libtbb.so.2")
JLLWrappers.@declare_library_product(libtbbmalloc, "libtbbmalloc.so.2")
JLLWrappers.@declare_library_product(libtbbmalloc_proxy, "libtbbmalloc_proxy.so.2")
JLLWrappers.@declare_executable_product(model_1)
JLLWrappers.@declare_executable_product(model_2)
JLLWrappers.@declare_executable_product(model_3)
JLLWrappers.@declare_executable_product(model_4)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libtbb,
        "lib/libtbb.so.2",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libtbbmalloc,
        "lib/libtbbmalloc.so.2",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libtbbmalloc_proxy,
        "lib/libtbbmalloc_proxy.so.2",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        model_1,
        "bin/model_1",
    )

    JLLWrappers.@init_executable_product(
        model_2,
        "bin/model_2",
    )

    JLLWrappers.@init_executable_product(
        model_3,
        "bin/model_3",
    )

    JLLWrappers.@init_executable_product(
        model_4,
        "bin/model_4",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Test

# development often happens in lockstep with other packages,
# so check-out the master branch of those packages.
using Pkg
if haskey(ENV, "GITLAB_CI")
  Pkg.add([PackageSpec(name = x; rev = "master")
           for x in ["CUDAapi", "GPUArrays", "CUDAnative", "NNlib", "CUDAdrv"]])
end

include("util.jl")

using Random
Random.seed!(1)

using CuArrays

using GPUArrays
import GPUArrays: allowscalar, @allowscalar

testf(f, xs...; kwargs...) = GPUArrays.TestSuite.compare(f, CuArray, xs...; kwargs...)

allowscalar(false)

@testset "CuArrays" begin

include("base.jl")
include("dnn.jl")
include("blas.jl")
include("sparse.jl")
include("solver.jl")
include("fft.jl")
include("rand.jl")
include("sparse_solver.jl")

CuArrays.pool_status()
CuArrays.pool_timings()

end


text:
using TrajectoryOptimization
using LinearAlgebra
using Plots
using Logging

#####################################
## Minimum Time - First-order hold ##
#####################################
# NOTES
#   - Nov. 28, 2018 -T
#   - Found that increasing regularization by +10 after failed fp works well
#   - It's key to not perform an outerloop update when fp fails but to increase reg. as described above, this means convergence for dJ must be > 0.0
#   - the solves are very sensitive to initial conditions, namely max_dt since the initial dt = max_dt/2 and the initial controls
#####################################

# Pendulum
model,obj = TrajectoryOptimization.Dynamics.pendulum!
n,m = model.n, model.m

u_bound = 5.
Q = Array(1e-3*Diagonal(I,n))
R = Array(1e-3*Diagonal(I,m))
Qf = Array(Diagonal(I,n)*0.0)
tf = obj.tf
x0 = obj.x0
xf = obj.xf

obj_uncon = LQRObjective(Q, R, Qf, tf, x0, xf)
obj_c = TrajectoryOptimization.ConstrainedObjective(obj_uncon, u_min=-u_bound, u_max=u_bound) # constrained objective
obj_min = TrajectoryOptimization.update_objective(obj_c, tf=:min)

solver_min = TrajectoryOptimization.Solver(model,obj_min,integration=:rk4,N=31)

solver_min.opts.verbose = false
solver_min.opts.use_static = false
solver_min.opts.max_dt = 0.15
solver_min.opts.min_dt = 1e-3
solver_min.opts.constraint_tolerance = 0.001 # 0.005
solver_min.opts.R_minimum_time = 10.0#15.0 #13.5 # 12.0
solver_min.opts.bp_reg_initial = 0
solver_min.opts.constraint_decrease_ratio = .25
solver_min.opts.penalty_scaling = 2.0
solver_min.opts.outer_loop_update_type = :individual
solver_min.opts.iterations = 1000
solver_min.opts.iterations_outerloop = 25 # 20

U = ones(m,solver_min.N-1)
@time results_min,stats_min = TrajectoryOptimization.solve(solver_min,U)

println("Pendulum tf (minimum time): $(TrajectoryOptimization.total_time(solver_min,results_min))")

plot(TrajectoryOptimization.to_array(results_min.X)[1:2,:]')
plot(TrajectoryOptimization.to_array(results_min.U)[1:2,:]')
plot(stats_min["cost"])

# # Cartpole
# model,obj = TrajectoryOptimization.Dynamics.cartpole_analytical
# n, m = model.n, model.m
# Q = Array(0.0*Diagonal(I,n))
# R = Array(1e-3*Diagonal(I,m))
# Qf = Array(Diagonal(I,n)*0.0)
# u_bound = 20.
#
# obj_uncon = LQRObjective(Q, R, Qf, obj.tf, obj.x0, obj.xf)
# obj_c = TrajectoryOptimization.ConstrainedObjective(obj_uncon, u_min=-u_bound, u_max=u_bound) # constrained objective
# obj_min = TrajectoryOptimization.update_objective(obj_c, tf=:min)
#
# opts = TrajectoryOptimization.SolverOptions()
# opts.verbose = false
# opts.max_dt = 0.1
# opts.min_dt = 1e-3
# opts.constraint_tolerance = 0.001
# opts.R_minimum_time = 10.0 #1000.0
# opts.penalty_initial = 1.0
# opts.bp_reg_initial = 0.0
# opts.constraint_decrease_ratio = .25
# opts.penalty_scaling = 2.0
# opts.outer_loop_update_type = :default
# opts.iterations = 500
# opts.iterations_outerloop = 50 # 20
# opts.bp_reg_fp = 10.0
# opts.iterations_linesearch = 20
#
# solver = TrajectoryOptimization.Solver(model,obj,integration=:rk4,N=31,opts=opts)
# U = zeros(model.m,solver.N-1)
# results,stats = TrajectoryOptimization.solve(solver,U)
# solver_min = TrajectoryOptimization.Solver(model,obj_min,integration=:rk4,N=31,opts=opts)
# @time results_min,stats_min = TrajectoryOptimization.solve(solver_min,U)
#
# println("Cartpole tf (minimum time) : $(TrajectoryOptimization.total_time(solver_min,results_min))")
# plot(TrajectoryOptimization.to_array(results_min.X)[1:2,:]')
# plot(TrajectoryOptimization.to_array(results_min.U)[1:2,1:solver_min.N-1]')
# plot(stats_min["cost"])
# max_violation(results_min)
#
#


text:
module CVProcessing

export open_file, save_file, defaultimage, resize_image, scale_image, crop_image

using Images, ImageView, ImageTransformations, FileIO, ImageIO, ImageMagick

open_file(path) = FileIO.load(path)

save_file(path, img) = FileIO.save(path, img)

defaultimage() = FileIO.load(joinpath(@__DIR__, "../resource/img", "lighthouse.png"))

resize_image(img, w, h) = imresize(img, (convert(Int, h), convert(Int, w)))

function scale_image(img, w, h, scale)
    height, width = size(img)

    scale_percentage_h = round(convert(Int, h) / height + scale, digits = 2)
    scale_percentage_w = round(convert(Int, w) / width + scale, digits = 2)

    new_size = trunc.(Int, (height * scale_percentage_h, width * scale_percentage_w))
    imresize(img, new_size)
end

function crop_image(img, topLeftX, topLeftY, bottomRightX, bottomRightY)
    img[
        convert(Int, topLeftY):convert(Int, bottomRightY),
        convert(Int, topLeftX):convert(Int, bottomRightX),
    ]
end

end

max_stars_repo_name: bemcho/ComputerVisionWithJuliaBook

max_issues_repo_name: bemcho/ComputerVisionWithJuliaBook

max_forks_repo_name: bemcho/ComputerVisionWithJuliaBook


text:
using Extremes, Dates, DataFrames, Distributions, Plots, StatsBase, EmpiricalCDFs, CSV
pyplot()

d1 = Pareto(1, 1)
d2 = Normal(5, 1)

df = DataFrame(CSV.File("data.csv"))
sort!(df, [:date])
dropmissing!(df, :stakeValue)
filter!(row -> !(row.stakeValue == 0.0),  df)
# filter!(:stakeValue => x -> !any(f -> f(x), (ismissing, isnothing, isnan)), df)

plot(df[:stakeValue])

sv = df[:stakeValue]

empcdf = EmpiricalCDF()
append!(empcdf, sv)
sort!(empcdf)

mean(empcdf) / mad(empcdf)

empcdf(sv)
data(empcdf)

xs = rand(d1, 10_000)
dccdf = ccdf(d1, xs)

plot(xs)

plot(xs, dccdf, seriestype=:scatter, xaxis=:log, yaxis=:log)
# xlabel!("ccdf(pareto, xs)")
# ylabel!("xs")

plot(sv, empcdf(sv), seriestype=:scatter, xaxis=:log, yaxis=:log)


text:
"Holds the tableau of a Hamiltonian Partitioned Additive Runge-Kutta methods."
const TableauHPARK = AbstractTableauSPARK{:hpark}

"Parameters for right-hand side function of Hamiltonian Partitioned Additive Runge-Kutta methods."
const ParametersHPARK = AbstractParametersSPARK{:hpark}

@doc raw"""
Partitioned Additive Runge-Kutta integrator for Hamiltonian systems subject
to Dirac constraints *EXPERIMENTAL*.

This integrator solves the following system of equations for the internal stages,
```math
\begin{aligned}
Q_{n,i} &= q_{n} + h \sum \limits_{j=1}^{s} a_{ij} V_{n,j} + h \sum \limits_{j=1}^{r} \alpha_{ij} U_{n,j} , & i &= 1, ..., s , \\
P_{n,i} &= p_{n} + h \sum \limits_{j=1}^{s} a_{ij} F_{n,j} + h \sum \limits_{j=1}^{r} \alpha_{ij} G_{n,j} , & i &= 1, ..., s , \\
\tilde{Q}_{n,i} &= q_{n} + h \sum \limits_{j=1}^{s} \tilde{a}_{ij} V_{n,j} + h \sum \limits_{j=1}^{r} \tilde{\alpha}_{ij} U_{n,j} , & i &= 1, ..., r , \\
\tilde{P}_{n,i} &= p_{n} + h \sum \limits_{j=1}^{s} \tilde{a}_{ij} F_{n,j} + h \sum \limits_{j=1}^{r} \tilde{\alpha}_{ij} G_{n,j} , & i &= 1, ..., r , \\
\tilde{\Phi}_{n,i} &= 0 , & i &= 1, ..., r ,
\end{aligned}
```
with definitions
```math
\begin{aligned}
V_{n,i} &= \hphantom{-} \frac{\partial H}{\partial p} (Q_{n,i}, P_{n,i}) , & i &= 1, ..., s , \\
F_{n,i} &=           -  \frac{\partial H}{\partial q} (Q_{n,i}, P_{n,i}) , & i &= 1, ..., s , \\
U_{n,i} &= \hphantom{-} \frac{\partial \phi}{\partial p} (\tilde{Q}_{n,i}, \tilde{P}_{n,i})^{T} \Lambda_{n,i} , & i &= 1, ..., r , \\
G_{n,i} &=           -  \frac{\partial \phi}{\partial q} (\tilde{Q}_{n,i}, \tilde{P}_{n,i})^{T} \Lambda_{n,i} , & i &= 1, ..., r , \\
\tilde{\Phi}_{n,i} &= \phi(\tilde{Q}_{n,i}, \tilde{P}_{n,i}) , & i &= 1, ..., r ,
\end{aligned}
```
and update rule
```math
\begin{aligned}
q_{n+1} &= q_{n} + h \sum \limits_{i=1}^{s} b_{i} V_{n,i} + h \sum \limits_{i=1}^{r} \beta_{i} U_{n,i} , \\
p_{n+1} &= p_{n} + h \sum \limits_{i=1}^{s} b_{i} F_{n,i} + h \sum \limits_{i=1}^{r} \beta_{i} G_{n,i} .
\end{aligned}
```
"""
struct IntegratorHPARK{DT, TT, D, S, R, PT <: ParametersHPARK{DT,TT,D,S,R},
                                        ST <: NonlinearSolver{DT},
                                        IT <: InitialGuessPODE{TT}} <: AbstractIntegratorHSPARK{DT,TT,D,S,R}
    params::PT
    solver::ST
    iguess::IT
    caches::CacheDict{PT}

    function IntegratorHPARK(params::ParametersHPARK{DT,TT,D,S,R}, solver::ST, iguess::IT, caches) where {DT,TT,D,S,R,ST,IT}
        new{DT, TT, D, S, R, typeof(params), ST, IT}(params, solver, iguess, caches)
    end

    function IntegratorHPARK{DT,D}(equations::NamedTuple, tableau::TableauHPARK{TT}, Δt::TT) where {DT,TT,D}
        # get number of stages
        S = tableau.s
        R = tableau.r

        N = 2*D*S + 3*D*R

        # create params
        params = ParametersHPARK{DT,D}(equations, tableau, Δt)

        # create cache dict
        caches = CacheDict(params)

        # create solver
        solver = create_nonlinear_solver(DT, N, params, caches)

        # create initial guess
        iguess = InitialGuessPODE(get_config(:ig_extrapolation), equations[:v̄], equations[:f̄], Δt)

        # create integrator
        IntegratorHPARK(params, solver, iguess, caches)
    end

    function IntegratorHPARK(equation::Union{PDAE{DT}, HDAE{DT}}, tableau::TableauHPARK, Δt; kwargs...) where {DT}
        IntegratorHPARK{DT, ndims(equation)}(get_functions(equation), tableau, Δt; kwargs...)
    end
end

GeometricBase.nconstraints(::IntegratorHPARK{DT,TT,D}) where {DT,TT,D} = D

function compute_stages!(x::Vector{ST}, cache::IntegratorCacheSPARK{ST,D,S,R},
                                        params::ParametersHPARK{DT,TT,D,S,R}) where {ST,DT,TT,D,S,R}
    local tqᵢ::TT
    local tpᵢ::TT
    local tλᵢ::TT

    for i in 1:S
        for k in 1:D
            # copy x to Y, Z
            cache.Yi[i][k] = x[2*(D*(i-1)+k-1)+1]
            cache.Zi[i][k] = x[2*(D*(i-1)+k-1)+2]

            # compute Q and P
            cache.Qi[i][k] = params.q[k] + params.Δt * cache.Yi[i][k]
            cache.Pi[i][k] = params.p[k] + params.Δt * cache.Zi[i][k]
        end

        # compute f(X)
        tqᵢ = params.t + params.Δt * params.tab.q.c[i]
        tpᵢ = params.t + params.Δt * params.tab.p.c[i]
        params.equs[:v](tqᵢ, cache.Qi[i], cache.Pi[i], cache.Vi[i])
        params.equs[:f](tpᵢ, cache.Qi[i], cache.Pi[i], cache.Fi[i])
    end

    for i in 1:R
        for k in 1:D
            # copy y to Y, Z and Λ
            cache.Yp[i][k] = x[2*D*S+3*(D*(i-1)+k-1)+1]
            cache.Zp[i][k] = x[2*D*S+3*(D*(i-1)+k-1)+2]
            cache.Λp[i][k] = x[2*D*S+3*(D*(i-1)+k-1)+3]

            # compute Q and V
            cache.Qp[i][k] = params.q[k] + params.Δt * cache.Yp[i][k]
            cache.Pp[i][k] = params.p[k] + params.Δt * cache.Zp[i][k]
        end

        # compute f(X)
        tλᵢ = params.t + params.Δt * params.tab.λ.c[i]
        params.equs[:u](tλᵢ, cache.Qp[i], cache.Pp[i], cache.Λp[i], cache.Up[i])
        params.equs[:g](tλᵢ, cache.Qp[i], cache.Pp[i], cache.Λp[i], cache.Gp[i])
        params.equs[:ϕ](tλᵢ, cache.Qp[i], cache.Pp[i], cache.Φp[i])
    end
end

"Compute stages of variational partitioned additive Runge-Kutta methods."
function Integrators.function_stages!(y::Vector{ST}, b::Vector{ST}, params::ParametersHPARK{DT,TT,D,S,R},
                                      caches::CacheDict) where {ST,DT,TT,D,S,R}

    # get cache for internal stages
    cache = caches[ST]

    compute_stages!(y, cache, params)

    # compute b = - [(Y-AV-AU), (Z-AF-AG)]
    for i in 1:S
        for k in 1:D
            b[2*(D*(i-1)+k-1)+1] = - cache.Yi[i][k]
            b[2*(D*(i-1)+k-1)+2] = - cache.Zi[i][k]
            for j in 1:S
                b[2*(D*(i-1)+k-1)+1] += params.tab.q.a[i,j] * cache.Vi[j][k]
                b[2*(D*(i-1)+k-1)+2] += params.tab.p.a[i,j] * cache.Fi[j][k]
            end
            for j in 1:R
                b[2*(D*(i-1)+k-1)+1] += params.tab.q.α[i,j] * cache.Up[j][k]
                b[2*(D*(i-1)+k-1)+2] += params.tab.p.α[i,j] * cache.Gp[j][k]
            end
        end
    end

    # compute b = - [(Y-AV-AU), (Z-AF-AG), Φ]
    for i in 1:R
        for k in 1:D
            b[2*D*S+3*(D*(i-1)+k-1)+1] = - cache.Yp[i][k]
            b[2*D*S+3*(D*(i-1)+k-1)+2] = - cache.Zp[i][k]
            b[2*D*S+3*(D*(i-1)+k-1)+3] = - cache.Φp[i][k]
            for j in 1:S
                b[2*D*S+3*(D*(i-1)+k-1)+1] += params.tab.q̃.a[i,j] * cache.Vi[j][k]
                b[2*D*S+3*(D*(i-1)+k-1)+2] += params.tab.p̃.a[i,j] * cache.Fi[j][k]
            end
            for j in 1:R
                b[2*D*S+3*(D*(i-1)+k-1)+1] += params.tab.q̃.α[i,j] * cache.Up[j][k]
                b[2*D*S+3*(D*(i-1)+k-1)+2] += params.tab.p̃.α[i,j] * cache.Gp[j][k]
            end
        end
    end

    # compute b = - [Λ₁-λ]
    if params.tab.λ.c[1] == 0
        for k in 1:D
            b[2*D*S+3*(k-1)+3] = - cache.Λp[1][k] + params.λ[k]
        end
    end
end

function update_solution!(int::IntegratorHPARK{DT,TT}, sol::AtomicSolutionPDAE{DT,TT},
                          cache::IntegratorCacheSPARK{DT}=int.caches[DT]) where {DT,TT}
    # compute final update
    update_solution!(sol.q, sol.q̃, cache.Vi, int.params.tab.q.b, timestep(int))
    update_solution!(sol.p, sol.p̃, cache.Fi, int.params.tab.p.b, timestep(int))

    # compute projection
    update_solution!(sol.q, sol.q̃, cache.Up, int.params.tab.q.β, timestep(int))
    update_solution!(sol.p, sol.p̃, cache.Gp, int.params.tab.p.β, timestep(int))
    # update_multiplier!(sol.λ, cache.Λp, int.params.tab.λ.b)
end


text:
using POMDPs
using Random
using POMDPSimulators
using RLInterface
using BeliefUpdaters
using POMDPPolicies
using POMDPTesting

@testset "single crosswalk" begin
    params = CrosswalkParams(obstacles_visible = true)
    env = CrosswalkEnv(params)
    pomdp = SingleOCPOMDP(env = env,
                        collision_cost=-1.5)
    trans_prob_consistency_check(pomdp)

    obs_prob_consistency_check(pomdp)
end

@testset "crosswalk" begin
    rng = MersenneTwister(1)

    pomdp = OCPOMDP(ΔT = 0.5, p_birth = 0.3, max_peds = 1)

    env = KMarkovEnvironment(pomdp, k=4)

    up = KMarkovUpdater(4)

    policy = RandomPolicy(pomdp, rng=rng)

    hr = HistoryRecorder(rng=rng, max_steps = 100)
    s0 = initialstate(pomdp, rng)
    initial_observation = initialobs(pomdp, s0, rng)
    initial_obs_vec = fill(initial_observation, 4)
    hist = simulate(hr, pomdp, policy, up, initial_obs_vec, s0)

    @test n_steps(hist) > 1
end

max_stars_repo_name: MaximeBouton/AutomotivePOMDPs

max_issues_repo_name: MaximeBouton/AutomotivePOMDPs

max_forks_repo_name: MaximeBouton/AutomotivePOMDPs


text:
export ruleSPChanceConstraintOutG

const default_atol = 1e-4

"""
Computes normalizing constant and central moments of a truncated Gaussian distribution.
"""
function truncatedGaussianMoments(m::Float64, V::Float64, a::Float64, b::Float64)
    # Normalize parameters
    sigma = sqrt(clamp(V, tiny, huge))
    alpha = (a - m)/sigma
    beta = (b - m)/sigma

    # Compute normalizing constant
    Phi_alpha = normcdf(alpha)
    Phi_beta = normcdf(beta)
    Z = Phi_beta - Phi_alpha

    # Compute moments if truncation is valid
    if Z < tiny
        # Invalid region; return undefined mean and variance of truncated distribution
        Z    = 0.0
        m_tr = 0.0
        V_tr = 0.0
    else
        phi_alpha = normpdf(alpha)
        phi_beta = normpdf(beta)
        alpha_phi_alpha = clamp(alpha, -huge, huge)*phi_alpha
        beta_phi_beta = clamp(beta, -huge, huge)*phi_beta

        m_tr = m + (phi_alpha - phi_beta)/Z*sigma
        V_tr = V*( 1.0 + (alpha_phi_alpha - beta_phi_beta)/Z - ((phi_alpha - phi_beta)/Z)^2 )
    end

    return (Z, m_tr, V_tr)
end

#-------------
# Update rules
#-------------

function ruleSPChanceConstraintOutG(msg_out::Message{<:Gaussian, Univariate}, G::Tuple, epsilon::Float64; atol=default_atol)
    (m_bw, V_bw) = unsafeMeanCov(msg_out.dist)
    (xi_bw, W_bw) = unsafeWeightedMeanPrecision(msg_out.dist)

    # Default forward message is vague; only corrected when constraint is active
    xi_fw = 0.0
    W_fw = tiny

    # Correct marginal belief if constraint is active (1-D only)
    (min_G, max_G) = G # Extract minimum and maximum of safe region G
    (Phi_G, m_G, V_G) = truncatedGaussianMoments(m_bw, V_bw, min_G, max_G) # Compute statistics (and normalizing constant) of q in G
    if epsilon <= 1.0 - Phi_G # If constraint is active
        # Initialize statistics of uncorrected belief
        m_tilde = m_bw
        V_tilde = V_bw
        for i = 1:100 # Iterate at most this many times
            (Phi_lG, m_lG, V_lG) = truncatedGaussianMoments(m_tilde, V_tilde, -Inf, min_G) # Statistics for q in region left of G
            (Phi_rG, m_rG, V_rG) = truncatedGaussianMoments(m_tilde, V_tilde, max_G, Inf) # Statistics for q in region right of G

            # Compute moments of non-G region as a mixture of left and right truncations
            Phi_nG = Phi_lG + Phi_rG
            m_nG = Phi_lG/Phi_nG*m_lG + Phi_rG/Phi_nG*m_rG
            V_nG = Phi_lG/Phi_nG*(V_lG + m_lG^2) + Phi_rG/Phi_nG*(V_rG + m_rG^2) - m_nG^2

            # Compute moments of corrected belief as a mixture of G and non-G regions
            m_tilde = (1.0-epsilon)*m_G + epsilon*m_nG
            V_tilde = (1.0-epsilon)*(V_G + m_G^2) + epsilon*(V_nG + m_nG^2) - m_tilde^2

            # Re-compute statistics (and normalizing constant) of corrected belief
            (Phi_G, m_G, V_G) = truncatedGaussianMoments(m_tilde, V_tilde, min_G, max_G)

            if (1.0 - Phi_G) < (1.0 + atol)*epsilon
                break # Break the loop if the belief is sufficiently corrected
            end
        end

        # Convert moments of corrected belief to canonical form
        W_tilde = cholinv(V_tilde)
        xi_tilde = W_tilde*m_tilde

        # Compute canonical parameters of forward message
        xi_fw = xi_tilde - xi_bw
        W_fw  = W_tilde - W_bw
    end

    return Message(Univariate, GaussianWeightedMeanPrecision, xi=xi_fw, w=W_fw)
end

#--------------------------
# Custom inbounds collector
#--------------------------

function collectSumProductNodeInbounds(node::ChanceConstraint, entry::ScheduleEntry)
    inbounds = Any[]

    interface_to_schedule_entry = current_inference_algorithm.interface_to_schedule_entry
    for node_interface in node.interfaces
        inbound_interface = ultimatePartner(node_interface)
        if isa(inbound_interface.node, Clamp)
            # Hard-code outbound message of constant node in schedule
            push!(inbounds, assembleClamp!(inbound_interface.node, Message))
        else
            # Collect message from previous result
            push!(inbounds, interface_to_schedule_entry[inbound_interface])
        end
    end

    # Push function and value to calling signature; these need to be defined in the scope of the user
    push!(inbounds, Dict{Symbol, Any}(:G => node.G,
                                      :keyword => false))
    push!(inbounds, Dict{Symbol, Any}(:epsilon => node.epsilon,
                                      :keyword => false))
    if node.atol != nothing
        push!(inbounds, Dict{Symbol, Any}(:atol => node.atol,
                                          :keyword => true))
    end

    return inbounds
end


text:
"""
DNRM2 returns the euclidean norm of a vector via the function
name, so that

   DNRM2 := sqrt( x'*x )

-- This version written on 25-October-1982.
   Modified on 14-October-1993 to inline the call to DLASSQ.
   Sven Hammarling, Nag Ltd.
"""
function tron_dnrm2(n,x,incx)
    one = 1.0
    zero = 0.0

    if n < 1 || incx < 1
        xnorm = zero
    elseif n == 1
        xnorm = abs(x[1])
    else
        scale = zero
        ssq = one

        # The following loop is equivalent to this call to the LAPACK
        # auxiliary routine:
        # CALL DLASSQ( N, X, INCX, SCALE, SSQ )

        for ix=1:incx:1+(n-1)*incx
            if x[ix] != zero
                absxi = abs(x[ix])
                if scale < absxi
                    ssq = one + ssq*(scale/absxi)^2
                    scale = absxi
                else
                    ssq = ssq + (absxi/scale)^2
                end
            end
        end
        xnorm = scale*sqrt(ssq)
    end

    return xnorm
end

if isequal(BLAS_LIBRARY, :OpenBlas)
    dnrm2(n,x,incx) = BLAS.nrm2(n, x, incx)
    dnrm2(n, x::CuArray, incx) = CUBLAS.nrm2(n, x)
else
    dnrm2(n,x,incx) = tron_dnrm2(n, x, incx)
end

@inline function dnrm2(n::Int,x::CuDeviceArray{Float64,1},incx::Int)
    tx = threadIdx().x

    v = 0.0
    if tx <= n  # No check on ty so that each warp has v.
        @inbounds v = x[tx]*x[tx]
    end

    # shfl_down_sync() will automatically sync threads in a warp.

    offset = 16
    while offset > 0
        v += CUDA.shfl_down_sync(0xffffffff, v, offset)
        offset >>= 1
    end
    v = sqrt(v)
    v = CUDA.shfl_sync(0xffffffff, v, 1)

    return v
end


text:
export EulerSymp
export step!

"""
    EulerSymp(arguments;Niter,implicit,realdata)

Symplectic Euler solver.
The implicit Euler method is first used on one equation,
then the explicit Euler method is used on the second one.
The implicit equation is solved via Neumann iteration

Construct an object of type `TimeSolver` to be used in `Problem(model, initial, param; solver::TimeSolver)`

Arguments can be either
0. an object of type `AbstractModel`;
1. an `Array` of size `(N,2)` where `N` is the number of collocation points;
2. a `NamedTuple` containing a key `N`.

The keyword argument `Niter` (optional, defaut value = 10) determines the number of steps in the Neumann iteration solver of the implicit step.
The keyword argument `implicit` (optional, defaut value = 1) determines which equation is implicit (must be `1` or `2`).
The keyword argument `realdata` is optional, and determines whether pre-allocated vectors are real- or complex-valued.
By default, they are either determined by the model or the type of the array in case `0.` and `1.`, complex-valued in case `2.`.

"""
struct EulerSymp <: TimeSolver

    U1 :: Array
    U2   :: Array
    Niter:: Int
    implicit:: Int
    label :: String
    info :: String

    function EulerSymp( U :: Array; Niter = 10, implicit = 1, realdata=nothing )
        U1 = copy(U[:,1])
        U2 = copy(U[:,1])
        if realdata==true
            U1 = real.(U1);U2 = real.(U2)
        end
        if realdata==false
            U1 = complex.(U1);U2 = complex.(U2)
        end
        if implicit!=1 && implicit!=2
            @warn "the keyword `implicit` must be 1 or 2. solve! will not work."
        end
        info = "Symplectic Euler time solver: equation $implicit is solved first via \
        the implicit Euler step (using the Neumann expansion with $Niter iterations) \
        and then equation $(3-implicit) is solved via the explicit Euler step."
        label = "symplectic Euler"
        new( U1, U2, Niter, implicit, label, info)
    end

    function EulerSymp( model :: AbstractModel; Niter = 10, implicit = 1, realdata=nothing )
        U=model.mapto(Init(x->0*x,x->0*x))
        EulerSymp( U; Niter = Niter, realdata=realdata, implicit=implicit)
    end
    function EulerSymp( param::NamedTuple; Niter = 10, implicit = 1, realdata=nothing )
        EulerSymp( zeros(Complex{Float64}, (param.N,2)) ; Niter = Niter, realdata=realdata,implicit=implicit)
    end
end

function step!(solver :: EulerSymp,
                model :: AbstractModel,
                U  ,
                dt )

    solver.U1 .= copy(U[:,1])
    solver.U2 .= copy(U[:,2])

    if solver.implicit == 1
        for i=1:solver.Niter
            model.f1!( solver.U1, solver.U2 )
            solver.U1 .= U[:,1] + dt * solver.U1
        end
        U[:,1] .= solver.U1
        model.f2!( solver.U1, solver.U2 )
        U[:,2] .+= dt * solver.U2
    elseif solver.implicit == 2
        for i=1:solver.Niter
            model.f2!( solver.U1, solver.U2 )
            solver.U2 .= U[:,2] + dt * solver.U2
        end
        U[:,2] .= solver.U2
        model.f1!( solver.U1, solver.U2 )
        U[:,1] .+= dt * solver.U1
    else
        error("when defining `EulerSymp`, the keyword `implicit` must be either 1 or 2.")
    end
end


text:
import Base:+,*,-,\,^,getindex,setindex!,show,print

# ****************************************************************
# Square Block Diagonal Matrix
# ****************************************************************

export Block, size, block_idx, broadcastf, full,
    copy, getindex, setindex!, +, -, *, \, inv, square;

type Block <: AbstractMatrix{Real}

  Blocks::Array{Any}

  Block(size::Int) = new(Array(Any,size))
  Block(Blk::Array{Any}) = new(Blk)
  Block(Blk::Array) = new(convert(Array{Any}, Blk))

end

function Base.size(A::Block)
  if length(A.Blocks) == 0; return (0,0); end
  n = sum([size(B,1) for B in A.Blocks])
  return (n,n)
end

Base.size(A::Block, i::Integer)    = (i == 1 || i == 2) ? size(A)[1] : 1
getindex(A::Block, i::Integer)     = A.Blocks[i]
setindex!(A::Block, B, i::Integer) = begin; A.Blocks[i] = B; end

function block_idx(A::Block)

  #block_sizes = [size(B,1) for B in A.Blocks]

  k = size(A.Blocks,1)
  IColl = Array(UnitRange, k)

  cum_count = 1
  for i = 1:k
    blk_size = size(A.Blocks[i],1);
    IColl[i] = cum_count:(cum_count + blk_size - 1)
    cum_count += blk_size
  end

  return IColl

end

function blockIter(A::Block)

  k = size(A.Blocks,1)
  cum_count = 1
  for i = 1:k
    blk_size = size(A.Blocks[i],1)
    produce(cum_count:(cum_count + blk_size - 1))
    cum_count += blk_size
  end

end

function broadcastf(op::Function, A::Block)

  B = copy(A)
  for i = 1:length(A.Blocks)
    B[i] = op(A[i])
  end
  return B

end

function broadcastf(op::Function, A::Block, B::Block)

  C = copy(A)
  for i = 1:length(A.Blocks)
    C[i] = op(A[i], B[i])
  end
  return C

end

function broadcastf(op::Function, A::Block, x::Vector)

  y = similar(x)
  i = 1
  @inbounds for I = Task( () -> blockIter(A) )
    xI = view(x,I);
    y[I] = op(A.Blocks[i], xI)
    i += 1;
  end
  return y;

end

function broadcastf(op::Function, A::Block, X::Matrix)

  Y = similar(X)
  i = 1
  @inbounds for I = Task( () -> blockIter(A) )
    XI = view(X,I,:);
    Y[I,:] = op(A.Blocks[i],XI)
    i += 1;
  end
  return Y

end

function Base.sparse(A::Block)

  I₊, J₊, V₊ = Int[], Int[], Float64[]
  @inbounds for (I,Blk) = zip(block_idx(A), A.Blocks)
    Aᵢ = sparse(Blk)
    rows = rowvals(Aᵢ)
    vals = nonzeros(Aᵢ)
    m, n = size(Aᵢ)
    for i = 1:n
       for j in nzrange(Aᵢ, i)
          row = rows[j]; val = vals[j]
          push!(J₊, i + I[1] - 1); push!(I₊, row + I[1] - 1); push!(V₊, val)
       end
    end
  end
  return sparse(I₊,J₊,V₊);

end

function Base.full(A::Block)

  O = zeros(size(A))
  for (I,Blk) = zip(block_idx(A), A.Blocks)
    O[I,I] = full(Blk)
  end
  return O;

end

*(A::Block, X::Array{Float64,2}) = broadcastf(*,A,X)
Base.Ac_mul_B(A::Block, X::Array{Float64,2}) = broadcastf(Ac_mul_B, A, X)

*(A::Block, X::Vector) = broadcastf(*,A,X)
Base.Ac_mul_B(A::Block, X::Vector) = broadcastf(Ac_mul_B, A, X)

Base.copy(A::Block)        = Block(copy(A.Blocks))
Base.deepcopy(A::Block)    = Block(deepcopy(A.Blocks))
+(A::Block, B::Block)      = Block(A.Blocks + B.Blocks)
-(A::Block, B::Block)      = A + (-B)
Base.inv(A::Block)         = broadcastf(inv, A)
-(A::Block)                = broadcastf(-, A)
Base.ctranspose(A::Block)  = broadcastf(ctranspose, A)
*(A::Block, B::Block)      = broadcastf(*, A, B)
Base.Ac_mul_B(A::Block, B::Block) = broadcastf(Ac_mul_B, A, B)

ViewTypes = Union{SubArray}
VectorTypes = Union{Matrix, Vector, ViewTypes}

# Extra functions for dealing with views and stuff

function +(A::Diagonal, B::Block)
  i = 1
  B0 = Block(length(B.Blocks))
  for I = block_idx(B)
    dI = A.diag[I];
    B0[i] = B[i] + Diagonal(dI);
    i = i + 1
  end
  return B0
end

^(A::Block,n::Integer) = broadcastf(x -> ^(x,n), A);


text:
using MathProgBase

function _setconstrB!(m::MathProgBase.AbstractConicModel, b, K)
    error("Not supported")
end

function getLPconstrbounds(bs, Ks)
    sumlen = sum(map(length, bs))
    for i in 1:length(bs)
        @assert BitSet(1:length(bs[i])) == Compat.reduce(∪, map(c -> BitSet(c[2]), Ks[i]), init=BitSet())
    end
    lb = Vector{Float64}(undef, sumlen)
    ub = Vector{Float64}(undef, sumlen)
    offset = 0
    for i in 1:length(bs)
        b = bs[i]
        for (cone, idx) in Ks[i]
            #if !(cone in [:Zero, :NonPos, :NonNeg])
            #  error("This cone is not supported")
            #end
            offidx = offset .+ idx
            if cone == :Zero || cone == :NonPos
                lb[offidx] .= b[idx]
            else
                lb[offidx] .= -Inf
            end
            if cone == :Zero || cone == :NonNeg
                ub[offidx] .= b[idx]
            else
                ub[offidx] .= Inf
            end
        end
        offset += length(b)
    end
    lb, ub
end

function _setconstrB!(m::MathProgBase.AbstractLinearQuadraticModel, bs, Ks)
    lb, ub = getLPconstrbounds(bs, Ks)
    MathProgBase.setconstrLB!(m, lb)
    MathProgBase.setconstrUB!(m, ub)
end

function _addconstr!(m::MathProgBase.AbstractConicModel, idx, a, β, cone)
    error("Not supported")
end

function _addconstr!(m::MathProgBase.AbstractLinearQuadraticModel, idx, a, β, cone)
    lb = -Inf
    ub = Inf
    if cone == :Zero || cone == :NonPos
        lb = β
    end
    if cone == :Zero || cone == :NonNeg
        ub = β
    end
    MathProgBase.addconstr!(m, idx, Vector(a), lb, ub)
end

function _load!(model::MathProgBase.AbstractConicModel, c, A, bs, Ks, C)
    MathProgBase.loadproblem!(model, c, A, Compat.reduce(vcat, bs, init=Float64[]), Compat.reduce(vcat, Ks, init=[]), C)
end

function _load!(model::MathProgBase.AbstractLinearQuadraticModel, c, A, bs, Ks, C)
    lb, ub = getLPconstrbounds(bs, Ks)
    l  = Vector{Float64}(undef, size(A, 2))
    u  = Vector{Float64}(undef, size(A, 2))
    @assert BitSet(1:size(A, 2)) == Compat.reduce(∪, map(c -> BitSet(c[2]), C), init=BitSet())
    for (cone, idx) in C
        if !(cone  in [:Free, :NonPos, :NonNeg])
            error("This cone is not supported")
        end
        if cone == :Free || cone == :NonPos
            l[idx] .= -Inf
        else
            l[idx] .= 0
        end
        if cone == :Free || cone == :NonNeg
            u[idx] .= Inf
        else
            u[idx] .= 0
        end
    end
    # TODO sparse objective triggers a bug in Gurobi solver, report it
    MathProgBase.loadproblem!(model, A, l, u, Vector(c), lb, ub, :Min)
end

function _getdual(model::MathProgBase.AbstractConicModel)
    -MathProgBase.getdual(model)
end
function _getdual(model::MathProgBase.AbstractLinearQuadraticModel)
    if MathProgBase.status(model) == :Infeasible
        MathProgBase.getinfeasibilityray(model)
    else
        MathProgBase.getconstrduals(model)
    end
end
function _getsolution(model)
    MathProgBase.getsolution(model)
end
function _getunboundedray(model::MathProgBase.AbstractConicModel)
    error("Unbounded Ray retrieval is unsupported for conic model")
end
function _getunboundedray(model::MathProgBase.AbstractLinearQuadraticModel)
    MathProgBase.getunboundedray(model)
end
function _getobjval(model)
    # We never know, a solver could return something else
    # when it is :Unbounded or :Infeasible, thinking that
    # objval should only be called when the status is :Optimal
    if MathProgBase.status(model) == :Unbounded
        -Inf
    elseif MathProgBase.status(model) == :Infeasible
        Inf
    else
        MathProgBase.getobjval(model)
    end
end


text:
using Documenter
using Metatheory

using Metatheory.EGraphs
using Metatheory.Rules
using Metatheory.Library

for m ∈ [Metatheory]
    for i ∈ propertynames(m)
       xxx = getproperty(m, i)
       println(xxx)
    end
 end

makedocs(
    modules = [Metatheory],
    sitename = "Metatheory.jl",
    pages = [
        "index.md"
        "theories.md"
        "egraphs.md"
        "analysis.md"
        "extraction.md"
        "schedulers.md"
        "classic.md"
        "options.md"
    ])

deploydocs(repo = "github.com/0x0f0f0f/Metatheory.jl.git")


text:
# Test RepresentativeStateTable and DiagonalizationSector with spin flip
# symmetry
#
# Also tests entanglement entropy stuff, expectation value of operators in the
# momentum basis, and use of sitefilter to restrict the hamiltonian to a
# subsystem.
let
    L = 8
    hs = SpinHalfHilbertSpace(ChainLattice([L]))
    seed_state!(hs, N_up=div(L, 2))
    apply_hamiltonian = spin_half_hamiltonian(J1=1)
    rst = RepresentativeStateTable(hs, apply_hamiltonian, [spinflip_symmetry])
    @test rst.state_info_v[findfirst(isequal([1,1,1,0,1,0,0,0]), hs.indexer)].representative_index == rst.state_info_v[findfirst(isequal([0,1,1,1,0,1,0,0]), hs.indexer)].representative_index
    @test rst.state_info_v[findfirst(isequal([1,1,1,0,1,0,0,0]), hs.indexer)].representative_index != rst.state_info_v[findfirst(isequal([0,1,1,1,0,0,1,0]), hs.indexer)].representative_index

    full_ham = operator_matrix(hs, apply_hamiltonian)

    my_operator = (apply_SpSm, 1, 3)
    full_op = operator_matrix(hs, my_operator...)

    # FIXME: we should really test this using a Hamiltonian with more terms
    # (i.e. not just J1)
    tracer3 = Tracer(hs, 1:3)
    lattice_A_3 = ChainLattice([3], diagm(0 => [0]))
    subsystem3_H_A = [operator_matrix(SpinHalfHilbertSpace(lattice_A_3, ts.indexer_A), apply_hamiltonian) for ts in tracer3.sectors]

    @test diagsizes(Tracer(hs, 1:4)) == Dict{Int,Int}(1=>2,4=>2,6=>1)

    processed_length = 0
    for k_idx in eachmomentumindex(hs.lattice)
        for spinflip_idx in 0:1
            diagsect = DiagonalizationSector(rst, 1, k_idx, [spinflip_idx])
            processed_length += length(diagsect)
            reduced_ham = Matrix(construct_reduced_hamiltonian(diagsect))
            reduced_op = construct_reduced_operator(diagsect, my_operator...)
            reduced_H_A_3 = construct_reduced_operator(diagsect, apply_hamiltonian, siteidx -> siteidx in 1:3)
            fact = eigen(Hermitian((reduced_ham + reduced_ham') / 2))
            evals, evecs = fact.values, fact.vectors
            for i in 1:length(evals)
                eval = evals[i]
                evec = evecs[:,i]
                ψ = get_full_psi(diagsect, evec)

                @test eigenstate_badness(full_ham, eval, ψ) ≈ 0 atol=1e-8
                @test eigenstate_badness(hs, apply_hamiltonian, eval, ψ) ≈ 0 atol=1e-8
                check_eigenstate(full_ham, eval, ψ)
                @test_throws CheckEigenstateError check_eigenstate(full_ham, eval + 1, ψ)

                @test dot(ψ, full_op * ψ) ≈ dot(evec, reduced_op * evec) atol=1e-10

                let
                    # FIXME: could use a higher level abstraction for this pattern!
                    exv1 = mapreduce(+, 1:length(tracer3.sectors)) do sect
                        ρ_A = ExactDiag.construct_ρ_A_block(tracer3.sectors[sect], ψ)
                        dot(ρ_A', subsystem3_H_A[sect]) # trace product
                    end
                    exv2 = dot(evec, reduced_H_A_3 * evec)
                    @test exv1 ≈ exv2 atol=1e-10
                end

                if i == 1
                    for L_A in 0:div(L, 2)
                        ent_cut1 = entanglement_entropy(Tracer(hs, 1:L_A), ψ)
                        ent_cut2 = entanglement_entropy(Tracer(hs, 1:L-L_A), ψ)
                        ent_cut1_dm = entanglement_entropy(Tracer(hs, 1:L_A), ψ * ψ')
                        # FIXME: test against known results
                        @test ent_cut1 ≈ ent_cut2 atol=1e-8
                        @test ent_cut1 ≈ ent_cut1_dm atol=1e-12
                    end
                end
            end
            if k_idx == 1 && spinflip_idx == 0
                @test evals[1] ≈ -3.651093408937176
            end
        end
    end
    @test processed_length == length(hs.indexer)
end


text:
## Simulate quartets and aminoacid sequences as in Zou2019
## Using PAML, see notebook.md
## Note that this simulations is hard-coded for quartets only
## (cannot be extended to any size of tree)
## Claudia March 2020

using PhyloNetworks, Random, Distributions, Flux, HDF5
include("functions-zou2019.jl")

## Input
rseed = 03011058
L = 1550 ## sequence length
ratealpha = 0.1580288 ## see notebook.md
model = 3 ##for control file PAML
modeldatfile = "dayhoff.dat"
blL = 0.02 ##lower bound unif for BL
blU = 1.02 ##upper bound unif for BL
##long branch attraction case -------------------
lba = true
b = 0.1 ##[0.1,1.0]
rab = 2 ##a=[2b,40b]
rcb = 0.01 ##c=[0.01b,b]
## ----------------------------------------------
nrep = 5000
onehot = false ## convert seqs to one-hot encoding?

if length(ARGS) > 0
    rseed = parse(Int,ARGS[1])
    nrep = parse(Int,ARGS[2])
    if length(ARGS) > 2
        b = parse(Float64,ARGS[3])
        rab = parse(Float64,ARGS[4])
        rcb = parse(Float64,ARGS[5])
    end
    if length(ARGS) > 5
        onehot = convert(Bool,parse(Int,ARGS[6])) ## ARGS[6] must be 0/1
    end
end

Random.seed!(rseed);
seeds = sample(1:5555555555,nrep)
makeOdd!(seeds) ## we need odd seed for PAML

if onehot
    labels = zeros(nrep)
    matrices = zeros(L)
end

outfilelab = string("labels",rseed,".in")
outfileseq = string("sequences",rseed,".in")
if lba
    outfilelab = string("labels",rseed,"-",b,"-",rab,"-",rcb,".in")
    outfileseq = string("sequences",rseed,"-",b,"-",rab,"-",rcb,".in")
end
f = open(outfilelab,"w")

for i in 1:nrep
    println("=====================================================")
    @show i
    app = i == 1 ? false : true
    tree,ind = sampleRootedMetricQuartet(blL,blU, seeds[i], lba=lba, b=b, rab=rab, rcb=rcb)

    if onehot
        labels[i] = ind
    else
        global f
        write(f,string(ind))
        write(f,"\n")
    end
    namectl = string("rep-",i,".dat")
    createCtlFile(namectl, tree, seeds[i], L, ratealpha, model, modeldatfile)
    run(`./evolver 7 MCaa.dat`)
    if onehot
        run(`cp mc.paml rep-$i.paml`)
        mat = convert2onehot("mc.paml",L)
        global matrices
        matrices = hcat(matrices,mat)
    else
        global app
        writeSequence2File("mc.paml",L,outfileseq,append=app)
        run(`rm mc.paml`)
        run(`rm rep-$i.dat`)
    end
end
close(f)

if onehot
    matrices = matrices'
    matrices = matrices[2:end,:] ## 80*nrep by L: each replicate has a 80xL matrix, all stacked

    h5write("labels.h5","labels",labels)
    h5write("matrices.h5","matrices",matrices)
end

max_stars_repo_name: Forgotten/nn-phylogenetics

max_issues_repo_name: Forgotten/nn-phylogenetics

max_forks_repo_name: Forgotten/nn-phylogenetics