text:
"""
    EXP3 Implementation
"""

type EXP3 <: BanditAlgorithmBase
    noOfArms::Int64
    noOfSteps::Int64
    lastPlayedArm::Int64

    γ::Float64
    wVec::Vector
    pDist::Categorical

    function EXP3( noOfArms::Integer, γ::Real )
        new( noOfArms,
             0,
             0,
             γ,
             ones(noOfArms),
             Categorical(1/noOfArms*ones(noOfArms))
        )
    end
end

function get_arm_index( agent::EXP3 )
    agent.lastPlayedArm = rand( agent.pDist )
    return agent.lastPlayedArm
end

function update_reward!( agent::EXP3, r::Real )
    agent.noOfSteps     = agent.noOfSteps + 1;
    # Calculate estimated reward
    r_est = r/agent.pDist.p[agent.lastPlayedArm]

    # Update weight of arm
    agent.wVec[agent.lastPlayedArm] = agent.wVec[agent.lastPlayedArm] * exp(agent.γ*r_est/agent.noOfArms)

    # Calculate new probabilties
    p = (1-agent.γ) * agent.wVec/sum(agent.wVec) + agent.γ/agent.noOfArms

    # Make it a distrubution
    agent.pDist = Categorical( p )

    nothing
end

function reset!( agent::EXP3 )
    agent.noOfSteps     = 0
    agent.lastPlayedArm = 0
    agent.wVec          = ones( agent.noOfArms )
    agent.pDist         = Categorical( 1/agent.noOfArms*ones(agent.noOfArms) )

    nothing
end

function info_str( agent::EXP3, latex::Bool )
    if latex
        return @sprintf( "EXP3 (\$\\gamma = %4.3f\$)", agent.γ )
    else
        return @sprintf( "EXP3 (γ = %4.3f)", agent.γ )
    end
end

"""
    EXP3.1 Implementation
    Based on : Fig.2 in Auer, P., Bianchi, N. C., Freund, Y., & E.Schapire, R. (2002). The Non-Stochastic Multi-Armed bandit problem. SIAM Journal of Computing, 22, 322–331. http://doi.org/10.1109/CDC.1983.269708
"""

type EXP31 <: BanditAlgorithmBase
    G_hat::Vector
    r::Int64
    g_r::Float64

    _EXP3::EXP3

    function EXP31( noOfArms )
        new( zeros( noOfArms ),
             0,
             (noOfArms*log(noOfArms))/(e-1),
             EXP3( noOfArms, 1.00 )  # Initial value for γ will simplify to 1.00
        )
    end
end

function get_arm_index( agent::EXP31 )
    get_arm_index( agent._EXP3 )
end

function update_reward!( agent::EXP31, r::Real )
    # Update G_hat for pulled arm
    agent.G_hat[agent._EXP3.lastPlayedArm] += r/agent._EXP3.pDist.p[agent._EXP3.lastPlayedArm]
    # println( "    G_hat = ", maximum(agent.G_hat), "    LHS = ", (agent.g_r - agent._EXP3.noOfArms/agent._EXP3.γ) )
    if maximum(agent.G_hat) > (agent.g_r - agent._EXP3.noOfArms/agent._EXP3.γ)
        # Calculate new g_r and reset EXP3 with new γ
        agent.r     += 1
        agent.g_r   = (agent._EXP3.noOfArms*log(agent._EXP3.noOfArms))/(e-1) * (4^agent.r)
        # println( "        g_r = ", agent.g_r )
        # println( "        New γ = ", √((agent._EXP3.noOfArms*log(agent._EXP3.noOfArms))/((e-1)*agent.g_r )) )
        agent._EXP3 = EXP3( agent._EXP3.noOfArms,
                            min(
                                1.00,
                                √((agent._EXP3.noOfArms*log(agent._EXP3.noOfArms))/((e-1)*agent.g_r ))
                            )
                        )
    else
        # Update reward to EXP3
        update_reward!( agent._EXP3, r )
    end

    nothing
end

function reset!( agent::EXP31 )
    reset!( agent._EXP3 )
    agent.G_hat = zeros( agent._EXP3.noOfArms)
    agent.r     = 0

    nothing
end

function info_str( agent::EXP31, latex::Bool )
    return @sprintf( "EXP3.1" )
end

"""
    REXP3 Implementation
    Based on Besbes, O., Gur, Y., & Zeevi, A. (2014). Stochastic Multi-Armed-Bandit Problem with Non-stationary Rewards. Advances in Neural Information Processing Systems, 2, 1–9.
"""

type REXP3 <: BanditAlgorithmBase
    noOfSteps::Int64    # Timesteps passed
    j::Int64            # Batch index
    Δ::Int64            # Batch size

    _EXP3::EXP3

    function REXP3( noOfArms::Integer, γ::Real, Δ::Integer )
        new( 0,
             1,
             Δ,
             EXP3(noOfArms,γ)
        )
    end
end

function get_arm_index( agent::REXP3 )
    get_arm_index( agent._EXP3 )
end

function update_reward!( agent::REXP3, r::Real )
    agent.noOfSteps = agent.noOfSteps + 1;
    update_reward!( agent._EXP3, r );
    # Reset if necessary
    if agent._EXP3.noOfSteps >= agent.Δ
        reset!( agent._EXP3 );
        agent.j = agent.j + 1;
    end

    nothing
end

function reset!( agent::REXP3 )
    reset!( agent._EXP3 )
    agent.noOfSteps     = 0;
    agent.j             = 1;

    nothing
end

function info_str( agent::REXP3, latex::Bool )
    if latex
        return @sprintf( "REXP3 (\$\\gamma = %4.3f, \\Delta = %d\$)", agent._EXP3.γ, agent.Δ )
    else
        return @sprintf( "REXP3 (γ = %4.3f, Δ = %d)", agent._EXP3.γ, agent.Δ )
    end
end

"""
    EXP3-IX Implementation
    Based on Neu, G. (2015). Explore no more: improved high-probability regret bounds for non-stochastic bandits. In Advances in Neural Information Processing Systems (pp. 1–9).
"""

type EXP3IX <: BanditAlgorithmBase
    noOfArms::Int64
    noOfSteps::Int64
    lastPlayedArm::Int64

    η::Float64
    γ::Float64
    wVec::Vector
    pDist::Categorical

    function EXP3IX( noOfArms::Integer, η::Real, γ::Real )
        new( noOfArms,
             0,
             0,
             η,
             γ,
             ones(noOfArms),
             Categorical(1/noOfArms*ones(noOfArms))
        )
    end
end

function get_arm_index( agent::EXP3IX )
    agent.lastPlayedArm = rand( agent.pDist )
    return agent.lastPlayedArm
end

function update_reward!( agent::EXP3IX, r::Real )
    agent.noOfSteps     = agent.noOfSteps + 1;

    # Calculate loss
    l   = 1 - r

    # Calculate estimated reward
    l_est = l/(agent.pDist.p[agent.lastPlayedArm]+agent.γ)

    # Update weight of arm
    agent.wVec[agent.lastPlayedArm] = agent.wVec[agent.lastPlayedArm] * exp(-agent.η*l_est)

    # Calculate new probabilties
    p = agent.wVec/sum(agent.wVec)

    # Make it a distrubution
    agent.pDist = Categorical( p )

    nothing
end

function reset!( agent::EXP3IX )
    agent.noOfSteps     = 0
    agent.lastPlayedArm = 0
    agent.wVec          = ones( agent.noOfArms )
    agent.pDist         = Categorical( 1/agent.noOfArms*ones(agent.noOfArms) )

    nothing
end

function info_str( agent::EXP3IX, latex::Bool )
    if latex
        return @sprintf( "EXP3-IX (\$\\eta = %4.3f, \\gamma = %4.3f\$)", agent.η, agent.γ )
    else
        return @sprintf( "EXP3-IX (η = %4.3f, γ = %4.3f)", agent.η, agent.γ )
    end
end


text:
@testset "Blackbody T=$T" for T in [2000, 4000, 6000]
    wave = range(1e3, 5e4, length = 1000)
    b = 2.897771955185172e7

    bb = @inferred blackbody(wave, T)
    @test typeof(bb) <: Spectra.Spectrum
    @test bb.T == T
    @test bb.wave[argmax(bb)] ≈ b / T rtol = 0.01

    wave *= u"angstrom"
    T *= u"K"
    bb = @inferred blackbody(wave, T)
    @test typeof(bb) <: Spectra.Spectrum
    @test unit(bb)[2] == u"W/m^2/angstrom"
    @test bb.T == T
    @test bb.wave[argmax(bb)] ≈ b * u"angstrom*K" / T rtol = 0.01
end

ext: jl


text:
"""
# Earth

This module defines seismic models for the Earth.
"""
module Earth

using ..SeisModels

"""
	ALL_MODELS

List of Earth model instances as `Symbol`s.

If adding new inbuilt Earth models, append the name of the new model to this list
"""
const ALL_MODELS = (:AK135, :EK137, :IASP91, :PREM, :PREM_NOOCEAN, :STW105)

include("ak135.jl")
include("ek137.jl")
include("iasp91.jl")
include("prem.jl")
include("stw105.jl")

end # module


text:
import MixedModels: LinearMixedModel,
                    setθ!,
                    updateL!
import RCall: rcopy,
              RClass,
              rcopytype,
              reval,
              S4Sxp,
              sexp,
              protect,
              unprotect,
              sexpclass,
              @rput,
              @rget,
              @R_str

# if RCall is available, then so is DataFrames
import DataFrames: DataFrame
import Tables: ColumnTable

# from R
# note that weights are not extracted
# TODO: document weights issue and warn
function rcopy(::Type{LinearMixedModel}, s::Ptr{S4Sxp})

    # these try blocks should probably be changed to an examination of the indices
    # this only extracts the name within the call, not the actual weights
    try
        wts = rcopy(s[:call][:weights])
        @error "weights are not supported"
    catch err
        if !isa(err, BoundsError) # this is the error we were expecting
            rethrow(err)
        end
        # no weights defined, we continue on our way
    end

    try
        contrasts = rcopy(s[:call][:contrasts])
        @error "Contrasts must be specified in the dataframe, not the lmer() call"
    catch err
        if !isa(err, BoundsError) # this is the error we were expecting
            rethrow(err)
        end
        # no extra contrasts defined, we continue on our way
    end

    # for some reason this doesn't always give a formula with lmerTest
    #f = rcopy(s[:call][:formula])
    f = rcopy(R"as.formula($(s)@call$formula)")
    data = rcopy(s[:frame])
    contrasts = get_r_contrasts(s[:frame])
    θ = rcopyarray(s[:theta])
    reml = rcopy(s[:devcomp][:dims][:REML]) ≠ 0

    m = LinearMixedModel(f, data, contrasts=contrasts)

    if length(θ) != length(m.θ)
        @error """You're probably using || in R with a categorical variable,
                  whose translation is currently unsupported with MixedModels 3.0."""
        throw(ArgumentError("Parameter vectors in R and Julia are different sizes."))
    end

    θ = _reorder_theta_from_lme4(θ, m)

    m.optsum.REML = reml
    m.optsum.feval = rcopy(s[:optinfo][:feval])
    # I'm wondering if this be filled in from the Julia side
    m.optsum.final = rcopyarray(s[:optinfo][:val])
    m.optsum.optimizer = Symbol("$(rcopy(s[:optinfo][:optimizer])) (lme4)")
    m.optsum.returnvalue = rcopy(s[:optinfo][:conv][:opt]) == 0 ? :FAILURE : :SUCCESS
    m.optsum.fmin = reml ? rcopy(s[:devcomp][:cmp][:REML]) : rcopy(s[:devcomp][:cmp][:dev])
    updateL!(setθ!(m, θ))
end

rcopytype(::Type{RClass{:lmerMod}}, s::Ptr{S4Sxp}) = LinearMixedModel

# add lmerTest::lmer and afex::lmer_alt support
rcopytype(::Type{RClass{:lmerModLmerTest}}, s::Ptr{S4Sxp}) = LinearMixedModel

# TODO: fix some conversions -- Julia->R->Julia roundtrip currently due to
#        ERROR: REvalError: Error in function (x, value, pos = -1, envir = as.environment(pos), inherits = FALSE,  :
#          SET_VECTOR_ELT() can only be applied to a 'list', not a 'character'
function sexp(::Type{RClass{:lmerMod}}, x::Tuple{LinearMixedModel{T}, DataFrame}) where T
    m, tbl = x
    if !isempty(m.sqrtwts)
        @error "weights are not currently supported"
    end

    m.optsum.feval > 0 || throw(ArgumentError("Model must be fitted"))

    jellyme4_data = tbl
    formula = convert_julia_to_r(m.formula)

    θ = m.θ
    rsteps = 1
    REML = m.optsum.REML ? "TRUE" : "FALSE"
    jellyme4_theta = _reorder_theta_to_lme4(m)
    fval = m.optsum.fmin
    feval = m.optsum.feval
    conv = m.optsum.returnvalue == :SUCCESS ? 0 : 1
    optimizer = String(m.optsum.optimizer)
    message = "fit with MixedModels.jl"
    @rput jellyme4_data
    @rput jellyme4_theta

    set_r_contrasts!("jellyme4_data", m.formula)

    r = """
    jellyme4_mod <- $LMER(formula = $(formula),
                           data=jellyme4_data,
                           REML=$(REML),
                           control=lme4::lmerControl(optimizer="nloptwrap",
                                                     optCtrl=list(maxeval=$(rsteps)),
                                    calc.derivs=FALSE,
                                    check.nobs.vs.nRE= "warning"),
                            start=list(theta=jellyme4_theta))
     jellyme4_mod@optinfo\$feval <- $(feval)
     jellyme4_mod@optinfo\$message <- "$(message)"
     jellyme4_mod@optinfo\$optimizer <- "$(optimizer)"
     jellyme4_mod
    """
    r = reval(r)
    r = protect(sexp(r))
    unprotect(1)
    r
end

sexpclass(x::Tuple{LinearMixedModel{T}, DataFrame}) where T = LMER in ("afex::lmer_alt", "lmer_alt") ? RClass{:lmerModLmerTest} : RClass{:lmerMod}
sexp(::Type{RClass{:lmerModLmerTest}}, x::Tuple{LinearMixedModel{T}, DataFrame}) where T = sexp(RClass{:lmerMod}, x)

# generalize to ColumnTable, which is what MixedModels actually requires
function sexp(ss::Type{RClass{:lmerMod}}, x::Tuple{LinearMixedModel{T}, ColumnTable}) where T
    m, t  = x
    sexp(ss, (m, DataFrame(t)))
end

sexpclass(x::Tuple{LinearMixedModel{T}, ColumnTable}) where T = RClass{:lmerMod}


text:
"""
    localtest.jl -- Feel free to modify anything you like in this file!

You can call this file as a script at the command line, or `include` it from within a julia session.

Example command-line usage:

    julia --color=yes localtest.jl                  # Evaluate all simple functions the default number of times (500):
    julia --color=yes localtest.jl 2000             # Evaluate all simple functions a custom number of times (2000 in this case):
    julia --color=yes localtest.jl simple1          # Evaluate only simple1 the default number of times:
    julia --color=yes localtest.jl simple3 2000     # Evaluate only simple3 a custom number of times:
"""

# Include the other relevant files:
include(joinpath("project2_jl", "helpers.jl"))
include(joinpath("project2_jl", "project2.jl"))

if length(ARGS) == 0
    # Default
    K = 500
    probnames = sort(collect(keys(PROBS)))
elseif length(ARGS) == 1
    # One input argument can be either evaluation count or function name, so
    # have to test for both. If it's not an integer it's got to be a function name
    K = tryparse(Int, ARGS[1])
    if K === nothing
        K = 500
        probnames = [ARGS[1]]
    else
        probnames = sort(collect(keys(PROBS)))
    end
elseif length(ARGS) == 2
    # Two inputs can come in at any order, so try both
    K = tryparse(Int, ARGS[1])
    if K === nothing
        K = tryparse(Int, ARGS[2])
        if K === nothing
            throw(ArgumentError("Can't detect which input is intended to be the evaluation count. Make sure it's an integer"))
        end
        probnames = [ARGS[1]]
    else
        probnames = [ARGS[2]]
    end
else
    throw(ArgumentError("Too many command-line inputs to localtest.jl"))
end

printstyled("\nTesting $K times\n\n", bold = true)
for nm in probnames
    try
        scores, n_evals, optima = main(nm, K, optimize)

        if maximum(n_evals) > PROBS[nm].n
            @warn "number of exaluations exceeded. Got $(maximum(n_evals)) on $nm."
        end

        # Check which optima are feasible
        feasible_optima = map(optima) do x
            all(<=(0), PROBS[nm].c(x))
        end

        # 95% of optima must be feasible!
        pass = sum(feasible_optima) >= (0.95*K)

        if pass
            printstyled("Pass: optimize returns a feasible solution on $(sum(feasible_optima))/$K random seeds.\n", color = :green)
        else
            printstyled("Fail: optimize returns a feasible solution on $(sum(feasible_optima))/$K random seeds.\n", color = :red)
        end

    catch e
        println("ERROR IN PROBLEM $nm")
        showerror(stdout, e, catch_backtrace())
    end
end


text:
import LaplacianOpt as LOpt
using JuMP
using CPLEX
using Gurobi
#using GLPK

include("solver.jl")

result_all = zeros(10,3)
lopt_optimizer = get_gurobi()

for i=1:1
    #-------------------------------#
    #      User-defined inputs      #
    #-------------------------------#
    params = Dict{String, Any}(
        "num_nodes" => 8,
        "instance" => i,
        "eigen_cuts_full" => true,
        "soc_linearized_cuts" => true
        )

    #-------------------------------------------#
    #      Optimization model and solution      #
    #-------------------------------------------#
    result = LOpt.run_LOpt_model(params, lopt_optimizer)

    result_all[i,1] = i
    result_all[i,2] = result["objective"]
    result_all[i,3] = result["solve_time"]
end


text:
using OrdinaryDiffEq, OpenQuantumTools
using Plots, StatsBase

H = DenseHamiltonian([(s)->1.0], [-σz], unit=:ħ)
u0 = PauliVec[1][1]

coupling = ConstantCouplings(["Z"], unit=:ħ)
# number of fluctuators
num = 10
# The values of b created here are in angular frequency units:
bvec = 0.01 * ones(num)
# γᵢ
γvec = log_uniform(0.01, 1, num)
# create the fluctuator coupling interaction
fluctuator_ensemble = EnsembleFluctuator(bvec, γvec);
interaction_fluctuator = Interaction(coupling, fluctuator_ensemble)
# create the Ohmic coupling interaction
ohmic_bath = Ohmic(1e-4, 4, 16)
interaction_ohmic = Interaction(coupling, ohmic_bath)
# merge these two bath objects into `InteractionSet`
interactions = InteractionSet(interaction_fluctuator, interaction_ohmic)
annealing = Annealing(H, u0, interactions=interactions)

tf = 100
prob = build_ensembles(annealing, tf, :ame, ω_hint = range(-2, 2, length=100))
t_list = range(0,tf,length=200)
sol = solve(prob, Tsit5(), EnsembleSerial(), trajectories=1000, reltol=1e-6, saveat=t_list)

dataset = []
st(s, so) = normalize(so(s, continuity=:left))
for s in t_list
    push!(dataset, [real(st(s, so)' * σx * st(s, so)) for so in sol])
end

pop_mean = []
pop_rmse = []
for data in dataset
    p_mean = sum(data)/1000
    p_rmse = sqrt(sum((x)->(x-p_mean)^2, data))/1000
    push!(pop_mean, p_mean)
    push!(pop_rmse, 2*p_rmse)
end

a_list = range(0,tf,length=400)
annealing_ame = Annealing(H, u0, coupling=coupling, bath=ohmic_bath)
sol_ame = solve_ame(annealing_ame, tf, alg=Tsit5(), ω_hint = range(-2, 2, length=100),
    reltol=1e-6, saveat=a_list)
ame_x = [real(tr(u*σx)) for u in sol_ame.u]

plot(a_list, ame_x, label="Ohmic", linewidth=2)
plot!(t_list, pop_mean, ribbon=pop_rmse, label="Hyrbrid", linewidth=2)
xlabel!("s (ns)")
ylabel!("<X>")

using HOQSTTutorials
HOQSTTutorials.tutorial_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])


text:
function ispositive(cumulant::Array{ComplexF64,4},nx::Int,ny::Int,Λ::Int)
    twopoint = reshape(cumulant,(2*ny-1)*(nx-Λ),(2*ny-1)*(nx-Λ))
    D = eigvals(twopoint)
    !any(x->x<0.0,D)
end

function positivity!(cumulant::Array{ComplexF64,4},nx::Int,ny::Int,Λ::Int)
    twopoint = reshape(cumulant,(2*ny-1)*(nx-Λ),(2*ny-1)*(nx-Λ))
    D,V = eigen(twopoint)
    # @info "Removing following eignvalues from second cumulant ", D[isless.(D,0)]
    @info "Largest (abs) eigenvalue removed from second cumulant = ", maximum(abs.(D[isless.(D,0)]))
    Dpos = max.(D,0.0)
    twopoint = V*diagm(Dpos)*inv(V)
    # optimise further:
    # mul!(twopoint,V,lmul!(diagm(Dpos),inv(V)))
    cumulant = reshape(twopoint,2*ny-1,nx-Λ,2*ny-1,nx-Λ)
end

function positivity!(cumulant::Array{ComplexF64,3},temp::Array{ComplexF64,4},nx::Int,ny::Int)
    temp .= zero(ComplexF64)
    @inbounds for m1=1:nx-1
        @inbounds for n1=-(ny-1):ny-1
            @inbounds for n2=-(ny-1):ny-1
                 temp[n2+ny,m1,n1+ny,m1] = cumulant[n2+ny,n1+ny,m1]
            end
        end
    end
    twopoint = reshape(temp,(2*ny-1)*(nx-1),(2*ny-1)*(nx-1))
    D,V = eigen(twopoint)
    D = real.(D)
    # @info D
    # @info "Removing following eignvalues from second cumulant ", D[isless.(D,0)]
    @info "Largest (abs) eigenvalue removed from second cumulant = ", maximum(abs.(D[isless.(D,0)]))
    Dpos = max.(D,0.0)
    twopoint = V*diagm(Dpos)*inv(V)
    # optimise further:
    # mul!(twopoint,V,lmul!(diagm(Dpos),inv(V)))
    temp = reshape(twopoint,2*ny-1,nx-1,2*ny-1,nx-1)
    @inbounds for m1=1:nx-1
        @inbounds for n1=-(ny-1):ny-1
            @inbounds for n2=-(ny-1):ny-1
                 cumulant[n2+ny,n1+ny,m1] = temp[n2+ny,m1,n1+ny,m1]
            end
        end
    end

end

function truncatecumulant!(d::AbstractDomain,cc::SecondCumulant{T},temp::SecondCumulant{T}) where T
    for m=1:d.nx-1
        D,V = eigen(cc[:,:,m])
        D[1:2d.nx-2] .= zero(T)
        # @show D
        # Dmax = [real(e) >= maximum(real.(D)) ? e : 0.0 for e in D]
        temp[:,:,m] = V*diagm(D)*inv(V)
        # @show Dmax
    end
    cc .= temp
    nothing
end


text:
"""
    KLocalOperator

A representation of a local operator, acting on certain `sites`, each with
`hilb_dims` dimension. The local operator is written in matrix form in this
basis as `mat`. For every row of `mat` there are several non-zero values contained
in `mel`, and to each of those, `to_change` contains the sites that must change the basis
new value, contained in `new_value`
"""
struct KLocalOperator{H<:AbstractHilbert,T,SV,Vel,Vti,Vtc,Vtv,OC} <: AbsLinearOp{T}
    hilb::H

    # list of sites on which this operator acts
    sites::SV

    # local hilbert space dimensions for each site
    hilb_dims::SV

    # Operator (as a matrix) represented
    mat::Matrix{T}

    # Matrix elements, per row
    mel::Vector{Vel}

    # List of indices that should be changed
    to_change::Vector{Vtc}

    # List of values of that should be changed
    new_values::Vector{Vtv}

    # Connections
    op_conns::Vector{OC}

    # List of state indices corresponding to mel
    # This is redundant, and equivalent to specifying to_change and new_values
    new_indices::Vector{Vti}
end

#KLocalOperatorRow(sites::AbstractVector, hilb_dims::AbstractVector, operator) = begin
#    hilb = DiscreteHilbert(hilb_dims)
#    KLocalOperatorRow(hilb, eltype(operator), sites, hilb_dims, operator)
#end

"""
    KLocalOperatorRow(sites::Vector, hilb_dims::Vector, operator)

Creates a KLocalOperator where connections are stored by row for the operator
    `operator` acting on `sites` each with `hilb_dims` local dimension.
"""
KLocalOperatorRow(hilb::AbstractHilbert, sites::AbstractVector, operator) =
    KLocalOperatorRow(eltype(operator), hilb, sites, shape(hilb)[sites], operator)

function KLocalOperatorRow(T::Type{<:Number}, hilb::AbstractHilbert, sites::AbstractVector, hilb_dims::AbstractVector, operator)
    # TODO Generalize to arbistrary spaces and not only uniform
    st  = state(T, hilb)
    st1 = state(T, hilb)

    mel         = Vector{Vector{T}}()
    new_indices = Vector{Vector{Int}}()
    to_change   = Vector{Vector{Vector{eltype(Int)}}}()
    new_values  = Vector{Vector{Vector{eltype(st)}}}()
    op_conns    = Vector{OpConnection{Vector{T},
                        Vector{eltype(Int)}, Vector{eltype(st)}}}()
    #SCT = StateChanges{eltype(to_change_els), eltype(new_values_els)}

    for (r, row) = enumerate(eachrow(operator))
        # create the containers for this row
        mel_els        = T[]
        mel_indices    = Int[]
        to_change_els  = Vector{Vector{Int}}()
        new_values_els = Vector{Vector{T}}()
        conns_els      = eltype(op_conns)()

        #if abs(row[r]) > 10e-6
            push!(mel_els, row[r])
            push!(mel_indices, r)
            push!(to_change_els, Int[])
            push!(new_values_els, eltype(st)[])
            push!(conns_els, (row[r], Int[], eltype(st)[]))
        #end

        set!(st, hilb, r)
        for (c, val) = enumerate(row)
            r == c && continue
            abs(val) < 10e-6 && continue

            set!(st1, hilb, c)
            cngd = Int[]
            nwvls = eltype(st)[]
            for (i, (v, vn)) = enumerate(zip(st, st1))
                if v != vn
                    push!(cngd, sites[i])
                    push!(nwvls, vn)
                end
            end
            push!(mel_els, val)
            push!(mel_indices, c)
            push!(to_change_els, cngd)
            push!(new_values_els, nwvls)
            push!(conns_els, (val, cngd, nwvls))
        end

        push!(mel, mel_els)
        push!(new_indices, mel_indices)
        push!(to_change, to_change_els)
        push!(new_values, new_values_els)
        push!(op_conns, conns_els)
    end
    KLocalOperator(hilb, sites, hilb_dims, convert(Matrix{T}, operator |> collect),
                        mel, to_change, new_values, op_conns, new_indices)
end

KLocalOperator(op::KLocalOperator, mat::AbstractMatrix) =
    KLocalOperatorRow(basis(op), copy(sites(op)), mat)

## Accessors
"""
    sites(op::KLocalOperator)

Returns the vector of `Int` labelling the lattice sites on which this operator
acts, in no particular order.
"""
sites(op::KLocalOperator) = op.sites

hilb_dims(op::KLocalOperator) = op.hilb_dims

QuantumOpticsBase.basis(op::KLocalOperator) = op.hilb

operators(op::KLocalOperator) = (op,)

densemat(op::KLocalOperator) = op.mat

conn_type(top::Type{KLocalOperator{H,T,SV,Vel,Vti,Vtc,Vtv,OC}}) where {H, T, SV, Vel, Vti, Vtc, Vtv, OC} =
    OpConnection{Vel, eltype(Vtc), eltype(Vtv)}
conn_type(op::KLocalOperator{H,T,SV,Vel,Vti,Vtc,Vtv,OC}) where {H, T, SV, Vel, Vti, Vtc, Vtv, OC} =
    OpConnection{Vel, eltype(Vtc), eltype(Vtv)}

eye_conn_type(top::Type{KLocalOperator{H,T,SV,Vel,Vti,Vtc,Vtv,OC}}) where {H, T, SV, Vel, Vti, Vtc, Vtv, OC} =
    OpConnectionIdentity{eltype(Vtc), eltype(Vtv)}

# Copy
duplicate(op::KLocalOperator) = copy(op)
Base.copy(op::KLocalOperator) =
    KLocalOperator(basis(op), deepcopy(op.sites), deepcopy(op.hilb_dims), deepcopy(op.mat),
                    deepcopy(op.mel), deepcopy(op.to_change),
                    deepcopy(op.new_values), deepcopy(op.op_conns),
                    deepcopy(op.new_indices))

##
function _row_valdiff!(conn::OpConnection, op::KLocalOperator, v::AState)
    @boundscheck checkbounds(v, nsites(basis(op)))

    # Find row index
    @inbounds r = local_index(v, basis(op), sites(op))

    append!(conn, op.op_conns[r])
end

function row_valdiff_index!(conn::OpConnectionIndex, op::KLocalOperator, v::AState)
    @boundscheck checkbounds(v, nsites(basis(op)))

    # Find row index
    r = local_index(v, basis(op), sites(op))

    # Return values
    mel = op.mel[r]
    ids = op.new_indices[r]

    append!(conn, (mel, ids))
end

function map_connections(fun::Function, op::KLocalOperator, v::AState)
    @boundscheck checkbounds(v, nsites(basis(op)))

    @inbounds r = local_index(v, basis(op), sites(op))
    for (mel, changes) = op.op_conns[r]
        fun(mel, changes, v)
    end
    return nothing
end

function accumulate_connections!(acc::AbstractAccumulator, op::KLocalOperator, v::AState)
    @boundscheck checkbounds(v, nsites(basis(op)))

    # If it is a doublestate, we are probably computing Operator x densitymatrix,
    # so we only iterate along the column of v
    if v isa ADoubleState
        @inbounds r = local_index(col(v), basis(op), sites(op))
    else
        @inbounds r = local_index(v, basis(op), sites(op))
    end

    for (mel,changes)=op.op_conns[r]
        acc(mel, changes, v)
    end

    return acc
end

Base.:(==)(l::KLocalOperator, r::KLocalOperator) =
    l.mat == r.mat && l.hilb_dims == r.hilb_dims

Base.isapprox(l::KLocalOperator, r::KLocalOperator; kwargs...) =
    l.hilb_dims == r.hilb_dims && isapprox(l.mat, r.mat; kwargs...)

Base.:-(op::KLocalOperator) = KLocalOperator(op, -op.mat)

_add_samesite(op_l::KLocalOperator, op_r::KLocalOperator) = _add_samesite!(duplicate(op_l), op_r)

function _add_samesite!(op_l::KLocalOperator, op_r::KLocalOperator)
    @assert op_l.sites == op_r.sites
    @assert length(op_l.mel) == length(op_r.mel)

    op_l.mat .+= op_r.mat

    n = length(op_l.mel)
    for (i, (mel_vals_l, mel_vals_r)) = enumerate(zip(op_l.mel, op_r.mel))

        for (mel, ids, tch, nv) = zip(op_r.mel[i], op_r.new_indices[i],
                                      op_r.to_change[i], op_r.new_values[i])
            found = 0
            # if mel_val_l is empty, then for sure there is nothing in here
            if !isempty(mel_vals_l)
                ids_tch = findall(isequal(tch), op_l.to_change[i])
                for id=ids_tch
                    if op_l.new_values[i][id] == nv
                        found = id
                        break
                    end
                end
            end

            if found == 0
                push!(op_l.mel[i], mel)
                push!(op_l.new_indices[i], ids)
                push!(op_l.to_change[i], tch)
                push!(op_l.new_values[i], tch)
            else
                op_l.mel[i][found] += mel
            end
        end
    end

    for (i, (connsᵢ_l , connsᵢ_r)) = enumerate(zip(op_l.op_conns, op_r.op_conns))
        add!(connsᵢ_l, connsᵢ_r)
    end
    return op_l
end

Base.transpose(op::KLocalOperator) =
    KLocalOperator(op, collect(transpose(op.mat)))

function Base.conj!(op::KLocalOperator)
    conj!(op.mat)
    map(conj!, op.op_conns)
    for el=op.mel
        conj!(el)
    end
    return op
end

Base.conj(op::KLocalOperator) = conj!(duplicate(op))
Base.adjoint(op::KLocalOperator) = conj(transpose(op))

Base.:*(a::Number, b::KLocalOperator) =
    _op_alpha_prod(b,a)
_op_alpha_prod(op::KLocalOperator, a::Number) =
    KLocalOperator(op, a*op.mat)

function Base.:*(opl::KLocalOperator, opr::KLocalOperator)
    @assert shape(basis(opl)) == shape(basis(opr))

    if sites(opl) == sites(opr)
        return KLocalOperator(opl, opl.mat * opr.mat)
    else
        disjoint = true
        for s=sites(opr)
            if s ∈ sites(opl)
                disjoint = false
                break
            end
        end
        if disjoint
            _kop_kop_disjoint_prod(opl,opr)
        else
            _kop_kop_joint_prod(opl, opr)
        end
    end
end

function _kop_kop_disjoint_prod(opl::KLocalOperator, opr::KLocalOperator)
    sl = sites(opl); sr = sites(opr)
    if length(sl) == 1 && length(sr) == 1
        # it's commutative
        if sl[1] > sr[1]
            _op =opl
            opl = opr
            opr = _op
        end
        sl = first(sites(opl))
        sr = first(sites(opr))

        # inverted also in QuantumOptics... who knows why
        mat = kron(opr.mat, opl.mat)
        return KLocalOperatorRow(basis(opl), [sl, sr], mat)
    else
        sites_new = sort(vcat(sites(opl), sites(opr)))
        ids_l = [findfirst(i .==sites_new) for i=sites(opl)]
        ids_r = [findfirst(i .==sites_new) for i=sites(opr)]
        throw("Tensor product between two operators on disjoint support but on
        more than 1 site each must still be implemneted. Sites in question are $(sites(opl)) and $(sites(opr)).")
    end
end

function _kop_kop_joint_prod(opl::KLocalOperator, opr::KLocalOperator)
    sl = sites(opl); sr = sites(opr)
    if length(sl) == 1 || length(sr) == 1
        reversed = false
        if length(sl) == 1
            _op = opl
            opl = opr
            opr = _op
            reversed = true
        end
        # opl has many dims, opr only 1
        sr = first(sites(opr))
        r_index = findfirst(sr .== sl)
        hdim_r = first(hilb_dims(opr))

        matrices = [Matrix(I, d, d) for d=hilb_dims(opl)]
        matrices[r_index] = opr.mat

        # inverted also in QuantumOptics... who knows why
        mat_r = kron(reverse(matrices)...)
        prod_mat = reversed ? mat_r*opl.mat : opl.mat*mat_r

        return KLocalOperator(opl, prod_mat)
    else
        sites_new = sort(vcat(sites(opl), sites(opr)))
        ids_l = [findfirst(i .==sites_new) for i=sites(opl)]
        ids_r = [findfirst(i .==sites_new) for i=sites(opr)]
        throw("Tensor product between two operators on overlapping support but
        both on more than 1 site must still be implemneted.")
    end
end

function permutesystems(a::AbstractMatrix, h_dims::Vector, perm::Vector{Int})
    #@assert length(a.basis_l.bases) == length(a.basis_r.bases) == length(perm)
    #@assert isperm(perm)
    data = reshape(a, [h_dims; h_dims]...)
    data = permutedims(data, [perm; perm .+ length(perm)])
    data = reshape(data, prod(h_dims), prod(h_dims))
    return data
end

Base.eltype(::T) where {T<:KLocalOperator} = eltype(T)
Base.eltype(::Type{KLocalOperator{H,T,SV,Vel,Vti,Vtc,Vtv,OC}}) where {H,T,SV,Vel,Vti,Vtc,Vtv,OC} =
    eltype(eltype(Vel))

Base.show(io::IO, op::KLocalOperator) = begin
    T    = eltype(op)
    s    = sites(op)
    h    = basis(op)
    dims = hilb_dims(op)
    mat  = densemat(op)

    print(io, "KLocalOperatorRow($T, $h, $s, $dims, $mat)")
end

Base.show(io::IO, m::MIME"text/plain", op::KLocalOperator) = begin
    T    = eltype(op)
    s    = sites(op)
    h    = basis(op)
    dims = hilb_dims(op)
    mat  = densemat(op)

    print(io, "KLocalOperator($T)\n  Hilb: $h\n  sites: $s  (size: $dims)\n")
    Base.print_array(io, mat)
end


text:
# load packages

using GeoPhyInv
using Statistics

# create simple (almost) homogeneous acoustic model

model=J.Gallery.Seismic(:acou_homo1)
J.Models.Seismic_addon!(model, randn_perc=0.01)

# a simple acquisition geometry

acqgeom = GeoPhyInv.Gallery.Geom(model.mgrid,:xwell);


text:
# taken from tutorial of Makie.jl

using Makie

function sactter_withcolors()
    x = [1, 2, 3, 4]
    y = [4, 3, 2, 1]
    colors = [0,2,4,6]
    scene = scatter(x, y, color = colors)
end

function scatter_increase_its_size()
    x = [1, 2, 3, 4]
    y = [4, 3, 2, 1]
    markersize = [1,3,2,4] ./ 10
    scene = scatter(x, y, markersize = markersize)
end

function scatter_mixin()
    x = [1, 2, 3, 4]
    y = [4, 3, 2, 1]
    markersize = [1,3,2,4] ./ 10
    colors = [0,2,4,6]
    scatter(x, y, color=colors, markersize=markersize)
end

function main()
    scatter_mixin()
end

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
import .FixedSizeArrays: @fixed_vector

abstract type AbstractDistanceField end
abstract type AbstractUnsignedDistanceField <: AbstractDistanceField end
abstract type AbstractSignedDistanceField <: AbstractDistanceField end

"""
Abstract to categorize geometry primitives of dimensionality `N` and
the numeric element type `T`.
"""
abstract type AbstractGeometry{N, T} end
abstract type AbstractMesh{VertT, FaceT}  end # <: AbstractGeometry
abstract type GeometryPrimitive{N, T} <: AbstractGeometry{N, T} end

"""
Abstract to classify Simplices. The convention for N starts at 1, which means
a Simplex has 1 point. A 2-simplex has 2 points, and so forth. This convention
is not the same as most mathematical texts.
"""
abstract type AbstractSimplex{S, T} <: StaticVector{S, T} end

"""
A `Simplex` is a generalization of an N-dimensional tetrahedra and can be thought
of as a minimal convex set containing the specified points.

* A 0-simplex is a point.
* A 1-simplex is a line segment.
* A 2-simplex is a triangle.
* A 3-simplex is a tetrahedron.

Note that this datatype is offset by one compared to the traditional
mathematical terminology. So a one-simplex is represented as `Simplex{2,T}`.
This is for a simpler implementation.

It applies to infinite dimensions. The structure of this type is designed
to allow embedding in higher-order spaces by parameterizing on `T`.
"""
struct Simplex{S, T} <: AbstractSimplex{S, T}
    data::NTuple{S, T}
    Simplex{S, T}(x::NTuple{S, T}) where {S, T} = new{S, T}(x)
    Simplex{S, T}(x::NTuple{S}) where {S, T} = new{S, T}(StaticArrays.convert_ntuple(T, x))
end

@fixed_vector Vec StaticVector

@fixed_vector Point StaticVector

@fixed_vector Normal StaticVector

@fixed_vector TextureCoordinate StaticVector

@fixed_vector Face StaticVector

"""
A Face is typically used when constructing subtypes of `AbstractMesh` where
the `Face` should not reproduce the vertices, but rather index into them.
Face is parameterized by:

* `N` - The number of vertices in the face.
* `T` - The type of the indices in the face, a subtype of Integer.

"""
Face

"""
OffsetInteger type mainly for indexing.
* `O` - The offset relative to Julia arrays. This helps reduce copying when
communicating with 0-indexed systems such ad OpenGL.
"""
struct OffsetInteger{O, T <: Integer} <: Integer
    i::T

    OffsetInteger{O, T}(x::Integer) where {O, T <: Integer} = new{O, T}(T(x + O))
end

raw(x::OffsetInteger) = x.i
raw(x::Integer) = x

"""
A `HyperRectangle` is a generalization of a rectangle into N-dimensions.
Formally it is the cartesian product of intervals, which is represented by the
`origin` and `width` fields, whose indices correspond to each of the `N` axes.
"""
struct HyperRectangle{N, T} <: GeometryPrimitive{N, T}
    origin::Vec{N, T}
    widths::Vec{N, T}
end

struct HyperCube{N, T} <: GeometryPrimitive{N, T}
    origin::Vec{N, T}
    width::T
end

"""
A `HyperSphere` is a generalization of a sphere into N-dimensions.
A `center` and radius, `r`, must be specified.
"""
struct HyperSphere{N, T} <: GeometryPrimitive{N, T}
    center::Point{N, T}
    r::T
end

struct SimpleRectangle{T} <: GeometryPrimitive{2, T}
    x::T
    y::T
    w::T
    h::T
end

# TODO remove before 0.2.0 tag
const Rectangle = SimpleRectangle

"""
A rectangle in 3D space.
"""
struct Quad{T} <: GeometryPrimitive{3, T}
    downleft::Vec{3, T}
    width   ::Vec{3, T}
    height  ::Vec{3, T}
end

struct Pyramid{T} <: GeometryPrimitive{3, T}
    middle::Point{3, T}
    length::T
    width ::T
end

struct Particle{N, T} <: GeometryPrimitive{N, T}
    position::Point{N, T}
    velocity::Vec{N, T}
end

"""
A `SignedDistanceField` is a uniform sampling of an implicit function.
The `bounds` field corresponds to the sampling space intervals on each axis.
The `data` field represents the value at each point whose exact location
can be rationalized from `bounds`.
The type is parameterized by:

* `N` - The dimensionality of the sampling space.
* `SpaceT` - the type of the space where we will uniformly sample.
* `FieldT` - the type resulting from evaluation of the implicit function.

Note that decoupling the space and field types is useful since geometry can
be formulated with integers and distances can be measured with floating points.
"""
mutable struct SignedDistanceField{N,SpaceT,FieldT} <: AbstractSignedDistanceField
    bounds::HyperRectangle{N,SpaceT}
    data::Array{FieldT,N}
end

"""
The `HomogenousMesh` type describes a polygonal mesh that is useful for
computation on the CPU or on the GPU.
All vectors must have the same length or must be empty, besides the face vector
Type can be Void or a value, this way we can create many combinations from this
one mesh type.
This is not perfect, but helps to reduce a type explosion (imagine defining
every attribute combination as a new type).
"""
struct HomogenousMesh{VertT, FaceT, NormalT, TexCoordT, ColorT, AttribT, AttribIDT} <: AbstractMesh{VertT, FaceT}
    vertices            ::Vector{VertT}
    faces               ::Vector{FaceT}
    normals             ::Vector{NormalT}
    texturecoordinates  ::Vector{TexCoordT}
    color               ::ColorT
    attributes          ::AttribT
    attribute_id        ::Vector{AttribIDT}
end

function HomogenousMesh(
        vertices::AbstractVector{VertT},
        faces::AbstractVector{FaceT},
        normals::AbstractVector{NormalT},
        texturecoordinates::AbstractVector{TexCoordT},
        color::ColorT, attributes::AttribT,
        attribute_id::AbstractVector{AttribIDT}
    ) where {VertT, FaceT, NormalT, TexCoordT, ColorT, AttribT, AttribIDT}

    HomogenousMesh{VertT, FaceT, NormalT, TexCoordT, ColorT, AttribT, AttribIDT}(
        convert(Vector{VertT}, vertices),
        convert(Vector{FaceT}, faces),
        convert(Vector{NormalT}, normals),
        convert(Vector{TexCoordT}, texturecoordinates),
        color,
        attributes,
        convert(Vector{AttribIDT}, attribute_id),
    )
end
"""
AbstractFlexibleGeometry{T}

AbstractFlexibleGeometry refers to shapes, which are somewhat mutable.
"""
abstract type AbstractFlexibleGeometry{T} end
const AFG = AbstractFlexibleGeometry

"""
FlexibleConvexHull{T}

Represents the convex hull of finitely many points. The number of points is not fixed.
"""
struct FlexibleConvexHull{T} <: AFG{T}
    vertices::Vector{T}
end

"""
FlexibleSimplex{T}

Represents a Simplex whos number of vertices is not fixed.
"""
struct FlexibleSimplex{T} <: AFG{T}
    vertices::Vector{T}
end

"""
A `Cylinder` is a 2D rectangle or a 3D cylinder defined by its origin point,
its extremity and a radius. `origin`, `extremity` and `r`, must be specified.
"""
struct Cylinder{N,T<: AbstractFloat} <: GeometryPrimitive{N,T}
    origin::Point{N,T}
    extremity::Point{N,T}
    r::T
end

"""
AbstractConvexHull

Groups all geometry types, that can be described as the convex hull of finitely
many points.
"""
const AbstractConvexHull = Union{
    Simplex, FlexibleConvexHull, FlexibleSimplex,
    HyperCube, HyperRectangle
} # should we parametrize ACH by the type of points T?


text:
#using DelimitedFiles
#using StaticArrays

ssCO = readdlm((@__DIR__)*"/CO_shielding_functions/shield.03.5.69-557-36.dat"); #CO linewidth = 0.3, T_ex = 5
const Nbin_nco = 46
const Nbin_nh2 = 41
const Nco = SVector{Nbin_nco,Float64}(ssCO[10:55])
const Nh2 = SVector{Nbin_nh2,Float64}(ssCO[58:98])
const istart = 100 #start line number for 12CO

function get_fco(is)
    a=[]
    for i in 0:4
        a = vcat(a, ssCO[is+i,:])
    end
    a = a[1:Nbin_nco]
end

ssCOtable = zeros(Nbin_nco, Nbin_nh2)
for j in 1:Nbin_nh2
    jj = istart+(j-1)*5
    ssCOtable[:, j] .= get_fco(jj)
end

#=
fCOss = interpolate((Nco, Nh2), ssCOtable, Gridded(Linear()));
function fCOselfshield(Nco,Nh2)
    Nh2 = Nh2 >= 1e23 ? 1e23 : Nh2
    Nco = Nco >= 1e19 ? 1e19 : Nco
    Nh2 = Nh2 <= 1e10 ? 1e10 : Nh2
    Nco = Nco <= 1e10 ? 1e10 : Nco
    return fCOss(Nco, Nh2)
end
=#

include("interpolation.jl")
const par_fCOss = InterpRange{Float64}(10., 19., 1/0.2, 15., 23., 1/0.2)
#const tblCO = SMatrix{Nbin_nco,Nbin_nh2}(ssCOtable) #size of SArray should not be larger than 100
const tblCO = Matrix{Float64}(ssCOtable)
function fCOselfshield(NCO,NH2)
    x, y = log10(NCO), log10(NH2)
    interp(x,y, tblCO, par_fCOss)
end

function fH2selfshield(Nh2)
    x = Nh2 * 2e-15
    b5 = 2.0
    return 0.965 / (1.0 + x/b5)^2 + 0.035 / sqrt(1 + x) * exp(-8.5e-4 * sqrt(1 + x))
end

function fCselfshield(NC,NH2)
    rH2 = 2.8e-22 * NH2
    return exp( -1.6e-17 * NC ) * exp(-rH2) / (1.0 + rH2)
end


text:
# function that fulfills all incidents in the exchange queue during
# the current minute "mnt" if possible
function fulfill_exchange_queue!(incidents::DataFrame,
                                 sim_data::DataFrame,
                                 ressource_flow::Array{Int64,3},
                                 drivingtime::Array{Float64,2},
                                 exchange_queue::Array{Union{Missing,Int64},1},
                                 candidates::Vector{Union{Missing,Float64}},
                                 locations::Vector{Int64},
                                 exchange_reserve::Vector{Int64},
                                 card_districts::Int64,
                                 mnt::Int64,
                                 dispatch::String,
                                 max_drive::Float64,
                                 p::Int64)
    for j = 1:length(exchange_queue)
        # check whether the queue holds candidates for exchange
        if ismissing(exchange_queue[j]) == false
            # clear the list of current candidates in the queue
            current_case = exchange_queue[j]
            candidate = false
            candidates .= missing
            # check for all districts within an acceptable perimeter
            # how long the driving time would take during the current
            # traffic flow
            @simd for i = 1:card_districts
                if ressource_flow[mnt,i,1] > exchange_reserve[p] &&
                    drivingtime[locations[i],incidents[current_case,:location]] <= max_drive
                    drive = dispatch_drivingtime(drivingtime[locations[i], incidents[current_case,:location]], traffic, mnt)
                    # save the location as exchange candidate and weight it by
                    # the current number of cars available for dispatch
                    if drive < max_drive
                        candidates[i] = ressource_flow[mnt,i,1]/drive
                        candidate = true
                    end
                end
            end
            # determine the best candidate from the candidate list (if available)
            if candidate == true
                i = findmin(skipmissing(candidates))[2]
                k = 1
                fastest_time = 0
                # dispatch one car from the exchange location
                incident_dispatch!(sim_data::DataFrame,
                                   incidents::DataFrame,
                                   ressource_flow::Array{Int64,3},
                                   drivingtime::Array{Float64,2},
                                   traffic::Array{Float64,2},
                                   locations::Vector{Int64},
                                   current_case::Int64,
                                   mnt::Int64,
                                   i::Int64,
                                   k::Int64,
                                   dispatch::String,
                                   fastest_time::Int64)
            end
        else
            break
        end
    end
end


text:
using Bench
using Test

@testset "Bench.jl" begin
    # Write your own tests here.
end


text:
# Magic Square Generator
# ======================
#
# # Reference
#
# - Moler, Cleve B. Experiments with MATLAB. Society for Industrial and Applied Mathematics, 2011.

# The JuliaLang community is famous for naming their functions as something like
# `verylongfunctionnamethatistotallyunreadable` and abusing unicode characters all over the source.
# A more appropriate name for this is probably `is🎩`, but let's stick with the old style for now.
function is_magic(m::AbstractMatrix{<:Integer})
    n_r, n_c = size(m)
    # Is square?
    if n_r != n_c
        return false
    end

    n = n_r
    magic_const = (1 + n ^ 2) * n ÷ 2

    # Are all rows equal to magic constant?
    if !all(sum(m, dims=2) .== magic_const)
        return false
    end
    # Are all columns equal to magic constant?
    if !all(sum(m, dims=1) .== magic_const)
        return false
    end
    # And the two diagonals?
    if sum([m[i, i] for i in 1:n]) != magic_const
        return false
    end
    if sum([m[p...] for p in zip(1:n, n:-1:1)]) != magic_const
        return false
    end

    true
end

# Odd ordered magic square generator
#
# It is easy to see the correctness of this algorithm. For example, compare the elements in each
# column with the one to the northeast of itself, there are `n - 1` elements that are 1 less while
# exactly 1 that is `n - 1` larger, which makes the two columns add up to the same value.

"Generates a magic square of an odd order `n` using Siamese method (De la Loubère method)."
function odd_magic(n::Integer)
    @assert isodd(n)

    m = zeros(typeof(n), n, n)
    r, c = 1, div(n + 1, 2)

    for i in 1 : n*n
        m[r, c] = i

        # No need to do bounds checking, as only the last one will fall out of the board.
        p = r + 1, c
        r, c = mod1(r - 1, n), mod1(c + 1, n)
        if m[r, c] != 0
            r, c = p
        end
    end

    m
end

# A vectorized version of `odd_magic`.
function odd_magic_vec(n::Integer)
    @assert isodd(n)
    i = 1:n
    offset = div(n - 3, 2)
    A = @. mod(i + i' + offset, n)
    B = @. mod(i + i' * 2 - 2, n)
    @. n * A + B + 1
end

# TODO: Even ordered magic square generator
function singly_even_magic(n::Integer)
    error("unimplemented")
end

function doubly_even_magic(n::Integer)
    error("unimplemented")
end

function magic(n::Integer)
    if isodd(n)
        odd_magic(n)
    else
        half_n = div(n, 2)
        if isodd(half_n)
            singly_even_magic(n)
        else
            doubly_even_magic(n)
        end
    end
end

# Unit tests
using Test

@testset "Odd magic square tests" begin
    orders = 1:2:21
    for n in orders
        m = odd_magic(n)
        @test m == odd_magic_vec(n)
        @test is_magic(m)
    end
end

max_stars_repo_name: gyk/TrivialSolutions

max_issues_repo_name: gyk/TrivialSolutions

max_forks_repo_name: gyk/TrivialSolutions


text:
# Laplace distribution

import StatsFuns
export Laplace

@parameterized Laplace(μ,σ) ≪ (1/2) * Lebesgue(ℝ)

@kwstruct Laplace()

function logdensity(d::Laplace{()} , x)
    return -abs(x)
end

Base.rand(rng::AbstractRNG, μ::Laplace{()}) = rand(rng, Dists.Laplace())

≪(::Laplace, ::Lebesgue{X}) where X <: Real = true

TV.as(::Laplace) = asℝ

@μσ_methods Laplace()
@half Laplace()


text:
module ISPC

include("ast_tools.jl")
include("toolchain.jl")
export load_ispc, ispc_native, ispc_llvm
include("runtime.jl")
include("basefuncs.jl")
include("codegen.jl")
include("transform.jl")
include("extract.jl")
include("macros.jl")
export @foreach, @unmasked, @coherent, @ispc, @kernel

end # module


text:
module ConfidenceWeighted

export SoftConfidenceWeighted

include("SoftConfidenceWeighted.jl")

end # module


text:
module Graph6

using SimpleTraits
using Graphs
using Graphs: AbstractGraphFormat

import Graphs: loadgraph, loadgraphs, savegraph

export Graph6Format

struct Graph6Format <: AbstractGraphFormat end

function _bv2int(x::BitVector)
  @assert(length(x) <= 8 * sizeof(Int))
  acc = 0
  for i = 1:length(x)
    acc = acc << 1 + x[i]
  end
  return acc
end

function _int2bv(n::Int, k::Int)
  bitstr = lstrip(bitstring(n), '0')
  l = length(bitstr)
  padding = k - l
  bv = falses(k)
  for i = 1:l
    bv[padding + i] = (bitstr[i] == '1')
  end
  return bv
end

function _g6_R(_x::BitVector)::Vector{UInt8}
  k = length(_x)
  padding = cld(k, 6) * 6 - k
  x = vcat(_x, falses(padding))
  nbytes = div(length(x), 6)
  bytevec = Vector{UInt8}(undef, nbytes)   # uninitialized data!
  for i = 1:nbytes
    xslice  = x[((i - 1) * 6 + 1):(i * 6)]

    intslice = 0
    for bit in xslice
      intslice = (intslice << 1) + bit
    end
    intslice += 63
    bytevec[i] = intslice
  end
  return UInt8.(bytevec)
end

_g6_R(n::Int, k::Int) = _g6_R(_int2bv(n, k))

function _g6_Rp(bytevec::Vector{UInt8})
  nbytes = length(bytevec)
  x = BitVector()
  for byte in bytevec
    bits = _int2bv(byte - 63, 6)
    x = vcat(x, bits)
  end
  return x
end

function _g6_N(x::Integer)::Vector{UInt8}
  if (x < 0) || (x > 68719476735) error("x must satisfy 0 <= x <= 68719476735")
  elseif (x <= 62) nvec = [x + 63]
  elseif (x <= 258047)
    nvec = vcat([0x7e], _g6_R(x, 18))
  else
    nvec = vcat([0x7e; 0x7e], _g6_R(x, 36))
  end
  return UInt8.(nvec)
end

function _g6_Np(N::Vector{UInt8})
  if N[1] < 0x7e return (Int(N[1] - 63), N[2:end])
  elseif N[2] < 0x7e return (_bv2int(_g6_Rp(N[2:4])), N[5:end])
  else return(_bv2int(_g6_Rp(N[3:8])), N[9:end])
  end
end

"""
    _graphToG6String(g)

Given a graph `g`, create the corresponding Graph6 string.
"""
function _graphToG6String(g::Graphs.Graph)
  A = adjacency_matrix(g, Bool)
  n = nv(g)
  nbits = div(n * (n - 1), 2)
  x  = BitVector(undef, nbits)

  ind = 0
  for col = 2:n, row = 1:(col - 1)
    ind += 1
    x[ind] = A[row, col]
  end
  return join([">>graph6<<", String(_g6_N(n)), String(_g6_R(x))])
end

function _g6StringToGraph(s::AbstractString)
  if startswith(s, ">>graph6<<")
    s = s[11:end]
  end
  V = Vector{UInt8}(s)
  (nv, rest) = _g6_Np(V)
  bitvec = _g6_Rp(rest)

  g = Graphs.Graph(nv)
  n = 0
  for i in 2:nv, j in 1:(i - 1)
    n += 1
    if bitvec[n]
      add_edge!(g, j, i)
    end
  end
  return g
end

function loadgraph6_mult(io::IO)
  n = 0
  graphdict = Dict{String,Graphs.Graph}()
  while !eof(io)
    n += 1
    line = strip(chomp(readline(io)))
    gname = "graph$n"
    if length(line) > 0
        g = _g6StringToGraph(line)
        graphdict[gname] = g
    end
  end
  return graphdict
end

"""
    loadgraph6(io, gname="graph")

Read a graph from IO stream `io` in the [Graph6](http://users.cecs.anu.edu.au/%7Ebdm/data/formats.txt)
format. Return the graph.
"""
loadgraph6(io::IO, gname::String="graph") = loadgraph6_mult(io)[gname]

"""
    savegraph6(io, g, gname="graph")

Write a graph `g` to IO stream `io` in the [Graph6](http://users.cecs.anu.edu.au/%7Ebdm/data/formats.txt)
format. Return 1 (number of graphs written).
"""
function savegraph6 end
@traitfn function savegraph6(io::IO, g::::(!Graphs.IsDirected), gname::String = "graph")
  str = _graphToG6String(g)
  println(io, str)
  return 1
end

function savegraph6_mult(io::IO, graphs::Dict)
  ng = 0
  sortkeys = sort(collect(keys(graphs)))
  for gname in sortkeys
    ng += savegraph6(io, graphs[gname], gname)
  end
  return ng
end

loadgraph(io::IO, gname::String, ::Graph6Format) = loadgraph6(io, gname)
loadgraphs(io::IO, ::Graph6Format) = loadgraph6_mult(io)
savegraph(io::IO, g::AbstractGraph, gname::String, ::Graph6Format) = savegraph6(io, g, gname)
savegraph(io::IO, d::Dict, ::Graph6Format) = savegraph6_mult(io, d)

end # module


text:
__precompile__()

module GAFramework

using Random
using FileIO

export GAModel, ga!,
fitness, genauxga, crossover!, mutation!, selection, randcreature,
printfitness, savecreature,
GAState, GAIterable, loadgastate, savegastate,
RouletteWheelSelection, TournamentSelection

"""
To create a GA with a specific GAModel, import this module,
make a GAModel with the following interface functions:
    fitness (has default)
    genauxga  (has default)
    crossover! (no default)
    mutation! (has identity function as default)
    selection (has default)
    randcreature (no default)
    printfitness (has default)
    savecreature (has default)
    stopcondition (has default)
 """
abstract type GAModel end

const DEFAULT_GASTATE_PARAMS = Dict(:mutation_rate=>0.1,
    :print_fitness_iter=>1,
    :save_creature_iter=>0,
    :save_state_iter=>0,
    :file_name_prefix=>"gamodel")

include("ga.jl")

"""
Fitness function.
    fitness(x) is maximized always
    To minimize x.objvalue, dispatch fitness(x) to -x.objvalue for your Creature
    Recommended to make this either x.objvalue to maximize or -x.objvalue to minimize

    Since fitness(x) used for selecting the fittest creature, elites, and parents,
    all the computationally expensive part of calculating the fitness value should
    be implemented in the randcreature method.
"""
fitness(x) = x.objvalue

"""
    genauxga(model::GAModel) :: GAModel auxiliary structure

Given model GAM <: GAModel, generate auxiliary scratch space for calculating fitness scores
    model = GAM(G1,G2)
    aux = genauxga(model)
The purpose is to not allocate memory every time you calculate fitness for a new creature.
"""
genauxga(model::GAModel) = nothing

"""
    crossover!(z, x,y, model::GAModel, aux, rng)

Crosses over x and y to create a child. Optionally use space in z as a
scratch space or to create the child. aux is more scratch space. rng is random number generator.
    model = GAM(G1,G2)
    aux = genauxga(model)
    x = randcreature(model,aux)
    y = randcreature(model,aux)
    z = randcreature(model,aux)
    child = crossover(z,x,y,model,aux,rng)
"""
crossover!(z, x, y, st::GAState, aux, rng::AbstractRNG) =
    error("crossover not implemented for $(typeof(z)) and $(typeof(st.model))")

"""
    Mutates a incoming creature and outputs mutated creature
"""
mutation!(creature, st::GAState, aux, rng::AbstractRNG) = creature

"""
    selection(pop::Vector, n::Integer, rng)

    Generate a vector of n tuples (i,j) where i and j are
    indices into pop, and where pop[i] and pop[j] are the
    selected parents.
    Uses binary tournament selection by default.
"""
selection(pop::Vector, n::Integer, st::GAState, rng::AbstractRNG) =
    selection(TournamentSelection(2), pop, n, st, rng)

"""
    randcreature(model::GAModel, aux)

    Create a random instance of a creature, given a GAModel.
    There is always a creature associated with a GAModel
"""
randcreature(model::GAModel, aux, rng::AbstractRNG) =
    error("randcreature not implemented for $(typeof(model))")

"""
Logging
    * saves state every save_state_iter iterations to file
        - restart using state = loadgastate(filename) & ga!(state)
    * outputs creature every save_creature_iter iterations to file
    * prints fitness value every print_fitness_iter iterations to screen

    print the fitness of fittest creature every n iteration
        print_fitness_iter::Int
    save the fittest creature to file every n iteration
        save_creature_iter::Int
    save the entire state of the GA (i.e. this struct) to file every n iteration
        save_state_iter::Int
    prefix for the files to be save
        file_name_prefix::AbstractString
"""
printfitness(curgen::Int, x::Any) =
    println("curgen: $(curgen) fitness: $(fitness(x))")
savecreature(file_name_prefix::AbstractString, curgen::Int, x) =
    save("$(file_name_prefix)_creature_$(curgen).bson", "creature", x)

stopcondition(st::GAState) = st.curgen > st.ngen

"""
    Saves ga state to file
"""
function savegastate(file_name_prefix::AbstractString, curgen::Integer, state::GAState)
    filename = "$(file_name_prefix)_state_$(curgen).bson"
    println("Saving state to file $filename")
    save(filename,"state", state)
end

function loadgastate(filename::AbstractString)
    println("Load state from file $filename")
    load(filename, "state")
end

include("selections.jl")

include("models/coordinatega.jl")
include("models/permga.jl")
include("models/netalga.jl")
include("models/magnaga.jl")

function logiteration(st::GAState)
    creature = st.pop[1]
    gp = st.params
    file_name_prefix = gp[:file_name_prefix]
    if gp[:print_fitness_iter] > 0 && mod(st.curgen, gp[:print_fitness_iter]) == 0
        printfitness(st.curgen, creature)
    end
    if gp[:save_creature_iter] > 0 && mod(st.curgen, gp[:save_creature_iter]) == 0
        savecreature(file_name_prefix, st.curgen, creature)
    end
    if gp[:save_state_iter] > 0 && mod(st.curgen, gp[:save_state_iter]) == 0
        filename = "$(file_name_prefix)_state_$(st.curgen).bson"
        save(filename, "state", st)
    end
end

end


text:
@doc raw"""
    g = first_order_system(K, τ)

construct a first-order transfer function with gain `K` and time constant `τ`:

$$g(s)=\frac{K}{\tau s+1}$$

# example
```jldoctest
K = 1.0
τ = 3.0
g = first_order_system(K, τ)
# output
    1.0
-----------
3.0*s + 1.0
```

# returns
* `g::TransferFunction`: the first order transfer function. well, (0, 1) order.
"""
first_order_system(K::Float64, τ::Float64) = K / (τ * s + 1)

@doc raw"""
    g = second_order_system(K, τ, ξ)

construct a second-order transfer function with gain `K`, time constant `τ`, and damping coefficient `ξ`:

$$g(s)=\frac{K}{\tau^2 s^2 + 2\tau \xi s +1}$$

# example
```jldoctest
K = 1.0
τ = 2.0
ξ = 0.1
g = second_order_system(K, τ, ξ)
# output
         1.0
---------------------
4.0*s^2 + 0.4*s + 1.0
```

# returns
* `g::TransferFunction`: the second order transfer function. well, (0, 2) order.
"""
second_order_system(K::Float64, τ::Float64, ξ::Float64) = K / (τ^2 * s^2 + 2 * τ * ξ * s + 1)

@doc raw"""
    τ = time_constant(g)

compute the time constant τ of an order (0, 1) or order (0, 2) transfer function.

order (0, 1) representation:

$$g(s)=\frac{K}{\tau s+1}$$

order (0, 2) representation:

$$g(s)=\frac{K}{\tau^2 s^2 + 2\tau \xi s +1}$$

# returns
`τ::Float64`: the time constant.

# examples

```jldoctest
g = 4 / (6 * s + 2)
time_constant(g)
# output
3.0

g = 1.0 / (8 * s^2 + 0.8 * s + 2)
time_constant(g)
# output
2.0
```
"""
function time_constant(g::TransferFunction)
    if system_order(g) == (0, 1)
        # xx / (a s + b) ==> τ = a / b
        return g.denominator[1] / g.denominator[0]
    elseif system_order(g) == (0, 2)
        # xx / (a s^2 + b s + c) ==> τ = sqrt(a/c)
        return sqrt(g.denominator[2] / g.denominator[0])
    else
        error("`time_constant` only supports order (0, 1) and (0, 2) transfer functions.")
    end
end

@doc raw"""
    ξ = damping_coefficient(g)

compute the damping coefficient ξ of an order (0, 2) transfer function.

order (0, 2) representation:

$$g(s)=\frac{K}{\tau^2 s^2 + 2\tau \xi s +1}$$

# returns
`ξ::Float64`: the damping coefficient

# examples
```jldoctest
g = 1.0 / (8 * s^2 + 0.8 * s + 2)
damping_coefficient(g)
# output
0.1
```
"""
function damping_coefficient(g::TransferFunction)
    if system_order(g) != (0, 2)
        error("`damping_coefficient` only pertains to (0, 2) transfer functions.")
    end
    τ = time_constant(g)
    # xx / (a s² + b s + c) then 2τξ = b / c
    τξ2 = g.denominator[1] / g.denominator[0]
    return τξ2 / (2 * τ)
end

max_stars_repo_name: SimonEnsemble/Controlz

max_issues_repo_name: SimonEnsemble/Controlz

max_forks_repo_name: SimonEnsemble/Controlz


text:
############################################################################
#
#  main.jl
#
#  Purpose:    Main routine for Quadrotor simulation.
#
#  Author:     Edward Daughtery
#
#  Date:       05 August 2016
#
#
############################################################################
using PyPlot
using Aerospace
include("Quad_Types.jl")
using Quad_Types
include("Software.jl")
include("Quadrotor.jl")
include("Quad_Forces.jl")

# Closing all Figures
CloseAll();

# Simulation Control
time = 0.0;
dt = 0.001;
tfinal = 10.0;

## Number of Sample points
nn = (tfinal - time)/dt + 1;
NN = UInt(nn);

TT = zeros(NN);
ALT = zeros(NN);
VMAG = zeros(NN);
GAM = zeros(NN);
AZ = zeros(NN);
ROLL = zeros(NN);
PITCH = zeros(NN);
HEADING = zeros(NN);
W1 = zeros(NN);
W2 = zeros(NN);
W3 = zeros(NN);
W4 = zeros(NN);

## Position Initialization
latitude = 0.0*d2r;
longitude = 0.0*d2r;
altitude = 0.0*d2r;
position = InitPos( time, latitude, longitude, altitude );

## Velocity Initialization
vmag = 0.0;
gamma = 0.0*d2r;
azimuth = 0.0*d2r;
velocity = InitVel( vmag, gamma, azimuth, position );

# Attitude Initialization
roll = 0.0*d2r;
pitch = 0.0*d2r;
heading = 0.0*d2r;
euler = InitEuler(roll, pitch, heading);

# Angle Rate Intialization
wb_bn_body = [0.0;0.0;0.0];
angle_rates = InitAngleRates(wb_bn_body, position, euler);

# Quadrotor Initialization
quad1 = QuadInit( time, dt, position, velocity, euler, angle_rates );

for i = 1:NN
    TT[i] = time;
    ROLL[i] = quad1.airframe.euler.roll*r2d;
    PITCH[i] = quad1.airframe.euler.pitch*r2d;
    HEADING[i] = quad1.airframe.euler.heading*r2d;
    ALT[i] = quad1.airframe.position.altitude;
    W1[i] = quad1.motors.motor1.position;
    W2[i] = quad1.motors.motor2.position;
    W3[i] = quad1.motors.motor3.position;
    W4[i] = quad1.motors.motor4.position;

    quad1 = QuadUpdate(quad1);
    time = time + dt;

end

figure()
plot(TT, ALT)
grid
xlabel("Time (s)")
ylabel("Altitude (m)")

figure()
plot(TT, ROLL)
grid
xlabel("Time (s)")
ylabel("Roll (deg)")

figure()
plot(TT, PITCH)
grid
xlabel("Time (s)")
ylabel("Pitch (deg)")

figure()
plot(TT, HEADING)
grid
xlabel("Time (s)")
ylabel("Heading (deg)")

figure()
plot(TT, W1, TT, W2, TT, W3, TT, W4)
grid
xlabel("Time (s)")
ylabel("Motor Speed (rad/s)")
legend(["Motor 1", "Motor 2", "Motor 3", "Motor 4"]);


text:
module MochaClassificationWrapper

using Mocha
using Compat

export MWrap, forwardbackward,evaluateTestNN, initMochaBackend, getpredprob

"""
`initMochaBackend(use_gpu::Bool)` returns a Mocha Backend for use in MochaClassification and MochaRegression.
"""
function initMochaBackend(use_gpu::Bool)
  if use_gpu == false
    backend = CPUBackend()
  else
    backend = GPUBackend()
  end
  init(backend)
  backend
end

"""
Internal wrapper type for interaction with Mocha.
"""
type MWrap
    net          ::Net
    nparams      ::Int64
    paramRanges  ::Array{Int64,2}
    layersWParam ::Array{Int64,1}
    # Constructor
    function MWrap(
      data_layer::MemoryDataLayer,
      barebonesNN::Function,
      name::AbstractString,
      accuracylayer::Bool,
      predproblayer::Bool,
      bck::Backend)

        (common_layers,pred_prob_layer, storage_layer, loss_layer,acc_layer) = barebonesNN()

        if accuracylayer
            net = Net(name,bck,[data_layer, common_layers..., loss_layer, acc_layer])
        elseif predproblayer
            net = Net(name,bck,[data_layer,common_layers...,pred_prob_layer,storage_layer])
        else
            net = Net(name,bck,[data_layer, common_layers..., loss_layer])
        end
        init(net)
        (npara,paramRanges,layersWParam) = parameterisation(net)
        new(net,npara,paramRanges,layersWParam)
    end
end

"""
Calculates negative loss and gradient over minibatch
"""
function forwardbackward(mw::MWrap,para::Array{Float64,1}; regu_coef = 0.)

    mw = setparams!(mw,para)
    val = forward(mw.net, regu_coef)
    backward(mw.net, regu_coef)
    grad = getgrad(mw)
    return (-val, -grad)
end

"""
Evaluate loss over the whole dataset (normalized by batchsize)
"""
function forward_epoch_loss(net::Net)
    ll = 0.;
    epoch = get_epoch(net)
    count_batches = 0;
    while get_epoch(net) == epoch
        count_batches += 1
        result_forward = forward(net, 0.0)
        ll -= result_forward
    end
    return ll
end

"""
Evaluate outputs (eg predictive probabilities) over whole dataset
"""
function forward_epoch_outputs(net::Net)
    epoch = get_epoch(net)
    while get_epoch(net) == epoch
        result_forward = forward(net, 0.0)
    end
    #read out outputs
    result =  hcat(net.states[end].outputs[1]...)
    #reset network
    Mocha.reset_outputs(net.states[end])
    return  result

end

"""
Evaluate accuracy and loglikelihood over whole dataset.
"""
function evaluateTestNN(mw::MWrap, para::Array{Float64,1}, batch_size::Int64)

    mw = setparams!(mw,para)
    Mocha.reset_statistics(mw.net)
    # sum of normalised batch losses.
    ll = forward_epoch_loss(mw.net)
    # scale up for full loss.
    ll = ll * batch_size
    return (mw.net.states[end].accuracy, ll)
end

"""
Evaluate predictive probabilities on test set.
"""
function getpredprobs(mw::MWrap, para::Array{Float64,1}, batch_size::Int64)
    setparams!(mw,para)
    Mocha.reset_statistics(mw.net)
    predprobs = forward_epoch_outputs(mw.net)

    return predprobs
end

"""
Sets weights to para
"""
function setparams!(mw::MWrap,para::Array{Float64,1})
  r = 1
 for i in mw.layersWParam
      for j = 1:length(mw.net.states[i].parameters)
           copy!( mw.net.states[i].parameters[j].blob,
                   para[mw.paramRanges[1,r]:mw.paramRanges[2,r]] )
           r+=1
      end
 end
 mw
end

"""
Get number of weights
"""
function getnparams(mw::MWrap)
  mw.nparams
end

"""
Get weights
"""
function getparams(mw::MWrap)
  para::Array{Float64,1} = zeros(mw.nparams)
  r = 1
  for i in mw.layersWParam
      for j = 1:length(mw.net.states[i].parameters)
         x = zeros(mw.paramRanges[2,r]-mw.paramRanges[1,r]+1)
         copy!(x, mw.net.states[i].parameters[j].blob)
         para[mw.paramRanges[1,r]:mw.paramRanges[2,r]] = x
         r+=1
      end
  end
  para
end

"""
Get gradient. This only reads out the gradient and does not trigger a forward backward step.
"""
function getgrad(mw::MWrap)
  grad::Array{Float64,1} = zeros(mw.nparams)
  r=1
  for i in mw.layersWParam
      for j=1:length(mw.net.states[i].parameters)
          p = unsafe_wrap(Array{Float64},pointer(grad, mw.paramRanges[1,r] ), mw.paramRanges[2,r]-mw.paramRanges[1,r]+1)
          copy!(p, mw.net.states[i].parameters[j].gradient)
          r+=1
      end
  end
  grad
end

"""
Xavier initilization
"""
function init_xavier(mw::MWrap)
  for i  = 2:(length(mw.net.states)-1)
    if Mocha.has_param(mw.net.layers[i])
      #for j = 1:length(net.states[i].parameters)
      weights_blob = mw.net.states[i].parameters[1].blob
      bias_blob = mw.net.states[i].parameters[2].blob

      # Initialize weights
      # NOTE: Here is where I make a change from the Mocha code to match the equations in the Glorot and Bengio paper
      fan_in = get_fea_size(weights_blob) + size(weights_blob)[end]
      scale = sqrt(6.0 / fan_in)
      init_val = rand(eltype(weights_blob), size(weights_blob)) * 2scale - scale
      copy!(mw.net.states[i].parameters[1].blob, init_val)

      # Initialize biases
      fill!(mw.net.states[i].parameters[2].blob, 0)
    end
  end
end

"""
Initialization based on simple fanin
"""
function init_simple_fanin(mw::MWrap)
  for i  = 2:(length(mw.net.states)-1)
    if Mocha.has_param(mw.net.layers[i])
      #for j = 1:length(net.states[i].parameters)
      weights_blob = mw.net.states[i].parameters[1].blob
      bias_blob = mw.net.states[i].parameters[2].blob

      # Initialize weights
      # NOTE: Here is where I make a change from the Mocha code to match the equations in the Glorot and Bengio paper
      fan_in = get_fea_size(weights_blob)
      scale = sqrt(1.0 / fan_in)
      init_val = rand(eltype(weights_blob), size(weights_blob)) * 2scale - scale
      copy!(mw.net.states[i].parameters[1].blob, init_val)

      # Initialize biases
      fill!(mw.net.states[i].parameters[2].blob, 0)
    end
  end
end

"""
Gaussian initialization for weights, zeros for biases.
"""
function init_gaussian(mw::MWrap, initvar::Float64)
  for i  = 2:(length(mw.net.states)-1)
    if Mocha.has_param(mw.net.layers[i])
      #for j = 1:length(net.states[i].parameters)
      weights_blob = mw.net.states[i].parameters[1].blob
      bias_blob = mw.net.states[i].parameters[2].blob

      # Initialize weights
      init_val = randn(size(weights_blob)) * sqrt(initvar)
      copy!(mw.net.states[i].parameters[1].blob, init_val)

      # Initialize biases
      fill!(mw.net.states[i].parameters[2].blob, 0)
    end
  end
end

"""
Uniform initialization for weights, zeros for biases.
"""
function init_uniform(mw::MWrap, initvar::Float64)
  for i = 2:(length(mw.net.states)-1)
    if Mocha.has_param(mw.net.layers[i])
      #for j = 1:length(net.states[i].parameters)
      weights_blob = mw.net.states[i].parameters[1].blob
      bias_blob = mw.net.states[i].parameters[2].blob

      # Initialize weights
      b = sqrt(12 * initvar)
      init_val = rand(eltype(weights_blob), size(weights_blob)) * b - 0.5 * b
      copy!(mw.net.states[i].parameters[1].blob, init_val)

      # Initialize biases
      fill!(mw.net.states[i].parameters[2].blob, 0)
    end
  end
end

"""
Matches indices to weight/bias blocks. Called during MWrap construction.
"""
function parameterisation(net::Net)
    istart = 1
    iend   = 1
    layersWParam=Int64[]
    paramRangesI=Int64[]
    for i  = 2:(length(net.states)-1)
        if Mocha.has_param(net.layers[i])
            push!(layersWParam,i)
            for j = 1:length(net.states[i].parameters)
                iend   = istart+length(net.states[i].parameters[j].blob)-1
                push!(paramRangesI,istart)
                push!(paramRangesI,iend)
                istart = iend+1
            end
        end
    end
    nparams = paramRangesI[end]
    paramRanges = reshape(paramRangesI,(2,div(length(paramRangesI),2)))
    return(nparams,paramRanges,layersWParam)
end
end


text:
River = Symbol
Area = Int
Plant = Symbol

mutable struct HydroPlantData{T <: AbstractFloat}
    M₀::T                # Initial reservoir contents
    M̄::T                 # Maximum reservoir capacities
    H̄::T                 # Maximal production
    Q̄::T                 # Maximal discharge
    S̱::T                 # Minimum spillage
    Rq::Int              # Discharge flow time in minutes
    Rs::Int              # Spillage flow time in minutes
    Mt::Int              # Maintenance time in whole hours
    Q̃::T                 # Yearly mean flow of each plant
    V::T                 # Local inflow

    function HydroPlantData(M₀::AbstractFloat,
                            M̄::AbstractFloat,
                            H̄::AbstractFloat,
                            Q̄::AbstractFloat,
                            S̱::AbstractFloat,
                            Rq::Integer,
                            Rs::Integer,
                            Mt::Integer,
                            Q̃::AbstractFloat,
                            V::AbstractFloat)
        T = promote_type(typeof(M₀),typeof(M̄),typeof(H̄),eltype(Q̄),typeof(S̱),typeof(Q̃),typeof(V),Float32)
        return new{T}(M₀,M̄,H̄,Q̄,S̱,Rq,Rs,Mt,Q̃,V)
    end
end
Base.eltype(::HydroPlantData{T}) where T <: AbstractFloat = T

struct HydroPlantCollection{T <: AbstractFloat, S}
    # Plants
    # ========================================================
    plants::Vector{Plant}                   # All possible plant instances
    rivers::Dict{River,Vector{Plant}}       # Plants sorted according to river
    areas::Dict{Area,Vector{Plant}}         # Plants sorted according to price area
    # Parameters
    # ========================================================
    plantdata::Dict{Plant,HydroPlantData{T}}     # Data for each plant
    segmenter::Segmenter{S}                      # Segmenter
    Qd::Dict{Plant,Vector{Plant}}                # All discharge outlets located downstream (including itself)
    Qu::Dict{Plant,Vector{Plant}}                # Discharge outlet(s) located directly upstream
    Sd::Dict{Plant,Vector{Plant}}                # Spillage outlets located downstream (including itself)
    Su::Dict{Plant,Vector{Plant}}                # Spillage outlet(s) located directly upstream

    function HydroPlantCollection(segmenter::Segmenter{S},
                                  plants::Dict{Plant,HydroPlantData{T}},
                                  Qd::Dict{Plant,Vector{Plant}},
                                  Qu::Dict{Plant,Vector{Plant}},
                                  Sd::Dict{Plant,Vector{Plant}},
                                  Su::Dict{Plant,Vector{Plant}}) where {T <: AbstractFloat,S}
        return new{T,S}(collect(keys(plants)),
                        Dict{River,Vector{Plant}}(),
                        Dict{Area,Vector{Plant}}(),
                        plants,
                        segmenter,
                        Qd,
                        Qu,
                        Sd,
                        Su)
    end

    function HydroPlantCollection(segmenter::Segmenter{S},
                                  plantnames::Vector{String},
                                  Qlinks::Vector{Int},
                                  Slinks::Vector{Int},
                                  H̄::AbstractVector,
                                  Q̄::AbstractVector,
                                  S̱::AbstractVector,
                                  M̄::AbstractVector,
                                  Q̃::AbstractVector,
                                  Rq::Vector{Int},
                                  Rs::Vector{Int},
                                  Mt::Vector{Int},
                                  rivers::Vector{String},
                                  areas::Vector{Area}) where S
        T = promote_type(eltype(H̄),eltype(Q̄),eltype(S̱),eltype(M̄),eltype(Q̃),Float32)
        modeldata = new{T,S}(Vector{Plant}(),
                             Dict{River,Vector{Plant}}(),
                             Dict{Area,Vector{Plant}}(),
                             Dict{Plant,HydroPlantData{T}}(),
                             segmenter,
                             Dict{Plant,Vector{Plant}}(),
                             Dict{Plant,Vector{Plant}}(),
                             Dict{Plant,Vector{Plant}}(),
                             Dict{Plant,Vector{Plant}}(),)
        define_model_parameters(modeldata,plantnames,Qlinks,Slinks,convert(Vector{T},H̄),convert(Vector{T},Q̄),convert(Vector{T},S̱),convert(Vector{T},M̄),convert(Vector{T},Q̃),Rq,Rs,Mt,rivers,areas)
        return modeldata
    end

    function HydroPlantCollection(::Type{T}, segmenter::Segmenter{S}, plantfilename::String) where {T <: AbstractFloat, S}
        modeldata = new{T,S}(Vector{Plant}(),
                             Dict{River,Vector{Plant}}(),
                             Dict{Area,Vector{Plant}}(),
                             Dict{Plant,HydroPlantData{T}}(),
                             segmenter,
                             Dict{Plant,Vector{Plant}}(),
                             Dict{Plant,Vector{Plant}}(),
                             Dict{Plant,Vector{Plant}}(),
                             Dict{Plant,Vector{Plant}}(),)
        define_model_parameters(modeldata, plantfilename)
        return modeldata
    end
end
Base.eltype(::HydroPlantCollection{T}) where T <: AbstractFloat = T

function show(io::IO, collection::HydroPlantCollection)
    if get(io, :multiline, false)
        print(io,"HydroPlantCollection")
    else
        println(io,"Collection of Hydropower Plants")
        if !isempty(collection.rivers)
            println(io,"Rivers:")
            for river in keys(collection.rivers)
                println(io,string(river))
            end
        end
    end
end

function HydroPlantCollection(plantnames::Vector{String},
                              Qlinks::Vector{Int},
                              Slinks::Vector{Int},
                              H̄::AbstractVector,
                              Q̄::AbstractVector,
                              S̱::AbstractVector,
                              M̄::AbstractVector,
                              Q̃::AbstractVector,
                              Rq::Vector{Int},
                              Rs::Vector{Int},
                              Mt::Vector{Int},
                              rivers::Vector{String},
                              areas::Vector{Area})
    return HydroPlantCollection(Segmenter{2}(), plantnames, Qlinks, Slinks, H̄, Q̄, S̱, M̄, Q̃, Rq, Rs, Mt, rivers, areas)
end

function HydroPlantCollection(plant_filename::String)
    return HydroPlantCollection(Float64, Segmenter{2}(), plant_filename)
end

Base.getindex(collection::HydroPlantCollection, plant::Plant) = collection.plantdata[plant]

function Q̄(collection::HydroPlantCollection, plant::Plant, s::Integer)
    Q̄s, μs = segment(collection.segmenter, collection[plant].Q̄, collection[plant].H̄)
    return Q̄s[s]
end
function μ(collection::HydroPlantCollection, plant::Plant, s::Integer)
    Q̄s, μs = segment(collection.segmenter, collection[plant].Q̄, collection[plant].H̄)
    return μs[s]
end
function %(collection::HydroPlantCollection, s::Integer)
    return segment_percentage(collection.segmenter, s)
end

function define_plants!(collection::HydroPlantCollection, plantnames::Vector{String})
    for plantname in plantnames
        push!(collection.plants,Plant(filter(x->!isspace(x),plantname)))
    end
end

function define_rivers!(collection::HydroPlantCollection, rivernames::Vector{String})
    for (p,rivername) in enumerate(rivernames)
        plant = collection.plants[p]
        river = River(filter(x->!isspace(x),rivername))
        if !haskey(collection.rivers,river)
            collection.rivers[river] = Plant[]
        end
        push!(collection.rivers[river],plant)
    end
end

function define_areas!(collection::HydroPlantCollection, areas::Vector{Area})
    for (p,area) in enumerate(areas)
        plant = collection.plants[p]
        if !haskey(collection.areas,area)
            collection.areas[area] = Plant[]
        end
        push!(collection.areas[area],plant)
    end
end

function define_plant_topology!(links::Dict{Plant,Plant},
                                downstream_plants::Dict{Plant,Vector{Plant}},
                                upstream_plants::Dict{Plant,Vector{Plant}})
    linker = (downstream_plants, current, links ) -> begin
        if current == :NoLink
            return
        end
        push!(downstream_plants, current)
        linker(downstream_plants, links[current], links)
    end

    for (plant,link) in links
        if link == :NoLink
            continue
        end
        push!(downstream_plants[plant],plant)
        linker(downstream_plants[plant],link,links)
        push!(upstream_plants[link],plant)
    end
end

function calculate_inflow(plant::Plant,w::Dict{Plant,<:AbstractFloat},upstream_plants::Vector{Plant})
    V = w[plant]
    if !isempty(upstream_plants)
        V -= sum(w[p] for p in upstream_plants)
    end
    return V
end

function define_model_parameters(collection::HydroPlantCollection{T,S},
                                 plantnames::Vector{String},
                                 Qlinks::Vector{Int},
                                 Slinks::Vector{Int},
                                 H̄::Vector{T},
                                 Q̄::Vector{T},
                                 S̱::Vector{T},
                                 M̄::Vector{T},
                                 Q̃::Vector{T},
                                 Rq::Vector{Int},
                                 Rs::Vector{Int},
                                 Mt::Vector{Int},
                                 rivers::Vector{String},
                                 areas::Vector{Area}) where {T <: AbstractFloat, S}
    Vsf = 0.2278           # Scale factor for local inflow (0.2278)
    δ = 0.363              # Initial reservoir content factor (0.363)
    M_end = 0.89           # Target water level as factor of M_0 (0.89)

    # Define available plants
    define_plants!(collection, plantnames)

    # Define sortings of plants
    define_rivers!(collection,rivers)
    define_areas!(collection,areas)

    for p in collection.plants
        collection.Qd[p] = Plant[]
        collection.Qu[p] = Plant[]
        collection.Sd[p] = Plant[]
        collection.Su[p] = Plant[]
    end

    # Define plant topologies
    links_as_plants(links) = Dict([(collection.plants[p],(l == 0 ? :NoLink : collection.plants[l])) for (p,l) in enumerate(links)])
    define_plant_topology!(links_as_plants(Qlinks), collection.Qd, collection.Qu)
    define_plant_topology!(links_as_plants(Slinks), collection.Sd, collection.Su)

    w = Dict(zip(collection.plants,Vsf*Q̃))

    for i in 1:length(plantnames)
        p = collection.plants[i]
        V = calculate_inflow(p,w,collection.Qu[p])
        collection.plantdata[p] = HydroPlantData(δ*M̄[i],
                                                 M̄[i],
                                                 H̄[i],
                                                 Q̄[i],
                                                 S̱[i],
                                                 Rq[i],
                                                 Rs[i],
                                                 Mt[i],
                                                 Vsf*Q̃[i],
                                                 V)
    end
end

function define_model_parameters(collection::HydroPlantCollection{T}, plantdata::Matrix) where T <: AbstractFloat
    @assert size(plantdata,2) == 13 "Invalid plant data format"

    define_model_parameters(collection,
                            convert(Vector{String},plantdata[:,1]),
                            convert(Vector{Int},plantdata[:,2]),
                            convert(Vector{Int},plantdata[:,3]),
                            convert(Vector{T},plantdata[:,4]),
                            convert(Vector{T},plantdata[:,5]),
                            convert(Vector{T},plantdata[:,6]),
                            convert(Vector{T},plantdata[:,7]),
                            convert(Vector{T},plantdata[:,8]),
                            convert(Vector{Int},plantdata[:,9]),
                            convert(Vector{Int},plantdata[:,10]),
                            convert(Vector{Int},plantdata[:,11]),
                            convert(Vector{String},plantdata[:,12]),
                            convert(Vector{Area},plantdata[:,13]))
end

function define_model_parameters(collection::HydroPlantCollection, plantfilename::String)
    plantdata = readdlm(plantfilename, ',')
    define_model_parameters(collection,plantdata[2:end,:])
end

function plants_in_river(collection::HydroPlantCollection,river::River)
    if river == :All
        return collection.plants
    end
    if !haskey(collection.rivers,river)
        warn(string("Invalid river name: ", river))
        return Plant[]
    end
    return collection.rivers[river]
end

function plants_in_river(collection::HydroPlantCollection,rivers::Vector{River})
    plants = Plant[]
    for river in rivers
        append!(plants,plants_in_river(collection,river))
    end
    return plants
end

function plants_in_area(collection::HydroPlantCollection,area::Area)
    if area == 0
        return collection.plants
    end
    if !haskey(collection.areas,area)
        warn(string("Invalid area: ",area))
        return Plant[]
    end
    return collection.areas[area]
end

function plants_in_area(collection::HydroPlantCollection,areas::Vector{Area})
    plants = Plant[]
    for area in areas
        append!(plants,plants_in_area(collection,area))
    end
    return plants
end

function plants_in_areas_and_rivers(collection::HydroPlantCollection,areas::Vector{Area},rivers::Vector{River})
    return plants_in_river(collection,rivers) ∩ plants_in_area(collection,areas)
end


text:
module NonLinearDynamicsCourse

    using Plots
    using LinearAlgebra
    using DifferentialEquations
    using ForwardDiff
    using IntervalRootFinding
    using StaticArrays

    export  flux1D,
            potential1D,
            flux2d_vectorfield,
            flux2d_nullclines,
            flux2d_animated,
            flux2d_manifolds,
            flux2d_forced,
            classification_linear,
            phase_portrait,
            attractor_basin,
            poincare_forced,
            poincare_forced_zoom,
            recurrence_plot,
            saddle_orbit2D,
            saddle_manifolds_forced,
            butterfly

    include("DNL_utils.jl")

end


text:
using Documenter, TableWidgets

makedocs(
    format = :html,
    sitename = "TableWidgets",
    authors = "Pietro Vertechi",
    pages = [
        "Introduction" => "index.md",
        "API reference" => "api_reference.md",
    ]
)

deploydocs(
    repo = "github.com/piever/TableWidgets.jl.git",
    target = "build",
    julia  = "1.0",
    osname = "linux",
    deps   = nothing,
    make   = nothing
)