text:
"""
    gradient_optimizer(obj_fun, x0; kwargs...)

Performs a very rudimentary version of gradient descent/ascent. Follows basic
pattern of NLopt package: `obj_fun` is a function which accepts arguments `x`, `g`; it
returns the objective function at `x` and sets the gradient `g` in place. The
value `x0` is an initial guess. Default behavior is to descend along the
gradient to minimize the function, setting keyword argument `objective` to
`:max` changes this behavior (default is `:min`). Returns a tuple
`(F_opt, x_opt, stop)` where `stop` is an array of the stopping criteria met.

List of keywords and default values:
 * `objective = :min` minimize `obj_fun`, setting to `:max` will maximize
 * `verbose = 0` Options are 0 (no output), 1 (some output), 2 (all the output)
 * `print_eval = 100` In maximum verbosity, waits this many evaluations before printing
 * `lr = 1.0` A "learning rate", the gradient is scaled by this much.
 * `adjust_lr = lr > 1` If true, every time a smaller value of `obj_fun` is reached, `lr` is halved.

Stopping criteria (except vanishing gradient) can be set with keyword arguments.
 * `maxeval = 1000` stop after evaluating `obj_fun` this many times
 * `xtol_abs = 0.0` stop when `|x_cur - x_prev| < xtol_abs`
 * `xtol_rel = 0.0` stop when `|x_cur - x_prev|/|x_cur| < xtol_rel`
 * `ftol_abs = 0.0` stop when `|F_cur - F_prev| < ftol_abs`
 * `ftol_rel = 0.0` stop when `|F_cur - F_prev|/|F_cur| < ftol_rel`

"""
function gradient_optimizer(obj_fun, x0; objective=:min,
    verbose::Int=0, print_eval::Int=100,
    lr::Float64=1.0, adjust_lr::Bool=(lr > 1.0),
    maxeval::Int=1000,
    xtol_abs::Float64=0.0, xtol_rel::Float64=0.0,
    ftol_abs::Float64=0.0, ftol_rel::Float64=0.0,
    kwargs...)

    # println("The type of obj_fun is $(typeof(obj_fun))")
    # obj_fun = (x,g) -> obj_fun(x,g)
    # println("Now the type of obj_fun is $(typeof(obj_fun))")
    # bookkeeping set up
    n_vars = length(x0)
    dkwargs = Dict(kwargs)
    obj_sign = (objective == :min ? -1 : 1)
    stopping_criteria = [:maxeval_reached,
                        :opt_reached,
                        :xtol_abs_reached,
                        :xtol_rel_reached,
                        :ftol_abs_reached,
                        :ftol_rel_reached]
    stop = falses(stopping_criteria)

    # gradient algorithm set up
    g_cur = zeros(x0); g_prev = zeros(x0); g_opt = zeros(x0)
    x_cur = zeros(x0); x_prev = zeros(x0); x_opt = zeros(x0)
    F_cur = 0.0; F_prev = 0.0; F_opt = -obj_sign * Inf;
    dF_cur = 0.0; dx_cur = 0.0; dx_next = 0.0
    eval = 0; eval_opt = 0
    x_prev[:] = x0[:]

    whoopscount = 0

    if verbose > 0
        println("gradient_optimizer: $(obj_sign > 0 ? "maximizing" : "minimizing") function of $n_vars variable$(n_vars > 1 ? "s" : "")")
    end
    if verbose > 1
        println("\tlr           : $lr ($(adjust_lr ? "will adjust" : "no adjust"))")
        println("\tmaxeval      : $maxeval")
        println("\txtol abs/rel : $xtol_abs / $xtol_rel")
        println("\tftol abs/rel : $ftol_abs / $ftol_rel")
        println()
        println("\teval\tF(x_cur)\t\tdF\t|dx_next|\tlr")
    end
    while stop.chunks[1] == 0
        # get current value of x_cur from x_prev and gradient
        x_cur[:] = x_prev[:] + obj_sign * lr * g_prev[:]
        # evaluate at new x_cur
        F_cur = obj_fun(x_cur, g_cur)
        dF_cur = F_cur - F_prev
        dx_next = lr * norm(g_cur) # size of next step

        # update local extremum
        if F_cur * obj_sign > F_opt * obj_sign
            # if we're more extreme than current best, update
            g_opt[:] = g_cur[:]
            x_opt[:] = x_cur[:]
            F_opt = F_cur
            eval_opt = eval
        end

        # print if necessary
        if verbose > 1 && mod(eval, print_eval) == 0
            println("\t$eval\t$F_cur\t$dF_cur\t$dx_next\t$lr$(eval_opt == eval ? " *" : "")")
        end

        # update the learning rate if we overshot
        if adjust_lr && (dF_cur * obj_sign < 0)
            whoopscount += 1
            lr /= 2
        end

        # increment counter; eval is now the index of the next step.
        eval += 1

        #push cur into prev slot
        x_prev[:] = x_cur[:]
        g_prev[:] = g_cur[:]
        F_prev = F_cur

        # Check stopping conditions
        # 1 maxeval
        # 2 opt reached (dF = 0)
        # 3 xtol_abs reached (lr * norm(g_cur) < xtol_abs)
        # 4 xtol_rel reached (lr * norm(g_cur) / norm(x_cur) < xtol_rel)
        # 5 ftol_abs reached (dF_cur < ftol_abs)
        # 6 ftol_rel reached (|dF_cur / F_cur| < ftol_rel )
        Fstep_rel = F_cur == 0 ? Inf : abs(dF_cur / F_cur)
        xstep_rel = norm(x_cur) == 0 ? Inf : abs(dx_next / norm(x_cur))
        stop = stop | [eval > maxeval,
                        dF_cur == 0.0,
                        dx_next < xtol_abs,
                        xstep_rel < xtol_rel,
                        abs(dF_cur) < ftol_abs,
                        Fstep_rel < ftol_rel]
    end
    if verbose > 0
        println("gradient_optimizer: Stopping criteria reached: $(stopping_criteria[stop])")
        println("gradient_optimizer: optimal value found: $F_opt")
    end
    if verbose > 1
        println("gradient_optimizer: opt. value reached at evaluation $eval_opt / $maxeval")
        println("gradient_optimizer: Learning rate adjusted $whoopscount time$(whoopscount != 1 ? "s" : "")")
        println("gradient_optimizer: final learning rate: $lr")
    end

    return F_opt, x_opt, stopping_criteria[stop]
end

"""
    brute_force_optimizer(obj_fun, d; objective=:min, kwargs...)

Perform a brute force search over a `d`-dimensional cube for the optimal value
of `obj_fun`. Returns `F_opt, x_opt, [:brute_force_search]` to mirror output
from `gradient_optimizer`

Some kwargs that can be passed:
 * `ranges::Vector{Iterable}` A vector of `d` iterables explicitly defining the search space
 * `save_fun::Array{Float64,d}` a matrix that will be modified in-place with the function values over the whole cube. Size must match the `length` of each iterable in `ranges`

WARNING: Currently only accepts functions which take exactly 3 parameters,
because that's all I need it for. And let's be real, you shouldn't be doing this
for large `d` anyway. (The further reason I am not generalizing is that I'll be
using this on a machine running Julia v0.5, so no Itertools package, and don't
feel like writing my own product iterator.)

Current functionality is limited: Must pass `ranges = [iter1, iter2, iter3]` as
a kwarg, this will define the cube in which search takes place. Each `iter` can
be any iterator that is indexed by a range of the form `1:N` and returns `Real`
numbers (and reports `length(iter) = N`)

"""
function brute_force_optimizer(obj_fun, d::Integer; verbose=0,
    objective=:min, kwargs...)
    if d == 3
        return brute_force_optimizer3(obj_fun, d; verbose=verbose, objective=objective, kwargs...)
    else
        error("brute_force_optimizer: Can't actually handle dimensionality other than 3 right now, sorry.")
        #MAYBEDO expand to any dimensionality
    end
end
function brute_force_optimizer3(obj_fun, d; verbose=0,
    objective=:min, save_fun::Array{Float64}=[], kwargs...)

    d = 3

    dkwargs = Dict(kwargs)
    obj_sign = (objective == :min) ? -1 : 1
    # lb = length(lb) < d ? lb * ones(d) : [lb...]
    # ub = lengtH(ub) < d ? ub * ones(d) : [ub...]

    #MAYBEDO implement lots of options to allow automated creation of ranges
    #from, e.g. lb and ub.

    # for now, make life easy and assume explicit ranges have been passed as
    # `ranges` kwarg.
    x_opt = zeros(d)
    F_opt = -obj_sign * Inf
    F_cur = 0.0
    ranges = pop!(dkwargs, :ranges)
    if verbose > 0
        println("brute_force_optimizer: $(obj_sign > 0 ? "maximizing" : "minimizing") function of $d variables")
        for k = 1:d
            println("\trange $k: $(length(ranges[k])) options")
        end
        if length(save_fun) > 0
            println("\tWriting function values to matrix")
        end
    end
    for (i1,x1) in enumerate(ranges[1])
        for (i2,x2) in enumerate(ranges[2])
            for (i3,x3) in enumerate(ranges[3])
                # score[i,j,k] = obj_fun([x1,x2,x3])
                F_cur = obj_fun([x1,x2,x3])
                if F_cur * obj_sign > F_opt * obj_sign
                    F_opt = F_cur
                    x_opt = [x1,x2,x3]
                end
                if length(save_fun) > 0
                    save_fun[i1,i2,i3] = F_cur
                end
            end
        end
    end
    return F_opt, x_opt, [:brute_force_search]
end


text:
module CUnion

using MacroTools

abstract type AbstractUnion end

function unionfields end

export AbstractUnion,unionfields

function reinterpret_cast(::Type{T}, x) where T
    s = max(sizeof(x), sizeof(T))
    a = Vector{UInt8}(undef, s)
    p = pointer(a)
    unsafe_store!(Ptr{typeof(x)}(p), x)
    unsafe_load(Ptr{T}(p))
end

function Base.getproperty(u::T, name::Symbol) where T <: AbstractUnion
    fields = unionfields(T)
    i = findfirst(((n, t),)->n == name, fields)
    if !isnothing(i)
        ft = fields[i].second
        isprimitivetype(ft) && sizeof(ft) == sizeof(T) && return reinterpret(ft, u)
        return reinterpret_cast(ft, u)
    end
    getfield(u, name)
end

function unionof(::Type{U}, x::T) where {U <: AbstractUnion,T}
    types = tuple(unique(last.(unionfields(U)))...)
    ret = findfirst(t->T == t, types)
    if isnothing(ret)
        error("$(typeof(x)) does not match any types of fields: $types")
    end
    if isprimitivetype(T) && sizeof(T) == sizeof(U)
        reinterpret(U, x)
    else
        reinterpret_cast(U, x)
    end
end

walkctx(x, inner, outer, ctx) = outer(x, ctx)
walkctx(x::Expr, inner, outer, ctx) = outer(Expr(x.head, map(inner, x.args)...), ctx)
postwalkctx(f, x, ctx) = walkctx(x, x->postwalkctx(f, x, isexpr(x, :struct) ? "$ctx#$(x.args[2])" : ctx), f, ctx)

function structflatten(lines, basety)
    defs = []
    result = postwalkctx(lines, "#$basety") do ex, ctx
        if @capture(ex,struct fname_Symbol lines__ end)
            fty = Symbol(ctx)
            def = Expr(:struct, false, fty, Expr(:block, lines...))
            push!(defs, def)
            :($fname::$fty)
        else
            ex
        end
    end
    return (Expr(:block, defs...) |> MacroTools.flatten, result)
end

macro union(ex)
    @capture(ex,
    struct T_Symbol
        rawbody_
    end
    ) || error("expected a struct definition without parameters: $ex")
    (defs, body) = structflatten(rawbody, T)
    @capture(body,lines__)
    members = map(lines) do line
        @capture(line,name_::type_) || error("expected a typed field definition: $line")
        (name, type)
    end |> x->(x...,)
    membersex = Expr(:tuple, map(members) do (name, type)
        :($(QuoteNode(name)) => $(type |> esc))
    end...)
    names = first.(members)
    if length(names) != length(unique(names))
        error("duplicate field name")
    end
    types = last.(members)
    typesex = Expr(:tuple, map(types) do t
        t |> esc
    end...)
    modulesym(x::Symbol) = Expr(:., Symbol(@__MODULE__) |> esc, QuoteNode(x))
    quote
        $(defs |> esc)
        maxsize = maximum($typesex) do t
            isbitstype(t) || error("not bits type: $t")
            sizeof(t)
        end * 8
        primitive type $T <: AbstractUnion maxsize end
        $(modulesym(:unionfields))(::Type{$(T |> esc)}) = $membersex
        $(T |> esc)($(:x |> esc)) = unionof($(T |> esc), $(:x |> esc))
        nothing
    end |> x->MacroTools.postwalk(unblock ∘ rmlines, x)
end

export @union

end # module

ext: jl


text:
using CVChannel
using Test
using DelimitedFiles

"""
This obtains data about the Werner-Holevo channel when tensored
with the Siddhu channel.
"""

print("\nThis script shows that there is non-multiplicativity of")
print("cvPPT for the Werner-Holevo channel tensored with the Siddhu channel.")

@testset "Werner-Holevo with Siddhu is Non-Multiplicative" begin
    println("Now we actually verify non-multiplicativity.")
    λ_range = [0:0.025:0.325;]
    s_range = [0:0.1:0.5;]
    λ_ctr, s_ctr = 1,1
    data_table = zeros(length(s_range),length(λ_range))
    for λ in λ_range
        println("---Now scanning λ=",λ,"---")
        A,B,g,a = generalWHLPConstraints(1,3,λ*ones(1))
        cvWH, v1 =wernerHolevoCVPPT(1,3,A,B,g,a)
        whChan(X) = wernerHolevoChannel(X,λ)
        wh_chan = Choi(whChan,3,3)
        for s in s_range
            println("Now scanning s=",s)
            target_val = 2*cvWH
            sidchan(X) = siddhuChannel(X,s)
            sid_chan= Choi(sidchan,3,3)

            par_cv, = pptCV(parChoi(wh_chan, sid_chan), :dual)

            non_mult = par_cv - target_val
            isapprox(non_mult,0,atol=3e-6) ?
                data_table[s_ctr,λ_ctr] = 0 :
                    data_table[s_ctr,λ_ctr] = non_mult

            s_ctr += 1
            end
        λ_ctr += 1
        s_ctr = 1
    end

    println("\n Here is the data.")
    show(stdout, "text/plain", data_table)
    println("Note that the non-multiplicativity is the same for all s.")
    println("Also note that as λ increases, the non-multiplicativity decreases.")

    header, s_label = hcat("s|λ:",λ_range'), s_range
    data_to_save = vcat(header,hcat(s_label, data_table))

    info_vec = Vector{Union{Nothing,String}}(nothing, size(data_to_save)[1])
    info_vec[1] = "INFO:"
    info_vec[2] = "Generated by werner-holevo-with-siddhu.jl"
    info_vec[4] = "λ_range ="*string(λ_range)
    info_vec[5] = "s_range="*string(s_range)
    data_to_save = hcat(data_to_save,info_vec)

    println("Please name the file you'd like to write the results to: \n")
    file_name = readline()
    file_to_open = string(file_name,".csv")
    writedlm(file_to_open, data_to_save, ',')

    #This is sufficient for checking that there is non-multiplicativity somewhere
    @test data_table[1:1:1]>[0.1]
end

max_forks_repo_name: ChitambarLab/cv-channel


text:
"A fixed-length vector whose data is stored in a Blob."
struct BlobVector{T} <: AbstractArray{T, 1}
    data::Blob{T}
    length::Int64
end

Base.@propagate_inbounds function get_address(blob::BlobVector{T}, i::Int)::Blob{T} where T
    @boundscheck begin
        (0 < i <= blob.length) || throw(BoundsError(blob, i))
    end
    blob.data + (i-1)*self_size(T)
end

# blob interface

function Base.size(blob::Blob{BlobVector})
    (blob.length[],)
end

function Base.IndexStyle(_::Type{Blob{BlobVector{T}}}) where T
    Base.IndexLinear()
end

Base.@propagate_inbounds function Base.getindex(blob::Blob{BlobVector{T}}, i::Int)::Blob{T} where T
    get_address(blob[], i)
end

# array interface

function Base.size(blob::BlobVector)
    (blob.length,)
end

function Base.IndexStyle(_::Type{BlobVector{T}}) where T
    Base.IndexLinear()
end

Base.@propagate_inbounds function Base.getindex(blob::BlobVector{T}, i::Int)::T where T
    get_address(blob, i)[]
end

Base.@propagate_inbounds function Base.setindex!(blob::BlobVector{T}, v, i::Int)::T where T
    get_address(blob, i)[] = v
end

# copying, with correct handling of overlapping regions
# TODO use memcopy
function Base.copy!(dest::BlobVector{T}, doff::Int, src::BlobVector{T}, soff::Int, n::Int) where T
    if doff < soff
        for i in 0:n-1 dest[doff+i] = src[soff+i] end
    else
        for i in n-1:-1:0 dest[doff+i] = src[soff+i] end
    end
end

# iterate interface

@inline function Base.iterate(blob::BlobVector, i=1)
    (i % UInt) - 1 < length(blob) ? (@inbounds blob[i], i + 1) : nothing
end


text:
function solve(
    i::UniformMatroidInstance{Float64, Maximise},
    ::DynamicProgramming,
)
    # V[µ, δ]: µ items, up to δ.
    V = Matrix{Float64}(undef, i.m, dimension(i))
    S = Dict{Tuple{Int, Int}, Vector{Int}}()

    # Initialise: µ == 1, just take the best element among [δ+1, d]; δ == d, take nothing.
    V[1, dimension(i)] = 0.0
    S[1, dimension(i)] = Int[]

    for δ in (dimension(i) - 1):-1:1
        v, x = findmax(i.rewards[(δ + 1):dimension(i)])
        V[1, δ] = v
        S[1, δ] = [x]
    end

    for µ in 2:(i.m)
        V[µ, dimension(i)] = 0.0
        S[µ, dimension(i)] = Int[]
    end

    # Dynamic part.
    for µ in 2:(i.m)
        for δ in (dimension(i) - 1):-1:0
            obj_idx = δ + 1
            take_δ = i.rewards[obj_idx] + V[µ - 1, δ + 1]
            dont_take_δ = V[µ, δ + 1]

            if take_δ > dont_take_δ
                V[µ, δ + 1] = take_δ
                S[µ, δ + 1] = vcat(δ, S[µ - 1, δ + 1])
            else
                V[µ, δ + 1] = dont_take_δ
                S[µ, δ + 1] = S[µ, δ + 1]
            end
        end
    end

    return UniformMatroidSolution(i, S[i.m, 1] .+ 1)
end


text:
export AbstractComparer, Comparer, fn, reset!

abstract AbstractComparer

immutable Comparer <: AbstractComparer
    fn::Function
    reset::Nullable{Function}
end

Comparer() = Comparer((x, y) -> <(x, y), Nullable{Function}())

Comparer(fn::Function) = Comparer(fn, Nullable{Function}())

Comparer(ac::AbstractComparer) = Comparer(fn(ac), () -> reset!(ac))

fn(c::Comparer) = c.fn

reset!(c::Comparer) = !isnull(c.reset) && c.reset()

max_stars_repo_name: unbornchikken/julia-ml-proto

max_issues_repo_name: unbornchikken/julia-ml-proto

max_forks_repo_name: unbornchikken/julia-ml-proto


text:
# This is quasi derived from the Rosenbrock function
# whereby F(x,p) is almost the gradient of the Rosenbrock function
# except the Hessian matroix is made non-symmetric by slightly altering
# F(x,p)

# State function ≈ gradient of Rosenbrock.
# With p = [a, b] the parameters defining the constants
# in Rosenbrocks function, and the slight alteration to make the
# Hessian matrix non-symmetric.
# Specifically, if r(x,y) denotes the rosenbrock function,
#     F₁([x, y], p) = ∂r/∂x
#     F₂([x, y], p) = 0.5 * ∂r/∂y
F(x,p) = [
    -2 * (p[1] - x[1]) - 4 * p[2] * (x[2] - x[1]^2) * x[1]
    p[2] * (x[2] - x[1]^2)
]

# Jacobian function of F w.r.t. p
∇ₓF(x,p) = ForwardDiff.jacobian(x -> F(x,p), x)

# Define mismatch function f(x,p) and its derivative ∇ₓf(x,p)
# (Note ∇ₓF and ∇ₓf are required by the F1 method)
function state_mismatch(x)
    δ(x) = x .- 1
    return 0.5δ(x)'δ(x)
end
function parameter_mismatch(p)
    δ(p) = log.(p)
    return 0.5δ(p)'δ(p)
end
f(x,p) = state_mismatch(x) + parameter_mismatch(p)
∇ₓf(x,p) = ForwardDiff.jacobian(x -> [f(x,p)], x)

# Initialize the starting state and parameters
x₀ = rand(2)
p₀ = rand(2)
# Initialize the memory cache for storing reusable objects
mem = F1Method.initialize_mem(x₀, p₀)

# Define the functions via the F1 method
F1_objective(p) = F1Method.objective(f, F, ∇ₓF, mem, p, MyAlg())
F1_gradient(p) = F1Method.gradient(f, F, ∇ₓf, ∇ₓF, mem, p, MyAlg())
F1_Hessian(p) = F1Method.hessian(f, F, ∇ₓf, ∇ₓF, mem, p, MyAlg())

# Define the exact solution and objective for comparison
exact_solution(p) = [p[1], p[1]^2]
exact_objective(p) = f(exact_solution(p), p)

# Run the tests
@testset "objective" begin
    @test exact_objective(p₀) ≈ F1_objective(p₀)
end
@testset "gradient" begin
    @test ForwardDiff.jacobian(p -> [exact_objective(p)], p₀) ≈ F1_gradient(p₀)
end
@testset "Hessian" begin
    @test ForwardDiff.hessian(p -> exact_objective(p), p₀) ≈ F1_Hessian(p₀)
end