text:
"""
    guesswork_upper_bound(
        p::AbstractVector{T},
        ρBs::AbstractVector{<:AbstractMatrix};
        make_solver,
        c::AbstractVector = T.(1:length(p)),
        max_retries = 50,
        max_time = Inf,
        num_constraints = Inf,
        verbose::Bool = false,
        num_steps_per_SA_run::Integer = length(p)^2 * 500,
    ) where {T<:Number} -> NamedTuple

Computes an upper bound to the guesswork problem associated to the c-q state
specified by `p` and `ρBs`, as in [`guesswork`](@ref). A custom cost vector `c`
may be optionally passed. If the keyword argument `verbose` is set to true,
information is printed about each iteration of the algorithm.

The keyword argument `make_solver` is required, and must pass a *function that
creates a solver instances*. For example, instead of passing `SCSSolver()`, pass
`() -> SCSSolver()`. This is needed because the algorithm used in
`guesswork_upper_bound` solves a sequence of SDPs, not just one.

The algorithm has three termination criteria which are controlled by keyword
arguments. The algorithm stops when any of the following occur:

* `max_retries` simulated annealing attempts fail to find a violated constraint.
* `num_constraints` constraints have been added to the dual SDP
* The total runtime of the algorithm is projected to exceed `max_time` on the next iteration.

By default, `max_retries` is set to 50, while `num_constraints` and `max_time`
are set to infinity.

Lastly, the keyword argument `num_steps_per_SA_run` controls the runtime of the
simulated annealing algorithm. Increase `num_steps_per_SA_run` to search longer
for a violated constraint within a given simulated annealing run.
"""
function guesswork_upper_bound(
    p::AbstractVector{T},
    ρBs::AbstractVector{<:AbstractMatrix};
    make_solver,
    c::AbstractVector = T.(1:length(p)),
    max_retries = 50,
    max_time = Inf,
    num_constraints = Inf,
    verbose::Bool = false,
    num_steps_per_SA_run::Integer = length(p)^2 * 500,
    mutate! = rand_rev!,
    debug = false,
) where {T<:Number}

    length(p) == length(ρBs) ||
    throw(ArgumentError("Length of prior and vector of side information must match J"))
    J = length(p)
    dB = size(ρBs[1], 1)
    all(ρB -> size(ρB) == (dB, dB), ρBs) ||
    throw(ArgumentError("All side-information states must be square matrices of the same size"))

    if max_retries == Inf && max_time == Inf && num_constraints == Inf
        throw(ArgumentError("All three termination criteria (`max_retries`, `max_time`, and `num_constraints`) are infinite; algorithm would never terminate."))
    end

    if (length(c) < length(p))
        throw(ArgumentError("Need `length(c) >= length(p)`."))
    end

    constraints::Vector{Vector{Int}} = []
    convex_problem, convex_Y, convex_R = make_problem(p, ρBs, c, dB, constraints, T)
    solve_time::Float64 = 0.0
    convex_Y.value = Matrix{Complex{T}}(maximum(c)*I, dB, dB)
    upper_bound::T = 0.0

    # We optimize over the inverse of the permutations we are interested in
    # that way we don't need to take the inverse until the end.
    Rg_inv! = let ρBs_subnormalized = [p[x] * ρBs[x] for x = 1:J], c = c
        (R, g_inv) -> begin
            for x = 1:J
                R .+= ρBs_subnormalized[x] .* c[g_inv[x]]
            end
            R
        end
    end

    # Initialize some variables
    π = collect(1:J)
    scratch_π = copy(π)
    R_scratch::Matrix{Complex{T}} = zeros(Complex{T}, dB, dB)

    total_time_so_far::Float64 = 0.0

    # run the main loop
    # verbose is passed as a "value type" so logging statements can be compiled out when it is set to false.
    upper_bound = _loop(
        upper_bound,
        R_scratch,
        constraints,
        π,
        scratch_π,
        convex_R,
        convex_problem,
        convex_Y,
        Rg_inv!,
        num_steps_per_SA_run,
        total_time_so_far,
        max_retries,
        max_time,
        num_constraints,
        mutate!,
        make_solver,
        T,
        Val(verbose),
    )

    input_data = (J = J, K = J, c = c, p = p, ρBs = ρBs, dB = dB)

    return (
        optval = upper_bound,
        status = convex_problem.status,
        Y = evaluate(convex_Y),
        povm_outcomes = invperm.(constraints),
        input_data...,
    )
end

# We put the inner loop behind a function barrier
function _loop(
    upper_bound::T,
    R_scratch,
    constraints,
    π,
    scratch_π,
    convex_R,
    convex_problem,
    convex_Y,
    Rg_inv!,
    num_steps_per_SA_run,
    total_time_so_far,
    max_retries,
    max_time,
    num_constraints,
    mutate!,
    make_solver,
    ::Type{T},
    ::Val{verbose},
) where {T,verbose}

    while length(constraints) < num_constraints
        f = let Y::Matrix{Complex{T}} = evaluate(convex_Y)
            π -> begin
                R_scratch .= -1 .* Y
                eigmin(Hermitian(Rg_inv!(R_scratch, π)))::T
            end
        end

        # Choose `fval` as something `> 0` to enter the loop.
        fval::T = 1.0
        n_tries::Int = 0
        SA_time = @elapsed while fval > 0 && n_tries < max_retries
            n_tries += 1
            # initialize `π` to hold a random permutation, then use the simulated annealing algorithm to try to minimize `f`, starting from this permutation. Repeat until we find a violated constraint, or we've tried enough times.
            randperm!(π)
            fval = SA!(π, scratch_π, f, mutate!, num_steps_per_SA_run)
        end
        total_time_so_far += SA_time

        if fval < 0
            # add the permutation to the list of constraints and resolve the problem.
            new_constraint = copy(π)

            push!(constraints, new_constraint)
            update_problem!(convex_problem, convex_Y, new_constraint, convex_R)

            solve_time = @elapsed solve!(convex_problem, make_solver(); verbose = verbose)
            upper_bound::T = convex_problem.optval

            if verbose
                _log(
                    length(constraints),
                    fval,
                    upper_bound,
                    SA_time,
                    solve_time,
                    total_time_so_far,
                )
            end
        else
            # We've reached `max_retries` attempts to find a violated constraint.
            verbose && _log_retries(max_retries)
            break
        end

        # if we do another loop, it will take a round of SA plus a solve, so let's stop if
        # that will exceed the max time.
        projected_time = total_time_so_far + 1.1*(SA_time + solve_time)
        if projected_time > max_time
            if verbose
                _log_time(total_time_so_far, projected_time, max_time)
            end
            break
        end

    end
    return upper_bound
end

@noinline function _log(
    lconstraints,
    violation,
    upper_bound,
    SA_time,
    solve_time,
    total_time_so_far,
)
    @info "Adding constraint $(lconstraints) " violation upper_bound SA_time solve_time total_time_so_far
end

@noinline function _log_retries(max_retries)
    @info "Unable to find violated constraint after repeating ($max_retries) simulated annealing runs. Change `max_retries` to adjust this behavior."
end

@noinline function _log_time(total_time_so_far, projected_time, max_time)
    if total_time_so_far < max_time
        @info "Next iteration is projected to exceed maximum time $(max_time); terminating here. Change `max_time` to adjust this behavior." total_time_so_far projected_time
    else
        @info "Reached maximum time $(max_time); terminating here. Change `max_time` to adjust this behavior." total_time_so_far
    end
end

function make_problem(p, ρBs, c, dB, constraints, ::Type{T}) where {T}
    R = g -> sum(c[N(invperm(g), x)] * p[x] * ρBs[x] for x in eachindex(p, ρBs))
    Y = ComplexVariable(dB, dB)
    povm_outcomes = constraints
    Y = ComplexVariable(dB, dB)
    constraints = Convex.Constraint[R(outcome) ⪰ Y for outcome in povm_outcomes]
    push!(constraints, Y' == Y)
    objective = real(tr(Y))
    return compatible_problem(Convex.maximize, objective, constraints, T), Y, R
end

function update_problem!(problem, Y, new_constraint, R)
    push!(problem.constraints, R(new_constraint) ⪰ Y)
end

function get_dual_and_val(
    problem,
    Y,
    make_solver,
    ::Val{dB},
    ::Val{verbose},
    ::Type{T},
) where {dB,T,verbose}
    solve!(problem, make_solver(); verbose = verbose)
    return problem.optval::T, SMatrix{dB,dB,Complex{T}}(evaluate(Y))
end

"""
SA!(π, scratch_π, f, mutate! = rand_rev!)

Try to minimize `f` over permutations of 1:`length(π)` via a simulated annealing algorithm, returning the value of `f` on the final choice of `π`. Modifies `π` to hold the optimal permutation at the end of the search, and uses `scratch_π` to hold proposal permutations. The optional argument `mutate!` should be a function

    mutate!(π_to_mutate, scratch_π, f, J) -> nothing

which modifies `π_to_mutate` in place, can use `scratch_π` to hold temporary permutations, and has access to `f` and `J`, the length of `π_to_mutate`. By default, `rand_rev!` is chosen, which is very fast and requires no evaluations of `f`. Other choices of `mutate!` include  `swap_all!` and `two_opt!` (uses ~ `n^2` evaluations of `f`), or a random choice at each step, via e.g.

    mutate! = (π_to_mutate, scratch_π, f, J) -> begin
        if rand() < .01
            two_opt!(π_to_mutate, scratch_π, f, J)
        else
            rand_rev!(π_to_mutate, scratch_π, f, J)
        end
    end
"""
function SA!(π, scratch_π, f, mutate!, num_steps)
    J = length(π)
    init_temp = exp(8)
    final_temp = exp(-6.5)
    cool_rate = (final_temp / init_temp)^(1 / (num_steps - 1))
    # divide by cool_rate so when we first multiply we get init_temp
    temp = init_temp / cool_rate
    proposed_π = scratch_π

    current_cost = f(π)
    for _ = 1:num_steps
        temp *= cool_rate
        mutate!(proposed_π, π, f, J)
        proposed_cost = f(proposed_π)
        @fastmath accept = proposed_cost < current_cost ? true :
            rand() < exp((current_cost - proposed_cost) / temp)
        if accept
            current_cost = proposed_cost
            π .= proposed_π
        else
            proposed_π .= π
        end
    end
    # try all possible swaps at the end just to try to really have a optimal value.
    fval = two_opt!(π, scratch_π, f, J)
    return fval
end

"""
    rand_rev!(π, scratch_π = nothing, f = nothing, n = length(π))

Modifies `π` by randomly reversing a section of the permutation. Has the arguments `scratch_π` and `f` so it can be used in `SA!`.
"""
@inline function rand_rev!(π, scratch_π = nothing, f = nothing, n = length(π))
    i, j = rand(1:n), rand(1:n)
    if i > j
        i, j = j, i
    end
    reverse!(π, i, j)
end

"""
    swap_all!(π, scratch_π, f, n = length(π))

Modifies `π` in place to find a lower value of `f` by swapping each of the `n*(n-1)/2` pairs of entries and choosing the optimal one.
"""
function swap_all!(π, scratch_π, f, n = length(π))
    @inline swap!(π, i, j) = π[i], π[j] = π[j], π[i]
    _mutate_two!(π, scratch_π, f, swap!, n)
end

"""
    two_opt!(π, scratch_π, f, n = length(π))

Modifies `π` in place to find a lower value of `f` by reversing each possible ordered section of  `π`  by trying all `n*(n-1)/2` possible pairs of endpoints and choosing the optimal one.
"""
function two_opt!(π, scratch_π, f, n = length(π))
    _mutate_two!(π, scratch_π, f, reverse!, n)
end

@inline function _mutate_two!(π, scratch_π, f, mut!, n)
    best_fval_so_far = Inf
    best_π_so_far = scratch_π
    for (i, j) in combinations(1:n, 2)
        mut!(π, i, j)
        fval = f(π)
        if fval < best_fval_so_far
            best_fval_so_far = fval
            best_π_so_far .= π
        end
        mut!(π, i, j)
    end
    π .= best_π_so_far
    return best_fval_so_far
end


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

"""
    Arrow.Struct

An `ArrowVector` where each element is a "struct" of some kind with ordered, named fields, like a `NamedTuple{names, types}` or regular julia `struct`.
"""
struct Struct{T, S} <: ArrowVector{T}
    validity::ValidityBitmap
    data::S # Tuple of ArrowVector
    ℓ::Int
    metadata::Union{Nothing, Dict{String, String}}
end

Base.size(s::Struct) = (s.ℓ,)

@propagate_inbounds function Base.getindex(s::Struct{T}, i::Integer) where {T}
    @boundscheck checkbounds(s, i)
    NT = Base.nonmissingtype(T)
    if ArrowTypes.structtype(NT) === ArrowTypes.NAMEDTUPLE
        if NT !== T
            return s.validity[i] ? NT(ntuple(j->s.data[j][i], fieldcount(NT))) : missing
        else
            return NT(ntuple(j->s.data[j][i], fieldcount(NT)))
        end
    elseif ArrowTypes.structtype(NT) === ArrowTypes.STRUCT
        if NT !== T
            return s.validity[i] ? NT(ntuple(j->s.data[j][i], fieldcount(NT))...) : missing
        else
            return NT(ntuple(j->s.data[j][i], fieldcount(NT))...)
        end
    end
end

@propagate_inbounds function Base.setindex!(s::Struct{T}, v::T, i::Integer) where {T}
    @boundscheck checkbounds(s, i)
    if v === missing
        @inbounds s.validity[i] = false
    else
        NT = Base.nonmissingtype(T)
        N = fieldcount(NT)
        foreach(1:N) do j
            @inbounds s.data[j][i] = getfield(v, j)
        end
    end
    return v
end

struct ToStruct{T, i, A} <: AbstractVector{T}
    data::A # eltype is NamedTuple or some struct
end

ToStruct(x::A, j::Integer) where {A} = ToStruct{fieldtype(Base.nonmissingtype(eltype(A)), j), j, A}(x)

Base.IndexStyle(::Type{<:ToStruct}) = Base.IndexLinear()
Base.size(x::ToStruct) = (length(x.data),)

Base.@propagate_inbounds function Base.getindex(A::ToStruct{T, j}, i::Integer) where {T, j}
    @boundscheck checkbounds(A, i)
    @inbounds x = A.data[i]
    return x === missing ? ArrowTypes.default(T) : getfield(x, j)
end

arrowvector(::StructType, x::Struct, i, nl, fi, de, ded, meta; kw...) = x

function arrowvector(::StructType, x, i, nl, fi, de, ded, meta; kw...)
    len = length(x)
    validity = ValidityBitmap(x)
    T = Base.nonmissingtype(eltype(x))
    data = Tuple(arrowvector(ToStruct(x, j), i, nl + 1, j, de, ded, nothing; kw...) for j = 1:fieldcount(T))
    return Struct{eltype(x), typeof(data)}(validity, data, len, meta)
end

function compress(Z::Meta.CompressionType, comp, x::A) where {A <: Struct}
    len = length(x)
    nc = nullcount(x)
    validity = compress(Z, comp, x.validity)
    buffers = [validity]
    children = Compressed[]
    for y in x.data
        push!(children, compress(Z, comp, y))
    end
    return Compressed{Z, A}(x, buffers, len, nc, children)
end

function makenodesbuffers!(col::Struct{T}, fieldnodes, fieldbuffers, bufferoffset, alignment) where {T}
    len = length(col)
    nc = nullcount(col)
    push!(fieldnodes, FieldNode(len, nc))
    @debug 1 "made field node: nodeidx = $(length(fieldnodes)), col = $(typeof(col)), len = $(fieldnodes[end].length), nc = $(fieldnodes[end].null_count)"
    # validity bitmap
    blen = nc == 0 ? 0 : bitpackedbytes(len, alignment)
    push!(fieldbuffers, Buffer(bufferoffset, blen))
    @debug 1 "made field buffer: bufferidx = $(length(fieldbuffers)), offset = $(fieldbuffers[end].offset), len = $(fieldbuffers[end].length), padded = $(padding(fieldbuffers[end].length, alignment))"
    bufferoffset += blen
    for child in col.data
        bufferoffset = makenodesbuffers!(child, fieldnodes, fieldbuffers, bufferoffset, alignment)
    end
    return bufferoffset
end

function writebuffer(io, col::Struct, alignment)
    @debug 1 "writebuffer: col = $(typeof(col))"
    @debug 2 col
    writebitmap(io, col, alignment)
    # write values arrays
    for child in col.data
        writebuffer(io, child, alignment)
    end
    return
end

ext: jl

max_stars_repo_name: sweb/arrow

max_issues_repo_name: sweb/arrow

max_forks_repo_name: sweb/arrow


text:
set_error(bs::ByteStream)::Nothing = (bs.error = true; nothing)

Base.error(bs::ByteStream)::Bool = bs.error

end_input!(bs::ByteStream)::Nothing = (bs.has_eof = true; nothing)

input_ended(bs::ByteStream)::Bool = bs.has_eof

buffer_empty(bs::ByteStream)::Bool = isempty(bs.buf)

Base.eof(bs::ByteStream)::Bool = bs.has_eof && isempty(bs.buf)

function peek_out(bs::ByteStream, len)
    @inbounds @simd for i in 1:len
        bs.str_arr[i] = bs.buf[i]
    end
    return unsafe_string(pointer(bs.str_arr), len)
end

function pop_output!(bs::ByteStream, len)
    bs.buf.first = mod1(bs.buf.first + len, bs.buf.capacity)
    bs.buf.length -= len
    bs.bytes_read_count += len
    nothing
end

remain_cap(bs::ByteStream)::Int = bs.buf.capacity - length(bs.buf)

buffer_size(bs::ByteStream)::Int = length(bs.buf)

function write!(bs::ByteStream, data::AbstractString)::Int
    isempty(data) && return 0
    bytes_to_write = min(length(data), remain_cap(bs))
    @inbounds append!(bs.buf, @view(codeunits(data)[1:bytes_to_write]))
    bs.bytes_written_count += bytes_to_write
    return bytes_to_write
end

function Base.read!(bs::ByteStream, len::Int=typemax(Int))
    read_len = min(len, buffer_size(bs))
    s =  peek_out(bs, read_len)
    pop_output!(bs, read_len)
    return s
end

max_stars_repo_name: AquaIndigo/JLSponge

max_issues_repo_name: AquaIndigo/JLSponge

max_forks_repo_name: AquaIndigo/JLSponge


text:
function plot_mesh!(ax, m; d0 = m.d, vertex_score = true)
    ax.clear()
    ax.set_aspect("equal")
    ax.axis("off")
    ax.set_xlim(-1.1, 1.1)
    ax.set_ylim(-1.1, 1.1)

    ax.tripcolor(
        m.p[:, 1],
        m.p[:, 2],
        Array(m.t .- 1),
        0 * m.t[:, 1],
        cmap = "Set3",
        edgecolors = "k",
        linewidth = 1,
        facecolor = "g",
    )

    dd = m.d - d0
    tpars = Dict(
        :color => "w",
        #:fontfamily=>"sans-serif",
        :horizontalalignment => "center",
        :verticalalignment => "center",
        :fontsize => 12,
        :fontweight => "bold",
    )

    if vertex_score
        ax.scatter(m.p[dd.<0, 1], m.p[dd.<0, 2], 450, color = "r")
        ax.scatter(m.p[dd.>0, 1], m.p[dd.>0, 2], 450, color = "m")
        # ax.scatter(m.p[dd.==0, 1], m.p[dd.==0, 2], 450, color = "b")
        for i = 1:size(m.p, 1)
            if dd[i] != 0
                txt = string(dd[i])
                if dd[i] > 0
                    txt = "+" * txt
                end
                ax.text(m.p[i, 1], m.p[i, 2], txt; tpars...)
            end
        end
    end
end

function plot_mesh(m; d0 = m.d, vertex_score = true, figsize = 10, filename = "")
    fig, ax = PyPlot.subplots(figsize = (figsize, figsize))
    plot_mesh!(ax, m, d0 = d0, vertex_score = vertex_score)

    if length(filename) > 0
        fig.savefig(filename)
    end

    return fig, ax
end


text:
using Test, YaoBase, BitBasis, YaoArrayRegister, YaoBlocks, LinearAlgebra

A(i, j) = control(i, j => shift(2π / (1 << (i - j + 1))))
B(n, i) = chain(n, i == j ? put(i => H) : A(j, i) for j in i:n)
qft(n) = chain(B(n, i) for i in 1:n)

struct QFT{N,T} <: PrimitiveBlock{N,T} end

QFT(::Type{T}, n::Int) where {T} = QFT{n,T}()
QFT(n::Int) = QFT(ComplexF64, n)
circuit(::QFT{N}) where {N} = qft(N)
YaoBlocks.mat(x::QFT) = mat(circuit(x))
YaoBlocks.print_block(io::IO, x::QFT{N}) where {N} = print(io, "QFT($N)")

using FFTW, LinearAlgebra

function YaoBlocks.apply!(r::ArrayReg, x::QFT)
    α = sqrt(length(statevec(r)))
    invorder!(r)
    lmul!(α, ifft!(statevec(r)))
    return r
end

r = rand_state(5)
r1 = r |> copy |> QFT(5)
r2 = r |> copy |> circuit(QFT(5))

@test r1 ≈ r2

Hadamards(n) = repeat(H, 1:n)
ControlU(n, m, U) = chain(n + m, control(k, n+1:n+m => matblock(U^(2^(k - 1)))) for k in 1:n)
PE(n, m, U) = chain(
    n + m, # total number of the qubits
    concentrate(Hadamards(n), 1:n), # apply H in local scope
    ControlU(n, m, U),
    concentrate(QFT(n)', 1:n),
)

using LinearAlgebra

N, M = 3, 5
P = eigen(rand_unitary(1 << M)).vectors
θ = Int(0b110) / 1 << N
phases = rand(1 << M)
phases[bit"010"] = θ
U = P * Diagonal(exp.(2π * im * phases)) * P'

psi = P[:, 3]

λ = exp(2π * im * θ)
@test isapprox(U * psi, λ * psi; atol = 1e-8)

r = join(ArrayReg(psi), zero_state(N))
r |> PE(N, M, U)

p, a = findmax(probs(partial_tr(r, N+1:N+M)))
@test p ≈ 1.0
@test breflect(Int(0b110); nbits = 3) == a - 1


text:
#=
# 3D Poisson, Dirichlet bc
=#
if !@isdefined(Femshop)
    include("../Femshop.jl");
    using .Femshop
end
init_femshop("poisson3d");

# Try making an optional log
@useLog("poisson3dlog")

n = 10;
ord = 3;

# Set up the configuration (order doesn't matter)
@domain(3, SQUARE, UNIFORM_GRID)    # dimension, geometry, decomposition
@solver(CG)                         # DG, CG, etc.
@functionSpace(LEGENDRE, ord)         # function, order (or use testFunction and trialFunction)
@nodes(LOBATTO)                     # elemental node arrangement

# Specify the problem
@mesh(HEXMESH, n)                   # .msh file or generate our own

@variable(u)                        # same as @variable(u, SCALAR)

@testSymbol(v)                    # sets the symbol for a test function

@boundary(u, 1, DIRICHLET, 0)

# Write the weak form
@coefficient(f, "-3*pi*pi*sin(pi*x)*sin(pi*y)*sin(pi*z)")
@weakForm(u, "-dot(grad(u), grad(v)) - f*v")

solve(u);

# exact solution is sin(pi*x)*sin(pi*y)*sin(pi*z)
# check error
maxerr = 0;
exact(x,y,z) = sin(pi*x)*sin(pi*y)*sin(pi*z);

for i=1:size(Femshop.grid_data.allnodes,2)
    x = Femshop.grid_data.allnodes[1,i];
    y = Femshop.grid_data.allnodes[2,i];
    z = Femshop.grid_data.allnodes[3,i];
    err = abs(u.values[i] - exact(x,y,z));
    global maxerr;
    maxerr = max(err,maxerr);
end
println("max error = "*string(maxerr));

# solution is stored in the variable's "values"
# using Plots
# pyplot();
# N = n*ord+1;
# half = Int(round(N/2));
# range = (N*N*half+1):(N*N*(half+1));
# display(plot(Femshop.grid_data.allnodes[1,range], Femshop.grid_data.allnodes[2,range], u.values[range], st=:surface))

# check
log_dump_config();
log_dump_prob();

@finalize()


text:
rmps = readdir(joinpath(@__DIR__, "../dat/rmp"))

for (solver, ext) in [("cplex", :cpx), ("gurobi", :grb), ("mosek", :msk) , ("tulip", :tlp), ("tulip_uba", :tlpu)]
    for frmp in rmps
        println("julia --project --sysimage=JuliaLP.so src/RMP/$(solver).jl dat/rmp/$(frmp) > logs/rmp/$(solver)/$(frmp).$(ext) 2>&1")
    end
end

max_stars_repo_name: mtanneau/LPBenchmarks

max_issues_repo_name: mtanneau/LPBenchmarks

max_forks_repo_name: mtanneau/LPBenchmarks


text:
using MAB

alg = Algorithms.SWUCB( 3, 5 )

println( Algorithms.info_str(alg,false) )
for i = 1:10
    idx = Algorithms.getArmIndex( alg )
    r   = rand()
    Algorithms.updateReward!( alg, r )
    println( "i = $i" )
    println( "\tnoOfSteps       : ", alg.noOfSteps )
    println( "\tlastPlayedArm   : ", alg.lastPlayedArm )
    println( "\tcummReward      : ", alg.cummReward )
    println( "\tcount           : ", alg.count )
    println( "\tswCount         : ", alg.swCount )
    println( "\tarmsInWindow    : ", alg.armsInWindow )
    println( "\trewardsInWindow : ", alg.rewardsInWindow )
    println( "\tindices         : ", alg.indices)
end
Algorithms.reset!( alg )
println( Algorithms.info_str(alg,false) )
for i = 1:10
    idx = Algorithms.getArmIndex( alg )
    r   = rand()
    Algorithms.updateReward!( alg, r )
    println( "i = $i" )
    println( "\tnoOfSteps       : ", alg.noOfSteps )
    println( "\tlastPlayedArm   : ", alg.lastPlayedArm )
    println( "\tcummReward      : ", alg.cummReward )
    println( "\tcount           : ", alg.count )
    println( "\tswCount         : ", alg.swCount )
    println( "\tarmsInWindow    : ", alg.armsInWindow )
    println( "\trewardsInWindow : ", alg.rewardsInWindow )
    println( "\tindices         : ", alg.indices)
end
Algorithms.reset!( alg )


text:
# Benders Algorithm for Mixed Integer Linear Programming
# To Solve the News Boy Problem
# Edward J. Xu
# Feb 28th, 2019
using JuMP
using GLPKMathProgInterface
# Benders Decomposition algorithm
demand = [12, 14, 16, 18, 20, 22, 24, 26, 28, 30] # Demand of newspapers in each scenario
lengthS = length(demand)
prob = [0.05, 0.10, 0.10, 0.10, 0.15, 0.15, 0.10, 0.10, 0.10, 0.05 ] # probability of scenario
c = 20  # purchase price
p = 70  # selling price
h = 10  # scrap value
y = 20
let
    mas = Model(solver = GLPKSolverMIP())
    @variable(mas, q)
    @variable(mas, 0 <= y <= 40, Int)
    # add master-problem objective
    @objective(mas, Min, (c - h) * y + q)
    boundUp = Inf
    boundLow = - Inf
    epsilon = 0.01
    # definition of sub variables
    uBar = ones(10)
    # initial value of master variables
    yBar = 0
    #
    objMas = 0
    objSub = 0
    timesIteration = 1
    while (boundUp - boundLow > epsilon && timesIteration < 10)
        # sub-problem
        let
            sub = Model(solver = GLPKSolverLP())
            # create full sub-problem model, using the value of master variables yBar
            @variable(sub, u[1: 10] >= 0)
            @objective(sub, Max, - u)
            consX_1 = @constraint(sub, - (u_1 + u_2 + u_3 + u_4 + u_5 + u_6 + u_7 + u_8 + u_9 + u_10) <= prob[1] * (h - p))
            consX_2 = @constraint(sub, - (u_1 + u_2 + u_3 + u_4 + u_5 + u_6 + u_7 + u_8 + u_9 + u_10) <= prob[2] * (h - p))
            consX_3 = @constraint(sub, - (u_1 + u_2 + u_3 + u_4 + u_5 + u_6 + u_7 + u_8 + u_9 + u_10) <= prob[3] * (h - p))
            consX_4 = @constraint(sub, - (u_1 + u_2 + u_3 + u_4 + u_5 + u_6 + u_7 + u_8 + u_9 + u_10) <= prob[4] * (h - p))
            consX_5 = @constraint(sub, - (u_1 + u_2 + u_3 + u_4 + u_5 + u_6 + u_7 + u_8 + u_9 + u_10) <= prob[5] * (h - p))
            consX_6 = @constraint(sub, - (u_1 + u_2 + u_3 + u_4 + u_5 + u_6 + u_7 + u_8 + u_9 + u_10) <= prob[6] * (h - p))
            consX_7 = @constraint(sub, - (u_1 + u_2 + u_3 + u_4 + u_5 + u_6 + u_7 + u_8 + u_9 + u_10) <= prob[7] * (h - p))
            consX_8 = @constraint(sub, - (u_1 + u_2 + u_3 + u_4 + u_5 + u_6 + u_7 + u_8 + u_9 + u_10) <= prob[8] * (h - p))
            consX_9 = @constraint(sub, - (u_1 + u_2 + u_3 + u_4 + u_5 + u_6 + u_7 + u_8 + u_9 + u_10) <= prob[9] * (h - p))
            consX_10 = @constraint(sub, - (u_1 + u_2 + u_3 + u_4 + u_5 + u_6 + u_7 + u_8 + u_9 + u_10) <= prob[10] * (h - p))
            solve(sub)
            objSub = getobjectivevalue(sub)
            u1Bar = getvalue(u_1)
            u2Bar = getvalue(u_2)
            u3Bar = getvalue(u_3)
            u4Bar = getvalue(u_4)
            u5Bar = getvalue(u_5)
            u6Bar = getvalue(u_6)
            u7Bar = getvalue(u_7)
            u8Bar = getvalue(u_8)
            u9Bar = getvalue(u_9)
            u10Bar = getvalue(u_10)
        end
        boundUp = min(boundUp, objSub + (c - h) * yBar)
        # master-problem only add the feasibility constraint
        @constraint(mas, sum((- demand .- y).* [u1Bar, u2Bar, u3Bar, u4Bar, u5Bar, u6Bar, u7Bar, u8Bar, u9Bar, u10Bar]) <= q)
        solve(mas)
        objMas = getobjectivevalue(mas)
        yBar = getvalue(y)
        boundLow = max(objMas, boundLow)
        println("timesIteration: $(timesIteration)\t boundUp: $(boundUp)\t boundLow: $(boundLow)")
        timesIteration += 1
    end
end
println("Correct Ending")


text:
module MyTools

using LinearAlgebra

include("vec.jl")
export Vec, Vec2, Vec3, Vec4
include("pauli.jl")
export σ0, σx, σy, σz, σv, ⊗

module Atom
    using Printf
    include("atom.jl")
    export Element, NoneElement, is_none, ElementConstructType,
        getZ, atomic_number, symbol, chinese, english, latin, pinyin,
        find_element_with_Z,
        find_element_with_symbol,
        find_element_with_chinese,
        find_element_with_english,
        find_element_with_latin,
        find_element_with_pinyin,
        find_element,
        show_element_table
end # module AtomName

module Nucleus
    using ..Atom
    import ..Atom:getZ
    include("nucleus.jl")
    export Isotope,
        getN, getA,
        SingleParticleOrbit,
        JOrbit, MOrbit,
        name,
        NuclearShell,
        s_shell, p_shell, sd_shell, pf_shell,
        j_orbits, m_orbits,
        jsize, msize, psize, nsize,
        max_pj, max_nj,
        partition,
        m_config_size,
        HO_shell, HO_orbits,
        ValenceSpace,
        valence,
        p_space, sd_space, pf_space
end

module Ini
    include("inifile.jl")
    export IniFile, IniItem, IniSection,
       isgood, error_msg, saveAs,
       getValue, getString, getBool, getInt, getFloat
end # module Ini

module Tool
    include("tool.jl")
    export GPA, get_id_check_code, check_idcard
end # module Tool

end # module MyTools


text:
module ShapML

using Distributed
using DataFrames
using Random

include("shap_sample.jl")  # Load _shap_sample().
include("aggregate.jl")  # Load _aggregate().
include("predict.jl")  # Load _predict().

export shap

"""
    shap(explain::DataFrame,
         reference::Union{DataFrame, Nothing} = nothing,
         model,
         predict_function::Function,
         target_features::Union{Vector, Nothing} = nothing,
         sample_size::Integer = 60,
         parallel::Symbol = [:none, :samples, :features, :both],
         seed::Integer = 1,
         precision::Union{Integer, Nothing} = nothing,
         chunk::Bool = true,
         reconcile_instance::Bool = false
         )

Compute stochastic feature-level Shapley values for any ML model.

# Arguments
- `explain::DataFrame`: A DataFrame of model features with 1 or more instances to be explained using Shapley values.
- `reference`: Optional. A DataFrame with the same format as `explain` which serves as a reference group against which the Shapley value deviations from `explain` are compared (i.e., the model intercept).
- `model`: A trained ML model that is passed into `predict_function`.
- `predict_function`: A wrapper function that takes 2 required positional arguments–(1) the trained model from `model` and (2) a DataFrame of instances with the same format as `explain`. The function should return a 1-column DataFrame of model predictions; the column name does not matter.
- `target_features`: Optional. An `Array{String, 1}` of model features that is a subset of feature names in `explain` for which Shapley values will be computed. For high-dimensional models, selecting a subset of features may dramatically speed up computation time. The default behavior is to return Shapley values for all instances and features in `explain`.
- `sample_size::Integer`: The number of Monte Carlo samples used to compute the stochastic Shapley values for each feature.
- `parallel::Union{Symbol, Nothing}`: One of [:none, :samples, :features, :both]. Whether to perform the calculation serially (:none) or in parallel over Monte Carlo samples (:samples) with `pmap()` and/or multi-threaded over target features (:features) with @threads or :both.
- `seed::Integer`: A number passed to `Random.seed!()` to get reproducible results.
- `precision::Union{Integer, Nothing}`: The number of digits to `round()` results in the ouput (to reduce the size of the returned DataFrame).
- `chunk::Bool`: Default `true`. Increases speed on data with many instances and/or features. Calls the `predict()` function once per sample in `sample_size` instead of once per call to `ShapML.shap()`.
- `reconcile_instance`: EXPERIMENTAL. For each instance in `explain`, the stochastic feature-level Shapley values are adjusted so that their sum equals the model prediction. The adjustments are based on feature-level sampling variances and are typically small compared to the model prediction.

# Return
- A `size(explain, 1)` * `length(target_features)` row by 6 column DataFrame.
    + `index`: An instance in `explain`.
    + `feature_name`: Model feature.
    + `feature_value`: Feature value.
    + `shap_effect`: The average Shapley value across Monte Carlo samples.
    + `shap_effect_sd`: The standard deviation of Shapley values across Monte Carlo samples.
    + `intercept`: The average model prediction from `explain` or `reference`.
"""
function shap(;explain::DataFrame,
              reference::Union{DataFrame, Nothing} = nothing,
              model,
              predict_function::Function,
              target_features::Union{Vector, Nothing} = nothing,
              sample_size::Integer = 60,
              parallel::Union{Symbol, Nothing} = nothing,
              seed::Integer = 1,
              precision::Union{Integer, Nothing} = nothing,
              chunk::Bool = true,
              reconcile_instance::Bool = false
              )

    feature_names = String.(names(explain))
    feature_names_symbol = Symbol.(feature_names)

    if (target_features === nothing)

        target_features = copy(feature_names)  # Default is to explain with all features.

    else

        if !all(isa.(target_features, String))
            error(""""target_features" should be an array of feature names of type "String".""")
        end

        if !all(map(x -> any(x .== target_features), target_features))
            error("""One or more "target_features" is not in String.(names(explain)).""")
        end
    end
    #----------------------------------------------------------------------------
    n_instances_explain = size(explain, 1)
    n_features = size(explain, 2)
    n_target_features = length(target_features)
    #----------------------------------------------------------------------------
    if (reference === nothing)  # Default is to explain all instances in 'explain' without a specific reference group.

        reference = copy(explain)

    else

        if names(explain) != names(reference)
            error(""""explain" and "reference" should have the same model features and no outcome column.""")
        end
    end

    n_instances = size(reference, 1)
    #--------------------------------------------------------------------------
    # Parallel computation setup; the type of parallelization, if any, depends on
    # the 'parallel' argument.
    if (parallel === nothing)
        parallel = :none
    end

    if !any(parallel .== [:none, :samples, :features, :both])
         error(""""parallel" should be one of [:none, :samples, :features, :both].""")
    end
    #--------------------------------------------------------------------------
    # Create a vector of random seeds to get reproducible results with both
    # serial and parallel computations. This is not the perfect solution because there
    # could potentially be correlation between the seeds, but the effect on randomness,
    # if any, will be small. To-do: Pass in seed generator objects.
    Random.seed!(seed)
    seeds = abs.(rand(Int, sample_size))
    #--------------------------------------------------------------------------
    # Main Shapley value computation from _shap_sample(). This code is either
    # run serially or in parallel.
    if any(parallel .== [:none, :features])

        data_predict = _shap_sample(explain,
                                    reference,
                                    n_instances,
                                    n_instances_explain,
                                    n_features,
                                    n_target_features,
                                    target_features,
                                    feature_names,
                                    feature_names_symbol,
                                    sample_size,
                                    parallel,
                                    seeds,
                                    chunk,
                                    model,  # chunk = true.
                                    predict_function,  # chunk = true.
                                    precision  # chunk = true.
                                    )

    elseif any(parallel .== [:samples, :both])

        data_predict = pmap(_i -> _shap_sample(explain,
                                               reference,
                                               n_instances,
                                               n_instances_explain,
                                               n_features,
                                               n_target_features,
                                               target_features,
                                               feature_names,
                                               feature_names_symbol,
                                               sample_size,
                                               parallel,
                                               seeds[_i],
                                               chunk,
                                               model,  # chunk = true.
                                               predict_function,  # chunk = true.
                                               precision  # chunk = true.
                                               ), 1:sample_size)
    end  # End Shapley value Monte Carlo calculation.
    #--------------------------------------------------------------------------
    # Put all Frankenstein instances from all instances passed in 'explain' into
    # a single DataFrame for the user-defined predict() function.
    data_predict = vcat(data_predict...)

    if any(parallel .== [:samples, :both])
        data_predict = vcat(data_predict...)
        data_predict.sample = repeat(1:sample_size, inner = n_instances_explain * n_target_features * 2)
    end

    data_shap = _predict(reference = reference,  # input arg.
                         data_predict = data_predict,  # Calculated.
                         model = model,  # input arg.
                         predict_function = predict_function,  # input arg.
                         n_features = n_features,  # Calculated.
                         n_target_features = n_target_features,  # Calculated.
                         n_instances_explain = n_instances_explain,  # Calculated.
                         sample_size = sample_size,  # input arg.
                         precision = precision,  # input arg.
                         chunk = chunk,  # input arg.
                         reconcile_instance = reconcile_instance,  # input arg.
                         explain = explain  # input arg; needed if reconcile_instance = true.
                         )
    #--------------------------------------------------------------------------
    # Melt the input 'explain' DataFrame for merging the model features to the Shapley values.
    data_merge = DataFrames.stack(explain, feature_names_symbol)
    rename!(data_merge, Dict(:variable => "feature_name", :value => "feature_value"))
    data_merge.feature_name = String.(data_merge.feature_name)  # Coerce for merging.

    data_merge.index = repeat(1:n_instances_explain, n_features)  # The merge index for each instance.

    # Each instance in explain has one Shapley value per instance in a long DataFrame format.
    data_shap = leftjoin(data_shap, data_merge, on = [:index, :feature_name])

    # Re-order columns for easier reading.
    DataFrames.select!(data_shap, [:index, :feature_name, :feature_value, :shap_effect, :shap_effect_sd, :intercept])

    return data_shap

end  # End shap().
end  # End module.


text:
@doc raw"""
    bollingerbands(ta::TimeArray, ma=20, width=2.0)

Bollinger Bands

# Formula

```math
\begin{align*}
  \text{Up}   & = \text{SMA} + \text{width} \times \sigma(P) \\
  \text{Mean} & = \text{SMA} \\
  \text{Down} & = \text{SMA} - \text{width} \times \sigma(P)
\end{align*}
```
"""
function bollingerbands(ta::TimeArray, ma::Integer = 20, width::AbstractFloat = 2.0)
  tama   = sma(ta, ma)
  upband = tama .+ moving(nanstd, ta, ma) .* width .* sqrt((ma - 1) / ma) # take out Bessel correction, per algorithm
  dnband = tama .- moving(nanstd, ta, ma) .* width .* sqrt((ma - 1) / ma)
  bands  =  merge(upband, dnband)
  merge(bands, tama, colnames = [:up, :down, :mean])
end

@doc raw"""
    donchianchannels(ohlc::TimeArray, n = 20; h = :High, l = :Low)

Donchian Channels

# Formula

```math
\begin{align*}
  \text{Up}   & = \max(P_1^\text{High}, \dots, P_t^\text{High}) \\
  \text{Mid}  & = \frac{\text{Up} + \text{Down}}{2} \\
  \text{Down} & = \min(P_1^\text{Low}, \dots, P_t^\text{Low})
\end{align*}
```

# References

- [TradingView Wiki]
  (https://www.tradingview.com/wiki/Donchian_Channels_(DC))
"""
function donchianchannels(ohlc::TimeArray, n::Integer = 20; h = :High, l = :Low)
  up   = rename(moving(nanmax, ohlc[h], n), :up)
  down = rename(moving(nanmin, ohlc[l], n), :down)
  mid  = rename((up .+ down) ./ 2, :mid)
  merge(up, merge(mid, down))
end

@doc raw"""
    truerange(ohlc::TimeArray; h = :High, l = :Low, c = :Close)

True Range

# Formula

```math
\text{TR} = \max (P_t^\text{High}, P_{t-1}^\text{Close}) -
            \min (P_t^\text{Low},  P_{t-1}^\text{Close})
```
"""
function truerange(ohlc::TimeArray; h = :High, l = :Low, c = :Close)
  highs    = merge(ohlc[h], lag(ohlc[c]))
  lows     = merge(ohlc[l], lag(ohlc[c]))
  truehigh = TimeArray(timestamp(highs), maximum(values(highs), dims=2), [:hi], meta(highs))
  truelow  = TimeArray(timestamp(lows),  minimum(values(lows), dims=2),  [:lo], meta(lows))
  rename(truehigh .- truelow, :tr)
end

@doc raw"""
    atr(ohlc::TimeArray, n = 14; h= :High, l= :Low, c = :Close)

Average True Range

It's the exponential moving average of [`truerange`](@ref)

# Formula

```math
\text{ATR} = \text{EMA}(\text{TR}, n)
```
"""
function atr(ohlc::TimeArray, n::Integer = 14; h = :High, l = :Low, c = :Close)
  # atr was invented by Wilder, so only his ema is currently supported
  res = ema(truerange(ohlc; h = h, l = l, c = c), n, wilder = true)
  TimeArray(timestamp(res), values(res), [:atr], meta(ohlc))
end

@doc raw"""
    keltnerbands(ohlc, n = 20, w = 2; h = :High, l = :Low, c = :Close)

Keltner Channels

Linda Bradford Raschke introduced the newer version of Keltner Channels
in the 1980s. We implement the newer version.

# Formula

```math
\begin{align*}
  \text{Up}   & = \text{Mid} + w \times \text{ATR}(n) \\
  \text{Mid}  & = \text{EMA}(P^{typical}, n) \\
  \text{Down} & = \text{Mid} - w \times \text{ATR}(n)
\end{align*}
```

# References

- [StockCharts]
  (http://stockcharts.com/school/doku.php?id=chart_school:technical_indicators:keltner_channels)

- [Wikipedia]
  (https://en.wikipedia.org/wiki/Keltner_channel)
"""
function keltnerbands(ohlc::TimeArray, n::Integer = 20, w::AbstractFloat = 2.0;
                      h = :High, l = :Low, c = :Close)
  kma = rename(ema(typical(ohlc, h=h, l=l, c=c), n), :kma)
  rng = atr(ohlc, n, h=h, l=l, c=c)

  kup = rename(kma .+ (2 .* rng), :kup)
  kdn = rename(kma .- (2 .* rng), :kdn)

  merge(kup, merge(kma, kdn))
end

@doc raw"""
    chaikinvolatility(ta, n = 10, p = 10; h = High, l = :Low)

Chaikin Volatility

# Arguments

- `n` is the smooth period
- `p` is the previous period

# Formula

```math
\text{Chaikin Vola} =
  \frac{\text{EMA}(P^\text{High}_t - P^\text{Low}_t, n) - \text{EMA}(P^\text{High}_{t-p} - P^\text{Low}_{t-p}, n)}
       {\text{EMA}(P^\text{High}_{t-p} - P^\text{Low}_{t-p}, n)}
  \times 100
```

# References

- [IncredibleCharts]
  (https://www.incrediblecharts.com/indicators/chaikin_volatility.php)
"""
function chaikinvolatility(ta::TimeArray, n::Integer = 10, p::Integer = 10;
                           h = :High, l = :Low)
  rng = ema(ta[h] .- ta[l], n)
  prev = lag(rng, p)
  rename(@.((rng - prev) / prev * 100), :chaikinvolatility)
end


text:
abstract type DataSet end
abstract type ObjectClass; end
Base.print(io::IO, oc::ObjectClass) = print(io, labels(typeof(oc))[oc.class])

struct ValidationImage
    set::Type{<:DataSet}
    idx::Int
    img
    ground_truth::ObjectClass
end

struct ValData{T<:DataSet} <: AbstractVector{ValidationImage}
    set::T
end
valimgs(set::DataSet) = ValData(set)

struct TestData{T<:DataSet} <: AbstractVector{Any}
    set::T
end
testimgs(set::DataSet) = TestData(set)

struct TrainingImage
    set::Type{<:DataSet}
    idx::Int
    img
    ground_truth::ObjectClass
end

struct TrainData{T<:DataSet} <: AbstractVector{TrainingImage}
    set::T
end
trainimgs(set::DataSet) = TrainData(set)


text:
using PairwiseAlignment
using Base.Test

let
    seq = "ACGT"
    gseq = GappedSequence(seq)
    @test isa(gseq, GappedSequence{ASCIIString})
    @test length(gseq) == 4

    gseq = GappedSequence(seq, 2)
    @test length(gseq) == 0
    @test string(gseq) == ""
    @test n_chars(gseq) == 0
    @test n_gaps(gseq) == 0
    @test leading_chars(gseq) == leading_gaps(gseq) == 0
    @test trailing_chars(gseq) == trailing_gaps(gseq) == 0
    @test gapmap(gseq) == BitVector([])
    push_chars!(gseq, 1)
    @test length(gseq) == 1
    @test string(gseq) == "C"
    @test n_chars(gseq) == 1
    @test n_gaps(gseq) == 0
    @test leading_chars(gseq) == trailing_chars(gseq) == 1
    @test leading_gaps(gseq) == trailing_gaps(gseq) == 0
    @test gapmap(gseq) == BitVector([false])
    push_gaps!(gseq, 1)
    @test length(gseq) == 2
    @test string(gseq) == "C-"
    @test n_chars(gseq) == 1
    @test n_gaps(gseq) == 1
    @test leading_chars(gseq) == 1
    @test leading_gaps(gseq) == 0
    @test trailing_chars(gseq) == 0
    @test trailing_gaps(gseq) == 1
    @test gapmap(gseq) == BitVector([false, true])

    counts = [2, 2, 2, 0]
    gseq = GappedSequence(seq, counts)
    @test length(gseq) == 6
    @test string(gseq) == "AC--GT"
    @test leading_chars(gseq) == trailing_chars(gseq) == 2
    @test leading_gaps(gseq) == trailing_gaps(gseq) == 0
    @test PairwiseAlignment.counts(gseq) == counts
    @test gapmap(gseq) == BitVector([false, false, true, true, false, false])

    counts = [2, 2]
    gseq = GappedSequence(seq, counts)
    @test length(gseq) == 4
    @test string(gseq) == "AC--"
    @test n_chars(gseq) == 2
    @test n_gaps(gseq) == 2
    @test leading_chars(gseq) == 2
    @test leading_gaps(gseq) == 0
    @test trailing_chars(gseq) == 0
    @test trailing_gaps(gseq) == 2
    @test PairwiseAlignment.counts(gseq) == counts

    counts = [0, 2, 2, 2, 2, 2]
    gseq = GappedSequence(seq, counts)
    @test length(gseq) == 10
    @test string(gseq) == "--AC--GT--"
    @test n_chars(gseq) == 4
    @test n_gaps(gseq) == 6
    @test leading_chars(gseq) == trailing_chars(gseq) == 0
    @test leading_gaps(gseq) == leading_gaps(gseq) == 2
    @test PairwiseAlignment.counts(gseq) == counts

    counts = [1, 1, 1, 1]
    gseq = GappedSequence(seq, 3, counts)
    @test length(gseq) == 4
    @test string(gseq) == "G-T-"

    counts = [2, 0]
    gseq = GappedSequence(seq, counts)
    @test length(gseq) == 2
    @test string(gseq) == "AC"
    push_chars!(gseq, 1)
    @test length(gseq) == 3
    @test string(gseq) == "ACG"
    append!(gseq, [1, 2])
    @test length(gseq) == 6
    @test string(gseq) == "ACGT--"
end

let
    # a substitution matrix
    submat = SubstitutionMatrix(eye(5, 5))
    @test submat[0,0] == 1.0
    @test submat[0,1] == 0.0
    submat[0,1] = -1.0
    @test submat[0,1] == -1.0
    # models
    affine_gap = AffineGapScoreModel(submat, gap_open_penalty=2, gap_extend_penalty=1)
    @test isa(affine_gap, AffineGapScoreModel{Float64})
    cost = CostModel(submat, insertion_cost=1, deletion_cost=2)
    @test isa(cost, CostModel{Float64})
end

immutable SimpleSubstMatrix{T} <: AbstractSubstitutionMatrix{T}
    matching_score::T
    mismatching_score::T
end
Base.getindex(m::SimpleSubstMatrix, x, y) = ifelse(x == y, m.matching_score, m.mismatching_score)

rmspaces(s) = replace(s, ' ', "")

let
    # Global Alignment
    affine_gap = AffineGapScoreModel(
        SimpleSubstMatrix(0, -6),
        gap_open_penalty=5,
        gap_extend_penalty=2
    )

    a = ""
    b = "A"
    aln1 = pairalign(GlobalAlignment(), a, b, affine_gap)
    aln2 = pairalign(GlobalAlignment(), b, a, affine_gap)
    @test aln1.score == aln2.score == -7

    a = ""
    b = "AC"
    aln1 = pairalign(GlobalAlignment(), a, b, affine_gap)
    aln2 = pairalign(GlobalAlignment(), b, a, affine_gap)
    @test aln1.score == aln2.score == -9

    a = "A"
    b = "A"
    aln1 = pairalign(GlobalAlignment(), a, b, affine_gap)
    aln2 = pairalign(GlobalAlignment(), b, a, affine_gap)
    @test aln1.score == aln2.score == 0

    a = "ACGT"
    b = "ACGT"
    aln1 = pairalign(GlobalAlignment(), a, b, affine_gap)
    aln2 = pairalign(GlobalAlignment(), b, a, affine_gap)
    @test aln1.score == aln2.score == 0

    a = "ACGT"
    b = "ATGT"
    aln1 = pairalign(GlobalAlignment(), a, b, affine_gap)
    aln2 = pairalign(GlobalAlignment(), b, a, affine_gap)
    @test aln1.score == aln2.score == -6

    a = "ACGT"
    b = "TCGA"
    aln1 = pairalign(GlobalAlignment(), a, b, affine_gap)
    aln2 = pairalign(GlobalAlignment(), b, a, affine_gap)
    @test aln1.score == aln2.score == -12

    a = "ACGT"
    b = "ATTT"
    aln1 = pairalign(GlobalAlignment(), a, b, affine_gap)
    aln2 = pairalign(GlobalAlignment(), b, a, affine_gap)
    @test aln1.score == aln2.score == -12

    a = "ACGT"
    b = "ACG"
    aln1 = pairalign(GlobalAlignment(), a, b, affine_gap)
    aln2 = pairalign(GlobalAlignment(), b, a, affine_gap)
    @test aln1.score == aln2.score == -7

    a = "ACGT"
    b = "AC"
    aln1 = pairalign(GlobalAlignment(), a, b, affine_gap)
    aln2 = pairalign(GlobalAlignment(), b, a, affine_gap)
    @test aln1.score == aln2.score == -9

    a = "ACGT"
    b = "A  T" |> rmspaces
    aln1 = pairalign(GlobalAlignment(), a, b, affine_gap)
    aln2 = pairalign(GlobalAlignment(), b, a, affine_gap)
    @test aln1.score == aln2.score == -9

    a = "GC   T" |> rmspaces
    b = "GCACGT"
    aln1 = pairalign(GlobalAlignment(), a, b, affine_gap)
    aln2 = pairalign(GlobalAlignment(), b, a, affine_gap)
    @test aln1.score == aln2.score == -11
    @test string(aln1[1]) == "GC---T"
    @test string(aln1[2]) == "GCACGT"

    a = "CGTATCCAACGG     TTGTGGGA " |> rmspaces
    #     ||||| || ||      |||| ||
    b = " GTATCAAATGGACCCACTGTGTGAG" |> rmspaces
    aln1 = pairalign(GlobalAlignment(), a, b, affine_gap)
    aln2 = pairalign(GlobalAlignment(), b, a, affine_gap)
    @test aln1.score == aln2.score == -7 + -6 + -6 + -15 + -6 + -6 + -7
end

let
    # Banded Global Alignment
    affine_gap = AffineGapScoreModel(
        SimpleSubstMatrix(0, -6),
        gap_open_penalty=5,
        gap_extend_penalty=2
    )

    a = ""
    b = "A"
    aln1 = pairalign(GlobalAlignment(), a, b, affine_gap, banded=true, lower=1, upper=1)
    aln2 = pairalign(GlobalAlignment(), b, a, affine_gap, banded=true, lower=1, upper=1)
    @test aln1.score == aln2.score == -7

    a = ""
    b = "AC"
    aln1 = pairalign(GlobalAlignment(), a, b, affine_gap, banded=true, lower=2, upper=2)
    aln2 = pairalign(GlobalAlignment(), b, a, affine_gap, banded=true, lower=2, upper=2)
    @test aln1.score == aln2.score == -9

    a = "A"
    b = "A"
    aln1 = pairalign(GlobalAlignment(), a, b, affine_gap, banded=true, lower=1, upper=1)
    aln2 = pairalign(GlobalAlignment(), b, a, affine_gap, banded=true, lower=1, upper=1)
    @test aln1.score == aln2.score == 0

    a = "ACGT"
    b = "ACGT"
    aln1 = pairalign(GlobalAlignment(), a, b, affine_gap, banded=true, lower=1, upper=1)
    aln2 = pairalign(GlobalAlignment(), b, a, affine_gap, banded=true, lower=1, upper=1)
    @test aln1.score == aln2.score == 0

    a = "ACGT"
    b = "ATGT"
    aln1 = pairalign(GlobalAlignment(), a, b, affine_gap, banded=true, lower=1, upper=1)
    aln2 = pairalign(GlobalAlignment(), b, a, affine_gap, banded=true, lower=1, upper=1)
    @test aln1.score == aln2.score == -6

    a = "ACGT"
    b = "TCGA"
    aln1 = pairalign(GlobalAlignment(), a, b, affine_gap, banded=true, lower=1, upper=1)
    aln2 = pairalign(GlobalAlignment(), b, a, affine_gap, banded=true, lower=1, upper=1)
    @test aln1.score == aln2.score == -12

    a = "ACGT"
    b = "ATTT"
    aln1 = pairalign(GlobalAlignment(), a, b, affine_gap, banded=true, lower=1, upper=1)
    aln2 = pairalign(GlobalAlignment(), b, a, affine_gap, banded=true, lower=1, upper=1)
    @test aln1.score == aln2.score == -12

    a = "ACGT"
    b = "ACG"
    aln1 = pairalign(GlobalAlignment(), a, b, affine_gap, banded=true, lower=1, upper=1)
    aln2 = pairalign(GlobalAlignment(), b, a, affine_gap, banded=true, lower=1, upper=1)
    @test aln1.score == aln2.score == -7

    a = "ACGT"
    b = "AC"
    aln1 = pairalign(GlobalAlignment(), a, b, affine_gap, banded=true, lower=2, upper=2)
    aln2 = pairalign(GlobalAlignment(), b, a, affine_gap, banded=true, lower=2, upper=2)
    @test aln1.score == aln2.score == -9

    a = "ACGT"
    b = "A  T" |> rmspaces
    aln1 = pairalign(GlobalAlignment(), a, b, affine_gap, banded=true, lower=2, upper=2)
    aln2 = pairalign(GlobalAlignment(), b, a, affine_gap, banded=true, lower=2, upper=2)
    @test aln1.score == aln2.score == -9

    a = "GC   T" |> rmspaces
    b = "GCACGT"
    aln1 = pairalign(GlobalAlignment(), a, b, affine_gap, banded=true, lower=3, upper=3)
    aln2 = pairalign(GlobalAlignment(), b, a, affine_gap, banded=true, lower=3, upper=3)
    @test aln1.score == aln2.score == -11
    @test string(aln1[1]) == "GC---T"
    @test string(aln1[2]) == "GCACGT"

    a = "TATGGCCCG TAG   TTT" |> rmspaces
    b = "TA    CCGATAGGGGTTA" |> rmspaces
    aln1_1 = pairalign(GlobalAlignment(), a, b, affine_gap, banded=true, lower=1, upper=1)
    aln2_1 = pairalign(GlobalAlignment(), b, a, affine_gap, banded=true, lower=1, upper=1)
    aln1_3 = pairalign(GlobalAlignment(), a, b, affine_gap, banded=true, lower=3, upper=3)
    aln2_3 = pairalign(GlobalAlignment(), b, a, affine_gap, banded=true, lower=3, upper=3)
    aln1_5 = pairalign(GlobalAlignment(), a, b, affine_gap, banded=true, lower=5, upper=5)
    aln2_5 = pairalign(GlobalAlignment(), b, a, affine_gap, banded=true, lower=5, upper=5)
    #@show aln1_1 aln1_3 aln1_5
    @test aln1_1.score == aln2_1.score < aln1_3.score == aln1_3.score < aln1_5.score == aln2_5.score == -13 + -7 + -11 + -6
end

let
    affine_gap = AffineGapScoreModel(
        SimpleSubstMatrix(0, -6),
        gap_open_penalty=5,
        gap_extend_penalty=3
    )
    # same length
    for (a, b) in [("A", "A"), ("ACGT", "ACGT"), ("ACGT", "ACGG")]
        aln1 = pairalign(GlobalAlignment(), a, b, affine_gap, banded=false, score_only=false)
        aln2 = pairalign(GlobalAlignment(), a, b, affine_gap, banded=false, score_only=true)
        aln3 = pairalign(GlobalAlignment(), a, b, affine_gap, banded=true,  score_only=false)
        aln4 = pairalign(GlobalAlignment(), a, b, affine_gap, banded=true,  score_only=true)
        @test aln1.score == aln2.score == aln3.score == aln4.score
    end
    for (a, b) in [("A", "AT"), ("A", "ATT"), ("ACGT", "ACGTT"), ("ACGT", "ACGGT"), ("ACGT", "ACCCGT")]
        # b is longer
        aln1 = pairalign(GlobalAlignment(), a, b, affine_gap, banded=false, score_only=false)
        aln2 = pairalign(GlobalAlignment(), a, b, affine_gap, banded=false, score_only=true)
        aln3 = pairalign(GlobalAlignment(), a, b, affine_gap, banded=true,  lower=2, upper=2, score_only=false)
        aln4 = pairalign(GlobalAlignment(), a, b, affine_gap, banded=true,  lower=2, upper=2, score_only=true)
        @test aln1.score == aln2.score == aln3.score == aln4.score

        # a is longer
        a, b = b, a
        aln1 = pairalign(GlobalAlignment(), a, b, affine_gap, banded=false, score_only=false)
        aln2 = pairalign(GlobalAlignment(), a, b, affine_gap, banded=false, score_only=true)
        aln3 = pairalign(GlobalAlignment(), a, b, affine_gap, banded=true,  lower=2, upper=2, score_only=false)
        aln4 = pairalign(GlobalAlignment(), a, b, affine_gap, banded=true,  lower=2, upper=2, score_only=true)
        @test aln1.score == aln2.score == aln3.score == aln4.score
    end
end

let
    # Local Alignment
    affine_gap = AffineGapScoreModel(
        SimpleSubstMatrix(3, -1),
        gap_open_penalty=2,
        gap_extend_penalty=1
    )

    # equivalent sequences
    a = "A"
    b = "A"
    aln = pairalign(LocalAlignment(), a, b, affine_gap)
    @test aln.score == 3
    a = "ACGT"
    b = "ACGT"
    aln = pairalign(LocalAlignment(), a, b, affine_gap)
    @test aln.score == 12

    # substitution(s)
    a = "ACGT"
    b = "ATGT"
    aln = pairalign(LocalAlignment(), a, b, affine_gap)
    @test aln.score == 3 * 3 - 1 * 1
    a = "ACGT"
    b = "ATTT"
    aln = pairalign(LocalAlignment(), a, b, affine_gap)
    @test aln.score == 3 * 2 - 1 * 2

    # b is shorter than a
    a = "ACGT"
    b = "ACG"
    aln = pairalign(LocalAlignment(), a, b, affine_gap)
    @test aln.score == 9
    a = "ACGT"
    b = "AC"
    aln = pairalign(LocalAlignment(), a, b, affine_gap)
    @test aln.score == 6

    a =   "TTAGC   TAGT" |> rmspaces
    b = "GATTAGCACGTAGTA"
    aln = pairalign(LocalAlignment(), a, b, affine_gap)
    @test string(aln[1]) == "TTAGC---TAGT"
    @test string(aln[2]) == "TTAGCACGTAGT"
    @test aln.score == 3 * 5 + (-2 + -1 * 3) + 3 * 4

    # partial approximate matching
    a =     "ACGTAG"
    b = "AGATACGTAGGGAGA"
    aln = pairalign(LocalAlignment(), a, b, affine_gap)
    @test aln.score == 18
    a =     "ACG  GGG" |> rmspaces
    b = "AGATACGTAGGGAGA"
    aln = pairalign(LocalAlignment(), a, b, affine_gap)
    @test aln.score == 3 * 6 - (2 + 1 + 1)
end

let
    # Semi-global Alignment
    affinegap = AffineGapScoreModel(
        SimpleSubstMatrix(0, -6),
        gap_open_penalty=5,
        gap_extend_penalty=2
    )

    a = "A"
    b = "A"
    aln = pairalign(SemiGlobalAlignment(), a, b, affinegap)
    @test aln.score == 0
    aln = pairalign(SemiGlobalAlignment(), a, b, affinegap, score_only=true)
    @test aln.score == 0

    a =        "ACGT"
    b = "GATHAGCACGTAGTA"
    aln = pairalign(SemiGlobalAlignment(), a, b, affinegap)
    @test aln.score == 0

    a =   "TTAGC   TAGT" |> rmspaces
    b = "GATTAGCACGTAGTA"
    aln = pairalign(SemiGlobalAlignment(), a, b, affinegap)
    @test string(aln[1]) == "TTAGC---TAGT"
    @test string(aln[2]) == "TTAGCACGTAGT"
    @test aln.score == -11
end

let
    # Hamming Distance
    a = "ACGT"
    b = "ACGT"
    aln = pairalign(HammingDistance(), a, b)
    @test aln.score == 0
    a = "ACGT"
    b = "ACCT"
    aln = pairalign(HammingDistance(), a, b)
    @test aln.score == 1
    a = "ACGT"
    b = "ACG"
    @test_throws ErrorException pairalign(HammingDistance(), a, b)
end

let
    # Levenshtein Distance
    a = "ACGT"
    b = "ACGT"
    aln = pairalign(LevenshteinDistance(), a, b)
    @test aln.score == 0
    a = "ACGT"
    b = "ACCT"
    aln = pairalign(LevenshteinDistance(), a, b)
    @test aln.score == 1
    a = "ACGT"
    b = "ACGGA"
    aln = pairalign(LevenshteinDistance(), a, b)
    @test aln.score == 2
    aln = pairalign(LevenshteinDistance(), b, a)
    @test aln.score == 2
end

immutable SimpleCostMatrix{T} <: AbstractSubstitutionMatrix{T}
    mismatching_cost::T
end

Base.getindex{T}(m::SimpleCostMatrix{T}, x, y) = ifelse(x == y, 0, m.mismatching_cost)

let
    # Edit Distance
    costmodel = CostModel(
        SimpleCostMatrix(3),
        insertion_cost=4,
        deletion_cost=5
    )

    # equivalent sequences
    a = "ACGT"
    b = "ACGT"
    aln = pairalign(EditDistance(), a, b, costmodel)
    @test aln.score == 0
    aln = pairalign(EditDistance(), a, b, costmodel, distance_only=true)
    @test aln.score == 0

    # substitution
    a = "ACGT"
    b = "ACCT"
    aln = pairalign(EditDistance(), a, b, costmodel)
    @test aln.score == 3

    # indel
    a = "ACG"
    b = "ACGT"
    aln = pairalign(EditDistance(), a, b, costmodel)
    @test aln.score == 4
    aln = pairalign(EditDistance(), b, a, costmodel)
    @test aln.score == 5
end


text:
using GoogleAnalytics
using Test

printstyled("Running tests:\n", color=:blue)

include("authentication.jl")
include("data_formatting.jl")
include("api_comms.jl")


text:
# Convert between different parameterizations.

module Transform

# TODO: don't import Model; transformations should operate on
# generic ParamSets
using ..Model
using ..SensitiveFloats
import ..Log

export DataTransform, ParamBounds, ParamBox, SimplexBox,
       get_mp_transform, enforce_bounds!

################################
# Elementary functions.

"""
Unconstrain x in the unit interval to lie in R.
"""
function inv_logit{NumType <: Number}(x::NumType)
    @assert(x >= 0)
    @assert(x <= 1)
    -log(1.0 / x - 1)
end

function inv_logit{NumType <: Number}(x::Array{NumType})
    @assert(all(x .>= 0))
    @assert(all(x .<= 1))
    -log(1.0 ./ x - 1)
end

"""
Convert x in R to lie in the unit interval.
"""
function logit{NumType <: Number}(x::NumType)
    1.0 / (1.0 + exp(-x))
end

function logit{NumType <: Number}(x::Array{NumType})
    1.0 ./ (1.0 + exp(-x))
end

"""
Convert an (n - 1)-vector of real numbers to an n-vector on the simplex, where
the last entry implicitly has the untransformed value 1.
"""
function constrain_to_simplex{NumType <: Number}(x::Vector{NumType})
    if any(x .== Inf)
        z = NumType[ x_entry .== Inf ? one(NumType) : zero(NumType) for x_entry in x]
        z ./ sum(z)
        push!(z, 0)
        return(z)
    else
        z = exp(x)
        z_sum = sum(z) + 1
        z ./= z_sum
        push!(z, 1 / z_sum)
        return(z)
    end
end

"""
Convert an n-vector on the simplex to an (n - 1)-vector in R^{n -1}.  Entries
are expressed relative to the last element.
"""
function unconstrain_simplex{NumType <: Number}(z::Vector{NumType})
    n = length(z)
    NumType[ log(z[i]) - log(z[n]) for i = 1:(n - 1)]
end

################################
# The transforms for Celeste.

immutable ParamBox
    lower_bound::Float64
    upper_bound::Float64
    scale::Float64

    function ParamBox(lower_bound, upper_bound, scale)
        @assert lower_bound > -Inf # Not supported
        @assert scale > 0.0
        @assert lower_bound < upper_bound
        new(lower_bound, upper_bound, scale)
    end
end

immutable SimplexBox
    lower_bound::Float64
    scale::Float64
    n::Int

    function SimplexBox(lower_bound, scale, n)
        @assert n >= 2
        @assert 0.0 <= lower_bound < 1 / n
        new(lower_bound, scale, n)
    end
end

# The vector of transform parameters for a Symbol.
typealias ParamBounds Dict{Symbol, Union{Vector{ParamBox}, Vector{SimplexBox}}}

###############################################
# Functions for a "free transform".

function unbox_parameter{NumType <: Number}(param::NumType, param_box::ParamBox)
    lower_bound = param_box.lower_bound
    upper_bound = param_box.upper_bound
    scale = param_box.scale

    positive_constraint = (upper_bound == Inf)

    # exp and the logit functions handle infinities correctly, so
    # parameters can equal the bounds.
    @assert(lower_bound .<= param .<= upper_bound,
                    string("unbox_parameter: param outside bounds: ",
                                 "$param ($lower_bound, $upper_bound)"))

    if positive_constraint
        return log(param - lower_bound) * scale
    else
        param_bounded = (param - lower_bound) / (upper_bound - lower_bound)
        return inv_logit(param_bounded) * scale
    end
end

function box_parameter{NumType <: Number}(
            free_param::NumType, param_box::ParamBox)
    lower_bound = param_box.lower_bound
    upper_bound = param_box.upper_bound
    scale = param_box.scale

    positive_constraint = (upper_bound == Inf)
    if positive_constraint
        return (exp(free_param / scale) + lower_bound)
    else
        return(
            logit(free_param / scale) * (upper_bound - lower_bound) + lower_bound)
    end
end

"""
Convert an unconstrained (n-1)-vector to a simplicial n-vector, z, such that
  - sum(z) = 1
  - z >= simplex_box.lower_bound
See notes for a derivation and reasoning.
"""
function simplexify_parameter{NumType <: Number}(
            free_param::Vector{NumType}, simplex_box::SimplexBox)
    n = simplex_box.n
    lower_bound = simplex_box.lower_bound
    scale = simplex_box.scale

    @assert length(free_param) == (n - 1)

    # Broadcasting doesn't work with DualNumbers and Floats. :(
    # z_sim is on an unconstrained simplex.
    z_sim = constrain_to_simplex(NumType[ p / scale for p in free_param ])
    param = NumType[ (1 - n * lower_bound) * p + lower_bound for p in z_sim ]

    param
end

"""
Invert the transformation simplexify_parameter() by converting an n-vector
on a simplex to R^{n - 1}.
"""
function unsimplexify_parameter{NumType <: Number}(
            param::Vector{NumType}, simplex_box::SimplexBox)
    n = simplex_box.n
    lower_bound = simplex_box.lower_bound
    scale = simplex_box.scale

    @assert length(param) == n
    @assert all(param .>= lower_bound)
    @assert(abs(sum(param) - 1) < 1e-14, abs(sum(param) - 1))

    # z_sim is on an unconstrained simplex.
    # Broadcasting doesn't work with DualNumbers and Floats. :(
    z_sim = NumType[ (p - lower_bound) / (1 - n * lower_bound) for p in param ]
    free_param = NumType[ p * scale for p in unconstrain_simplex(z_sim) ]

    free_param
end

##################
# Derivatives

"""
Return derivatives of an unscaled transform from free parameters to a simplex.

Args:
  - z_sim: A vector in a simplex (NB: the function returns derivatives of the
           function f(unconstrained) = simplex, but it the answer is expressed
           in terms of the output of the function.)

Returns:
  - jacobian: n by (n -1) matrix of derivatives of the simplex output (in rows)
              wrt the free parameters (in columns)
  - hessian_vec: An n-length vector of the hessian of each simplex output
                 parameter with respect to the (n-1) free input parameters.
"""
function simplex_derivatives{NumType <: Number}(z_sim::Vector{NumType})
    n = length(z_sim)
    hessian_vec = Array(Array{NumType}, n)

    for i = 1:n
        hessian_vec[i] = Array(NumType, n - 1, n - 1)
        for j=1:(n - 1), k=1:(n - 1)
            if j != k
                if (j == i)
                    hessian_vec[i][j, k] = -z_sim[i] * z_sim[k] * (1 - 2 * z_sim[i])
                elseif (k == i)
                    hessian_vec[i][j, k] = -z_sim[i] * z_sim[j] * (1 - 2 * z_sim[i])
                else
                    hessian_vec[i][j, k] = 2 * z_sim[i] * z_sim[j] * z_sim[k]
                end
            else # j == k
                if i == j # All equal
                    hessian_vec[i][j, k] = z_sim[i] * (1 - z_sim[j]) * (1 - 2 * z_sim[k])
                else # j == k, but both are different from i
                    hessian_vec[i][j, k] = - z_sim[i] * z_sim[j] * (1 - 2 * z_sim[j])
                end
            end
        end
    end

    jacobian =
        NumType[ z_sim[i] * (i == j) - z_sim[i] * z_sim[j] for i=1:n, j=1:(n - 1) ]

    jacobian, hessian_vec
end

"""
Return the derivative and hessian of a simplex transform given the constrained
parameters.

Args:
  - param: The constrained parameter (NB: the derivatives are expressed
               as a function of the constrained parameterd despite being
                     the derivative of the function unconstrained -> constrained)
    - simplex_box: A box simplex constraint
"""
function box_simplex_derivatives{NumType <: Number}(
    param::Vector{NumType}, simplex_box::SimplexBox)
    lower_bound = simplex_box.lower_bound
  scale = simplex_box.scale
  n = simplex_box.n

  @assert length(param) == n

  # z_sim is on an unconstrained simplex.
  # Broadcasting doesn't work with DualNumbers and Floats. :(
  z_sim = NumType[ (p - lower_bound) / (1 - n * lower_bound) for p in param ]

    jacobian, hessian_vec = simplex_derivatives(z_sim)

    for i in 1:n
        hessian_vec[i] *= (scale ^ 2) * (1 - n * lower_bound)
    end
    jacobian *= scale * (1 - n * lower_bound)
  jacobian, hessian_vec
end

"""
Return the derivative and hessian of a box transform given the constrained
parameters.

Args:
  - param: The constrained parameter (NB: the derivatives are expressed
               as a function of the constrained parameterd despite being
                     the derivative of the function unconstrained -> constrained)
    - param_box: A box constraint
"""
function box_derivatives{NumType <: Number}(param::NumType, param_box::ParamBox)
    lower_bound = param_box.lower_bound
  upper_bound = param_box.upper_bound
  scale = param_box.scale

    if upper_bound == Inf
        centered_param = param - lower_bound
        return scale * centered_param, scale ^ 2 * centered_param
    else
        param_range = upper_bound - lower_bound
        centered_param = (param - lower_bound) / param_range
        derivative = param_range * centered_param * (1 - centered_param) / scale
        return derivative, derivative * (1 - 2 * centered_param) / scale
    end
end

"""
A datatype containing derivatives of the transform from free to constrained
parameters.

Members:
  dparam_dfree: The Jacobian of the transformation
                contrainted_param = f(free param)
  d2param_dfree2: A vector of hessians.  Each element is the Hessian of one
                  component of the aforementioned f()
"""
type TransformDerivatives{NumType <: Number}
    dparam_dfree::Matrix{NumType}
    d2param_dfree2::Array{Matrix{NumType}}
    Sa::Int

    # TODO: use sparse matrices?
    function TransformDerivatives(Sa::Int)
        dparam_dfree =
            zeros(NumType,
                        Sa * length(CanonicalParams), Sa * length(UnconstrainedParams))
        d2param_dfree2 = Array(Matrix{NumType}, Sa * length(CanonicalParams))
        for i in 1:(Sa * length(CanonicalParams))
            d2param_dfree2[i] =
                zeros(NumType,
                            Sa * length(UnconstrainedParams), Sa * length(UnconstrainedParams))
        end
        new(dparam_dfree, d2param_dfree2, Sa)
    end
end

"""
Populate a TransformDerivatives object in place.

Args:
  - bounds: A vector containing one ParamBounds for each active source in the
            same order as active_sources.
  - transform_derivatives: TransformDerivatives to be populated.

Returns:
  Update transform_derivatives in place.
"""
function get_transform_derivatives!{NumType <: Number}(
                            vp::VariationalParams{NumType},
                            active_sources::Vector{Int},
                            bounds::Vector{ParamBounds},
                            transform_derivatives::TransformDerivatives)

  @assert transform_derivatives.Sa == length(active_sources)

  for param in fieldnames(ids), sa = 1:length(active_sources)
    s = active_sources[sa]
      constraint_vec = bounds[sa][param]

      if isa(constraint_vec[1], ParamBox) # It is a box constraint
          @assert(length(constraint_vec) == length(getfield(ids_free, param)) ==
        length(getfield(ids, param)))

          # Get each components' derivatives one by one.
          for ind = 1:length(constraint_vec)
              @assert isa(constraint_vec[ind], ParamBox)
              vp_ind = getfield(ids, param)[ind]
              vp_free_ind = getfield(ids_free, param)[ind]

              jac, hess = box_derivatives(vp[s][vp_ind], constraint_vec[ind])

              vp_sf_ind = length(CanonicalParams) * (sa - 1) + vp_ind
              vp_free_sf_ind = length(UnconstrainedParams) * (sa - 1) + vp_free_ind

              transform_derivatives.dparam_dfree[vp_sf_ind, vp_free_sf_ind] = jac
              transform_derivatives.d2param_dfree2[
                  vp_sf_ind][vp_free_sf_ind, vp_free_sf_ind] = hess
          end
      else # It is a simplex constraint

          # If a param is not a box constraint, it must have all simplex constraints.
          @assert all([ isa(constraint, SimplexBox)  for constraint in constraint_vec])

          param_size = size(getfield(ids, param))
          @assert length(constraint_vec) == param_size[2]
          for col=1:(param_size[2])
              vp_free_ind = getfield(ids_free, param)[:, col]
              vp_ind = getfield(ids, param)[:, col]
              vp_sf_ind = length(CanonicalParams) * (sa - 1) + vp_ind
              vp_free_sf_ind = length(UnconstrainedParams) * (sa - 1) + vp_free_ind

              jac, hess = Transform.box_simplex_derivatives(
                  vp[s][vp_ind], constraint_vec[col])

              transform_derivatives.dparam_dfree[vp_sf_ind, vp_free_sf_ind] = jac
              for row in 1:(param_size[1])
                  transform_derivatives.d2param_dfree2[
                      vp_sf_ind[row]][vp_free_sf_ind, vp_free_sf_ind] = hess[row]
              end
          end
      end
  end
end

function get_transform_derivatives{NumType <: Number}(
                                vp::VariationalParams{NumType},
                                active_sources::Vector{Int},
                                bounds::Vector{ParamBounds})
  transform_derivatives =
    TransformDerivatives{Float64}(length(active_sources))
  get_transform_derivatives!(vp, active_sources, bounds, transform_derivatives)
  transform_derivatives
end

######################
# Functions to take actual parameter vectors.

"""
Convert between variational parameter vectors and unconstrained parameters.

Args:
  - vp: The vector of contrainted variationl parameters
  - vp_free: The vector of uncontrainted parameters for optimization
  - bounds: ParamBounds describing the transformation
  - to_unconstrained: If true, converts vp to vp_free.  If false, converts
                      vp_free to vp.

Returns:
  If to_unconstrained is true, updates vp_free in place.
  If to_unconstrained is false, updates vp in place.
"""
function perform_transform!{NumType <: Number}(
            vp::Vector{NumType}, vp_free::Vector{NumType}, bounds::ParamBounds,
            to_unconstrained::Bool)
    for (param, constraint_vec) in bounds
        is_box = isa(bounds[param], Array{ParamBox})
        if is_box
            # Apply a box constraint to each parameter.
            @assert(length(getfield(ids, param)) == length(getfield(ids_free, param)) ==
                length(bounds[param]))
            for ind in 1:length(getfield(ids, param))
                constraint = constraint_vec[ind]
                free_ind = getfield(ids_free, param)[ind]
                vp_ind = getfield(ids, param)[ind]
                to_unconstrained ?
                    vp_free[free_ind] = unbox_parameter(vp[vp_ind], constraint):
                    vp[vp_ind] = box_parameter(vp_free[free_ind], constraint)
            end
        else
            # Apply a simplex constraint to each parameter.
            @assert isa(bounds[param], Array{SimplexBox})
            param_size = size(getfield(ids, param))

            # Each column is a simplex and should have its own
            # simplicial constraints.
            @assert length(bounds[param]) == param_size[2]
            for col in 1:(param_size[2])
                free_ind = getfield(ids_free, param)[:, col]
                vp_ind = getfield(ids, param)[:, col]
                constraint = constraint_vec[col]
                to_unconstrained ?
                    vp_free[free_ind] = unsimplexify_parameter(vp[vp_ind], constraint):
                    vp[vp_ind] = simplexify_parameter(vp_free[free_ind], constraint)
            end
        end
    end
end

function free_to_vp!{NumType <: Number}(
        vp_free::Vector{NumType}, vp::Vector{NumType}, bounds::ParamBounds)
    perform_transform!(vp, vp_free, bounds, false)
end

function vp_to_free!{NumType <: Number}(
        vp::Vector{NumType}, vp_free::Vector{NumType}, bounds::ParamBounds)
    perform_transform!(vp, vp_free, bounds, true)
end

#####################
# Conversion to and from variational parameter vectors and arrays.

"""
Transform VariationalParams to an array.

Args:
  - vp = variational parameters
  - omitted_ids = ids in ParamIndex

There is probably no use for this function, since you'll only be passing
trasformations to the optimizer, but I'll include it for completeness.
"""
function free_vp_to_array{NumType <: Number}(vp::FreeVariationalParams{NumType},
                                             omitted_ids::Vector{Int})

    left_ids = setdiff(1:length(UnconstrainedParams), omitted_ids)
    new_P = length(left_ids)
    S = length(vp)
    x_new = zeros(NumType, new_P, S)

    for p1 in 1:length(left_ids), s=1:S
        p0 = left_ids[p1]
        x_new[p1, s] = vp[s][p0]
    end

    x_new
end

"""
Transform a parameter vector to variational parameters in place.

Args:
 - xs: A (param x sources) matrix created from free variational parameters.
 - vp_free: Free variational parameters.  Only the ids not in omitted_ids
            will be updated.
 - omitted_ids: Ids to omit (from ids_free)

Returns:
 - Update vp_free in place.
"""
function array_to_free_vp!{NumType <: Number}(
    xs::Matrix{NumType}, vp_free::FreeVariationalParams{NumType},
    omitted_ids::Vector{Int})

    left_ids = setdiff(1:length(UnconstrainedParams), omitted_ids)
    P = length(left_ids)
    S = length(vp_free)
    @assert size(xs) == (P, S)

    for s in 1:S, p1 in 1:P
        p0 = left_ids[p1]
        vp_free[s][p0] = xs[p1, s]
    end
end

#########################
# Define the exported variables.

"""
Functions to move between a single source's variational parameters and a
transformation of the data for optimization.

to_vp: A function that takes transformed parameters and returns
             variational parameters
from_vp: A function that takes variational parameters and returned
                 transformed parameters
to_vp!: A function that takes (transformed paramters, variational parameters)
                and updates the variational parameters in place
from_vp!: A function that takes (variational paramters, transformed parameters)
                    and updates the transformed parameters in place
...
transform_sensitive_float: A function that takes (sensitive float, model
    parameters) where the sensitive float contains partial derivatives with
    respect to the variational parameters and returns a sensitive float with total
    derivatives with respect to the transformed parameters.
bounds: The bounds for each parameter and each object in ElboArgs.
active_sources: The sources that are being optimized.    Only these sources'
    parameters are transformed into the parameter vector.
"""
type DataTransform
    to_vp::Function
    from_vp::Function
    to_vp!::Function
    from_vp!::Function
    vp_to_array::Function
    array_to_vp!::Function
    transform_sensitive_float::Function
    bounds::Vector{ParamBounds}
    active_sources::Vector{Int}
    active_S::Int
    S::Int
end

# TODO: Maybe this should be initialized with ElboArgs with optional
# custom bounds.    Or maybe it should be part of ElboArgs with one transform
# per celestial object rather than a single object containing an array of
# transforms.
function DataTransform(bounds::Vector{ParamBounds};
                       active_sources=collect(1:length(bounds)),
                       S=length(bounds))
    @assert length(bounds) == length(active_sources)
    @assert maximum(active_sources) <= S
    active_S = length(active_sources)

    function from_vp!{NumType <: Number}(
            vp::VariationalParams{NumType}, vp_free::VariationalParams{NumType})
        S = length(vp)
        @assert length(vp_free) == active_S
        for si=1:active_S
            s = active_sources[si]
            vp_to_free!(vp[s], vp_free[si], bounds[si])
        end
    end

    function from_vp{NumType <: Number}(vp::VariationalParams{NumType})
            vp_free = Array{NumType, 1}[
                zeros(NumType, length(ids_free)) for si = 1:active_S]
            from_vp!(vp, vp_free)
            vp_free
    end

    function to_vp!{NumType <: Number}(
            vp_free::FreeVariationalParams{NumType}, vp::VariationalParams{NumType})
        @assert length(vp_free) == active_S
        for si=1:active_S
            s = active_sources[si]
            free_to_vp!(vp_free[si], vp[s], bounds[si])
        end
    end

    function to_vp{NumType <: Number}(vp_free::FreeVariationalParams{NumType})
            @assert(active_S == S,
                    string("to_vp is not supported when active_sources is a ",
                           "strict subset of all sources."))
            vp = Array{NumType, 1}[
                zeros(NumType, length(CanonicalParams)) for s = 1:S]
            to_vp!(vp_free, vp)
            vp
    end

    function vp_to_array{NumType <: Number}(vp::VariationalParams{NumType},
                                            omitted_ids::Vector{Int})
                                            vp_trans = from_vp(vp)
        free_vp_to_array(vp_trans, omitted_ids)
    end

    function array_to_vp!{NumType <: Number}(xs::Matrix{NumType},
                                             vp::VariationalParams{NumType},
                                             omitted_ids::Vector{Int})
        # This needs to update vp in place so that variables in omitted_ids
        # stay at their original values.
        vp_trans = from_vp(vp)
        array_to_free_vp!(xs, vp_trans, omitted_ids)
        to_vp!(vp_trans, vp)
    end

    # Given a sensitive float with derivatives with respect to all the
    # constrained parameters, calculate derivatives with respect to
    # the unconstrained parameters.
    #
    # Note that all the other functions in DeterministicVI calculated derivatives with
    # respect to the constrained parameterization.
    function transform_sensitive_float{NumType <: Number}(
                                    sf::SensitiveFloat,
                                    vp::VariationalParams{NumType},
                                    active_sources::Vector{Int})
        @assert size(sf.d) == (length(CanonicalParams), length(active_sources))
        @assert length(active_sources) == active_S

        transform_derivatives = get_transform_derivatives(vp, active_sources, bounds)
        sf_free = zero_sensitive_float(UnconstrainedParams, NumType, active_S)

        sf_d_vec = sf.d[:]
        sf_free.v[1] = sf.v[1]
        sf_free.d =
        reshape(transform_derivatives.dparam_dfree' * sf_d_vec,
                        length(UnconstrainedParams), active_S)

        sf_free.h =
                transform_derivatives.dparam_dfree' *
        sf.h * transform_derivatives.dparam_dfree
        for ind in 1:length(sf_d_vec)
                sf_free.h += transform_derivatives.d2param_dfree2[ind] * sf_d_vec[ind]
        end
        # Ensure exact symmetry, which is necessary for some numerical linear
        # algebra routines.
        sf_free.h = 0.5 * (sf_free.h' + sf_free.h)
        sf_free
    end

    DataTransform(to_vp, from_vp, to_vp!, from_vp!, vp_to_array, array_to_vp!,
                  transform_sensitive_float, bounds, active_sources, active_S, S)
end

function get_mp_transform{NumType <: Number}(
                          vp::VariationalParams{NumType},
                          active_sources::Vector{Int};
                          loc_width::Float64=1.5e-3)
    bounds = Array(ParamBounds, length(active_sources))

    # Note that, for numerical reasons, the bounds must be on the scale
    # of reasonably meaningful changes.
    for si in 1:length(active_sources)
        s = active_sources[si]
        bounds[si] = ParamBounds()
        bounds[si][:u] = Array(ParamBox, 2)
        u = vp[s][ids.u]
        for axis in 1:2
            bounds[si][:u][axis] =
                ParamBox(u[axis] - loc_width, u[axis] + loc_width, 1.0)
        end
        bounds[si][:r1] = Array(ParamBox, Ia)
        bounds[si][:r2] = Array(ParamBox, Ia)
        for i in 1:Ia
            bounds[si][:r1][i] = ParamBox(-1.0, 10., 1.0)
            bounds[si][:r2][i] = ParamBox(1e-4, 0.1, 1.0)
        end
        bounds[si][:c1] = Array(ParamBox, 4 * Ia)
        bounds[si][:c2] = Array(ParamBox, 4 * Ia)
        for ind in 1:length(ids.c1)
            bounds[si][:c1][ind] = ParamBox(-10., 10., 1.0)
            bounds[si][:c2][ind] = ParamBox(1e-4, 1., 1.0)
        end
        bounds[si][:e_dev] = ParamBox[ ParamBox(1e-2, 1 - 1e-2, 1.0) ]
        bounds[si][:e_axis] = ParamBox[ ParamBox(1e-2, 1 - 1e-2, 1.0) ]
        bounds[si][:e_angle] = ParamBox[ ParamBox(-10.0, 10.0, 1.0) ]
        bounds[si][:e_scale] = ParamBox[ ParamBox(0.1, 70., 1.0) ]

        const simplex_min = 0.005
        bounds[si][:a] = Array(SimplexBox, 1)
        bounds[si][:a][1] = SimplexBox(simplex_min, 1.0, 2)

        bounds[si][:k] = Array(SimplexBox, D)
        for d in 1:D
            bounds[si][:k][d] = SimplexBox(simplex_min, 1.0, 2)
        end
    end

    DataTransform(bounds, active_sources=active_sources, S=length(vp))
end

"""
Put the variational parameters within the bounds of the transform.
Returns:
  Updates vp in place.
"""
function enforce_bounds!{NumType <: Number}(
                        vp::VariationalParams{NumType},
                        active_sources::Vector{Int},
                        transform::DataTransform)
    for sa=1:transform.active_S, (param, constraint_vec) in transform.bounds[sa]
        s = active_sources[sa]
        is_box = isa(constraint_vec, Array{ParamBox})
        if is_box
            # Box parameters.
            for ind in 1:length(getfield(ids, param))
                constraint = constraint_vec[ind]
                if !(constraint.lower_bound <=
                         vp[s][getfield(ids, param)[ind]] <=
                         constraint.upper_bound)
                    # Don't set the value to exactly the lower bound to avoid Inf
                    diff = constraint.upper_bound - constraint.lower_bound
                    epsilon = diff == Inf ? 1e-12: diff * 1e-12
                    vp[s][getfield(ids, param)[ind]] =
                        min(vp[s][getfield(ids, param)[ind]], constraint.upper_bound - epsilon)
                    vp[s][getfield(ids, param)[ind]] =
                        max(vp[s][getfield(ids, param)[ind]], constraint.lower_bound + epsilon)
                end
            end
        else
            # Simplex parameters.
            param_size = size(getfield(ids, param))
            for col in 1:param_size[2]
                constraint = constraint_vec[col]
                param_sum = zero(NumType)
                for row in 1:param_size[1]
                    if !(constraint.lower_bound <= vp[s][getfield(ids, param)[row, col]] <= 1.0)
                        # Don't set the value to exactly the lower bound to avoid Inf
                        epsilon = (1.0 - constraint.lower_bound) * 1e-12
                        vp[s][getfield(ids, param)[row, col]] =
                            min(vp[s][getfield(ids, param)[row, col]], 1.0 - epsilon)
                        vp[s][getfield(ids, param)[row, col]] =
                            max(vp[s][getfield(ids, param)[row, col]],
                                    constraint.lower_bound + epsilon)
                    end
                    param_sum += vp[s][getfield(ids, param)[row, col]]
                end
                if param_sum != 1.0
                    # Normalize in a way that maintains the lower bounds
                    rescale =
                      (1 - constraint.n * constraint.lower_bound) /
                      (param_sum - constraint.n * constraint.lower_bound)
                    vp[s][getfield(ids, param)[:, col]] =
                      constraint.lower_bound +
                      rescale * (vp[s][getfield(ids, param)[:, col]] -
                                 constraint.lower_bound)
                end
            end
        end
    end
end

end


text:
#  ------------------------------------------------------------------------------------------
#  Copyright (c) Microsoft Corporation. All rights reserved.
#  Licensed under the MIT License (MIT). See LICENSE in the repo root for license information.
#  ------------------------------------------------------------------------------------------

create_batch_client(credentials::Nothing) = nothing

create_blob_client(credentials::Nothing) = nothing

create_queue_client(credentials::Nothing) = nothing

function create_clients(credentials::Nothing; batch::Bool=false, blob::Bool=false, queue::Bool=false)

    # List of all clients
    clients = Array{Any}(undef, 0)

    batch ? (batch_client = create_batch_client(credentials)) : batch_client = nothing
    blob ? (blob_client = create_blob_client(credentials)) : blob_client = nothing
    queue ? (queue_client = create_queue_client(credentials)) : queue_client = nothing

    return [Dict(
        "batch_client" => batch_client,
        "blob_client" => blob_client,
        "queue_client" => queue_client
        )]
end

###################################################################################################
# Blob stuff

# Create containers given a list of container names
create_blob_containers(blob_client::Nothing, container_name_list::Array{String, 1}) = nothing

create_batch_output_file(blob_client::Nothing, storage_account_name, container_name, filename) = nothing

###################################################################################################
# Batch stuff

# Create pool
function create_pool(batch_service_client::Nothing, pool_id, pool_vm_size, pool_node_count, node_os_publisher,
    node_os_offer, node_os_sku; image_resource_id=nothing, resource_files=nothing,
    enable_auto_scale=false, auto_scale_formula=nothing, auto_scale_evaluation_interval_minutes=nothing)

    return nothing
end

# Create pool and resource file
function create_pool_and_resource_file(clients::Dict{String,Nothing}, pool_id, pool_vm_size, pool_node_count,
    node_os_publisher, node_os_offer, node_os_sku, file_name, container_name; image_resource_id=nothing,
    enable_auto_scale=false, auto_scale_formula=nothing, auto_scale_evaluation_interval_minutes=15)

    return nothing
end

# Enable auto scaling for batch pool
enable_auto_scale(batch_client::Nothing, pool_id, auto_scale_formula; auto_scale_evaluation_interval_minutes=5) = nothing

# Resize pool
resize_pool(batch_client::Nothing, pool_id, target_dedicated_nodes, target_low_priority_nodes;
    resize_timeout_minutes=nothing, node_deallocation_option=nothing, pool_resize_options=nothing) = nothing

# Create resource file from file
create_batch_resource_from_file(blob_client::Nothing, container, file) = [nothing]

# Create resource file from existing blob
create_batch_resource_from_blob(blob_client::Nothing, container, blob) = [nothing]

# Create resource file from bytes
create_batch_resource_from_bytes(blob_client::Nothing, container, blob_name, blob) = [nothing]

# Create batch job
create_batch_job(batch_client::Nothing, job_id, pool_id; uses_task_dependencies=false, priority=0) = nothing

# Wait for all tasks to complete
wait_for_tasks_to_complete(batch_service_client::Nothing, job_id, timeout, verbose=true, num_restart=0) = true

# Wait for specified task to complete
wait_for_task_to_complete(batch_service_client::Nothing, job_id, task_id, timeout, verbose=true) = true

function wait_for_one_task_from_multi_pool(batch_service_client::Array{Nothing,1},
        job_id, task_id_list, timedelta_minutes, verbose=true)
    idx = randperm(length(task_id_list))[1]
    task_id = task_id_list[idx]["taskname"]    # return random index from task list
    pool_id = task_id_list[idx]["pool"]
    return task_id, pool_id
end

function wait_for_one_task_from_multi_jobs(batch_service_client::Array{Nothing,1},
        job_id_list, task_id_list, timedelta_minutes, verbose=true)
    idx = randperm(length(task_id_list))
    task_id = task_id_list[idx]["taskname"]    # return random index from task list
    pool_id = task_id_list[idx]["pool"]
    return task_id, pool_id
end


text:
using Documenter, LiquidsNEqProcesses

makedocs(;
    modules=[LiquidsNEqProcesses],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://gitlab.com/lanimfe/LiquidsNEqProcesses.jl/blob/{commit}{path}#L{line}",
    sitename="LiquidsNEqProcesses.jl",
    authors="Pablo Zubieta, Jesús B. Zepeda",
    assets=String[],
)

deploydocs(;
    repo="gitlab.com/lanimfe/LiquidsNEqProcesses.jl",
)


text:
# There are two derrivatives of interest in the following functions, variables
# differentiated with respect to time (t) are expressed as dot, and variables
# differentiated with respect to ficticious time (s) are expressed as prime
#
# For a variable V:
#   dV/dt = V_dot
#   dV/ds = V_prime
using StaticArrays
using LinearAlgebra
function L_fx(u)

    return @SMatrix [u[1] -u[2] -u[3]  u[4];
                     u[2]  u[1] -u[4] -u[3];
                     u[3]  u[4]  u[1]  u[2];
                     u[4] -u[3]  u[2] -u[1]]

end

function uprime_from_xdot(xdot,u)
    # given a u and xdot, return u prime
    # xdot = float(xdot)
    # u = float(u)

    u1 = u[1]
    u2 = u[2]
    u3 = u[3]
    u4 = u[4]

    xd1 = xdot[1]
    xd2 = xdot[2]
    xd3 = xdot[3]

    up1 = (1/2)*( u1*xd1  +  u2*xd2  +  u3*xd3  )
    up2 = (1/2)*(-u2*xd1  +  u1*xd2  +  u4*xd3  )
    up3 = (1/2)*(-u3*xd1  -  u4*xd2  +  u1*xd3  )
    up4 = (1/2)*( u4*xd1  -  u3*xd2  +  u2*xd3  )

    uprime = [up1;up2;up3;up4]

    return uprime

    # return SVector((1/2)*( u[1]*xd[1]  +  u[2]*xd[2]  +  u[3]*xd3  ))

end

function u_from_x(x)
    # get KS u from x, u1 is always 1
    x = float(x)
    r = norm(x)

    u1 = .1
    u4 = sqrt(.5*(r+x[1]) - u1^2)
    u2 = (x[2]*u1 + x[3]*u4)/(r+x[1])
    u3 = (x[3]*u1 - x[2]*u4)/(r+x[1])

    u = [u1;u2;u3;u4]
    return u
end

function x_from_u(u)
    # gets x from u

    return SVector( u[1]^2 - u[2]^2 - u[3]^2 + u[4]^2,
                     2*(u[1]*u[2] - u[3]*u[4]),
                     2*(u[1]*u[3] + u[2]*u[4]))

end

function xdot_from_u(u,u_prime)
    up = u_prime;
    u1 = u[1];u2 = u[2];u3 = u[3];u4 = u[4]
    up1 = up[1];up2 = up[2];up3 = up[3];up4 = up[4]
    r = dot(u,u)
    x1d = (2/r)*(u1*up1 - u2*up2 - u3*up3 + u4*up4)
    x2d = (2/r)*(u2*up1 + u1*up2 - u4*up3 - u3*up4)
    x3d = (2/r)*(u3*up1 + u4*up2 + u1*up3 + u2*up4)

    return [x1d;x2d;x3d]
end

function get_time_transfer(X,dt,tscale)
    """get the time in days for a run. This allows us to remove t from state"""
    # X = states(solver)
    t = 0.0
    t_hist = zeros(length(X))
    for i = 1:(length(X)-1)

        p0 = X[i][1:4]
        R0 = dot(p0,p0)

        p1 = X[i+1][1:4]
        R1 = dot(p1,p1)

        t += (R0 + R1)*dt/2
        t_hist[i+1] = copy(t)
    end
    t_days = t*tscale/(24*3600)
    t_hist = t_hist*tscale/(24*3600)
    return t_days, t_hist
end
function discretize_t(t_hist,n)
    """Take the time history vec, and discretize it even further"""
    # new time vec
    new_t = [0.0]
    for i = 1:(length(t_hist)-1)

        # time spacing betwen two knot points
        δt = t_hist[i+1] - t_hist[i]

        # add intermediate points between the knot points
        for j = 1:n
            push!(new_t,δt*j/n + t_hist[i])
        end
    end
    return new_t
end
function thrust_angle_from_rtn(u_rtn)
    u_rtn = normalize(u_rtn)
    α = atan(u_rtn[1],u_rtn[2])
    β = atan(u_rtn[3],norm(u_rtn[1:2]))
    return [α; β]
end
function get_rtn_in_out(U_real,r_eci_hist,v_eci_hist)
    """Get rtn thrust angles from eci thrust vector"""
    return [thrust_angle_from_rtn(rECItoRTN(
    [r_eci_hist[:,i];v_eci_hist[:,i]])*U_real[:,i] ) for i = 1:size(U_real,2)]
end

max_stars_repo_name: RoboticExplorationLab/KSLowThrust

max_issues_repo_name: RoboticExplorationLab/KSLowThrust

max_forks_repo_name: RoboticExplorationLab/KSLowThrust


text:
export ModelKitHomotopy

"""
    ModelKitHomotopy(H::Homotopy, parameters = nothing)

Construct a homotopy from the given homotopy `H` with the given `parameters` fixed.
The difference to `H` is that this compiles a straight line programm for the fast
evaluation of `H` and that `ModelKitHomotopy <: AbstractHomotopy`.
"""
struct ModelKitHomotopy{S,P<:Union{Nothing,AbstractVector}} <: AbstractHomotopy
    homotopy::ModelKit.CompiledHomotopy{S}
    parameters::P
end

ModelKitHomotopy(H::ModelKit.Homotopy, parameters = nothing) =
    ModelKitHomotopy(ModelKit.compile(H), parameters)
ModelKitHomotopy(H::ModelKit.CompiledHomotopy, parameters = nothing) =
    ModelKitHomotopy(H, parameters)

Base.size(H::ModelKitHomotopy) = size(H.homotopy)

evaluate!(u, H::ModelKitHomotopy, x, t) =
    ModelKit.evaluate!(u, H.homotopy, x, t, H.parameters)
evaluate_and_jacobian!(u, U, H::ModelKitHomotopy, x, t) =
    ModelKit.evaluate_and_jacobian!(u, U, H.homotopy, x, t, H.parameters)
function taylor!(u, v::Val, H::ModelKitHomotopy, x, t)
    ModelKit.taylor!(u, v, H.homotopy, x, t, H.parameters)
    u
end


text:
# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.4.2
#   kernelspec:
#     display_name: Julia 1.4.0
#     language: julia
#     name: julia-1.4
# ---

readdir()

using JLD2, FileIO, CSV, Plots
using FFTW

frame = CSV.read("frame.csv");
tmp1=frame[!,:KineticEnergy];
ftmp1=fft(tmp1)./50000;
ftmp1[1]=0.;
plot(abs.(ftmp1[1:1000]),xlabel="freq",ylabel="spectrum kinetic energy",legend=false)
#savefig("Kinetic_spectrum_dirac.png")
tmp2=frame[!,:PotentialEnergyE2];
ftmp2=fft(tmp2)./50000;
ftmp2[1]=0.;
#@show kmax=findmax(abs.(ftmp)[1:200])
#@show 77*2*3.1415/50.
plot(abs.(ftmp2[1:1000])./50000,xlabel="freq",ylabel="spectrum Ex energy",legend=false)
#savefig("Ex_spectrum_dirac.png")
tmp3=frame[!,:Kineticspin];
ftmp3=fft(tmp3)./50000;
ftmp3[1]=0.;
#@show kmax=findmax(abs.(ftmp)[1:200])
#@show 77*2*3.1415/50.
plot(abs.(ftmp3[1:1000]),xlabel="freq",ylabel="spectrum Spin energy",legend=false)
#savefig("Spin_spectrum_dirac.png")
tmp1=frame[!,:Momentum1];
ftmp1=fft(tmp1)./50000;
ftmp1[1]=0.;
#@show kmax=findmax(abs.(ftmp)[1:200])
#@show 77*2*3.1415/50.
plot(abs.(ftmp1[1:1000]),xlabel="freq",ylabel="Ex energy",legend=false)
#savefig("Ex_spectrum_dirac.png")
tmp1=frame[!,:PotentialEnergyE2];
ftmp1=fft(tmp1)./50000;
ftmp1[1]=0.;
#@show kmax=findmax(abs.(ftmp)[1:200])
#@show 77*2*3.1415/50.
plot(abs.(ftmp1[1:1000]),xlabel="freq",ylabel="Ey energy",legend=false)
#savefig("Ey_spectrum_dirac.png")
tmp1=frame[!,:PotentialEnergyE3];
ftmp1=fft(tmp1)./50000;
ftmp1[1]=0.;
#@show kmax=findmax(abs.(ftmp)[1:200])
#@show 77*2*3.1415/50.
plot(abs.(ftmp1[1:1000]),xlabel="freq",ylabel="Ez energy",legend=false)
#savefig("Ez_spectrum_dirac.png")
tmp1=frame[!,:PotentialEnergyB2];
ftmp1=fft(tmp1)./50000;
ftmp1[1]=0.;
#@show kmax=findmax(abs.(ftmp)[1:200])
#@show 77*2*3.1415/50.
plot(abs.(ftmp1[1:1000]),xlabel="freq",ylabel="By energy",legend=false)
#savefig("By_spectrum_dirac.png")
tmp1=frame[!,:PotentialEnergyB3];
ftmp1=fft(tmp1)./50000;
ftmp1[1]=0.;
#@show kmax=findmax(abs.(ftmp)[1:200])
#@show 77*2*3.1415/50.
plot(abs.(ftmp1[1:1000]),xlabel="freq",ylabel="Bz energy",legend=false)
#savefig("Bz_spectrum_dirac.png")
# +
#Kineticspin,Momentum1,Momentum2,Momentum3,Momentum4,
#Momentum5,Momentum6,Momentum7,PotentialEnergyE1,
#PotentialEnergyE2,PotentialEnergyE3,PotentialEnergyB2,PotentialEnergyB3
# -