text:
"""
    allequal(x::Array)
    allequal(x::Tuple)
Return the equality of all elements of the array
"""
function allequal(x::Array; comparison=isequal)::Bool
    all(comparison.(first(x),  x))
end

allequal(x::Vector{Bool})::Bool = !any(x) || all(x)

allequal(x::NTuple{N,Bool}) where {N} = !any(x) || all(x)

function allequal(x::Tuple; comparison=isequal)::Bool
    all(comparison.(first(x),  x))
end

function allequal(x...; comparison=isequal)::Bool
    all(comparison.(first(x),  x))
end

"""
    average(a::T, b::T) where {T<:AbstractFloat}

Return the average of `a` and `b`
"""
average(a::T, b::T) where {T<:AbstractFloat} = (a + b) / 2

"""
    cpad(s::AbstractString, n::Integer)

Return a string of length `n` containing `s` in the center (more-or-less).
"""
cpad(s::String, n::Integer) = rpad(lpad(s, (n + textwidth(s)) >> 1), n)

"""
densify(S::SparseMatrix, threshold=0.25)

Convert sparse `S` to `Diagonal` if `S` is diagonal or to `Array(S)` if
the proportion of nonzeros exceeds `threshold`.
"""
function densify(A::SparseMatrixCSC, threshold::Real = 0.25)
    m, n = size(A)
    if m == n && isdiag(A)  # convert diagonal sparse to Diagonal
        # the diagonal is always dense (otherwise rank deficit)
        # so make sure it's stored as such
        Diagonal(Vector(diag(A)))
    elseif nnz(A) / (m * n) ≤ threshold
        A
    else
        Array(A)
    end
end
densify(A::AbstractMatrix, threshold::Real = 0.3) = A

densify(A::SparseVector, threshold::Real = 0.3) = Vector(A)
densify(A::Diagonal{T,SparseVector{T,Ti}}, threshold::Real = 0.3) where {T,Ti} =
    Diagonal(Vector(A.diag))

"""
    RaggedArray{T,I}

A "ragged" array structure consisting of values and indices

# Fields
- `vals`: a `Vector{T}` containing the values
- `inds`: a `Vector{I}` containing the indices

For this application a `RaggedArray` is used only in its `sum!` method.
"""
struct RaggedArray{T,I}
    vals::Vector{T}
    inds::Vector{I}
end

function Base.sum!(s::AbstractVector{T}, a::RaggedArray{T}) where {T}
    for (v, i) in zip(a.vals, a.inds)
        s[i] += v
    end
    s
end

function rownormalize!(A::AbstractMatrix)
    for r in eachrow(A)
        # all zeros arise in zerocorr situations
        if !iszero(r)
            normalize!(r)
        end
    end
    A
end

"""
    ltriindprs

A row-major order `Vector{NTuple{2,Int}}` of indices in the strict lower triangle.
"""
const ltriindprs = NTuple{2,Int}[]

function checkindprsk(k::Integer)
    kchoose2 = (k * (k - 1)) >> 1
    if length(ltriindprs) < kchoose2
        sizehint!(empty!(ltriindprs), kchoose2)
        for i in 1:k, j in 1:(i-1)
            push!(ltriindprs, (i,j))
        end
    end
    ltriindprs
end

"""
    replicate(f::Function, n::Integer; use_threads=false)

Return a vector of the values of `n` calls to `f()` - used in simulations where the value of `f` is stochastic.

Note that if `f()` is not thread-safe or depends on a non thread-safe RNG,
    then you must set `use_threads=false`. Also note that ordering of replications
    is not guaranteed when `use_threads=true`, although the replications are not
    otherwise affected for thread-safe `f()`.
"""
function replicate(f::Function, n::Integer; use_threads=false)
    if use_threads
        # no macro version yet: https://github.com/timholy/ProgressMeter.jl/issues/143
        p = Progress(n)
        # get the type
        rr = f()
        next!(p)
        # pre-allocate
        results = [rr for _ in Base.OneTo(n)]
        Threads.@threads for idx = 2:n
            results[idx] = f()
            next!(p)
        end
    else
        results = @showprogress [f() for _ in Base.OneTo(n)]
    end
    results
end

cacheddatasets = Dict{String,Any}()
"""
    dataset(nm)

Return the data frame of test data set named `nm`, which can be a `String` or `Symbol`
"""
function dataset(nm::AbstractString)
    get!(cacheddatasets, nm) do
        path = joinpath(TestData, nm * ".feather")
        if !isfile(path)
            throw(ArgumentError(
                "Dataset \"$nm\" is not available.\nUse MixedModels.datasets() for available names."))
        end
        Feather.read(path)
    end
end
dataset(nm::Symbol) = dataset(string(nm))

"""
    datasets()

Return a vector of names of the available test data sets
"""
datasets() = first.(Base.Filesystem.splitext.(filter(Base.Fix2(endswith, ".feather"), readdir(TestData))))

"""
    PCA{T<:AbstractFloat}

Principal Components Analysis

## Fields

* `covcorr` covariance or correlation matrix
* `sv` singular value decomposition
* `corr` is this a correlation matrix
"""
struct PCA{T<:AbstractFloat}
    covcor::Symmetric{T,Matrix{T}}
    sv::SVD{T,T,Matrix{T}}
    corr::Bool
end

"""
    PCA(::AbstractMatrix; corr::Bool=true)
    PCA(::ReMat; corr::Bool=true)
    PCA(::LinearMixedModel; corr::Bool=true)

Constructs a [`MixedModels.PCA`](@ref]) object from a covariance matrix.

For `LinearMixedModel`, a named tuple of PCA on each of the random-effects terms is returned.

If `corr=true`, then the covariance is first standardized to the correlation scale.
"""

function PCA(covfac::AbstractMatrix; corr::Bool=true)
    covf = corr ? rownormalize!(copy(covfac)) : copy(covfac)
    PCA(Symmetric(covf*covf', :L), svd(covf), corr)
end

function Base.getproperty(pca::PCA, s::Symbol)
    if s == :cumvar
        cumvv = cumsum(abs2.(pca.sv.S))
        cumvv ./ last(cumvv)
    elseif s == :loadings
        pca.sv.U
    else
        getfield(pca, s)
    end
end

Base.propertynames(pca::PCA, private = false) = (
    :covcor,
    :sv,
    :corr,
    :cumvar,
    :loadings,
#    :rotation,
)

Base.show(pca::PCA;
        ndigitsmat=2, ndigitsvec=2, ndigitscum=4,
        covcor=true, loadings=true, variances=false, stddevs=false) =
        Base.show(Base.stdout, pca,
                    ndigitsmat=ndigitsmat,
                    ndigitsvec=ndigitsvec,
                    ndigitscum=ndigitscum,
                    covcor=covcor,
                    loadings=loadings,
                    variances=variances,
                    stddevs=stddevs)

function Base.show(io::IO, pca::PCA;
        ndigitsmat=2, ndigitsvec=2, ndigitscum=4,
        covcor=true, loadings=true, variances=false, stddevs=false)
    println(io)
    if covcor
        println(io,
                "Principal components based on ",
                pca.corr ? "correlation" : "(relative) covariance",
                " matrix")
        # only display the lower triangle of symmetric matrix
        Base.print_matrix(io, round.(LowerTriangular(pca.covcor), digits=ndigitsmat))
        println(io)
    end
    if stddevs
        println(io, "Standard deviations:")
        sv = pca.sv
        show(io, round.(sv.S, digits=ndigitsvec))
        println(io)
    end
    if variances
        println(io, "Variances:")
        vv = abs2.(sv.S)
        show(io, round.(vv, digits=ndigitsvec))
        println(io)
    end
    println(io, "Normalized cumulative variances:")
    show(io, round.(pca.cumvar, digits=ndigitscum))
    println(io)
    if loadings
        println(io, "Component loadings")
        Base.print_matrix(io, round.(pca.loadings, digits=ndigitsmat))
    end

    nothing
end


text:
#####
##### Custom Threading
#####

"""
    threaded(f, A::AbstractArray, args...; [iterations = 1], [nthreads = Threads.nthreads()])

Split `A` into `nthreads` separate, sequential views.
For each view `v`, call `f(v, args...)` on a distinct thread.
This inner call is performed `iterations` times.

**Note**: It is required that `mod(length(A), nthreads) == 0` so that `A` can be evenly
broken into views.
"""
function threaded(f, A::AbstractArray, args...; iterations = 1, nthreads = Threads.nthreads())
    @assert iszero(mod(length(A), nthreads))

    step = div(length(A), nthreads)
    Threads.@threads for i in 1:nthreads
        start = step * (i-1) + 1
        stop = step * i
        x = view(A, start:stop)

        # Run the inner loop
        for j in 1:iterations
            f(x, args...)
        end
    end
    return nothing
end

ext: jl


text:
###
# UCS Oracle
###
# At the moment only supports cutting planes
import JuMP.UnsetSolver  #remove this dependence

export UCSOracle

type UCSOracle <: AbstractOracle
    eps_kappa::Float64
    Gamma1::Float64
    Gamma2::Float64

    # Cutting plane algorithm
    muhat::Vector{Float64}
    covbar::Array{Float64, 2}
    C::Array{Float64, 2}  #C'C = covbar
    cut_tol::Float64  ##defaults to 1e-6
    cut_model::Model
    cut_vars::Vector{Variable}
    unbounded_support::Bool   #True if no support bounds on any uncertainties

    # Other options
    debug_printcut::Bool
end

function UCSOracle(muhat, covhat, Gamma1, Gamma2, eps_)
    covbar = covhat + Gamma2 * eye(size(covhat)...)
    C = chol(covbar, :U)   #C'C = covbar
    UCSOracle(  eps_, Gamma1, Gamma2, muhat, covbar, C,
                1e-6, Model(), Variable[], true, false)
end

function UCSOracle(data, eps_, alpha1, alpha2; numBoots=10000)
    muhat  = vec(mean(data, 1))
    covhat = cov(data)
    Gamma1 = boot_mu(data, alpha1, numBoots)
    Gamma2 = boot_sigma(data, alpha2, numBoots)
    UCSOracle(muhat, covhat, Gamma1, Gamma2, eps_)
end

#preferred interface
UCSOracle(data, eps_, alpha; numBoots=10000) = UCSOracle(data, eps_, alpha/2, alpha/2, numBoots=numBoots)

#the supp fcn when support = Rd
function suppFcnUCSRd(xs, muhat, Gamma1, covbar, eps_k, cut_sense=:Max)
    toggle = 1.
    if cut_sense == :Min
        xs = copy(-xs)
        toggle = -1.
    end
    norm_x = norm(xs)
    sig_term = sqrt(xs' * covbar * xs)[1]
    ustar =  muhat + Gamma1/norm_x * xs
    ustar += kappa(eps_k) * sig_term / norm_x / norm_x * xs
    dot(ustar, xs)*toggle, ustar
end

function suppFcn(xs, w::UCSOracle, cut_sense)
    if w.unbounded_support
        return suppFcnUCSRd(xs, w.muhat, w.Gamma1, w.covbar, w.eps_kappa, cut_sense)
    else
        setObjective(w.cut_model, cut_sense, sum([xs[ix] * w.cut_vars[ix] for ix=1:length(w.cut_vars)]))

        cut_solve_status = solve(w.cut_model, suppress_warnings=true)
        cut_solve_status != :Optimal && error("Cutting plane problem failed: $cut_solve_status")
        ustar = getValue(w.cut_vars)
        zstar = getObjectiveValue(w.cut_model)
        return zstar, ustar
    end
end

function setup(w::UCSOracle, rm::Model, prefs)
    # Extract preferences we care about
    w.debug_printcut = get(prefs, :debug_printcut, false)
    w.cut_tol        = get(prefs, :cut_tol, w.cut_tol)

    rd = JuMPeR.getRobust(rm)
    w.cut_model.solver   = isa(rd.cutsolver, UnsetSolver) ? rm.solver : rd.cutsolver

    d = size(w.covbar, 1)
    @assert (rd.numUncs == d) "Num Uncertainties doesn't match columns in data"
    @assert (length(rd.uncertaintyset) == 0) #does not support additional cnsts on unctertainties for now
    #w.cut_model.colCat   = rd.uncCat  #only supports continuous variables for now

    #if there are no bounds, can use the simpler support function
    w.unbounded_support = true
    for i = 1:d
        if (rd.uncLower[i] > -Inf) || (rd.uncUpper[i] < Inf )
            w.unbounded_support=false
            break
        end
    end

    #else build an SOCP for separation
    if ! w.unbounded_support
        @defVar(w.cut_model, rd.uncLower[i] <= us[i=1:d] <= rd.uncUpper[i])
        @defVar(w.cut_model, z1[1:d])
        @defVar(w.cut_model, z2[1:d])
        for i = 1:d
            setName(us[i], rd.uncNames[i])
            addConstraint(w.cut_model, us[i] == w.muhat[i] + z1[i] + dot(w.C[:, i], z2))
        end

        addConstraint(w.cut_model, dot(z1, z1) <= w.Gamma1^2 )
        addConstraint(w.cut_model, dot(z2, z2) <= kappa(w.eps_kappa)^2)
        w.cut_vars = us[:]
    end
end


text:
function BigReal(arg0::BigDecimal)
    return BigReal((BigDecimal,), arg0)
end

function BigReal(arg0::BigInteger)
    return BigReal((BigInteger,), arg0)
end

function BigReal(arg0::BigInteger, arg1::MathContext)
    return BigReal((BigInteger, MathContext), arg0, arg1)
end

function BigReal(arg0::BigInteger, arg1::jint)
    return BigReal((BigInteger, jint), arg0, arg1)
end

function BigReal(arg0::BigInteger, arg1::jint, arg2::MathContext)
    return BigReal((BigInteger, jint, MathContext), arg0, arg1, arg2)
end

function BigReal(arg0::JString)
    return BigReal((JString,), arg0)
end

function BigReal(arg0::JString, arg1::MathContext)
    return BigReal((JString, MathContext), arg0, arg1)
end

function BigReal(arg0::Vector{jchar})
    return BigReal((Vector{jchar},), arg0)
end

function BigReal(arg0::Vector{jchar}, arg1::MathContext)
    return BigReal((Vector{jchar}, MathContext), arg0, arg1)
end

function BigReal(arg0::Vector{jchar}, arg1::jint, arg2::jint)
    return BigReal((Vector{jchar}, jint, jint), arg0, arg1, arg2)
end

function BigReal(arg0::Vector{jchar}, arg1::jint, arg2::jint, arg3::MathContext)
    return BigReal((Vector{jchar}, jint, jint, MathContext), arg0, arg1, arg2, arg3)
end

function BigReal(arg0::jdouble)
    return BigReal((jdouble,), arg0)
end

function BigReal(arg0::jdouble, arg1::MathContext)
    return BigReal((jdouble, MathContext), arg0, arg1)
end

function BigReal(arg0::jint)
    return BigReal((jint,), arg0)
end

function BigReal(arg0::jint, arg1::MathContext)
    return BigReal((jint, MathContext), arg0, arg1)
end

function BigReal(arg0::jlong)
    return BigReal((jlong,), arg0)
end

function BigReal(arg0::jlong, arg1::MathContext)
    return BigReal((jlong, MathContext), arg0, arg1)
end

function add(obj::BigReal, arg0::BigReal)
    return jcall(obj, "add", BigReal, (BigReal,), arg0)
end

function big_decimal_value(obj::BigReal)
    return jcall(obj, "bigDecimalValue", BigDecimal, ())
end

function compare_to(obj::BigReal, arg0::BigReal)
    return jcall(obj, "compareTo", jint, (BigReal,), arg0)
end

function divide(obj::BigReal, arg0::BigReal)
    return jcall(obj, "divide", BigReal, (BigReal,), arg0)
end

function double_value(obj::BigReal)
    return jcall(obj, "doubleValue", jdouble, ())
end

function equals(obj::BigReal, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_field(obj::BigReal)
    return jcall(obj, "getField", Field, ())
end

function get_rounding_mode(obj::BigReal)
    return jcall(obj, "getRoundingMode", RoundingMode, ())
end

function get_scale(obj::BigReal)
    return jcall(obj, "getScale", jint, ())
end

function hash_code(obj::BigReal)
    return jcall(obj, "hashCode", jint, ())
end

function multiply(obj::BigReal, arg0::BigReal)
    return jcall(obj, "multiply", BigReal, (BigReal,), arg0)
end

function multiply(obj::BigReal, arg0::jint)
    return jcall(obj, "multiply", BigReal, (jint,), arg0)
end

function negate(obj::BigReal)
    return jcall(obj, "negate", BigReal, ())
end

function reciprocal(obj::BigReal)
    return jcall(obj, "reciprocal", BigReal, ())
end

function set_rounding_mode(obj::BigReal, arg0::RoundingMode)
    return jcall(obj, "setRoundingMode", void, (RoundingMode,), arg0)
end

function set_scale(obj::BigReal, arg0::jint)
    return jcall(obj, "setScale", void, (jint,), arg0)
end

function subtract(obj::BigReal, arg0::BigReal)
    return jcall(obj, "subtract", BigReal, (BigReal,), arg0)
end


text:
using Parameters:@with_kw

"""
Stores information of a self-play game.
"""
mutable struct GameHistory
    observation_history::Array{Float32,4}
    action_history::Vector{Int}
    reward_history::Vector{Float32}
    to_play_history::Vector{Int}
    child_visits::Matrix{Float32}
    root_values::Vector{Float32}
    reanalysed_predicted_root_values
    priorities
    game_priority
end

@with_kw struct Config
    seed::Int = 1337
    observation_shape::Tuple{Int,Int,Int} # Dimensions of the game observation, must be 3D (WHC). For a 1D array, please reshape it to (1, 1, length of array)
    action_space
	players
	stacked_observations::Int # Number of previous observations and previous actions to add to the current observation
    muzero_player::Int = 1 # Turn Muzero begins to play (1: MuZero plays first, 2: MuZero plays second)
    opponent::String = "expert" # Hard coded agent that MuZero faces to assess its progress in multiplayer games. It doesn't influence training. None, "random" or "expert" if implemented in the Game class
	intermediate_rewards::Bool = false

    num_workers::Int # Number of simultaneous threads/workers self-playing to feed the replay buffer
    selfplay_on_gpu::Bool = false
    max_moves::Int # Game dependent
    temperature_threshold::Union{Int,Nothing} = nothing # Number of moves before dropping the temperature given by visit_softmax_temperature_fn to 0 (ie selecting the best action). If None, visit_softmax_temperature_fn is used every time

    dirichlet_α::Float32 = 0.25
    exploration_ϵ::Float32 = 0.25
    pb_c_base::Int = 19652
    pb_c_init::Float32 = 1.25
    discount::Float32 = 0.997 # Chronological discount of the reward
    num_iters::Int # Number of future moves self-simulated

    replay_buffer_size::Int = 10000 # Number of self-play games to keep in the replay buffer
    num_unroll_steps::Int # Number of game moves to keep for every batch element
    td_steps::Int # Number of steps in the future to take into account for calculating the target value
    PER::Bool # Prioritized Replay (See paper appendix Training), select in priority the elements in the replay buffer which are unexpected for the network
    PER_alpha::Int = 1 # How much prioritization is used, 0 corresponding to the uniform case, paper suggests 1

    results_path::String = mkpath("./results")
    networks_path::String = mkpath("./networks")
    training_steps::Int # Total number of training steps (ie weights update according to a batch)
    batch_size::Int # Number of parts of games to train on at each training step
    checkpoint_interval::Int = 10 # Number of training steps before using the model for self-playing
    value_loss_weight::Float32 = 0.25 # Scale the value loss to avoid overfitting of the value function, paper recommends 0.25 (See paper appendix Reanalyze)
end

function update_remote_counter!(remote_counter::RemoteChannel, update_count::Int)
	old_value=take!(remote_counter)
	put!(remote_counter, old_value + update_count)
	return nothing
end

@with_kw struct FeedForwardHP
	width_hidden::Int
	depth_representation::Int
	depth_prediction::Int
	depth_dynamics::Int
	depth_policy::Int
	depth_value::Int
	depth_reward::Int
	depth_state_head::Int
	use_batch_norm::Bool = false
	batch_norm_momentum::Float32 = 0.6f0
	hidden_state_size::Int #should be equal to prod(observation_shape)
	reward_activation
end

@with_kw mutable struct ResNetHP
	num_blocks::Int # Number of blocks in the ResNet
	depth_representation::Int
	num_filters::Int
	conv_kernel_size::Tuple{Int,Int}
	num_second_head_filters::Int = 2
	num_first_head_filters::Int = 1
	batch_norm_momentum::Float32 = 0.6f0
	downsample::Bool = false
	hidden_state_size::Int
	representation_output_size
	depth_policy::Int
	depth_value::Int
end


text:
Pkg.clone("https://github.com/Keno/AbstractTrees.jl")
Pkg.clone("https://github.com/JuliaComputing/Deprecations.jl")
Pkg.clone("https://github.com/JuliaComputing/FemtoCleaner.jl")

Pkg.update()

using FemtoCleaner
FemtoCleaner.cleanrepo("C:\\Users\\bernhard.konig\\Documents\\ASync\\home\\Code\\Julia\\DecisionTrees.jl"; show_diff=true, delete_local=false)


text:
#
# Methods for orbit states and state vectors.
#

Calculations.keplerian(state::CartesianStateVector, μ) = keplerian(position(state), velocity(state), μ)
Calculations.keplerian(orbit::CartesianR2BPOrbit) = keplerian(state(orbit), massparameter(system(orbit)))
Calculations.keplerian(orbit::KeplerianR2BPOrbit) = orbit

Calculations.cartesian(orbit::CartesianR2BPOrbit) = orbit
Calculations.cartesian(orbit::KeplerianR2BPOrbit) = cartesian(state(orbit), massparameter(system(orbit)))
Calculations.cartesian(state::KeplerianState, μ) = cartesian(
    eccentricity(state),
    semimajor_axis(state),
    inclination(state),
    RAAN(state),
    argument_of_periapsis(state),
    true_anomoly(state),
    μ
)

Base.position(state::CartesianStateVector) = States.get_r(state) * lengthunit(state)
Base.position(orbit::KeplerianR2BPOrbit) = position(CartesianState(cartesian(state(orbit), massparameter(state(orbit)))))
Base.position(orbit::Orbit) = position(state(orbit))

Calculations.distance(state::CartesianStateVector) = norm(position(state))
Calculations.distance(orbit::Orbit) = norm(position(orbit))

velocity(state::CartesianStateVector) = States.get_v(state) * velocityunit(state)
velocity(orbit::Orbit) = velocity(state(orbit))

Calculations.speed(state::CartesianStateVector) = norm(velocity(state))
Calculations.speed(orbit::KeplerianR2BPOrbit) = speed(semi_parameter(semimajor_axis(state(orbit)), eccentricity(state(orbit))), semimajor_axis(state(orbit)), massparameter(system(orbit)))
Calculations.speed(orbit::Orbit) = norm(velocity(state(orbit)))

Calculations.eccentricity(state::KeplerianState) = States.get_e(state)
Calculations.eccentricity(orbit::KeplerianR2BPOrbit) = eccentricity(state(orbit))
Calculations.eccentricity(orbit::CartesianR2BPOrbit) = eccentricity(position(state(orbit)), velocity(state(orbit)), massparameter(system(orbit)))

Calculations.semimajor_axis(state::KeplerianState) = States.get_a(state) * lengthunit(state)
Calculations.semimajor_axis(orbit::CartesianR2BPOrbit) = semimajor_axis(distance(state(orbit)), speed(state(orbit)), massparameter(system(orbit)))
Calculations.semimajor_axis(orbit::KeplerianR2BPOrbit) = semimajor_axis(state(orbit))

inclination(state::KeplerianState) = States.get_i(state) * angularunit(state)
inclination(orbit::KeplerianR2BPOrbit) = inclination(state(orbit))

RAAN(state::KeplerianState) = States.get_Ω(state) * angularunit(state)
RAAN(orbit::KeplerianR2BPOrbit) = RAAN(state(orbit))

argument_of_periapsis(state::KeplerianState) = States.get_ω(state) * angularunit(state)
argument_of_periapsis(orbit::KeplerianR2BPOrbit) = argument_of_periapsis(state(orbit))

Calculations.true_anomoly(state::KeplerianState) = States.get_ν(state) * angularunit(state)
Calculations.true_anomoly(orbit::KeplerianR2BPOrbit) = true_anomoly(state(orbit))
Calculations.true_anomoly(orbit::CartesianR2BPOrbit) = true_anomoly(distance(orbit), specific_angular_momentum(orbit), eccentricity(orbit), massparameter(system(orbit)))

Calculations.conic(orbit::R2BPOrbit) = conic(eccentricity(orbit))

Calculations.specific_angular_momentum_vector(orbit::CartesianR2BPOrbit) = specific_angular_momentum_vector(position(state(orbit)), velocity(state(orbit)))

Calculations.specific_angular_momentum(orbit::CartesianR2BPOrbit) = specific_angular_momentum(position(state(orbit)), velocity(state(orbit)))

Calculations.specific_energy(orbit::CartesianR2BPOrbit) = specific_energy(distance(state(orbit)), speed(state(orbit)), massparameter(system(orbit)))
Calculations.specific_energy(orbit::KeplerianR2BPOrbit) = specific_energy(semimajor_axis(orbit), massparameter(system(orbit)))

Calculations.C3(orbit::R2BPOrbit) = C3(distance(orbit), speed(orbit), massparameter(system(orbit)))

Calculations.v_infinity(orbit::R2BPOrbit) = v_infinity(distance(orbit), speed(orbit), massparameter(system(orbit)))

Calculations.eccentricity_vector(orbit::CartesianR2BPOrbit) = eccentricity_vector(position(state(orbit)), velocity(state(orbit)), massparameter(system(orbit)))

Calculations.semi_parameter(orbit::R2BPOrbit) = semi_parameter(semimajor_axis(orbit), eccentricity(orbit))

Calculations.periapsis_radius(orbit::R2BPOrbit) = periapsis_radius(semimajor_axis(orbit), eccentricity(orbit))

Calculations.apoapsis_radius(orbit::R2BPOrbit) = apoapsis_radius(semimajor_axis(orbit), eccentricity(orbit))

periapsis_velocity(orbit::R2BPOrbit) = speed(periapsis_radius(orbit), semimajor_axis(orbit), massparameter(system(orbit)))

apoapsis_velocity(orbit::R2BPOrbit) = speed(apoapsis_radius(orbit), semimajor_axis(orbit), massparameter(system(orbit)))

Calculations.period(orbit::R2BPOrbit) = conic(orbit) ∈ (Circular, Elliptical) ? period(semimajor_axis(orbit), massparameter(system(orbit))) : conic(orbit) == Parabolic ? Inf * timeunit(orbit) : NaN * timeunit(orbit)

Calculations.mean_motion(orbit::R2BPOrbit) = mean_motion(semimajor_axis(orbit), massparameter(system(orbit)))

function mean_motion_vector(orbit::R2BPOrbit)
#   î = SVector{3, Float64}([1, 0, 0])
#   ĵ = SVector{3, Float64}([0, 1, 0])
    k̂ = SVector{3}(0, 0, 1)
    return k̂ × specific_angular_momentum_vector(orbit)
end

Calculations.time_since_periapsis(orbit::R2BPOrbit) = time_since_periapsis(mean_motion(orbit), eccentricity(orbit), eccentric_anomoly(orbit))

Calculations.specific_potential_energy(orbit::CartesianR2BPOrbit) = specific_potential_energy(distance(orbit.state), massparameter(system(orbit)))

Calculations.jacobi_constant(orbit::CR3BPOrbit) = jacobi_constant(States.get_r(state(orbit)), States.get_v(state(orbit)), massparameter(system(orbit)))

Calculations.potential_energy(orbit::CR3BPOrbit) = potential_energy(States.get_r(state(orbit)), massparameter(system(orbit)))

distance_to_primary(orbit::CR3BPOrbit) = norm(States.get_r(state(orbit)) .- SVector(-massparameter(orbit), 0, 0))

distance_to_secondary(orbit::CR3BPOrbit) = norm(States.get_r(state(orbit)) .- SVector(1-massparameter(orbit), 0, 0))


text:
using GridWorld

g1 = Grid(3, 4, (x,y)->rand(1:9))


text:
function alterMini(w::Array{Float64,1}, shot::Shot, fidMtx::FidMtx, pos::Array{Int64,2}, mu::Float64, lambda::Float64; outer_it=20, fistait=10, CGit=5, tmax=0.5)
    nz = fidMtx.nz; nx = fidMtx.nx
    ntw = length(w)
    for k = 1 : outer_it
        println("outer loop: $k")
        dis = fista(shot, w, mu, lambda, fidMtx, pos_rec, tmax, fistait)
        (w, obj_w) = CG_w(ntw, shot, fidMtx, pos, dis, maxit=CGit, tmax=tmax)
    end
    return dis, w
end

function alterCG(w::Array{Float64,1}, shot::Shot, fidMtx::FidMtx, pos::Array{Int64,2}; outer_it=20, inner_it=5, tmax=0.5)
    nz = fidMtx.nz; nx = fidMtx.nx
    ntw = length(w)
    for k = 1 : outer_it
        println("outer loop: $k")
        (dis, obj_dis) = CG_dis(shot, fidMtx, pos, w, maxit=inner_it, tmax=tmax)
        (w  , obj_w  ) = CG_w(ntw, shot, fidMtx, pos, dis, maxit=inner_it, tmax=tmax)
    end
    return dis, w
end

function CG_w(ntw::Int64, shot::Shot, fidMtx::FidMtx, pos::Array{Int64,2}, dis::Array{Float64,1}; maxit=10, mu=0.1, tmax=0.5) # the default initial guess is zeros
    obj = zeros(maxit+1)
    x = zeros(ntw)
    r = copy(shot)
    s = MultiStepAdjoint(ntw, dis, r, fidMtx)
    p = copy(s)
    gamma = dot(s, s)
    obj[1] = gamma
    println("iteration 0, objective $gamma")
    for k = 1 : maxit
        q = MultiStepForward(pos, p, dis, fidMtx, tmax=tmax)
        delta = dot(vec(q.d), vec(q.d)) + mu*dot(p,p)
        alpha = gamma / delta
        x     = x + alpha*p
        r.d   = r.d - alpha*(q.d)
        s     = MultiStepAdjoint(ntw, dis, r, fidMtx) - mu*x
        gamma0= gamma
        gamma = dot(s,s)
        println("iteration $k, objective $gamma")
        obj[k+1] = gamma
        beta  = gamma / gamma0
        p     = s + beta*p
    end
    return x, obj
end

function CG_dis(shot::Shot, fidMtx::FidMtx, pos::Array{Int64,2}, w::Array{Float64,1}; maxit=10, mu=0.1, tmax=0.5) # the default initial guess is zeros
    obj = zeros(maxit+1)
    N = fidMtx.nz * fidMtx.nx
    x = zeros(N)
    r = copy(shot)
    s = MultiStepAdjoint(w, r, fidMtx)
    p = copy(s)
    gamma = dot(s, s)
    obj[1] = gamma
    println("iteration 0, objective $gamma")
    for k = 1 : maxit
        q = MultiStepForward(pos, w, p, fidMtx, tmax=tmax)
        delta = dot(vec(q.d), vec(q.d)) + mu*dot(p,p)
        alpha = gamma / delta
        x     = x + alpha*p
        r.d   = r.d - alpha*(q.d)
        s     = MultiStepAdjoint(w, r, fidMtx) - mu*x
        gamma0= gamma
        gamma = dot(s,s)
        println("iteration $k, objective $gamma")
        beta  = gamma / gamma0
        obj[k+1] = gamma
        p     = s + beta*p
    end
    return x, obj
end

max_stars_repo_name: WneleiGao/AcousticWaveLab

max_issues_repo_name: WneleiGao/AcousticWaveLab

max_forks_repo_name: WneleiGao/AcousticWaveLab


text:
load_file1 = joinpath("src", "algorithms", "smithwaterman.jl")
source_dir = pwd()

include(joinpath(source_dir, "../src", "algorithms", "smithwaterman.jl"));
include(joinpath(source_dir, "../src", "utils", "parse.jl"));
include(joinpath(source_dir, "../src", "utils", "benchmarking.jl"));

# include(normpath(joinpath(Base.source_dir(), "..", "src", "utils", "benchmarking.jl")))
# include(normpath(joinpath(Base.source_dir(),"..",load_file)))

using Base.Test
using Compat
import Compat: UTF8String, ASCIIString, view

function test_tpr()

    scores = [1,2,3,4,5,6,7,8,9,10]
    thresh = find_tp_threshold(scores, .7)

    @test thresh < 4

    scores = [.1, .2, .3, .4, .5, .6, .7, 8., .9, .10]
    tresh = find_tp_threshold(scores, .7)

    @test thresh < .4

    return

end

function test_sw_score()

    seq1 = "PRTEINS";
    seq2 = "PRTWPSEIN";

    smat = parse_parse_score_matrix(joinpath(source_dir, "../scoring", "BLOSUM50"))

    score = sw_align(seq1, seq2, smat, 1, 1, "test.txt"; write_align = false);

    @test score = 37.0

    score = sw_align(seq1, seq2, smat, 5, 3, "test.txt"; write_align = false);

    @test score = 29.0

    return

end

println("Running tests:")
test_tpr();


text:
module WFLearning

using Compat
using Yao, Kernels, Utils, Yao.Zoo, GradOptim, UnicodePlots, JLD2

################################################################################
#                           Define Circuit

using Yao, Yao.Blocks, Yao.LuxurySparse, Yao.Zoo
import Yao.Blocks: blocks, apply!, print_block
import Base: gradient
import Kernels: loss

"""
    There are many born machines, but this one is mine.
"""
struct OhMyBM{N, NL, CT, T, Basis <: AbstractBlock} <: CompositeBlock{N, T}
    circuit::CT
    basis::Basis

    function OhMyBM{N, NL}(pairs) where {N, NL}
        circuit = diff_circuit(N, NL, pairs)
        basis = rot_basis(N)
        new{N, NL, typeof(circuit), datatype(circuit), typeof(basis)}(circuit, basis)
    end
end

# Basic Components
entangler(pairs) = chain(control(ctrl, target=>X) for (ctrl, target) in pairs)
layer(x::Symbol) = layer(Val(x))
layer(::Val{:first}) = rollrepeat(chain(Rx(0.0), Rz(0.0)))
layer(::Val{:last}) = rollrepeat(chain(Rz(0.0), Rx(0.0)))
layer(::Val{:mid}) = rollrepeat(chain(Rz(0.0), Rx(0.0), Rz(0.0)))

# This is just some shortcuts
apply_zero(c) = apply!(zero_state(nqubits(c)), c)
initialize!(c) = dispatch!(c, 2pi * rand(nparameters(c)))
change_basis!(mine) = dispatch!(mine.basis, vec(randpolar(nqubits(mine))))

# Right, let's forward some APIs
function apply!(r::AbstractRegister, mine::OhMyBM)
    apply!(r, mine.circuit)
    apply!(r, mine.basis)
end

blocks(mine::OhMyBM) = blocks(mine.circuit)

function print_block(io::IO, qcbm::OhMyBM)
    printstyled(io, "OhMyBM"; bold=true, color=:red)
end

################################################################################
#                               Gradients
function gradient(c::OhMyBM, rots::Vector, prob, kernel, ptrain)
    map(gate->gradient(c, gate, prob, kernel, ptrain), rots)
end

function gradient(c::OhMyBM, gate, prob, kernel, ptrain)
    dispatch!(+, gate, pi / 2)
    prob_pos = probs(apply_zero(c))

    dispatch!(-, gate, pi)
    prob_neg = probs(apply_zero(c))

    dispatch!(+, gate, pi / 2)

    grad_pos = Kernels.expect(kernel, prob, prob_pos) - Kernels.expect(kernel, prob, prob_neg)
    grad_neg = Kernels.expect(kernel, ptrain, prob_pos) - Kernels.expect(kernel, ptrain, prob_neg)
    grad_pos - grad_neg
end

loss(c::OhMyBM, kernel, ptrain) = Kernels.loss(probs(apply_zero(c)), kernel, ptrain)

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

function train!(c::OhMyBM, psi, optim; learning_rate=0.1, epochs=200, nbatch=10)
    initialize!(c)
    kernel = Kernels.RBFKernel(nqubits(c), [0.25], false)
    history = Float64[]

    itr = Iterators.filter(x->(hasparameter(x) && isprimitive(x)), BlockTreeIterator(:DFS, c))
    rots = collect(itr)
    batch_grad = zeros(nparameters(c))
    for i = 1:epochs
        fill!(batch_grad, 0)
        curr_loss = 0

        for m = 1:nbatch
            change_basis!(c)
            ptrain = probs(apply!(copy(psi), c.basis))
            prob = probs(apply_zero(c))
            curr_loss += loss(c, kernel, ptrain)
            grad = gradient(c, rots, prob, kernel, ptrain)
            batch_grad += grad
        end
        batch_grad ./= nbatch
        curr_loss /= nbatch
        push!(history, curr_loss)
        println("epoch ", i, " loss = ", curr_loss)

        params = collect(parameters(c))
        update!(params, batch_grad, optim)
        dispatch!(c, params)
    end
    history
end

struct TrainData
    n::Int
    nlayers::Int
    nbatch::Int
    epochs::Int
    learning_rate::Float64
    history::Vector{Float64}
    fidelity::Float64
end

function task(filename, n, nlayers, nbatch, epochs, learning_rate)
    # set up machine
    bm = OhMyBM{n, nlayers}([(i%n + 1)=>((i+1)%n + 1) for i = 1:n])
    initialize!(bm)
    # set up optimizer
    optim = Adam(lr=0.1);

    # set target state
    # r = register(bit"0"^n) + register(bit"1"^n)
    r = rand_state(n)
    normalize!(r)
    his = train!(bm, r, optim, nbatch=nbatch, epochs=epochs, learning_rate=learning_rate)
    fs = fidelity(apply!(zero_state(nqubits(bm)), bm.circuit), r)
    data = TrainData(n, nlayers, nbatch, epochs, learning_rate, his, fs[1])
end

function ptask(filename, n, nlayers, nbatch, epochs, learning_rate; nsamples=100)
    data = pmap(x->task(filename, n, nlayers, nbatch, epochs, learning_rate), 1:nsamples)

    jldopen("data/$filename.jld2", "a+") do file
        file["$n/$nlayers/$nbatch/$epochs/$learning_rate"] = data
    end
end

export task, ptask

end

max_stars_repo_name: Roger-luo/QCBM

max_issues_repo_name: Roger-luo/QCBM

max_forks_repo_name: Roger-luo/QCBM


text:
"""
    AbstractΔSizeStrategy

Abstract base type for strategies for how to change the size.

See [`Δsize!`](@ref).
"""
abstract type AbstractΔSizeStrategy end

"""
    AbstractJuMPΔSizeStrategy <: AbstractΔSizeStrategy

Abstract type for strategies to change or align the sizes of vertices using JuMP.
"""
abstract type AbstractJuMPΔSizeStrategy <: AbstractΔSizeStrategy end

"""
    DecoratingJuMPΔSizeStrategy <: AbstractJuMPΔSizeStrategy

Abstract type for [`AbstractJuMPΔSizeStrategies`](@ref AbstractJuMPΔSizeStrategy) which wants to delegate some parts of the problem formulation to another strategy.

More concretely: If `s` is a [`DecoratingJuMPΔSizeStrategy`](@ref) then [`base(s)`](@ref base(s::DecoratingJuMPΔSizeStrategy)) will be used unless explicitly stated through dispatch.
"""
abstract type DecoratingJuMPΔSizeStrategy <: AbstractJuMPΔSizeStrategy end

"""
    base(s::DecoratingJuMPΔSizeStrategy)

Return the strategy wrapped by `s`.
"""
base(s::DecoratingJuMPΔSizeStrategy) = s.base # Just make a wild guess rather than throwing an error

"""
    ΔSizeFailError <: Exception

Size change could not be solved.
"""
struct ΔSizeFailError <: Exception
    msg::String
end
Base.showerror(io::IO, e::ΔSizeFailError) = print(io, "ΔSizeFailError: ", e.msg)

"""
    ThrowΔSizeFailError <: AbstractJuMPΔSizeStrategy
    ThrowΔSizeFailError(msg::String)

Throws an `ΔSizeFailError` with message `msg`.
"""
struct ThrowΔSizeFailError{F} <: AbstractJuMPΔSizeStrategy
    msgfun::F
end
ThrowΔSizeFailError(msg::AbstractString) = ThrowΔSizeFailError(vs -> msg)
ThrowΔSizeFailError() = ThrowΔSizeFailError(vs -> string("Could not change size of vertices ", join(nameorrepr.(vs), ", ", " and "),"!"))

add_participants!(::ThrowΔSizeFailError, vs=AbstractVertex[]) = vs

"""
    ΔSizeFailNoOp <: AbstractJuMPΔSizeStrategy
    ΔSizeFailNoOp()

Does not perform any action.
"""
struct ΔSizeFailNoOp <: AbstractJuMPΔSizeStrategy end

add_participants!(::ΔSizeFailNoOp, vs=AbstractVertex[]) = vs

"""
    LogΔSizeExec <: AbstractJuMPΔSizeStrategy
    LogΔSizeExec(msgfun)
    LogΔSizeExec(msgfun, level::Logging.LogLevel)
    LogΔSizeExec(msgfun, level::Logging.LogLevel, andthen::AbstractJuMPΔSizeStrategy)

Logs `msgfun(v)` at log level `level`, then executes [`AbstractJuMPΔSizeStrategy`](@ref) `andthen` for vertex `v`.
"""
struct LogΔSizeExec{F, S} <: DecoratingJuMPΔSizeStrategy
    msgfun::F
    level::LogLevel
    andthen::S
end
LogΔSizeExec(msgfun, level=Logging.Info) = LogΔSizeExec(msgfun, level, ΔSizeFailNoOp())
LogΔSizeExec(msg::String, level::LogLevel=Logging.Info, andthen=ΔSizeFailNoOp()) = LogΔSizeExec(v -> msg, level, andthen)
LogSelectionFallback(nextstr, andthen; level=Logging.Warn) = LogΔSizeExec(v -> "Size change for vertex $(nameorrepr(v)) failed! $nextstr", level, andthen)
base(s::LogΔSizeExec) = s.andthen
fallback(s::LogΔSizeExec) = fallback(base(s))

add_participants!(s::LogΔSizeExec, vs=AbstractVertex[]) = add_participants!(s.andthen, vs)

"""
    DefaultJuMPΔSizeStrategy <: AbstractJuMPΔSizeStrategy
    DefaultJuMPΔSizeStrategy()

Default strategy which applies default constraints or objectives to a model.

Intended use is for strategies which apply special constraints or objectives to a subset of the vertices.

Also useful when configuring [`DecoratingJuMPΔSizeStrategies`](@ref DecoratingJuMPΔSizeStrategy).

### Examples

Use to get add special constraints for a subset of vertices:
```julia
function NaiveNASlib.addsomeconstraint!(s::MySpecialStrategy, v::AbstractVertex, data)
    if v === s.targetvertex
        # Add some special constraints here
    else
        # Use the default
        NaiveNASlib.addsomeconstraint!(DefaultJuMPΔSizeStrategy(), v::AbstractVertex, data)
    end
end
```
Use with a [`DecoratingJuMPΔSizeStrategy`](@ref):
```jldoctest
julia> using NaiveNASlib, NaiveNASlib.Advanced

julia> WithUtilityFun(defaultutility, DefaultJuMPΔSizeStrategy())
WithUtilityFun{typeof(defaultutility), DefaultJuMPΔSizeStrategy}(NaiveNASlib.defaultutility, DefaultJuMPΔSizeStrategy())
```
"""
struct DefaultJuMPΔSizeStrategy <: AbstractJuMPΔSizeStrategy end
fallback(::DefaultJuMPΔSizeStrategy) = default_sizefallback("Solver failed!")

"""
    TimeLimitΔSizeStrategy{S} <: DecoratingJuMPΔSizeStrategy
    TimeLimitΔSizeStrategy(limit::Number)
    TimeLimitΔSizeStrategy(limit::Number, base::S)

Sets the time limit for the solver to `limit`. Use strategy `base` for all other aspects.
"""
struct TimeLimitΔSizeStrategy{S,F} <: DecoratingJuMPΔSizeStrategy
    limit::Float64
    base::S
    fallback::F
end
TimeLimitΔSizeStrategy(limit::Number, base=DefaultJuMPΔSizeStrategy(); fallback=default_sizefallback("Solver failed!")) = TimeLimitΔSizeStrategy(limit, base, fallback)
TimeLimitΔSizeStrategy(limit::Number, base::S, fallback::F) where {S,F} = TimeLimitΔSizeStrategy{S,F}(Float64(limit), base, fallback)
base(s::TimeLimitΔSizeStrategy) = s.base
fallback(s::TimeLimitΔSizeStrategy) = s.fallback
add_participants!(s::TimeLimitΔSizeStrategy, vs=AbstractVertex[]) = add_participants!(base(s), vs)

"""
    TimeOutAction{S,A,F} <: DecoratingJuMPΔSizeStrategy
    TimeOutAction(action::A, base::S, fallback::F)
    TimeOutAction(base; fallback)

Calls `action(model)` if JuMP model `model` has status `MOI.TIME_LIMIT` after optimization stops. Use strategy `base` for all other aspects.

Default action is to display a warning and then apply `fallback` (default [`ΔSizeFailNoOp`](@ref)).
"""
struct TimeOutAction{S,A,F} <: DecoratingJuMPΔSizeStrategy
    action::A
    base::S
    fallback::F
end
base(s::TimeOutAction) = s.base
function TimeOutAction(;action=default_timeout, base::AbstractJuMPΔSizeStrategy=DefaultJuMPΔSizeStrategy(), fallback=default_sizefallback("Solver failed!"))
    TimeOutAction(action, base, fallback)
end
fallback(s::TimeOutAction) = s.fallback
function default_timeout(args...)
    @warn "Solver timed out"
    return false
end
add_participants!(s::TimeOutAction, vs=AbstractVertex[]) = add_participants!(base(s), vs)

struct Exact end
struct Relaxed end

"""
    ΔNout <: AbstractJuMPΔSizeStrategy
    ΔNout{T, V, F}(Δs::Dict{V, Int}, fallback::F)

Strategy for changing nout of vertices.

For each key-value pair in `v, Δ` in `Δs`, change [`nout(v)`](@ref) by `Δ`, i.e new size is `nout(v) + Δ`.

If `T == Exact`, size change will be added as a constraint to the model which means that the operation will fail
if it is not possible to change `nout(v)` by exactly `Δ`.

If `T == Relaxed`, size change will be added as an objective to the model which means that `nout(v)` might
not change by exactly `Δ`. In addition, a constraint that `nout(v)` must change is also added.
"""
struct ΔNout{T, V, F} <: AbstractJuMPΔSizeStrategy
    Δs::Dict{V, Int}
    fallback::F
end
fallback(s::ΔNout) = s.fallback

add_participants!(s::ΔNout, vs=AbstractVertex[]) = append!(vs, filter(v -> v ∉ vs, all_in_Δsize_graph(keys(s.Δs), Output())))

generic_Δnin_docstring_examples(fname::String; kwargs...) = generic_Δnin_docstring_examples(args ->  "$fname($args)"; kwargs...)
generic_Δnin_docstring_examples(fgen; footer = "```", kwargs...) = generic_Δnout_docstring_examples(fgen; footer="", kwargs...) * """

julia> v3 = conc(affine(iv, 4), affine(iv, 5), affine(iv, 6); dims=2);

julia> v4 = conc(affine(iv, 7), affine(iv, 8); dims=2);

julia> $(fgen("v3, 3, missing, 2"));

julia> $(fgen("v3 => (3, missing, 2), v4 => (-2, 0)"));
""" * footer

generic_Δnout_docstring_examples(fname::String; kwargs...) = generic_Δnout_docstring_examples(args -> "$fname($args)"; kwargs...)
generic_Δnout_docstring_examples(fgen; header="```julia-repl", footer = "```") = """
### Examples
$header
julia> using NaiveNASlib, NaiveNASlib.Extend, NaiveNASlib.Advanced

julia> module TestNNLib
           using NaiveNASlib, NaiveNASlib.Extend
           mutable struct Affine{T} W::Matrix{T} end
           NaiveNASlib.nout(a::Affine) = size(a.W, 1)
           NaiveNASlib.nin(a::Affine) = [size(a.W, 2)]
           NaiveNASlib.Δsize!(a::Affine, ins::AbstractVector, outs::AbstractVector) = a.W = parselect(a.W, 2=>ins[1], 1=>outs)
           affine(in, outsize) = absorbvertex(Affine(ones(outsize, nout(in))), in)
           export affine
       end;

julia> using .TestNNLib;

julia> iv = affine(inputvertex("in", 3), 6);

julia> v1 = affine(iv, 4);

julia> v2 = affine(iv, 5);

julia> $(fgen("v1, 3"));

julia> $(fgen("v1 => -3, v2 => -2"));

julia> $(fgen("Dict(v1 => 3, v2 => 2)"));
""" * footer

"""
    ΔNout{T}(args...;fallback::F)

Return a `ΔNout{T, V, F}` created from `args` where `args` can be either a `Dict` or arguments used to create a `Dict`
(e.g. a set of `Pair`s or a generator).

Also helps creating a `Dict{V, Int}` from e.g. a `Dict{Any, Any}` as a convenience.

`T` is typically either `Exact` or `Relaxed`.

$(generic_Δnout_docstring_examples(args -> "ΔNout{NaiveNASlib.Exact}($args; fallback=ΔSizeFailNoOp())"; header="```jldoctest"))
"""
ΔNout{T}(v::AbstractVertex, Δ::Integer; fallback) where T = ΔNout{T}(v=>Int(Δ); fallback)
ΔNout{T}(args...; fallback) where T = ΔNout{T}(Dict(args...); fallback)
ΔNout{T}(Δs::AbstractDict; fallback) where T = ΔNout{T}(Dict(k=>Int(v) for (k,v) in Δs); fallback)
ΔNout{T}(Δs::AbstractDict{V, Int}; fallback::F) where {T,V,F} = ΔNout{T,V,F}(Dict(Δs), fallback)

"""
    ΔNoutExact(args...; fallback)

Return a `ΔNout{Exact}` with `fallback` set to give a warning and then try again with [`ΔNoutRelaxed(args...)`](@ref).

Accepts either a `Dict` or arguments used to create a `Dict` (e.g. a set of `Pair`s or a generator).

$(generic_Δnout_docstring_examples("ΔNoutExact"; header="```jldoctest"))
"""
ΔNoutExact(args...; fallback=default_noutfallback(ΔNoutRelaxed,"nout", args)) = ΔNout{Exact}(args...;fallback)

"""
    ΔNoutRelaxed(args...;fallback)

Return a `ΔNout{Relaxed}` with `fallback` set to [`ΔSizeFailNoOp`](@ref).

Accepts either a `Dict` or arguments used to create a `Dict` (e.g. a set of `Pair`s or a generator).

$(generic_Δnout_docstring_examples("ΔNoutRelaxed"; header="```jldoctest"))
"""
ΔNoutRelaxed(args...;fallback=default_noutfallback("nout", args)) = ΔNout{Relaxed}(args...;fallback)

const Maybe{T} = Union{T, Missing}

Δnout_err_info(v, Δ::Union{<:Maybe{Int}, <:Pair}...) = "$(nameorrepr(v)) by $(join(first.(Δ), ", ", " and "))"
Δnout_err_info(v, Δ::Tuple) = Δnout_err_info(v, Δ...)
Δnout_err_info(ps::Pair...) = join(map(p ->Δnout_err_info(p...), ps), ", ", " and ")
Δnout_err_info(d::AbstractDict) = join((Δnout_err_info(k, v) for (k,v) in d), ", ", " and ")

function default_noutfallback(nextfallback, dirstr, args)
    msgfun = v -> string("Could not change ", dirstr, " of ", Δnout_err_info(args...), "! Relaxing constraints...")
    return LogΔSizeExec(msgfun, Logging.Warn, nextfallback(args...))
end

default_noutfallback(dirstr, args) = LogΔSizeExec(vs -> string("Could not change ", dirstr, " of ", Δnout_err_info(args...), "! Attempt aborted!"), Logging.Warn)
default_sizefallback(msg) = LogΔSizeExec(msg, Logging.Warn)

"""
    ΔNinExact(args...; fallback)

Return a `ΔNout{Exact}` configured to set `nout` of the inputs to vertices in `args` to the given sizes.

For vertices with more than one input, the size change must be expressed as a tuple with one element per input.
Use `missing` to indicate that no special treatment is needed for an input.

Accepts either a `Dict` or arguments used to create a `Dict` (e.g. a set of `Pair`s or a generator).

By default, `fallback` is set to give a warning and then try again with [`ΔNinRelaxed(args...)`](@ref).

$(generic_Δnin_docstring_examples("ΔNinExact"; header="```jldoctest"))
"""
ΔNinExact(args...; fallback=default_noutfallback(ΔNinRelaxed, "nin", args)) = ΔNout{Exact}(Δnin2Δnout(args...); fallback)

"""
    ΔNinRelaxed(args...; fallback)

Same as [`ΔNinExact`](@ref) except `Exact` is replaced by `Relaxed` and `fallback` set to [`ΔSizeFailNoOp`](@ref) by default.

$(generic_Δnin_docstring_examples("ΔNinRelaxed"; header="```jldoctest"))
"""
ΔNinRelaxed(args...; fallback=default_noutfallback("nout", args)) = ΔNout{Relaxed}(Δnin2Δnout(args...); fallback)

Δnin2Δnout(d::AbstractDict) = reduce(merge!, (Δnin2Δnout(k,v) for (k,v) in d); init=Dict())
Δnin2Δnout(ps::Pair{<:AbstractVertex}...) = mapreduce(p -> Δnin2Δnout(p...), merge!, ps; init=Dict())
Δnin2Δnout(v::AbstractVertex, Δs...) = Δnin2Δnout(v, Δs)
Δnin2Δnout(v::AbstractVertex, Δs::Pair{<:Tuple, <:Union{Relaxed, Exact}}) = Dict(k => (v => last(Δs)) for (k,v) in Δnin2Δnout(v, first(Δs)))
function Δnin2Δnout(v::AbstractVertex, Δs::Tuple)
    @assert length(Δs) == length(inputs(v)) "Must supply same number of Δs as v has inputs! Got $Δs for $(nameorrepr(v)) with $(length(inputs(v))) inputs.
    $(length(Δs) < length(inputs(v)) ? "Tip: Use missing to indicate that no special requirement shall be used for an input" : "")"

    inds = findall(!ismissing, Δs)
    return Dict(inputs(v)[i] => Δs[i] for i in inds)
end

"""
    ΔNoutMix{VE, VR, F} <: AbstractJuMPΔSizeStrategy
    ΔNoutMix(exact::ΔNout{Exact, VE, F}, relax::ΔNout{Relaxed, VR, F}, fallback::F)

Strategy for changing nout of vertices.

Applies both `exact` and `relax`, thereby allowing for having both hard and relaxed size constraints for different vertices.

Can be conveniently created by [`ΔNout`](@ref) or [`ΔNin`](@ref).
"""
struct ΔNoutMix{VE, VR, F} <: AbstractJuMPΔSizeStrategy
    exact::ΔNout{Exact, VE, F}
    relax::ΔNout{Relaxed, VR, F}
    fallback::F
end
fallback(s::ΔNoutMix) = s.fallback

function add_participants!(s::ΔNoutMix, vs=AbstractVertex[])
    add_participants!(s.exact, vs)
    add_participants!(s.relax, vs)
    return vs
end

"""
    relaxed(Δ)

Return `Δ => Relaxed()` which indicates that `Δ` shall be relaxed when changing size.

See [`Δnout!`](@ref) and [`Δnin!`](@ref).
"""
relaxed(Δ::Integer) = Δ => Relaxed()
relaxed(Δs::Tuple{Vararg{<:Maybe{Int}}}) = Δs => Relaxed()
relaxed(Δs...) = Δs => Relaxed()

"""
    ΔNout(args...;[fallback])

Splits `args` into relaxed and exact size changes and creates the appropriate strategy (one of `ΔNout{Exact}, ΔNout{Relaxed} or ΔNoutMix`).

Fallback strategy may be supplied through the `fallback` keyword. Defaults to `ΔNoutRelaxed` for all given vertices.

$(generic_Δnout_docstring_examples("ΔNout"; header="```jldoctest", footer=""))
julia> ΔNout(v1, relaxed(2));

julia> ΔNout(v1 => relaxed(2), v2 => -1);
```
"""
function ΔNout(args...;kwargs...)
    exact, relaxed = split_exact_relaxed(args)
    return ΔNout(exact, relaxed; kwargs...)
end

"""
    ΔNin(args...;[fallback])

Splits `args` into relaxed and exact size changes and creates the appropriate strategy (one of `ΔNout{Exact}, ΔNout{Relaxed} or ΔNoutMix`).

Fallback strategy may be supplied through the `fallback` keyword. Defaults to `ΔNinRelaxed` for all given vertices.

$(generic_Δnin_docstring_examples("ΔNin"; header="```jldoctest", footer=""))
julia> ΔNin(v3, relaxed(3), missing, 2);

julia> ΔNin(v3 => (relaxed(3), missing, 2), v4 => relaxed(-2, 0));
```
"""
function ΔNin(args...;fallback=nothing)
    exact,relaxed = split_exact_relaxed(Δnin2Δnout(args...)) # To support mixed cases, e.g. Δnin!(v => (2, relaxed(4)))
    fallback = fallback === nothing ? default_noutfallback((args...) -> ΔNoutRelaxed(merge(exact, relaxed); fallback=default_noutfallback("nin", args)), "nin", args) : fallback
    return ΔNout(exact, relaxed; fallback)
end

function ΔNout(exact::Dict, relaxed::Dict; fallback = default_noutfallback(ΔNoutRelaxed, "nout", merge(exact,relaxed)))
    isempty(relaxed) && return ΔNoutExact(exact;fallback)
    isempty(exact) && return ΔNoutRelaxed(relaxed; fallback=_relaxed_fallback(relaxed, fallback))
    return ΔNoutMix(ΔNoutExact(exact;fallback), ΔNoutRelaxed(relaxed;fallback), fallback)
end

_relaxed_fallback(Δs, fb) = _is_same_fallback(Δs, fb) ? fallback(fb) : fb

_is_same_fallback(Δs, fb::LogΔSizeExec) = _is_same_fallback(Δs, base(fb))
_is_same_fallback(Δs, fb::ΔNout{Relaxed}) = fb.Δs == Δs
_is_same_fallback(Δs, fb) = false

split_exact_relaxed(p::Tuple{AbstractVertex, Any}) = split_exact_relaxed(tuple(Pair(p...)))
split_exact_relaxed(d::Tuple{Dict}) = split_exact_relaxed(first(d))
function split_exact_relaxed(args::Union{AbstractDict, Tuple{Vararg{Pair}}})
    exact=Dict()
    relaxed=Dict()
    for (k, v) in args
        if !isa(v, Pair) || v isa Pair{<:Any, Exact}
            exact[k] = first(v) # first works for numbers too
        elseif v isa Pair{<:Any, Relaxed}
            relaxed[k] = first(v)
        else
            throw(ArgumentError(string("Unhandled value type: ", v)))
        end
    end
    return exact, relaxed
end

"""
    AlignNinToNout{S, F} <: DecoratingJuMPΔSizeStrategy
    AlignNinToNout(;vstrat::S, fallback::F)
    AlignNinToNout(vstrat::S, fallback::F)

Adds variables and constraints for `nin(vi) == nout.(inputs(vi))`.

Generally useful when doing structural changes such are removing/adding vertices or edges.

If it fails, the operation will be retried with the `fallback` strategy (default [`ΔSizeFailNoOp`](@ref)).
"""
struct AlignNinToNout{S, F} <: DecoratingJuMPΔSizeStrategy
    nindict::Dict{AbstractVertex, Vector{JuMP.VariableRef}}
    vstrat::S
    fallback::F
end
AlignNinToNout(;vstrat=DefaultJuMPΔSizeStrategy(), fallback=default_sizefallback("Failed to align Nin to Nout!!")) = AlignNinToNout(vstrat, fallback)
AlignNinToNout(vstrat, fallback) = AlignNinToNout(Dict{AbstractVertex, Vector{JuMP.VariableRef}}(), vstrat, fallback)
fallback(s::AlignNinToNout) = s.fallback
base(s::AlignNinToNout) = s.vstrat

function add_participants!(s::AlignNinToNout, vs=AbstractVertex[])
    add_participants!(base(s), vs)
    append!(vs, filter(v -> v ∉ vs, all_in_Δsize_graph(keys(s.nindict), Both())))
end

"""
    AlignNinToNoutVertices{V1,V2,S,F} <: AbstractJuMPΔSizeStrategy
    AlignNinToNoutVertices(vin::V1, vout::V2, inds; vstrat::S=AlignNinToNout(), fallback::F=ΔSizeFailNoOp())
    AlignNinToNoutVertices(vin::V1, vout::V2, inds, vstrat::S, fallback::F)

Same as [`AlignNinToNout`](@ref) with an additional constraint that [`nin(s.vin)`](@ref)`[s.ininds] ==` [`nout(s.vout)`](@ref)
where `s` is a [`AlignNinToNoutVertices`](@ref).

Useful in the context of adding edges to vertices to align sizes before the edge has been added.

If it fails, the operation will be retried with the `fallback` strategy (default [`ΔSizeFailNoOp`](@ref)).
"""
struct AlignNinToNoutVertices{V1,V2,S,F} <: DecoratingJuMPΔSizeStrategy
    vin::V1
    vout::V2
    ininds::Vector{Int}
    vstrat::S
    fallback::F
end
AlignNinToNoutVertices(vin, vout, inds; vstrat=AlignNinToNout(), fallback=failtoalign(vin, vout)) = AlignNinToNoutVertices(vin, vout, inds, vstrat, fallback)
AlignNinToNoutVertices(vin, vout, inds, vstrat, fallback) = AlignNinToNoutVertices(vin, vout, collect(Int, inds), vstrat, fallback)

fallback(s::AlignNinToNoutVertices) = s.fallback
base(s::AlignNinToNoutVertices) = s.vstrat

function add_participants!(s::AlignNinToNoutVertices, vs=AbstractVertex[])
    add_participants!(base(s), vs)
    append!(v -> v ∉ vs, vcat(all_in_Δsize_graph(s.vin, Input()), all_in_Δsize_graph(s.vout, Output())))
end

failtoalign(vin, vout) = LogΔSizeExec(vs -> string("Could not align nout of ", nameorrepr(vin), " to nin of ", nameorrepr(vout), "!!"), Logging.Warn)

"""
    TruncateInIndsToValid{S} <: AbstractΔSizeStrategy
    TruncateInIndsToValid()
    TruncateInIndsToValid(s::S)

Ensures that all selected input indices are within range of existing input indices after applying `s`
(default [`DefaultJuMPΔSizeStrategy`](@ref)).

Not needed in normal cases, but certain structural mutations (e.g create_edge!) may cause this to happen
due to how constraints are (not) created when original sizes do not align in conjunction with how result of
selection is interpreted.

An example of when it is needed is when adding an edge from vertex `vi` to an invariant vertex `vo` where
[`nout(vi)`](@ref) `>` [`nout(vo)`](@ref). In this case it is expected that the result of the optimization
is that the indices `1:nout(vi)` of `vi` shall be kept. However, this will propagate to `vo` which will be
instructed to keep indices it does not have. With this strategy, all indices which are larger than `nout(vo)`
will be replaced by `-1` (which indicates that new parameters shall be created)

While this may be considered a flaw in the output selection procedure, it is rare enough so that in most cases
when it happens it is the result of a user error or lower level bug. Therefore this strategy is left optional
to be used only in cases when mismatches are expected.

Note that in most normal cases, this has no effect since vertices capable of getting new input edges generally
don't have parameters.
"""
struct TruncateInIndsToValid{S} <: AbstractΔSizeStrategy
    strategy::S
end
TruncateInIndsToValid() = TruncateInIndsToValid(DefaultJuMPΔSizeStrategy())
base(s::TruncateInIndsToValid) = s.strategy

add_participants!(s::TruncateInIndsToValid, vs=AbstractVertex[]) = add_participants!(base(s), vs)

"""
    WithUtilityFun{F, S} <: AbstractΔSizeStrategy
    WithUtilityFun(utilityfun::F, strategy::S)

Applies neuron indices selection with `strategy` and using `utilityfun` to compute the utility of neurons indices.

Note that `utilityfun` will override any utility function supplied in function call. Thus it is possible use
`WithUtilityFun` to change utility function e.g. when switching to a fallback strategy.
"""
struct WithUtilityFun{F, S} <: AbstractΔSizeStrategy
    utilityfun::F
    strategy::S
end
WithUtilityFun(f) = s -> WithUtilityFun(f ,s)
base(s::WithUtilityFun) = s.strategy
fallback(s::WithUtilityFun) = fallback(s.strategy)

add_participants!(s::WithUtilityFun, vs=AbstractVertex[]) = add_participants!(base(s), vs)

"""
    UpperInsertBound{V<:AbstractDict{<:AbstractVertex, <:Integer}, S, F} <: DecoratingJuMPΔSizeStrategy
    UpperInsertBound(vbounds; strategy, fallback)
    UpperInsertBound(vs, bound::Integer; strategy, fallback)

Sets upper bound for how many new neurons to insert of each vertex `v` in `vbounds` to `vbounds[v]` and then applies neuron indices
selection with `strategy` (default [`DefaultJuMPΔSizeStrategy`](@ref)).

Alternatively, sets the upper bound of each vertex in `vs` to `bound`. Applies `bound` to all vertices if `vs` is empty.

Mostly useful to minimize solver time in cases when one does not want any new neurons as disallowing this allows skipping a couple of constraints.

If it fails, the operation will be retried with the `fallback` strategy (default [`ΔSizeFailNoOp`](@ref)).
"""
struct UpperInsertBound{VB, S, F} <: DecoratingJuMPΔSizeStrategy
    vbounds::VB
    strategy::S
    fallback::F
end

function UpperInsertBound(vs::AbstractVector{<:AbstractVertex}, bound::Integer=0; kwargs...)
    isempty(vs) ? UpperInsertBound(bound; kwargs...) : UpperInsertBound(Dict(v => bound for v in vs); kwargs...)
end
function UpperInsertBound(vbounds::Union{Integer, AbstractDict}; strategy=DefaultJuMPΔSizeStrategy(), fallback=default_sizefallback("Solver failed!"))
    UpperInsertBound(vbounds, strategy, fallback)
end
base(s::UpperInsertBound) = s.strategy
fallback(s::UpperInsertBound) = s.fallback
add_participants!(s::UpperInsertBound, vs=AbstractVertex[]) = add_participants!(base(s), vs)

"""
    AbstractAfterΔSizeStrategy <: DecoratingJuMPΔSizeStrategy

Abstract base type for strategies which perform actions after size has changed (e.g validation and logging).
"""
abstract type AbstractAfterΔSizeStrategy <: DecoratingJuMPΔSizeStrategy end

"""
    WithKWargs{KW, S} <: AbstractAfterΔSizeStrategy
    WithKwargs(strategy; kwargs...)

Adds keyword argument `kwargs` to any calls to [`Δsize!`](@ref Δsize!(::Any, ::AbstractVector, ::AbstractVector)) after sizes have been determined.
"""
struct WithKwargs{KW, S} <: AbstractAfterΔSizeStrategy
    kwargs::KW
    strategy::S
end
WithKwargs(strategy; kwargs...) = WithKwargs(kwargs, strategy)
base(s::WithKwargs) = s.strategy
fallback(s::WithKwargs) = fallback(s.strategy)

add_participants!(s::WithKwargs, vs=AbstractVertex[]) = add_participants!(base(s), vs)

"""
    AfterΔSizeCallback{F, S} <: AbstractAfterΔSizeStrategy
    AfterΔSizeCallback(cbfun::F, basestrat::S=ThrowΔSizeFailError())

Calls `cbfun(v, Δ, isnout)` for all vertices which change size after having been asked to change their sizes as a result of `basestrat`.
"""
struct AfterΔSizeCallback{F, S} <: AbstractAfterΔSizeStrategy
    cbfun::F
    base::S
end
AfterΔSizeCallback(cbfun) = AfterΔSizeCallback(cbfun, ThrowΔSizeFailError())
base(s::AfterΔSizeCallback) = s.base
add_participants!(s::AfterΔSizeCallback, vs=AbstractVertex[]) = add_participants!(base(s), vs)
fallback(s::AfterΔSizeCallback) = fallback(base(s))

"""
    logafterΔsize(printfun=nameorrepr;level=Logging.Info, base=DefaultJuMPΔSizeStrategy())

Return an [`AfterΔSizeCallback`](@ref) configured to log size changes with log level `level`.

For a given vertex `v`, `printfun(v)` will be used in the logged string.

Strategy `base` will be used to change sizes (e.g if [`Δsize!(logafterΔsize(base))`](@ref) is called).
"""
logafterΔsize(printfun=nameorrepr;level=Logging.Info, base=DefaultJuMPΔSizeStrategy()) = AfterΔSizeCallback(base) do v, Δ, dirlabel, ischanged
    ischanged || return
    dstr, Δstr = if dirlabel === :nout
                        "nout", compressed_string(Δ)
                    else
                        "nin", join(compressed_string.(Δ), ", ", " and ")
                    end
    @logmsg level "Change $dstr of $(printfun(v)) by $Δstr"
end

"""
    validateafterΔsize(printfun=nameorrepr; base=DefaultJuMPΔSizeStrategy())

Return an [`AfterΔSizeCallback`](@ref) configured to validate that sizes (nin and nout) are consistent
after a size change and throw a `ΔSizeFailError` if validation fails.

For a given vertex `v`, `printfun(v)` will be used in the error message should the size validation fail.

Strategy `base` will be used to change sizes (e.g if [`Δsize!(validateafterΔsize(base))`](@ref) is called).
"""
validateafterΔsize(printfun=nameorrepr; base=DefaultJuMPΔSizeStrategy()) = AfterΔSizeCallback(base) do v, Δ, dirlabel, ischanged
    dirlabel === :nout && return validate_Δnout(printfun, v, Δ)
    dirlabel === :nin && return validate_Δnin(printfun, v, Δ)
end

function validate_Δnin(pf, v::AbstractVertex, Δ)
    length(Δ) == length(inputs(v)) || throw(ArgumentError("Length of Δ must be equal to number of inputs for $(pf(v))! length(Δ) = $(length(Δ)), length(inputs(v)) = $(length(inputs(v)))"))
    nout.(inputs(v)) == nin(v) || throw(ΔSizeFailError("Nin change of $(compressed_string.(Δ)) to $(pf(v)) did not result in expected size! Expected: $(nout.(inputs(v))), actual: $(nin(v))"))
end

function validate_Δnout(pf, v::AbstractVertex, Δ)
    nin_of_outputs = unique(mapreduce(vi -> nin(vi)[inputs(vi) .== v], vcat, outputs(v), init=nout(v)))
    nin_of_outputs == [nout(v)] || throw(ΔSizeFailError("Nout change of $(compressed_string(Δ)) to $(pf(v)) resulted in size mismatch! Nin of outputs: $nin_of_outputs, nout of this: $([nout(v)])"))
end


text:
module Data

import ..Flux

export CMUDict, cmudict

deps(path...) = joinpath(@__DIR__, "..", "..", "deps", path...)

function __init__()
  mkpath(deps())
end

include("mnist.jl")
export MNIST

include("fashion-mnist.jl")
export FashionMNIST

include("cmudict.jl")
using .CMUDict

include("tree.jl")
include("sentiment.jl")
using .Sentiment

end


text:
function next_time_period(TIME::Timer, QC::Array{QuayCrane, 1}, CTS::Constants)
    min_unit = 5
    TIME.period += min_unit
    l = Array{Int, 1}()
    #update all QC atributtes
    for qc in QC
        #update QC position
        if qc.status == "moving"
            if  qc.next_bay < qc.current_bay
                qc.current_bay -= min_unit/CTS.tt
            elseif qc.next_bay > qc.current_bay
                qc.current_bay += min_unit/CTS.tt
            end
            if abs(qc.current_bay-round(qc.current_bay)) < 1e-10
                qc.current_bay = round(qc.current_bay)
            end
        end
        #update QC time left and status
        if qc.time_left > min_unit
            qc.time_left -= min_unit
        elseif qc.time_left == min_unit
            if length(qc.task_buffer) > 0
                if qc.status == "loading"
                    if qc.current_bay == qc.task_buffer[1].b
                        qc.time_left = 2*qc.task_buffer[1].t
                        qc.status = "loading"
                        deleteat!(qc.task_buffer, 1)
                    else
                        qc.time_left = travel_time(qc.current_bay, qc.task_buffer[1].b, CTS)
                        qc.next_bay = qc.task_buffer[1].b
                        qc.status = "moving"
                    end
                elseif qc.status == "moving" || qc.status == "waiting to load"
                    if qc.task_buffer[1].t ==0
                        qc.status = "idle"
                    else
                        qc.status = "loading"
                    end
                    qc.time_left = 2*qc.task_buffer[1].t
                    deleteat!(qc.task_buffer, 1)
                elseif qc.status == "waiting to move"
                    qc.time_left = travel_time(qc.current_bay, qc.task_buffer[1].b, CTS)
                    qc.next_bay = qc.task_buffer[1].b
                    qc.status = "moving"
                end
            else
                qc.time_left = 0
                qc.status = "idle"
                qc.task_buffer = Array{LTask, 1}()
            end
        elseif qc.time_left == 0
            qc.time_left = 0
            qc.status = "idle"
        end
        if qc.status == "idle"
            push!(l, qc.q)
        end
    end
    TIME.available_cranes = l
end

function travel_time(current_bay::Number, target_bay::Int, CTS::Constants)
    return(round(abs(current_bay - target_bay)*CTS.tt))
end

function get_qc_last_time(q::Int, LS::LoadingSequence)
    for t in reverse(LS.order)
        if t.qc == q
            return(t.start_time + 2*t.task.t)
        end
    end
end

function get_qc_start_time(q::Int, LS::LoadingSequence)
    for t in LS.order
        if t.qc == q
            return(t.start_time)
        end
    end
end

function get_qc_last_bay(q::Int, LS::LoadingSequence)
    for t in reverse(LS.order)
        if t.qc == q
            return(t.task.b)
        end
    end
    return(0)
end
function get_qc_first_bay(q::Int, LS::LoadingSequence)
    for t in LS.order
        if t.qc == q
            return(t.task.b)
        end
    end
end

function update_quay_crane(TIME::Timer, QC::Array{QuayCrane, 1}, task::LTask, q::Int, CTS::Constants)
    if abs(QC[q].current_bay-round(QC[q].current_bay)) < 1e-10
        QC[q].current_bay = round(QC[q].current_bay)
    end
    if task.b == QC[q].current_bay
        QC[q].time_left = 2*task.t
        QC[q].status = "loading"
        QC[q].next_bay = QC[q].current_bay
    else
        QC[q].time_left = travel_time(QC[q].current_bay, task.b, CTS)
        QC[q].status = "moving"
        QC[q].next_bay = task.b
        if task.c != 0 && task.p != 0
            push!(QC[q].task_buffer, task)
        end
    end
end

function update_load_seq(TIME::Timer, LS::LoadingSequence, QC::Array{QuayCrane, 1}, task::LTask, q::Int, CTS::Constants)
    LS.len += 1
    LS.tasks_left -= 1
    start_time = TIME.period + travel_time(QC[q].current_bay, task.b, CTS)
    push!(LS.order, (task=task, start_time=start_time, qc=q))
    push!(LS.filled_pos, task.p)
    push!(LS.loaded_cont,task.c)
end

function update_wait_quay_crane(TIME::Timer, QC::Array{QuayCrane, 1}, task::LTask, q::Int, move_q::Int, move_bay::Int, CTS::Constants)
    if abs(QC[q].current_bay-round(QC[q].current_bay)) < 1e-10
        QC[q].current_bay = round(QC[q].current_bay)
    end
    if task.b == QC[q].current_bay
        QC[q].time_left = travel_time(QC[move_q].current_bay, QC[move_q].next_bay, CTS)
        QC[q].status = "waiting to load"
        QC[q].next_bay = QC[q].current_bay
        push!(QC[q].task_buffer, task)
    elseif travel_time(QC[q].current_bay, task.b, CTS) >= travel_time(QC[move_q].current_bay, QC[move_q].next_bay, CTS)
        QC[q].time_left = travel_time(QC[q].current_bay, task.b, CTS)
        QC[q].status = "moving"
        QC[q].next_bay = task.b
        if task.t != 0 && task.p !=0
            push!(QC[q].task_buffer, task)
        end
    else
        QC[q].time_left = travel_time(QC[move_q].current_bay, QC[move_q].next_bay, CTS) - travel_time(QC[q].current_bay, task.b, CTS)
        QC[q].status = "waiting to move"
        QC[q].next_bay = task.b
        if task.t != 0 && task.p !=0
            push!(QC[q].task_buffer, task)
        end
    end
end

function update_wait_load_seq(TIME::Timer, LS::LoadingSequence, QC::Array{QuayCrane, 1}, task::LTask, q::Int, move_q::Int, move_bay::Int, CTS::Constants)
    LS.len += 1
    LS.tasks_left -= 1
    start_time = TIME.period + max(travel_time(QC[q].current_bay, task.b, CTS), travel_time(QC[move_q].current_bay, QC[move_q].next_bay, CTS))
    push!(LS.order, (task=task, start_time=start_time, qc=q))
    push!(LS.filled_pos, task.p)
    push!(LS.loaded_cont,task.c)
end

function update_dummy_load_seq(TIME::Timer, LS::LoadingSequence, QC::Array{QuayCrane, 1}, move_q::Int, move_bay::Int, CTS::Constants)
    start_time = TIME.period + travel_time(QC[move_q].current_bay, move_bay, CTS)
    push!(LS.order, (task = LTask(0, move_bay, 0, 0), start_time=start_time, qc=move_q))
end

function add_task(task::LTask, q::Int, TIME::Timer, LS::LoadingSequence, QC::Array{QuayCrane, 1}, CTS::Constants)
    update_load_seq(TIME, LS, QC, task, q, CTS)
    update_quay_crane(TIME, QC, task, q, CTS)
    deleteat!(TIME.available_cranes, findall(x->x==q, TIME.available_cranes))
end

function add_task_move(task::LTask, q::Int, move_q::Int, move_bay::Int, move_status::String, LS::LoadingSequence, TIME::Timer, QC::Array{QuayCrane, 1}, CTS::Constants)
    if move_status == "idle"
        update_quay_crane(TIME, QC, LTask(0, move_bay, 0, 0), move_q, CTS)
        update_dummy_load_seq(TIME, LS, QC, move_q, move_bay, CTS)
        deleteat!(TIME.available_cranes, findall(x->x==move_q, TIME.available_cranes))
        update_wait_quay_crane(TIME, QC, task, q, move_q, move_bay, CTS)
        update_wait_load_seq(TIME, LS, QC, task, q, move_q, move_bay, CTS)
        deleteat!(TIME.available_cranes, findall(x->x==q, TIME.available_cranes))
    elseif move_status == "moving"
        update_wait_quay_crane(TIME, QC, task, q, move_q, move_bay, CTS)
        update_wait_load_seq(TIME, LS, QC, task, q, move_q, move_bay, CTS)
        deleteat!(TIME.available_cranes, findall(x->x==q, TIME.available_cranes))
    end
end

function add_move(move_q::Int, move_bay::Int, LS::LoadingSequence, TIME::Timer, QC::Array{QuayCrane, 1}, CTS::Constants)
    update_quay_crane(TIME, QC, LTask(0, move_bay, 0, 0), move_q, CTS)
    update_dummy_load_seq(TIME, LS, QC, move_q, move_bay, CTS)
    deleteat!(TIME.available_cranes, findall(x->x==move_q, TIME.available_cranes))
end

max_stars_repo_licenses: Zlib MIT

max_issues_repo_licenses: Zlib MIT

max_forks_repo_licenses: Zlib MIT


text:
"""

In this example, we plot the temperature, read from the log file, as function
of simulation time

"""

using M3GTools

println(" Loading Plots... ")
using Plots
ENV["GKSwstype"]="nul" # This supresses the need of a display while plotting

# Initialize simulation data and path to VMD executable
println(" Reading simulation data ... ")
mysim = M3GTools.init(psf="../simulation_files/structure.psf",
                  log="../simulation_files/example.log")

println(" Plotting... ")
plot(mysim.log.time,mysim.log.temperature,label="example.log")
plot!(xlabel="time / ns")
plot!(ylabel="temperature / K")

savefig("temperature.pdf")
println(" Created figure: temperature.pdf")


text:
@testset "Arithmetic" begin
    x = rand(Complex{Double64})
    y = rand(Double64)

    @test x + y isa Complex{Double64}
    @test y + x isa Complex{Double64}
    @test x - y isa Complex{Double64}
    @test y - x isa Complex{Double64}
    @test x * y isa Complex{Double64}
    @test y * x isa Complex{Double64}
    @test x / y isa Complex{Double64}
    @test y / x isa Complex{Double64}
    @test x ^ y isa Complex{Double64}
    @test y ^ x isa Complex{Double64}
    @test x ^ x isa Complex{Double64}
    @test y ^ y isa Double64

    x, y = rand(2)
    @test add2(x, y) == (Double64(x) + Double64(y))
    @test sub2(x, y) == (Double64(x) - Double64(y))
    @test mul2(x, y) == (Double64(x) * Double64(y))
    @test div2(x, y) == (Double64(x) / Double64(y))

    @test x ⊕ y == (Double64(x) + Double64(y))
    @test x ⊖ y == (Double64(x) - Double64(y))
    @test x ⊗ y == (Double64(x) * Double64(y))
    @test x ⊘ y == (Double64(x) / Double64(y))
end

@testset "div,rem,.." begin
    a = 17.0
    b =  5.0
    da = Double64(a)
    db = Double64(b)

    @test div(da, db) == div(a, b)
    @test cld(da, db) == cld(a, b)
    @test fld(da, db) == fld(a, b)
    @test rem(da, db) == rem(a, b)
    @test mod(da, db) == mod(a, b)
    @test divrem(da, db) == divrem(a, b)

    @test fldmod(da, db) == fldmod(a, b)
end

@testset "Trig functions" begin
    x = rand(Double64)
    sinx, cosx = sincos(x)
    @test sinx ≈ sin(x)
    @test cosx ≈ cos(x)
end

@testset "promotion" begin
    @test df64"123" - df64"123" * df64"0.1" * df64"10.0" == 0.0
    @test df64"123" - 123 * df64"0.1" * df64"10.0" == 0.0
end

@testset "modpi" begin
    @test iszero(DoubleFloats.mod1pi(Double64(pi)))
    @test iszero(DoubleFloats.mod2pi(2*Double64(pi)))
    @test iszero(DoubleFloats.modhalfpi(0.5*Double64(pi)))
    @test iszero(DoubleFloats.modqrtrpi(0.25*Double64(pi)))

    @test iszero(DoubleFloats.mod1pi(Double32(pi)))
    @test iszero(DoubleFloats.mod2pi(2*Double32(pi)))
    @test iszero(DoubleFloats.modhalfpi(0.5*Double32(pi)))
    @test iszero(DoubleFloats.modqrtrpi(0.25*Double32(pi)))

    @test 0 < DoubleFloats.mod1pi(Double64(pi)*10) < eps(Double64(pi)*10)
    @test 0 < DoubleFloats.mod2pi(Double64(pi)*10) < eps(2*Double64(pi)*10)
    @test 0 < DoubleFloats.modhalfpi(Double64(pi)*10) < eps(0.5*Double64(pi)*10)
    @test 0 < DoubleFloats.modqrtrpi(0.25*Double64(pi)*10) < eps(0.25*Double64(pi)*10)

    @test 0 < DoubleFloats.mod1pi(Double32(pi)*10) < eps(Double32(pi)*10)
    @test 0 < DoubleFloats.mod2pi(Double32(pi)*10) < eps(2*Double32(pi)*10)
    @test 0 < DoubleFloats.modhalfpi(Double32(pi)*10) < eps(0.5*Double32(pi)*10)
    @test 0 < DoubleFloats.modqrtrpi(0.25*Double32(pi)*10) < eps(0.25*Double32(pi)*10)
end

@testset "rempi" begin
    @test iszero(DoubleFloats.rem1pi(Double64(pi)))
    @test iszero(DoubleFloats.rem2pi(2*Double64(pi)))
    @test iszero(DoubleFloats.remhalfpi(0.5*Double64(pi)))
    @test iszero(DoubleFloats.remqrtrpi(0.25*Double64(pi)))

    @test iszero(DoubleFloats.rem1pi(Double32(pi)))
    @test iszero(DoubleFloats.rem2pi(2*Double32(pi)))
    @test iszero(DoubleFloats.remhalfpi(0.5*Double32(pi)))
    @test iszero(DoubleFloats.remqrtrpi(0.25*Double32(pi)))
end

@testset "hypot" begin
    @test DoubleFloats.hypot(Double64(1), Double64(1))^2 == 2
    @test sum(DoubleFloats.normalize(Double64(1), Double64(1)))^2 == 2
end


text:
x = 100

println("x is ", x)

x*5

println("x is ", x)

x = x*5

println("x is ", x)

y = 200

x,y = y,x

println("x is ", x, " y is ", y)

अंशुल= 101

println("Value stored in symbols (Anshul) from language Hindi ", अंशुल)

# String
langname = "Julia"
println(typeof(langname))

println(typemax(Int32))
println(typemin(Int32))

x = Int16(10000)
println(x*x)

# largest possible value of 16-bit Int: 32767
x = typemax(Int16)
println(x)
# 32767

# Julia will maintain variable types even if the max value # is exceeding
# -32768
println(x + Int16(1))

# Float64
x=2.99e8
println(typeof(x))
# Replace e by f to generate Float32
y=2.99f8
println(typeof(y))

# true
println(x==y)

# This happens when that particular number doesn’t have the expected
# floating point representation.
x = 1.1; y = 0.1;
println(x+y)

# expressions
x = 4; y = 11;
println(3x + 4y + 93)
# 149

# this creates an empty array
empty_array = Float64[]
println(empty_array)
# 0-element Array{Float64,1}

# the function push!() adds values to the array
push!(empty_array,1.1)
println(empty_array)
# 1-element Array{Float64,1}:
# 1.1

# adding more values
push!(empty_array,2.2,3.3)
println(empty_array)
# 3-element Array{Float64,1}:
# 1.1
# 2.2
# 3.3

# constructing a matrix
X = [1 1 2; 3 5 8; 13 21 34]
println(X)
# 3×3 Array{Int64,2}:
# 1 1 2
# 3 5 8
# 13 21 34

# a 3x2 matrix
A = [2 4; 8 16; 32 64]
println(A)
# 3×2 Array{Int64,2}:
# 2 4
# 8 16
# 32 64

# We can change the shape of a matrix using the reshape() function:
println(reshape(A,2,3))
# 2×3 Array{Int64,2}:
# 2 32 16
# 8 4 64

# multidimensional array
multiA = rand(3,3,3)
println(multiA)
# 3×3×3 Array{Float64,3}:
# .........

# access the values just as we accessed values in arrays and matrices, using the index:
println(multiA[1,3,2])

# NA datavalues in dataset
using DataArrays
x = DataArray([1.1, 2.2, 3.3, 4.4, 5.5, 6.6])
println(x)

x[1] = NA
println(x)

# operations on DataArrays
println(mean(x[2:6]))

# Type aliases of Vector (a onedimensional
# Array type) and Matrix (a two-dimensional Array type) are DataVector
# and DataMatrix, provided by DataArray.
dvector = data([10,20,30,40,50])
println(dvector)
dmatrix = data([10 20 30; 40 50 60])
println(dmatrix)

# Pkg.add("DataFrames")
using DataFrames
df = DataFrame(Name = ["Julia", "Python"], Version = [0.5, 3.6])
println(df)

max_stars_repo_name: PacktPublishing/Learning-Julia

max_issues_repo_name: PacktPublishing/Learning-Julia

max_forks_repo_name: PacktPublishing/Learning-Julia


text:
@testset "Laplacian" begin
  points = P2[(0,0), (1,0), (0,1), (1,1), (0.5,0.5)]
  connec = connect.([(1,2,5),(2,4,5),(4,3,5),(3,1,5)], Triangle)
  mesh = SimpleMesh(points, connec)
  L = laplacematrix(mesh, weights=:uniform, normalize=true)
  @test L == [
     -1 1/3 1/3   0 1/3
    1/3  -1   0 1/3 1/3
    1/3   0  -1 1/3 1/3
      0 1/3 1/3  -1 1/3
    1/4 1/4 1/4 1/4  -1
  ]
end


text:
# # Global stress objective example
#-
#md # !!! tip
#md #     This example is also available as a Jupyter notebook:
#md #     [`global_stress.ipynb`](@__NBVIEWER_ROOT_URL__/examples/global_stress.ipynb)
#-
# ## Commented Program
#
# What follows is a program spliced with comments.
#md # The full program, without comments, can be found in the next [section](@ref global-stress-plain-program).

using TopOpt, LinearAlgebra

# ### Define the problem
E = 1.0 # Young’s modulus
v = 0.3 # Poisson’s ratio
f = 1.0 # downward force
rmin = 3.0 # filter radius

problems = Any[
    PointLoadCantilever(Val{:Linear}, (60, 20), (1.0, 1.0), E, v, f),
    HalfMBB(Val{:Linear}, (60, 20), (1.0, 1.0), E, v, f),
]
problem_names = [
    "Cantilever beam",
    "Half MBB beam",
    "L-beam",
    "Tie-beam",
]

i = 1
println(problem_names[i])
problem = problems[i]

# ### Parameter settings
V = 0.5 # volume fraction
xmin = 0.001 # minimum density
steps = 40 # maximum number of penalty steps, delta_p0 = 0.1
convcriteria = Nonconvex.KKTCriteria()
penalty = TopOpt.PowerPenalty(1.0)

# ### Define a finite element solver
solver = FEASolver(
    Displacement, Direct, problem, xmin = xmin, penalty = penalty,
)

# ### Define **stress** objective
# Notice that gradient is derived automatically by automatic differentiation (Zygote.jl)!
stress = TopOpt.MicroVonMisesStress(solver)
filter = if problem isa TopOptProblems.TieBeam
    identity
else
    DensityFilter(solver, rmin = rmin)
end
volfrac = TopOpt.Volume(problem, solver)

obj = x -> volfrac(filter(x))
constr = x -> norm(stress(filter(x)), 5) - 1.0
options = Nonconvex.MMAOptions(
    maxiter=2000, tol = Nonconvex.Tolerance(kkt = 1e-4),
)

# ### Define subproblem optimizer
x0 = fill(1.0, length(solver.vars))
optimizer = Optimizer(
    obj, constr, x0, Nonconvex.MMA87(),
    options = options, convcriteria = convcriteria,
)

# ### Define continuation SIMP optimizer
simp = SIMP(optimizer, solver, 3.0)

# ### Solve
result = simp(x0)

@show result.convstate
@show optimizer.workspace.iter
@show result.objval

# ### (Optional) Visualize the result using Makie.jl
# Need to run `using Pkg; Pkg.add(Makie)` first
# ```julia
# using TopOpt.TopOptProblems.Visualization: visualize
# fig = visualize(problem; topology = result.topology)
# Makie.display(fig)
# ```

#md # ## [Plain Program](@id global-stress-plain-program)
#md #
#md # Below follows a version of the program without any comments.
#md # The file is also available here: [global-stress.jl](global_stress.jl)
#md #
#md # ```julia
#md # @__CODE__
#md # ```