text:
using RecipesBase, Measures
export IexMcmcSampler, AuxiliaryMcmc, PosteriorMcmcOutput

struct AuxiliaryMcmc{S,T}
    mcmc::S  # Mcmc move (leave general since will also use for graph models)
    data::Vector{T} # To store auxiliary samples
    init_at_prev::Bool
end

struct IexMcmcSampler{T<:InvMcmcMove,S<:AuxiliaryMcmc}
    move::T    # Move for mode update
    ε::Float64 # Neighborhood for gamma proposal
    aux::S
    desired_samples::Int
    burn_in::Int
    lag::Int
    K::Int
    pointers::InteractionSequence{Int}
    suff_stats::Vector{Float64}
    γ_counts::Vector{Int}
    function IexMcmcSampler(
        move::T, aux_mcmc::InvMcmcSampler;
        ε::Float64=0.15,
        desired_samples::Int=1000, burn_in::Int=0, lag::Int=1,
        K::Int=100, aux_init_at_prev::Bool=false
        ) where {T<:InvMcmcMove}
        pointers = [Int[] for i in 1:(2K)]
        suff_stats = [0.0,0.0]
        aux_data = [Vector{Int}[]]
        γ_counts = [0,0]
        aux = AuxiliaryMcmc(aux_mcmc, aux_data, aux_init_at_prev)
        new{T,typeof(aux)}(
            move, ε, aux,
            desired_samples, burn_in, lag,
            K, pointers,
            suff_stats,
            γ_counts
            )
    end
end

Base.show(io::IO, x::IexMcmcSampler) = print(io, typeof(x))

struct PosteriorMcmcOutput{T<:Union{SisPosterior,SimPosterior}}
    S_sample::Vector{InteractionSequence{Int}}
    γ_sample::Vector{Float64}
    suff_stat_trace::Vector{Float64}
    posterior::T
end

Base.show(io::IO, x::PosteriorMcmcOutput) = print(io, typeof(x))

@recipe function f(
    output::PosteriorMcmcOutput,
    S_true::InteractionSequence{Int}
    )

    S_sample = output.S_sample
    γ_sample = output.γ_sample
    layout := (2,1)
    legend --> false
    xguide --> "Index"
    yguide --> ["Distance from True Mode" "γ"]
    size --> (800, 600)
    margin --> 5mm
    y1 = map(x->output.posterior.dist(S_true,x), S_sample)
    y2 = γ_sample
    hcat(y1,y2)
end

acceptance_prob(mcmc::IexMcmcSampler) = (mode=acceptance_prob(mcmc.move),γ=mcmc.γ_counts[1]/mcmc.γ_counts[2])

function eval_accept_prob(
    S_curr::InteractionSequence{Int},
    S_prop::InteractionSequence{Int},
    γ_curr::Float64,
    aux::AuxiliaryMcmc,
    posterior::T,
    log_ratio::Float64,
    suff_stats::Vector{Float64}
    ) where {T<:Union{SisPosterior,SimPosterior}}

    dist, V, K_inner, K_outer = (
        posterior.dist, posterior.V,
        posterior.K_inner, posterior.K_outer
    )
    mode_prior, γ_prior = (
        posterior.S_prior.mode, posterior.S_prior.γ
    )

    # Infer aux model from mode prior
    aux_model = typeof(posterior.S_prior)(
        S_prop, γ_curr,
        dist, V,
        K_inner, K_outer
    )

    if aux.init_at_prev
        @inbounds tmp = deepcopy(aux.data[end])
        draw_sample!(aux.data, aux.mcmc, aux_model, init=tmp)
    else
        draw_sample!(aux.data, aux.mcmc, aux_model)
    end

    aux_log_lik_ratio = -γ_curr * (
        sum(dist(x, S_curr) for x in aux.data)
        - sum(dist(x, S_prop) for x in aux.data)
    )

    suff_stats[2] = sum(dist(x, S_prop) for x in posterior.data)

    log_lik_ratio = -γ_curr * (
        suff_stats[2] - suff_stats[1] # (prop-curr)
    )

    log_prior_ratio = -γ_prior * (
        dist(S_prop, mode_prior) - dist(S_curr, mode_prior)
    )

    log_α = log_lik_ratio + log_prior_ratio + aux_log_lik_ratio + log_ratio

    # Log acceptance probability
    if typeof(posterior)==SisPosterior
        return log_α
    else
        log_multinom_term = log_multinomial_ratio(S_curr, S_prop)
        return log_α + log_multinom_term
    end

end

function accept_reject_mode!(
    S_curr::InteractionSequence{Int},
    S_prop::InteractionSequence{Int},
    γ_curr::Float64,
    pointers::InteractionSequence{Int},
    move::InvMcmcMove,
    aux::AuxiliaryMcmc,
    posterior::S,
    suff_stats::Vector{Float64}
    ) where {S<:Union{SisPosterior,SimPosterior}}

    move.counts[2] += 1

    # Do imcmc proposal (log_ratio gets passed to eval_accept_prob()...)
    log_ratio = prop_sample!(S_curr, S_prop, move, pointers, posterior.V)

    # Adjust for dimension bounds (reject if outside of them)
    if any(!(1 ≤ length(x) ≤ posterior.K_inner.u) for x in S_prop)
        log_α = -Inf
    elseif !(1 ≤ length(S_prop) ≤ posterior.K_outer.u)
        log_α = -Inf
    else
        log_α = eval_accept_prob(
            S_curr, S_prop, γ_curr,
            aux,
            posterior,
            log_ratio,
            suff_stats
        )
    end

    if log(rand()) < log_α
        # We accept!
        move.counts[1] += 1
        suff_stats[1] = suff_stats[2]  # Update sufficient stat
        enact_accept!(S_curr, S_prop, pointers, move)
    else
        # We reject!
        enact_reject!(S_curr, S_prop, pointers, move)
    end

end

function accept_reject_mode!(
    S_curr::InteractionSequence{Int},
    S_prop::InteractionSequence{Int},
    γ_curr::Float64,
    pointers::InteractionSequence{Int},
    mix_move::InvMcmcMixtureMove,
    aux::AuxiliaryMcmc,
    posterior::T,
    suff_stats::Vector{Float64}
    ) where {T<:Union{SisPosterior,SimPosterior}}

    # Sample move
    p = mix_move.p  # Mixture proportions
    β = rand()      # Random unif(0,1)
    z, i = (0.0, 0) # To store cumulative prob and index
    for prob in p
        if z > β
            break
        end
        i += 1
        z += prob
    end

    # Select ith move
    move = mix_move.moves[i]
    # Do accept-reject for the move
    accept_reject_mode!(
        S_curr, S_prop, γ_curr,
        pointers,
        move,
        aux,
        posterior,
        suff_stats
    )
end

function accept_reject_mode!(
    S_curr::InteractionSequence{Int},
    S_prop::InteractionSequence{Int},
    γ_curr::Float64,
    pointers::InteractionSequence{Int},
    mcmc::IexMcmcSampler,
    posterior::T
    ) where {T<:Union{SisPosterior,SimPosterior}}
    accept_reject_mode!(
        S_curr, S_prop, γ_curr,
        pointers,
        mcmc.move,
        mcmc.aux,
        posterior,
        mcmc.suff_stats
    )
end

function initialise!(
    aux::AuxiliaryMcmc{S,T},
    S_init::T, # Initial value for aux chain
    γ_init::Float64, # Used to set-up auxiliary model
    S_prior::V,
    n::Int # Sample size
    ) where {T,S,V<:Union{SIS,SIM}}

    if length(aux.data) ≥ n
        resize!(aux.data, n)
    else
        append!(aux.data, [similar(S_init) for i in 1:(n-length(aux.data))])
    end
    aux_model = similar(
        S_prior,
        S_init,
        γ_init
    )
    draw_sample!(aux.data, aux.mcmc, aux_model)

end

function draw_sample_mode!(
    sample_out::Union{InteractionSequenceSample{Int}, SubArray},
    mcmc::IexMcmcSampler,
    posterior::T,
    γ_fixed::Float64;
    burn_in::Int=mcmc.burn_in,
    lag::Int=mcmc.lag,
    S_init::Vector{Path{Int}}=sample_frechet_mean(posterior.data, posterior.dist)[1],
    loading_bar::Bool=true
    ) where {T<:Union{SisPosterior,SimPosterior}}

    if loading_bar
        iter = Progress(
            length(sample_out)*lag + burn_in, # How many iters
            1,  # At which granularity to update loading bar
            "Chain for γ = $(γ_fixed) and n = $(posterior.sample_size) (mode conditional)....")  # Loading bar. Minimum update interval: 1 second
    end

    pointers = mcmc.pointers
    S_curr, S_prop = intialise_states!(pointers, S_init)

    sample_count = 1    # Keeps which sample to be stored we are working to get
    i = 0               # Keeps track all samples (included lags and burn_ins)
    reset_counts!(mcmc.move) # Reset counts for mcmc move (tracks acceptances)

    # Initialise aux data
    initialise!(mcmc.aux, S_curr, γ_fixed, posterior.S_prior, posterior.sample_size)

    # Initialise sufficient statistic
    dist, data = (posterior.dist, posterior.data)
    suff_stats = mcmc.suff_stats
    suff_stats[1] = sum(dist(S_curr, x) for x in data)

    # Storage for sample statistic
    sample_suff_stat = zeros(length(sample_out))

    while sample_count ≤ length(sample_out)
        i += 1
        # Store value
        if (i > burn_in) & (((i-1) % lag)==0)
            @inbounds sample_out[sample_count] = deepcopy(S_curr)
            @inbounds sample_suff_stat[sample_count] = suff_stats[1]
            sample_count += 1
        end
        accept_reject_mode!(
            S_curr, S_prop,
            γ_fixed, pointers,
            mcmc, posterior
        )
        # println(S_prop)
        if loading_bar
            next!(iter)
        end
    end
    return_states!(S_curr, S_prop, pointers)
    return sample_suff_stat
end

function draw_sample_mode(
    mcmc::IexMcmcSampler,
    posterior::T,
    γ_fixed::Float64;
    desired_samples::Int=mcmc.desired_samples,
    args...
    ) where {T<:Union{SisPosterior,SimPosterior}}

    sample_out = [[Int[]] for i in 1:desired_samples]
    sample_suff_stat = draw_sample_mode!(
        sample_out,
        mcmc, posterior, γ_fixed,
        args...
    )

    return sample_out, sample_suff_stat

end

function (mcmc::IexMcmcSampler)(
    posterior::T,
    γ_fixed::Float64;
    args...
    ) where {T<:Union{SisPosterior,SimPosterior}}

    S_sample, suff_stat_trace = draw_sample_mode(
        mcmc, posterior, γ_fixed;
        args...
    )
    γ_sample = fill(γ_fixed, length(S_sample))

    return PosteriorMcmcOutput(
        S_sample, γ_sample,
        suff_stat_trace,
        posterior
    )
end

function accept_reject_gamma(
    γ_curr::Float64,
    S_curr::InteractionSequence{Int},
    mcmc::IexMcmcSampler,
    posterior::T
    ) where {T<:Union{SisPosterior,SimPosterior}}

    dist, V, K_inner, K_outer = (
        posterior.dist, posterior.V,
        posterior.K_inner, posterior.K_outer
    )
    aux = mcmc.aux
    ε, suff_stats = (
        mcmc.ε,
        mcmc.suff_stats
    )

    γ_prop = rand_reflect(γ_curr, ε, 0.0, Inf)

    aux_model = typeof(posterior.S_prior)(
        S_curr, γ_prop,
        dist, V,
        K_inner, K_outer
    )

    if aux.init_at_prev
        @inbounds tmp = deepcopy(aux.data[end])
        draw_sample!(aux.data, aux.mcmc, aux_model, init=tmp)
    else
        draw_sample!(aux.data, aux.mcmc, aux_model)
    end

    # Accept reject
    log_lik_ratio = (γ_curr - γ_prop) * suff_stats[1]
    aux_log_lik_ratio = (γ_prop - γ_curr) * sum(dist(x,S_curr) for x in aux.data)

    log_α = (
        logpdf(posterior.γ_prior, γ_prop)
        - logpdf(posterior.γ_prior, γ_curr)
        + log_lik_ratio + aux_log_lik_ratio
    )
    if log(rand()) < log_α
        return γ_prop, 1
    else
        return γ_curr, 0
    end
end

function draw_sample_gamma!(
    sample_out::Union{Vector{Float64}, SubArray},
    mcmc::IexMcmcSampler,
    posterior::T,
    S_fixed::InteractionSequence{Int};
    burn_in::Int=mcmc.burn_in,
    lag::Int=mcmc.lag,
    γ_init::Float64=4.0,
    loading_bar::Bool=true
    ) where {T<:Union{SisPosterior,SimPosterior}}

    if loading_bar
        iter = Progress(
            length(sample_out)*lag + burn_in, # How many iters
            1,  # At which granularity to update loading bar
            "Chain for n = $(posterior.sample_size) (disperison conditional)....")  # Loading bar. Minimum update interval: 1 second
    end

    γ_curr = γ_init
    S_curr = deepcopy(S_fixed)

    sample_count = 1    # Keeps which sample to be stored we are working to get
    i = 0               # Keeps track all samples (included lags and burn_ins)

    # Initialise aux data
    initialise!(mcmc.aux, S_curr, γ_curr, posterior.S_prior, posterior.sample_size)

    # Initialise sufficient statistic
    dist, data = (posterior.dist, posterior.data)
    suff_stats = mcmc.suff_stats
    suff_stats[1] = sum(dist(S_curr, x) for x in data)

    acc_count = 0
    i = 0 # Counter for iterations
    sample_count = 1  # Which sample we are working to get

    while sample_count ≤ length(sample_out)
        i += 1
        # Store value
        if (i > burn_in) & (((i-1) % lag)==0)
            sample_out[sample_count] = γ_curr
            sample_count += 1
        end

        γ_curr, was_acc = accept_reject_gamma(
            γ_curr, S_curr,
            mcmc, posterior
        )
        acc_count += was_acc

        if loading_bar
            next!(iter)
        end

    end
    mcmc.γ_counts[1:2] = [acc_count, i]
    return
end

function draw_sample_gamma(
    mcmc::IexMcmcSampler,
    posterior::T,
    S_fixed::InteractionSequence{Int};
    desired_samples::Int=mcmc.desired_samples,
    args...
    ) where {T<:Union{SisPosterior,SimPosterior}}

    sample_out = zeros(desired_samples)
    draw_sample_gamma!(
        sample_out,
        mcmc, posterior,
        S_fixed;
        args...
        )
    return sample_out

end

function (mcmc::IexMcmcSampler)(
    posterior::T,
    S_fixed::InteractionSequence{Int};
    args...
    ) where {T<:Union{SisPosterior,SimPosterior}}

    γ_sample = draw_sample_gamma(
        mcmc, posterior, S_fixed;
        args...
    )
    S_sample = [deepcopy(S_fixed) for i in eachindex(γ_sample)]
    suff_stat = sum(posterior.dist(x,S_fixed) for x in posterior.data)
    suff_stat_trace = fill(suff_stat, length(γ_sample))
    return PosteriorMcmcOutput(
        S_sample, γ_sample,
        suff_stat_trace,
        posterior
    )
end

function draw_sample!(
    sample_out_S::T,
    sample_out_gamma::S,
    mcmc::IexMcmcSampler,
    posterior::V;
    burn_in::Int=mcmc.burn_in,
    lag::Int=mcmc.lag,
    S_init::Vector{Path{Int}}=sample_frechet_mean(posterior.data, posterior.dist)[1],
    γ_init::Float64=5.0,
    loading_bar::Bool=true
    ) where {T<:Union{InteractionSequenceSample{Int},SubArray},S<:Union{Vector{Float64},SubArray},V<:Union{SisPosterior,SimPosterior}}

    if loading_bar
        iter = Progress(
            length(sample_out_S)*lag + burn_in, # How many iters
            1,  # At which granularity to update loading bar
            "Chain for n = $(posterior.sample_size) (joint)....")  # Loading bar. Minimum update interval: 1 second
    end

    pointers = mcmc.pointers
    S_curr, S_prop = intialise_states!(pointers, S_init)
    γ_curr = γ_init

    sample_count = 1    # Keeps which sample to be stored we are working to get
    i = 0               # Keeps track all samples (included lags and burn_ins)
    γ_acc_count = 0
    reset_counts!(mcmc.move) # Reset counts for mcmc move (tracks acceptances)

    # Initialise aux data
    initialise!(mcmc.aux, S_curr, γ_curr, posterior.S_prior, posterior.sample_size)

    # Initialise sufficient statistic
    dist, data = (posterior.dist, posterior.data)
    suff_stats = mcmc.suff_stats
    suff_stats[1] = sum(dist(S_curr, x) for x in data)

    # Storage for sample statistic
    sample_suff_stat = zeros(length(sample_out_S))

    while sample_count ≤ length(sample_out_S)
        i += 1
        # Store value
        if (i > burn_in) & (((i-1) % lag)==0)
            @inbounds sample_out_S[sample_count] = deepcopy(S_curr)
            @inbounds sample_out_gamma[sample_count] = γ_curr
            @inbounds sample_suff_stat[sample_count] = suff_stats[1]
            sample_count += 1
        end
        # Accept-reject move
        accept_reject_mode!(
            S_curr, S_prop,
            γ_curr, pointers,
            mcmc, posterior
        )
        # Accept-reject dispersion
        γ_curr, was_acc = accept_reject_gamma(
            γ_curr, S_curr,
            mcmc,
            posterior
        )
        # println(S_curr)
        γ_acc_count += was_acc
        # println(S_prop)
        if loading_bar
            next!(iter)
        end
    end
    mcmc.γ_counts[1:2] = [γ_acc_count, i]
    return_states!(S_curr, S_prop, pointers)
    return sample_suff_stat
end

function draw_sample(
    mcmc::T,
    posterior::S;
    desired_samples::Int=mcmc.desired_samples,
    args...
    ) where {T<:IexMcmcSampler,S<:Union{SisPosterior,SimPosterior}}

    sample_out_S = InteractionSequenceSample{Int}(undef, desired_samples)
    sample_out_gamma = zeros(desired_samples)
    suff_stat_sample = draw_sample!(
        sample_out_S,
        sample_out_gamma,
        mcmc,
        posterior;
        args...
        )
    return sample_out_S, sample_out_gamma, suff_stat_sample
end

function (mcmc::T where {T<:IexMcmcSampler})(
    posterior::S;
    args...
    ) where {S<:Union{SisPosterior,SimPosterior}}
    S_sample, γ_sample, suff_stat_trace = draw_sample(
        mcmc, posterior;
        args...
    )
    return PosteriorMcmcOutput(
        S_sample, γ_sample,
        suff_stat_trace,
        posterior
    )
end


text:
#=
    The file contains the implementation of the FGLM algorithm,
    an algorithm for ordering conversion
=#

# does not insert duplicates
function insertsorted!(array, item, rev=false)
    pos = searchsortedfirst(array, item, rev=rev)
    if !isempty(array) && pos <= length(array) && array[pos] == item
        return array
    end
    insert!(array, pos, item)
    return array
end

# TODO: revise this criterion
function insert_nexts!(nextmonomials, monom)
    xs = gens(parent(monom))
    for new_monom in monom .* xs
        insertsorted!(nextmonomials, new_monom)
    end
end

"""
    Represents the given polynomial as an element from Kⁿ
    vector space with respect to monomials occuring in monoms
"""
function decompose_poly(f, monoms)
    v = zeros(base_ring(f), length(monoms))
    for (j, c, m) in zip(1:length(f), coefficients(f), monomials(f))
        ind = findfirst(c -> c == m, monoms)
        v[ind] = c
    end
    v
end

"""
    Represents the given vector from Kⁿ as a polynomial
    in the given monomials
"""
function compose_poly(v, monoms)
    f = zero(parent(first(monoms)))
    for (i, m) in enumerate(monoms)
        f += m * v[i]
    end
    f
end

struct RrefBasis{C}
    rref_vectors::Vector{Vector{C}} # vectors stored in rref form
    monoms::Vector{MPoly{C}}   # monomials corresponding to columns of rref
    new_monoms::Vector{MPoly{C}} # monomials from the new basis
end

RrefBasis{C}() where {C} = RrefBasis(Vector{Vector{C}}(), Vector{MPoly{C}}(), Vector{MPoly{C}}())

function update_basis!(basis::RrefBasis, newpoly)
    ground = base_ring(newpoly)

    n = length(basis.monoms)
    for new_monom in monomials(newpoly)
        (new_monom in basis.monoms) && continue
        push!(basis.monoms, new_monom)
    end

    for vec in basis.rref_vectors
        append!(vec, zeros(ground, length(basis.monoms) - n))
    end
end

function solve_in_rref(rref_vectors, vector)
    coeffs = zeros(parent(first(vector)), length(rref_vectors))
    for (i, basis_vector) in enumerate(rref_vectors)
        pivot = findfirst(!iszero, basis_vector)
        iszero(vector[pivot]) && continue
        coeffs[i] = vector[pivot]
        vector .-= coeffs[i] .* basis_vector
    end
    exists = iszero(vector)
    coeffs, exists, vector
end

function linear_relation!(basis::RrefBasis, poly)
    vector = decompose_poly(poly, basis.monoms)
    coeffs, exists, vector = solve_in_rref(basis.rref_vectors, vector)
    coeffs, exists, vector
end

function update_rref!(basis::RrefBasis, new_monom, new_vector, coeffs)
    push!(basis.new_monoms, new_monom)
    pivot = findfirst(!iszero, new_vector)

    n_coeff = new_vector[pivot]
    new_vector .//= n_coeff

    for (i, c) in enumerate(coeffs)
        basis.new_monoms[end] -= c * basis.new_monoms[i]
    end
    basis.new_monoms[end] *= inv(n_coeff)

    es = basis.rref_vectors
    for i in 1:length(es)
        iszero(es[i][pivot]) && continue
        basis.new_monoms[i] -= es[i][pivot] * basis.new_monoms[end]
        es[i] .-= es[i][pivot] .* new_vector
    end
    push!(basis.rref_vectors, new_vector)

end

function linear_combination(basis::RrefBasis, coeffs)
    sum(c*x  for (c, x) in zip(coeffs, basis.new_monoms))
end

"""
    Converts the basis `G` into the :lex ordering
"""
function fglm(ring, gbexps, gbcoeffs, rng)

    # TODO #
    RR, _ = PolynomialRing(GF(Int(ring.ch)), ["x$i" for i in 1:length(gbcoeffs)], ordering=:degrevlex)
    println(typeof(RR))
    G = export_basis(RR, gbexps, gbcoeffs)
    ###

    origR  = parent(first(G))
    origxs = gens(origR)
    ground = base_ring(origR)
    T      = elem_type(origR)

    @assert ordering(origR) == :degrevlex
    @assert characteristic(ground) > 0

    # @info "converting from $(ordering(origR)) to :lex with fglm"

    newR, _ = PolynomialRing(ground, string.(origxs), ordering=:lex)

    newG      = T[]            #  lex
    staircase = T[]            #  orderless
    nextmonomials = T[newR(1)] # lex

    i = 0
    basis = RrefBasis{elem_type(ground)}()

    while !isempty(nextmonomials)
        i += 1
        if i > 10000
            @error "Something probably went wrong in FGLM"
            break
        end
        if i % 10 == 0
            # @debug "fglm iter $i" size(newG)
        end

        m_lex = popfirst!(nextmonomials)
        m_deg = change_base_ring(ground, m_lex, parent=origR)

        # do we even need this ?
        if any(x -> is_term_divisible(m_lex, x), staircase)
            continue
        end

        # we convert to origR to ensure the output lives there
        nf = normal_form(m_deg, G)

        # pad vectors with zeros and add new monoms
        update_basis!(basis, nf)

        # try to find linear relation between vectors of basis
        # resulting to nf
        λ, exists, reduced = linear_relation!(basis, nf)

        if exists
            poly = m_lex - linear_combination(basis, λ)
            push!(newG, poly)
            push!(staircase, m_lex)
        else
            update_rref!(basis, m_lex, reduced, λ)
            insert_nexts!(nextmonomials, m_lex)
        end
    end

    _, newexps, newcoeffs = convert_to_internal(newG)

    newexps, newcoeffs
end

ext: jl


text:
function file_to_array(from::IO, t::DataType)
	seekend(from)
	eof = position(from)
	seekstart(from)
	buffer = Array{t, 1}(undef, fld(eof, sizeof(t)))
	read!(from, buffer)
	return buffer
end

function is_file_sorted(from::IO, t::DataType)
	seekend(from)
	eof = position(from)
	seekstart(from)
	buffer = Array{t, 1}(undef, 1048576)
	sorted = true
	while eof - position(from) >= (1048576 * sizeof(t))
		read!(from, buffer)
		sorted &= is_sorted(buffer)
	end
	if eof - position(from) > 0
		buffer2 = Array{t, 1}(undef, fld(eof - position(from), sizeof(t)))
		read!(from, buffer2)
		sorted &= is_sorted(buffer2)
	end
	return sorted
end

function generate_random_file(to::IO, t::DataType, size::Int)
	seekstart(to)
	while size - position(to) >= 1048576
		write(to, rand(t, 1048576))
	end
	if size - position(to) > 0
		write(to, rand(t, size - position(to)))
	end
end

function generate_random_file(to::IO, t::DataType, size::Int, range)
	seekstart(to)
	while size - position(to) >= 1048576
		buffer::Array{t, 1} = rand(range, 1048576)
		write(to, buffer)
	end
	if size - position(to) > 0
		buffer2::Array{t, 1} = rand(range, size - position(to))
		write(to, buffer2)
	end
end

function merge_into_file(from::IO, to::IO, data::Array, t::DataType)
	seekend(from)
	eof::Int = position(from)
	seekstart(from)
	inbuffer::Array{t, 1} = Array{t, 1}(undef, 1048576)
	outbuffer::Array{t, 1} = Array{t, 1}(undef, 1048576)

	filebuffer_max = fld(eof, sizeof(t)) > 1048576 ? 1048576 : fld(eof, sizeof(t))

	ptrdata::Int32 = 1
	ptrfile::Int32 = 1048576
	optr::Int32 = 1

	moreindata = true
	moreinfile = true

	while moreindata || moreinfile
		if moreinfile && (ptrfile > filebuffer_max)
			if eof - position(from) >= 1048576 * sizeof(t)
				read!(from, inbuffer)
				ptrfile = 1
			elseif eof - position(from) > 0
				filebuffer_max = fld(eof - position(from), 4)
				temp_buffer = Array{t, 1}(undef, fld(eof - position(from), 4))
				read!(from, temp_buffer)
				for i in 1:filebuffer_max
					inbuffer[i] = temp_buffer[i]
				end
				ptrfile = 1
			else
				moreinfile = false
			end
		elseif ptrfile > filebuffer_max
			moreinfile = false
		end

		if ptrdata > length(data)
			moreindata = false
		end

		if moreindata && moreinfile
			if data[ptrdata] <= inbuffer[ptrfile]
				outbuffer[optr] = data[ptrdata]
				ptrdata += 1
				optr += 1
			else
				outbuffer[optr] = inbuffer[ptrfile]
				ptrfile += 1
				optr += 1
			end
		elseif moreinfile
			outbuffer[optr] = inbuffer[ptrfile]
			ptrfile += 1
			optr += 1
		end

		if optr >= 1048576
			write(to, outbuffer)
			optr = 1
		end
	end
	if optr != 1
		temp_buffer = Array{t, 1}(undef, optr - 1)
		for i in 1:(optr - 1)
			temp_buffer[i] = outbuffer[i]
		end
		write(to, temp_buffer)
	end
end

# function simple_filesort(from::IO, to::IO, t::DataType)

# end

# function filesort(from::IO, to::IO, t::DataType, batch_size=1048576)
# 	seekend(from)
# 	eof = position(from)
# 	seekstart(from)

# 	temp_path = ""

# 	buffer = Array{t, 1}(undef, batch_size)
# 	while eof - position(from) >= (batch_size * sizeof(t))
# 		read!(from, buffer)
# 		quicksort!(buffer)

# 	end
# end


text:
using Documenter
using Plots
using CapacityExpansion

makedocs(sitename="CapacityExpansion.jl",
    authors = "Elias Kuepper, and Holger Teichgraeber",
    pages = [
        "Introduction" => "index.md",
        "Workflow" => "workflow.md",
        "Data" => ["preparing_clust_data.md", "preparing_opt_cep_data.md", "csv_structure.md", "README_GER_18.md", "README_GER_1.md", "README_CA_14.md", "README_CA_1.md", "README_TX_1.md"],
        "Optimization" => ["opt_cep.md","opt_cep_examples.md"],
        ],
    format = Documenter.HTML(assets = [
        "assets/clust_for_opt_text.svg",
        "assets/opt_cep.svg",
        "assets/workflow.svg",
        "assets/GER_1.svg",
        "assets/GER_18.svg",
        "assets/CA_1.svg",
        "assets/CA_14.svg",
        "assets/opt_cep_cap_plot",
        "assets/preparing_clust_data_load",
        "assets/preparing_clust_data_agg"])
        )

deploydocs(repo = "github.com/YoungFaithful/CapacityExpansion.jl.git", devbranch="dev")


text:
struct TeXParseError <: Exception
    msg::String
    stack::Stack
    position::Int
    data
end

function Base.showerror(io::IO, e::TeXParseError)
    println(io, "TeXParseError: ",  e.msg)
    show_state(io, e.stack, e.position, e.data)
end

function show_state(io::IO, stack, position, data)
    if position > lastindex(data)
        println(io, "after the end of the data")
        println(io, data)
    else
        # Compute the index of the char from the byte index
        position_to_id = zeros(Int, lastindex(data))
        id = 0
        for p in 1:lastindex(data)
            if isvalid(data, p)
                id += 1
            end
            position_to_id[p] = id
        end
        p = position_to_id[position]
        println(io, "at position ", position, " (string index ", p, ")")
        println(io, data)
        println(io, " " ^ (p - 1), "^")
    end

    println(io, "with stack (length $(length(stack))):")
    for (k, level) in enumerate(stack)
        println(io, "[$k] ", level)
    end
end

show_state(stack, position, data) = show_state(stdout, stack, position, data)

# Super and subscript
super = re"\^"
super.actions[:exit] = [:end_command_builder, :setup_decorated, :begin_super]

sub = re"_"
sub.actions[:exit] = [:end_command_builder, :setup_decorated, :begin_sub]

# Groups
lbrace = re"{"
lbrace.actions[:exit] = [:end_command_builder, :begin_group]

rbrace = re"}"
rbrace.actions[:exit] = [:end_command_builder, :end_group, :end_token]

# Commands
bslash = re"\\"
bslash.actions[:exit] = [:end_command_builder, :begin_command_builder]

command_char = re"[A-Za-z]"
command_char.actions[:exit] = [:push_char, :end_token]

# Characters
space = re" "
special_char = lbrace | rbrace | bslash | super | sub | command_char | space
other_char = re"." \ special_char
other_char.actions[:enter] = [:end_command_builder]
other_char.actions[:exit] = [:push_char, :end_token]

mathexpr = re.rep(special_char | other_char)
mathexpr.actions[:all] = [:show_debug_info]
mathexpr.actions[:exit] = [:end_command_builder]

machine = Automa.compile(mathexpr)

current(stack) = first(stack)
current_head(stack) = head(current(stack))
push_to_current!(stack, arg) = push!(current(stack).args, arg)

number_of_arguments = Dict(
    :frac => 2,
    :sqrt => 1,
    :accent => 1,
    :mathfont => 1
)

function has_all_arguments(texexpr)
    required = get(number_of_arguments, head(texexpr), 0)
    return required <= length(texexpr.args)
end

const require_token = [
    :subscript,
    :superscript,
    :left_delimiter,
    :right_delimiter
]

function requirement(stack)
    !has_all_arguments(current(stack)) && return :argument
    current_head(stack) in require_token && return :token
    return :none
end

function end_token!(stack)
    requirement(stack) != :token && return
    token = pop!(stack)

    if current_head(stack) in [:decorated, :underover, :integral]
        decorated = pop!(stack)
        id = head(token) == :subscript ? 2 : 3

        !isnothing(decorated.args[id]) && throw(
            TeXParseError("multiple subscripts or superscripts", stack, p, data))
        decorated.args[id] = first(token.args)
        push_to_current!(stack, decorated)
    elseif head(token) == :left_delimiter
        push_to_current!(stack, token)
    elseif head(token) == :right_delimiter
        push_to_current!(stack, token)
        delimited = pop!(stack)

        # Simplify by constructing a single :delimited expr
        left = delimited.args[1]
        right = delimited.args[end]
        args = delimited.args[2:end-1]
        content = length(args) == 1 ? first(args) : TeXExpr(:group, args)
        simplified = TeXExpr(:delimited, [left.args[1], content, right.args[1]])
        push_to_current!(stack, simplified)
    end
end

actions = Dict(
    :show_debug_info => quote
        if showdebug
            @info "New transition"
            show_state(stack, p, data)
        end
    end,
    :begin_group => quote
        push!(stack, TeXExpr(:group))
    end,
    :end_group => quote
        current_head(stack) != :group && throw(
            TeXParseError("Unexpected '}' at position $(p-1)", stack, p, data))
        group = pop!(stack)

        # Remove nestedness for group with a single element
        if length(group.args) == 1
            group = first(group.args)
        end

        if requirement(stack) == :argument
            push_to_current!(stack, group)

            if has_all_arguments(current(stack))
                command = pop!(stack)
                push_to_current!(stack, command)
            end
        else
            push_to_current!(stack, group)
        end
    end,
    :push_char => quote
        if isvalid(data, p-1)
            char = data[prevind(data, p)]
            push_to_current!(stack, get_symbol_expr(char))
        end
    end,
    :end_token => quote
        if isvalid(data, p-1)
            end_token!(stack)
        end
    end,
    :begin_sub => quote
        push!(stack, TeXExpr(:subscript))
    end,
    :begin_super => quote
        push!(stack, TeXExpr(:superscript))
    end,
    :setup_decorated => quote
        core = pop!(current(stack).args)
        if head(core) ∉ [:decorated, :integral, :underover]
            push!(stack, TeXExpr(:decorated, Any[core, nothing, nothing]))
        else
            push!(stack, core)
        end
    end,
    :begin_command_builder => quote
        push!(stack, TeXExpr(:command_builder))
    end,
    :end_command_builder => quote
        if current_head(stack) == :command_builder
            command_builder = pop!(stack)
            command_name = String(Char.(command_builder.args))

            if command_name == "frac"
                push!(stack, TeXExpr(:frac))
            elseif command_name == "sqrt"
                push!(stack, TeXExpr(:sqrt))
            elseif command_name == "left"
                push!(stack, TeXExpr(:delimited))
                push!(stack, TeXExpr(:left_delimiter))
            elseif command_name == "right"
                current_head(stack) != :delimited && throw(
                    TeXParseError("unexpected '\\right' at position $(p-1)",
                    stack, p, data))
                push!(stack, TeXExpr(:right_delimiter))
            elseif is_supported_command(command_name)
                push_to_current!(stack, get_command_expr(command_name))
                end_token!(stack)
            else
                throw(
                    TeXParseError("unsupported command \\$command_name",
                    stack, p, data))
            end
        end
    end
)

context = Automa.CodeGenContext()
@eval function texparse(data ; showdebug=false)
    $(Automa.generate_init_code(context, machine))
    p_end = p_eof = lastindex(data)

    # Needed to avoid problem with multi bytes unicode chars
    stack = Stack{Any}()
    push!(stack, TeXExpr(:expr))

    try
        $(Automa.generate_exec_code(context, machine, actions))
    catch
        throw(TeXParseError("unexpected error while parsing", stack, p, data))
    end

    if length(stack) > 1
        err = TeXParseError(
            "end of string reached with unfinished $(current(stack).head)",
            stack, p_eof, data)
        throw(err)
    end

    return current(stack)
end


text:
using FileIO
using DataFrames
##
function failure_idx(count::T, rewarded) where {T <: Number}
    if rewarded
        count = 0.0
    else
        count + 1.0
    end
end
##
filename = "/home/beatriz/mainen.flipping.5ht@gmail.com/Flipping/Datasets/Stimulations/DRN_Opto_Flipping/pokesDRN_Opto_Flipping.csv"
data = FileIO.load(filename) |> DataFrame
data[!,:LastPoke] .= false
data[!,:FailuresIdx] .= 0.0
by(data,[:Session,:Streak]) do dd
    dd.LastPoke[end] = true
    dd[:,:FailuresIdx] = Base.accumulate(failure_idx, occursin.("true",dd.Reward);init=0.0)
end
checkpoint = by(data,:Session) do d
    d.PokeIn[1] == d.PokeIn[2]
end
any(checkpoint[:,2])
println(describe(data))
filtered = filter(t -> t[:Protocol].!="45/15", data)
##
"""
SCATTER OF MEAN ANS SD PER PROTOCOL AND BARRIER
"""
##
function scatter_mu_sd(df,var)
    μ = mean(select(df,var))
    σ = stderr(select(df,var))
    dict = OrderedDict(μ => σ)
end


text:
using MacroTools
import ..LambdaCalculus.LocallyNameless

export @lambda, @λ

function meta_convert1(expr, Γ::NamingContext)
    if isexpr(expr, :$)
        @assert length(expr.args) == 1
        return esc(expr.args[1])
    end

    @match expr begin
        (f_(arg1_, args__)) => mapfoldl(e -> meta_convert1(e, Γ),
                                        (e, arg) -> :(App($e, $arg)),
                                        [f, arg1, args...])
        (arg_Symbol -> body_) => :(Abs($(meta_convert1(body, pushfirst(Γ, arg)))))
        (name_Symbol) => begin
            index = findfirst(isequal(name), Γ)
            if index === nothing
                :(FVar($(Meta.quot(name))))
            else
                :(BVar($index))
            end
        end
        other_ => error("unhandled expression: $other")
    end
end

meta_convert(expr, Γ::NamingContext) = meta_convert1(MacroTools.prewalk(unblock ∘ rmlines, expr), Γ)

"Convert a (well-formed) Julia expression to a `Term`."
macro lambda(expr)
    meta_convert(expr, NamingContext())
end

"Convert a (well-formed) Julia expression to a `Term`."
macro λ(expr)
    meta_convert(expr, NamingContext())
end


text:
abstract type AbstractVecArray{T<:AbstractVec,N} <: AbstractArray{T,N} end

@inline Base.size(u::AbstractVecArray) =
    size(u.x)

@inline Base.length(u::AbstractVecArray) =
    length(u.x)

@inline Base.linearindices(u::AbstractVecArray) =
    linearindices(u.x)

Base.IndexStyle(::Type{<:AbstractVecArray}) =
    IndexLinear()

struct VecArray{T<:Number,N} <: AbstractVecArray{Vec{T},N}
    x::Array{T,N}
    y::Array{T,N}
    z::Array{T,N}

    function VecArray{T,N}(x,y,z) where {T,N}
        (size(x) == size(y) == size(z)) || throw(ArgumentError("Arrays must have the same dimensions"))
        return new{T,N}(x,y,z)
    end

end

VecArray(x,y,z) =
    VecArray{eltype(x),ndims(x)}(x,y,z)

@inline function Base.getindex(u::VecArray{T,N},I) where {T,N}
    x = u.x
    y = u.y
    z = u.z
    @boundscheck checkbounds(x,I)
    @inbounds begin
        xv = x[I]
        yv = y[I]
        zv = z[I]
    end
    return Vec{T}(xv,yv,zv)
end

@inline function Base.setindex!(u::VecArray{T,N},v::Vec,I) where {T,N}
    x = u.x
    y = u.y
    z = u.z
    @boundscheck checkbounds(x,I)
    @inbounds begin
        x[I] = v.x
        y[I] = v.y
        z[I] = v.z
    end
    return u
end

struct Vec2DArray{T<:Number,N} <: AbstractVecArray{Vec2D{T},N}
    x::Array{T,N}
    y::Array{T,N}

    function Vec2DArray{T,N}(x,y) where {T,N}
        (size(x) == size(y)) || throw(ArgumentError("Arrays must have the same dimensions"))
        return new{T,N}(x,y)
    end

end

Vec2DArray(x::AbstractArray,y) =
    Vec2DArray{eltype(x),ndims(x)}(x,y)

Vec2DArray{T}(L::Integer) where T = Vec2DArray{T,1}(zeros(T,L),zeros(T,L))

Base.similar(a::Vec2DArray{T}) where {T} = Vec2DArray{T}(length(a))

@inline function Base.getindex(u::Vec2DArray{T,N},I) where {T,N}
    x = u.x
    y = u.y
    @boundscheck checkbounds(x,I)
    @inbounds begin
        xv = x[I]
        yv = y[I]
    end
    return Vec2D{T}(xv,yv)
end

@inline function Base.setindex!(u::Vec2DArray{T,N},v::Vec2D,I) where {T,N}
    x = u.x
    y = u.y
    @boundscheck checkbounds(x,I)
    @inbounds begin
        x[I] = v.x
        y[I] = v.y
    end
    return u
end

struct Vec1DArray{T<:Number,N} <: AbstractVecArray{Vec1D{T},N}
    x::Array{T,N}

    function Vec1DArray{T,N}(x) where {T,N}
        return new{T,N}(x)
    end

end

Vec1DArray(x::AbstractArray) =
    Vec1DArray{eltype(x),ndims(x)}(x)

@inline function Base.getindex(u::Vec1DArray{T,N},I) where {T,N}
    x = u.x
    @boundscheck checkbounds(x,I)
    @inbounds begin
        xv = x[I]
    end
    return Vec1D{T}(xv)
end

@inline function Base.setindex!(u::Vec1DArray{T,N},v::Vec1D,I) where {T,N}
    x = u.x
    @boundscheck checkbounds(x,I)
    @inbounds begin
        x[I] = v.x
    end
    return u
end


text:
### Utility Functions needed for the construction of the Power System, mostly used for consistency checking ####

## Time Series Length ##

function timeseriescheckload(loads::Array{T}) where {T<:ElectricLoad}
    t = length(loads[1].scalingfactor)
    for l in loads
        if t == length(l.scalingfactor)
            continue
        else
            @error "Inconsistent load scaling factor time series length"
        end
    end
    return t
end

function timeserieschecksources(generators::Array{T}, expected) where {T<:Generator}
    for g in generators
        actual = length(g.scalingfactor)
        if expected == actual
            continue
        else
            @error "Inconsistent generation scaling factor time series length for $(g.name)" expected actual
        end
    end
end

## Check that all the buses have a type defintion ##

function buscheck(buses::Array{Bus})
    for b in buses
        if b.bustype == nothing
            @warn "Bus/Nodes data does not contain information to build an a network"
        end
    end

    check_ascending_order([x.number for x in buses], "Bus")
end

## Slack Bus Definition ##

function slackbuscheck(buses::Array{Bus})
    slack = -9
    for b in buses
        if b.bustype == "SF"
            slack = b.number
        end
    end
    if slack == -9
        @error "Model doesn't contain a slack bus"
    end
end

### PV Bus Check ###

function pvbuscheck(buses::Array{Bus}, generators::Array{T}) where {T<:Generator}
    pv_list = -1*ones(Int64, length(generators))
    for (ix,g) in enumerate(generators)
        g.bus.bustype == "PV" ? pv_list[ix] = g.bus.number : continue
    end

    for b in buses
        if b.bustype == "PV"
            b.number in pv_list ? continue : 0 #@warn "The bus ", b.number, " is declared as PV without a generator connected to it"
        else
            continue
        end
    end
end

# TODO: Check for islanded Buses

# check for minimum timediff
function minimumtimestep(loads::Array{T})where {T<:ElectricLoad}
    if length(loads[1].scalingfactor) > 1
        timeseries = loads[1].scalingfactor
        n = length(timeseries)-1
        ts = []
        for i in 1:n
            push!(ts,TimeSeries.timestamp(timeseries)[n+1]-TimeSeries.timestamp(timeseries)[n])
        end
        return minimum(ts)
    else
        ts =Dates.Dates.Minute(1)
        return ts
    end
end

# convert generator ramp rates to a consistent denominator
function convertramp(ramplimits::Union{NamedTuple{(:up, :down),Tuple{Float64,Float64}},Nothing}, ts::Dates.TimePeriod)
    if isa(ramplimits,NamedTuple)
        hr = typeof(ts)(Dates.Dates.Minute(1))
        scaling  = hr/ts
        up = ramplimits.up/scaling
        down = ramplimits.down/scaling
        R = (up = up, down = down)
        return R
    else
        return ramplimits
    end
end

# check for valid ramp limits
function checkramp(generators::Array{T}, ts::Dates.TimePeriod) where {T<:Generator}
    for (ix,g) in enumerate(generators)
        if isa(g,ThermalDispatch)
            R = convertramp(g.tech.ramplimits,ts)
            generators[ix] = ThermalDispatch(deepcopy(g.name),deepcopy(g.available),deepcopy(g.bus),
                                            TechThermal(deepcopy(g.tech.activepower),deepcopy(g.tech.activepowerlimits),
                                                        deepcopy(g.tech.reactivepower),deepcopy(g.tech.reactivepowerlimits),
                                                        R,deepcopy(g.tech.timelimits)),
                                            deepcopy(g.econ)
                                            )
            if isa(g.tech.ramplimits, NamedTuple)
                if g.tech.ramplimits.up >= (g.tech.activepowerlimits.max - g.tech.activepowerlimits.min)
                    @warn "The generator $(g.name) has a nonbinding ramp up limit."
                end
                if g.tech.ramplimits.down >= (g.tech.activepowerlimits.max - g.tech.activepowerlimits.min)
                    @warn "The generator $(g.name) has a nonbinding ramp down limit."
                end
            else
                @info "Ramp defined as nothing for $(g.name)"
            end
        elseif isa(g,ThermalGenSeason)
            R = convertramp(g.tech.ramplimits,ts)
            generators[ix] = ThermalGenSeason(deepcopy(g.name),deepcopy(g.available),deepcopy(g.bus),
                                            TechThermal(deepcopy(g.tech.activepower),deepcopy(g.tech.activepowerlimits),
                                                        deepcopy(g.tech.reactivepower),deepcopy(g.tech.reactivepowerlimits),
                                                        R,deepcopy(g.tech.timelimits)),
                                            deepcopy(g.econ),
                                            deepcopy(g.scalingfactor)
                                            )
            if isa(g.tech.ramplimits, NamedTuple)
                if g.tech.ramplimits.up >= (g.tech.activepowerlimits.max - g.tech.activepowerlimits.min)
                    @warn "The generator $(g.name) has a nonbinding ramp up limit."
                end
                if g.tech.ramplimits.down >= (g.tech.activepowerlimits.max - g.tech.activepowerlimits.min)
                    @warn "The generator $(g.name) has a nonbinding ramp down limit."
                end
            else
                @info "Ramp defined as nothing for $(g.name)"
            end
        elseif isa(g,HydroCurtailment)
            R = convertramp(g.tech.ramplimits,ts)
            generators[ix] = HydroCurtailment(deepcopy(g.name),deepcopy(g.available),deepcopy(g.bus),
                                            TechHydro(deepcopy(g.tech.installedcapacity),deepcopy(g.tech.activepower),deepcopy(g.tech.activepowerlimits),
                                                    deepcopy(g.tech.reactivepower),deepcopy(g.tech.reactivepowerlimits),
                                                    R,deepcopy(g.tech.timelimits)),
                                            deepcopy(g.econ.curtailpenalty),
                                            deepcopy(g.scalingfactor)
                                            )
            if isa(g.tech.ramplimits, NamedTuple)
                if g.tech.ramplimits.up >= (g.tech.activepowerlimits.max - g.tech.activepowerlimits.min)
                    @info "The generator $(g.name) has a nonbinding ramp up limit." maxlog=PS_MAX_LOG
                end
                if g.tech.ramplimits.down >= (g.tech.activepowerlimits.max - g.tech.activepowerlimits.min)
                    @info "The generator $(g.name) has a nonbinding ramp down limit." maxlog=PS_MAX_LOG
                end
            else
                @info "Ramp defined as nothing for $(g.name)"
            end
        elseif isa(g,HydroStorage)
            R = convertramp(g.tech.ramplimits,ts)
            generators[ix] = HydroStorage(deepcopy(g.name),deepcopy(g.available),deepcopy(g.bus),
                                            TechHydro(deepcopy(g.tech.installedcapacity),deepcopy(g.tech.activepower),deepcopy(g.tech.activepowerlimits),
                                                    deepcopy(g.tech.reactivepower),deepcopy(g.tech.reactivepowerlimits),
                                                    R,deepcopy(g.tech.timelimits)),
                                            deepcopy(g.econ),
                                            deepcopy(g.storagecapacity),
                                            deepcopy(g.scalingfactor)
                                            )
            if isa(g.tech.ramplimits, NamedTuple)
                if g.tech.ramplimits.up >= (g.tech.activepowerlimits.max - g.tech.activepowerlimits.min)
                    @info "The generator $(g.name) has a nonbinding ramp up limit." maxlog=PS_MAX_LOG
                end
                if g.tech.ramplimits.down >= (g.tech.activepowerlimits.max - g.tech.activepowerlimits.min)
                    @info "The generator $(g.name) has a nonbinding ramp down limit." maxlog=PS_MAX_LOG
                end
            else
                @info "Ramp defined as nothing for $(g.name)"
            end
        end
    end
    return generators
end


text:
function toy_rnn()
    batch_size  = 4
    input_size  = 1
    output_size = 1
    hidden_size = 10

    function toy_data(n)
        x = rand(Float32,  input_size,  n*batch_size)
        y = zeros(Float32, output_size, n*batch_size)
        for i in 1:n
            for t in 1:batch_size
                y[i*batch_size] += x[1,(i-1)*batch_size+t]
            end
        end
        (x, y)
    end

    function toy_net(xtrn, ytrn, xtst, ytst, iter)
        net = NNN.Network(batch_size)
        NNN.push!(net,
                  NNN.SimpleRecurrentLayer(batch_size, input_size, hidden_size;
                                           activation=NNN.LeakyRelu()),
                  NNN.SGD(;ε=0.05),
                  NNN.SGD(;ε=0.05))
        NNN.push!(net,
                  NNN.OutputLayer(batch_size, hidden_size, output_size;
                                  activation=NNN.IdentityAct(),
                                  loss=NNN.rnn_quadratic_cost),
                  NNN.SGD(;ε=0.05))

        NNN.train(xtrn, ytrn, net; iter=iter, batch_size=batch_size)
        errors = NNN.test(xtst, ytst, net; batch_size=batch_size)

        (net, errors)
    end

    xtrn, ytrn  = toy_data(1000);
    xtst, ytst  = toy_data(100);
    net, errors = toy_net(xtrn, ytrn, xtst, ytst, 5)
    (net, errors, xtst, ytst)
end

net, errors, xtst, ytst = toy_rnn()

facts("Simple RNN") do
    @fact sum(abs(NNN.forward!(net, xtst[:,1:4])[4] - ytst[:,4])) => less_than(1.0)
end


text:
println("==============flag! MOI_wrapper>MOI_wrapper.jl======================================")
"""
MOI_wrapper.jl defines the Juniper.Optimizer struct
with all mandatory MOI functions overloaded
"""

"""
MOI functions, sets and, other type definitions
"""
# indices
println("==================flag!=========MOI_wrapper_indices=======================")
const VI = MOI.VariableIndex
const CI = MOI.ConstraintIndex

# sets
const BOUNDS = Union{
    MOI.EqualTo{Float64},
    MOI.GreaterThan{Float64},
    MOI.LessThan{Float64},
}

const VAR_TYPES = Union{
    MOI.ZeroOne,
    MOI.Integer
}

# other MOI types
const AFF_TERM = MOI.ScalarAffineTerm{Float64}
const QUAD_TERM = MOI.ScalarQuadraticTerm{Float64}

"""
Variable information struct definition
"""
mutable struct VariableInfo
    lower_bound::Float64  # May be -Inf even if has_lower_bound == true
    has_lower_bound::Bool # false implies lower_bound == -Inf
    upper_bound::Float64  # May be Inf even if has_upper_bound == true
    has_upper_bound::Bool # false implies upper_bound == Inf
    is_fixed::Bool        # Implies lower_bound == upper_bound and !has_lower_bound and !has_upper_bound
    is_binary::Bool       # Implies lower_bound == 0, upper_bound == 1 and is MOI.ZeroOne
    is_integer::Bool      # Implies variable is MOI.Integer
    start::Real           # Primal start
    name::String
end
VariableInfo() = VariableInfo(-Inf, false, Inf, false, false, false, false, 0.0, "")

"""
Optimizer struct
"""
mutable struct Optimizer <: MOI.AbstractOptimizer
    inner::Union{JuniperProblem, Nothing}
    variable_info::Vector{VariableInfo}
    nlp_data::MOI.NLPBlockData
    sense::MOI.OptimizationSense
    objective::Union{SVF, SAF, SQF, Nothing}
    linear_le_constraints::Vector{Tuple{SAF, MOI.LessThan{Float64}}}
    linear_ge_constraints::Vector{Tuple{SAF, MOI.GreaterThan{Float64}}}
    linear_eq_constraints::Vector{Tuple{SAF, MOI.EqualTo{Float64}}}
    quadratic_le_constraints::Vector{Tuple{SQF, MOI.LessThan{Float64}}}
    quadratic_ge_constraints::Vector{Tuple{SQF, MOI.GreaterThan{Float64}}}
    quadratic_eq_constraints::Vector{Tuple{SQF, MOI.EqualTo{Float64}}}
    options::SolverOptions
end
println("==============flag! MOI_wrapper>MOI_wrapper.jl==>Optimizer(get, supports)===================================")
MOI.get(::Optimizer, ::MOI.SolverName) = "Juniper"
MOI.supports(::Optimizer, ::MOI.Silent) = true
MOI.supports(::Optimizer, ::MOI.NumberOfThreads) = true
MOI.supports(::Optimizer, ::MOI.TimeLimitSec) = true
MOI.supports(::Optimizer, ::MOI.RawParameter) = true

function MOI.set(model::Optimizer, ::MOI.Silent, value::Bool)
    if value === true
        model.options.log_levels = []
    end
    model.options.silent = value
    return
end

function MOI.set(model::Optimizer, ::MOI.NumberOfThreads, value::Union{Nothing,Int})
    println("==============flag! MOI_wrapper>MOI_wrapper.jl==>MOI.set(::MOI.NumberOfThreads)===================================")
    if value === nothing
        model.options.processors = 1
    else
        model.options.processors = value
    end
    return
end

function MOI.set(model::Optimizer, ::MOI.TimeLimitSec, value::Union{Nothing,Float64})
    println("==============flag! MOI_wrapper>MOI_wrapper.jl==>MOI.set(::MOI.TimeLimitSec, value::Union{Nothing,Float64})===================================")
    if value === nothing
        model.options.time_limit = Inf
    else
        model.options.time_limit = value
    end
    return
end

function MOI.set(model::Optimizer, p::MOI.RawParameter, value)
    println("==============flag! MOI_wrapper>MOI_wrapper.jl==>MOI.set(p::MOI.RawParameter)===================================")
    p_symbol = Symbol(p.name)
    if in(p_symbol, fieldnames(SolverOptions))
        type_of_param = fieldtype(SolverOptions, p_symbol)
        if hasmethod(convert, (Type{type_of_param}, typeof(value)))
            passed_checks = true
            if p_symbol == :traverse_strategy
                if !(value in [:BFS, :DFS, :DBFS])
                    passed_checks = false
                    @error "Traverse strategy $(value) is not supported. Use one of `[:BFS, :DFS, :DBFS]`."
                end
            end
            if p_symbol == :branch_strategy
                if !(value in [:StrongPseudoCost, :PseudoCost, :Reliability, :MostInfeasible])
                    passed_checks = false
                    @error "Branch strategy $(value) is not supported. Use one of `[:StrongPseudoCost, :PseudoCost, :Reliability, :MostInfeasible]`."
                end
            end
            passed_checks && setfield!(model.options, p_symbol, convert(type_of_param, value))
        else
            @error "The option $(p.name) has a different type ($(type_of_param))"
        end
    else
        @error "The option $(p.name) doesn't exist."
    end
    return
end

# Returns the number of solutions that can be retrieved
println("==============flag! MOI_wrapper>MOI_wrapper.jl==>Optimizer(get, supports2)===================================")
MOI.get(model::Optimizer, ::MOI.ResultCount) = length(model.inner.solutions)

MOI.get(model::Optimizer, ::MOI.NumberOfThreads) = model.options.processors

MOI.get(model::Optimizer, ::MOI.TimeLimitSec) = model.options.time_limit

MOI.get(model::Optimizer, ::MOI.Silent) = model.options.silent

function MOI.get(model::Optimizer, p::MOI.RawParameter)
    if in(p.name, fieldnames(SolverOptions))
        return getfield(model.options, p.name)
    end
    @error "The option $(p.name) doesn't exist."
end

"""
EmptyNLPEvaluator struct and associated functions
"""
struct EmptyNLPEvaluator <: MOI.AbstractNLPEvaluator end
MOI.features_available(::EmptyNLPEvaluator) = [:Grad, :Jac, :Hess]
MOI.initialize(::EmptyNLPEvaluator, features) = nothing
MOI.eval_objective(::EmptyNLPEvaluator, x) = NaN
function MOI.eval_constraint(::EmptyNLPEvaluator, g, x)
    @assert length(g) == 0
    return
end
function MOI.eval_objective_gradient(::EmptyNLPEvaluator, g, x)
    fill!(g, 0.0)
    return
end
MOI.jacobian_structure(::EmptyNLPEvaluator) = Tuple{Int64,Int64}[]
MOI.hessian_lagrangian_structure(::EmptyNLPEvaluator) = Tuple{Int64,Int64}[]
function MOI.eval_constraint_jacobian(::EmptyNLPEvaluator, J, x)
    @assert length(J) == 0
    return
end
function MOI.eval_hessian_lagrangian(::EmptyNLPEvaluator, H, x, σ, μ)
    @assert length(H) == 0
    return
end
empty_nlp_data() = MOI.NLPBlockData([], EmptyNLPEvaluator(), false)

function register(s::Symbol, dimension::Integer, f::Function; autodiff::Bool=false)
    return RegisteredFunction(s, dimension, f, nothing, nothing, autodiff)
end

function register(s::Symbol, dimension::Integer, f::Function, gradf::Function; autodiff::Bool=false)
    return RegisteredFunction(s, dimension, f, gradf, nothing, autodiff)
end

function register(s::Symbol, dimension::Integer, f::Function, gradf::Function, grad2f::Function; autodiff::Bool=false)
    return RegisteredFunction(s, dimension, f, gradf, grad2f, autodiff)
end

"""
Optimizer struct constructor
"""
function Optimizer(;options...)

    solver_options = combine_options(options)

    return Optimizer(
    nothing,
    [],
    empty_nlp_data(),
    MOI.FEASIBILITY_SENSE,
    nothing,
    [], [], [], # linear constraints
    [], [], [], # quadratic constraints
    solver_options)
end

function Optimizer(options::Vector{Pair{String,Any}})
    symbol_options = Dict{Symbol, Any}()
    for option in options
        symbol_options[Symbol(option.first)] = option.second
    end
    return Optimizer(symbol_options)
end

function Optimizer(options::Dict{Symbol,Any})

    solver_options = combine_options(options)

    return Optimizer(
    nothing,
    [],
    empty_nlp_data(),
    MOI.FEASIBILITY_SENSE,
    nothing,
    [], [], [], # linear constraints
    [], [], [], # quadratic constraints
    solver_options)
end

"""
Printing the optimizer
"""
function Base.show(io::IO, model::Optimizer)
    println("A Juniper MathOptInterface model with backend")
    return
end

"""
Copy constructor for the optimizer
"""
MOIU.supports_default_copy_to(model::Optimizer, copy_names::Bool) = true
function MOI.copy_to(model::Optimizer, src::MOI.ModelLike; kws...)
    return MOI.Utilities.automatic_copy_to(model, src; kws...)
end

"""
``MOI.is_empty(model::Optimizer)`` overload for Alpine.Optimizer
"""
function MOI.is_empty(model::Optimizer)
    return isempty(model.variable_info) &&
        model.nlp_data.evaluator isa EmptyNLPEvaluator &&
        model.sense == MOI.FEASIBILITY_SENSE &&
        isempty(model.linear_le_constraints) &&
        isempty(model.linear_ge_constraints) &&
        isempty(model.linear_eq_constraints) &&
        isempty(model.quadratic_le_constraints) &&
        isempty(model.quadratic_ge_constraints) &&
        isempty(model.quadratic_eq_constraints)
end

"""
``MOI.empty!(model::Optimizer)`` overload for Alpine.Optimizer
"""
function MOI.empty!(model::Optimizer)
    model.inner = nothing
    empty!(model.variable_info)
    model.nlp_data = empty_nlp_data()
    model.sense = MOI.FEASIBILITY_SENSE
    model.objective = nothing
    empty!(model.linear_le_constraints)
    empty!(model.linear_ge_constraints)
    empty!(model.linear_eq_constraints)
    empty!(model.quadratic_le_constraints)
    empty!(model.quadratic_ge_constraints)
    empty!(model.quadratic_eq_constraints)
end

println("==============flag! MOI_wrapper>MOI_wrapper.jl==>ordering of constraints provided to Juniper.jl===================================")
"""
ordering of constraints provided to Juniper.jl
"""
linear_le_offset(model::Optimizer) = 0
linear_ge_offset(model::Optimizer) = length(model.linear_le_constraints)
linear_eq_offset(model::Optimizer) = linear_ge_offset(model) + length(model.linear_ge_constraints)
quadratic_le_offset(model::Optimizer) = linear_eq_offset(model) + length(model.linear_eq_constraints)
quadratic_ge_offset(model::Optimizer) = quadratic_le_offset(model) + length(model.quadratic_le_constraints)
quadratic_eq_offset(model::Optimizer) = quadratic_ge_offset(model) + length(model.quadratic_ge_constraints)
nlp_constraint_offset(model::Optimizer) = quadratic_eq_offset(model) + length(model.quadratic_eq_constraints)

function info_array_of_variables(variable_info::Vector{VariableInfo}, attr::Symbol)
    len_var_info = length(variable_info)
    type_dict = get_type_dict(variable_info[1])
    result = Array{type_dict[attr], 1}(undef, len_var_info)
    for i = 1:len_var_info
        result[i] = getfield(variable_info[i], attr)
        # if type is binary then set bounds correctly
        if result[i] < 0 && attr == :lower_bound && getfield(variable_info[i], :is_binary)
            result[i] = 0
        end
        if result[i] > 1 && attr == :upper_bound && getfield(variable_info[i], :is_binary)
            result[i] = 1
        end
    end
    return result
end

function replace_solution!(m::JuniperProblem, tree::BnBTreeObj)
    status_dict = Dict{Symbol,MOI.TerminationStatusCode}()
    status_dict[:Time] = MOI.TIME_LIMIT
    status_dict[:Infeasible] = !tree.global_solver ? MOI.LOCALLY_INFEASIBLE : MOI.INFEASIBLE
    status_dict[:MipGap] = MOI.OBJECTIVE_LIMIT
    status_dict[:BestObjStop] = MOI.OBJECTIVE_LIMIT
    status_dict[:EnoughSolutions] = MOI.SOLUTION_LIMIT

    if isdefined(tree, :incumbent)
        incumbent = tree.incumbent
        m.objval = incumbent.objval
        m.solution = incumbent.solution
    end
    if tree.limit == :None
        if tree.incumbent.only_almost
            if !tree.global_solver
                m.status = MOI.ALMOST_LOCALLY_SOLVED
            else
                m.status = MOI.ALMOST_OPTIMAL
            end
        else
            if !tree.global_solver
                m.status = MOI.LOCALLY_SOLVED
            else
                m.status = MOI.OPTIMAL
            end
        end
    else
        m.status = status_dict[tree.limit]
    end
    m.best_bound = tree.best_bound
end

"""
``MOI.optimize!()`` for Juniper
"""
function MOI.optimize!(model::Optimizer)
    println("==============flag! MOI_wrapper>MOI_wrapper.jl==>Optimizer!===================================")
    #println(model.inner)
    # feasibility_pump only if mip solver exists
    if model.options.feasibility_pump && model.options.mip_solver === nothing
        model.options.feasibility_pump = false
    end

    println("=======================0.5=====================")
    Random.seed!(1)
    MOI.initialize(model.nlp_data.evaluator, [:ExprGraph])

    if ~isa(model.nlp_data.evaluator, EmptyNLPEvaluator)

    else
        @info "no explicit NLP constraints or objective provided using @NLconstraint or @NLobjective macros"
    end

    println("=======================1=====================")
    # fill JuniperProblem
    #println(model.inner)
    model.inner = JuniperProblem()

    @views jp = model.inner

    init_juniper_problem!(jp, model)
    println("=========after init_juniper_program=========jp num_l_constr is ===========",jp.num_l_constr)
    println("=======================2=====================")
    ps = jp.options.log_levels

    jp.debugDict = Dict{Any,Any}()

    (:All in ps || :AllOptions in ps) && print_options(jp;all=true)
    (:Options in ps) && print_options(jp;all=false)

    if !jp.options.fixed_gain_mu && jp.obj_sense == :Max
        jp.options.gain_mu = 1-jp.options.gain_mu
    end
    println("=======================3=====================")
    nw = nworkers()
    if nw < jp.options.processors
        jp.options.processors = nw
        @warn "Julia was started with less processors than you defined in your options. Start julia with: `julia -p "*string(jp.options.processors)*"`"
    end
    # set incumbent to nothing might be updated using start values or the feasibility_pump
    incumbent = nothing
    println("=======================4=====================")
    create_root_model!(model, jp)
    (:All in ps || :Info in ps) && print_info(jp)
    # fix primal start to check if we have an incumbent
    println("=======================5=====================")
    fix_primal_start!(jp)

    incumbent = solve_root_incumbent_model(jp)
    println("=======================6=====================")

    unfix_primal_start!(jp)
    relax_start_time = time()
    restarts = solve_root_model!(jp)
    jp.relaxation_time = time()-relax_start_time

    (:All in ps || :Info in ps) && println("Status of relaxation: ", jp.relaxation_status)
    jp.soltime = time()-jp.start_time

    jp.options.debug && debug_fill_basic(jp.debugDict,jp,restarts)

    # if infeasible or unbounded => return
    if !state_is_optimal(jp.relaxation_status; allow_almost=jp.options.allow_almost_solved)
        jp.status = jp.relaxation_status
        if jp.options.debug && jp.options.debug_write
            write(jp.options.debug_file_path, JSON.json(jp.debugDict))
        end
        return
    end

    (:All in ps || :Info in ps || :Timing in ps) && println("Time for relaxation: ", jp.soltime)

    backend                = JuMP.backend(jp.model)
    jp.relaxation_objval   = JuMP.objective_value(jp.model)
    jp.relaxation_solution = JuMP.value.(jp.x)

    jp.options.debug && debug_objective(jp.debugDict,jp)
    # TODO free model for Knitro

    (:All in ps || :Info in ps || :Timing in ps) && println("Relaxation Obj: ", jp.relaxation_objval)

    println("=======================7=====================")
    fpump(model,jp)
    #write(jp.options.debug_file_path, JSON.json(jp.debugDict))
    #return jp
end

function Juniper_result(model::Optimizer)
    println("==============flag! MOI_wrapper>MOI_wrapper.jl==>Juniper_result!===================================")
    #println(model.inner)
    # feasibility_pump only if mip solver exists
    if model.options.feasibility_pump && model.options.mip_solver === nothing
        model.options.feasibility_pump = false
    end

    println("=======================0=====================")
    Random.seed!(1)
    MOI.initialize(model.nlp_data.evaluator, [:ExprGraph])

    if ~isa(model.nlp_data.evaluator, EmptyNLPEvaluator)

    else
        @info "no explicit NLP constraints or objective provided using @NLconstraint or @NLobjective macros"
    end

    println("=======================1=====================")
    # fill JuniperProblem
    #println(model.inner)
    model.inner = JuniperProblem()

    @views jp = model.inner

    init_juniper_problem!(jp, model)
    println("=========after init_juniper_program=========jp num_l_constr is ===========",jp.num_l_constr)
    println("=======================2=====================")
    println(jp)
    return jp
end

getnsolutions(m::JuniperProblem) = m.nsolutions
getsolutions(m::JuniperProblem) = m.solutions
getnbranches(m::JuniperProblem) = m.nbranches

include("variables.jl")
include("constraints.jl")
include("objective.jl")
include("results.jl")
include("nlp.jl")


text:
#=
Created on Thursday 23 July 2020
Last update: -

@author: Michiel Stock
michielfmstock@gmail.com

Illustration of Particle Swarm Optimization.
=#

# first we define a particle with a position `x` and a velocity `v`

using STMO.TestFuns

fun = TestFuns.ackley
x1lims = (-10, 10)
x2lims = (-10, 10)

mutable struct Particle{T}
    x::T
    v::T
    x_best::T
end

Particle(x) = Particle(x, zero(x), x)

init_population(n_particles, lims...) = [Particle([(u - l) * rand() + l for (l, u) in lims]) for i in 1:n_particles]

"""
    particle_swarm_optimization!(f, population, k_max;
            w=1, c1=1, c2=1, tracker=nothing)

Performs Particle Swarm Optimization to minimize a function `f`. Give an initial
vector of particles (type `Particle`) and the `k_max`, the number of iterations.

Optionally set hyperparameters `w`, `c1` and `c2` (default value of 1).
"""
function particle_swarm_optimization!(f, population::Vector{Particle{T}} where T, k_max;
        w=1, c1=1, c2=1, tracker=nothing)
    # find best point
    y_best, x_best = minimum((((f(part.x), part.x)) for part in population))
    for k in 1:k_max
        # update population
        for particle in population
            r1, r2 = rand(2)
            particle.v .= w * particle.v + c1 * r1 * (particle.x_best .- particle.x) .+
                c2 * r2 * (x_best .- particle.x)
            particle.x .+= particle.v
            fx = f(particle.x)
            # update current best
            fx < f(x_best) && (particle.x_best .= particle.x)
            # update global best
            if fx < y_best
                y_best = fx
                x_best .= particle.x
            end
        end
        tracker isa Nothing || tracker(population)
    end
    return y_best, x_best
end

population = init_population(50, x1lims, x2lims)
particle_swarm_optimization!(fun, population, 100, w=1)

function pso_animation(f, population, k_max;
        w=1, c1=1, c2=1)
    # find best point
    y_best, x_best = minimum((((f(part.x), part.x)) for part in population))
    objective = [y_best]
    # determine xvals and yvals depending on the spread of the particles
    x1min = minimum((part.x[1] for part in population))
    x1max = maximum((part.x[1] for part in population))
    x2min = minimum((part.x[2] for part in population))
    x2max = maximum((part.x[2] for part in population))
    x1steps = (x1max - x1min) / 200
    x2steps = (x2max - x2min) / 200
    anim = @animate for k in 1:k_max

        swarmplot = heatmap(x1min:x1steps:x1max, x2min:x2steps:x2max,
                        (x,y)->f((x,y)), color=:speed)
        xlims!((x1min, x1max))
        ylims!((x2min, x2max))
        # update population
        for particle in population
            r1, r2 = rand(2)
            particle.v .= w * particle.v + c1 * r1 * (particle.x_best .- particle.x) .+
                c2 * r2 * (x_best .- particle.x)
            particle.x .+= particle.v
            fx = f(particle.x)
            # update current best
            fx < f(x_best) && (particle.x_best .= particle.x)
            # update global best
            if fx < y_best
                y_best = fx
                x_best .= particle.x
            end
            scatter!(swarmplot, [particle.x[1]], [particle.x[2]],
                        color=:red, label="")
        end
        push!(objective, y_best)
        pobj = plot(0:k, objective, label="objective")
        xlabel!("iteration")
        ylabel!("best objective")
        plot(swarmplot, pobj)
    end
    return anim
end

population = init_population(50, x1lims, x2lims)
pso_no_momentum = pso_animation(fun, population, 100, w=1, c1=0.9, c2=0.9)
gif(pso_no_momentum, "figures/pso_no_momentum.gif", fps=4)

population = init_population(50, x1lims, x2lims)
pso_momentum = pso_animation(fun, population, 100, w=0.8, c1=0.9, c2=0.9)
gif(pso_momentum, "figures/pso_momentum.gif", fps=4)


text:
module TestResidualJacobian
include("preamble.jl")
include("sym_diff.jl")

using StaticArrays: SVector

using .SymDiff: test_residual_jacobian
using Bifurcations.BifurcationsBase: SaddleNodeCont, HopfCont
using Bifurcations.Codim2: NormalizingAS, BackReferencingAS
using Bifurcations.Examples: examples

EXAMPLES = [
    (name, ex) for (name, ex) in examples()
    if try
        ex.make_codim2_prob
        true
    catch
        false
    end
]

@testset "residual_jacobian $name" for (name, ex) in EXAMPLES
    @testset "$(nameof(ck)) $(nameof(ASType)) $(utype)" for
        utype in [:SVector, :Vector],
        ASType in [BackReferencingAS, NormalizingAS],
        ck in [SaddleNodeCont, HopfCont]

        if utype == :SVector
            u0 = SVector(ex.u0...)
        elseif utype == :Vector
            u0 = Vector(ex.u0)
        else
            error("utype: $utype")
        end

        v0 = copy(u0)
        w0 = 0
        @assert ck in (SaddleNodeCont, HopfCont)
        if ck == HopfCont
            v0 = v0 .+ 1im
            w0 = 1.0
        end

        prob = ex.make_codim2_prob(
            u0 = u0,
            v0 = v0,
            w0 = w0,
            augmented_system = ASType(),
        )
        test_residual_jacobian(prob)
    end
end

end  # module


text:
abstract type AbstractImportanceSampling end
abstract type UnweightedIS <: AbstractImportanceSampling end
abstract type WeightedIS <: AbstractImportanceSampling end

struct IS <: UnweightedIS end
struct PDIS <: UnweightedIS end
struct WIS <: WeightedIS end
struct WPDIS <: WeightedIS end
struct CWPDIS <: WeightedIS end

include("bandit_is.jl")
include("trajectory_is.jl")

function estimate_returns(ism::AbstractImportanceSampling, π, H)
    return @. estimate_return((ism,), (π,), H)
end

function estimate_returns!(G, ism::AbstractImportanceSampling, π, H)
    @. G = estimate_return((ism,), (π,), H)
end

function estimate_returns_withentropy(ism::AbstractImportanceSampling, π, H, α)
    return @. estimate_return_withentropy((ism,), (π,), H, α)
end

function estimate_returns_withentropy!(G, ism::AbstractImportanceSampling, π, H, α)
    @. G = estimate_return_withentropy((ism,), (π,), H, α)
end

function estimate_returns(ism::WIS, π, H)
    N = length(H)
    ρtot = 0.0
    rets = estimate_return.((ism,), (π,), H)
    G = [r[1] for r in rets]

    for i in 1:N
        ρi = rets[i][2]
        ρtot += ρi
    end
    G ./ ρtot
end

function estimate_returns_withentropy(ism::WIS, π, H, α)
    N = length(H)
    ρtot = 0.0
    rets = estimate_return_withentropy.((ism,), (π,), H, (α,))
    G = [r[1] for r in rets]

    for i in 1:N
        ρi = rets[i][2]
        ρtot += ρi
    end
    G / ρtot
end

function estimate_returns!(G, ism::WIS, logpdf_fn, H)
    N = length(H)
    ρtot = 0.0
    for i in 1:N
        Gi, ρi = estimate_return(ism, logpdf_fn, H[i])
        G[i] = Gi
        ρtot += ρi
    end
    @. G ./= ρtot
end

function estimate_returns_withentropy!(G, ism::WIS, π, H, α)
    N = length(H)
    ρtot = 0.0
    for i in 1:N
        Gi, ρi = estimate_return_withentropy(ism, π, H[i], α)
        G[i] = Gi
        ρtot += ρi
    end
    @. G ./= ρtot
end


text:
using LinearAlgebra, Test
using Combinatorics

using MultivariatePolynomials
const MP = MultivariatePolynomials

function div_number_test()
    Mod.@polyvar x
    @test div(2x, 2) == x
    @test div(6x + 9x^2, 3) == 2x + 3x^2
    @test div(6x + 9x^2, 4) == x + 2x^2
    if VERSION >= v"1.6"
        @test div(10x^3, 4, RoundUp) == 3x^3
        @test div(6x + 9x^2, 4, RoundUp) == 2x + 3x^2
    end
end

function gcd_lcm_test()
    Mod.@polyvar x y z
    @test gcd(x*y, y) == y
    @test lcm(x*y, y) == x*y
    @test gcd(x*y^2*z, y*z^3) == y*z
    @test lcm(x*y^2*z, y*z^3) == x*y^2*z^3
    @test gcd(x^2*y^7*z^3, x^4*y^5*z^2) == x^2*y^5*z^2
    @test lcm(x^2*y^7*z^3, x^4*y^5*z^2) == x^4*y^7*z^3
end

# Taken from
# Ideals, Varieties, and Algorithms
# Cox, Little and O'Shea, Fourth edition
# They have been adapted to the grlex ordering
function divides_test()
    Mod.@polyvar x y z
    @test divides(leadingmonomial(x+y), x) # leadingmonomial(x+y) will be x^1*y^0 -> tricky test !
    @test !divides(leadingmonomial(x^2+y), x)
    @test divides(x*y, x^2*y)
    @test divides(x*y, x*y^2)
    @test divides(y*z, x*y*z)
    @test !divides(y*z, x*z)
    @test !divides(x^2*y, x*y^2)
end

function leading_term_test()
    Mod.@polyvar x y z
    @test @inferred(monic(2x^2 + 4y + 2)) == x^2 + 2y + 1
    # See page 60
    p = 4x*y^2*z + 4z^2 + 7x^2*z^2 - 5x^3
    @test leadingcoefficient(p) == 7
    @test leadingmonomial(p) == x^2*z^2
    @test leadingterm(p) == 7x^2*z^2
end

function divrem_test()
    Mod.@polyvar x y
    @test (@inferred div(x*y^2 + 1, x*y + 1)) == y
    @test (@inferred rem(x*y^2 + 1, x*y + 1)) == -y + 1
    @test (@inferred div(x*y^2 + x, y)) == x*y
    @test (@inferred rem(x*y^2 + x, y)) == x
    @test (@inferred rem(x^4 + x^3 + (1+1e-10)*x^2 + 1, x^2 + x + 1)) == 1
    @test (@inferred rem(x^4 + x^3 + (1+1e-10)*x^2 + 1, x^2 + x + 1; ztol=1e-11)) ≈ -((1+1e-10)-1)x + 1
end

function testdiv(p, ps)
    q, r = @inferred divrem(p, ps)
    @test p == dot(q, ps) + r
    q, r
end

function multi_div_test()
    Mod.@polyvar x y
    # Example 1
    q, r = testdiv(x*y^2 + 1, [x*y + 1, y + 1])
    @test q == [y, -1]
    @test r == 2
    # Example 2
    q, r = testdiv(x^2*y + x*y^2 + y^2, [x*y - 1, y^2 - 1])
    @test q == [x + y, 1]
    @test r == x + y + 1
    # Example 4
    q, r = testdiv(x^2*y + x*y^2 + y^2, [y^2 - 1, x*y - 1])
    @test q == [x + 1, x]
    @test r == 2x + 1
    # Example 5
    q, r = testdiv(x*y^2 - x, [x*y - 1, y^2 - 1])
    @test q == [y, 0]
    @test r == -x + y
    q, r = testdiv(x*y^2 - x, [y^2 - 1, x*y - 1])
    @test q == [x, 0]
    @test r == 0

    @testset "Issue DynamicPolynomials#62" begin
        p = x^2 + x + 1
        q = rem(p, [x^2-y])
        @test q == x + y + 1
    end

    @test (@inferred rem(x^2*y + (1+1e-10)*x*y + 1, [x^2 + x, y + 1])) == 1
    @test (@inferred rem(x^2*y + (1+1e-10)*x*y + 1, [x^2 + x, y + 1]; ztol=1e-11)) == -((1+1e-10)-1)x + 1
end

function test_gcd_unit(expected, p1, p2, algo)
    g = @inferred gcd(p1, p2, algo)
    # it does not make sense, in general, to speak of "the" greatest common
    # divisor of u and v; there is a set of greatest common divisors, each
    # one being a unit multiple of the others [Knu14, p. 424] so `expected` and
    # `-expected` are both accepted.
    @test g == expected || g == -expected
end

function test_gcdx_unit(expected, p1, p2, algo)
    test_gcd_unit(expected, p1, p2, algo)
    a, b, g = @inferred gcdx(p1, p2, algo)
    # it does not make sense, in general, to speak of "the" greatest common
    # divisor of u and v; there is a set of greatest common divisors, each
    # one being a unit multiple of the others [Knu14, p. 424] so `expected` and
    # `-expected` are both accepted.
    @test iszero(MP.pseudo_rem(g, expected, algo)[2])
    @test a * p1 + b * p2 == g
end
function _test_gcdx_unit(expected, p1, p2, algo)
    test_gcdx_unit(expected, p1, p2, algo)
    test_gcdx_unit(expected, p2, p1, algo)
end

function univariate_gcd_test(algo=GeneralizedEuclideanAlgorithm())
    Mod.@polyvar x
    test_gcdx_unit(x + 1, x^2 - 1, x^2 + 2x + 1, algo)
    test_gcdx_unit(x + 1, x^2 + 2x + 1, x^2 - 1, algo)
    test_gcdx_unit(x + 1, x^2 - 1, x + 1, algo)
    test_gcdx_unit(x + 1, x + 1, x^2 - 1, algo)
    test_gcdx_unit(x + 1, x - x, x + 1, algo)
    test_gcdx_unit(x + 1, x + 1, x - x, algo)
    test_gcdx_unit(x - x + 1, x + 1, x + 2, algo)
    test_gcdx_unit(x - x + 1, x + 1, 2x + 1, algo)
    test_gcdx_unit(x - x + 1, x - 1, 2x^2 - 2x - 2, algo)
    @test       0  == @inferred gcd(x - x, x^2 - x^2, algo)
    @test       0  == @inferred gcd(x^2 - x^2, x - x, algo)
end

function _mult_test(a::Number, b)
    @test iszero(maxdegree(b))
end
function _mult_test(a::Number, b::Number)
    @test iszero(rem(a, b))
    @test iszero(rem(b, a))
end
function _mult_test(a, b)
    @test iszero(rem(a, b))
    @test iszero(rem(b, a))
end
function mult_test(expected, a, b, algo)
    g = @inferred MP._simplifier(a, b, algo)
    @test g isa promote_type(polynomialtype(a), polynomialtype(b))
    _mult_test(expected, g)
end
function mult_test(expected, a::Number, b, algo)
    g = @inferred MP._simplifier(a, b, algo)
    @test g isa promote_type(typeof(a), MP.coefficienttype(b))
    _mult_test(expected, g)
end
function mult_test(expected, a, b::Number, algo)
    g = @inferred MP._simplifier(a, b, algo)
    @test g isa promote_type(MP.coefficienttype(a), typeof(b))
    _mult_test(expected, g)
end
function sym_test(a, b, g, algo)
    mult_test(g, a, b, algo)
    mult_test(g, b, a, algo)
end
function triple_test(a, b, c, algo)
    sym_test(a * c, b * c, gcd(a, b, algo) * c, algo)
    sym_test(b * a, c * a, gcd(b, c, algo) * a, algo)
    sym_test(c * b, a * b, gcd(c, a, algo) * b, algo)
end

function multivariate_gcd_test(::Type{T}, algo=GeneralizedEuclideanAlgorithm()) where {T}
    Mod.@polyvar x y z
    o = one(T)
    zr = zero(T)
    sym_test(o, o * x, o, algo)
    sym_test(o * x, zr * x, o * x, algo)
    sym_test(o * x + o, zr * x, o * x + o, algo)
    # Inspired from https://github.com/JuliaAlgebra/MultivariatePolynomials.jl/issues/160
    f1 = o * x * y + o * x
    f2 = o * y^2
    f3 = o * x
    sym_test(f1, f2, 1, algo)
    sym_test(f2, f3, 1, algo)
    sym_test(f3, f1, x, algo)
    triple_test(f1, f2, f3, algo)

    @testset "Issue #173" begin
        p1 = o*x*y + x
        p2 = x^2
        sym_test(p1, p2, x, algo)
    end

    p1 = o*z - z
    p2 = z
    @test gcd(p1, p2, algo) == z
    p1 = o*y - y
    p2 = z
    @test gcd(p1, p2, algo) == z
    test_relatively_prime(p1, p2, algo) = test_gcd_unit(o, p1, p2, algo)
    test_relatively_prime(2o*y*z - o, y - z, algo)
    test_relatively_prime(2o*y^2*z - y, y - z, algo)
    test_relatively_prime(2o*y^2*z - y, y^3*z - y - z, algo)
    test_relatively_prime(2o*y^3 + 2*y^2*z - y, y^3 + y^3*z - y - z, algo)
    test_relatively_prime(z^4 + 2o*y^3 + 2*y^2*z - y, y*z^4 + y^3 + y^3*z - y - z, algo)
    test_relatively_prime(
        y*z^3 + z^4 + 2o*y^3 + 2*y^2*z - y,
        y^2*z^3 + y*z^4 + y^3 + y^3*z - y - z,
        algo,
    )
    if T != Int || (algo != GeneralizedEuclideanAlgorithm(false, false) && algo != GeneralizedEuclideanAlgorithm(true, false))
        test_relatively_prime(
            -3o*y^2*z^3 - 3*y^4 + y*z^3 + z^4 + 2*y^3 + 2*y^2*z - y,
            3o*y^3*z^3 - 2*y^5 + y^2*z^3 + y*z^4 + y^3 + y^3*z - y - z,
            algo,
        )
    end
    test_relatively_prime(
        -z^6 - 3o*y^2*z^3 - 3*y^4 + y*z^3 + z^4 + 2*y^3 + 2*y^2*z - y,
        -y*z^6 - 3o*y^3*z^3 - 2*y^5 + y^2*z^3 + y*z^4 + y^3 + y^3*z - y - z,
        algo,
    )
    test_relatively_prime(
        -z^6 - 3o*y^2*z^3 - 3*y^4 + y*z^3 + z^4 + 2*y^3 + 2*y^2*z - y,
        -y^2*z^6 - 3o*y^4*z^3 - 2*y^6 + y^3*z^3 + y^2*z^4 + y^5 + y^4*z - y^2 - y*z,
        algo,
    )

    a = (o * x + o * y^2) * (o * z^3 + o * y^2 + o * x)
    b = (o * x + o * y + o * z) * (o * x^2 + o * y)
    c = (o * x + o * y + o * z) * (o * z^3 + o * y^2 + o * x)
    if T != Int || (algo != GeneralizedEuclideanAlgorithm(false, false) && algo != GeneralizedEuclideanAlgorithm(true, false))
        sym_test(a, b, 1, algo)
    end
    sym_test(b, c, x + y + z, algo)
    sym_test(c, a, z^3 + y^2 + x, algo)
    if (T != Int || (algo != GeneralizedEuclideanAlgorithm(false, false) && algo != GeneralizedEuclideanAlgorithm(true, false) && algo != GeneralizedEuclideanAlgorithm(false, true))) &&
        (T != Float64 || (algo != GeneralizedEuclideanAlgorithm(false, true) && algo != GeneralizedEuclideanAlgorithm(true, true)))
        triple_test(a, b, c, algo)
    end
end

function lcm_test()
    Mod.@polyvar x
    l = @inferred lcm(x^2 - 1, x^2 + 2x + 1)
    exp = -(x-1) * (x+1)^2
    @test l == exp || l == -exp
    l = @inferred lcm(x^2 + 2x + 1, x^2 - 1)
    @test l == exp || l == -exp
    @test x^2 - 1 == @inferred lcm(x^2 - 1, x - 1)
    @test x^2 - 1 == @inferred lcm(x - 1, x^2 - 1)
    @test       0 == @inferred lcm(x - x, x + 1)
    @test       0 == @inferred lcm(x + 1, x - x)
    @test       0 == @inferred lcm(x - x, x^2 - x^2)
    @test       0 == @inferred lcm(x^2 - x^2, x - x)
end

function deflation_test()
    Mod.@polyvar a b c
    function _test(p, eshift, edefl)
        shift, defl = MP.deflation(p)
        @test shift == eshift
        @test defl == edefl
        q = MP.deflate(p, shift, defl)
        @test p == MP.inflate(q, shift, defl)
    end
    for (x, y, z) in permutations((a, b, c))
        _test(x + x^2 + y * x, x, x * y)
        p = x^2 + x^2 * y^3 + y^6 * x^4
        _test(p, x^2, x^2 * y^3)
        @test p == MP.deflate(p, z^0, z^0)
        _test(p, x^2, x^2 * y^3)
    end
end

function extracted_variable_test()
    Mod.@polyvar x y z
    function _test(p1, p2, w1, w2)
        i1, i2, n = MP._extracted_variable(p1, p2)
        v(p, i) = iszero(i) ? nothing : variables(p)[i]
        if string(Mod) == "DynamicPolynomials"
            alloc_test(() -> MP._extracted_variable(p1, p2), 0)
        end
        @test v(p1, i1) == w1
        @test v(p2, i2) == w2
        i2, i1, n = MP._extracted_variable(p2, p1)
        if string(Mod) == "DynamicPolynomials"
            alloc_test(() -> MP._extracted_variable(p2, p1), 0)
        end
        @test v(p1, i1) == w1
        @test v(p2, i2) == w2
    end
    _test(x - x, y - y, nothing, nothing)
    _test(x + y + z, x - z, y, nothing)
    _test(x + y + z, x - y, z, nothing)
    _test(x + y + z, y - z, x, nothing)
    _test(x^4 + y^5 + z^6, x^3 + y^2 + z, z, z)
    _test(x^6 + y^5 + z^4, x + y^2 + z^3, x, x)
    _test(x^6 + y + z^4 - y, x + y + z^3 - y, x, x)
end

@testset "Division" begin
    @testset "div by number" begin
        div_number_test()
    end
    @testset "GCD and LCM" begin
        gcd_lcm_test()
    end
    @testset "Divides" begin
        divides_test()
    end
    @testset "Leading function" begin
        leading_term_test()
    end
    @testset "Division examples" begin
        divrem_test()
    end
    @testset "Division by multiple polynomials examples" begin
        multi_div_test()
    end
    @testset "Univariate gcd primitive_rem=$primitive_rem" for primitive_rem in [false, true]
        @testset "skip_last=$skip_last" for skip_last in [false, true]
            univariate_gcd_test(GeneralizedEuclideanAlgorithm(primitive_rem, skip_last))
        end
    end
    @testset "Multivariate gcd $T" for T in [Int, BigInt, Rational{BigInt}, Float64]
        if T != Rational{BigInt} || VERSION >= v"1.6"
            # `gcd` for `Rational{BigInt}` got defined at some point between v1.0 and v1.6
            @testset "primitive_rem=$primitive_rem" for primitive_rem in [false, true]
                @testset "skip_last=$skip_last" for skip_last in [false, true]
                    multivariate_gcd_test(T, GeneralizedEuclideanAlgorithm(primitive_rem, skip_last))
                end
            end
        end
    end
    @testset "lcm" begin
        lcm_test()
    end
end


text:
# Automation and the Future of Work
# Michele Fornino and Andrea Manera

#---------------------------------------------------#
#					  BASE model 					#
#---------------------------------------------------#

"""
	FM19base_params{F<:Real, I<:Integer}

	Composite type containing all parameters of the Base model in Fornino Manera (2019).

"""
@with_kw struct FM19base_params{F<:Real, I<:Integer} @deftype F

	# Parameters
	θ = 0.5
	w = 0.5
	pR = 1.5
	Γ = 0.5
	m = 0.2
	δ = log(1 + 1/12)
	ρ = log(1 + 0.04)
	ψR = 25.0
	θz = 0.1
	σz = 0.4
	zd = 1.0
	p = 1.0

	# Derived Parameters
	Ω = (1.0 - Γ)/Γ * w - m
	Rmax =  1.0 ./ (δ .* ψR) .* (Ω ./ (ρ .+ δ) .- pR)

	# Grids Parameters
	NR::I = 200
	Nz::I = 200
	zGrid_lower_quantile = 1e-4
	zGrid_upper_quantile = 1 - 1e-4
	RGrid_lower = 0.0
	RGrid_upper = Rmax

	# Assertions
	@assert (θ > 0) & (θ < 1)
	@assert w > 0
	@assert pR > 0
	@assert (Γ > 0) & (Γ < 1)
	@assert m > 0
	@assert δ > 0
	@assert ρ > 0
	@assert ψR > 0
	@assert θz > 0
	@assert σz > 0
	@assert zd > 0
	@assert p > 0
	@assert Ω > 0 "Labor savings Ω cannot be negative."

	# Check values for grids make sense
	@assert NR > 1
	@assert Nz > 0
	@assert zGrid_upper_quantile > zGrid_lower_quantile
	@assert zGrid_lower_quantile > 0
	@assert zGrid_upper_quantile < 1
	@assert RGrid_lower >= 0

	# Enforce definitions (prevents accidental creation of object with wrong derived parameters)
	@assert isequal(Ω, (1.0 - Γ)/Γ * w - m) "Ω is a derived parameter. Please modify the underlying parameters if you wish to change it."
	@assert isequal(Rmax, 1.0 ./ (δ .* ψR) .* (Ω ./ (ρ .+ δ) .- pR)) "Rmax is a derived parameter. Please modify the underlying parameters if you wish to change it."

end

"""
	function HamiltonJacobiBellmanEquations.HJBEquation(params::FM19base_params)

Overloaded constructor for input of type FM19base_params.

"""
function HamiltonJacobiBellmanEquations.HJBEquation(params::FM19base_params)

	# Use UnPack.jl to localize all members of params
	@unpack_FM19base_params params

	# Profit Function
	function Π(R::T, z::T)::T where T<:Real
		Rb = 1.0/ (1.0 - Γ) * ((p * z * θ * Γ) / w) ^ (1.0 / (1.0 - θ))
		Rh = 1.0/ (1.0 - Γ) * ((p * z * θ * (1.0 - Γ)) / m) ^ (1.0 / (1.0 - θ))
		if R ≤ Rb
			static_profit = (1.0 - θ) * p * z * ((1.0 - Γ) * Rb) ^ θ + Ω * R
		elseif Rb ≤ R ≤ Rh
			static_profit = p * z * ((1.0 - Γ) * R) ^ θ - m * R
		else
			static_profit = p * z * ((1.0 - Γ) * Rh) ^ θ - m * Rh
		end
		return static_profit
	end

	# Returns Function
	function flow_returns(state::AbstractVector{T}, I::Union{T, AbstractVector{T}})::T where T<:Real
		return Π(state[1], state[2]) - pR * I[1] - ψR / 2 * I[1] ^ 2
	end

	# Policy Function in term of marginal value function
	function optimal_controls(state::AbstractVector{T}, ∂V∂R::Union{T, AbstractVector{T}})::T where T<:Real
		return 1 / ψR * (∂V∂R[1] - pR)
	end

	# State equations
	function state_equations(state::AbstractVector{T}, I::Union{T, AbstractVector{T}})::T where T<:Real
		return I[1] - δ * state[1]
	end

	# State constraint
	function state_constraints(state::AbstractVector{T})::T where T<:Real
		return ψR * δ * state[1] + pR
	end

	# Tuple of endogenous state grids
	RGrid = Vector(range(RGrid_lower, RGrid_upper, length = NR))
	state_grids_endog = (RGrid, )

	# Exogenous shock
	exog_shock = ContinuousTimeMarkovChain(params)

	# Allocate hjb object
	return HJBEquation(flow_returns,
					   state_equations,
					   optimal_controls,
					   state_constraints,
					   state_grids_endog,
					   exog_shock,
					   ρ)

end

"""
	function MarkovChains.ContinuousTimeMarkovChain(params::FM19base_params)

Overloaded constructor for input of type FM19base_params.

"""
function MarkovChains.ContinuousTimeMarkovChain(params::FM19base_params)

	@unpack_FM19base_params params

	# Construct exog_shock ContinuousTimeMarkovChain object (if problem is stochastic).
	if Nz > 1

		# Obtain discretize diffusion process object as a MarkovChains.ContinuousTimeMarkovChain.
		μ = z -> -θz .* (z .- zd)
		σ = z -> sqrt(2 * θz * σz^2 / zd) .* sqrt.(z)
		dp = ItoDiffusionProcess(μ, σ)

		# Theoretical benchmark for CIR stationary distribution, used for zGrid.
		shape = zd^2 / σz^2
		rate = zd / shape
		Fz_theory = Gamma(shape, rate)
		zGrid_lower = quantile(Fz_theory, zGrid_lower_quantile)
		zGrid_upper = quantile(Fz_theory, zGrid_upper_quantile)

		# Allocate Grid
		zGrid = Vector(range(zGrid_lower, zGrid_upper, length = Nz))

		# Exogenous Shock: ContinuousTimeMarkovChain object
		exog_shock = ContinuousTimeMarkovChain(dp, zGrid)

	elseif Nz == 1

		exog_shock = ContinuousTimeMarkovChain(reshape([0], 1, 1), [zd])

	end

	return exog_shock

end

function HamiltonJacobiBellmanEquations.solve(params::FM19base_params)

	return solve(HJBEquation(params))

end

# # IDEA?
# function plot(params::FM19base_params, sol::HJBSolution)

# 	# PLOTTING FOR 2-dim problem (1 end + 1 exo)

# 	# IDEA IS TO MAKE MULTIPLE FIGURES
# 	# FIG 1:
# 	# [subplot surf/contour of value function ]
# 	# [subplot slice x] ||| [subplot slice z  ]

# 	# FIG 2:
# 	# same as above but one for each policy c_j

# 	# FIG 3:
# 	# stationary distribution (countour)

# 	# PLOTTING FOR N-dim problem (N-1 end + 1 exo)
# 	# THINK ABOUT HOW TO MAKE IT EFFECTIVELY
# end

# # IDEA?
# function summarize(sol::HJBSolution)

# 	# OBTAIN STATISTICS REGARDING SOLUTION OBJECTS (MUST BE GENERIC FOR ANY PROBLEM)

# end


text:
using FastClosures

mutable struct BFGSData{T}
    M :: Matrix{T}
    scaling::Bool
    Ax :: Vector{T}  # pre allocated for mul5
end

function BFGSData(
    M₀ :: Matrix{T};
    scaling :: Bool = true ) where {T}

    n,  = size(M₀) # add checks for square symmetric matrix

    BFGSData{T}(M₀, scaling, Vector{T}(undef, n))
end

"A type for BFGS matrix updates."
mutable struct BFGSOperator{T, N <: Integer, F, Ft, Fct} <: AbstractLinearOperator{T}
    nrow::N
    ncol::N
    symmetric::Bool
    hermitian::Bool
    prod!::F    # apply the operator to a vector
    tprod!::Ft    # apply the transpose operator to a vector
    ctprod!::Fct   # apply the transpose conjugate operator to a vector
    #inverse::Bool
    data::BFGSData{T}
    nprod::N
    ntprod::N
    nctprod::N
end

BFGSOperator{T}(
    nrow :: N,
    ncol :: N,
    symmetric :: Bool,
    hermitian :: Bool,
    prod!::F,
    tprod!::Ft,
    ctprod!::Fct,
    data::BFGSData{T},
) where {T, N, F, Ft, Fct} = BFGSOperator{T, N, F, Ft, Fct}(
    nrow,
    ncol,
    symmetric,
    hermitian,
    prod!,
    tprod!,
    ctprod!,
    data,
    0,
    0,
    0,
)

import LinearOperators:has_args5, use_prod5!, isallocated5
has_args5(op::BFGSOperator) = true
use_prod5!(op::BFGSOperator) = true
isallocated5(op::BFGSOperator) = true

"""
    InverseBFGSOperator(M₀, n [; scaling=true])
    InverseBFGSOperator(n, [; scaling=true])
Construct a BFGS approximation in inverse form.
"""
function InverseBFGSOperator(M :: Matrix{T}; kwargs...) where {T <: Real}
    kwargs = Dict(kwargs)
    bfgs_data = BFGSData(M; kwargs...)
    n,  = size(M)

    function bfgs_multiply(res::AbstractVector,
                           data::BFGSData,
                           x::AbstractArray,
                           αm,
                           βm::T2,
                           ) where T2

        q = data.Ax  # pre allocated
        q .= data.M * x

        # mul5 stuff
        if βm == zero(T2)
            res .= αm .* q
        else
            res .= αm .* q .+ βm .* res
        end

        return res
    end

    prod! = @closure (res, x, α, β) -> bfgs_multiply(res, bfgs_data, x, α, β)
    return BFGSOperator{T}(n, n, true, true, prod!, prod!, prod!, bfgs_data)
end

function InverseBFGSOperator(T, n :: N; kwargs...) where {N <: Integer}
    Eye = Matrix{T}(I,n,n)
    InverseBFGSOperator(Eye; kwargs...)
end

#function Matrix(Op::InverseBFGSOperator)
#    return Op.data.M
#end

max_stars_repo_name: vepiteski/LSDescent

max_issues_repo_name: vepiteski/LSDescent

max_forks_repo_name: vepiteski/LSDescent


text:
export name, connected, disconnected, disconnect, url, set_url!,
       push_url, set_push_url!, fetch_refspecs, push_refspecs,
       add_fetch!, add_push!, clear_refspecs!, save!, rename!,
       update_tips!

function check_valid_url(url::AbstractString)
    if !bool(ccall((:git_remote_valid_url, libgit2), Cint, (Ptr{UInt8},), url))
        throw(ArgumentError("Invalid URL : $url"))
    end
    return true
end

GitRemote(r::GitRepo, url::AbstractString) = begin
    check_valid_url(url)
    remote_ptr = Ptr{Void}[0]
    @check ccall((:git_remote_create_anonymous, libgit2), Cint,
                 (Ptr{Ptr{Void}}, Ptr{Void}, Ptr{UInt8}, Ptr{Void}), remote_ptr, r, url, C_NULL)
    return GitRemote(remote_ptr[1])
end

Base.show(io::IO, r::GitRemote) = begin
    write(io, "GitRemote(\"$(name(r))\" @ $(url(r)))")
end

function name(r::GitRemote)
    name_ptr = ccall((:git_remote_name, libgit2), Ptr{UInt8}, (Ptr{Void},), r)
    return name_ptr != C_NULL ? bytestring(name_ptr) : nothing
end

connected(r::GitRemote) = bool(ccall((:git_remote_connected, libgit2), Cint, (Ptr{Void},), r))
disconnected(r::GitRemote) = !connected(r)

function disconnect(r::GitRemote)
    @check ccall((:git_remote_disconnect, libgit2), Cint, (Ptr{Void},), r)
    return r
end

# direction can be :fetch or :push
Base.connect(r::GitRemote, direction::Symbol) = begin
    local dir::Cint
    dir = direction == :fetch ? GitConst.DIRECTION_FETCH :
          direction == :push  ? GitConst.DIRECTION_PUSH  :
          throw(ArgumentError("direction can be :fetch or :push, got :$direction"))
    @check ccall((:git_remote_connect, libgit2), Cint, (Ptr{Void}, Cint), r, dir)
    return r
end

Base.connect(f::Function, r::GitRemote, direction::Symbol) = begin
    connect(r, direction)
    try
        f(r)
    finally
        disconnect(r)
    end
end

function url(r::GitRemote)
    url_ptr = ccall((:git_remote_url, libgit2), Ptr{UInt8}, (Ptr{Void},), r)
    return url_ptr != C_NULL ? bytestring(url_ptr) : nothing
end

function set_url!(r::GitRemote, url::AbstractString)
    check_valid_url(url)
    @check ccall((:git_remote_set_url, libgit2), Cint, (Ptr{Void}, Ptr{UInt8}), r, url)
    return r
end

function push_url(r::GitRemote)
    url_ptr = ccall((:git_remote_pushurl, libgit2), Ptr{UInt8}, (Ptr{Void},), r)
    return url_ptr != C_NULL ? bytestring(url_ptr) : nothing
end

function set_push_url!(r::GitRemote, url::AbstractString)
    check_valid_url(url)
    @check ccall((:git_remote_set_pushurl, libgit2), Cint, (Ptr{Void}, Ptr{UInt8}), r, url)
    return r
end

function fetch_refspecs(r::GitRemote)
    sa_ptr = [StrArrayStruct()]
    @check ccall((:git_remote_get_fetch_refspecs, libgit2), Cint,
                  (Ptr{StrArrayStruct}, Ptr{Void}), sa_ptr, r)
    sa = sa_ptr[1]
    refs = Array(UTF8String, sa.count)
    for i in 1:sa.count
        refs[i] = utf8(bytestring(unsafe_load(sa.strings, i)))
    end
    return refs
end

function push_refspecs(r::GitRemote)
    sa_ptr = [StrArrayStruct()]
    @check ccall((:git_remote_get_push_refspecs, libgit2), Cint,
                  (Ptr{StrArrayStruct}, Ptr{Void}), sa_ptr, r)
    sa = sa_ptr[1]
    refs = Array(UTF8String, sa.count)
    for i in 1:sa.count
        refs[i] = utf8(bytestring(unsafe_load(sa.strings, i)))
    end
    return refs
end

function add_push!(r::GitRemote, ref::AbstractString)
    @check ccall((:git_remote_add_push, libgit2), Cint, (Ptr{Void}, Ptr{UInt8}), r, ref)
    return r
end

function add_fetch!(r::GitRemote, ref::AbstractString)
    @check ccall((:git_remote_add_fetch, libgit2), Cint, (Ptr{Void}, Ptr{UInt8}), r, ref)
    return r
end

function clear_refspecs!(r::GitRemote)
    ccall((:git_remote_clear_refspecs, libgit2), Void, (Ptr{Void},), r)
    return r
end

function save!(r::GitRemote)
    @check ccall((:git_remote_save, libgit2), Cint, (Ptr{Void},), r)
    return r
end

function rename!(r::GitRemote, newname::AbstractString)
    sa_ptr = [StrArrayStruct()]
    @check ccall((:git_remote_rename, libgit2), Cint,
                  (Ptr{StrArrayStruct}, Ptr{Void}, Ptr{UInt8}), sa_ptr, r, newname)
    sa = sa_ptr[1]
    errs = UTF8String[]
    for i=1:sa.count
        push!(errs, utf8(bytestring(unsafe_load(sa.strings, i))))
    end
    return errs
end

type RemoteHead
    islocal::Bool
    id::Oid
    lid::MaybeOid
    name::ByteString
end

RemoteHead(struct::RemoteHeadStruct) =  begin
    id  = Oid(string(struct.id))
    lid = Oid(string(struct.lid))
    name = struct.name != C_NULL ? "" : bytestring(struct.name)
    return RemoteHead(bool(struct.islocal), id, iszero(lid) ? nothing : lid, name)
end

ls(r::GitRemote) = begin
    nheads = Csize_t[0]
    head_ptrs = Ptr{Ptr{RemoteHeadStruct}}[0]
    @check ccall((:git_remote_ls, libgit2), Cint,
                 (Ptr{Ptr{Ptr{RemoteHeadStruct}}}, Ptr{Csize_t}, Ptr{Void}), head_ptrs, nheads, r)
    head_ptr = head_ptrs[1]
    remote_heads = RemoteHead[]
    for i in 1:nheads[1]
        struct = unsafe_load(unsafe_load(head_ptr, i))::RemoteHeadStruct
        push!(remote_heads, RemoteHead(struct))
    end
    return remote_heads
end

Base.download(r::GitRemote) = begin
    @check ccall((:git_remote_download, libgit2), Cint, (Ptr{Void},), r)
    return r
end

function update_tips!(r::GitRemote,
                     sig::MaybeSignature=nothing,
                     logmsg::MaybeString=nothing)
    @check ccall((:git_remote_update_tips, libgit2), Cint,
                 (Ptr{Void}, Ptr{SignatureStruct}, Ptr{UInt8}),
                 r, sig != nothing ? sig : C_NULL, logmsg != nothing ? logmsg : C_NULL)
    return r
end

function fetch(r::GitRemote, sig::MaybeSignature=nothing, logmsg::MaybeString=nothing)
    @check ccall((:git_remote_fetch, libgit2), Cint,
                 (Ptr{Void}, Ptr{SignatureStruct}, Ptr{UInt8}),
                 r, sig != nothing ? sig : C_NULL, logmsg != nothing ? logmsg : C_NULL)
    return r
end

function default_branch(r::GitRemote)
    buf_ptr = [BufferStruct()]
    @check ccall((:git_remote_default_branch, libgit2), Cint,
                 (Ptr{BufferStruct}, Ptr{Void}), buf_ptr, r)
    name = utf8(bytestring(buf_ptr[1]))
    ccall((:git_buf_free, libgit2), Void, (Ptr{BufferStruct},), buf_ptr)
    return name
end