text:
using LinearAlgebra
import Base.copy

abstract type Algorithm end
abstract type GASolution end

#= Problem struct =#
"""
    Problem

Expects `eval_function` to have the form `f(s.x)` where `s` is a `Solution` and
to return a `Vector{Float64}`. If `has_constraint` is true, expects the final
value of the returned objective vector to be zero if solution is feasible, or
a positive Float64 indicating size of the constraint violation if infeasible.
"""
struct Problem
    nobjs::Int64
    maximize_objectives::Vector{Bool} # false -> minimize, true -> maximize
    nvars::Int64
    var_types::Vector{MOGA_Type}
    eval_fn::Function
    has_constraint::Bool
end

"no constraints constructor"
function Problem(nobjs::Int64, maximize_objectives::Vector{Bool},
                 nvars::Int64, var_types::Vector{T},
                 eval_fn::F) where {T <: MOGA_Type, F <: Function}
    Problem(nobjs, maximize_objectives, nvars, var_types,
            eval_fn, false)
end

#= Solutions =#
"""
    Solution

Type for candidate solutions, currently designed for NSGAII and NSGAIII.
"""
mutable struct Solution <: GASolution
    problem::Problem
    x::Vector{Float64}
    evaluated::Bool
    objectives::Vector{Float64}
    crowding_distance::Float64
    rank::Int64
    generation::Int64
    feasible::Bool
    constraint_violation::Float64
end

"no constraints constructor"
function Solution(problem::Problem,
                  x::Vector{Float64},
                  evaluated::Bool,
                  objectives::Vector{Float64},
                  crowding_distance::Float64,
                  rank::Int64;
                  generation::Int64=0)
    Solution(problem, x, evaluated, objectives, crowding_distance, rank,
             generation, true, 0.0)
end

"construct a candidate `Solutions` with a specific genome"
function Solution(problem::Problem,
                  x::Vector{Float64};
                  generation::Int64=0)
    Solution(problem, x, false, zeros(problem.nobjs), 0.0, 0,
        generation, true, 0)
end

"allocate a new solution that copies argument"
function copy(s::Solution)
    return Solution(
        s.problem,
        copy(s.x),
        s.evaluated,
        copy(s.objectives),
        s.crowding_distance,
        s.rank,
        s.generation,
        s.feasible,
        s.constraint_violation
    )
end

"Copies second argument solution into first argument"
function copyinto!(s_to::Solution, s_from::Solution)
    s_to.problem = s_from.problem
    s_to.x .= s_from.x
    s_to.evaluated = s_from.evaluated
    s_to.objectives .= s_from.objectives
    s_to.crowding_distance = s_from.crowding_distance
    s_to.rank = s_from.rank
    s_to.generation = s_from.generation
    s_to.feasible = s_from.feasible
    s_to.constraint_violation = s_from.constraint_violation
end

"""
    evaluate!(s::Solution)

Calls `s.problem.eval_fn`, saves into `s.objectives`, sets `s.evaluated` to
true. If `s.problem.has_constraint` is `true`, the last value in the objective
score array is treated as the magnitude of constraint violation, and the point
is considered feasible if the constraint violation is 0.
"""
function evaluate!(s::Solution)
    if !s.evaluated
        result = s.problem.eval_fn(s.x)
        if s.problem.has_constraint
            s.feasible = result[end] == 0.0
            s.constraint_violation = result[end]
            s.objectives = result[1:end-1]
        else
            s.objectives = result
        end
        s.evaluated = true
    end
    return s.objectives
end

# Archives
mutable struct Archive
    dominance::Function
    solutions::Vector{Solution}
end

function insert_solution!(archive::Archive, solution::Solution)
    flags = [archive.dominance(solution, arch_sol) for
             arch_sol in archive.solutions]
    dominates = [x > 0 for x in flags]
    nondominated = [x == 0 for x in flags]

    if any(dominates)
        return false
    else
        archive.solutions = archive.solutions[nondominated]
        push!(archive.solutions, solution)
        return true
    end
end

function insert_solutions!(archive::Archive, solutions::Vector{Solution})
    for s in solutions
        insert_solution!(archive, s)
    end
end

# Helper functions
"""
    compare_pareto_dominance(s1, s2)

returns -1 if `s1` dominates `s2`, 1 if `s2` dominates `s1`, and 0 otherwise
"""
function compare_pareto_dominance(s1, s2)
    prob = s1.problem
    # THIS IS WHERE TO INSERT CONSTRAINT CHECK

    dom1 = false
    dom2 = false

    for i in 1:prob.nobjs
        o1 = s1.objectives[i]
        o2 = s2.objectives[i]

        if prob.maximize_objectives[i] # == MAXIMIZE
            o1 = -1 * o1
            o2 = -1 * o2
        end

        if o1 < o2
            dom1 = true
            if dom2
                return 0
            end
        elseif o1 > o2
            dom2 = true
            if dom1
                return 0
            end
        end
    end

    if dom1 == dom2
        return 0
    elseif dom1
        return -1
    else
        return 1
    end

end

"""
    nondominated(solutions::Vector{Solution})

Returns all nondominated elements of `solutions`. This requires constructing
an `Archive` and inserting all solutions, so I'm not sure it's particularly fast.
"""
function nondominated(solutions::Vector{Solution})
    archive = Archive()
    insert_solutions!(archive, solutions)
    return archive.solutions
end

"""
    nondominated_sort(solutions::Vector{Solution})

Implements the nondominated sort for NSGAII. This function assigns front rank
and crowding distance scores to each solution in `solutions`. Note: better
function name might be `assign_rank_and_crowding!`, since it does not return
a sorted list.
"""
function nondominated_sort(solutions::Vector{Solution})
    rank::Int64 = 1

    for s in solutions
        s.rank = -1
    end

    while length(solutions) > 0
        archive = Archive(compare_pareto_dominance, Vector{Solution}())
        insert_solutions!(archive, solutions)

        for s in archive.solutions
            s.rank = rank
        end

        crowding_distance(archive.solutions)
        solutions = [s for s in solutions if s.rank == -1]
        rank += 1
    end
end

"""
    nondominated_cmp(x::Solution, y::Solution)

Returns -1 if x dominates y, 1 if y dominates x, and 0 otherwise. In the case
that both x and y are feasible, dominance is determined first by front rank,
and then by crowding distance. If only one solution is feasible, it automatically
dominates. If both are infeasbile, the solution with the smaller constraint
violation dominates.
"""
function nondominated_cmp(x::Solution, y::Solution)
    if x.feasible && y.feasible
        if x.rank == y.rank
            if -x.crowding_distance < -y.crowding_distance
                return -1
            elseif -x.crowding_distance > -y.crowding_distance
                return 1
            else
                return 0
            end
        else
            if x.rank < y.rank
                return -1
            elseif x.rank > y.rank
                return 1
            else
                return 0
            end
        end
    elseif x.feasible && !y.feasible
        return -1
    elseif !x.feasible && y.feasible
        return 1
    else
        if x.constraint_violation < y.constraint_violation
            return -1
        elseif y.constraint_violation < x.constraint_violation
            return 1
        else
            return 0
        end
    end
end

"""
    crowding_distance(solutions::Vector{Solution})

Assigns crowding distance, as defined in Deb 2002, to each element of
`solutions`.
"""
function crowding_distance(solutions::Vector{Solution})
    for s in solutions
        s.crowding_distance = 0.0
    end

    nobjs = length(solutions[1].objectives)

    if length(solutions) < 3
        for s in solutions
            s.crowding_distance = Inf
        end
    else
        # main case
        for i in 1:nobjs
            sorted_solutions = sort!(solutions, by=s->s.objectives[i])
            min_value = sorted_solutions[1].objectives[i]
            max_value = sorted_solutions[end].objectives[i]

            sorted_solutions[1].crowding_distance += Inf
            sorted_solutions[end].crowding_distance += Inf

            for j in 2:length(solutions)-1
                if max_value - min_value < eps(Float64)
                    sorted_solutions[j].crowding_distance = Inf
                else
                    diff = sorted_solutions[j+1].objectives[i] -
                           sorted_solutions[j-1].objectives[i]
                    sorted_solutions[j].crowding_distance += diff /
                                                (max_value - min_value)
                end
            end
        end
    end
end

"""
    nondominated_truncate(pop::Vector{Solution}, num::Int64,
                          dominance::Function=nondominated_cmp)

Implements nondominated truncate used in NSGAII. `num` gives the desired number
of individuals in the truncated population. The `dominance` function is used to
determine sort order, so dominance(x, y) should return -1 if x dominates, 1 if
y dominates, and 0 if neither dominates.
"""
function nondominated_truncate(pop::Vector{Solution}, num::Int64,
                               dominance::Function=nondominated_cmp)
    # maybe more efficient to not sort the whole thing but pull out by
    # rank until reach the marginal rank class. Not sure. This is easy.
    k(c1, c2) = dominance(c1, c2) < 0
    sorted = sort(pop, lt=k)    # TODO: sort in place (sort!)?
    return sorted[1:num]
end

"""
    reference_point_truncate(
        pop::Vector{Solution}, N::Int64, refpts::Array{Float64,2})

Implements reference point truncation algorithm used in NSGAIII. `N` is the
desired population size, and `refpts` is an array containing the reference
points.
"""
function reference_point_truncate(pop::Vector{Solution}, N::Int64,
              refpts::Array{Float64, 2})

    M = length(pop[1].objectives)

    # find cutoff dominance rank
    ranks = [p.rank for p in pop]
    rankcounts = zeros(Int64, maximum(ranks))
    for r in ranks
        rankcounts[r] += 1
    end
    marginalrank = findfirst(x -> x >= N, cumsum(rankcounts))
    @show marginalrank

    St = [i for i=1:length(pop) if pop[i].rank <= marginalrank]

    if length(St) == N
        return [copy(pop[i]) for i in St]
    end
    # else:
    Pnext = [i for i in St if pop[i].rank < marginalrank]
    marginalfront = [i for i in St if pop[i].rank == marginalrank]

    ## Find ideal point: minimal value for each objective across all solutions
    zmin = zeros(M)
    for j = 1:M
        zmin[j] = minimum([pop[i].objectives[j] for i in St])
    end
    @show zmin

    ## FIND EXTREME VALUES
    ## the goal of this operation is to find, for each axis j, the solution that
    ## minimizes the value of non-j values. This isn't exactly the same as
    ## maximizing j if the solution space has a weird shape.
    ## TODO: I think this is right, but I'm not clear about the max vs mins here...
    expts = [-1 for j = 1:M]
    for j = 1:M  # looping through objectives/axes
        w = 0.00001 * ones(M)
        w[j] = 1.0
        maxval = -Inf
        for i in St
            val = maximum((pop[i].objectives) ./ w)
            if val > maxval
                expts[j] = i
                maxval = val
            end
        end
    end

    ## FIND INTERCEPT VALUES
    ## ENTER: Have found one extreme point per objective. These define a hyper-plane.
    ## EXIT: Find the intersections of the hyper-plane with each objective axis.
    ##    These will be used to rescale the objective scores for these axes.
    ##
    b = ones(M)
    Z = transpose(reduce(hcat, [pop[i].objectives for i in expts]))
    # try Z = factorize(Z)
    if rank(Z) == M
        a = 1.0 ./ (Z\b)
    else
        maxvals = zeros(M)
        for j=1:M
            maxvals[j] = maximum(pop[i].objectives[j] for i in St)
        end
        a = 1.0 ./ maxvals
    end

    # fp = [p.objectives - zmin for p in St] # ::Array{Array{Float64, 1}, 1}
    # fn = [f ./ a for f in fp] # ::Array{Array{Float64, 1}, 1}

    @show length(Pnext)
    @show length(marginalfront)

    fp = Array{Float64,1}[(pop[i].objectives - zmin) ./ a for i in Pnext]
    fm = Array{Float64,1}[(pop[i].objectives - zmin) ./ a for i in marginalfront]

    @show typeof(fp)
    @show typeof(fm)

    ## ASSOCIATE SOLUTIONS WITH REFERENCE POINTS
    ## Do some linear algebra to find distances between the objective values
    ## and lines passing through each reference point.
    fp_ref_index, fp_ref_distance = associate_points(fp, refpts)
    fm_ref_index, fm_ref_distance = associate_points(fm, refpts)
    fm_index_dict = Dict(i => d for (i, d) in zip(marginalfront, fm_ref_index))
    fm_dist_dict = Dict(i => d for (i, d) in zip(marginalfront, fm_ref_distance))

    ## NICHING PROCEDURE
    ## select points from the marginal front to include in P_{t+1}
    ## Prefer points associated with reference points with fewest associated
    ## solutions.

    K = N - length(Pnext)
    k = 1
    refpt_deck = Set(1:length(refpts))
    candidate_deck = Set(marginalfront) # set of indices into pop

    ρ = Dict(rp => 0 for rp in refpt_deck)
    for i in fp_ref_index
        ρ[i] += 1
    end

    while k <= K
        # find all reference points with minimal # of associated solutions from Pnext

        ρmin = minimum(values(ρ))
        # @show(ρmin)
        Jmin = [r for r in refpt_deck if ρ[r] == ρmin]
        # @show(Jmin)  # set of ref pts with minimal representation
        jbar = rand(Jmin)
        # find the set of all solutions in the marginal front associated with this ref point
        Ij = [i for i in candidate_deck if fm_index_dict[i] == jbar]
        if length(Ij) > 0
            # pick a point
            if ρ[jbar] == 0
                # pick the closest point
                distances = [fm_dist_dict[i] for i in Ij]
                chosen = Ij[argmin(distances)]
            else
                # choose any of the points
                chosen = rand(Ij)
            end
            push!(Pnext, chosen)
            ρ[jbar] += 1
            pop!(candidate_deck, chosen)
            k += 1
        else
            # case with no solutions in marginal front associated with this
            # reference point. Remove it from the pool for this iteration.
            pop!(refpt_deck, jbar)
            pop!(ρ, jbar)
        end
    end

    newpop = [copy(pop[i]) for i in Pnext] ## TODO: copies necessary?
    return newpop
end

max_stars_repo_name: phawthorne/echidna

max_issues_repo_name: phawthorne/echidna

max_forks_repo_name: phawthorne/echidna


text:
@testset "Gradients" begin
    inter = LennardJones(force_units=NoUnits, energy_units=NoUnits)
    box_size = SVector(5.0, 5.0, 5.0)
    a1, a2 = Atom(σ=0.3, ϵ=0.5), Atom(σ=0.3, ϵ=0.5)

    function force_direct(dist)
        c1 = SVector(1.0, 1.0, 1.0)
        c2 = SVector(dist + 1.0, 1.0, 1.0)
        vec = vector(c1, c2, box_size)
        F = force(inter, vec, c1, c2, a1, a2, box_size)
        return F[1]
    end

    function force_grad(dist)
        grad = gradient(dist) do dist
            c1 = SVector(1.0, 1.0, 1.0)
            c2 = SVector(dist + 1.0, 1.0, 1.0)
            vec = vector(c1, c2, box_size)
            potential_energy(inter, vec, c1, c2, a1, a2, box_size)
        end
        return -grad[1]
    end

    dists = collect(0.2:0.01:1.2)
    forces_direct = force_direct.(dists)
    forces_grad = force_grad.(dists)
    @test all(isapprox.(forces_direct, forces_grad))

    abs2_vec(x) = abs2.(x)

    # Function is strange in order to work with gradients on the GPU
    function mean_min_separation(coords, box_size)
        n_atoms = length(coords)
        coords_rep = repeat(reshape(coords, n_atoms, 1), 1, n_atoms)
        vec2arg(c1, c2) = vector(c1, c2, box_size)
        diffs = vec2arg.(coords_rep, permutedims(coords_rep, (2, 1)))
        disps = Array(sum.(abs2_vec.(diffs)))
        disps_diag = disps .+ Diagonal(100 * ones(typeof(box_size[1]), n_atoms))
        return mean(sqrt.(minimum(disps_diag; dims=1)))
    end

    function test_grad(gpu::Bool, forward::Bool, f32::Bool, pis::Bool,
                        sis::Bool, implicit_solvent::Bool)
        n_atoms = 50
        n_steps = 100
        atom_mass = f32 ? 10.0f0 : 10.0
        box_size = f32 ? SVector(3.0f0, 3.0f0, 3.0f0) : SVector(3.0, 3.0, 3.0)
        temp = f32 ? 1.0f0 : 1.0
        simulator = VelocityVerlet(
            dt=f32 ? 0.001f0 : 0.001,
            coupling=RescaleThermostat(temp),
        )
        coords = place_atoms(n_atoms, box_size, f32 ? 0.6f0 : 0.6)
        velocities = [velocity(atom_mass, temp) for i in 1:n_atoms]
        coords_dual = [ForwardDiff.Dual.(x, f32 ? 0.0f0 : 0.0) for x in coords]
        velocities_dual = [ForwardDiff.Dual.(x, f32 ? 0.0f0 : 0.0) for x in velocities]
        nb_cutoff = f32 ? 1.2f0 : 1.2
        lj = LennardJones(
            cutoff=DistanceCutoff(nb_cutoff),
            nl_only=true,
            force_units=NoUnits,
            energy_units=NoUnits,
        )
        crf = CoulombReactionField(
            dist_cutoff=nb_cutoff,
            solvent_dielectric=f32 ? Float32(Molly.solventdielectric) : Molly.solventdielectric,
            nl_only=true,
            coulomb_const=f32 ? Float32(ustrip(Molly.coulombconst)) : ustrip(Molly.coulombconst),
            force_units=NoUnits,
            energy_units=NoUnits,
        )
        pairwise_inters = pis ? (lj, crf) : ()
        bond_is, bond_js = collect(1:(n_atoms ÷ 2)), collect((1 + n_atoms ÷ 2):n_atoms)
        bond_dists = [norm(vector(Array(coords)[i], Array(coords)[i + n_atoms ÷ 2], box_size)) for i in 1:(n_atoms ÷ 2)]
        angles_inner = [HarmonicAngle(th0=f32 ? 2.0f0 : 2.0, cth=f32 ? 10.0f0 : 10.0) for i in 1:15]
        angles = InteractionList3Atoms(
            collect(1:15),
            collect(16:30),
            collect(31:45),
            repeat([""], 15),
            gpu ? cu(angles_inner) : angles_inner,
        )
        torsions_inner = [PeriodicTorsion(
                periodicities=[1, 2, 3],
                phases=f32 ? [1.0f0, 0.0f0, -1.0f0] : [1.0, 0.0, -1.0],
                ks=f32 ? [10.0f0, 5.0f0, 8.0f0] : [10.0, 5.0, 8.0],
                n_terms=6,
            ) for i in 1:10]
        torsions = InteractionList4Atoms(
            collect(1:10),
            collect(11:20),
            collect(21:30),
            collect(31:40),
            repeat([""], 10),
            gpu ? cu(torsions_inner) : torsions_inner,
        )
        atoms_setup = [Atom(charge=f32 ? 0.0f0 : 0.0, σ=f32 ? 0.0f0 : 0.0) for i in 1:n_atoms]
        imp = ImplicitSolventOBC(
            gpu ? cu(atoms_setup) : atoms_setup,
            [AtomData(element="O") for i in 1:n_atoms],
            InteractionList2Atoms(bond_is, bond_js, [""], nothing);
            use_OBC2=true,
        )
        general_inters = implicit_solvent ? (imp,) : ()
        neighbor_finder = DistanceVecNeighborFinder(
            nb_matrix=gpu ? cu(trues(n_atoms, n_atoms)) : trues(n_atoms, n_atoms),
            n_steps=10,
            dist_cutoff=f32 ? 1.5f0 : 1.5,
        )

        function loss(σ, kb)
            if f32
                atoms = [Atom(i, i % 2 == 0 ? -0.02f0 : 0.02f0, atom_mass, σ, 0.2f0, false) for i in 1:n_atoms]
            else
                atoms = [Atom(i, i % 2 == 0 ? -0.02 : 0.02, atom_mass, σ, 0.2, false) for i in 1:n_atoms]
            end

            bonds_inner = [HarmonicBond(bond_dists[i], kb) for i in 1:(n_atoms ÷ 2)]
            bonds = InteractionList2Atoms(
                bond_is,
                bond_js,
                repeat([""], length(bonds_inner)),
                gpu ? cu(bonds_inner) : bonds_inner,
            )
            cs = deepcopy(forward ? coords_dual : coords)
            vs = deepcopy(forward ? velocities_dual : velocities)

            s = System(
                atoms=gpu ? cu(atoms) : atoms,
                pairwise_inters=pairwise_inters,
                specific_inter_lists=sis ? (bonds, angles, torsions) : (),
                general_inters=general_inters,
                coords=gpu ? cu(cs) : cs,
                velocities=gpu ? cu(vs) : vs,
                box_size=box_size,
                neighbor_finder=neighbor_finder,
                gpu_diff_safe=true,
                force_units=NoUnits,
                energy_units=NoUnits,
            )

            simulate!(s, simulator, n_steps)

            return mean_min_separation(s.coords, box_size)
        end

        return loss
    end

    runs = [
        ("cpu"                 , [false, false, false, true , true , false], 0.1 , 0.25),
        ("cpu forward"         , [false, true , false, true , true , false], 0.01, 0.01),
        ("cpu f32"             , [false, false, true , true , true , false], 0.2 , 10.0),
        ("cpu nospecific"      , [false, false, false, true , false, false], 0.1 , 0.0 ),
        ("cpu nopairwise"      , [false, false, false, false, true , false], 0.0 , 0.25),
        ("cpu implicit"        , [false, false, false, true , true , true ], 0.1 , 0.25),
        ("cpu implicit forward", [false, true , false, true , true , true ], 0.01, 0.01),
    ]
    if run_gpu_tests
        push!(runs, ("gpu"                 , [true , false, false, true , true , false], 0.25, 20.0))
        push!(runs, ("gpu forward"         , [true , true , false, true , true , false], 0.01, 0.01))
        push!(runs, ("gpu f32"             , [true , false, true , true , true , false], 0.5 , 50.0))
        push!(runs, ("gpu nospecific"      , [true , false, false, true , false, false], 0.25, 0.0 ))
        push!(runs, ("gpu nopairwise"      , [true , false, false, false, true , false], 0.0 , 10.0))
        push!(runs, ("gpu implicit"        , [true , false, false, true , true , true ], 0.25, 20.0))
        push!(runs, ("gpu implicit forward", [true , true , false, true , true , true ], 0.01, 0.01))
    end

    for (name, args, tol_σ, tol_kb) in runs
        forward, f32 = args[2], args[3]
        σ = f32 ? 0.4f0 : 0.4
        kb = f32 ? 100.0f0 : 100.0
        f = test_grad(args...)
        if forward
            # Run once to setup
            grad_zygote = (
                gradient((σ, kb) -> Zygote.forwarddiff(σ  -> f(σ, kb), σ ), σ, kb)[1],
                gradient((σ, kb) -> Zygote.forwarddiff(kb -> f(σ, kb), kb), σ, kb)[2],
            )
            grad_zygote = (
                gradient((σ, kb) -> Zygote.forwarddiff(σ  -> f(σ, kb), σ ), σ, kb)[1],
                gradient((σ, kb) -> Zygote.forwarddiff(kb -> f(σ, kb), kb), σ, kb)[2],
            )
        else
            # Run once to setup
            grad_zygote = gradient(f, σ, kb)
            grad_zygote = gradient(f, σ, kb)
        end
        grad_fd = (
            central_fdm(6, 1)(σ  -> ForwardDiff.value(f(σ, kb)), σ ),
            central_fdm(6, 1)(kb -> ForwardDiff.value(f(σ, kb)), kb),
        )
        for (prefix, gzy, gfd, tol) in zip(("σ", "kb"), grad_zygote, grad_fd, (tol_σ, tol_kb))
            if abs(gfd) < 1e-13
                @test isnothing(gzy) || abs(gzy) < 1e-12
            else
                frac_diff = abs(gzy - gfd) / abs(gfd)
                @info "$(rpad(name, 20)) - $(rpad(prefix, 2)) - FD $gfd, Zygote $gzy, fractional difference $frac_diff"
                @test frac_diff < tol
            end
        end
    end
end

ext: jl


text:
function link(j::Int64,
              q::Int64,
              k::Int64,
              g::Int64,
              L::Array{Float64,2},
              X::Array{Int64,2})
    out = 0.0
    if q != j && X[q,g] == 1
        out = L[j,q] * X[j,k] * X[q,g]
    end
    return out::Float64
end

function INL(j::Int64,
            w1::Int64,
            L::Array{Float64,2},
            X::Array{Int64,2})
    out = 0.0
    for q = 1:size(L,1)
        out += link(j::Int64,
                    q::Int64,
                    w1::Int64,
                    w1::Int64,
                    L::Array{Float64,2},
                    X::Array{Int64,2})
    end
    return out::Float64
end

function OUTL(j::Int64,
              w1::Int64,
              w2::Int64,
              L::Array{Float64,2},
              X::Array{Int64,2})
    out = 0.0
    for q = 1:size(L,1)
        out += link(j::Int64,
                    q::Int64,
                    w1::Int64,
                    w2::Int64,
                    L::Array{Float64,2},
                    X::Array{Int64,2})
    end
    return out::Float64
end

function DL(j::Int64,
            w1::Int64,
            w2::Int64,
            L::Array{Float64,2},
            X::Array{Int64,2})
    OUTL(j,w1,w2,L,X) - INL(j,w1,L,X)
end

function PAYOFF(j::Int64,
                q::Int64,
                w1::Int64,
                w2::Int64,
                L::Array{Float64,2},
                X::Array{Int64,2})
    DL(j,w1,w2,L,X) + DL(q,w2,w1,L,X) - 2*L[j,q]
end

function LINKS(trans::Array{Int64,2},
               ov::Array{Float64,1})
    L = Array{Float64,2}(undef,size(trans,2),size(trans,2)) .= 0
    for j = 1:size(trans,2)
        for q = 1:size(trans,2)
            if j == q
                L[j,q] = 0
            else
                for i = 1:size(trans,1)
                    L[j,q] += ov[i] * trans[i,j] * trans[i,q]
                end
            end
        end
    end
    return L::Array{Float64,2}
end

function LINKADJUST(trans::Array{Int64,2})
    ov = Array{Float64,1}(undef,size(trans,1))
    ov .= 0
    for i = 1:size(trans,1)
        if sum(trans[i,:]) > 1
            ov[i] = sum(trans[i,:])-1/(factorial(big(sum(trans[i,:])))/factorial(2)*factorial(big(sum(trans[i,:])-2)))
        end
    end
    return ov::Array{Float64,1}
end

function LWbest(L::Array{Float64,2})
    LW_best = 0
    for j = 1:size(L,2)
        for q = j+1:size(L,2)
            if j < q
                LW_best += L[j,q]
            end
        end
    end
    return LW_best::Int64
end

function LW(L::Array{Float64,2},
            X::Array{Int64,2})
    out = 0
    for k = 1:(size(X,2))
        for j = 1:size(X,1)
            for q = j+1:size(X,1)
                if j < q
                    out += L[j,q] * X[j,k] * X[q,k]
                end
            end
        end
    end
    return out::Float64
end

function STRATEGY1(X::Array{Int64,2},
                   m::Int64,
                   L::Array{Float64,2},
                   capacity::Array{Int64,1},
                   stop::Int64)
    dl_arr = Array{Float64,3}(undef,size(X,1),size(X,1),size(X,2)) .= 0
    for i in 1:size(X,1)
        if X[i,m] == 1
            for g = 1:size(X,2)
                if sum(X[:,g]) < capacity[g] && m!=g
                    dl_arr[i,m,g] = DL(i,m,g,L,X)
                end
            end
        end
    end
    if findmax(dl_arr)[1] > 0
        X[getindex(findmax(dl_arr)[2],1),getindex(findmax(dl_arr)[2],2)] = 0
        X[getindex(findmax(dl_arr)[2],1),getindex(findmax(dl_arr)[2],3)] = 1
        stop = 0
    end
    return X::Array{Int64,2}
end

function STRATEGY2(X::Array{Int64,2},
                   m::Int64,
                   L::Array{Float64,2},
                   capacity::Array{Int64,1},
                   stop::Int64)
    pay_arr = Array{Float64,4}(undef,size(X,1),size(X,1),size(capacity,1),size(capacity,1)) .= 0
    for i in 1:size(X,1)
        if X[i,m] == 1
            for j in 1:size(X,1)
                for g = 1:size(X,2)
                    if m != g && X[j,g] == 1
                        pay_arr[i,j,m,g] = PAYOFF(i,j,m,g,L,X)
                    end
                end
            end
        end
    end
    if findmax(pay_arr)[1] > 0
        X[getindex(findmax(pay_arr)[2],1),getindex(findmax(pay_arr)[2],3)] = 0
        X[getindex(findmax(pay_arr)[2],1),getindex(findmax(pay_arr)[2],4)] = 1
        X[getindex(findmax(pay_arr)[2],2),getindex(findmax(pay_arr)[2],4)] = 0
        X[getindex(findmax(pay_arr)[2],2),getindex(findmax(pay_arr)[2],3)] = 1
        stop = 0
    end
    return X::Array{Int64,2}
end

max_stars_repo_name: beyondatlas/splitdeliveries

max_issues_repo_name: beyondatlas/splitdeliveries

max_forks_repo_name: beyondatlas/splitdeliveries


text:
using Base: Float64
using AlphaZero
using ProgressLogging

include("../mu_game_wrapper.jl")
include("../network.jl")
include("../trace.jl")
include("../play.jl")
#include("memory.jl")
include("../training.jl")

# experiment = Examples.experiments["tictactoe"]
# session = Session(experiment)
# env = session.env

# oracle = Network.copy(env.bestnn, on_gpu=env.params.self_play.sim.use_gpu, test_mode=true)
# player = MctsPlayer(env.gspec, oracle, env.params.self_play.mcts)

# game = GI.init(env.gspec)

# function AlphaZero.think(p::MctsPlayer, game)
#     hidden_state = Representation()(GI.current_state(game))
#     mugame = MuGameEnvWrapper(game,
#         CachedOracle(Dynamics(gspec=env.gspec),Dict{Tuple{typeof(hidden_state), Int},Tuple{Float64, typeof(hidden_state)}}()),
#         Representation(),
#         hidden_state,
#         true,
#         GI.white_playing(game),
#         0.)

#     if isnothing(p.timeout) # Fixed number of MCTS simulations
#         MCTS.explore!(p.mcts, mugame, p.niters)
#     else # Run simulations until timeout
#         start = time()
#         while time() - start < p.timeout
#         MCTS.explore!(p.mcts, mugame, p.niters)
#         end
#     end
#     return MCTS.policy(p.mcts, mugame)
# end

# actions, π_target = think(player, game) #overload this function, perhaps MuMCTSplayer

gspec = Examples.games["tictactoe"]

representationhp = RepresentationHP(width=200, depth=4, hiddenstate_shape=32)
representation_oracle = RepresentationNetwork(gspec, representationhp)

dynamicshp = DynamicsHP(hiddenstate_shape=32, width=200, depth_common=4)
dynamics_oracle = DynamicsNetwork(gspec, dynamicshp)

predictionhp = PredictionHP(hiddenstate_shape=32, width=200, depth_common=4)
prediction_oracle = PredictionNetwork(gspec, predictionhp)

nns = (f=prediction_oracle, g=dynamics_oracle, h=representation_oracle)

import LinearAlgebra
LinearAlgebra.BLAS.set_num_threads(1)

n=10
sim=SimParams(
    num_games=10,
    num_workers=n,
    batch_size=n,
    use_gpu=false,
    reset_every=4,
    flip_probability=0.,
    alternate_colors=false)

  mcts = MctsParams(
    num_iters_per_turn=400,
    cpuct=1.0,
    temperature=ConstSchedule(1.0),
    dirichlet_noise_ϵ=0.2,
    dirichlet_noise_α=1.0)

  self_play = (; sim, mcts)

  #memory = CircularBuffer{Trace{GI.state_type(gspec)}}(1024)

  μparams = MuParams(self_play, 1, 100)

  env = MuEnv(gspec, μparams, nns)

  t = @timed self_play_step!(env)
  # push!(res, (n, t))

n_workers = 1:16
res2 = []
for n in n_workers
  sim=SimParams(
    num_games=100,
    num_workers=n,
    batch_size=n,
    use_gpu=false,
    reset_every=4,
    flip_probability=0.,
    alternate_colors=false)

  mcts = MctsParams(
    num_iters_per_turn=400,
    cpuct=1.0,
    temperature=ConstSchedule(1.0),
    dirichlet_noise_ϵ=0.2,
    dirichlet_noise_α=1.0)

  self_play = (; sim, mcts)

  #memory = CircularBuffer{Trace{GI.state_type(gspec)}}(1024)

  μparams = MuParams(self_play, 1, 100)

  env = MuEnv(gspec, μparams, nns)

  t = @timed self_play_step!(env)
  push!(res2, (n,t))
  @info res2

end

using Plots
# bar(first.(res), last.(res), xlabel="num_workers", ylabel="time [s]", label="num_games=100")
bar(first.(res2), last.(res2), xlabel="num_workers", ylabel="time [s]", label="num_games=100")


text:
"""
    stress(Δ, Λ, W = ones(size(Δ)), iterable=eachindex(Δ))

Compute the stress: Σ wij (Δij -  Λij)^2
"""
function stress(Δ, Λ, W = ones(size(Δ)))
    n = size(Δ, 1)
    s = 0.0
    for idx in eachindex(Δ)
        s += W[idx] * (Δ[idx] - Λ[idx])^2
    end
    s
end

function log_output(k, σ, verbose)
    verbose == false && return
    if k == 1
        println("Iteration\tStress")
        println("---------\t------")
    end
    @printf("%3d\t      %-6f\n", k, round(σ[k], digits = 6))
end


text:
using TimeZones: Transition
using Dates: Hour, Second, UTM, @dateformat_str

@testset "ZonedDateTime" begin
    utc = FixedTimeZone("UTC", 0, 0)
    warsaw = first(compile("Europe/Warsaw", tzdata["europe"]))

    @testset "dateformat parsing" begin
        @testset "successful parsing: $f" for f in (parse, tryparse)
            # Make sure all dateformat codes parse correctly
            # yYmuUdHMSseEzZ and yyyymmdd
            zdt = ZonedDateTime(1, 2, 3, 4, 5, 6, 7, utc)
            # Test y, u, d, H, M, S, s, Z
            p_zdt = f(
                ZonedDateTime,
                "Feb 3 1, 4:5:6.007 UTC",
                dateformat"u d y, H:M:S.s Z",
            )
            @test zdt == p_zdt

            # Test m, e, Y, z
            p_zdt = f(
                ZonedDateTime,
                "2 mon 3 1, 4:5:6.007+00:00",
                dateformat"m e d Y, H:M:S.s+z",
            )
            @test zdt == p_zdt

            # Test E, U
            p_zdt = f(
                ZonedDateTime,
                "February Monday 3 1 4:5:6.007 UTC",
                dateformat"U E d y H:M:S.s Z",
            )
            @test zdt == p_zdt

            # Test yyyymmdd
            p_zdt = f(
                ZonedDateTime,
                "00010203 4:5:6.007 UTC",
                dateformat"yyyymmdd H:M:S.s Z",
            )
            @test zdt == p_zdt
        end

        @testset "failed parsing: parse" begin
            @test_throws ArgumentError parse(
                ZonedDateTime,
                "2015-07-29 11:12:13.456 FakeTZ",
                dateformat"yyyy-mm-dd HH:MM:SS.sss Z",
            )

            @test_throws ArgumentError parse(
                ZonedDateTime,
                "2015-07-29 11:12:13.456",
                dateformat"yyyy-mm-dd HH:MM:SS.sss Z",
            )
        end

        @testset "failed parsing: tryparse" begin
            @test tryparse(
                ZonedDateTime,
                "2015-07-29 11:12:13.456 FakeTZ",
                dateformat"yyyy-mm-dd HH:MM:SS.sss Z",
            ) === nothing

            @test tryparse(
                ZonedDateTime,
                "2015-07-29 11:12:13.456",
                dateformat"yyyy-mm-dd HH:MM:SS.sss Z",
            ) === nothing
        end
    end

    @testset "standard time" begin
        local_dt = DateTime(1916, 2, 1, 0)
        utc_dt = DateTime(1916, 1, 31, 23)

        # Disambiguating parameters ignored when there is no ambiguity.
        @test ZonedDateTime(local_dt, warsaw).zone.name == "CET"
        @test ZonedDateTime(local_dt, warsaw, 0).zone.name == "CET"
        @test ZonedDateTime(local_dt, warsaw, 1).zone.name == "CET"
        @test ZonedDateTime(local_dt, warsaw, 2).zone.name == "CET"
        @test ZonedDateTime(local_dt, warsaw, true).zone.name == "CET"
        @test ZonedDateTime(local_dt, warsaw, false).zone.name == "CET"
        @test ZonedDateTime(utc_dt, warsaw, from_utc=true).zone.name == "CET"

        @test ZonedDateTime(local_dt, warsaw).utc_datetime == utc_dt
        @test ZonedDateTime(local_dt, warsaw, 0).utc_datetime == utc_dt
        @test ZonedDateTime(local_dt, warsaw, 1).utc_datetime == utc_dt
        @test ZonedDateTime(local_dt, warsaw, 2).utc_datetime == utc_dt
        @test ZonedDateTime(local_dt, warsaw, true).utc_datetime == utc_dt
        @test ZonedDateTime(local_dt, warsaw, false).utc_datetime == utc_dt
        @test ZonedDateTime(utc_dt, warsaw, from_utc=true).utc_datetime == utc_dt
    end

    @testset "daylight saving time" begin
        local_dt = DateTime(1916, 6, 1, 0)
        utc_dt = DateTime(1916, 5, 31, 22)

        # Disambiguating parameters ignored when there is no ambiguity.
        @test ZonedDateTime(local_dt, warsaw).zone.name == "CEST"
        @test ZonedDateTime(local_dt, warsaw, 0).zone.name == "CEST"
        @test ZonedDateTime(local_dt, warsaw, 1).zone.name == "CEST"
        @test ZonedDateTime(local_dt, warsaw, 2).zone.name == "CEST"
        @test ZonedDateTime(local_dt, warsaw, true).zone.name == "CEST"
        @test ZonedDateTime(local_dt, warsaw, false).zone.name == "CEST"
        @test ZonedDateTime(utc_dt, warsaw, from_utc=true).zone.name == "CEST"

        @test ZonedDateTime(local_dt, warsaw).utc_datetime == utc_dt
        @test ZonedDateTime(local_dt, warsaw, 0).utc_datetime == utc_dt
        @test ZonedDateTime(local_dt, warsaw, 1).utc_datetime == utc_dt
        @test ZonedDateTime(local_dt, warsaw, 2).utc_datetime == utc_dt
        @test ZonedDateTime(local_dt, warsaw, true).utc_datetime == utc_dt
        @test ZonedDateTime(local_dt, warsaw, false).utc_datetime == utc_dt
        @test ZonedDateTime(utc_dt, warsaw, from_utc=true).utc_datetime == utc_dt
    end

    @testset "spring-forward" begin
        local_dts = (
            DateTime(1916,4,30,22),
            DateTime(1916,4,30,23),
            DateTime(1916,5,1,0),
        )
        utc_dts = (
            DateTime(1916,4,30,21),
            DateTime(1916,4,30,22),
        )
        @test_throws NonExistentTimeError ZonedDateTime(local_dts[2], warsaw)
        @test_throws NonExistentTimeError ZonedDateTime(local_dts[2], warsaw, 0)
        @test_throws NonExistentTimeError ZonedDateTime(local_dts[2], warsaw, 1)
        @test_throws NonExistentTimeError ZonedDateTime(local_dts[2], warsaw, 2)
        @test_throws NonExistentTimeError ZonedDateTime(local_dts[2], warsaw, true)
        @test_throws NonExistentTimeError ZonedDateTime(local_dts[2], warsaw, false)

        @test ZonedDateTime(local_dts[1], warsaw).zone.name == "CET"
        @test ZonedDateTime(local_dts[3], warsaw).zone.name == "CEST"
        @test ZonedDateTime(utc_dts[1], warsaw, from_utc=true).zone.name == "CET"
        @test ZonedDateTime(utc_dts[2], warsaw, from_utc=true).zone.name == "CEST"

        @test ZonedDateTime(local_dts[1], warsaw).utc_datetime == utc_dts[1]
        @test ZonedDateTime(local_dts[3], warsaw).utc_datetime == utc_dts[2]
        @test ZonedDateTime(utc_dts[1], warsaw, from_utc=true).utc_datetime == utc_dts[1]
        @test ZonedDateTime(utc_dts[2], warsaw, from_utc=true).utc_datetime == utc_dts[2]
    end

    @testset "fall-back" begin
        local_dt = DateTime(1916, 10, 1, 0)
        utc_dts = (DateTime(1916, 9, 30, 22), DateTime(1916, 9, 30, 23))

        @test_throws AmbiguousTimeError ZonedDateTime(local_dt, warsaw)
        @test_throws AmbiguousTimeError ZonedDateTime(local_dt, warsaw, 0)

        @test ZonedDateTime(local_dt, warsaw, 1).zone.name == "CEST"
        @test ZonedDateTime(local_dt, warsaw, 2).zone.name == "CET"
        @test ZonedDateTime(local_dt, warsaw, true).zone.name == "CEST"
        @test ZonedDateTime(local_dt, warsaw, false).zone.name == "CET"
        @test ZonedDateTime(utc_dts[1], warsaw, from_utc=true).zone.name == "CEST"
        @test ZonedDateTime(utc_dts[2], warsaw, from_utc=true).zone.name == "CET"

        @test ZonedDateTime(local_dt, warsaw, 1).utc_datetime == utc_dts[1]
        @test ZonedDateTime(local_dt, warsaw, 2).utc_datetime == utc_dts[2]
        @test ZonedDateTime(local_dt, warsaw, true).utc_datetime == utc_dts[1]
        @test ZonedDateTime(local_dt, warsaw, false).utc_datetime == utc_dts[2]
        @test ZonedDateTime(utc_dts[1], warsaw, from_utc=true).utc_datetime == utc_dts[1]
        @test ZonedDateTime(utc_dts[2], warsaw, from_utc=true).utc_datetime == utc_dts[2]
    end

    @testset "standard offset reduced" begin
        # The zone's standard offset was changed from +2 to +1 creating an ambiguous hour
        local_dt = DateTime(1922,5,31,23)
        utc_dts = (DateTime(1922, 5, 31, 21), DateTime(1922, 5, 31, 22))
        @test_throws AmbiguousTimeError ZonedDateTime(local_dt, warsaw)

        @test ZonedDateTime(local_dt, warsaw, 1).zone.name == "EET"
        @test ZonedDateTime(local_dt, warsaw, 2).zone.name == "CET"
        @test_throws AmbiguousTimeError ZonedDateTime(local_dt, warsaw, true)
        @test_throws AmbiguousTimeError ZonedDateTime(local_dt, warsaw, false)
        @test ZonedDateTime(utc_dts[1], warsaw, from_utc=true).zone.name == "EET"
        @test ZonedDateTime(utc_dts[2], warsaw, from_utc=true).zone.name == "CET"

        @test ZonedDateTime(local_dt, warsaw, 1).utc_datetime == utc_dts[1]
        @test ZonedDateTime(local_dt, warsaw, 2).utc_datetime == utc_dts[2]
        @test ZonedDateTime(utc_dts[1], warsaw, from_utc=true).utc_datetime == utc_dts[1]
        @test ZonedDateTime(utc_dts[2], warsaw, from_utc=true).utc_datetime == utc_dts[2]
    end

    @testset "2-hour daylight saving time offset" begin
        # Check behaviour when the "save" offset is larger than an hour.
        paris = first(compile("Europe/Paris", tzdata["europe"]))

        @test ZonedDateTime(DateTime(1945,4,2,1), paris).zone == FixedTimeZone("WEST", 0, 3600)
        @test_throws NonExistentTimeError ZonedDateTime(DateTime(1945,4,2,2), paris)
        @test ZonedDateTime(DateTime(1945,4,2,3), paris).zone == FixedTimeZone("WEMT", 0, 7200)

        @test_throws AmbiguousTimeError ZonedDateTime(DateTime(1945,9,16,2), paris)
        @test ZonedDateTime(DateTime(1945,9,16,2), paris, 1).zone == FixedTimeZone("WEMT", 0, 7200)
        @test ZonedDateTime(DateTime(1945,9,16,2), paris, 2).zone == FixedTimeZone("CET", 3600, 0)

        # Ensure that dates are continuous when both a UTC offset and the DST offset change.
        @test ZonedDateTime(DateTime(1945,9,16,1), paris).utc_datetime == DateTime(1945,9,15,23)
        @test ZonedDateTime(DateTime(1945,9,16,2), paris, 1).utc_datetime == DateTime(1945,9,16,0)
        @test ZonedDateTime(DateTime(1945,9,16,2), paris, 2).utc_datetime == DateTime(1945,9,16,1)
        @test ZonedDateTime(DateTime(1945,9,16,3), paris).utc_datetime == DateTime(1945,9,16,2)
    end

    @testset "multiple hour transitions" begin
        # Transitions changes that exceed an hour. Results in having two sequential
        # non-existent hour and two sequential ambiguous hours.
        t = VariableTimeZone("Testing", [
            Transition(DateTime(1800,1,1), FixedTimeZone("TST",0,0)),
            Transition(DateTime(1950,4,1), FixedTimeZone("TDT",0,7200)),
            Transition(DateTime(1950,9,1), FixedTimeZone("TST",0,0)),
        ])

        # A "spring forward" where 2 hours are skipped.
        @test ZonedDateTime(DateTime(1950,3,31,23), t).zone == FixedTimeZone("TST",0,0)
        @test_throws NonExistentTimeError ZonedDateTime(DateTime(1950,4,1,0), t)
        @test_throws NonExistentTimeError ZonedDateTime(DateTime(1950,4,1,1), t)
        @test ZonedDateTime(DateTime(1950,4,1,2), t).zone == FixedTimeZone("TDT",0,7200)

        # A "fall back" where 2 hours are duplicated. Never appears to occur in reality.
        @test ZonedDateTime(DateTime(1950,8,31,23), t).utc_datetime == DateTime(1950,8,31,21)  # TDT

        # First occurrences of duplicated hours.
        @test ZonedDateTime(DateTime(1950,9,1,0), t, 1).utc_datetime == DateTime(1950,8,31,22) # TST
        @test ZonedDateTime(DateTime(1950,9,1,1), t, 1).utc_datetime == DateTime(1950,8,31,23) # TST

        # Second occurrences of duplicated hours.
        @test ZonedDateTime(DateTime(1950,9,1,0), t, 2).utc_datetime == DateTime(1950,9,1,0)   # TDT
        @test ZonedDateTime(DateTime(1950,9,1,1), t, 2).utc_datetime == DateTime(1950,9,1,1)   # TDT

        @test ZonedDateTime(DateTime(1950,9,1,2), t).utc_datetime == DateTime(1950,9,1,2)      # TDT
    end

    @testset "highly ambiguous hour" begin
        # Ambiguous local DateTime that has more than 2 solutions. Never occurs in reality.
        t = VariableTimeZone("Testing", [
            Transition(DateTime(1800,1,1), FixedTimeZone("TST",0,0)),
            Transition(DateTime(1960,4,1), FixedTimeZone("TDT",0,7200)),
            Transition(DateTime(1960,8,31,23), FixedTimeZone("TXT",0,3600)),
            Transition(DateTime(1960,9,1), FixedTimeZone("TST",0,0)),
        ])

        @test ZonedDateTime(DateTime(1960,8,31,23), t).utc_datetime == DateTime(1960,8,31,21)  # TDT
        @test ZonedDateTime(DateTime(1960,9,1,0), t, 1).utc_datetime == DateTime(1960,8,31,22) # TDT
        @test ZonedDateTime(DateTime(1960,9,1,0), t, 2).utc_datetime == DateTime(1960,8,31,23) # TXT
        @test ZonedDateTime(DateTime(1960,9,1,0), t, 3).utc_datetime == DateTime(1960,9,1,0)   # TST
        @test ZonedDateTime(DateTime(1960,9,1,1), t).utc_datetime == DateTime(1960,9,1,1)      # TST

        @test_throws AmbiguousTimeError ZonedDateTime(DateTime(1960,9,1,0), t, true)
        @test_throws AmbiguousTimeError ZonedDateTime(DateTime(1960,9,1,0), t, false)
    end

    @testset "skip an entire day" begin
        # Significant offset change: -11:00 -> 13:00.
        apia = first(compile("Pacific/Apia", tzdata["australasia"]))

        # Skips an entire day.
        @test ZonedDateTime(DateTime(2011,12,29,23),apia).utc_datetime == DateTime(2011,12,30,9)
        @test_throws NonExistentTimeError ZonedDateTime(DateTime(2011,12,30,0),apia)
        @test_throws NonExistentTimeError ZonedDateTime(DateTime(2011,12,30,23),apia)
        @test ZonedDateTime(DateTime(2011,12,31,0),apia).utc_datetime == DateTime(2011,12,30,10)
    end

    @testset "redundant transitions" begin
        # Redundant transitions should be ignored.
        # Note: that this can occur in reality if the TZ database parse has a Zone that ends
        # at the same time a Rule starts. When this occurs the duplicates always in standard
        # time with the same abbreviation.
        zone = Dict{AbstractString,FixedTimeZone}()
        zone["DTST"] = FixedTimeZone("DTST", 0, 0)
        zone["DTDT-1"] = FixedTimeZone("DTDT-1", 0, 3600)
        zone["DTDT-2"] = FixedTimeZone("DTDT-2", 0, 3600)

        dup = VariableTimeZone("DuplicateTest", [
            Transition(DateTime(1800,1,1), zone["DTST"])
            Transition(DateTime(1935,4,1), zone["DTDT-1"])  # Ignored
            Transition(DateTime(1935,4,1), zone["DTDT-2"])
            Transition(DateTime(1935,9,1), zone["DTST"])
        ])

        # Make sure that the duplicated hour only doesn't contain an additional entry.
        @test_throws AmbiguousTimeError ZonedDateTime(DateTime(1935,9,1), dup)
        @test ZonedDateTime(DateTime(1935,9,1), dup, 1).zone.name == "DTDT-2"
        @test ZonedDateTime(DateTime(1935,9,1), dup, 2).zone.name == "DTST"
        @test_throws BoundsError ZonedDateTime(DateTime(1935,9,1), dup, 3)

        # Ensure that DTDT-1 is completely ignored.
        @test_throws NonExistentTimeError ZonedDateTime(DateTime(1935,4,1), dup)
        @test ZonedDateTime(DateTime(1935,4,1,1), dup).zone.name == "DTDT-2"
        @test ZonedDateTime(DateTime(1935,8,31,23), dup).zone.name == "DTDT-2"
    end

    @testset "equality" begin
        # Check equality between ZonedDateTimes
        apia = first(compile("Pacific/Apia", tzdata["australasia"]))

        spring_utc = ZonedDateTime(DateTime(2010, 5, 1, 12), utc)
        spring_apia = ZonedDateTime(DateTime(2010, 5, 1, 1), apia)

        @test spring_utc.zone == FixedTimeZone("UTC", 0, 0)
        @test spring_apia.zone == FixedTimeZone("SST", -39600, 0)
        @test spring_utc == spring_apia
        @test spring_utc !== spring_apia
        @test isequal(spring_utc, spring_apia)
        @test hash(spring_utc) == hash(spring_apia)
        @test astimezone(spring_utc, apia) === spring_apia  # Since ZonedDateTime is immutable
        @test astimezone(spring_apia, utc) === spring_utc
        @test isequal(astimezone(spring_utc, apia), spring_apia)
        @test hash(astimezone(spring_utc, apia)) == hash(spring_apia)

        fall_utc = ZonedDateTime(DateTime(2010, 10, 1, 12), utc)
        fall_apia = ZonedDateTime(DateTime(2010, 10, 1, 2), apia)

        @test fall_utc.zone == FixedTimeZone("UTC", 0, 0)
        @test fall_apia.zone == FixedTimeZone("SDT", -39600, 3600)
        @test fall_utc == fall_apia
        @test fall_utc !== fall_apia
        @test !(fall_utc < fall_apia)
        @test !(fall_utc > fall_apia)
        @test isequal(fall_utc, fall_apia)
        @test hash(fall_utc) == hash(fall_apia)
        @test astimezone(fall_utc, apia) === fall_apia  # Since ZonedDateTime is immutable
        @test astimezone(fall_apia, utc) === fall_utc
        @test isequal(astimezone(fall_utc, apia), fall_apia)
        @test hash(astimezone(fall_utc, apia)) == hash(fall_apia)
    end

    @testset "broadcastable" begin
        # Validate that ZonedDateTime is treated as a scalar during broadcasting
        zdt = ZonedDateTime(2000, 1, 2, 3, utc)
        @test size(zdt .== zdt) == ()
    end

    @testset "deepcopy hash" begin
        # Issue #78
        x = ZonedDateTime(2017, 7, 6, 15, 44, 55, 28, warsaw)
        y = deepcopy(x)

        @test x == y
        @test x !== y
        @test !(x < y)
        @test !(x > y)
        @test isequal(x, y)
        @test hash(x) == hash(y)

        # Check that the ZonedDateTime and plain DateTimes don't hash to the same value.
        @test hash(x) != hash(y.utc_datetime)
    end

    @testset "multiple time zones" begin
        # ZonedDateTime constructor that takes any number of Period or TimeZone types
        @test_throws ArgumentError ZonedDateTime(FixedTimeZone("UTC", 0, 0), FixedTimeZone("TMW", 86400, 0))
    end

    @testset "cutoff" begin
        # The absolutely min DateTime you can create. Even smaller than `typemin(DateTime)`
        early_utc = ZonedDateTime(DateTime(UTM(typemin(Int64))), utc)

        # A FixedTimeZone is effective for all of time where as a VariableTimeZone has as
        # start.
        @test DateTime(early_utc, UTC) < warsaw.transitions[1].utc_datetime
        @test_throws NonExistentTimeError astimezone(early_utc, warsaw)
    end

    @testset "UnhandledTimeError" begin
        # VariableTimeZone with a cutoff set
        cutoff_tz = VariableTimeZone(
            "cutoff", [Transition(DateTime(1970, 1, 1), utc)], DateTime(1988, 5, 6),
        )

        ZonedDateTime(DateTime(1970, 1, 1), cutoff_tz)  # pre cutoff
        @test_throws UnhandledTimeError ZonedDateTime(DateTime(1988, 5, 6), cutoff_tz)  # on cutoff
        @test_throws UnhandledTimeError ZonedDateTime(DateTime(1989, 5, 7), cutoff_tz)
        @test_throws UnhandledTimeError ZonedDateTime(DateTime(1988, 5, 5), cutoff_tz) + Hour(24)

        zdt = ZonedDateTime(DateTime(2038, 3, 28), warsaw, from_utc=true)
        @test_throws UnhandledTimeError zdt + Hour(1)
    end

    @testset "no cutoff" begin
        # TimeZones that no longer have any transitions after the max_year shouldn't have a cutoff
        # eg. Asia/Hong_Kong, Pacific/Honolulu, Australia/Perth
        perth = first(compile("Australia/Perth", tzdata["australasia"]))
        zdt = ZonedDateTime(DateTime(2200, 1, 1), perth, from_utc=true)
    end

    @testset "convenience constructors" begin
        # Convenience constructors for making a DateTime on-the-fly
        digits = [2010, 1, 2, 3, 4, 5, 6]
        for i in eachindex(digits)
            @test ZonedDateTime(digits[1:i]..., warsaw) == ZonedDateTime(DateTime(digits[1:i]...), warsaw)
            @test ZonedDateTime(digits[1:i]..., utc) == ZonedDateTime(DateTime(digits[1:i]...), utc)
        end

        # Convenience constructor dealing with ambiguous time
        digits = [1916, 10, 1, 0, 2, 3, 4]  # Fall DST transition in Europe/Warsaw
        for i in eachindex(digits)
            expected = [
                ZonedDateTime(DateTime(digits[1:i]...), warsaw, 1)
                ZonedDateTime(DateTime(digits[1:i]...), warsaw, 2)
            ]

            if i > 1
                @test_throws AmbiguousTimeError ZonedDateTime(digits[1:i]..., warsaw)
            end

            @test ZonedDateTime(digits[1:i]..., warsaw, 1) == expected[1]
            @test ZonedDateTime(digits[1:i]..., warsaw, 2) == expected[2]
            @test ZonedDateTime(digits[1:i]..., warsaw, true) == expected[1]
            @test ZonedDateTime(digits[1:i]..., warsaw, false) == expected[2]
        end
    end

    @testset "promotion" begin
        @test_throws ErrorException promote_type(ZonedDateTime, Date)
        @test_throws ErrorException promote_type(ZonedDateTime, DateTime)
        @test_throws ErrorException promote_type(Date, ZonedDateTime)
        @test_throws ErrorException promote_type(DateTime, ZonedDateTime)
        @test promote_type(ZonedDateTime, ZonedDateTime) == ZonedDateTime

        # Issue #52
        dt = now()
        @test_throws ErrorException ZonedDateTime(dt, warsaw) > dt
    end

    @testset "extrema" begin
        @test typemin(ZonedDateTime) <= ZonedDateTime(typemin(DateTime), utc)
        @test typemax(ZonedDateTime) >= ZonedDateTime(typemax(DateTime), utc)
    end
end


text:
"""
    BipartiteView{T<:Real} <: AbstractGraph{Int64}

Represents a bipartite view of a hypergraph `h`.
Note this is a view - changes to the original hypergraph will be automatically reflected in the view.

**Constructors**

BipartiteView(::Hypergraph)

The bipartite view of a hypergraph is suitable for processing with the LightGraphs.jl package.
Several LightGraphs methods are provided for the compability.

"""
struct BipartiteView{T<:Real} <: AbstractGraph{Int}
    h::Hypergraph{T}
end

LightGraphs.nv(b::BipartiteView) = length(b.h.v2he)+length(b.h.he2v)

LightGraphs.vertices(b::BipartiteView) = Base.OneTo(LightGraphs.nv(b))

LightGraphs.ne(b::BipartiteView) = 2*sum(length.(b.h.v2he))

function LightGraphs.all_neighbors(b::BipartiteView, v::Integer)
    n1 = length(b.h.v2he)
    if v <= n1
      n1 .+ keys(b.h.v2he[v])
    else
      keys(b.h.he2v[v-n1])
    end
end

function LightGraphs.has_edge(b::BipartiteView, s, d)
    n1 = length(b.h.v2he)
    if s <= n1
        d > n1 && has_key(b.v2he[s], d - n1)
    else
        d <= n1 && has_key(b.h.he2v[s-n1], d)
    end
end

LightGraphs.outneighbors(b::BipartiteView, v::Integer) = LightGraphs.all_neighbors(b::BipartiteView, v)

LightGraphs.inneighbors(b::BipartiteView, v::Integer) = LightGraphs.all_neighbors(b::BipartiteView, v)

"""
    LightGraphs.SimpleGraph(b::BipartiteView)

Creates a `LightGraphs.SimpleGraph` representation of a `BipartiteView` b.

This creates a copy of the date. Note that the weights information is not stored
in the created `SimpleGraph`.
"""
function LightGraphs.SimpleGraph(b::BipartiteView)
    g = SimpleGraph(nv(b))
    for v in keys(b.h.v2he)
        for he in keys(b.h.v2he[v])
            add_edge!(g, v, length(b.h.v2he) + he)
        end
    end
    g
end

LightGraphs.is_directed(b::BipartiteView) = false

"""
    shortest_path(b::BipartiteView,source::Int, target::Int)

Finds a single shortest path in a graph `b` between vertices
`source` and `target`.
Note that if several paths of the same length exist, only one
will be returned.

"""
function shortest_path(b::BipartiteView,source::Int, target::Int)
    checkbounds(b.h.v2he, source)
    checkbounds(b.h.v2he, target)
    dj = dijkstra_shortest_paths(b, source)
    enumerate_paths(dj)[target][1:2:end]
end


text:
using MonteCarloMeasurements, ForwardDiff, Test
const FD = ForwardDiff

@testset "forwarddiff" begin
    c = 1 ± 0.1 # These are the uncertain parameters
    d = 1 ± 0.1 # These are the uncertain parameters
    # In the cost function below, we ensure that $cx+dy > 10 \; ∀ \; c,d ∈ P$ by looking at the worst case
    function cost(params)
        x,y = params
        -(3x+2y) + 10000sum(params .< 0) + 10000*(maximum(c*x+d*y) > 10)
    end

    params = [1., 2] # Initial guess
    paramsp = [1., 2] .± 0.001 # Initial guess
    @test cost(params) == -7    # Try the cost function

    @test FD.gradient(cost, params) == @unsafe FD.gradient(cost, paramsp)
    @test FD.gradient(sum, params) == FD.gradient(sum, paramsp)
    @test all(FD.gradient(prod, params) .≈ FD.gradient(prod, paramsp))

    unsafe_comparisons(false)
    @test FD.gradient(x -> params'x, params) == params
    @test FD.gradient(x -> paramsp'x, params) == paramsp
    @test FD.gradient(x -> params'x, paramsp) == params
    r = FD.gradient(x -> paramsp'x, paramsp)
    @test mean(mean(r[1])) ≈ params[1] atol=1e-2
    @test mean(mean(r[2])) ≈ params[2] atol=1e-2

    @test FD.jacobian(x -> params+x, params) == I
    @test FD.jacobian(x -> paramsp+x, params) == I
    @test FD.jacobian(x -> params+x, paramsp) == I

    @test FD.jacobian(x -> params-x, params) == -I
    @test FD.jacobian(x -> paramsp-x, params) == -I
    @test FD.jacobian(x -> params-x, paramsp) == -I

    @test FD.jacobian(x -> x-params, params) == I
    @test FD.jacobian(x -> x-paramsp, params) == I
    @test FD.jacobian(x -> x-params, paramsp) == I

    function strange(x,y)
        (x.^2)'*(y.^2)
    end
    ref = FD.gradient(x->strange(x,params), params)
    FD.gradient(x->strange(x,params), paramsp) != ref
    FD.gradient(x->strange(x,params), paramsp) ≈ ref
    FD.gradient(x->strange(x,paramsp), params) != ref
    FD.gradient(x->strange(x,paramsp), params) ≈ ref
    r = FD.gradient(x->strange(x,paramsp), paramsp) # maybe this is a bit overkill
    @test mean(mean(r[1])) ≈ ref[1] atol=1e-2
    @test mean(mean(r[2])) ≈ ref[2] atol=1e-2
    @test mean(mean(r[1])) != ref[1]
    @test mean(mean(r[2])) != ref[2]

    unsafe_comparisons(false)
end


text:
abstract type VTPRRuleModel <: MixingRule end

struct VTPRRule{γ} <: VTPRRuleModel
    components::Array{String,1}
    activity::γ
end

@registermodel VTPRRule
export VTPRRule
function VTPRRule(components::Vector{String}; activity = Wilson, userlocations::Vector{String}=String[],activity_userlocations::Vector{String}=String[], verbose::Bool=false)
    init_activity = activity(components;userlocations = activity_userlocations,verbose)

    model = VTPRRule(components, init_activity)
    return model
end

function mixing_rule(model::PRModel,V,T,z,mixing_model::VTPRRuleModel,α,a,b,c)
    n = sum(z)
    x = z./n
    invn2 = (one(n)/n)^2
    lnγ_res_ = lnγ_res(mixing_model.activity,1e5,T,z)
    g_E_res = sum(x[i]*R̄*T*lnγ_res_[i] for i ∈ @comps)
    b = Diagonal(b).diag
    b = ((b.^(3/4).+b'.^(3/4))/2).^(4/3)
    b̄ = dot(z,Symmetric(b),z) * invn2
    c̄ = dot(z,c)/n
    ā = b̄*R̄*T*(sum(x[i]*a[i,i]*α[i]/b[i,i]/(R̄*T) for i ∈ @comps)-1/0.53087*(g_E_res/(R̄*T)))
    return ā,b̄,c̄
end


text:
# This file is a part of SearchModels.jl

include("polynomial.jl")
include("rosenbrock.jl")


text:
function SpecialFunctions.besselj(ν::Arblib.ArbOrRef, z::Arblib.ArbSeries)
    deg = Arblib.degree(z)
    res = zero(z)

    x = z[0]
    tmp = zero(x)
    ai = zero(x)

    # ai = besselj(ν, x)
    Arblib.hypgeom_bessel_j!(ai, ν, x)

    res[0] = ai

    if deg >= 1
        # ai = (besselj(ν - 1, x) - besselj(ν + 1, x)) / 2

        # ai = besselj(ν - 1, x)
        Arblib.sub!(tmp, ν, 1)
        Arblib.hypgeom_bessel_j!(ai, tmp, x)

        # tmp = besselj(ν + 1, x)
        # ai += tmp
        Arblib.add!(tmp, ν, 1)
        Arblib.hypgeom_bessel_j!(tmp, tmp, x)
        Arblib.sub!(ai, ai, tmp)

        # ai /= 2
        Arblib.mul_2exp!(ai, ai, -1)

        res[1] = ai

        if deg >= 2
            # ai = (besselj(ν - 2, x) + besselj(ν + 2, x) - 2a0) / 8

            # ai = besselj(ν - 2, x)
            Arblib.sub!(tmp, ν, 2)
            Arblib.hypgeom_bessel_j!(ai, tmp, x)

            # tmp = besselj(ν + 2, x)
            # ai += tmp
            Arblib.add!(tmp, ν, 2)
            Arblib.hypgeom_bessel_j!(tmp, tmp, x)
            Arblib.add!(ai, ai, tmp)

            # ai -= 2a0
            Arblib.mul_2exp!(tmp, Arblib.ref(res, 0), 1)
            Arblib.sub!(ai, ai, tmp)

            # ai /= 8
            Arblib.mul_2exp!(ai, ai, -3)

            res[2] = ai

            if deg >= 3
                # ai = (besselj(ν - 3, x) - besselj(ν + 3, x) - 6a1) / 48

                # ai = besselj(ν - 3, x)
                Arblib.sub!(tmp, ν, 3)
                Arblib.hypgeom_bessel_j!(ai, tmp, x)

                # tmp = besselj(ν + 3, x)
                # ai += tmp
                Arblib.add!(tmp, ν, 3)
                Arblib.hypgeom_bessel_j!(tmp, tmp, x)
                Arblib.sub!(ai, ai, tmp)

                # ai -= 6a1
                Arblib.mul!(tmp, Arblib.ref(res, 1), 6)
                Arblib.sub!(ai, ai, tmp)

                # ai /= 48
                Arblib.div!(ai, ai, 48)

                res[3] = ai

                if deg >= 4
                    x2 = x^2
                    x2mν2 = x2 - ν^2
                    for i = 4:deg
                        k = i - 4

                        # tmp = 2x * res[k + 1]
                        Arblib.mul_2exp!(tmp, x, 1)
                        Arblib.mul!(tmp, tmp, Arblib.ref(res, k + 1))

                        # ai = res[k] + tmp
                        Arblib.add!(ai, Arblib.ref(res, k), tmp)

                        # tmp = (k^2 + 4k + 4 - ν^2 + x^2)
                        # ai += tmp * res[k + 2]
                        Arblib.add!(tmp, x2mν2, k^2 + 4k + 4)
                        Arblib.addmul!(ai, tmp, Arblib.ref(res, k + 2))

                        # tmp = (2k^2 + 11k + 15) * x
                        # ai += tmp * res[k + 3]
                        Arblib.mul!(tmp, x, 2k^2 + 11k + 15)
                        Arblib.addmul!(ai, tmp, Arblib.ref(res, k + 3))

                        # tmp = x^2 * (k^2 + 7k + 12)
                        # ai /= tmp
                        Arblib.mul!(tmp, x2, k^2 + 7k + 12)
                        Arblib.div!(ai, ai, tmp)

                        # ai = -ai
                        Arblib.neg!(ai, ai)

                        res[i] = ai
                    end
                end
            end
        end
    end

    # Compose the Taylor series for the Bessel function with that of z

    # First set the constant in z to 0 temporarily
    Arblib.zero!(Arblib.ref(z, 0))
    Arblib.compose_series!(res, res, z, deg + 1)
    Arblib.set!(Arblib.ref(z, 0), x)

    return res
end

SpecialFunctions.besselj0(z::ArbSeries) = besselj(zero(Arb), z)
SpecialFunctions.besselj1(z::ArbSeries) = besselj(one(Arb), z)


text:
module Selenium

using PyCall

"Sets up the Python selenium module."
function __init__()
    py"""
    import selenium.webdriver as webdriver
    """
    global webdriver = py"webdriver"
end

"Spawns a Firefox browser controlled by Selenium."
function spawn_driver(;timewait::Int=0, headless::Bool=false)
    opts = headless && begin
        opts = webdriver.FirefoxOptions()
        opts.set_headless()
	opts
    end
    if opts != false
        driver = webdriver.Firefox(firefox_options=opts)
    else
        driver = webdriver.Firefox()
    end
    if timewait > 0
	    driver.implicitly_wait(timewait)
    end
    return driver
end

export spawn_driver

end

max_stars_repo_name: stivenroytman/Selenium

max_issues_repo_name: stivenroytman/Selenium

max_forks_repo_name: stivenroytman/Selenium


text:
function normal_form_f4!(
            ring::PolyRing,
            basis::Basis{C},
            ht::MonomialHashtable,
            tobereduced::Basis{C}) where {C<:Coeff}

    matrix = initialize_matrix(ring, C)
    symbol_ht = initialize_secondary_hash_table(ht)

    select_tobereduced!(basis, tobereduced, matrix, symbol_ht, ht)

    symbolic_preprocessing!(basis, matrix, ht, symbol_ht)

    convert_hashes_to_columns!(matrix, symbol_ht)

    #dump(matrix, maxdepth=5)
    #@error "((()))"

    sort_matrix_rows_decreasing!(matrix) # for pivots,  AB part

    #=
    dump(matrix, maxdepth=5)
    @error "UPROWS"
    for (i, u) in enumerate(matrix.uprows)
        print(map(c -> symbol_ht.exponents[matrix.col2hash[c]], u))
        println(" ", basis.coeffs[matrix.up2coef[i]])
    end
    @error "LOWrows"
    for (i,u) in enumerate(matrix.lowrows)
        println(map(c -> symbol_ht.exponents[matrix.col2hash[c]], u))
        println(" ", tobereduced.coeffs[matrix.low2coef[i]])
    end

    @error "" Int(ring.ch)
    =#

    exact_sparse_linear_algebra_nf!(matrix, tobereduced, basis)

    # @warn "lll"
    # dump(matrix, maxdepth=5)

    convert_nf_rows_to_basis_elements!(matrix, tobereduced, ht, symbol_ht)

    tobereduced
end

function normal_form_f4(
                ring::PolyRing,
                basisexps::Vector{Vector{ExponentVector}},
                basiscoeffs::Vector{Vector{C}},
                tobereducedexps::Vector{Vector{ExponentVector}},
                tobereducedcfs::Vector{Vector{C}},
                rng::Rng,
                tablesize::Int=2^16) where {Rng, C<:Coeff}

    basis, ht = initialize_structures(ring, basisexps,
                                        basiscoeffs, rng, tablesize)

    tobereduced, ht = initialize_structures(ring, tobereducedexps,
                                        tobereducedcfs, rng, tablesize, ht)

    tobereduced = normal_form_f4!(ring, basis, ht, tobereduced)

    export_basis_data(tobereduced, ht)
end


text:
################################################################################
#
#	ABSTRACT TYPE
#
#   Site{L,D}
#   --> L is the label type
#   --> D is the dimension of embedding space
#
#   FILE CONTAINS
#       - abstract type definition
#       - interface definition
#       - interface testing
#
################################################################################

################################################################################
#
#   ABSTRACT TYPE DEFINITION
#
################################################################################
abstract type AbstractSite{L,D} end

# export the type
export AbstractSite

# DOCSTRING
"""
    abstract type AbstractSite{L,D}

The abstract site type that describes all site implementations. It is parametric in two types,
- `L` the type of the bond label
- `D` an integer number specifying the number of spatial dimensions in which the site is embedded

The abstract type is equipped with interface functions which are listed in the following. These interface functions have to be implemented by concrete subtypes / structs for the concrete type to be usable in all pre-implemented functions.
- `newSite` to create a new site object of a passed concrete site type
- `label` and `label!` to access the site label
- `point` and `point!` to access the site position in real space

The interface functions throw adequate errors when called but not implemented for concrete types.
See specific documentation for the individual functions.

# Examples

```julia-REPL
julia> AbstractSite{Symbol,3} <: AbstractSite{L,D} where{L,D}
true
```

"""
AbstractSite

"""
Number of spatial dimensions the site is embedded in.
"""
Base.ndims(::Type{S}) where {D,S<:AbstractSite{L,D} where L} = D

################################################################################
#
#	INTERFACING / ACCESSING SITES
#	(functions have to be overwritten by concrete types)
#
################################################################################

# default constructor interface
# used for creation of new sites
"""
    function newSite(
            :: Type{S},
            point   :: Vector{<:Real},
            label   :: L
        ) :: S where {L,S<:AbstractSite{L,D} where D}

Interface function for creation of new `AbstractSite` object of a passed type `S` from a passed `point` and `label`.
Returns a new site object of type `S`.

This function has to be overwritten by a concrete type,
otherwise it will throw an error when used in further functions.

# Examples

```julia-REPL
julia> s = newSite(Site{String,2}, [0,0], "site1")
...
```
"""
function newSite(
            :: Type{S},
            point   :: Vector{<:Real},
            label   :: L
        ) :: S where {L,S<:AbstractSite{L,D} where D}

    # print an error because implementation for concrete type is missing
    error(  "not implemented function 'newSite' for concrete site type " *
            string(S) * " with label type " * string(L) )
end

# export the newSite interface
export newSite

# get label
"""
    function label(
            s :: AbstractSite{L,D}
        ) :: L where {L,D}

Interface function for obtaining the label of a passed `AbstractSite` object `s`.
Returns an object of type `L`.

This function has to be overwritten by a concrete type,
otherwise it will throw an error when used in further functions.

# Examples

```julia-REPL
julia> s = newSite(Site{String,2}, [0,0], "site1")

julia> label(s)
"site1"
```
"""
function label(
            s :: AbstractSite{L,D}
        ) :: L where {L,D}

    # print an error because implementation for concrete type is missing
    error("not implemented interface function 'label' for site type " * string(typeof(s)))
end
# set label
"""
    function label!(
            s :: AbstractSite{L,D},
            l :: L
        ) where {L,D}

Interface function for setting the label of a passed `AbstractSite` object `s`
to a passed value `l` of label type `L`. Returns nothing.

This function has to be overwritten by a concrete type,
otherwise it will throw an error when used in further functions.

# Examples

```julia-REPL
julia> s = newSite(Site{String,2}, [0,0], "site1")

julia> label!(s, "mylabel")

julia> label(s)
"mylabel"
```
"""
function label!(
            s :: AbstractSite{L,D},
            l :: L
        ) where {L,D}

    # print an error because implementation for concrete type is missing
    error("not implemented interface function 'label!' for site type " * string(typeof(s)))
end

# export the label interface
export label, label!

# get point
"""
    function point(
            s :: AbstractSite{L,D}
        ) ::Vector{Float64} where {L,D}

Interface function for obtaining the position of a passed `AbstractSite` object `s`.
Returns an object of type `Vector{Float64}` with length `D`.

This function has to be overwritten by a concrete type,
otherwise it will throw an error when used in further functions.

# Examples

```julia-REPL
julia> s = newSite(Site{String,2}, [0,0], "site1")

julia> point(s)
[0.0, 0.0]
```
"""
function point(
            s :: AbstractSite{L,D}
        ) :: Vector{Float64} where {L,D}

    # print an error because implementation for concrete type is missing
    error("not implemented interface function 'point' for site type " * string(typeof(s)))
end
# set point
"""
    function point!(
            s :: AbstractSite{L,D},
            p :: Vector{<:Real}
        ) where {L,D}

Interface function for setting the position of a passed `AbstractSite` object `s`
to a passed value `p` of label type `Vector{<:Real}`. Returns nothing.

This function has to be overwritten by a concrete type,
otherwise it will throw an error when used in further functions.

# Examples

```julia-REPL
julia> s = newSite(Site{String,2}, [0,0], "site1")

julia> point!(s, [2.0,1.0])

julia> point(s)
[2.0,1.0]
```
"""
function point!(
            s :: AbstractSite{L,D},
            p :: Vector{<:Real}
        ) where {L,D}

    # print an error because implementation for concrete type is missing
    error("not implemented interface function 'point!' for site type " * string(typeof(s)))
end

# export the point interface
export point, point!

function saveSites(
        ss :: Vector{S},
        fn :: AbstractString,
        group :: AbstractString = "sites"
        ;
        append :: Bool = false
    ) where {LS,D,S<:AbstractSite{LS,D}}

    # print an error because implementation for concrete type is missing
    error("not implemented interface function 'saveSites' for site type " * string(S))
end

function loadSites(
        ::Type{S},
        fn :: AbstractString,
        group :: AbstractString = "sites"
    ) where {LS,D,S<:AbstractSite{LS,D}}

    # print an error because implementation for concrete type is missing
    error("not implemented interface function 'loadSites' for site type " * string(S))
end

export saveSites, loadSites

# SIMILAR FUNCTION (can be overwritten but does not have to be overwritten)

# without new parameters
function similar(
            s :: S
        ) :: S where {L,D,S<:AbstractSite{L,D}}

    # return a new site object
    return newSite(S, deepcopy(point(s)), deepcopy(label(s)))
end
# with new parameters
function similar(
            s :: S,
            p :: Vector{<:Real},
            l :: L
        ) :: S where {L,D,S<:AbstractSite{L,D}}

    # create a new site object
    s_new = similar(s)
    # set parameters
    point!(s_new, p)
    label!(s_new, l)
    # return the new object
    return s_new
end

# export the similar interface
export similar

################################################################################
#
#	TESTING THE INTERFACE OF SITES
#	(should never be overwritten by concrete types)
#
################################################################################

# TESTING THE SITE INTERFACE
function testInterface(
            ::Type{T}
        ) :: Bool where {T<:AbstractSite}

	# get the parameterless constructor
	S = Base.typename(T).wrapper

    # iterate over some standard points
    for p in Vector{Float64}[[1.0,], [1.0, 1.0], [0.0, 0.0, 0.0, 0.0]]
    # iterate over some standard labels
    for l in ["t", 1, 1.0]
        # create a new site
        s = newSite(S{typeof(l), length(p)}, p,l)
        # test the interface by getting label and point
		label(s)
		point(s)
        # set both label and point
        label!(s, l)
        point!(s, p)
    end
    end

	# return true to indicate the test passed
	return true
end

# export the test functions
export testInterface


text:
export
    AccelTurnrate,
    AccelDesang,
    LatLonAccel,
    AccelSteeringAngle

include("accel_turnrate.jl")
include("accel_desang.jl")
include("lat_lon_accel.jl")
include("accel_steering.jl")


text:
function calcbranchprops!(system::CVSystem,old=Dict("a"=>0),restart="no")
    if restart == "no"
        for i in 1:length(system.branches.ID)
            # reference area, m^2
            push!(system.branches.A0,[π*(system.branches.radiusincm[i]*cmTom)^2])
            # stiffness parameter, Pa/m
            push!(system.branches.beta,[(sqrt(π)*system.branches.thicknessincm[i]*
                cmTom*system.branches.YoungsModinMPa[i]*MPaToPa)/
                ((1-system.solverparams.nu^2)*system.branches.A0[i][end])])
            # Moens-Korteweg wave speed, m/s
            push!(system.branches.c0,[sqrt(0.5*system.branches.beta[i][end]/
                system.solverparams.rho)*system.branches.A0[i][end]^0.25])
        end
    elseif restart == "yes"
        for i = 1:length(system.branches.ID)
            push!(system.branches.A0,[old["A0"][i][end]]);
            push!(system.branches.beta,[old["beta"][i][end]]);
            push!(system.branches.c0,[old["c0"][i][end]]);
        end
    end
end

max_stars_repo_name: dcanuto/CardioModeling

max_issues_repo_name: dcanuto/CardioModeling

max_forks_repo_name: dcanuto/CardioModeling


text:
"""
    IdealPDController

This piece of software relies on a truth sensor and an ideal actuator to align
the spacecraft with the ICRF.

"""
function IdealPDController()

    function step(t, constants, state, inputs, outputs_bus,
                  inputs_bus) # for output only

        # Extract what we care about.
        ω_BI_B = outputs_bus["truth_sensor"].ω_BI_B
        q_BI   = outputs_bus["truth_sensor"].q_BI
        q_TI   = constants # Target quaternion

        # Run the actual algorithm that will fly.
        (f_B, τ_B) = pd_controller(q_TI, q_BI, ω_BI_B)

        # Write out our command for the ideal actuator.
        inputs_bus["ideal_actuator"] = IdealActuatorCommand(f_B, τ_B)

        return (state,      # Updated state (none)
                nothing,    # Outputs
                inputs_bus, # Updated inputs bus (this is not necessary)
                true)       # Active

    end

    DynamicalModel("ideal_pd_controller",
                   init = step,
                   update = step,
                   timing = ModelTiming(0.05),
                   constants = [0.; 0.; 0.; 1.]) # target quaternion

end

# This represents the "real" algorithm under development.
function pd_controller(q_TI, q_BI, ω_BI_B)
    theta, r_B = q2aa(qdiff(q_TI, q_BI))
    τ_B        = 0.2 * theta * r_B - 2 * ω_BI_B
    f_B        = [0.; 0.; 0.]
    return (f_B, τ_B)
end

"""
ReactionWheelController

This piece of software relies on a truth sensor and an ideal actuator to align
the spacecraft with the ICRF.

"""
function ReactionWheelController()

    function step(t, q_TI, state, inputs, outputs_bus,
                  inputs_bus) # for output only

        # Extract what we care about.
        q_BI   = outputs_bus["star_tracker"]
        ω_BI_B = outputs_bus["gyro"]
        ω_rw   = [outputs_bus["rw1"]; outputs_bus["rw2"]; outputs_bus["rw3"]]

        # Get the desired net torque on the body.
        (f_B, τ_B) = pd_controller(q_TI, q_BI, ω_BI_B)

        # Feed the dynamics forward to get the resulting torque necessary from the RW motors.
        I_b = 5.
        I_r = 0.25
        τ_br_B_cmd = (  τ_B
                      # + ω_BI_B × (I_b * ω_BI_B) # Zero (MOI are all the same)
                      + ω_BI_B × (I_r * ω_rw))

        # H_B = (  I_b * ω_BI_B
        #        + 0.25*[1.; 1/sqrt(2.); 1/sqrt(2.)] .* ([ω_rw[1]; 0.; 0.] + ω_BI_B)
        #        + 0.25*[1/sqrt(2.); 1.; 1/sqrt(2.)] .* ([0.; ω_rw[2]; 0.] + ω_BI_B)
        #        + 0.25*[1/sqrt(2.); 1/sqrt(2.); 1.] .* ([0.; 0.; ω_rw[3]] + ω_BI_B))
        # H_I = qrot(qinv(q_BI), H_B)
        # println("Angular momentum in inertial frame: ", H_I)

        # Write out our command for the ideal actuator.
        inputs_bus["rw1"] = τ_br_B_cmd[1]
        inputs_bus["rw2"] = τ_br_B_cmd[2]
        inputs_bus["rw3"] = τ_br_B_cmd[3]

        return (state,      # Updated state (none)
                nothing,    # Outputs
                inputs_bus, # Updated inputs bus (this is not necessary)
                true)       # Active

    end

    DynamicalModel("rw_controller",
                   init = step,
                   update = step,
                   timing = ModelTiming(0.05),
                   constants = [0.; 0.; 0.; 1.]) # target quaternion

end


text:
function solve(d, f)
    p, i = 1, 0
    while p <= length(d)
        v = d[p]
        d[p] = f(v)
        p += v
        i += 1
    end
    return i
end

d = readdlm("input.txt", Int)
println(solve(copy(d), x -> x + 1))
println(solve(d, x -> x + (x < 3 ? 1 : -1)))


text:
# Classification.jl

module Classification

# Classification related code here

end # module


text:
module LearnJulia

include("PersonalAssets.jl")
include("ExampleSubmodule.jl")

using .ExampleSubmodule
export anexamplefunction

"""
  helloworld()

A dreadfully important function which does a thing.
"""
function helloworld()
  print("hello world")
end

end

max_stars_repo_name: BradLyman/learn-julia

max_issues_repo_name: BradLyman/learn-julia

max_forks_repo_name: BradLyman/learn-julia


text:
#---------------------------------------------------------------------#
#This function computes the Initial & Exact Solutions.
#Written by F.X. Giraldo on April 19, 2019
#           Department of Applied Mathematics
#           Naval Postgraduate School
#           Monterey; CA 93943-5216
#---------------------------------------------------------------------#
function exact_solution(coord,Npoin,case,DFloat)

    #Set some constants
    c=2*π

    #Initialize
    qe=zeros(DFloat,Npoin)
    qeₓ=zeros(DFloat,Npoin)
    fe=zeros(DFloat,Npoin)

    #Generate Grid Points
    for I=1:Npoin
        x=coord[I]
        if (case == 1) #Gaussian
            qe[I]=sin(c*x)
            qeₓ[I]=-c*cos(c*x)
            fe[I]=-c^2*sin(c*x)
        end
    end #I
    return (qe,qeₓ,fe)
end


text:
using Test, ConventionalApp

#printsep(label, sep="-") = println("\n", label, "\n", repeat(sep, 80))
#printheader(label) = printsep(label, "=")

function show_testset_section()
	println()
	@info "SECTION: " * Test.get_testset().description * "\n" * repeat("=", 80)
end

function show_testset_description()
	@info "Testing: " * Test.get_testset().description
end

@testset "ConventionalApp tests" begin
	testfiles = ["test1.jl", "test2.jl"]

	for testfile in testfiles
#		include(testfile)
	end

	@test 1 == 1

end #testset

:Test_Complete

max_stars_repo_name: ma-laforge/JuliaApplications

max_issues_repo_name: ma-laforge/JuliaApplications

max_forks_repo_name: ma-laforge/JuliaApplications


text:
import Oceananigans.Grids: required_halo_size
using Oceananigans.Utils: prettysummary

struct ScalarDiffusivity{TD, F, N, K} <: AbstractScalarDiffusivity{TD, F}
    ν :: N
    κ :: K

    function ScalarDiffusivity{TD, F}(ν::N, κ::K) where {TD, F, N, K}
        return new{TD, F, N, K}(ν, κ)
    end
end

"""
    ScalarDiffusivity(time_disc=ExplicitTimeDiscretization,
                      formulation=ThreeDimensionalFormulation, FT=Float64;
                      ν=0, κ=0,
                      discrete_form = false)

Return `ScalarDiffusivity` with viscosity `ν` and tracer diffusivities `κ`
for each tracer field in `tracers`. If a single `κ` is provided, it is
applied to all tracers. Otherwise `κ` must be a `NamedTuple` with values
for every tracer individually.

`ν` and the fields of `κ` may be constants (converted to `FT`), arrays, fields or
    - functions of `(x, y, z, t)` if `discrete_form = false`
    - functions of `(i, j, k, grid, LX, LY, LZ)` with `LX`, `LY` and `LZ` are either `Face()` or `Center()` if
      `discrete_form = true`.
"""
function ScalarDiffusivity(time_disc=ExplicitTimeDiscretization(),
                           formulation=ThreeDimensionalFormulation(), FT=Float64;
                           ν=0, κ=0,
                           discrete_form = false)

    if formulation == HorizontalFormulation() && time_discretization == VerticallyImplicitTimeDiscretization()
        throw(ArgumentError("VerticallyImplicitTimeDiscretization is only supported for `HorizontalFormulation` or `ThreeDimensionalFormulation`"))
    end

    κ = convert_diffusivity(FT, κ; discrete_form)
    ν = convert_diffusivity(FT, ν; discrete_form)

    return ScalarDiffusivity{typeof(time_disc), typeof(formulation)}(ν, κ)
end

const VerticalScalarDiffusivity{TD} = ScalarDiffusivity{TD, VerticalFormulation} where TD
const HorizontalScalarDiffusivity{TD} = ScalarDiffusivity{TD, HorizontalFormulation} where TD

  VerticalScalarDiffusivity(time_disc=ExplicitTimeDiscretization(), FT::DataType=Float64; kwargs...) = ScalarDiffusivity(time_disc, VerticalFormulation(), FT; kwargs...)
HorizontalScalarDiffusivity(time_disc=ExplicitTimeDiscretization(), FT::DataType=Float64; kwargs...) = ScalarDiffusivity(time_disc, HorizontalFormulation(), FT; kwargs...)

# Aliases that allow specify the floating type, assuming that the discretization is Explicit in time
          ScalarDiffusivity(FT::DataType; kwargs...) = ScalarDiffusivity(ExplicitTimeDiscretization(), ThreeDimensionalFormulation(), FT; kwargs...)
  VerticalScalarDiffusivity(FT::DataType; kwargs...) = ScalarDiffusivity(ExplicitTimeDiscretization(), VerticalFormulation(), FT; kwargs...)
HorizontalScalarDiffusivity(FT::DataType; kwargs...) = ScalarDiffusivity(ExplicitTimeDiscretization(), HorizontalFormulation(), FT; kwargs...)

required_halo_size(closure::ScalarDiffusivity) = 1

function with_tracers(tracers, closure::ScalarDiffusivity{TD, F}) where {TD, F}
    κ = tracer_diffusivities(tracers, closure.κ)
    return ScalarDiffusivity{TD, F}(closure.ν, κ)
end

@inline viscosity(closure::ScalarDiffusivity, args...) = closure.ν
@inline diffusivity(closure::ScalarDiffusivity, ::Val{tracer_index}, args...) where tracer_index = closure.κ[tracer_index]

# z_viscosity and z_diffusivity fallback to viscosity and diffusivity
# for both ThreeDimensional and VerticalScalarDiffusivity
@inline z_viscosity(closure::HorizontalScalarDiffusivity, args...) = 0
@inline z_diffusivity(closure::HorizontalScalarDiffusivity, args...) = 0

calculate_diffusivities!(diffusivities, ::ScalarDiffusivity, args...) = nothing

# Note: we could compute ν and κ (if they are Field):
# function calculate_diffusivities!(diffusivities, closure::ScalarDiffusivity, args...)
#     compute!(viscosity(closure, diffusivities))
#     !isnothing(closure.κ) && Tuple(compute!(diffusivity(closure, Val(c), diffusivities) for c=1:length(closure.κ)))
#     return nothing
# end

function Base.summary(closure::ScalarDiffusivity)
    TD = summary(time_discretization(closure))
    prefix = replace(summary(formulation(closure)), "Formulation" => "")
    prefix === "ThreeDimensional" && (prefix = "")
    return string(prefix, "ScalarDiffusivity{$TD}(ν=", prettysummary(closure.ν), ", κ=", prettysummary(closure.κ), ")")
end

Base.show(io::IO, closure::ScalarDiffusivity) = print(io, summary(closure))


text:
mutable struct TraitValidationXMLElement <: MyXMLElement
    el::XMLOrNothing
    treeModel_el::TreeModelXMLElement
    traitLikelihood_el::TraitLikelihoodXMLElement
    mis_ind::Int
    obs_ind::Int
    standardize::Bool
    id::String

    function TraitValidationXMLElement(tm_el::TreeModelXMLElement,
                    tl_el::TraitLikelihoodXMLElement)

        @assert length(tm_el.node_traits) == 2

        mis_ind = 1
        obs_ind = 2

        return new(nothing, tm_el, tl_el, mis_ind, obs_ind, false, bn.TRAIT_VALIDATION)
    end
end

mutable struct FactorValidationXMLElement <: MyXMLElement
    el::XMLOrNothing
    ifl::IntegratedFactorsXMLElement
    tdl::TraitLikelihoodXMLElement
    tml::TreeModelXMLElement
    trait_name::String
    trait_param_name::String
    id::String
end

function FactorValidationXMLElement(ifl::IntegratedFactorsXMLElement,
                                    tdl::TraitLikelihoodXMLElement,
                                    tml::TreeModelXMLElement
                                    )

    @assert length(tml.param_names) == 2

    return FactorValidationXMLElement(nothing, ifl, tdl, tml,
                                      tml.node_traits[2], tml.param_names[2],
                                      "factorValidation")
end

mutable struct CrossValidationXMLElement <: MyXMLElement
    el::XMLOrNothing
    validationProvider_el::MyXMLElement
    log_sum::Bool
    validation_type::String
    id::String

    CrossValidationXMLElement(xml::MyXMLElement) =
            new(nothing, xml, false, bn.SQUARED_ERROR, bn.CROSS_VALIDATION)
end

function set_validation_type!(cv::CrossValidationXMLElement, t::String)
    cv.validation_type = t
end

function set_log_sum!(cv::CrossValidationXMLElement, b::Bool)
    cv.log_sum = b
end

function make_xml(cv_el::CrossValidationXMLElement)

    el = new_element(bn.CROSS_VALIDATION)
    set_attributes(el, [(bn.ID, cv_el.id),
                        (bn.LOG_SUM, string(cv_el.log_sum)),
                        (bn.TYPE, cv_el.validation_type)])

    make_xml(cv_el.validationProvider_el)
    add_child(el, cv_el.validationProvider_el.el)

    cv_el.el = el
    return el
end

function make_xml(tv_el::TraitValidationXMLElement)

    obs_ind = tv_el.obs_ind
    mis_ind = tv_el.mis_ind

    el = new_element(bn.TRAIT_VALIDATION_PROVIDER)
    attrs = [(bn.ID, bn.TRAIT_VALIDATION),
            (bn.TRAIT_NAME, tv_el.treeModel_el.node_traits[obs_ind]),
            (bn.INFERRED_TRAIT, tv_el.traitLikelihood_el.attrs[bn.TRAIT_NAME]),
            (bn.STANDARDIZE, string(tv_el.standardize))]
    set_attributes(el, attrs)

    make_xml(tv_el.traitLikelihood_el)
    add_ref_el(el, tv_el.traitLikelihood_el.el)

    tp_el = new_child(el, bn.TRAIT_PARAMETER)
    add_ref_el(tp_el, bn.PARAMETER, tv_el.treeModel_el.param_names[obs_ind])

    tv_el.el = el
    return el
end

function make_xml(fac_val::FactorValidationXMLElement)
    el = new_element(bn.FACTOR_VALIDATION)

    make_xml(fac_val.ifl)
    make_xml(fac_val.tdl)
    make_xml(fac_val.tml)

    add_ref_el(el, fac_val.ifl.el)
    add_ref_el(el, fac_val.tdl.el)
    trait_name = attribute(fac_val.tdl.el, bn.TRAIT_NAME)
    set_attribute(el, bn.ID, fac_val.id)
    set_attribute(el, bn.TRAIT_NAME, trait_name)

    trait_provider = new_child(el, bn.TRAIT_DATA_PROVIDER)
    add_ref_el(trait_provider, fac_val.tml.el)
    set_attribute(trait_provider, bn.TRAIT_NAME, fac_val.trait_name)
    tp_el = new_child(trait_provider, bn.TRAIT_PARAMETER)
    add_ref_el(tp_el, bn.PARAMETER, fac_val.trait_param_name)

    fac_val.el = el
end

function get_loggables(cv_el::CrossValidationXMLElement)
    make_xml(cv_el)
    return [cv_el.el]
end


text:
function MultiStartUnivariateOptimizer(arg0::UnivariateOptimizer, arg1::jint, arg2::RandomGenerator)
    return MultiStartUnivariateOptimizer((UnivariateOptimizer, jint, RandomGenerator), arg0, arg1, arg2)
end

function get_evaluations(obj::MultiStartUnivariateOptimizer)
    return jcall(obj, "getEvaluations", jint, ())
end

function get_optima(obj::MultiStartUnivariateOptimizer)
    return jcall(obj, "getOptima", Vector{UnivariatePointValuePair}, ())
end

function optimize(obj::MultiStartUnivariateOptimizer, arg0::Vector{OptimizationData})
    return jcall(obj, "optimize", UnivariatePointValuePair, (Vector{OptimizationData},), arg0)
end