text:
export Variable, @variables, @parameters

struct Variable <: Expression
    name::Symbol
    dependents::Vector{Variable}
    known::Bool
    Variable(name, dependents = Variable[]; known = false) =
        new(name, dependents, known)
end

struct Constant <: Expression
    value::Number
end
Base.get(c::Constant) = c.value

Base.iszero(ex::Expression) = isa(ex, Constant) && iszero(ex.value)
Base.isone(ex::Expression)  = isa(ex, Constant) && isone(ex.value)

# Variables use isequal for equality since == is an Operation
Base.isequal(x::Variable, y::Variable) = (x.name, x.known) == (y.name, y.known)
Base.isequal(::Variable, ::Number) = false
Base.isequal(::Number, ::Variable) = false
Base.isequal(::Variable, ::Constant) = false
Base.isequal(::Constant, ::Variable) = false
Base.isequal(c::Constant, n::Number) = c.value == n
Base.isequal(n::Number, c::Constant) = c.value == n
Base.isequal(a::Constant, b::Constant) = a.value == b.value

function Base.convert(::Type{Expr}, x::Variable)
    x.known               || return x.name
    isempty(x.dependents) && return x.name
    return :($(x.name)($(convert.(Expr, x.dependents)...)))
end
Base.convert(::Type{Expr}, c::Constant) = c.value

Base.show(io::IO, x::Variable) = print(io, x.name)

# Build variables more easily
function _parse_vars(macroname, known, x)
    ex = Expr(:block)
    var_names = Symbol[]
    # if parsing things in the form of
    # begin
    #     x
    #     y
    #     z
    # end
    x = flatten_expr!(x)
    for _var in x
        iscall = isa(_var, Expr) && _var.head == :call
        issym    = _var isa Symbol
        @assert iscall || issym "@$macroname expects a tuple of expressions (`@$macroname x y z(t)`)"

        if iscall
            dependents = :(Variable[$(_var.args[2:end]...)])
            var_name = _var.args[1]
        else
            dependents = Variable[]
            var_name = _var
        end

        push!(var_names, var_name)
        expr = :($var_name = $Variable($(Meta.quot(var_name)), $dependents; known = $known))
        push!(ex.args, expr)
    end
    push!(ex.args, build_expr(:tuple, var_names))
    return ex
end
macro variables(xs...)
    esc(_parse_vars(:Variable, false, xs))
end
macro parameters(xs...)
    esc(_parse_vars(:Param, true, xs))
end


text:
function mimo_data(seed, m, n)
    Random.seed!(seed)
    # Channel
    H = randn((m, n))
    # Gaussian noise
    v = randn((m, 1))
    # True signal
    s = rand([-1, 1], n)
    # Received signal
    sigma = .0001
    y = H * s + sigma * v
    L = [hcat(H' * H, -H' * y); hcat(-y' * H, y' * y)]
    return s, H, y, L
end

function mimo_eval(s, H, y, L, XX)
    x_hat = sign.(XX[1:end-1, end])
    rank = length([eig for eig in eigen(XX).values if eig > 1e-7])
    @show decode_error = sum(abs.(x_hat - s))
    @show rank
    @show norm(y - H * x_hat)
    @show norm(y - H * s)
    @show tr(L * XX)
    return nothing
end

ext: jl


text:
using DataGenerators

regex_generator(:SCWildcardGen, ".")

@testset "choose(String) using regex containing wilcard" begin

	gn = SCWildcardGen()

	@testset "emits different ASCII strings that match regex" begin
	    td = choose(gn)
	    @test typeof(td) == String
	    @test ismatch(r"^.$", td)
	    @mcheck_values_vary td
	end

end

regex_generator(:SCQuantifiersGen, "a?b+c*d{4}e{5,6}f{7,}g{8,8}")

@testset "choose(String) using regex containing quantifiers" begin

	gn = SCQuantifiersGen()

	@testset "emits different ASCII strings that match regex" begin
	    td = choose(gn)
	    @test typeof(td) == String
	    @test ismatch(r"^a?b+c*d{4}e{5,6}f{7,}g{8,8}$", td)
	    @mcheck_values_include count(x->x=='a',td) [0,1]
	    @mcheck_values_include count(x->x=='b',td) [1,2,3]
	    @mcheck_values_include count(x->x=='c',td) [0,1,2]
	    @mcheck_values_include count(x->x=='d',td) [4,]
	    @mcheck_values_include count(x->x=='e',td) [5,6]
	    @mcheck_values_include count(x->x=='f',td) [7,8,9]
	    @mcheck_values_are count(x->x=='g',td) [8]
	end

end

regex_generator(:SCAlternationGen, "foo|bar|baz")

@testset "choose(String) using regex containing alternation" begin

	gn = SCAlternationGen()

	@testset "emits different ASCII strings that match regex" begin
	    td = choose(gn)
	    @test typeof(td) == String
	    @test ismatch(r"^foo|bar|baz$", td)
	    @mcheck_values_are td ["foo","bar","baz"]
	end

end

regex_generator(:SCBracketsGen, "a[uvw][x-z0-3]b")

@testset "choose(String) using regex containing bracket" begin

	gn = SCBracketsGen()

	@testset "emits different ASCII strings that match regex" begin
	    td = choose(gn)
	    @test typeof(td) == String
	    @test ismatch(r"^a[uvw][x-z0-3]b$", td)
	    @mcheck_values_are td[2] ['u','v','w']
	    @mcheck_values_are td[3] ['x','y','z','0','1','2','3']
	end

end

regex_generator(:SCParenthesesGen, "a(bc|de)")

@testset "choose(String) using regex containing parentheses" begin

	gn = SCParenthesesGen()

	@testset "emits different ASCII strings that match regex" begin
	    td = choose(gn)
	    @test typeof(td) == String
	    @test ismatch(r"^a(bc+|de+)$", td)
	    @mcheck_values_are td[2] ['b','d']
	    @mcheck_values_are td[3] ['c','e']
	end

end

regex_generator(:SCClassesGen, "\\s\\S\\d\\D\\w\\W")

@testset "choose(String) using regex containing classes" begin

	gn = SCClassesGen()

	@testset "emits different ASCII strings that match regex" begin
	    td = choose(gn)
	    @test typeof(td) == String
	    @test ismatch(r"^\s\S\d\D\w\W$", td)
	    @mcheck_values_vary td
	end

end

regex_generator(:SCEscapesGen, "\\.\\[\\]\\|\\?\\+\\*\\\\")

@testset "choose(String) using regex that escapes metacharacters" begin

	gn = SCEscapesGen()

	@testset "emits different ASCII strings that match regex" begin
	    td = choose(gn)
	    @test typeof(td) == String
	    @test ismatch(r"^\.\[\]\|\?\+\*\\$", td)
	end

end


text:
module Hello

using App, Luxor

function greet()
  Drawing(1000, 1000, joinpath(App.config.server_document_root, "hello-world.png"))
  origin()
  background("black")
  sethue("red")
  fontsize(50)
  text("hello world")
  finish()
end

end

max_stars_repo_name: essenciary/genie-demo-hello-world

max_issues_repo_name: essenciary/genie-demo-hello-world

max_forks_repo_name: essenciary/genie-demo-hello-world


text: println(sum(parse.(Int64, split(readline()))))

max_stars_repo_name: sugitanishi/competitive-programming

max_issues_repo_name: sugitanishi/competitive-programming

max_forks_repo_name: sugitanishi/competitive-programming


text:
function step(q, v, body, physics, get_contacts)
    M = body.M
    Mi = body.M_inv
    g = physics.g
    Δ = physics.Δ
    dt = physics.dt

    J, ϕ = get_contacts(q, body)
    fₑ = M*g
    Jₙ = J[:,2]
    Jₜ = J[:,1]
    qₚ = [0., Δ, 0.]

    if ϕ >= Δ
        v¹ = v + dt*Mi*fₑ
        q¹ = q + dt*v¹
    else
        V, p = form_lcp(J, v, body, physics)
        fᵧ = solve_lcp(V, p)
        fₙ = fᵧ[1]; fₜ₁ = fᵧ[2]; fₜ₂ = fᵧ[3]
        v¹ = v + dt*Mi*(fₑ + fₙ*Jₙ - fₜ₁*Jₜ + fₜ₂*Jₜ)
        q¹ = q + dt*v¹+qₚ
    end
    return q¹, v¹
end

function simulate(q₀, v₀, body, physics, get_contacts)
    q = [q₀]
    v = [v₀]
    T = physics.T

    for i=1:T-1
        q¹, v¹ = step(q₀, v₀, body, physics, get_contacts)
        q₀ = q¹
        v₀ = v¹
        push!(q, q¹)
        push!(v, v¹)
    end
    q = hcat(q...)
    v = hcat(v...)

    return q, v
end

function simulate_walking(z0, walker, steps...)
    l = walker.l
    flag = 1
    if length(steps) == 0
        flag = 0
        steps = [1]
    end
    θ₁ = z0[1]
    xh = -1.2
    yh = l*cos(θ₁)
    xh_start = xh

    t0 = 0.0
    dt = 4.0
    t_ode = t0
    z_ode = reshape(cat(z0, xh, yh, dims=2), (1, 6))
    zplus = zeros(4)
    z_ref = [z0[3] z0[4] 0]

    for i=1:steps[1]
        if walker.use_controller
            walker.t₀ = t0
            walker.Θ₂ = z0[3]
            walker.Θ₂̇  = z0[4]
        end
        tspan = (t0, t0+dt)
        prob = ODEProblem(single_stance!, z0, tspan, walker)
        sol = DifferentialEquations.solve(prob, Tsit5(), callback=collision_cb)

        if walker.use_controller
            n = length(sol.t)
            θ₂ᵣ, θ₂̇ᵣ, θ₂̈ᵣ = zeros(n), zeros(n), zeros(n)
            for j = 1:length(n)
                (u, θ₂ᵣ[j], θ₂̇ᵣ[j], θ₂̈ᵣ[j]) = controller(sol.t[j], sol.u[j][1:4], walker)
            end
        end

        zplus = footstrike(sol.t[end], sol.u[end, :], walker)

        z0 = zplus
        t0 = sol.t[end]
        z_temp = reduce(vcat, sol.u)
        t_temp = sol.t
        xh_temp = xh_start .+ l*sin(z_temp[1,1]).-l*sin.(z_temp[:,1])
        yh_temp = l*cos.(z_temp[:,1])
        t_ode = vcat(t_ode, t_temp[2:end])
        z_ode = vcat(z_ode, hcat(vcat(z_temp[2:end-1, :],zplus), xh_temp[2:end], yh_temp[2:end]))

        if walker.use_controller
            z_ref = vcat(z_ref, transpose(vcat(transpose(θ₂ᵣ[2:end]), transpose(θ₂̇ᵣ[2:end]), transpose(θ₂̈ᵣ[2:end]))))
            θ₂ᵣ, θ₂̇ᵣ, θ₂̈ᵣ = [],[],[]
        end

        xh_start = xh_temp[end]

    end
    z = zplus[1:4]
    if flag == 1
        return (z_ode, t_ode)
    else
        return z
    end
end


text:
using Base.Test
using LinearAlgebraicRepresentation

@testset "Bounding boxes building test" begin
    V = [.56 .28; .84 .57; .35  1.0; .22  .43]
    @test LinearAlgebraicRepresentation.bbox(V) == ([.22 .28], [.84 1.0])
end

@testset "Bounding boxes containment test" begin
    bboxA = ([0. 0.], [1. 1.])
    bboxB = ([.5 .5], [1.5 1.5])
    bboxC = ([1. 1.], [1.25 1.25])
    bboxD = ([0 .75], [.25 1])
    bboxE = ([0 1.25], [.25 1.5])

    @test LinearAlgebraicRepresentation.bbox_contains(bboxA, bboxD)
    @test LinearAlgebraicRepresentation.bbox_contains(bboxB, bboxC)
    @test !LinearAlgebraicRepresentation.bbox_contains(bboxA, bboxB)
    @test !LinearAlgebraicRepresentation.bbox_contains(bboxA, bboxE)
end
@testset "Face area calculation test" begin
    V = Float64[2 1; 1 2; 0 0; 1 1; 2 0; 0 2]
    EV = spzeros(Int8, 6, 6)
    EV[1, [1, 4]] = [-1, 1]; EV[2, [2, 4]] = [-1, 1]
    EV[3, [2, 6]] = [-1, 1]; EV[4, [3, 6]] = [-1, 1]
    EV[5, [3, 5]] = [-1, 1]; EV[6, [1, 5]] = [-1, 1]
    FE = spzeros(Int8, 2, 6)
    FE[1, :] = [ 1 -1  1 -1  1 -1]
    FE[2, :] = [-1  1 -1  1 -1  1]

    @test LinearAlgebraicRepresentation.face_area(V, EV, FE[1,:]) == -LinearAlgebraicRepresentation.face_area(V, EV, FE[2,:])
end


text:
using Gridap
import Gridap: ∇

# Define manufactured functions
ufun(x) = x[1] + x[2]
ufun_grad(x) = VectorValue(1.0,1.0,0.0)
∇(::typeof(ufun)) = ufun_grad
bfun(x) = 0.0

# Define forms of the problem
a(v,u) = inner(∇(v), ∇(u))

# Define norms to measure the error
l2(u) = inner(u,u)
h1(u) = a(u,u) + l2(u)

function main(n)

  # Construct the discrete model
  println("model")
  @time model = CartesianDiscreteModel(partition=(n,n,n))

  # Construct the FEspace
  order = 1
  diritag = "boundary"
  println("fespace")
  @time fespace = ConformingFESpace(Float64,model,order,diritag)

  # Define test and trial spaces
  println("V")
  @time V = TestFESpace(fespace)

  println("U")
  @time U = TrialFESpace(fespace,ufun)

  # Define integration mesh and quadrature
  println("trian")
  @time trian = Triangulation(model)

  println("quad")
  @time quad = CellQuadrature(trian,degree=2)

  # Define the source term
  println("bfield")
  @time bfield = CellField(trian,bfun)

  # Define Assembler
  println("assem")
  @time assem = SparseMatrixAssembler(V,U)

  # Define the FEOperator
  b(v) = inner(v,bfield)
  println("op")
  @time op = LinearFEOperator(a,b,V,U,assem,trian,quad)

  # Define the FESolver
  ls = LUSolver()
  solver = LinearFESolver(ls)

  # Solve!
  println("uh")
  @time uh = solve(solver,op)

  # Define exact solution and error
  println("u")
  @time u = CellField(trian,ufun)

  println("e")
  @time e = u - uh

  # Compute errors
  println("el2")
  @time el2 = sqrt(sum( integrate(l2(e),trian,quad) ))

  println("eh1")
  @time eh1 = sqrt(sum( integrate(h1(e),trian,quad) ))

  @assert el2 < 1.e-8
  @assert eh1 < 1.e-8

end


text:
################################################################################

"""
    setup(v::Array{T, 2}) where <: Number

Initiate hidden Markov model object.

See also: [`process!`](@ref), [`HMM`](@ref), [`HMMParams`](@ref)
"""
function setup(v::Array{T, 2}) where T <: Number
  hmm = HMM([zeros(size(v, 2))], [zeros(size(v, 1) + 1)])
  reset!(hmm)
  return hmm
end

################################################################################

max_stars_repo_name: DanielRivasMD/HiddenMarkovModelReaders

max_issues_repo_name: DanielRivasMD/HiddenMarkovModelReaders

max_forks_repo_name: DanielRivasMD/HiddenMarkovModelReaders


text:
module Dec3Part1
slope = [3,1]

function get_level(pattern, level)
    return pattern[level]
end

function get_x_value(level, index)
    n_l = length(level)
    if index <= n_l
        return index
    else
        return index - n_l
    end
end
function get_next_point(current_point, slope, pattern)
    x_ref, y_ref = current_point + slope
    n_level = y_ref
    level = get_level(pattern, n_level)
    x = get_x_value(level, x_ref)
    y = y_ref
    return [x,y]
end

function get_marker_point(pattern, current_point)
    row, col = current_point
    return pattern[col][row]
end

function get_n_trees(pattern)
    global slope
    initial_point = [1,1]
    n_max = length(pattern)
    current_point = initial_point
    counter = 0
    while current_point[2] < n_max
        current_point = get_next_point(current_point, slope, pattern)
        marker = get_marker_point(pattern, current_point)
        if marker == '#'
            counter += 1
        end
    end
    return counter
end
function main(file_path::String)
    lines = readlines(file_path)
    return get_n_trees(lines)
end

end


text:
# TODO: when Optim supports sparse arrays, make a SparseMatrixCSC version of jacobianx

abstract type AbstractBarrierState end

# These are used not only for the current state, but also for the step and the gradient
struct BarrierStateVars{T}
    slack_x::Vector{T}     # values of slack variables for x
    slack_c::Vector{T}     # values of slack variables for c
    λx::Vector{T}          # λ for equality constraints on slack_x
    λc::Vector{T}          # λ for equality constraints on slack_c
    λxE::Vector{T}         # λ for equality constraints on x
    λcE::Vector{T}         # λ for linear/nonlinear equality constraints
end
# Note on λxE:
# We could just set equality-constrained variables to their
# constraint values at the beginning of optimization, but this
# might make the initial guess infeasible in terms of its
# inequality constraints. This would be a much bigger problem than
# not matching the equality constraints.  So we allow them to
# differ, and require that the algorithm can cope with it.

function BarrierStateVars{T}(bounds::ConstraintBounds) where T
    slack_x = Array{T}(undef, length(bounds.ineqx))
    slack_c = Array{T}(undef, length(bounds.ineqc))
    λx = similar(slack_x)
    λc = similar(slack_c)
    λxE = Array{T}(undef, length(bounds.eqx))
    λcE = Array{T}(undef, length(bounds.eqc))
    sv = BarrierStateVars{T}(slack_x, slack_c, λx, λc, λxE, λcE)
end
BarrierStateVars(bounds::ConstraintBounds{T}) where T = BarrierStateVars{T}(bounds)

function BarrierStateVars(bounds::ConstraintBounds{T}, x) where T
    sv = BarrierStateVars(bounds)
    setslack!(sv.slack_x, x, bounds.ineqx, bounds.σx, bounds.bx)
    sv
end
function BarrierStateVars(bounds::ConstraintBounds{T}, x, c) where T
    sv = BarrierStateVars(bounds)
    setslack!(sv.slack_x, x, bounds.ineqx, bounds.σx, bounds.bx)
    setslack!(sv.slack_c, c, bounds.ineqc, bounds.σc, bounds.bc)
    sv
end
function setslack!(slack, v, ineq, σ, b)
    for i = 1:length(ineq)
        dv = v[ineq[i]]-b[i]
        slack[i] = abs(σ[i]*dv)
    end
    slack
end

slack(bstate::BarrierStateVars) = [bstate.slack_x; bstate.slack_c]
lambdaI(bstate::BarrierStateVars) = [bstate.λx; bstate.λc]
lambdaE(bstate::BarrierStateVars) = [bstate.λxE; bstate.λcE] # TODO: Not used by IPNewton?
lambdaI(state::AbstractBarrierState) = lambdaI(state.bstate)
lambdaE(state::AbstractBarrierState) = lambdaE(state.bstate) # TODO: Not used by IPNewton?

Base.similar(bstate::BarrierStateVars) =
    BarrierStateVars(similar(bstate.slack_x),
                     similar(bstate.slack_c),
                     similar(bstate.λx),
                     similar(bstate.λc),
                     similar(bstate.λxE),
                     similar(bstate.λcE))

Base.copy(bstate::BarrierStateVars) =
    BarrierStateVars(copy(bstate.slack_x),
                     copy(bstate.slack_c),
                     copy(bstate.λx),
                     copy(bstate.λc),
                     copy(bstate.λxE),
                     copy(bstate.λcE))

function Base.fill!(b::BarrierStateVars, val)
    fill!(b.slack_x, val)
    fill!(b.slack_c, val)
    fill!(b.λx, val)
    fill!(b.λc, val)
    fill!(b.λxE, val)
    fill!(b.λcE, val)
    b
end

Base.convert(::Type{BarrierStateVars{T}}, bstate::BarrierStateVars) where T =
    BarrierStateVars(convert(Array{T}, bstate.slack_x),
                     convert(Array{T}, bstate.slack_c),
                     convert(Array{T}, bstate.λx),
                     convert(Array{T}, bstate.λc),
                     convert(Array{T}, bstate.λxE),
                     convert(Array{T}, bstate.λcE))

Base.isempty(bstate::BarrierStateVars) = isempty(bstate.slack_x) &
    isempty(bstate.slack_c) & isempty(bstate.λxE) & isempty(bstate.λcE)

Base.eltype(::Type{BarrierStateVars{T}}) where T = T
Base.eltype(sv::BarrierStateVars) = eltype(typeof(sv))

function Base.show(io::IO, b::BarrierStateVars)
    print(io, "BarrierStateVars{$(eltype(b))}:")
    for fn in (:slack_x, :slack_c, :λx, :λc, :λxE, :λcE)
        print(io, "\n  $fn: ")
        show(io, getfield(b, fn))
    end
end

Base.:(==)(v::BarrierStateVars, w::BarrierStateVars) =
    v.slack_x == w.slack_x &&
    v.slack_c == w.slack_c &&
    v.λx == w.λx &&
    v.λc == w.λc &&
    v.λxE == w.λxE &&
    v.λcE == w.λcE

const bsv_seed = sizeof(UInt) == 64 ? 0x145b788192d1cde3 : 0x766a2810
Base.hash(b::BarrierStateVars, u::UInt) =
    hash(b.λcE, hash(b.λxE, hash(b.λc, hash(b.λx, hash(b.slack_c, hash(b.slack_x, u+bsv_seed))))))

function dot(v::BarrierStateVars, w::BarrierStateVars)
    dot(v.slack_x,w.slack_x) +
        dot(v.slack_c, w.slack_c) +
        dot(v.λx, w.λx) +
        dot(v.λc, w.λc) +
        dot(v.λxE, w.λxE) +
        dot(v.λcE, w.λcE)
end

function norm(b::BarrierStateVars, p::Real)
    norm(b.slack_x, p) + norm(b.slack_c, p) +
        norm(b.λx, p) + norm(b.λc, p) +
        norm(b.λxE, p) + norm(b.λcE, p)
end

"""
    BarrierLineSearch{T}

Parameters for interior-point line search methods that use only the value
"""
struct BarrierLineSearch{T}
    c::Vector{T}                  # value of constraints-functions at trial point
    bstate::BarrierStateVars{T}   # trial point for slack and λ variables
end
Base.convert(::Type{BarrierLineSearch{T}}, bsl::BarrierLineSearch) where T =
    BarrierLineSearch(convert(Vector{T}, bsl.c),
                      convert(BarrierStateVars{T}, bsl.bstate))

"""
    BarrierLineSearchGrad{T}

Parameters for interior-point line search methods that exploit the slope.
"""
struct BarrierLineSearchGrad{T}
    c::Vector{T}                  # value of constraints-functions at trial point
    J::Matrix{T}                  # constraints-Jacobian at trial point
    bstate::BarrierStateVars{T}   # trial point for slack and λ variables
    bgrad::BarrierStateVars{T}    # trial point's gradient
end
Base.convert(::Type{BarrierLineSearchGrad{T}}, bsl::BarrierLineSearchGrad) where T =
    BarrierLineSearchGrad(convert(Vector{T}, bsl.c),
                          convert(Matrix{T}, bsl.J),
                          convert(BarrierStateVars{T}, bsl.bstate),
                          convert(BarrierStateVars{T}, bsl.bgrad))

function ls_update!(out::BarrierStateVars, base::BarrierStateVars, step::BarrierStateVars, αs::NTuple{4,Number})
    ls_update!(out.slack_x, base.slack_x, step.slack_x, αs[2])
    ls_update!(out.slack_c, base.slack_c, step.slack_c, αs[2])
    ls_update!(out.λx, base.λx, step.λx, αs[3])
    ls_update!(out.λc, base.λc, step.λc, αs[3])
    ls_update!(out.λxE, base.λxE, step.λxE, αs[4])
    ls_update!(out.λcE, base.λcE, step.λcE, αs[4])
    out
end
ls_update!(out::BarrierStateVars, base::BarrierStateVars, step::BarrierStateVars, αs::Tuple{Number,Number}) =
    ls_update!(out, base, step, (αs[1],αs[1],αs[2],αs[1]))
ls_update!(out::BarrierStateVars, base::BarrierStateVars, step::BarrierStateVars, α::Number) =
    ls_update!(out, base, step, (α,α,α,α))
ls_update!(out::BarrierStateVars, base::BarrierStateVars, step::BarrierStateVars, αs::AbstractVector) =
    ls_update!(out, base, step, αs[1]) # (αs...,))

function initial_convergence(d, state, method::ConstrainedOptimizer, initial_x, options)
    # TODO: Make sure state.bgrad has been evaluated at initial_x
    # state.bgrad normally comes from constraints.c!(..., initial_x) in initial_state
    gradient!(d, initial_x)
    stopped = !isfinite(value(d)) || any(!isfinite, gradient(d))
    norm(gradient(d), Inf) + norm(state.bgrad, Inf) < options.g_abstol, stopped
end

function optimize(f, g, lower::AbstractArray, upper::AbstractArray, initial_x::AbstractArray, method::ConstrainedOptimizer=IPNewton(),
    options::Options = Options(;default_options(method)...))
    d = TwiceDifferentiable(f, g, initial_x)
    optimize(d, lower, upper, initial_x, method, options)
end
function optimize(f, g, h, lower::AbstractArray, upper::AbstractArray, initial_x::AbstractArray, method::ConstrainedOptimizer=IPNewton(),
    options::Options = Options(;default_options(method)...))
    d = TwiceDifferentiable(f, g, h, initial_x)
    optimize(d, lower, upper, initial_x, method, options)
end
function optimize(d::TwiceDifferentiable, lower::AbstractArray, upper::AbstractArray, initial_x::AbstractArray,
    options::Options = Options(;default_options(IPNewton())...))
    optimize(d, lower, upper, initial_x, IPNewton(), options)
end

function optimize(d, lower::AbstractArray, upper::AbstractArray, initial_x::AbstractArray, method::ConstrainedOptimizer,
                  options::Options = Options(;default_options(method)...))
    twicediffed = d isa TwiceDifferentiable ? d : TwiceDifferentiable(d, initial_x)

    bounds = ConstraintBounds(lower, upper, [], [])
    constraints = TwiceDifferentiableConstraints(
            (c,x)->nothing, (J,x)->nothing, (H,x,λ)->nothing, bounds)

    state = initial_state(method, options, twicediffed, constraints, initial_x)

    optimize(twicediffed,
             constraints,
             initial_x,
             method,
             options,
             state)
end
function optimize(d::AbstractObjective, constraints::AbstractConstraints, initial_x::AbstractArray, method::ConstrainedOptimizer,
                  options::Options = Options(;default_options(method)...),
                  state = initial_state(method, options, d, constraints, initial_x))
    #== TODO:
    Let's try to unify this with the unconstrained `optimize` in Optim
    The only thing we'd have to deal with is to dispatch
    the univariate `optimize` to one with empty constraints::AbstractConstraints
    ==#

    t0 = time() # Initial time stamp used to control early stopping by options.time_limit
    _time = t0

    tr = OptimizationTrace{typeof(value(d)), typeof(method)}()
    tracing = options.store_trace || options.show_trace || options.extended_trace || options.callback !== nothing
    stopped, stopped_by_callback, stopped_by_time_limit = false, false, false
    f_limit_reached, g_limit_reached, h_limit_reached = false, false, false
    x_converged, f_converged, f_increased, counter_f_tol = false, false, false, 0

    g_converged, stopped = initial_convergence(d, state, method, initial_x, options)
    converged = g_converged

    # prepare iteration counter (used to make "initial state" trace entry)
    iteration = 0

    options.show_trace && print_header(method)
    trace!(tr, d, state, iteration, method, options, t0)

    while !converged && !stopped && iteration < options.iterations
        iteration += 1

        update_state!(d, constraints, state, method, options) && break # it returns true if it's forced by something in update! to stop (eg dx_dg == 0.0 in BFGS or linesearch errors)

        update_fg!(d, constraints, state, method)

        # TODO: Do we need to rethink f_increased for `ConstrainedOptimizer`s?
        x_converged, f_converged,
        g_converged, f_increased = assess_convergence(state, d, options)
        # With equality constraints, optimization is not necessarily
        # monotonic in the value of the function. If the function
        # change is approximately canceled by a change in the equality
        # violation, it's possible to spuriously satisfy the f_tol
        # criterion. Consequently, we require that the f_tol condition
        # be satisfied a certain number of times in a row before
        # declaring convergence.
        counter_f_tol = f_converged ? counter_f_tol+1 : 0
        converged = x_converged || g_converged || (counter_f_tol > options.successive_f_tol)

        # We don't use the Hessian for anything if we have declared convergence,
        # so we might as well not make the (expensive) update if converged == true
        !converged && update_h!(d, constraints, state, method)

        if tracing
            # update trace; callbacks can stop routine early by returning true
            stopped_by_callback = trace!(tr, d, state, iteration, method, options)
        end

        # Check time_limit; if none is provided it is NaN and the comparison
        # will always return false.
        _time = time()
        stopped_by_time_limit = _time-t0 > options.time_limit ? true : false
        f_limit_reached = options.f_calls_limit > 0 && f_calls(d) >= options.f_calls_limit ? true : false
        g_limit_reached = options.g_calls_limit > 0 && g_calls(d) >= options.g_calls_limit ? true : false
        h_limit_reached = options.h_calls_limit > 0 && h_calls(d) >= options.h_calls_limit ? true : false

        if (f_increased && !options.allow_f_increases) || stopped_by_callback ||
            stopped_by_time_limit || f_limit_reached || g_limit_reached || h_limit_reached
            stopped = true
        end
    end # while

    after_while!(d, constraints, state, method, options)

    # we can just check minimum, as we've earlier enforced same types/eltypes
    # in variables besides the option settings
    T = typeof(options.f_reltol)
    Tf = typeof(value(d))
    f_incr_pick = f_increased && !options.allow_f_increases
    return MultivariateOptimizationResults(method,
                                        initial_x,
                                        pick_best_x(f_incr_pick, state),
                                        pick_best_f(f_incr_pick, state, d),
                                        iteration,
                                        iteration == options.iterations,
                                        x_converged,
                                        T(options.x_abstol),
                                        T(options.x_reltol),
                                        x_abschange(state),
                                        x_relchange(state),
                                        f_converged,
                                        T(options.f_abstol),
                                        T(options.f_reltol),
                                        f_abschange(d, state),
                                        f_relchange(d, state),
                                        g_converged,
                                        T(options.g_abstol),
                                        g_residual(d),
                                        f_increased,
                                        tr,
                                        f_calls(d),
                                        g_calls(d),
                                        h_calls(d),
                                        nothing,
                                        options.time_limit,
                                        _time-t0,
                                        NamedTuple())
end

# Fallbacks (for methods that don't need these)
after_while!(d, constraints::AbstractConstraints, state, method, options) = nothing
update_h!(d, constraints::AbstractConstraints, state, method) = nothing

"""
    initialize_μ_λ!(state, bounds, μ0=:auto, β=0.01)
    initialize_μ_λ!(state, bounds, (Hobj,HcI), μ0=:auto, β=0.01)

Pick μ and λ to ensure that the equality constraints are satisfied
locally (at the current `state.x`), and that the initial gradient
including the barrier would be a descent direction for the problem
without the barrier (μ = 0). This ensures that the search isn't pushed
out of the basin of the user-supplied initial guess.

Upon entry, the objective function gradient, constraint values, and
constraint jacobian must be set in `state.g`, `state.c`, and `state.J`
respectively. If you also wish to ensure that the projection of
Hessian is minimally-perturbed along the initial gradient, supply the
hessian of the objective (`Hobj`) and

    HcI = ∑_i (σ_i/s_i)∇∇ c_{Ii}

for the constraints. This can be obtained as

    HcI = hessianI(state.x, constraints, 1./state.slack_c)

You can manually specify `μ` by supplying a numerical value for
`μ0`. Whether calculated algorithmically or specified manually, the
values of `λ` are set using the chosen `μ`.
"""
function initialize_μ_λ!(state, bounds::ConstraintBounds, Hinfo, μ0::Union{Symbol,Number}, β::Number=1//100)
    if nconstraints(bounds) == 0 && nconstraints_x(bounds) == 0
        state.μ = 0
        fill!(state.bstate, 0)
        return state
    end
    gf = state.g  # must be pre-set to ∇f
    # Calculate projection of ∇f into the subspace spanned by the
    # equality constraint Jacobian
    JE = jacobianE(state, bounds)
    # QRF = qrfact(JE)
    # Q = QRF[:Q]
    # PEg = Q'*(Q*gf)   # in the subspace of JE
    C = JE*JE'
    Cc = cholesky(Positive, C)
    Pperpg = gf-JE'*(Cc \ (JE*gf))   # in the nullspace of JE
    # Set μ
    JI = jacobianI(state, bounds)
    if μ0 == :auto
        # Calculate projections of the Lagrangian's gradient, and
        # possibly hessian, along (∇f)_⟂
        Dperp = dot(Pperpg, Pperpg)
        σ, s = sigma(bounds), slack(state)
        σdivs = σ./s
        Δg = JI'*σdivs
        PperpΔg = Δg - JE'*(Cc \ (JE*Δg))
        DI = dot(PperpΔg, PperpΔg)
        κperp, κI = hessian_projections(Hinfo, Pperpg, (JI*Pperpg)./s)
        # Calculate μ and λI
        μ = β * (κperp == 0 ? sqrt(Dperp/DI) : min(sqrt(Dperp/DI), abs(κperp/κI)))
        if !isfinite(μ)
            Δgtilde = JI'*(1 ./ s)
            PperpΔgtilde = Δgtilde - JE'*(Cc \ (JE*Δgtilde))
            DItilde = dot(PperpΔgtilde, PperpΔgtilde)
            μ = β*sqrt(Dperp/DItilde)
        end
        if !isfinite(μ) || μ == 0
            μ = one(μ)
        end
    else
        μ = convert(eltype(state.x), μ0)
    end
    state.μ = μ
    # Set λI
    @. state.bstate.λx = μ / state.bstate.slack_x
    @. state.bstate.λc = μ / state.bstate.slack_c
    # Calculate λE
    λI = lambdaI(state)
    ∇bI = gf - JI'*λI
#    qrregularize!(QRF)  # in case of any 0 eigenvalues
    λE = Cc \ (JE*∇bI) + (cbar(bounds) - cE(state, bounds))/μ
    k = unpack_vec!(state.bstate.λxE, λE, 0)
    k = unpack_vec!(state.bstate.λcE, λE, k)
    k == length(λE) || error("Something is wrong when initializing μ and λ.")
    state
end
function initialize_μ_λ!(state, bounds::ConstraintBounds, μ0::Union{Number,Symbol}, β::Number=1//100)
    initialize_μ_λ!(state, bounds, nothing, μ0, β)
end

function hessian_projections(Hinfo::Tuple{AbstractMatrix,AbstractMatrix}, Pperpg, y)
    κperp = dot(Hinfo[1]*Pperpg, Pperpg)
    κI = dot(Hinfo[2]*Pperpg, Pperpg) + dot(y,y)
    κperp, κI
end
hessian_projections(Hinfo::Nothing, Pperpg::AbstractVector{T}) where T = convert(T, Inf), zero(T)

function jacobianE(state, bounds::ConstraintBounds)
    J, x = state.constr_J, state.x
    JEx = jacobianx(J, bounds.eqx)
    JEc = view(J, bounds.eqc, :)
    JE = vcat(JEx, JEc)
end
jacobianE(state, constraints) = jacobianE(state, constraints.bounds)

function jacobianI(state, bounds::ConstraintBounds)
    J, x = state.constr_J, state.x
    JIx = jacobianx(J, bounds.ineqx)
    JIc = view(J, bounds.ineqc, :)
    JI = vcat(JIx, JIc)
end
jacobianI(state, constraints) = jacobianI(state, constraints.bounds)

# TODO: when Optim supports sparse arrays, make a SparseMatrixCSC version
function jacobianx(J::AbstractArray, indx)
    Jx = zeros(eltype(J), length(indx), size(J, 2))
    for (i,j) in enumerate(indx)
        Jx[i,j] = 1
    end
    Jx
end

function sigma(bounds::ConstraintBounds)
    [bounds.σx; bounds.σc]  # don't include σz
end
sigma(constraints) = sigma(constraints.bounds)

slack(state) = slack(state.bstate)

cbar(bounds::ConstraintBounds) = [bounds.valx; bounds.valc]
cbar(constraints) = cbar(constraints.bounds)
cE(state, bounds::ConstraintBounds) = [state.x[bounds.eqx]; state.constr_c[bounds.eqc]]

function hessianI!(h, x, constraints, λcI, μ)
    λ = userλ(λcI, constraints)
    constraints.h!(h, x, λ)
    h
end

"""
   hessianI(x, constraints, λcI, μ) -> h

Compute the hessian at `x` of the `λcI`-weighted sum of user-supplied
constraint functions for just the inequalities.  This also includes
contributions from any variables with bounds at 0, since those do not
cause introduction of a slack variable. Other (nonzero) box
constraints do not contribute to `h`, because the hessian of `x_i` is
zero. (They contribute indirectly via their slack variables.)
"""
hessianI(x, constraints, λcI, μ) =
    hessianI!(zeros(eltype(x), length(x), length(x)), x, constraints, λcI, μ)

"""
    userλ(λcI, bounds) -> λ

Accumulates `λcI` into a vector `λ` ordered as the user-supplied
constraint functions `c`. Upper and lower bounds are summed, weighted
by `σ`. The resulting λ includes an overall negative sign so that this
becomes the coefficient for the user-supplied hessian.

This is relevant only for the inequalities. If you want the λ for just
the equalities, you can use `λ[bounds.ceq] = λcE` for a zero-filled `λ`.
"""
function userλ(λcI, bounds::ConstraintBounds)
    ineqc, σc = bounds.ineqc, bounds.σc
    λ = zeros(eltype(bounds), nconstraints(bounds))
    for i = 1:length(ineqc)
        λ[ineqc[i]] -= λcI[i]*σc[i]
    end
    λ
end
userλ(λcI, constraints) = userλ(λcI, constraints.bounds)

## Computation of the Lagrangian and its gradient
# This is in a parametrization that is also useful during linesearch
# TODO: `lagrangian` does not seem to be used (IPNewton)?
function lagrangian(d, bounds::ConstraintBounds, x, c, bstate::BarrierStateVars, μ)
    f_x = NLSolversBase.value!(d, x)
    ev = equality_violation(bounds, x, c, bstate)
    L_xsλ = f_x + barrier_value(bounds, x, bstate, μ) + ev
    f_x, L_xsλ, ev
end

function lagrangian_fg!(gx, bgrad, d, bounds::ConstraintBounds, x, c, J, bstate::BarrierStateVars, μ)
    fill!(bgrad, 0)
    f_x, g_x = NLSolversBase.value_gradient!(d,x)
    gx .= g_x
    ev = equality_violation(bounds, x, c, bstate)
    L_xsλ = f_x + barrier_value(bounds, x, bstate, μ) + ev
    barrier_grad!(bgrad, bounds, x, bstate, μ)
    equality_grad!(gx, bgrad, bounds, x, c, J, bstate)
    f_x, L_xsλ, ev
end

# TODO: do we need lagrangian_vec? Maybe for automatic differentiation?
## Computation of Lagrangian and derivatives when passing all parameters as a single vector
function lagrangian_vec(p, d, bounds::ConstraintBounds, x, c::AbstractArray, bstate::BarrierStateVars, μ)
    unpack_vec!(x, bstate, p)
    f_x, L_xsλ, ev = lagrangian(d, bounds, x, c, bstate, μ)
    L_xsλ
end
function lagrangian_vec(p, d, bounds::ConstraintBounds, x, c::Function, bstate::BarrierStateVars, μ)
    # Use this version when using automatic differentiation
    unpack_vec!(x, bstate, p)
    f_x, L_xsλ, ev = lagrangian(d, bounds, x, c(x), bstate, μ)
    L_xsλ
end
function lagrangian_fgvec!(p, storage, gx, bgrad, d, bounds::ConstraintBounds, x, c, J, bstate::BarrierStateVars, μ)
    unpack_vec!(x, bstate, p)
    f_x, L_xsλ, ev = lagrangian_fg!(gx, bgrad, d, bounds, x, c, J, bstate, μ)
    pack_vec!(storage, gx, bgrad)
    L_xsλ
end

## for line searches that don't use the gradient along the line
function lagrangian_linefunc(αs, d, constraints, state)
    _lagrangian_linefunc(αs, d, constraints, state)[2]
end

function _lagrangian_linefunc(αs, d, constraints, state)
    b_ls, bounds = state.b_ls, constraints.bounds
    ls_update!(state.x_ls, state.x, state.s, alphax(αs))
    ls_update!(b_ls.bstate, state.bstate, state.bstep, αs)
    constraints.c!(b_ls.c, state.x_ls)
    lagrangian(d, constraints.bounds, state.x_ls, b_ls.c, b_ls.bstate, state.μ)
end
alphax(α::Number) = α
alphax(αs::Union{Tuple,AbstractVector}) = αs[1]

function lagrangian_linefunc!(α, d, constraints, state, method::IPOptimizer{typeof(backtrack_constrained)})
    # For backtrack_constrained, the last evaluation is the one we
    # keep, so it's safe to store the results in state
    state.f_x, state.L, state.ev = _lagrangian_linefunc(α, d, constraints, state)
    state.L
end
lagrangian_linefunc!(α, d, constraints, state, method) = lagrangian_linefunc(α, d, constraints, state)

## for line searches that do use the gradient along the line
function lagrangian_lineslope(αs, d, constraints, state)
    f_x, L, ev, slope = _lagrangian_lineslope(αs, d, constraints, state)
    L, slope
end

function _lagrangian_lineslope(αs, d, constraints, state)
    b_ls, bounds = state.b_ls, constraints.bounds
    bstep, bgrad = state.bstep, b_ls.bgrad
    ls_update!(state.x_ls, state.x, state.s, alphax(αs))
    ls_update!(b_ls.bstate, state.bstate, bstep, αs)
    constraints.c!(b_ls.c, state.x_ls)
    constraints.jacobian!(b_ls.J, state.x_ls)
    f_x, L, ev = lagrangian_fg!(state.g, bgrad, d, bounds, state.x_ls, b_ls.c, b_ls.J, b_ls.bstate, state.μ)
    slopeα = slopealpha(state.s, state.g, bstep, bgrad)
    f_x, L, ev, slopeα
end

function lagrangian_lineslope!(αs, d, constraints, state, method::IPOptimizer{typeof(backtrack_constrained_grad)})
    # For backtrack_constrained, the last evaluation is the one we
    # keep, so it's safe to store the results in state
    state.f_x, state.L, state.ev, slope = _lagrangian_lineslope(αs, d, constraints, state)
     state.L, slope
end
lagrangian_lineslope!(αs, d, constraints, state, method) = lagrangian_lineslope(αs, d, constraints, state)

slopealpha(sx, gx, bstep, bgrad) = dot(sx, gx) +
    dot(bstep.slack_x, bgrad.slack_x) + dot(bstep.slack_c, bgrad.slack_c) +
    dot(bstep.λx, bgrad.λx) + dot(bstep.λc, bgrad.λc) +
    dot(bstep.λxE, bgrad.λxE) + dot(bstep.λcE, bgrad.λcE)

function linesearch_anon(d, constraints, state, method::IPOptimizer{typeof(backtrack_constrained_grad)})
    αs->lagrangian_lineslope!(αs, d, constraints, state, method)
end
function linesearch_anon(d, constraints, state, method::IPOptimizer{typeof(backtrack_constrained)})
    αs->lagrangian_linefunc!(αs, d, constraints, state, method)
end

## Computation of Lagrangian terms: barrier penalty
"""
    barrier_value(constraints, state) -> val
    barrier_value(bounds, x, sx, sc, μ) -> val

Compute the value of the barrier penalty at the current `state`, or at
a position (`x`,`sx`,`sc`), where `x` is the current position, `sx`
are the coordinate slack variables, and `sc` are the linear/nonlinear
slack variables. `bounds` holds the parsed bounds.
"""
function barrier_value(bounds::ConstraintBounds, x, sx, sc, μ)
    # bμ is the coefficient of μ in the barrier penalty
    bμ = _bv(sx) +  # coords with other bounds
         _bv(sc)    # linear/nonlinear constr.
    μ*bμ
end
barrier_value(bounds::ConstraintBounds, x, bstate::BarrierStateVars, μ) =
    barrier_value(bounds, x, bstate.slack_x, bstate.slack_c, μ)
barrier_value(bounds::ConstraintBounds, state) =
    barrier_value(bounds, state.x, state.bstate.slack_x, state.bstate.slack_c, state.μ)
barrier_value(constraints::AbstractConstraints, state) =
    barrier_value(constraints.bounds, state)

# don'tcall this barrier_value because it lacks μ
function _bv(v, idx, σ) # TODO: Not used, delete? (IPNewton)
    ret = loginf(one(eltype(σ))*one(eltype(v)))
    for (i,iv) in enumerate(idx)
        ret += loginf(σ[i]*v[iv])
    end
    -ret
end

_bv(v) = isempty(v) ? loginf(one(eltype(v))) : -sum(loginf, v)

loginf(δ) = δ > 0 ? log(δ) : -oftype(δ, Inf)

"""
    barrier_grad!(bgrad, bounds, x, bstate, μ)
    barrier_grad!(gsx, gsc, bounds, x, sx, sc, μ)

Compute the gradient of the barrier penalty at (`x`,`sx`,`sc`), where
`x` is the current position, `sx` are the coordinate slack variables,
and `sc` are the linear/nonlinear slack
variables. `bounds::ConstraintBounds` holds the parsed bounds.

The result is *added* to `gsx`, and `gsc`, so these vectors
need to be initialized appropriately.
"""
function barrier_grad!(gsx, gsc, bounds::ConstraintBounds, x, sx, sc, μ)
    barrier_grad!(gsx, sx, μ)
    barrier_grad!(gsc, sc, μ)
    nothing
end
barrier_grad!(bgrad, bounds::ConstraintBounds, x, bstate, μ) =
    barrier_grad!(bgrad.slack_x, bgrad.slack_c, bounds, x, bstate.slack_x, bstate.slack_c, μ)

function barrier_grad!(out, v, μ)
    for i = 1:length(out)
        out[i] -= μ/v[i]
    end
    nothing
end

## Computation of Lagrangian terms: equality constraints penalty

"""
    equality_violation([f=identity], bounds, x, c, bstate) -> val
    equality_violation([f=identity], bounds, x, c, sx, sc, λx, λc, λxE, λcE) -> val

Compute the sum of `f(v_i)`, where `v_i = λ_i*(target - observed)`
measures the difference between the current state and the
equality-constrained state. `bounds::ConstraintBounds` holds the
parsed bounds. `x` is the current position, `sx` are the coordinate
slack variables, and `sc` are the linear/nonlinear slack
variables. `c` holds the values of the linear-nonlinear constraints,
and the λ arguments hold the Lagrange multipliers for `x`, `sx`, `sc`, and
`c` respectively.
"""
function equality_violation(f, bounds::ConstraintBounds, x, c, sx, sc, λx, λc, λxE, λcE)
    ev = equality_violation(f, sx, x, bounds.ineqx, bounds.σx, bounds.bx, λx) +
         equality_violation(f, sc, c, bounds.ineqc, bounds.σc, bounds.bc, λc) +
         equality_violation(f, x, bounds.valx, bounds.eqx, λxE) +
         equality_violation(f, c, bounds.valc, bounds.eqc, λcE)
end
equality_violation(bounds::ConstraintBounds, x, c, sx, sc, λx, λc, λxE, λcE) =
    equality_violation(identity, bounds, x, c, sx, sc, λx, λc, λxE, λcE)
function equality_violation(f, bounds::ConstraintBounds, x, c, bstate::BarrierStateVars)
    equality_violation(f, bounds, x, c, bstate.slack_x, bstate.slack_c,
                       bstate.λx, bstate.λc, bstate.λxE, bstate.λcE)
end
equality_violation(bounds::ConstraintBounds, x, c, bstate::BarrierStateVars) =
    equality_violation(identity, bounds, x, c, bstate)
equality_violation(f, bounds::ConstraintBounds, state::AbstractBarrierState) =
    equality_violation(f, bounds, state.x, state.constr_c, state.bstate)
equality_violation(bounds::ConstraintBounds, state::AbstractBarrierState) =
    equality_violation(identity, bounds, state)
equality_violation(f, constraints::AbstractConstraints, state::AbstractBarrierState) =
    equality_violation(f, constraints.bounds, state)
equality_violation(constraints::AbstractConstraints, state::AbstractBarrierState) =
    equality_violation(constraints.bounds, state)

# violations of s = σ*(v-b)
function equality_violation(f, s, v, ineq, σ, b, λ)
    ret = f(zero(eltype(λ))*(zero(eltype(s))-zero(eltype(σ))*(zero(eltype(v))-zero(eltype(b)))))
    for (i,iv) in enumerate(ineq)
        ret += f(λ[i]*(s[i] - σ[i]*(v[iv]-b[i])))
    end
    ret
end

# violations of v = target
function equality_violation(f, v, target, idx, λ)
    ret = f(zero(eltype(λ))*(zero(eltype(v))-zero(eltype(target))))
    for (i,iv) in enumerate(idx)
        ret += f(λ[i]*(target[i] - v[iv]))
    end
    ret
end

"""
    equality_grad!(gx, gbstate, bounds, x, c, J, bstate)

Compute the gradient of `equality_violation`, storing the result in `gx` (an array) and `gbstate::BarrierStateVars`.
"""
function equality_grad!(gx, gsx, gsc, gλx, gλc, gλxE, gλcE, bounds::ConstraintBounds, x, c, J, sx, sc, λx, λc, λxE, λcE)
    equality_grad_var!(gsx, gx, bounds.ineqx, bounds.σx, λx)
    equality_grad_var!(gsc, gx, bounds.ineqc, bounds.σc, λc, J)
    gx[bounds.eqx] .= gx[bounds.eqx] .- λxE
    equality_grad_var!(gx, bounds.eqc, λcE, J)
    equality_grad_λ!(gλx, sx, x, bounds.ineqx, bounds.σx, bounds.bx)
    equality_grad_λ!(gλc, sc, c, bounds.ineqc, bounds.σc, bounds.bc)
    equality_grad_λ!(gλxE, x, bounds.valx, bounds.eqx)
    equality_grad_λ!(gλcE, c, bounds.valc, bounds.eqc)
end
equality_grad!(gx, gb::BarrierStateVars, bounds::ConstraintBounds, x, c, J, b::BarrierStateVars) =
    equality_grad!(gx, gb.slack_x, gb.slack_c, gb.λx, gb.λc, gb.λxE, gb.λcE,
                   bounds, x, c, J,
                   b.slack_x, b.slack_c, b.λx, b.λc, b.λxE, b.λcE)

# violations of s = σ*(x-b)
function equality_grad_var!(gs, gx, ineq, σ, λ)
    for (i,ix) in enumerate(ineq)
        λi = λ[i]
        gs[i] += λi
        gx[ix] -= λi*σ[i]
    end
    nothing
end

function equality_grad_var!(gs, gx, ineq, σ, λ, J)
    @. gs = gs + λ
    if !isempty(ineq)
        gx .= gx .- view(J, ineq, :)'*(λ.*σ)
    end
    nothing
end

function equality_grad_λ!(gλ, s, v, ineq, σ, b)
    for (i,iv) in enumerate(ineq)
        gλ[i] += s[i] - σ[i]*(v[iv]-b[i])
    end
    nothing
end

# violations of v = target
function equality_grad_var!(gx, idx, λ, J)
    if !isempty(idx)
        gx .= gx .- view(J, idx, :)'*λ
    end
    nothing
end

function equality_grad_λ!(gλ, v, target, idx)
    for (i,iv) in enumerate(idx)
        gλ[i] += target[i] - v[iv]
    end
    nothing
end

"""
    isfeasible(constraints, state) -> Bool
    isfeasible(constraints, x, c) -> Bool
    isfeasible(constraints, x) -> Bool
    isfeasible(bounds, x, c) -> Bool

Return `true` if point `x` is feasible, given the `constraints` which
specify bounds `lx`, `ux`, `lc`, and `uc`. `x` is feasible if

    lx[i] <= x[i] <= ux[i]
    lc[i] <= c[i] <= uc[i]

for all possible `i`.
"""
function isfeasible(bounds::ConstraintBounds, x, c)
    isf = true
    for (i,j) in enumerate(bounds.eqx)
        isf &= x[j] == bounds.valx[i]
    end
    for (i,j) in enumerate(bounds.ineqx)
        isf &= bounds.σx[i]*(x[j] - bounds.bx[i]) >= 0
    end
    for (i,j) in enumerate(bounds.eqc)
        isf &= c[j] == bounds.valc[i]
    end
    for (i,j) in enumerate(bounds.ineqc)
        isf &= bounds.σc[i]*(c[j] - bounds.bc[i]) >= 0
    end
    isf
end
isfeasible(constraints, state::AbstractBarrierState) = isfeasible(constraints, state.x, state.constraints_c)
function isfeasible(constraints, x)
    # don't assume c! returns c (which means this is a little more awkward)
    c = Array{eltype(x)}(undef, constraints.bounds.nc)
    constraints.c!(c, x)
    isfeasible(constraints, x, c)
end
isfeasible(constraints::AbstractConstraints, x, c) = isfeasible(constraints.bounds, x, c)
isfeasible(constraints::Nothing, state::AbstractBarrierState) = true
isfeasible(constraints::Nothing, x) = true

"""
    isinterior(constraints, state) -> Bool
    isinterior(constraints, x, c) -> Bool
    isinterior(constraints, x) -> Bool
    isinterior(bounds, x, c) -> Bool

Return `true` if point `x` is on the interior of the allowed region,
given the `constraints` which specify bounds `lx`, `ux`, `lc`, and
`uc`. `x` is in the interior if

    lx[i] < x[i] < ux[i]
    lc[i] < c[i] < uc[i]

for all possible `i`.
"""
function isinterior(bounds::ConstraintBounds, x, c)
    isi = true
    for (i,j) in enumerate(bounds.ineqx)
        isi &= bounds.σx[i]*(x[j] - bounds.bx[i]) > 0
    end
    for (i,j) in enumerate(bounds.ineqc)
        isi &= bounds.σc[i]*(c[j] - bounds.bc[i]) > 0
    end
    isi
end
isinterior(constraints, state::AbstractBarrierState) = isinterior(constraints, state.x, state.constraints_c)
function isinterior(constraints, x)
    c = Array{eltype(x)}(undef, constraints.bounds.nc)
    constraints.c!(c, x)
    isinterior(constraints, x, c)
end
isinterior(constraints::AbstractConstraints, x, c) = isinterior(constraints.bounds, x, c)
isinterior(constraints::Nothing, state::AbstractBarrierState) = true
isinterior(constraints::Nothing, x) = true

## Utilities for representing total state as single vector
# TODO: Most of these seem to be unused (IPNewton)?
function pack_vec(x, b::BarrierStateVars)
    n = length(x)
    for fn in fieldnames(b)
        n += length(getfield(b, fn))
    end
    vec = Array{eltype(x)}(undef, n)
    pack_vec!(vec, x, b)
end

function pack_vec!(vec, x, b::BarrierStateVars)
    k = pack_vec!(vec, x, 0)
    for fn in fieldnames(b)
        k = pack_vec!(vec, getfield(b, fn), k)
    end
    k == length(vec) || throw(DimensionMismatch("vec should have length $k, got $(length(vec))"))
    vec
end
function pack_vec!(vec, x, k::Int)
    for i = 1:length(x)
        vec[k+=1] = x[i]
    end
    k
end
function unpack_vec!(x, b::BarrierStateVars, vec::Vector)
    k = unpack_vec!(x, vec, 0)
    for fn in fieldnames(b)
        k = unpack_vec!(getfield(b, fn), vec, k)
    end
    k == length(vec) || throw(DimensionMismatch("vec should have length $k, got $(length(vec))"))
    x, b
end
function unpack_vec!(x, vec::Vector, k::Int)
    for i = 1:length(x)
        x[i] = vec[k+=1]
    end
    k
end

## More utilities
function estimate_maxstep(αmax, x, s)
    for i = 1:length(s)
        si = s[i]
        if si < 0
            αmax = min(αmax, -x[i]/si)
        end
    end
    αmax
end

# TODO: This is not used anymore??
function qrregularize!(QRF)
    R = QRF[:R]
    for i = 1:size(R, 1)
        if R[i,i] == 0
            R[i,i] = 1
        end
    end
    QRF
end


text:
using NURBS.Utils
using NURBS.Bases

# NURBSBasis inner constructor
# ========================================================================

bs = BSplineBasis(range(0, stop=4, length=5), 3)

basis = NURBSBasis(bs, [1.1, 1.2, 1.3, 1.4, 1.5, 1.6], 1)
@test basis.bs == bs
@test basis.weights == [1.1, 1.2, 1.3, 1.4, 1.5, 1.6]
@test basis.deriv == 1

basis = NURBSBasis(bs)
@test basis.bs == bs
@test basis.weights == ones(6)
@test basis.deriv == 0

basis = NURBSBasis(bs, 2)
@test basis.bs == bs
@test basis.weights == ones(6)
@test basis.deriv == 2

basis = NURBSBasis(bs, range(4, stop=5, length=6))
@test basis.bs == bs
@test basis.weights == collect(range(4, stop=5, length=6))
@test basis.deriv == 0

@test_throws ArgumentError NURBSBasis(deriv(bs))
@test_throws ArgumentError NURBSBasis(bs, range(1, stop=2, length=5))
@test_throws ArgumentError NURBSBasis(bs, range(1, stop=2, length=7))
@test_throws ArgumentError NURBSBasis(bs, [-1.0; ones(5)])
@test_throws ArgumentError NURBSBasis(bs, 3)
@test_throws ArgumentError NURBSBasis(BSplineBasis([0, 1], 2), 2)

# NURBSBasis length, size and domain
# ========================================================================

for bs in [BSplineBasis(0, 1, 4, 4), BSplineBasis(0, 2, 5, 3)]
    nbs = NURBSBasis(bs)
    @test length(nbs) == length(bs)
    @test size(nbs) == size(bs)
    @test domain(nbs) == domain(bs)
end

# NURBSBasis evaluation
# ========================================================================

t = sqrt(2)/2
wts = [1, t, 1, t, 1, t, 1, t, 1]
kts = [0, pi/2, pi/2, pi, pi, 3pi/2, 3pi/2, 2pi]
coeffs = [1 0; 1 1; 0 1; -1 1; -1 0; -1 -1; 0 -1; 1 -1; 1 0]
bs = NURBSBasis(BSplineBasis(kts, 3), wts)

testeval_snn(bs, 0, [1, 0, 0], 1:3)
testeval_snn(bs, 0.5, [0.53237051599, 0.351554065456, 0.116075418554], 1:3)
testeval_snn(bs, 1.0, [0.1527439727714, 0.378440763761, 0.468815263468], 1:3)
testeval_snn(bs, 1.5, [0.0020838736033, 0.062440576626, 0.935475549771], 1:3)
testeval_snn(bs, 2.0, [0.597709490779, 0.317802576968, 0.0844879322529], 3:5)
testeval_snn(bs, 2.5, [0.1943378565224, 0.398037455731, 0.407624687747], 3:5)
testeval_snn(bs, 3.0, [0.00853541363591, 0.1218408648272, 0.869623721537], 3:5)
testeval_snn(bs, 3.5, [0.664247692256, 0.277702618453, 0.0580496892913], 5:7)
testeval_snn(bs, 4.0, [0.2406138000566, 0.41002634575, 0.349359854193], 5:7)
testeval_snn(bs, 4.5, [0.0196263352654, 0.1775220289045, 0.80285163583], 5:7)
testeval_snn(bs, 5.0, [0.73141253905, 0.2318427624256, 0.0367446985247], 7:9)
testeval_snn(bs, 5.5, [0.291244696239, 0.414209710077, 0.294545593684], 7:9)
testeval_snn(bs, 6.0, [0.0355815798651, 0.228799037625, 0.73561938251], 7:9)
testeval_snn(bs, 2pi, [0, 0, 1], 7:9)
norms = sum(bs(collect(range(0, stop=2pi/99, length=100)), coeffs) .^ 2, dims=2)
@test isapprox(norms, ones(100))

let bs = bs
    for j in 1:2
        for i in 0.5:0.5:6.0
            testeval_sdn(bs, i)
        end
        bs = deriv(bs)
    end
    @test_throws ArgumentError deriv(bs)
end


text:
#=
   Tourist with a short memory

   From
   Gunnar Blom, Lars Holst, Dennis Sandell:
   "Problems and Snapshots from the World of Probability"
   page 2, Problem 1.2

   A tourist wants to visits all the four cities, A, B, C, and D.
   When he is in A he then visit the other cities (B, C, or D)
   with equal probability, and the same for the other cities.
   Unfortunately, he don't remember which cities he has already visited
   so he always select the other 3 cities woth equal probability.

   How many cities will he visit until he had visited all four cities?

   The text state the "theoretical" expectation (geometry distribution) as:
   """
      E(N) = 1 + 1 + 3/2 + 3 = 13/2 [= 6.5]
   """

   Theoretical (numCities: 4): 6.5

   Summary Statistics
    parameters      mean       std   naive_se      mcse         ess      rhat   ess_per_sec
        Symbol   Float64   Float64    Float64   Float64     Float64   Float64       Float64

           len    6.4313    2.5985     0.0260    0.0293   6805.0938    1.0000     2741.7783

   Distributions of variable len (num:0)
   4.00000 =>    2315  (0.231500)
   5.00000 =>    2258  (0.225800)
   6.00000 =>    1749  (0.174900)
   7.00000 =>    1197  (0.119700)
   8.00000 =>     842  (0.084200)
   9.00000 =>     570  (0.057000)
   10.00000 =>     324  (0.032400)
   11.00000 =>     241  (0.024100)
   12.00000 =>     167  (0.016700)
   13.00000 =>     107  (0.010700)
   14.00000 =>      73  (0.007300)
   15.00000 =>      47  (0.004700)
   16.00000 =>      28  (0.002800)
   17.00000 =>      25  (0.002500)
   18.00000 =>      20  (0.002000)
   20.00000 =>       9  (0.000900)
   21.00000 =>       9  (0.000900)
   19.00000 =>       8  (0.000800)
   22.00000 =>       7  (0.000700)
   28.00000 =>       2  (0.000200)
   24.00000 =>       1  (0.000100)
   23.00000 =>       1  (0.000100)

  Cf ~/webppl/tourist_with_a_short_memory.wppl

=#
using Turing, StatsPlots, Distributions
include("jl_utils.jl")

@model function tourist_with_a_short_memory(numCities=4)
    cities = [i for i in 1:numCities]

    function visitCities(a)
        len = length(a)
        if length(unique(a)) == numCities
            return a
        elseif len == 0
            c = rand(cities)
            return visitCities(push!(a,c))
        else
            lastCity = a[end]
            nextCity = rand(setdiff(cities,[lastCity]))
            return visitCities(push!(a,nextCity))
        end
    end

    a = visitCities([])
    len ~ Dirac(length(a))

end

numCities = 4
theoretical_tourist = numCities -> 1 + (numCities-1)*sum([1/i for i in 1:numCities-1])
println("Theoretical (numCities: $numCities): ", theoretical_tourist(numCities))
model = tourist_with_a_short_memory(numCities)

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 10_000)
chns = sample(model, PG(5), 10_000)
# chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 10_000)

display(chns)
# display(plot(chns))

show_var_dist_pct(chns,:len)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
using POMCP
using POMDPs
using POMDPModels
using POMDPToolbox

problem = BabyPOMDP();
rng = MersenneTwister(1);

solver = POMCPSolver(rng=rng, tree_queries=5)
policy = solve(solver, problem)
up = updater(policy)

# setup
init_dist = initial_state_distribution(problem)
s = rand(rng, init_dist)
first_root_node = initialize_belief(up, init_dist)
typeof(first_root_node)

# plan and execute first action
a = action(policy, first_root_node)
(sp, o, r) = generate_sor(problem, s, a, rng)

# the updater simply chooses the next root node
second_root_node = update(up, first_root_node, a, o)
typeof(second_root_node)

# this new node contains particles representing the belief
second_root_node.B

# at the next step, POMCP uses the new root
action(policy, second_root_node);

# artificially simulate particle depletion
delete!(first_root_node.children[a].children, o)
@test_throws ErrorException update(up, first_root_node, a, o)

type UniformBabyReinvigorator <: ParticleReinvigorator end

function POMCP.reinvigorate!(pc::ParticleCollection,
        r::UniformBabyReinvigorator,
        old_node::BeliefNode, a::Bool, o::Bool)
    push!(pc, true)
    push!(pc, false)
    return pc
end

function POMCP.handle_unseen_observation(r::UniformBabyReinvigorator,
        old_node::BeliefNode, a::Bool, o::Bool)
    return ParticleCollection{Bool}([true, false])
end

up_with_reinvig = RootUpdater(UniformBabyReinvigorator())
# artificially simulate particle depletion
delete!(first_root_node.children[a].children, o)
update(up_with_reinvig, first_root_node, a, o)

init_dist = initial_state_distribution(problem)
action(policy, init_dist)

get(policy._tree_ref).children[true].children[false].B

solver = POMCPSolver(rng=rng, tree_queries=5,
                     node_belief_updater = VoidUpdater())
policy = solve(solver, problem)
a = action(policy, init_dist)

exact_updater = BabyBeliefUpdater(problem)
(sp, o, r) = generate_sor(problem, s, a, rng)
belief2 = update(exact_updater, init_dist, a, o)
a2 = action(policy, belief2)

solver = POMCPSolver(rng=rng, tree_queries=5,
                     node_belief_updater = exact_updater)
policy = solve(solver, problem)
a = action(policy, init_dist)
get(policy._tree_ref).children[true].children[false].B

POMCP.uses_states_from_planner(::BoolDistribution) = true
Base.push!(::BoolDistribution, s) = println("Received state $s from planner.")

solver = POMCPSolver(rng=rng, tree_queries=5,
                     node_belief_updater = exact_updater)
policy = solve(solver, problem)
a = action(policy, init_dist)


text:
using Talkon
using ConfigEnv
using MiniLoggers

MiniLogger(minlevel = MiniLoggers.Debug,
           errlevel = MiniLoggers.AboveMaxLevel,
           message_mode = :squash,
           format = "[{timestamp:func}] {level:func} {basename}:{line:cyan}{::yellow} {message}") |> global_logger

dotenv()

db = initialize(get(ENV, "TALKON_DATA_FILE", "varTEST.data"))
talka(db)


text:
## destruct.jl - deconstruction of structs
##
## TODO: now Espresso supports structures natively, remove this?

"Check if an object is of a struct type, i.e. not a number or array"
isstruct(::Type{T}) where T = !isbits(T) && !(T <: AbstractArray)
isstruct(obj) = !isbits(obj) && !isa(obj, AbstractArray)

field_values(m) = [getfield(m, f) for f in fieldnames(typeof(m))]
named_field_values(m) = [f => getfield(m, f) for f in fieldnames(typeof(m))]

"""
Replace all struct arguments by a list of their plain analogues.
Example:

    args = [:m, :x, :y]
    types = (Linear, Matrix{Float64}, Matrix{Float64})
    ex = :(sum((m.W * x .+ m.b) - y))

    destruct(args, types, ex)
    # ==>
    # ([:m_W, :m_b, :x, :y],
    #  :(sum((m_W * x .+ m_b) - y)),
    #  Dict(:(m.W) => :m_W, :(m.b) => :m_b))
"""
function destruct(args::Vector{Symbol}, types, ex)
    st = Dict()
    new_args = Symbol[]
    for (arg, T) in zip(args, types)
        if isstruct(T)
            for f in fieldnames(T)
                new_arg = Symbol("$(arg)_$(f)")
                push!(new_args, new_arg)
                f_ex = Expr(:., arg, QuoteNode(f))
                st[f_ex] = new_arg
            end
        else
            push!(new_args, arg)
        end
    end
    return new_args, subs(ex, st), st
end

function destruct_inputs(inputs)
    new_inputs = []
    for (k, v) in inputs
        if isstruct(v)
            for f in fieldnames(v)
                new_arg = Symbol("$(k)_$(f)")
                new_val = getfield(v, f)
                push!(new_inputs, new_arg => new_val)
            end
        else
            push!(new_inputs, k => v)
        end
    end
    return new_inputs
end


text:
using SafeTestsets

@time begin
    @time @safetestset "Base Utilities" begin
        include("utilities/base_util.jl")
    end
    @time @safetestset "Math Utilities" begin
        include("utilities/math_util.jl")
    end
    @time @safetestset "Multi-qubits Hamiltonian Construction" begin
        include("utilities/multi_qubits_construction.jl")
    end
    @time @safetestset "Development tools" begin
        include("dev_tools.jl")
    end
    @time @safetestset "Displays" begin
        include("utilities/display.jl")
    end
    @time @safetestset "Interpolations" begin
        include("interpolations.jl")
    end
    @time @safetestset "Dense Hamiltonian" begin
        include("hamiltonian/dense_hamiltonian.jl")
    end
    @time @safetestset "Sparse Hamiltonian" begin
        include("hamiltonian/sparse_hamiltonian.jl")
    end
    @time @safetestset "Adiabatic Frame Hamiltonian" begin
        include("hamiltonian/adiabatic_frame_hamiltonian.jl")
    end
    @time @safetestset "Interpolation Hamiltonian" begin
        include("hamiltonian/interp_hamiltonian.jl")
    end
    @time @safetestset "Custom Hamiltonian" begin
        include("hamiltonian/custom_hamiltonian.jl")
    end
    @time @safetestset "Coupling" begin
        include("coupling_bath_interaction/coupling.jl")
    end
    @time @safetestset "Bath" begin
        include("coupling_bath_interaction/bath.jl")
    end
    @time @safetestset "Interactions" begin
        include("coupling_bath_interaction/interaction.jl")
    end
    @time @safetestset "Davies and AME" begin
        include("opensys/davies.jl")
    end
    @time @safetestset "Redfield" begin
        include("opensys/redfield.jl")
    end
    @time @safetestset "Lindblad" begin
        include("opensys/lindblad.jl")
    end
    @time @safetestset "Stochastic" begin
        include("opensys/stochastic.jl")
    end
    @time @safetestset "Annealing/Evolution" begin
        include("annealing.jl")
    end
    @time @safetestset "Projections" begin
        include("projection.jl")
    end
end


text:
using Test, Self

@test indexeswhere([1,2,3,1], 1) == [1,4]


text:
module MCPhylo

using Reexport
@reexport using Distributions

#################### Imports ####################

using DelimitedFiles
using SpecialFunctions
using Serialization
using Distributed
using Printf: @sprintf
using LinearAlgebra
using Plots
@reexport using Plots
using StatsPlots
@reexport using StatsPlots
using RecipesBase
using StatsBase
using Zygote
using FiniteDiff
using Showoff: showoff
using Markdown
using DataFrames
using Random
using CSV
using ChainRules
using ChainRulesCore
using DataStructures
using ProgressMeter
using CUDA
if has_cuda()
  using GPUArrays
else
  @warn "The Julia CUDA library is installed, but no CUDA device detected.
         Computation is performed without CUDA functionality."
end

import Base: Matrix, names, summary, iterate
import Base.Threads.@spawn
import Compose: Context, context, cm, gridstack, inch, MeasureOrNumber, mm, pt, px
import LinearAlgebra: cholesky, dot, BlasFloat
import Statistics: cor
import Distributions:
       ## Generic Types
       Continuous, ContinuousUnivariateDistribution, Distribution,
       MatrixDistribution, Multivariate, MultivariateDistribution, PDiagMat,
       PDMat, ScalMat, Truncated, Univariate, UnivariateDistribution,
       ValueSupport,
       ## ContinuousUnivariateDistribution Types
       Arcsine, Beta, BetaPrime, Biweight, Cauchy, Chi, Chisq, Cosine,
       Epanechnikov, Erlang, Exponential, FDist, Frechet, Gamma, Gumbel,
       InverseGamma, InverseGaussian, Kolmogorov, KSDist, KSOneSided, Laplace,
       Levy, Logistic, LogNormal, NoncentralBeta, NoncentralChisq,
       NoncentralF, NoncentralT, Normal, NormalCanon, Pareto, Rayleigh,
       SymTriangularDist, TDist, TriangularDist, Triweight, Uniform, VonMises,
       Weibull,
       ## DiscreteUnivariateDistribution Types
       Bernoulli, Binomial, Categorical, DiscreteUniform, Geometric,
       Hypergeometric, NegativeBinomial, NoncentralHypergeometric, Pareto,
       PoissonBinomial, Skellam,
       ## MultivariateDistribution Types
       Dirichlet, Multinomial, MvNormal, MvNormalCanon, MvTDist,
       VonMisesFisher,
       ## MatrixDistribution Types
       InverseWishart, Wishart,
       ## Methods
       cdf, dim, gradlogpdf, insupport, isprobvec, logpdf, logpdf!, maximum,
       minimum, pdf, quantile, rand, sample!, support, length
using LightGraphs: DiGraph, add_edge!, outneighbors,
       topological_sort_by_dfs, vertices
import StatsBase: autocor, autocov, countmap, counts, describe, predict,
       quantile, sample, sem, summarystats

include("distributions/pdmats2.jl")
using .PDMats2
include("Tree/Node_Type.jl") # We need this to get the Node type in

#################### Types ####################

ElementOrVector{T} = Union{T, Vector{T}}

#################### Variate Types ####################

abstract type ScalarVariate <: Real end
abstract type ArrayVariate{N} <: DenseArray{Float64, N} end
abstract type TreeVariate <: AbstractNode end

const AbstractVariate = Union{ScalarVariate, ArrayVariate, TreeVariate}
const NumericalVariate = Union{ScalarVariate, ArrayVariate}
const VectorVariate = ArrayVariate{1}
const MatrixVariate = ArrayVariate{2}

#################### Distribution Types ####################

const DistributionStruct = Union{Distribution,
                                 Array{UnivariateDistribution},
                                 Array{MultivariateDistribution}}

#################### Dependent Types ####################

mutable struct ScalarLogical <: ScalarVariate
  value::Float64
  symbol::Symbol
  monitor::Vector{Int}
  eval::Function
  sources::Vector{Symbol}
  targets::Vector{Symbol}
end

mutable struct ArrayLogical{N} <: ArrayVariate{N}
  value::Array{Float64, N}
  symbol::Symbol
  monitor::Vector{Int}
  eval::Function
  sources::Vector{Symbol}
  targets::Vector{Symbol}
end

mutable struct TreeLogical{T} <: TreeVariate where T<:GeneralNode
  value::T
  symbol::Symbol
  monitor::Vector{Int}
  eval::Function
  sources::Vector{Symbol}
  targets::Vector{Symbol}
end

mutable struct ScalarStochastic <: ScalarVariate
  value::Float64
  symbol::Symbol
  monitor::Vector{Int}
  eval::Function
  sources::Vector{Symbol}
  targets::Vector{Symbol}
  distr::UnivariateDistribution
end

mutable struct ArrayStochastic{N} <: ArrayVariate{N}
  value::Union{CuArray{Float64, N}, Array{Float64, N}}
  symbol::Symbol
  monitor::Vector{Int}
  eval::Function
  sources::Vector{Symbol}
  targets::Vector{Symbol}
  distr::DistributionStruct
end

mutable struct TreeStochastic{T} <: TreeVariate where T<: GeneralNode
    value::T
    symbol::Symbol
    monitor::Vector{Int}
    eval::Function
    sources::Vector{Symbol}
    targets::Vector{Symbol}
    distr::DistributionStruct
end

const AbstractLogical = Union{ScalarLogical, ArrayLogical}
const AbstractStochastic = Union{ScalarStochastic, ArrayStochastic}
const AbstractTreeStochastic = Union{TreeLogical, TreeStochastic}
const AbstractDependent = Union{AbstractLogical, AbstractStochastic, AbstractTreeStochastic}

#################### Sampler Types ####################

mutable struct Sampler{T}
  params::Vector{Symbol}
  eval::Function
  tune::T
  targets::Vector{Symbol}
end

abstract type SamplerTune end

struct SamplerVariate{T<:SamplerTune} <: VectorVariate
  value::Union{Vector{Float64}, Vector{S}} where S<:GeneralNode
  tune::T

  function SamplerVariate{T}(x::AbstractVector, tune::T) where T<:SamplerTune
    v = new{T}(x, tune)
    validate(v)
  end

  function SamplerVariate{T}(x::AbstractVector, pargs...; kargs...) where T<:SamplerTune
    if !isa(x[1], GeneralNode)
      value = convert(Vector{Float64}, x)
    else
      mt = typeof(x[1])
      value = convert(Vector{mt}, x)
    end
    new{T}(value, T(value, pargs...; kargs...))
  end
end

#################### Model Types ####################

struct ModelGraph
  graph::DiGraph
  keys::Vector{Symbol}
end

struct ModelState
  value::Vector
  tune::Vector{Any}
end

mutable struct Model
  nodes::Dict{Symbol, Any}
  samplers::Vector{Sampler}
  states::Vector{ModelState}
  iter::Int
  burnin::Int
  hasinputs::Bool
  hasinits::Bool
  likelihood::Float64
end

#################### Chains Type ####################

abstract type AbstractChains end

struct Chains <: AbstractChains
  value::Array{Float64, 3}
  range::StepRange{Int, Int}
  names::Vector{AbstractString}
  chains::Vector{Int}
  trees::Array{AbstractString, 3}
  moves::Array{Int, 1}
  tree_names::Vector{AbstractString}
end

struct ModelChains <: AbstractChains
  value::Array{Float64, 3}
  range::StepRange{Int, Int}
  names::Vector{AbstractString}
  chains::Vector{Int}
  model::Model
  trees::Array{AbstractString, 3}
  moves::Array{Int, 1}
  tree_names::Vector{AbstractString}
end

#################### Includes ####################

include("progress.jl")
include("utils.jl")
include("variate.jl")

include("distributions/constructors.jl")
include("distributions/distributionstruct.jl")
include("distributions/extensions.jl")
include("distributions/pdmatdistribution.jl")
include("Likelihood/SubstitutionModels.jl")
include("distributions/Phylodist.jl")
include("distributions/TreeConstraints.jl")
include("distributions/transformdistribution.jl")

include("model/dependent.jl")
include("model/dependent_tree.jl")
include("model/graph.jl")
include("model/initialization.jl")
include("model/mcmc.jl")
include("model/model.jl")
include("model/simulation.jl")

include("output/chains.jl")
include("output/chainsummary.jl")
include("output/discretediag.jl")
include("output/fileio.jl")
include("output/gelmandiag.jl")
include("output/gewekediag.jl")
include("output/heideldiag.jl")
include("output/mcse.jl")
include("output/modelchains.jl")
include("output/modelstats.jl")
include("output/rafterydiag.jl")
include("output/stats.jl")
include("output/plot.jl")
include("output/tree_plot.jl")

include("samplers/sampler.jl")

include("samplers/abc.jl")
include("samplers/amm.jl")
include("samplers/amwg.jl")
include("samplers/bhmc.jl")
include("samplers/bia.jl")
include("samplers/bmc3.jl")
include("samplers/bmg.jl")
include("samplers/dgs.jl")
include("samplers/dmh.jl")
include("samplers/empirical.jl")
include("samplers/hmc.jl")
include("samplers/mala.jl")
include("samplers/miss.jl")
include("samplers/nuts.jl")
include("samplers/pnuts.jl")
include("samplers/rwm.jl")
include("samplers/rwmc.jl")
include("samplers/slice.jl")
include("samplers/slicesimplex.jl")

include("Tree/Tree_Basics.jl")
include("Tree/Converter.jl")
include("Tree/Tree_moves.jl")
include("Tree/Tree_Distance.jl")
include("Tree/Tree_Traversal.jl")
include("Tree/Tree_Search.jl")
include("Tree/Tree_Legacy.jl")
include("Tree/Tree_Clustering.jl")
include("Tree/Tree_Ladderizing.jl")
include("Tree/Tree_Pruning.jl")
include("Tree/Tree_Consensus.jl")
include("Tree/Tree_Statistics.jl")

include("Parser/Parser.jl")
include("Parser/ParseCSV.jl")
include("Parser/ParseNexus.jl")
include("Parser/ParseNewick.jl")

include("Utils/FileIO.jl")

include("Likelihood/LikelihoodCalculator_Node.jl")
include("Likelihood/Parsimony.jl")
include("Likelihood/Prior.jl")
include("Likelihood/Rates.jl")
#################### Exports ####################

export
  AbstractChains,
  AbstractDependent,
  AbstractLogical,
  AbstractStochastic,
  AbstractVariate,
  ArrayLogical,
  ArrayStochastic,
  ArrayVariate,
  TreeLogical,
  TreeVariate,
  TreeStochastic,
  GeneralNode,
  AbstractNode,
  Node,
  Chains,
  Logical,
  MatrixVariate,
  Model,
  ModelChains,
  Sampler,
  SamplerTune,
  SamplerVariate,
  ScalarLogical,
  ScalarStochastic,
  ScalarVariate,
  Stochastic,
  VectorVariate

export
  BDiagNormal,
  Flat,
  SymUniform,
  CompoundDirichlet,
  PhyloDist,
  MultiplePhyloDist,
  exponentialBL

export
  autocor,
  changerate,
  contourplot,
  cor,
  describe,
  dic,
  discretediag,
  draw,
  gelmandiag,
  gettune,
  gewekediag,
  gradlogpdf,
  gradlogpdf!,
  graph,
  graph2dot,
  heideldiag,
  hpd,
  insupport,
  invlink,
  invlogit,
  link,
  logit,
  logpdf,
  logpdf!,
  mcmc,
  mcse,
  plot,
  predict,
  quantile,
  rafterydiag,
  rand,
  readcoda,
  relist,
  relist!,
  sample!,
  setinits!,
  setinputs!,
  setmonitor!,
  setsamplers!,
  summarystats,
  unlist,
  update!

@static if VERSION >= v"1.0.0"
  export showall
end

export
  ABC,
  AMM, AMMVariate,
  AMWG, AMWGVariate,
  BHMC, BHMCVariate,
  BMC3, BMC3Variate,
  BMG, BMGVariate,
  DiscreteVariate,
  DGS, DGSVariate,
  HMC, HMCVariate,
  BIA, BIAVariate,
  MALA, MALAVariate,
  MISS,
  NUTS, NUTSVariate,
  RWM, RWMVariate,
  Slice, SliceMultivariate, SliceUnivariate,
  SliceSimplex, SliceSimplexVariate,
  DMH, DMHVariate
  PNUTS, PNUTSVariate,
  Empirical, EmpiricalVariate

export
  make_tree_with_data,
  make_tree_with_data_cu,
  to_file, drop_samples,
  to_df,
  to_covariance,
  to_distance_matrix,
  node_distance,
  get_path,
  cut_samples,
  NNI!, NNI,
  SPR!, SPR,
  slide!, slide,
  swing!, swing,
  RF, randomize!,
  BHV_bounds,
  get_branchlength_vector,
  set_branchlength_vector!,
  post_order,
  pre_order,
  level_order,
  add_child!,
  remove_child!,
  delete_node!,
  insert_node!,
  find_lca,
  find_by_binary,
  find_by_name,
  find_root,
  create_tree_from_leaves,
  create_tree_from_leaves_cu,
  newick,
  tree_length,
  tree_height,
  path_length,
  get_sister,
  get_leaves,
  check_leafsets,
  neighbor_joining,
  upgma,
  prune_tree!, prune_tree,
  ladderize_tree!, ladderize_tree,
  majority_consensus_tree,
  discrete_gamma_rates,
  Restriction, JC, GTR, freeK,
  ASDSF, parsimony, ParseNewick

export
  cm,
  inch,
  mm,
  pt,
  px

#################### Deprecated ####################

include("deprecated.jl")

end # module


text:
# Code for Video 1.7:

[1, 2, 3]
# 3-element Array{Int64,1}:
#  1
#  2
#  3
Array{Int64, 1} == Vector{Int64} # true

[1 2 3]
# 1x3 Array{Int64,2}:
#  1  2  3
Array{Int64, 2} == Matrix{Int64} # true

ma = [1 2 3; 4 5 6]
# 2x3 Array{Int64,2}:
#  1  2  3
#  4  5  6

size(ma, 1) # 2
size(ma, 2) # 3
ndims(ma) # 2
eltype(ma) # Int64

rand(3, 3)
# 3x3 Array{Float64,2}:
#  0.768931   0.504026  0.272933
#  0.0635822  0.160217  0.0623945
#  0.433924   0.64907   0.920496

rand(10:20, 5, 5)
# 5x5 Array{Int64,2}:
#  11  10  19  16  10
#  19  14  18  20  11
#  14  18  16  18  19
#  18  16  11  11  14
#  10  10  18  11  19

eye(3)
# 3x3 Array{Float64,2}:
#  1.0  0.0  0.0
#  0.0  1.0  0.0
#  0.0  0.0  1.0

plus = [n + m for n in 1:4, m in 1:4]
# 4x4 Array{Int64,2}:
#  2  3  4  5
#  3  4  5  6
#  4  5  6  7
#  5  6  7  8

ma = [1 2 3; 4 5 6; 7 8 9]
# 3x3 Array{Int64,2}:
#  1  2  3
#  4  5  6
#  7  8  9

ma[1, 1] # 1
ma[2, 3] # 6
ma[2, :] # 1x3 Array{Int64,2}: 4  5  6
ma[:, 3] # 3-element Array{Int64,1}:
 # 3
 # 6
 # 9

ma[1:2, 2:3]
# 2x2 Array{Int64,2}:
#  2  3
#  5  6

ma[2:end, 2:end]
# 2x2 Array{Int64,2}:
#  5  6
#  8  9

ma[2, 2] = 0 # 0
# ma
# 3x3 Array{Int64,2}:
#  1  2  3
#  4  0  6
#  7  8  9

ma[2:3, 3] = -1 # -1
# ma
# 3x3 Array{Int64,2}:
#  1  2   3
#  4  0  -1
#  7  8  -1

a = [1 2; 3 4]
# 2x2 Array{Int64,2}:
#  1  2
#  3  4

b = [5 6; 7 8]
# 2x2 Array{Int64,2}:
#  5  6
#  7  8

c = [a b]
# 2x4 Array{Int64,2}:
#  1  2  5  6
#  3  4  7  8

c = hcat(a, b)
# 2x4 Array{Int64,2}:
#  1  2  5  6
#  3  4  7  8

c = [a; b]
# 4x2 Array{Int64,2}:
#  1  2
#  3  4
#  5  6
#  7  8

c = vcat(a, b)
# 4x2 Array{Int64,2}:
#  1  2
#  3  4
#  5  6
#  7  8

reshape(a, 1, 4)
# 1x4 Array{Int64,2}:
#  1  3  2  4

a
# 2x2 Array{Int64,2}:
#  1  2
#  3  4

b
# 2x2 Array{Int64,2}:
#  5  6
#  7  8

a + b
# 2x2 Array{Int64,2}:
#   6   8
#  10  12

a * b
# 2x2 Array{Int64,2}:
#  19  22
#  43  50

a .* b
# 2x2 Array{Int64,2}:
#   5  12
#  21  32

a
# 2x2 Array{Int64,2}:
#  1  2
#  3  4

a'
# 2x2 Array{Int64,2}:
#  1  3
#  2  4

inv(a)
# 2x2 Array{Float64,2}:
#  -2.0   1.0
#   1.5  -0.5

a * inv(a)
# 2x2 Array{Float64,2}:
#  1.0          0.0
#  8.88178e-16  1.0

a - 2
# 2x2 Array{Int64,2}:
#  -1  0
#   1  2

a * 5
# 2x2 Array{Int64,2}:
#   5  10
#  15  20

max_stars_repo_path: Section 1/1_7.jl

max_stars_repo_name: lytemar/Julia-for-Data-Science-Video

max_issues_repo_path: Section 1/1_7.jl

max_issues_repo_name: lytemar/Julia-for-Data-Science-Video

max_forks_repo_path: Section 1/1_7.jl

max_forks_repo_name: lytemar/Julia-for-Data-Science-Video


text:
include("splinesSerial.jl")

knots = [0,0,0,1,1,2,2,3,3,4,4,4]
p = sqrt(2)/2.0
controls = [[-1,0,1], [-p,p,p], [0,1,1], [p,p,p],[1,0,1], [p,-p,p],
         [0,-1,1], [-p,-p,p], [-1,0,1]]
c1 = larBezier(S1, [[-1,0,0,1],[-1,0,1,1]])
c2 = larBezier(S1, [[-p,p,0,p],[-p,p,p,p]])
c3 = larBezier(S1, [[0,1,0,1],[0,1,1,1]])
c4 = larBezier(S1, [[p,p,0,p],[p,p,p,p]])
c5 = larBezier(S1, [[1,0,0,1],[1,0,1,1]])
c6 = larBezier(S1, [[p,-p,0,p],[p,-p,p,p]])
c7 = larBezier(S1, [[0,-1,0,1],[0,-1,1,1]])
c8 = larBezier(S1, [[-p,-p,0,p],[-p,-p,p,p]])
c9 = larBezier(S1, [[-1,0,0,1],[-1,0,1,1]])
controls = [c1,c2,c3,c4,c5,c6,c7,c8,c9]

tnurbs = t_nurbs(S2,2,knots,controls)
point = [0.7, 0.2]
result = tnurbs(point)


text:
using PyCall

@pyimport theano.sandbox.cuda as cuda

const long = pybuiltin("long")

function cudandarray(blob::CuTensorBlob)
  # Create a Theano CudaNdArray wrapper for the GPU blob data
  if eltype(blob) == Float64
    error("Only Float32 is supported for Theano GPU code")
  end
  ptr = long(blob.ptr.p)
  shape = reverse(blob.shape)
  stride = 1
  strides = Int[]
  for sz in blob.shape
    push!(strides, stride)
    stride *= sz
  end
  return cuda.from_gpu_pointer(ptr, blob.shape, strides, blob)
end

function forward(backend::GPUBackend, state::TheanoLossLayerState, inputs::Vector{Blob})
  data_type = eltype(inputs[1])
  n = length(inputs[1])
  # map Mocha input to Theano
  ins = map(cudandarray, inputs)

  # TODO factor out common code cpu/gpu

  # apply the forward function to get the loss (it's returned as a 0-d array)
  loss = state.forward(ins...)[]

  state.loss += loss*state.layer.weight

  # accumulate statistics
  state.loss_accum *= state.n_accum
  state.loss_accum += state.loss * n
  state.loss_accum /= state.n_accum + n

  state.n_accum += n
end

function backward(backend::GPUBackend, state::TheanoLossLayerState,
                  inputs::Vector{Blob}, diffs::Vector{Blob})
  ins = map(cudandarray, inputs)

  for (g, diff) in zip(state.gradients, diffs)
    if isa(diff, CuTensorBlob)
      # TODO Get theano to output to diffs directly without copying

      copy!(diff, g(ins...))
    end
  end
end

max_stars_repo_name: benmoran/MochaTheano

max_issues_repo_name: benmoran/MochaTheano

max_forks_repo_name: benmoran/MochaTheano


text:
## Show & convert as Nucleotide ##

const _to_bitscheme = zeros(UInt8,256)
_to_bitscheme[Int('A')] = 136
_to_bitscheme[Int('G')] = 72
_to_bitscheme[Int('C')] = 40
_to_bitscheme[Int('T')] = 24
_to_bitscheme[Int('U')] = 24
_to_bitscheme[Int('R')] = 192
_to_bitscheme[Int('M')] = 160
_to_bitscheme[Int('W')] = 144
_to_bitscheme[Int('S')] = 96
_to_bitscheme[Int('K')] = 80
_to_bitscheme[Int('Y')] = 48
_to_bitscheme[Int('V')] = 224
_to_bitscheme[Int('H')] = 176
_to_bitscheme[Int('D')] = 208
_to_bitscheme[Int('B')] = 112
_to_bitscheme[Int('N')] = 240
_to_bitscheme[Int('X')] = 240
_to_bitscheme[Int('-')] = 4
_to_bitscheme[Int('.')] = 4
_to_bitscheme[Int('?')] = 2
_to_bitscheme[Int('a')] = 136
_to_bitscheme[Int('g')] = 72
_to_bitscheme[Int('c')] = 40
_to_bitscheme[Int('t')] = 24
_to_bitscheme[Int('u')] = 24
_to_bitscheme[Int('r')] = 192
_to_bitscheme[Int('m')] = 160
_to_bitscheme[Int('w')] = 144
_to_bitscheme[Int('s')] = 96
_to_bitscheme[Int('k')] = 80
_to_bitscheme[Int('y')] = 48
_to_bitscheme[Int('v')] = 224
_to_bitscheme[Int('h')] = 176
_to_bitscheme[Int('d')] = 208
_to_bitscheme[Int('b')] = 112
_to_bitscheme[Int('n')] = 240
_to_bitscheme[Int('x')] = 240

const _to_uint8 = zeros(UInt8,256)
_to_uint8[136] = 'A'
_to_uint8[72] = 'G'
_to_uint8[40] = 'C'
_to_uint8[24] = 'T'
_to_uint8[192] = 'R'
_to_uint8[160] = 'M'
_to_uint8[144] = 'W'
_to_uint8[96] = 'S'
_to_uint8[80] = 'K'
_to_uint8[48] = 'Y'
_to_uint8[224] = 'V'
_to_uint8[176] = 'H'
_to_uint8[208] = 'D'
_to_uint8[112] = 'B'
_to_uint8[240] = 'N'
_to_uint8[4] = '-'
_to_uint8[2] = '?'

# Must be immutable
immutable Nucleotide8bit
  byte::UInt8
  Nucleotide8bit{T<:Union{Integer, Char}}(x::T) = new(_to_bitscheme[x])
end

## Conversions ##
# All conversions change the bit representation

convert{T<:Union{Integer, Char}}(::Type{Nucleotide8bit},x::T) = Nucleotide8bit(x)
convert{T<:Union{Integer, Char}}(::Type{T},x::Nucleotide8bit) = convert(T,_to_uint8[x.byte])
convert(::Type{UInt8},x::Nucleotide8bit) = _to_uint8[x.byte]

nucleotide(base::Nucleotide8bit) =  _to_uint8[base.byte]

nucleotide8bit{T<:Union{Integer, Char}}(value::T)   = convert(Nucleotide8bit,value)

nucleotide8bit{T<:Union{Integer, Char}}(vec::Vector{T})  = convert(Vector{Nucleotide8bit},vec)
nucleotide8bit{T<:Union{Integer, Char}}(mat::Matrix{T})  = convert(Matrix{Nucleotide8bit},mat)

nucleotide8bit(str::ASCIIString)  = convert(Vector{Nucleotide8bit},str)

bytestring(seq::Vector{Nucleotide8bit}) = bytestring(convert(Vector{UInt8},seq))

macro nt8_str(s);  :(convert(Vector{Nucleotide8bit}, @b_str($s) )); end

convert(::Type{ASCIIString}, seq::Vector{Nucleotide8bit}) = ASCIIString(convert(Vector{UInt8},seq))
convert(::Type{Vector{Nucleotide8bit}}, str::ASCIIString) = convert(Vector{Nucleotide8bit},str.data)

## Show ##

show(io::IO,x::Nucleotide8bit) = (write(io,_to_uint8[x.byte]); nothing)

## Bits representation ##

bits(x::Nucleotide8bit) = bits(x.byte)

## Comparisons ##
# By conversions to an Integer type

for fun in [:(==),:(!=),:(<),:(<=)]
  @eval $(fun){T<:Union{Integer, Char}}(x::Nucleotide8bit,y::T) = $(fun)(convert(T,x),y)
  @eval $(fun){T<:Union{Integer, Char}}(x::T,y::Nucleotide8bit) = $(fun)(x,convert(T,y))
end

## Bitwise Operations ##
# Using the byte field
# return Integer

~(x::Nucleotide8bit)              = ~x.byte

for fun in [:((&)),:(|),:(($))]
  @eval $(fun){T<:Union{Integer, Char}}(x::Nucleotide8bit,y::T) = $(fun)(x.byte,y)
  @eval $(fun){T<:Union{Integer, Char}}(x::T,y::Nucleotide8bit) = $(fun)(x,y.byte)
  @eval $(fun)(x::Nucleotide8bit,y::Nucleotide8bit) = $(fun)(x.byte,y.byte)
end

for fun in [:(<<),:(>>),:((>>>))]
  @eval $(fun)(x::Nucleotide8bit,y::Int32) = $(fun)(x.byte,y)
end

## mmap doesn't work

mmap_array{N}(::Type{Nucleotide8bit},dims::NTuple{N,Int64},s::IO,offset::Int64) = throw("Not implemented for Nucleotide8bit")

## Faster methods using the Bit Level Scheme

## 136
isadenine(base::Nucleotide8bit) = base.byte == 136
function isadenine{T<:AbstractArray{Nucleotide8bit}}(s::T)
  res = BitArray(size(s))
  for i in 1:length(s)
    res[i] = s[i].byte == 136
  end
  res
end
## 72
isguanine(base::Nucleotide8bit) = base.byte == 72
function isguanine{T<:AbstractArray{Nucleotide8bit}}(s::T)
  res = BitArray(size(s))
  for i in 1:length(s)
    res[i] = s[i].byte == 72
  end
  res
end
## 40
iscytosine(base::Nucleotide8bit)= base.byte == 40
function iscytosine{T<:AbstractArray{Nucleotide8bit}}(s::T)
  res = BitArray(size(s))
  for i in 1:length(s)
    res[i] = s[i].byte == 40
  end
  res
end
## 24
isthymine(base::Nucleotide8bit) = base.byte == 24
function isthymine{T<:AbstractArray{Nucleotide8bit}}(s::T)
  res = BitArray(size(s))
  for i in 1:length(s)
    res[i] = s[i].byte == 24
  end
  res
end
## & 55 == 0
ispurine(base::Nucleotide8bit)  = base & 55 == 0
function ispurine{T<:AbstractArray{Nucleotide8bit}}(s::T)
  res = BitArray(size(s))
  for i in 1:length(s)
    res[i] = s[i] & 55 == 0
  end
  res
end
## & 199 == 0
ispyrimidine(base::Nucleotide8bit) = base & 199 == 0
function ispyrimidine{T<:AbstractArray{Nucleotide8bit}}(s::T)
  res = BitArray(size(s))
  for i in 1:length(s)
    res[i] = s[i] & 199 == 0
  end
  res
end
## & 8 == 8 if a is known surely
isknown(base::Nucleotide8bit) = base & 8 == 8
function isknown{T<:AbstractArray{Nucleotide8bit}}(s::T)
  res = BitArray(size(s))
  for i in 1:length(s)
    res[i] = s[i] & 8 == 8
  end
  res
end

## (a & b) < 16 if a and b are different surely
!=(x::Nucleotide8bit,y::Nucleotide8bit) = x & y < 16
## isknown(a) && a == b if a and b are the same surely
==(x::Nucleotide8bit,y::Nucleotide8bit) = (x & 8 == 8) && (x.byte == y.byte)

function percentGC(seq::Vector{Nucleotide8bit})
  len = length(seq)
  sum = 0
  for elem in seq
    if elem.byte == 40 || elem.byte == 72 || elem.byte == 96
      sum += 1
    end
  end
  sum/len
end

# Hash as Nucleotide
hash(base::Nucleotide8bit) = hash(nucleotide(base))

==(S1::Nucleotide8bit, S2::Nucleotide8bit) = ==(S1.byte,S2.byte)

# Julia 0.2 compatibility
isequal(S1::Nucleotide8bit, S2::Nucleotide8bit) = isequal(S1.byte,S2.byte)

const _convert_to_nt8 = Nucleotide8bit[UInt8('A') UInt8('C');UInt8('T') UInt8('G')]

function convert(::Type{Nucleotide2bitBase},s::Nucleotide8bit)
  _convert_to_base2[s & 7]
end
function convert(::Type{Nucleotide8bit},s::Nucleotide2bitBase)
  _convert_to_nt8[s.b1+1,s.b2+1]
end

==(S1::Nucleotide8bit, S2::Nucleotide2bit) = ==(convert(T,S1),S2)
==(S1::Nucleotide2bit, S2::Nucleotide8bit) = ==(S1,convert(T,S2))

# Julia 0.2 compatibility
isequal{T<:Union{Integer,Nucleotide2bit}}(S1::Nucleotide8bit, S2::T) = isequal(convert(T,S1),S2)
isequal{T<:Union{Integer,Nucleotide2bit}}(S1::T, S2::Nucleotide8bit) = isequal(S1,convert(T,S2))

# Hash Vector{Nucleotide8bit} as Vector{Nucleotide}
hash(seq::Vector{Nucleotide8bit}) = hash(nucleotide(seq))

==(S1::Vector{Nucleotide8bit}, S2::Vector{Nucleotide}) = ==(nucleotide(S1),S2)
==(S1::Vector{Nucleotide}, S2::Vector{Nucleotide8bit}) = ==(S1,nucleotide(S2))

# Julia 0.2 compatibility
isequal(S1::Vector{Nucleotide8bit}, S2::Vector{Nucleotide}) = isequal(nucleotide(S1),S2)
isequal(S1::Vector{Nucleotide}, S2::Vector{Nucleotide8bit}) = isequal(S1,nucleotide(S2))


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

include("ZMQCM.jl")

start_worker(parse(Int,ARGS[1]), ARGS[2])

max_stars_repo_name: caodi/julia

max_issues_repo_name: caodi/julia

max_forks_repo_name: caodi/julia


text:
export Spinorama, evaluate, power_response, listening_window, early_reflections, project, project_vituix, orbit_evaluator, normlat

struct Spinorama{T<:AngularSpectra}
    hor::T
    ver::T
    step::Int # angle step
    lw_hor # listening window
    lw_ver
    pr_floor # power response
    pr_ceiling
    pr_front_wall
    pr_side_wall
    pr_rear_wall
    pr_rear_wall_standard

    function Spinorama(hor::T, ver::T;
        lw_hor=30,
        lw_ver=10,
        pr_floor=-40..(-20),
        pr_ceiling=40..60,
        pr_front_wall=-30..30,
        pr_side_wall=[-80..(-40), 40..80],
        pr_rear_wall=[-175..(-90), 90..180],
        pr_rear_wall_standard=[180, -90, 90],
        ) where T
        lw_hor >= 0 || error("hor < 0")
        lw_ver >= 0 || error("ver < 0")

        step_hor = angular_step(hor.angles)
        step_ver = angular_step(ver.angles)
        step_hor ≈ step_ver || error("horizontal step ($(step_hor)) != vetical step ($(step_ver))")

        new{T}(hor, ver, step_hor,
            lw_hor, lw_ver,
            pr_floor, pr_ceiling, pr_front_wall, pr_side_wall, pr_rear_wall, pr_rear_wall_standard)
    end
end

function Interpolations.interpolate(spinorama::Spinorama, template::T) where {U, T <: AbstractSampleArray{U}}
    Spinorama(
        interpolate(spinorama.hor, template),
        interpolate(spinorama.ver, template);
        lw_hor=spinorama.lw_hor,
        lw_ver=spinorama.lw_ver,
        pr_floor=spinorama.pr_floor,
        pr_ceiling=spinorama.pr_ceiling,
        pr_front_wall=spinorama.pr_front_wall,
        pr_side_wall=spinorama.pr_side_wall,
        pr_rear_wall=spinorama.pr_rear_wall,
        pr_rear_wall_standard=spinorama.pr_rear_wall_standard,
    )
end

function angular_step(angles)
    # expects sorted, normalized to (-180°, 180°] (see normdeg), non-duplicate angles
    # full 360° expected for now
    idx, α = angleidx(angles, 0.0)
    α ≈ 0 || error("missing zero angle!")
    collect(angles) == sort(angles) || error("angles not sorted!")
    step = 360.0 / length(angles) # expected angular step
    for i in 2:length(angles)
        diff = angles[i] - angles[i-1]
        diff ≈ step || error("irregular angle steps! $diff != $step")
    end
    diff = angles[end] + angles[1]
    diff ≈ step || error("irregular angle steps! $diff != $step")
    return step
end

function evaluate(spinorama::Spinorama)
    # note that spinorama results involve magnitude only

    ref = abs.(parent(spinorama.hor[0.0]))
    names!(ref, [:reference])

    pr = power_response(spinorama)
    lw = listening_window(spinorama)
    er = early_reflections(spinorama)

    spdi = lw[:, :listening_window] ./ pr[:, :power_response]
    names!(spdi, [:SPDI])

    erdi = lw[:, :listening_window] ./ er[:, :early_reflections]
    names!(erdi, [:ERDI])

    erdi_horizontal = lw[:, :listening_window] ./ er[:, :early_reflections_horizontal]
    names!(erdi_horizontal, [:ERDI_horizontal])

    erdi_vertical = lw[:, :listening_window] ./ er[:, :early_reflections_vertical]
    names!(erdi_vertical, [:ERDI_vertical])

    erdi_standard = lw[:, :listening_window] ./ er[:, :early_reflections_standard]
    names!(erdi_standard, [:ERDI_standard])

    hcat(ref, pr, lw, er, spdi, erdi, erdi_horizontal, erdi_vertical, erdi_standard)
end

function power_response(spinorama::Spinorama)
    # weights according to:
    # Tylka, Joseph G., and Edgar Y. Choueiri. "On the calculation of full and partial directivity indices." 3D Audio and Applied Acoustics Laboratory Technical Report (Princeton University Press, Princeton, NJ) (2014).
    function orbit_weights(orbit)
        angles = orbit.angles
        Δθ = deg2rad(angular_step(angles))

        weights = Float64[]
        for θ in angles
            if θ ≈ 0.0 || θ ≈ 180.0
                w = (1 - cos(Δθ/2))/4
            else
                θ = deg2rad(θ)
                w = abs(cos(θ + Δθ/2) - cos(θ - Δθ/2))/8
            end
            push!(weights, w)
        end
        weights
    end

    pr = sqrt.(+([sum(abs.(parent(orbit)).^2 .* reshape(orbit_weights(orbit), 1, :), dims=2)
        for orbit in [spinorama.hor, spinorama.ver]]...))
    names!(pr, [:power_response])
    pr
end

function listening_window(spinorama::Spinorama)
    # according to ANSI-CTA-2034-A
    hor, ver, step = spinorama.lw_hor, spinorama.lw_ver, spinorama.step

    hor_angle_idxs = [angleidx(spinorama.hor, α).idx for α in -hor:step:hor]
    ver_angle_idxs = [angleidx(spinorama.ver, α).idx for α in -ver:step:ver if α != 0]
    h = parent(spinorama.hor[:, hor_angle_idxs])
    v = parent(spinorama.ver[:, ver_angle_idxs])
    lw = sqrt.(mean(Float64.(abs.(hcat(h, v)).^2), dims=2)) # RMS
    names!(lw, [:listening_window])
    lw
end

function early_reflections(spinorama::Spinorama)
    # according to ANSI-CTA-2034-A
    # Devantier, Allan. "Characterizing the Amplitude Response of Loudspeaker Systems." Audio Engineering Society Convention 113. Audio Engineering Society, 2002.
    step = spinorama.step

    get_angle_idxs(orbit, α) = [angleidx(orbit, α).idx]

    get_angle_idxs(orbit, αs::AbstractVector) = vcat([get_angle_idxs(orbit, α) for α in αs]...)

    get_angle_idxs(orbit, αs::ClosedInterval) = get_angle_idxs(orbit, [α for α in infimum(αs):step:supremum(αs)])

    get_angles(orbit, αs) = parent(orbit)[:, get_angle_idxs(orbit, αs)]

    floor_ = get_angles(spinorama.ver, spinorama.pr_floor)
    ceiling_ = get_angles(spinorama.ver, spinorama.pr_ceiling)
    front_wall_ = get_angles(spinorama.hor, spinorama.pr_front_wall)
    side_wall_ = get_angles(spinorama.hor, spinorama.pr_side_wall)
    rear_wall_ = get_angles(spinorama.hor, spinorama.pr_rear_wall)
    rear_wall_standard_ = get_angles(spinorama.hor, spinorama.pr_rear_wall_standard)

    mse_mean(selection...) = mean(Float64.(abs.(hcat(selection...)).^2), dims=2)

    function sqrt_mean(name, mse_means...)
        er =  sqrt.(mean([mse_means...]))
        names!(er, [name])
        er
    end

    floor_mean = mse_mean(floor_)
    ceiling_mean = mse_mean(ceiling_)
    front_wall_mean = mse_mean(front_wall_)
    side_wall_mean = mse_mean(side_wall_)
    rear_wall_mean = mse_mean(rear_wall_)
    rear_wall_mean_standard = mse_mean(rear_wall_standard_)

    er = hcat(
        sqrt_mean(:early_reflections, floor_mean, ceiling_mean, front_wall_mean, side_wall_mean, rear_wall_mean),
        sqrt_mean(:early_reflections_horizontal, front_wall_mean, side_wall_mean, rear_wall_mean),
        sqrt_mean(:early_reflections_vertical, floor_mean, ceiling_mean),
        sqrt_mean(:early_reflections_ceiling, ceiling_mean),
        sqrt_mean(:early_reflections_floor, floor_mean),
        sqrt_mean(:early_reflections_front, front_wall_mean),
        sqrt_mean(:early_reflections_rear, rear_wall_mean),
        sqrt_mean(:early_reflections_side, side_wall_mean),
        sqrt_mean(:early_reflections_standard, floor_mean, ceiling_mean, front_wall_mean, side_wall_mean, rear_wall_mean_standard),
        sqrt_mean(:early_reflections_horizontal_standard, front_wall_mean, side_wall_mean, rear_wall_mean_standard),
        sqrt_mean(:early_reflections_rear_standard, rear_wall_mean_standard),
    )
    er
end

# see http://www.movable-type.co.uk/scripts/latlong-vectors.html
# According to Vituix x: left-right (speaker width), y: down-top (heigth), z: front-back (negative z ahead speaker)
# x -> z, y-> x, z -> y
#
#  y    z
#  |   /
#  |  /
#  | /
#  +----- x
# I use standard coordinates:
# x: back-front (positive x ahead speaker), y: left-right, z: down-top
#
#   z
#   |
#   |
#   |
#   +----- y
#  /
# x
function checklatlon(lat, lon)
    @assert lat >= -90 && lat <= 90
    @assert lon > -180 && lon <= 180
end

function tovect(lat::T, lon::U) where {T <: Real,U <: Real}
    checklatlon(lat, lon)
    [cosd(lat) * cosd(lon), cosd(lat) * sind(lon), sind(lat)]
end

function tolatlon(x::T, y::U, z::V) where {T <: Real, U <: Real, V <: Real}
    lat = atand(z, sqrt(x^2 + y^2))
    lon = atand(y, x >= 0 ? abs(x) : x) # check for -0.0
    checklatlon(lat, lon)
    lat, lon
end

function tolatlon(vec)
    length(vec) == 3 || throw(ArgumentError("vector length(vec) = $(length) != 3!"))
    tolatlon(vec[1], vec[2], vec[3])
end

# right-hand rotation matrices
Rx(α) = [
    1 0 0;
    0 cosd(α) -sind(α);
    0 sind(α) cosd(α);
]

Ry(α) = [
    cosd(α) 0 sind(α);
    0 1 0;
    -sind(α) 0 cosd(α);
]

Rz(α) = [
    cosd(α) -sind(α) 0;
    sind(α) cosd(α) 0 ;
    0 0 1;
]

function sort_closest_points(orbit_eval, lat, lon, k=2)
    # checklatlon(lat, lon)
    dst(lat2, lon2) = norm(tovect(lat, lon)-tovect(lat2, lon2))
    dists = [dst(p.latlon...) for p in orbit_eval]
    idxs = sortperm(dists)[1:k]
    orbit_eval[idxs], dists[idxs]
end

function orbit_evaluator(spinorama::Spinorama, ::Val{:horizontal})
    [(α=α, latlon=(0, α), f=() -> parent(spinorama.hor)[:, i]) for (i, α) in enumerate(angles(spinorama.hor))], 0.0
end

function normlat(α)
    # normalized latitude, longitude for the vertical orbit
    @assert -180 < α && α <= 180
    α > 90 && return 180 - α, 180
    α < -90 && return -180 - α, 180
    return α, 0
end
function orbit_evaluator(spinorama::Spinorama, ::Val{:vertical})
    [(α=α, latlon=normlat(α), f=() -> parent(spinorama.ver)[:, i]) for (i, α) in enumerate(angles(spinorama.ver))], 0.0
end

function orbit_evaluator(spinorama::Spinorama, lat, lon)
    function latlon_to_orbit_angle(lat, lon)
        # checklatlon(lat, lon)
        x, y, z = tovect(lat, lon)
        p = [0, y, z] # projection to n = [1, 0, 0] , i.e, Y-Z plane
        p = p[2] / norm(p) # = sum(p .* [0, 1, 0]) / norm(p)
        sign(z) * acosd(p)
    end
    β = latlon_to_orbit_angle(lat, lon)
    b = abs(β)

    w =  b >= 90 ? [b-90, 180-b] : [90-b, b]

    # change normal ordering of angles
    npts = nchannels(spinorama.hor)

    hidx(i) = -90 < b < 90 ? i : (i < npts ? npts - i : npts)
    vidx(i) = β < 0 ? (i < npts ? npts - i : npts) : i

    # hidx(i) = i
    # vidx(i) = (i < npts ? npts - i : npts)

    orbit = [(α=α, latlon=tolatlon(Rx(β) * tovect(0, α)),
            f=() -> mean(
                        hcat(
                            parent(spinorama.hor)[:, hidx(i)],
                            parent(spinorama.ver)[:, vidx(i)]
                        ), weights(w)
                    )
        )
        for (i, α) in enumerate(angles(spinorama.hor))]
    orbit, β
end

function interpolation_orbit_evaluator(spinorama::Spinorama, lat, lon)
    # checklatlon(lat, lon)
    lat ≈ 0 && return orbit_evaluator(spinorama, Val(:horizontal))
    lon ≈ 0 && return orbit_evaluator(spinorama, Val(:vertical))
    return orbit_evaluator(spinorama, lat, lon)
end

function project(spinorama::Spinorama, lat, lon)
    # TODO: don't have to construct whole orbit to find its two closest points
    orbit_eval, β = interpolation_orbit_evaluator(spinorama, lat, lon)
    # @show β
    selpts, dists = sort_closest_points(orbit_eval, lat, lon, 2)
    pts = hcat([o.f() for o in selpts]...)
    # pts = selpts[1].f()
    # pts = selpts[1].f()
    # return pts

    mean(pts, weights([dists[2], dists[1]]))
    # rms_mean(pts, weights([dists[2], dists[1]]))
end

function project(spinorama::Spinorama; lat=0.0, lon=0.0, ldist=2.0, H=nothing, x=0.0, y=0.0, z=0.0, r=0.0, t=0.0, c=344.0)
    # meters, degrees, speed of sound in m/s
    x, y, z, ldist = tom.((x, y, z, ldist))
    c = tomps(c)
    driver_center = [x, y, z]
    # @show driver_center

    target = tovect(lat, lon) .* ldist
    tvec = target - driver_center
    # @show tvec, norm(tvec)

    # driver_trans = Ry(t) * Rz(-r)
    # driver_invtrans = Rz(-r) * Ry(t)
    driver_invtrans = Ry(t) * Rz(-r) # OK

    tvec_invtrans = driver_invtrans * tvec
    # @show tvec_invtrans, norm(tvec_invtrans)

    lat_, lon_ = tolatlon(tvec_invtrans)
    # @show lat_, lon_

    tdist = norm(tvec)
    dratio = ldist/tdist
    # @show dratio

    response = project(spinorama, lat_, lon_) .* dratio
    delay_dist = (tdist - ldist)
    # @show delay_dist/c
    delay!(response, delay_dist / c)
    names!(response, [:response])

    if isnothing(H)
        return response
    else
        return H .* response
    end
end

# Vituix axes: x -> z, y-> x, z -> y
function project_vituix(spinorama::Spinorama; lat=0.0, lon=0.0, ldist=2.0, H=nothing, x=0.0, y=0.0, z=0.0, r=0.0, t=0.0, c=344.0)
    project(spinorama; lat=lat, lon=lon, ldist=ldist, H=H, x=-z, y=x, z=y, r=r, t=t, c=c)
end

# for now
const DriverDefinition = Dict{Symbol, Any}

function project(drivers::DriverDefinition...; lat=0.0, lon=0.0, ldist=2.0)
    Y = [project_vituix(d[:spinorama]; lat=lat, lon=lon, ldist=ldist, d[:transform]...) for d in drivers]
    sum(hcat(Y...); dims=2)
end

function project_vituix(drivers::DriverDefinition...; lat=0.0, lon=0.0, ldist=2.0)
    Y = [project_vituix(d[:spinorama]; lat=lat, lon=lon, ldist=ldist, d[:transform]...) for d in drivers]
    sum(hcat(Y...); dims=2)
end