text:
using FrankWolfe
using ProgressMeter
using Arpack
using Plots
using DoubleFloats
using ReverseDiff

# download movielens data
using ZipFile, DataFrames, CSV
import JSON

using Random

using Profile

using SparseArrays, LinearAlgebra

temp_zipfile = download("http://files.grouplens.org/datasets/movielens/ml-latest-small.zip")

# temp_zipfile = download("http://files.grouplens.org/datasets/movielens/ml-25m.zip")
#temp_zipfile = download("http://files.grouplens.org/datasets/movielens/ml-latest.zip")

zarchive = ZipFile.Reader(temp_zipfile)

movies_file = zarchive.files[findfirst(f -> occursin("movies", f.name), zarchive.files)]
movies_frame = CSV.read(movies_file, DataFrame)

ratings_file = zarchive.files[findfirst(f -> occursin("ratings", f.name), zarchive.files)]
ratings_frame = CSV.read(ratings_file, DataFrame)

# ratings_frame has columns user_id, movie_id
# we construct a new matrix with users as rows and all ratings as columns
# we use missing for non-present movies

users = unique(ratings_frame[:, :userId])
movies = unique(ratings_frame[:, :movieId])

@assert users == eachindex(users)
movies_revert = zeros(Int, maximum(movies))
for (idx, m) in enumerate(movies)
    movies_revert[m] = idx
end
movies_indices = [movies_revert[idx] for idx in ratings_frame[:, :movieId]]

const rating_matrix = sparse(
    ratings_frame[:, :userId],
    movies_indices,
    ratings_frame[:, :rating],
    length(users),
    length(movies),
)

missing_rate = 0.05

Random.seed!(42)

const missing_ratings = Tuple{Int,Int}[]
const present_ratings = Tuple{Int,Int}[]
let
    (I, J, V) = SparseArrays.findnz(rating_matrix)
    for idx in eachindex(I)
        if V[idx] > 0
            if rand() <= missing_rate
                push!(missing_ratings, (I[idx], J[idx]))
            else
                push!(present_ratings, (I[idx], J[idx]))
            end
        end
    end
end

function f(X)
    # note: we iterate over the rating_matrix indices,
    # since it is sparse unlike X
    r = 0.0
    for (i, j) in present_ratings
        r += 0.5 * (X[i, j] - rating_matrix[i, j])^2
    end
    return r
end

function grad!(storage, X)
    storage .= 0
    for (i, j) in present_ratings
        storage[i, j] = X[i, j] - rating_matrix[i, j]
    end
    return nothing
end

function test_loss(X)
    r = 0.0
    for (i, j) in missing_ratings
        r += 0.5 * (X[i, j] - rating_matrix[i, j])^2
    end
    return r
end

function project_nuclear_norm_ball(X; radius=1.0)
    U, sing_val, Vt = svd(X)
    if (sum(sing_val) <= radius)
        return X, -norm_estimation * U[:, 1] * Vt[:, 1]'
    end
    sing_val = FrankWolfe.projection_simplex_sort(sing_val, s=radius)
    return U * Diagonal(sing_val) * Vt', -norm_estimation * U[:, 1] * Vt[:, 1]'
end

#norm_estimation = 400 * Arpack.svds(rating_matrix, nsv=1, ritzvec=false)[1].S[1]
norm_estimation = 10 * Arpack.svds(rating_matrix, nsv=1, ritzvec=false)[1].S[1]

const lmo = FrankWolfe.NuclearNormLMO(norm_estimation)
const x0 = FrankWolfe.compute_extreme_point(lmo, zero(rating_matrix))
const k = 100

# benchmark the oracles
FrankWolfe.benchmark_oracles(
    f,
    (str, x) -> grad!(str, x),
    () -> randn(size(rating_matrix)),
    lmo;
    k=100,
)

gradient = spzeros(size(x0)...)
gradient_aux = spzeros(size(x0)...)

# pushes to the trajectory the first 5 elements of the trajectory and the test value at the current iterate
function build_callback(trajectory_arr)
    return function callback(state)
        return push!(trajectory_arr, (Tuple(state)[1:5]..., test_loss(state.x)))
    end
end

#Estimate the smoothness constant.
num_pairs = 1000
L_estimate = -Inf
for i in 1:num_pairs
    global L_estimate
    # computing random rank-one matrices
    u1 = rand(size(x0, 1))
    u1 ./= sum(u1)
    u1 .*= norm_estimation
    v1 = rand(size(x0, 2))
    v1 ./= sum(v1)
    x = FrankWolfe.RankOneMatrix(u1, v1)
    u2 = rand(size(x0, 1))
    u2 ./= sum(u2)
    u2 .*= norm_estimation
    v2 = rand(size(x0, 2))
    v2 ./= sum(v2)
    y = FrankWolfe.RankOneMatrix(u2, v2)
    grad!(gradient, x)
    grad!(gradient_aux, y)
    new_L = norm(gradient - gradient_aux) / norm(x - y)
    if new_L > L_estimate
        L_estimate = new_L
    end
end

# PGD steps

xgd = Matrix(x0)
function_values = Float64[]
timing_values = Float64[]
function_test_values = Float64[]

time_start = time_ns()
for _ in 1:k
    f_val = f(xgd)
    push!(function_values, f_val)
    push!(function_test_values, test_loss(xgd))
    push!(timing_values, (time_ns() - time_start) / 1e9)
    @info f_val
    grad!(gradient, xgd)
    xgd_new, vertex = project_nuclear_norm_ball(xgd - gradient / L_estimate, radius=norm_estimation)
    gamma, _ = FrankWolfe.backtrackingLS(f, gradient, xgd, xgd - xgd_new, 1.0)
    @. xgd -= gamma * (xgd - xgd_new)
end

trajectory_arr_fw = Vector{Tuple{Int64,Float64,Float64,Float64,Float64,Float64}}()
callback = build_callback(trajectory_arr_fw)
xfin, _, _, _, traj_data = FrankWolfe.frank_wolfe(
    f,
    grad!,
    lmo,
    x0;
    epsilon=1e-9,
    max_iteration=10 * k,
    print_iter=k / 10,
    verbose=true,
    linesearch_tol=1e-8,
    line_search=FrankWolfe.Adaptive(),
    emphasis=FrankWolfe.memory,
    gradient=gradient,
    callback=callback,
)

trajectory_arr_lazy = Vector{Tuple{Int64,Float64,Float64,Float64,Float64,Float64}}()
callback = build_callback(trajectory_arr_lazy)
xlazy, _, _, _, _ = FrankWolfe.lazified_conditional_gradient(
    f,
    grad!,
    lmo,
    x0;
    epsilon=1e-9,
    max_iteration=10 * k,
    print_iter=k / 10,
    verbose=true,
    linesearch_tol=1e-8,
    line_search=FrankWolfe.Adaptive(),
    emphasis=FrankWolfe.memory,
    gradient=gradient,
    callback=callback,
)

trajectory_arr_lazy_ref = Vector{Tuple{Int64,Float64,Float64,Float64,Float64,Float64}}()
callback = build_callback(trajectory_arr_lazy_ref)
xlazy, _, _, _, _ = FrankWolfe.lazified_conditional_gradient(
    f,
    grad!,
    lmo,
    x0;
    epsilon=1e-9,
    max_iteration=50 * k,
    print_iter=k / 10,
    verbose=true,
    linesearch_tol=1e-8,
    line_search=FrankWolfe.Adaptive(),
    emphasis=FrankWolfe.memory,
    gradient=gradient,
    callback=callback,
)

@info "Gdescent test loss: $(test_loss(xgd))"
@info "FW test loss: $(test_loss(xfin))"
@info "LCG test loss: $(test_loss(xlazy))"

fw_test_values = getindex.(trajectory_arr_fw, 6)
lazy_test_values = getindex.(trajectory_arr_lazy, 6)

open(joinpath(@__DIR__, "movielens_result.json"), "w") do f
    data = JSON.json((
        svals_gd=svdvals(xgd),
        svals_fw=svdvals(xfin),
        svals_lcg=svdvals(xlazy),
        fw_test_values=fw_test_values,
        lazy_test_values=lazy_test_values,
        trajectory_arr_fw=trajectory_arr_fw,
        trajectory_arr_lazy=trajectory_arr_lazy,
        function_values_gd=function_values,
        function_values_test_gd=function_test_values,
        timing_values_gd=timing_values,
        trajectory_arr_lazy_ref=trajectory_arr_lazy_ref,
    ))
    return write(f, data)
end

#Plot results w.r.t. iteration count
gr()
pit = plot(
    getindex.(trajectory_arr_fw, 1),
    getindex.(trajectory_arr_fw, 2),
    label="FW",
    xlabel="iterations",
    ylabel="Objective function",
    yaxis=:log,
    yguidefontsize=8,
    xguidefontsize=8,
    legendfontsize=8,
    legend=:bottomleft,
)
plot!(getindex.(trajectory_arr_lazy, 1), getindex.(trajectory_arr_lazy, 2), label="LCG")
plot!(eachindex(function_values), function_values, yaxis=:log, label="GD")
plot!(eachindex(function_test_values), function_test_values, label="GD_test")
plot!(getindex.(trajectory_arr_fw, 1), getindex.(trajectory_arr_fw, 6), label="FW_T")
plot!(getindex.(trajectory_arr_lazy, 1), getindex.(trajectory_arr_lazy, 6), label="LCG_T")
savefig(pit, "objective_func_vs_iteration.pdf")

#Plot results w.r.t. time
pit = plot(
    getindex.(trajectory_arr_fw, 5),
    getindex.(trajectory_arr_fw, 2),
    label="FW",
    ylabel="Objective function",
    yaxis=:log,
    xlabel="time (s)",
    yguidefontsize=8,
    xguidefontsize=8,
    legendfontsize=8,
    legend=:bottomleft,
)

plot!(getindex.(trajectory_arr_lazy, 5), getindex.(trajectory_arr_lazy, 2), label="LCG")
plot!(getindex.(trajectory_arr_lazy, 5), getindex.(trajectory_arr_lazy, 6), label="LCG_T")
plot!(getindex.(trajectory_arr_fw, 5), getindex.(trajectory_arr_fw, 6), label="FW_T")

plot!(timing_values, function_values, label="GD", yaxis=:log)
plot!(timing_values, function_test_values, label="GD_test")

savefig(pit, "objective_func_vs_time.pdf")


text:
using JuliaSetPkg
using Base.Test

x=Array(linspace(-1,1,200))
y=Array(linspace(-1,1,200))
R = z -> (z^2 + -0.4+0.6*im )

#test size
A=juliaSet(R,x,y)
nx,ny=size(A)
@test length(x)==nx
@test length(y)==ny

#test fixed point
R = z -> z
A=juliaSet(R,x,y,100)
@test A == zeros(200,200)+101

ext: jl

max_stars_repo_name: anjandn/JuliaSetPkg

max_issues_repo_name: anjandn/JuliaSetPkg

max_forks_repo_name: anjandn/JuliaSetPkg


text:
export  generateCareerProgression

"""
```
generateCareerProgression(
    mpSim::MPsim,
    idList::String... )
```
This function generates a report of the careers that the personnel members with IDs in `idList` had in the manpower simulation `mpSim`.

This function returns a `Dict{String, Tuple}` where the keys are the IDs that are in the simulation. The value associated with each IDs is a `Tuple{Float64, Dict}`, where the number is the age of the person at recruitment, and the dictionary consists of time stamps (`Float64` key) and a transition/source node/target node `NTuple{3, String}`.
"""
function generateCareerProgression( mpSim::MPsim,
    idList::String... )::Dict{String,Tuple{Float64,DataFrame}}

    result = Dict{String,Tuple{Float64,DataFrame}}()

    if now( mpSim ) == 0
        @warn "Simulation hasn't started yet, can't make report."
        return result
    end  # if now( mpSim ) == 0

    queryCmd = string( "SELECT `", mpSim.idKey, "`, ageAtRecruitment FROM `", mpSim.persDBname, "` WHERE `",
        mpSim.idKey, "` IN ('", join( idList, "', '" ), "')" )
    recruitmentAges = DataFrame( DBInterface.execute( mpSim.simDB, queryCmd ) )

    if isempty( recruitmentAges )
        return result
    end  # if isempty( recruitmentAges )

    idSymbol = Symbol( mpSim.idKey )
    queryCmd = string( "SELECT * FROM `", mpSim.transDBname, "` WHERE",
        "\n    `", mpSim.idKey, "` IN ('",
        join( recruitmentAges[:, idSymbol], "', '" ), "')",
        "\n    ORDER BY `", mpSim.idKey, "`, timeIndex" )
    careerPaths = DataFrame( DBInterface.execute( mpSim.simDB, queryCmd ) )

    for ii in 1:size( recruitmentAges, 1 )
        id = recruitmentAges[ii, idSymbol]
        careerPath = careerPaths[careerPaths[:, idSymbol] .== id, :]
        result[id] = (recruitmentAges[ii, :ageAtRecruitment], careerPath)
    end  # for ii in 1:size( recruitmentAges, 1 )

    return result

end  # generateCareerProgression( mpSim, idList )


text:
using JSON

mutable struct Body
    name::String
    mass::Float64
    x::Float64
    y::Float64
    z::Float64
    vx::Float64
    vy::Float64
    vz::Float64
end

function read_bodies(filename)
    bodies = Body[]
    open(filename, "r") do handle
        json = JSON.parse(handle)
        unitmass = if haskey(json, "unitmass")
            json["unitmass"]
        else
            1.0
        end
        for (key, value) in json
            if key != "unitmass"
                body = Body(key, unitmass * value["mass"],
                            value["position"]..., value["velocity"]...)
                push!(bodies, body)
            end
        end
    end
    bodies
end

function energy(bodies)
    E = 0.0
    for i in 1:length(bodies)
        E += kinetic_energy(bodies[i])
        for j in (i+1):length(bodies)
            E += potential_energy(bodies[i], bodies[j])
        end
    end
    E
end

kinetic_energy(body) = 0.5 * body.mass * (body.vx^2 + body.vy^2 + body.vz^2)

function potential_energy(a, b)
    dx = a.x - b.x
    dy = a.y - b.y
    dz = a.z - b.z
    -a.mass * b.mass / sqrt(dx^2 + dy^2 + dz^2)
end

function advance!(bodies, dt)
    for i in 1:length(bodies)
        a = bodies[i]
        for j in (i+1):length(bodies)
            b = bodies[j]

            dx = a.x - b.x
            dy = a.y - b.y
            dz = a.z - b.z
            d2 = dx^2 + dy^2 + dz^2
            ds = sqrt(d2)
            magnitude = dt / (d2*d2)
            magnitude *= ds

            a.vx -= b.mass * magnitude * dx
            a.vy -= b.mass * magnitude * dy
            a.vz -= b.mass * magnitude * dz

            b.vx += a.mass * magnitude * dx
            b.vy += a.mass * magnitude * dy
            b.vz += a.mass * magnitude * dz
        end
    end

    for b in bodies
        b.x += dt * b.vx
        b.y += dt * b.vy
        b.z += dt * b.vz
    end
end

function main(filename, time_steps)
    bodies = read_bodies(filename)

    for body in bodies
        println(body)
    end

    @printf "Initial Energy: %.9f\n" energy(bodies)
    duration = @elapsed for _ in 1:time_steps
        advance!(bodies, 0.01)
    end
    @printf "Final Energy:   %.9f\n" energy(bodies)
    @printf "Elapsed: %fs\n\n" duration
end

if length(ARGS) < 2
    error("usage: julia nbody.jl <json-file> <time-steps>")
end

const filename = ARGS[1]
const time_steps = parse(Int, ARGS[2])

main(filename, time_steps)

max_stars_repo_name: ELIFE-ASU/meet-julia

max_issues_repo_name: ELIFE-ASU/meet-julia

max_forks_repo_name: ELIFE-ASU/meet-julia


text:
using Distributed

"""
Run model for a single zone.

    run_model(zone, ts_func, callback_func)

Parameters
----------
zone : Zone AgComponent
ts_func : Function, defining actions for a time step.
        Must accept a Manager, zone, int, and Date/DateTime
        `ts_func(manager, zone, idx, dt)`

        All operations must be in-place as changes will not propagate.
pre : Function, defines additional actions for `zone` that occur
        at end of time step `dt_i`:

        `callback_func(zone, dt_i)`
post : Function, defines additional actions for `zone` that occur
                at end of time step `dt_i`:

                `callback_func(zone, dt_i)`

Returns
----------
DataFrame : results
"""
function run_model(zone::FarmZone, ts_func::Function; pre::Union{Function, Nothing}=nothing, post::Union{Function, Nothing}=nothing)

    time_sequence = zone.climate.time_steps

    @inbounds for (idx, dt_i) in enumerate(time_sequence)
        if !isnothing(pre)
            pre(zone, dt_i)
        end

        ts_func(zone.manager, zone, idx, dt_i)

        if !isnothing(post)
            post(zone, dt_i)
        end
    end

    return collect_results(zone)
end

"""Run timestep for all zones within a basin."""
function run_timestep!(basin::Basin, ts_func; pre::Union{Function, Nothing}=nothing, post::Union{Function, Nothing}=nothing)

    idx, dt_i = basin.current_ts

    for z in basin.zones
        if !isnothing(pre)
            pre(z, dt_i)
        end

        ts_func(z.manager, z, idx, dt_i)

        if !isnothing(post)
            post(z, dt_i)
        end
    end

    advance_timestep!(basin)
end

"""Run scenario for an entire basin."""
function run_model(basin::Basin, ts_func::Function; pre::Union{Function, Nothing}=nothing, post::Union{Function, Nothing}=nothing)
    for i in basin.climate.time_steps
        run_timestep!(basin, ts_func; pre=pre, post=post)
    end

    return collect_results(basin)
end

function advance_timestep!(basin)::Nothing
    idx::Int64, _ = basin.current_ts
    next_idx::Int64 = idx+1
    dt = nothing

    try
        dt = basin.climate.time_steps[next_idx]
    catch e
        if e isa BoundsError
            return
        end

        throw(e)
    end

    if isnothing(dt)
        throw(ArgumentError("datetime cannot be nothing"))
    end

    basin.current_ts = (idx=next_idx, dt=dt)

    return
end

"""
Run zone level scenarios for a given sample set.

Returns Dict with keys:
    "scenario_id/zone_results" = zone level results
    "scenario_id/field_results" = field level results
"""
function run_scenarios!(samples::DataFrame, zone::FarmZone, ts_func::Function;
                        pre::Union{Function, Nothing}=nothing, post::Union{Function, Nothing}=nothing)::Dict
    results = @sync @distributed (hcat) for row_id in 1:nrow(samples)
        tmp_z = deepcopy(zone)
        update_model!(tmp_z, samples[row_id, :])
        res = run_model(tmp_z, ts_func; pre=pre, post=post)

        # Return pair of scenario_id and results
        string(row_id), res
    end

    # Prep results into expected form
    transformed = Dict()
    for (idx, res) in results
        transformed["$(idx)/zone_results"] = res.zone_results
        transformed["$(idx)/field_results"] = res.field_results
    end

    return transformed
end

"""
Run basin level scenarios for a given sample set.

Returns Dict with keys:
    "scenario_id/zone_results" = zone level results for each zone in basin
    "scenario_id/field_results" = field level results for each zone in basin
"""
function run_scenarios!(samples::DataFrame, basin::Basin, ts_func::Function;
                        pre::Union{Function, Nothing}=nothing, post::Union{Function, Nothing}=nothing)::Dict
    results = @sync @distributed (hcat) for row_id in 1:nrow(samples)
        tmp_b = deepcopy(basin)
        update_model!(tmp_b, samples[row_id, :])
        res = run_model(tmp_b, ts_func; pre=pre, post=post)

        # Return pair of scenario_id and results
        string(row_id), res
    end

    # # Prep results into expected form
    transformed = Dict()
    for (idx, res) in results
        transformed["$(idx)/zone_results"] = Dict(
            k => v.zone_results for (k, v) in res
        )
        transformed["$(idx)/field_results"] = Dict(
            k => v.field_results for (k, v) in res
        )
    end

    return transformed
end

# todo
struct BasinResults
end

struct ZoneResults
end

# """
# Setup a model with a specific timestep function.
# """
# function setup_model(catchment::Union{FarmZone, Basin}, ts_func::Function; pre::Union{Function, Nothing}=nothing, post::Union{Function, Nothing}=nothing)
#     function run_model(catchment::Union{FarmZone, Basin})
#         time_sequence = catchment.climate.time_steps
#         @inbounds for (idx, dt_i) in enumerate(time_sequence)
#             if !isnothing(pre)
#                 pre(catchment, dt_i)
#             end

#             ts_func(catchment.manager, catchment, idx, dt_i)

#             if !isnothing(post)
#                 post(catchment, dt_i)
#             end
#         end

#         return collect_results(catchment)
#     end

#     return run_model
# end

max_stars_repo_name: ConnectedSystems/agtor-jl

max_issues_repo_name: ConnectedSystems/agtor-jl

max_forks_repo_name: ConnectedSystems/agtor-jl


text:
#############################################################################
# extreme_point.jl
# Compute the of the base polyhedron for an given ordering, i.e., the elements
# of the result is calculated in the order of the given ordering.
#############################################################################

export extreme_point

function extreme_point(F::Function, w::AbstractArray,
                       V = collect(1:length(w))::AbstractArray)
  @assert F([])[1] == 0 # f should be 0 at the empty set
  w = Int.(w)
  n = length(w)
  @assert length(V) == n
  x = zeros(n)
  for ii = 1:n
    x[V[w[ii]]] = F(V[w[1: ii]])[1] - F(V[w[1: ii-1]])[1]
  end
  return x
end

extreme_point(F::SubmodFunc, w::AbstractArray) = extreme_point(x -> evaluate(F, x), w)


text:
function bubblesort(arr)
    n = length(arr)

    for i=1:n, j=1:n-i
        if (arr[j] > arr[j+1])
            arr[j], arr[j+1] = arr[j+1], arr[j]
        end
    end
end

v = [144 89 4 9 95 12 86 25]
print("unbubble sorted array = ",v)

bubblesort(v)

print("\nbubble sorted array = ",v)

max_stars_repo_name: landrumb/julia

max_issues_repo_name: landrumb/julia

max_forks_repo_name: landrumb/julia


text:
function get_column_dimension(obj::AnyMatrix)
    return jcall(obj, "getColumnDimension", jint, ())
end

function get_row_dimension(obj::AnyMatrix)
    return jcall(obj, "getRowDimension", jint, ())
end

function is_square(obj::AnyMatrix)
    return jcall(obj, "isSquare", jboolean, ())
end


text:
include("splinesParallel.jl")

b1 = plarBezier(S1, [[0,1,0],[0,1,5]])
b2 = plarBezier(S1, [[0,0,0],[0,0,5]])
b3 = plarBezier(S1, [[1,0,0],[2,-1,2.5],[1,0,5]])
b4 = plarBezier(S1, [[1,1,0],[1,1,5]])
b5 = plarBezier(S1, [[0,1,0],[0,1,5]])
controls = [b1, b2, b3, b4, b5]
knots = [0,1,2,3,4,5,6,7]           # periodic B-spline
knots = [0,0,0,1,2,3,3,3]           # non-periodic B-spline
ptbspline = pt_bspline(S2, 2, knots, controls)
point = [0.7, 0.2]
result = ptbspline(point)


text:
using Distributed, Flux

addprocs(2)
@everywhere begin
  using DiffEqSensitivity, OrdinaryDiffEq, Test

  pa = [1.0]
  u0 = [3.0]
end

function model4()
  prob = ODEProblem((u, p, t) -> 1.01u .* p, u0, (0.0, 1.0), pa)

  function prob_func(prob, i, repeat)
    remake(prob, u0 = 0.5 .+ i/100 .* prob.u0)
  end

  ensemble_prob = EnsembleProblem(prob, prob_func = prob_func)
  sim = solve(ensemble_prob, Tsit5(), EnsembleDistributed(), saveat = 0.1, trajectories = 100)
end

# loss function
loss() = sum(abs2,1.0.-Array(model4()))

data = Iterators.repeated((), 10)

cb = function () # callback function to observe training
  @show loss()
end

pa = [1.0]
u0 = [3.0]
opt = Flux.ADAM(0.1)
println("Starting to train")
l1 = loss()
Flux.@epochs 10 Flux.train!(loss, Flux.params([pa,u0]), data, opt; cb = cb)
l2 = loss()
@test 10l2 < l1


text:
using Memoize

function testlcs()
	s1::String = "Please write comments if you find anything incorrect"
	s2::String = "or you want to share more information about the topic discussed above"
	@time a = lcs(s1, s2)
	@show a
end

function lcs(s1::String, s2::String)
	@memoize function rec(n1::Int, n2::Int)
		(n1 <= 0 || n2 <= 0) && return 0
		return max(	rec(n1 - 1, n2),
					rec(n1, n2 - 1),
					rec(n1 - 1, n2 - 1) +
					(s1[n1] == s2[n2] ? 1 : 0))
	end
	return rec(length(s1), length(s2))
end


text:
# Based on the introductory discussion on
# https://www.maa.org/press/periodicals/loci/joma/the-sir-model-for-spread-of-disease-the-differential-equation-model

## Uncomment the following lines if used as a standalone
# using DifferentialEquations
# using Plots

function SIR(du, u, p, t)
    s, i, r = u
    b, k = p
    du[1] = -b * s * i
    du[2] = b * s * i - k * i
    du[3] = k * i
end

# TODO: make the following a doctest example
    # #initian conditions
    # p = [0.5/7900000.0,0.33]
    # u0 = [7900000.0,10.0,0.0]
    # tspan = (0.0,140.0)

    # #solve
    # sir = ODEProblem(SIR,u0,tspan,p)
    # sol = solve(sir)

    # #plot
    # plot(sol)

max_stars_repo_name: ashwani-rathee/Julia

max_issues_repo_name: ashwani-rathee/Julia

max_forks_repo_name: ashwani-rathee/Julia


text:
function xchart(x)
  2x
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""PersistentVolumeClaimCondition contails details about state of pvc

    IoK8sApiCoreV1PersistentVolumeClaimCondition(;
        lastProbeTime=nothing,
        lastTransitionTime=nothing,
        message=nothing,
        reason=nothing,
        status=nothing,
        type=nothing,
    )

    - lastProbeTime::IoK8sApimachineryPkgApisMetaV1Time : Last time we probed the condition.
    - lastTransitionTime::IoK8sApimachineryPkgApisMetaV1Time : Last time the condition transitioned from one status to another.
    - message::String : Human-readable message indicating details about last transition.
    - reason::String : Unique, this should be a short, machine understandable string that gives the reason for condition's last transition. If it reports \"ResizeStarted\" that means the underlying persistent volume is being resized.
    - status::String
    - type::String
"""
mutable struct IoK8sApiCoreV1PersistentVolumeClaimCondition <: SwaggerModel
    lastProbeTime::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1Time } # spec name: lastProbeTime
    lastTransitionTime::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1Time } # spec name: lastTransitionTime
    message::Any # spec type: Union{ Nothing, String } # spec name: message
    reason::Any # spec type: Union{ Nothing, String } # spec name: reason
    status::Any # spec type: Union{ Nothing, String } # spec name: status
    type::Any # spec type: Union{ Nothing, String } # spec name: type

    function IoK8sApiCoreV1PersistentVolumeClaimCondition(;lastProbeTime=nothing, lastTransitionTime=nothing, message=nothing, reason=nothing, status=nothing, type=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1PersistentVolumeClaimCondition, Symbol("lastProbeTime"), lastProbeTime)
        setfield!(o, Symbol("lastProbeTime"), lastProbeTime)
        validate_property(IoK8sApiCoreV1PersistentVolumeClaimCondition, Symbol("lastTransitionTime"), lastTransitionTime)
        setfield!(o, Symbol("lastTransitionTime"), lastTransitionTime)
        validate_property(IoK8sApiCoreV1PersistentVolumeClaimCondition, Symbol("message"), message)
        setfield!(o, Symbol("message"), message)
        validate_property(IoK8sApiCoreV1PersistentVolumeClaimCondition, Symbol("reason"), reason)
        setfield!(o, Symbol("reason"), reason)
        validate_property(IoK8sApiCoreV1PersistentVolumeClaimCondition, Symbol("status"), status)
        setfield!(o, Symbol("status"), status)
        validate_property(IoK8sApiCoreV1PersistentVolumeClaimCondition, Symbol("type"), type)
        setfield!(o, Symbol("type"), type)
        o
    end
end # type IoK8sApiCoreV1PersistentVolumeClaimCondition

const _property_map_IoK8sApiCoreV1PersistentVolumeClaimCondition = Dict{Symbol,Symbol}(Symbol("lastProbeTime")=>Symbol("lastProbeTime"), Symbol("lastTransitionTime")=>Symbol("lastTransitionTime"), Symbol("message")=>Symbol("message"), Symbol("reason")=>Symbol("reason"), Symbol("status")=>Symbol("status"), Symbol("type")=>Symbol("type"))
const _property_types_IoK8sApiCoreV1PersistentVolumeClaimCondition = Dict{Symbol,String}(Symbol("lastProbeTime")=>"IoK8sApimachineryPkgApisMetaV1Time", Symbol("lastTransitionTime")=>"IoK8sApimachineryPkgApisMetaV1Time", Symbol("message")=>"String", Symbol("reason")=>"String", Symbol("status")=>"String", Symbol("type")=>"String")
Base.propertynames(::Type{ IoK8sApiCoreV1PersistentVolumeClaimCondition }) = collect(keys(_property_map_IoK8sApiCoreV1PersistentVolumeClaimCondition))
Swagger.property_type(::Type{ IoK8sApiCoreV1PersistentVolumeClaimCondition }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1PersistentVolumeClaimCondition[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1PersistentVolumeClaimCondition }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1PersistentVolumeClaimCondition[property_name]

function check_required(o::IoK8sApiCoreV1PersistentVolumeClaimCondition)
    (getproperty(o, Symbol("status")) === nothing) && (return false)
    (getproperty(o, Symbol("type")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1PersistentVolumeClaimCondition }, name::Symbol, val)
end


text:
using MPI
using CLIMA
using CLIMA.Mesh.Topologies
using CLIMA.Mesh.Grids
using CLIMA.DGBalanceLawDiscretizations
using CLIMA.DGBalanceLawDiscretizations.NumericalFluxes
using CLIMA.MPIStateArrays
using CLIMA.ODESolvers
using CLIMA.GenericCallbacks
using LinearAlgebra
using StaticArrays
using Logging, Printf, Dates
using CLIMA.VTK
using CLIMA.LinearSolvers
using CLIMA.GeneralizedConjugateResidualSolver

const γ_exact = 7 // 5 # FIXME: Remove this for some moist thermo approach

using CLIMA.MoistThermodynamics

using CLIMAParameters
using CLIMAParameters.Planet: R_d, cp_d, grav, cv_d, MSLP
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

const _nstate = 5
const _δρ, _ρu, _ρv, _ρw, _δρe = 1:_nstate
const _ρu⃗ = SVector(_ρu, _ρv, _ρw)
const statenames = ("δρ", "ρu", "ρv", "ρw", "δρe")

const _nauxstate = 6
const _a_ρ0, _a_ρe0, _a_ϕ, _a_ϕ_x, _a_ϕ_y, _a_ϕ_z = 1:_nauxstate
function auxiliary_state_initialization!(aux, x, y, z) #JK, dx, dy, dz)
    FT = eltype(aux)
    @inbounds begin
        ρ0, ρe0 = reference_ρ_ρe(x, y, z)
        _grav::FT = grav(param_set)
        aux[_a_ρ0] = ρ0
        aux[_a_ρe0] = ρe0
        aux[_a_ϕ] = y * _grav
        aux[_a_ϕ_x] = 0
        aux[_a_ϕ_y] = _grav
        aux[_a_ϕ_z] = 0
    end
end

@inline function pressure(Q, aux)
    @inbounds begin
        γ::eltype(Q) = γ_exact # FIXME: Remove this for some moist thermo approach
        δρ, δρe = Q[_δρ], Q[_δρe]
        ρu⃗ = SVector(Q[_ρu], Q[_ρv], Q[_ρw])
        ρ0, ρe0, ϕ = aux[_a_ρ0], aux[_a_ρe0], aux[_a_ϕ]
        ρ = ρ0 + δρ
        ρe = ρe0 + δρe
        ρinv = 1 / ρ
        (γ - 1) * (ρe - ρinv * (ρu⃗' * ρu⃗) / 2 - ϕ * ρ)
    end
end

@inline function wavespeed(n, Q, aux, t)
    γ::eltype(Q) = γ_exact # FIXME: Remove this for some moist thermo approach
    n⃗ = SVector(n)
    @inbounds begin
        P = pressure(Q, aux)
        δρ, δρe = Q[_δρ], Q[_δρe]
        ρu⃗ = SVector(Q[_ρu], Q[_ρv], Q[_ρw])
        ρ0, ρe0 = aux[_a_ρ0], aux[_a_ρe0]
        ρ = ρ0 + δρ
        ρe = ρe0 + δρe
        ρinv = 1 / ρ
        u⃗ = ρinv * ρu⃗
        abs(n⃗' * u⃗) + sqrt(ρinv * γ * P)
    end
end

@inline function euler_flux!(F, Q, _, aux, t)
    P = pressure(Q, aux)
    @inbounds begin
        δρ, δρe = Q[_δρ], Q[_δρe]
        ρu⃗ = SVector(Q[_ρu], Q[_ρv], Q[_ρw])
        ρ0, ρe0 = aux[_a_ρ0], aux[_a_ρe0]
        ρ = ρ0 + δρ
        ρe = ρe0 + δρe
        ρinv = 1 / ρ
        u⃗ = ρinv * ρu⃗

        F[:, _δρ] = ρu⃗
        F[:, _ρu⃗] = u⃗ * ρu⃗' + P * I
        F[:, _δρe] = u⃗ * (ρe + P)
    end
end

@inline function bcstate!(QP, QVP, auxP, nM, QM, QMP, auxM, bctype, t)
    if bctype == 1
        nofluxbc!(QP, QVP, auxP, nM, QM, QMP, auxM, bctype, t)
    else
        error("Unsupported boundary condition")
    end
end

@inline function nofluxbc!(QP, QVP, auxP, nM, QM, QMP, auxM, bctype, t)
    @inbounds begin
        ρu⃗M = SVector(QM[_ρu], QM[_ρv], QM[_ρw])
        n⃗M = SVector(nM)

        # No flux boundary conditions
        # No shear on walls (free-slip condition)
        ρunM = n⃗M' * ρu⃗M
        QP[_ρu⃗[1]] = ρu⃗M[1] - 2ρunM * n⃗M[1]
        QP[_ρu⃗[2]] = ρu⃗M[2] - 2ρunM * n⃗M[2]
        QP[_ρu⃗[3]] = ρu⃗M[3] - 2ρunM * n⃗M[3]
        QP[_δρ] = QM[_δρ]
        QP[_δρe] = QM[_δρe]
        nothing
    end
end

@inline function source!(S, Q, diffusive, aux, t)
    # Initialise the final block source term
    S .= -zero(eltype(Q))

    source_geopotential!(S, Q, aux, t)
end

@inline function source_geopotential!(S, Q, aux, t)
    @inbounds begin
        δρ = Q[_δρ]
        ρ0 = aux[_a_ρ0]
        ρ = ρ0 + δρ
        ∇ϕ = SVector(aux[_a_ϕ_x], aux[_a_ϕ_y], aux[_a_ϕ_z])
        S[_ρu⃗[1]] -= ρ * ∇ϕ[1]
        S[_ρu⃗[2]] -= ρ * ∇ϕ[2]
        S[_ρu⃗[3]] -= ρ * ∇ϕ[3]
    end
end

function reference_ρ_ρe(x, y, z)
    FT = eltype(x)
    _R_d::FT = R_d(param_set)
    _cp_d::FT = cp_d(param_set)
    _cv_d::FT = cv_d(param_set)
    _MSLP::FT = MSLP(param_set)
    _grav::FT = grav(param_set)
    # perturbation parameters for rising bubble

    θ0::FT = 303
    π_exner = 1 - _grav / (_cp_d * θ0) * y # exner pressure
    ρ = _MSLP / (_R_d * θ0) * (π_exner)^(_cv_d / _R_d) # density

    P = _MSLP * (_R_d * (ρ * θ0) / _MSLP)^(_cp_d / _cv_d) # pressure (absolute)
    T = P / (ρ * _R_d) # temperature
    u⃗ = SVector(-zero(FT), -zero(FT), -zero(FT))
    # energy definitions
    e_kin = u⃗' * u⃗ / 2
    e_pot = _grav * y
    e_int = _cv_d * T
    ρe = ρ * (e_int + e_kin + e_pot)
    ρ, ρe
end

# Initial Condition
function rising_bubble!(dim, Q, t, x, y, z, aux)
    FT = eltype(Q)
    _R_d::FT = R_d(param_set)
    _cp_d::FT = cp_d(param_set)
    _cv_d::FT = cv_d(param_set)
    _MSLP::FT = MSLP(param_set)
    _grav::FT = grav(param_set)
    # perturbation parameters for rising bubble

    r⃗ = SVector(x, y, z)
    r⃗_center = SVector(500, 260, 500)
    distance = norm(@view (r⃗ - r⃗_center)[1:dim])

    θ0::FT = 303
    θ_c::FT = 1 // 2
    Δθ::FT = -zero(FT)
    a::FT = 50
    s::FT = 100
    if distance <= a
        Δθ = θ_c
    elseif distance > a
        Δθ = θ_c * exp(-(distance - a)^2 / s^2)
    end
    θ = θ0 + Δθ # potential temperature
    π_exner = 1 - _grav / (_cp_d * θ) * y # exner pressure
    ρ = _MSLP / (_R_d * θ) * (π_exner)^(_cv_d / _R_d) # density

    P = _MSLP * (_R_d * (ρ * θ) / _MSLP)^(_cp_d / _cv_d) # pressure (absolute)
    T = P / (ρ * _R_d) # temperature
    u⃗ = SVector(-zero(FT), -zero(FT), -zero(FT))
    # energy definitions
    e_kin = u⃗' * u⃗ / 2
    e_pot = _grav * y
    e_int = _cv_d * T
    ρe = ρ * (e_int + e_kin + e_pot)
    ρu⃗ = ρ * u⃗

    @inbounds ρ0, ρe0 = aux[_a_ρ0], aux[_a_ρe0]
    @inbounds Q[_δρ], Q[_δρe] = ρ - ρ0, ρe - ρe0
    @inbounds Q[_ρu⃗] = ρu⃗
end

# {{{ Linearization
@inline function lin_eulerflux!(F, Q, _, aux, t)
    F .= -zero(eltype(Q))
    @inbounds begin
        FT = eltype(Q)
        γ::FT = γ_exact # FIXME: Remove this for some moist thermo approach

        δρ, δρe = Q[_δρ], Q[_δρe]
        ρu⃗ = SVector(Q[_ρu], Q[_ρv], Q[_ρw])

        ρ0, ρe0 = aux[_a_ρ0], aux[_a_ρe0]
        ϕ = aux[_a_ϕ]

        ρinv0 = 1 / ρ0
        e0 = ρinv0 * ρe0

        P0 = (γ - 1) * (ρe0 - ρ0 * ϕ)
        δP = (γ - 1) * (δρe - δρ * ϕ)

        p0 = ρinv0 * P0

        F[:, _δρ] = ρu⃗
        F[:, _ρu⃗] = δP * I + @SMatrix zeros(3, 3)
        F[:, _δρe] = (e0 + p0) * ρu⃗
    end
end

@inline function wavespeed_linear(n, Q, aux, t)
    @inbounds begin
        FT = eltype(Q)
        γ::FT = γ_exact # FIXME: Remove this for some moist thermo approach

        ρ0, ρe0, ϕ = aux[_a_ρ0], aux[_a_ρe0], aux[_a_ϕ]
        ρinv0 = 1 / ρ0
        P0 = (γ - 1) * (ρe0 - ρ0 * ϕ)

        sqrt(ρinv0 * γ * P0)
    end
end

@inline function lin_source!(S, Q, diffusive, aux, t)
    S .= 0
    lin_source_geopotential!(S, Q, diffusive, aux, t)
end
@inline function lin_source_geopotential!(S, Q, diffusive, aux, t)
    @inbounds begin
        δρ = Q[_δρ]
        ∇ϕ = SVector(aux[_a_ϕ_x], aux[_a_ϕ_y], aux[_a_ϕ_z])
        S[_ρu] -= δρ * ∇ϕ[1]
        S[_ρv] -= δρ * ∇ϕ[2]
        S[_ρw] -= δρ * ∇ϕ[3]
    end
end

@inline function lin_bcstate!(QP, QVP, auxP, nM, QM, QMP, auxM, bctype, t)
    if bctype == 1
        # this is already a linear boundary condition
        nofluxbc!(QP, QVP, auxP, nM, QM, QMP, auxM, bctype, t)
    else
        error("Unsupported boundary condition")
    end
end

# }}}

function run(
    mpicomm,
    ArrayType,
    dim,
    brickrange,
    periodicity,
    N,
    timeend,
    FT,
    dt,
    output_steps,
)

    topl = StackedBrickTopology(mpicomm, brickrange, periodicity = periodicity)

    grid = DiscontinuousSpectralElementGrid(
        topl,
        FloatType = FT,
        DeviceArray = ArrayType,
        polynomialorder = N,
    )

    numflux!(x...) = NumericalFluxes.rusanov!(x..., euler_flux!, wavespeed)
    numbcflux!(x...) = NumericalFluxes.rusanov_boundary_flux!(
        x...,
        euler_flux!,
        bcstate!,
        wavespeed,
    )

    # spacedisc = data needed for evaluating the right-hand side function
    spacedisc = DGBalanceLaw(
        grid = grid,
        length_state_vector = _nstate,
        flux! = euler_flux!,
        numerical_flux! = numflux!,
        numerical_boundary_flux! = numbcflux!,
        auxiliary_state_length = _nauxstate,
        auxiliary_state_initialization! = auxiliary_state_initialization!,
        source! = source!,
    )

    # This is a actual state/function that lives on the grid
    initialcondition(Q, x...) = rising_bubble!(dim, Q, FT(0), x...)
    Q = MPIStateArray(spacedisc, initialcondition)

    # {{{ Lineariztion Setup
    lin_numflux!(x...) = NumericalFluxes.rusanov!(
        x...,
        lin_eulerflux!,
        # (_...)->0) # central
        wavespeed_linear,
    )
    lin_numbcflux!(x...) = NumericalFluxes.rusanov_boundary_flux!(
        x...,
        lin_eulerflux!,
        lin_bcstate!,
        # (_...)->0) # central
        wavespeed_linear,
    )
    lin_spacedisc = DGBalanceLaw(
        grid = grid,
        length_state_vector = _nstate,
        flux! = lin_eulerflux!,
        numerical_flux! = lin_numflux!,
        numerical_boundary_flux! = lin_numbcflux!,
        auxiliary_state_length = _nauxstate,
        auxiliary_state_initialization! = auxiliary_state_initialization!,
        source! = lin_source!,
    )

    # NOTE: In order to get the same results on the CPU and GPU we force ourselves
    # to take the same number of iterations by setting at really high tolerance
    # specifying the number of restarts
    linearsolver = GeneralizedConjugateResidual(3, Q, rtol = 1e-8)

    timestepper = ARK548L2SA2KennedyCarpenter(
        spacedisc,
        lin_spacedisc,
        linearsolver,
        Q;
        dt = dt,
        t0 = 0,
    )
    # }}}

    eng0 = norm(Q)
    @info @sprintf """Starting
    norm(Q₀) = %.16e""" eng0

    # Set up the information callback
    starttime = Ref(now())
    cbinfo = GenericCallbacks.EveryXWallTimeSeconds(10, mpicomm) do (s = false)
        if s
            starttime[] = now()
        else
            energy = norm(Q)
            #globmean = global_mean(Q, _δρ)
            @info @sprintf(
                """Update
                simtime = %.16e
                runtime = %s
                norm(Q) = %.16e""",
                ODESolvers.gettime(timestepper),
                Dates.format(
                    convert(Dates.DateTime, Dates.now() - starttime[]),
                    Dates.dateformat"HH:MM:SS",
                ),
                energy
            )#, globmean)
        end
    end

    step = [0]
    vtkdir = "vtk_RTB_IMEX"
    mkpath(vtkdir)
    cbvtk =
        GenericCallbacks.EveryXSimulationSteps(output_steps) do (init = false)
            outprefix = @sprintf(
                "%s/RTB_%dD_mpirank%04d_step%04d",
                vtkdir,
                dim,
                MPI.Comm_rank(mpicomm),
                step[1]
            )
            @debug "doing VTK output" outprefix
            writevtk(outprefix, Q, spacedisc, statenames)
            pvtuprefix = @sprintf("RTB_%dD_step%04d", dim, step[1])
            prefixes = ntuple(
                i -> @sprintf(
                    "%s/RTB_%dD_mpirank%04d_step%04d",
                    vtkdir,
                    dim,
                    i - 1,
                    step[1]
                ),
                MPI.Comm_size(mpicomm),
            )
            writepvtu(pvtuprefix, prefixes, statenames)
            step[1] += 1
            nothing
        end

    solve!(Q, timestepper; timeend = timeend, callbacks = (cbinfo, cbvtk))

    # Print some end of the simulation information
    engf = norm(Q)
    @info @sprintf("""norm(QF) / norm(Q0) = %.16e""", engf / eng0)
    engf / eng0
end

using Test
let
    CLIMA.init()
    ArrayType = CLIMA.array_type()

    mpicomm = MPI.COMM_WORLD

    polynomialorder = 4
    FT = Float64

    expected_engf_eng0 = Dict()
    expected_engf_eng0[(Float64, 2)] = 1.5850821145834655e+00
    expected_engf_eng0[(Float64, 3)] = 1.4450465596993558e+00

    for dim in (2, 3)
        # Physical domain extents
        domain_start = (0, 0, 0)
        domain_end = (1000, 750, 1000)

        # Stable explicit time step
        Δxyz = MVector(25, 25, 25)
        dt = min(Δxyz...) / soundspeed_air(300.0) / polynomialorder
        dt *= dim == 2 ? 40 : 20

        output_time = 0.5
        output_steps = ceil(output_time / dt)

        #Get Ne from resolution
        Ls = MVector(domain_end .- domain_start)
        ratios = @. (Ls / Δxyz - 1) / polynomialorder
        Ne = ceil.(Int64, ratios)

        brickrange = ntuple(
            d -> range(
                domain_start[d],
                length = Ne[d] + 1,
                stop = domain_end[d],
            ),
            dim,
        )
        periodicity = ntuple(d -> false, dim)

        timeend = dim == 2 ? 10 : 1

        # only for printing
        if dim == 2
            Δxyz[3] = 0
            Ne[3] = 0
        end

        @info @sprintf """ ----------------------------------------------------"""
        @info @sprintf """   ______ _      _____ __  ________                  """
        @info @sprintf """  |  ____| |    |_   _|  ...  |  __  |               """
        @info @sprintf """  | |    | |      | | |   .   | |  | |               """
        @info @sprintf """  | |    | |      | | | |   | | |__| |               """
        @info @sprintf """  | |____| |____ _| |_| |   | | |  | |               """
        @info @sprintf """  | _____|______|_____|_|   |_|_|  |_|               """
        @info @sprintf """                                                     """
        @info @sprintf """ ----------------------------------------------------"""
        @info @sprintf """ Rising Bubble in %dD                                """ dim
        @info @sprintf """   Resolution:                                       """
        @info @sprintf """     (Δx, Δy, Δz)   = (%.2e, %.2e, %.2e)             """ Δxyz...
        @info @sprintf """     (Nex, Ney, Nez) = (%d, %d, %d)                  """ Ne...
        @info @sprintf """     dt         = %.2e                               """ dt
        @info @sprintf """ ----------------------------------------------------"""

        engf_eng0 = run(
            mpicomm,
            ArrayType,
            dim,
            brickrange,
            periodicity,
            polynomialorder,
            timeend,
            FT,
            dt,
            output_steps,
        )

        @test engf_eng0 ≈ expected_engf_eng0[FT, dim]
    end
end
nothing

max_stars_repo_name: yairchn/CLIMA

max_issues_repo_name: yairchn/CLIMA

max_forks_repo_name: yairchn/CLIMA


text:
export inverse_digits
inverse_digits(digits, base = 10) = sum([digits[k] * base^(k - 1) for k = 1:length(digits)])


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Pre-defined highlighters for the HTML backend.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export hl_cell

"""
    hl_cell(i::Number, j::Number, decoration::HTMLDecoration)

Highlight the cell `(i,j)` with the decoration `decoration` (see
`HTMLDecoration`).

    function hl_cell(cells::AbstractVector{NTuple(2,Int)}, decoration::HTMLDecoration)

Highlights all the cells in `cells` with the decoration `decoration` (see
`HTMLDecoration`).

# Remarks

Those functions return a `HTMLHighlighter` to be used with the HTML backend.

"""
hl_cell(i::Number, j::Number, decoration::HTMLDecoration) = HTMLHighlighter(
    f = (data,x,y)->begin
        return (x == i) && (y == j)
    end,
    decoration = decoration
)

hl_cell(cells::AbstractVector{NTuple{2,Int}}, decoration::HTMLDecoration) = HTMLHighlighter(
    f = (data,x,y)->begin
        return (x,y) ∈ cells
    end,
    decoration = decoration
)

"""
    function hl_col(i::Number, decoration::HTMLDecoration)

Highlight the entire column `i` with the decoration `decoration`.

    function hl_col(cols::AbstractVector{Int}, decoration::HTMLDecoration)

Highlights all the columns in `cols` with the decoration `decoration`.

# Remarks

Those functions return a `HTMLHighlighter` to be used with the HTML backend.

"""
hl_col(j::Number, decoration::HTMLDecoration) = HTMLHighlighter(
    f = (data,x,y)->begin
        return y == j
    end,
    decoration = decoration
)

hl_col(cols::AbstractVector{Int}, decoration::HTMLDecoration) = HTMLHighlighter(
    f = (data,x,y)->begin
        return y ∈ cols
    end,
    decoration = decoration
)

"""
    function hl_row(i::Number, decoration::HTMLDecoration)

Highlight the entire row `i` with the decoration `decoration`.

    function hl_row(rows::AbstractVector{Int}, decoration::HTMLDecoration)

Highlights all the rows in `rows` with the decoration `decoration`.

# Remarks

Those functions return a `HTMLHighlighter` to be used with the HTML backend.

"""
hl_row(i::Number, decoration::HTMLDecoration) = HTMLHighlighter(
    f = (data,x,y)->begin
        return x == i
    end,
    decoration = decoration
)

hl_row(rows::AbstractVector{Int}, decoration::HTMLDecoration) = HTMLHighlighter(
    f = (data,x,y)->begin
        return x ∈ rows
    end,
    decoration = decoration
)

"""
    hl_lt(n::Number, decoration::HTMLDecoration)

Highlight all elements that < `n` using the decoration `decoration`.

# Remarks

Those functions return a `HTMLHighlighter` to be used with the text backend.

"""
hl_lt(n::Number, decoration::HTMLDecoration) = HTMLHighlighter(
    f = (data,i,j)->begin
        if applicable(<,data[i,j],n) && data[i,j] < n
            return true
        else
            return false
        end
    end,
    decoration = decoration
)

"""
    hl_leq(n::Number, decoration::HTMLDecoration)

Highlight all elements that ≤ `n` using the decoration `decoration`.

# Remarks

Those functions return a `HTMLHighlighter` to be used with the text backend.

"""
hl_leq(n::Number, decoration::HTMLDecoration) = HTMLHighlighter(
    f = (data,i,j)->begin
        if applicable(≤,data[i,j],n) && data[i,j] ≤ n
            return true
        else
            return false
        end
    end,
    decoration = decoration
)

"""
    hl_gt(n::Number, decoration::HTMLDecoration)

Highlight all elements that > `n` using the decoration `decoration`.

# Remarks

Those functions return a `HTMLHighlighter` to be used with the text backend.

"""
hl_gt(n::Number, decoration::HTMLDecoration) = HTMLHighlighter(
    f = (data,i,j)->begin
        if applicable(>,data[i,j],n) && data[i,j] > n
            return true
        else
            return false
        end
    end,
    decoration = decoration
)

"""
    function hl_geq(n::Number, decoration::HTMLDecoration)

Highlight all elements that ≥ `n` using the decoration `decoration`.

# Remarks

Those functions return a `HTMLHighlighter` to be used with the text backend.

"""
hl_geq(n::Number, decoration::HTMLDecoration) = HTMLHighlighter(
    f = (data,i,j)->begin
        if applicable(≥,data[i,j],n) && data[i,j] ≥ n
            return true
        else
            return false
        end
    end,
    decoration = decoration
)

"""
    function hl_value(v::Any, decoration::HTMLDecoration)

Highlight all the values that matches `data[i,j] == v` using the decoration
`decoration`.

# Remarks

Those functions return a `HTMLHighlighter` to be used with the text backend.

"""
hl_value(v, decoration::HTMLDecoration) = HTMLHighlighter(
    f = (data,i,j)->data[i,j] == v,
    decoration = decoration
)


text:
using AffineArithmetic
using ValidatedNumerics

X = 1..3
X_aff = Affine(X)

Affine(0.0)

range(X_aff)

X_aff - X_aff

X_aff^2

range(X_aff^2)

X^2 - (X+X)
Y = X_aff^2 - (X_aff + X_aff)
range(Y)

o = Affine(1.0) # constant 1
X_aff - o
(X_aff - o)^2
(X_aff - o)^2 - o

X_aff^2 - (X_aff + X_aff)
(X_aff - o)^2 - o

X_aff^2 - (X_aff + X_aff) == (X_aff - o)^2 - o

range(X_aff^2 - (X_aff + X_aff))

X^2 - 2X
(X-1)^2 - 1   # true image, since only single use of interval

C = Affine(3.0, [4.0, 5.0])

range(C)

C = Affine(2., [1., -2, 3, 2])
D = Affine(1., [3., 0., -1.])

C + D

C - D

C = Affine(2., [2.])

C - C

D = Affine(1., [2.])

C = Affine(2., [1.])
C * C
range(C * C)
range(C)
range(C)^2

C * D
range(C+D)
range(C) + range(D)


text:
## using Float64
timer=TimerOutput("sparse pflap psnr reg"); # fix
method_kwargs =Dict(:timer => timer, :display => false, :max_repeats => 3, :match_source_histogram => false)
flow_est, source_reg, timer, results = test_registration_alg(sparse_pflap_psnr, img, imgw, flow, method_kwargs=method_kwargs, timer=timer);

## using Float16
timer=TimerOutput("sparse pflap psnr float 16 reg"); # fix
method_kwargs =Dict(:timer => timer, :display => false, :max_repeats => 3, :match_source_histogram => false)
flow_est, source_reg, timer, results = test_registration_alg(LAP_julia.sparse_pflap_psnr_float16, img, imgw, flow, method_kwargs=method_kwargs, timer=timer);

filter_half_size = 15
filter_size = filter_half_size*2 +1
img_size = (257,257)

sigma = Float16.((filter_half_size + 2) / 4)
centered_inds = centered(-filter_half_size:filter_half_size)
gaus = KernelFactors.gaussian(sigma, filter_size)
kern16 = ImageFiltering.kernelfactors((gaus,gaus))

img16 = rand(Float16, img_size)
holder16 = rand(Float16, img_size)

@code_warntype imfilter!(holder16, img16, kern16)

@benchmark imfilter!($holder16, $img16, $kern16)

# --------->> for some reason 16 is slow to filter.

##

sigma = Float64.((filter_half_size + 2) / 4)
centered_inds = centered(-filter_half_size:filter_half_size)
gaus = KernelFactors.gaussian(sigma, filter_size)
kern64 = ImageFiltering.kernelfactors((gaus,gaus))

img64 = rand(Float64, img_size)
holder64 = rand(Float64, img_size)

@code_warntype imfilter!(holder64, img64, kern64)

@benchmark imfilter!($holder64, $img64, $kern64)

using LinearAlgebra

A = [1 1; 2 2]

b = [1, 1]

A\b
qr(A) \ b
qr(A, Val(true)) \ b
lu(A) \ b
svd(A) \ b

A = Float16[0.0 0.0; 0.0 0.0]
b = Float16[-0.009766, 0.0]


text:
"""
    InfiniteOpt.add_measure_variable(model::JuMP.Model,
                                     var::InfiniteOpt.PointVariable,
                                     key::Val{:TransData}
                                     )::InfiniteOpt.GeneralVariableRef

Make a `PointVariableRef` and map it to the appropriate transcription variable
and return the `GeneralVariableRef`. This is an extension of
[`add_measure_variable`](@ref InfiniteOpt.add_measure_variable(::JuMP.Model,::Any,::Any))
for `TranscriptionOpt`.
"""
function InfiniteOpt.add_measure_variable(
    model::JuMP.Model,
    var::InfiniteOpt.PointVariable,
    key::Val{:TransData}
    )::InfiniteOpt.GeneralVariableRef
    # make negative index to not conflict with the InfiniteModel
    raw_index = transcription_data(model).last_point_index -= 1
    # make the reference and map it to a transcription variable
    ivref = var.infinite_variable_ref
    pvref = InfiniteOpt.GeneralVariableRef(JuMP.owner_model(ivref), raw_index,
                                           InfiniteOpt.PointVariableIndex)
    trans_var = lookup_by_support(model, ivref, var.parameter_values)
    transcription_data(model).finvar_mappings[pvref] = trans_var
    return pvref
end

"""
    InfiniteOpt.add_measure_variable(model::JuMP.Model,
                                     var::InfiniteOpt.SemiInfiniteVariable,
                                     key::Val{:TransData}
                                     )::InfiniteOpt.GeneralVariableRef

Make a `SemiInfiniteVariableRef` and add `var` to the transcription data
and return the `GeneralVariableRef`. This is an extension of
[`add_measure_variable`](@ref InfiniteOpt.add_measure_variable(::JuMP.Model,::Any,::Any))
for `TranscriptionOpt`. Note that `internal_semi_infinite_variable` is also extended
to be able to access the `var`.
"""
function InfiniteOpt.add_measure_variable(
    model::JuMP.Model,
    var::InfiniteOpt.SemiInfiniteVariable,
    key::Val{:TransData}
    )::InfiniteOpt.GeneralVariableRef
    # make negative index to not conflict with the InfiniteModel
    semi_infinite_vars = transcription_data(model).semi_infinite_vars
    raw_index = -1 * (length(semi_infinite_vars) + 1)
    # make the reference and map it to a transcription variable
    ivref = var.infinite_variable_ref
    rvref = InfiniteOpt.GeneralVariableRef(JuMP.owner_model(ivref), raw_index,
                                           InfiniteOpt.SemiInfiniteVariableIndex)
    push!(semi_infinite_vars, var)
    _set_semi_infinite_variable_mapping(model, var, rvref, InfiniteOpt._index_type(ivref))
    return rvref
end

"""
    InfiniteOpt.delete_semi_infinite_variable(model::JuMP.Model,
                                        vref::InfiniteOpt.SemiInfiniteVariableRef,
                                        key::Val{:TransData})::Nothing

This is an extension of
[`delete_semi_infinite_variable`](@ref InfiniteOpt.delete_semi_infinite_variable(::JuMP.Model, ::Any, ::Any))
for use in `TranscriptionOpt`. Here we do not delete semi-infinite variables once they
have been used since there is no performance gain for this paradigm and the
memory saving is small. Note this may change in the future.
"""
function InfiniteOpt.delete_semi_infinite_variable(
    model::JuMP.Model,
    rvref::InfiniteOpt.SemiInfiniteVariableRef,
    key::Val{:TransData}
    )::Nothing
    return
end


text:
const GLFW_LOCK_KEY_MODS = 0x00033004
const GLFW_MOD_CAPS_LOCK = 0x0010   # Caps Lock key is enabled
const GLFW_MOD_NUM_LOCK  = 0x0020   # Num Lock key is enabled

# keypad helper function
function keypadcmd(key::GLFW.Key, mods::Int32)::RobotCmd
    # Keypad keys have a different meaning if ctrl, alt, or both are pressed
    ctrl_or_alt = mods & (GLFW.MOD_CONTROL | GLFW.MOD_ALT) != 0
    ctrl_or_alt && return NO_COMMAND # ctrl, alt, or both are pressed

    # No other modifiers affect this set of keys
    key == GLFW.KEY_KP_SUBTRACT && return TOGGLE_TROT
    key == GLFW.KEY_KP_ADD      && return CYCLE_HOP

    # Only caps lock does not affect this set of keys
    unmodified = mods == 0 || mods == GLFW_MOD_CAPS_LOCK
    unmodified || return NO_COMMAND

    # The following comparisons assume that Num Lock is not toggled on
    key == GLFW.KEY_KP_7 && return INCREASE_HEIGHT    # Increase ride height
    key == GLFW.KEY_KP_8 && return PITCH_NOSE_UP      # Nose up, tail down
    key == GLFW.KEY_KP_9 && return INCREASE_VELOCITY  # Increase trot forward
    key == GLFW.KEY_KP_4 && return DECREASE_YAW       # Veer left
    key == GLFW.KEY_KP_6 && return INCREASE_YAW       # Veer right
    key == GLFW.KEY_KP_1 && return DECREASE_HEIGHT    # Decrease ride height
    key == GLFW.KEY_KP_2 && return PITCH_NOSE_DOWN    # Nose down, tail up
    key == GLFW.KEY_KP_3 && return DECREASE_VELOCITY  # Increase trot forward (reverse)
    key == GLFW.KEY_KP_0 && return ROLL_LEFT          # Roll body to the left
    key == GLFW.KEY_KP_DECIMAL && return ROLL_RIGHT   # Roll body to the right

    return NO_COMMAND
end

function keyboardcmd(key::GLFW.Key, mods::Int32)::RobotCmd
    # Ensure that key is in the shifted state (caps lock or shift key)
    keys_shifted = mods & GLFW.MOD_SHIFT > 0
    keys_shifted || return NO_COMMAND

    # Velocity PgUp / PgDn
    if     key == GLFW.KEY_PAGE_UP return INCREASE_VELOCITY
    elseif key == GLFW.KEY_PAGE_DOWN return DECREASE_VELOCITY

    # Height Home / End
    elseif key == GLFW.KEY_HOME return INCREASE_HEIGHT
    elseif key == GLFW.KEY_END return DECREASE_HEIGHT

    # Yaw left / right arrow
    elseif key == GLFW.KEY_LEFT return INCREASE_YAW
    elseif key == GLFW.KEY_RIGHT return DECREASE_YAW

    # Pitch up / down arrow
    elseif key == GLFW.KEY_UP return PITCH_NOSE_UP
    elseif key == GLFW.KEY_DOWN return PITCH_NOSE_DOWN

    # Roll left/right (Ins/Del)
    elseif key == GLFW.KEY_INSERT return ROLL_LEFT
    elseif key == GLFW.KEY_DELETE return ROLL_RIGHT

    # Toggle activate (`) / trot (-) / hop (+)
    elseif key == GLFW.KEY_1 return TOGGLE_ACTIVATION
    elseif key == GLFW.KEY_MINUS return TOGGLE_TROT
    elseif key == GLFW.KEY_EQUAL return CYCLE_HOP  # (+)

    # Turn left (<) / right (>)
    elseif key == GLFW.KEY_COMMA return TURN_LEFT   # Turn left while trotting in place
    elseif key == GLFW.KEY_PERIOD return TURN_RIGHT # Turn right while trotting in place

    end
end

in_range(key::Int, low::Int, high::Int) = key in low:high
is_keypad_robotcmd(key::GLFW.Key) = in_range(Int(key), Int(GLFW.KEY_KP_0), Int(GLFW.KEY_KP_ENTER))
is_keyboard_robotcmd(key::GLFW.Key) = in_range(Int(key), Int(GLFW.KEY_INSERT), Int(GLFW.KEY_END)) ||
            key == GLFW.KEY_1 || key == GLFW.KEY_MINUS || key == GLFW.KEY_EQUAL ||
            key == GLFW.KEY_COMMA || key == GLFW.KEY_PERIOD

is_robotcmd_turn(c::RobotCmd) = c == TURN_LEFT || c == TURN_RIGHT
# Repeat keys valid for velocity, height, roll, pitch, and yaw commands
is_robotcmd_real_valued(c::RobotCmd) = in_range(Int(c), Int(INCREASE_VELOCITY), Int(ROLL_RIGHT))

const keycmds = Dict{GLFW.Key, Function}(
    GLFW.KEY_F1=>(s)->begin  # help
        s.showhelp += 1
        if s.showhelp > 2 s.showhelp = 0 end
    end,
    GLFW.KEY_F2=>(s)->begin  # option
        s.showoption = !s.showoption
    end,
    GLFW.KEY_F3=>(s)->begin  # info
        s.showinfo = !s.showinfo
    end,
    GLFW.KEY_F4=>(s)->begin  # depth
        s.showdepth = !s.showdepth
    end,
    GLFW.KEY_F5=>(s)->begin  # toggle full screen
        s.showfullscreen = !s.showfullscreen
        s.showfullscreen ? GLFW.MaximizeWindow(s.window) : GLFW.RestoreWindow(s.window)
    end,
    #GLFW.KEY_F6=>(s)->begin  # stereo
    #   s.stereo = s.scn.stereo == MJCore.mjSTEREO_NONE ? mjSTEREO_QUADBUFFERED : MJCore.mjSTEREO_NONE
    #   s.scn[].stereo
    #end,
    GLFW.KEY_F7=>(s)->begin  # sensor figure
        s.showsensor = !s.showsensor
    end,
    GLFW.KEY_F8=>(s)->begin  # profiler
        s.showprofiler = !s.showprofiler
    end,
    GLFW.KEY_ENTER=>(s)->begin  # slow motion
        s.slowmotion = !s.slowmotion
        s.slowmotion ? println("Slow Motion Mode!") : println("Normal Speed Mode!")
    end,
    GLFW.KEY_SPACE=>(s)->begin  # pause
        s.paused = !s.paused
        s.paused ? println("Paused") : println("Running")
    end,
    GLFW.KEY_PAGE_UP=>(s)->begin    # previous keyreset
        s.keyreset = min(s.m.nkey - 1, s.keyreset + 1)
    end,
    GLFW.KEY_PAGE_DOWN=>(s)->begin  # next keyreset
        s.keyreset = max(-1, s.keyreset - 1)
    end,
    # continue with reset
    GLFW.KEY_BACKSPACE=>(s)->begin  # reset
        mj_resetData(s.m.m, s.d.d)
        if s.keyreset >= 0 && s.keyreset < s.m.nkey
            s.time = s.m.key_time[s.keyreset+1]
            s.d.qpos[:] = s.m.key_qpos[:,s.keyreset+1]
            s.d.qvel[:] = s.m.key_qvel[:,s.keyreset+1]
            s.d.act[:]  = s.m.key_act[:,s.keyreset+1]
        end
        mj_forward(s.m, s.d)
        #profilerupdate()
        sensorupdate(s)
    end,
    GLFW.KEY_RIGHT=>(s)->begin  # step forward
        if s.paused
            mj_step(s.m, s.d)
            #profilerupdate()
            sensorupdate(s)
        end
    end,
    GLFW.KEY_LEFT=>(s)->begin  # step back
    #    if s.paused
    #       dt = s.m.opt.timestep
    #       s.m.opt.timestep = -dt
    #       #cleartimers(s.d)
    #       mj_step(s.m, s.d)
    #       s.m.opt.timestep = dt
    #       #profilerupdate()
    #       sensorupdate(s)
    #    end
    end,
    GLFW.KEY_DOWN=>(s)->begin  # step forward 100
        if s.paused
            #cleartimers(d)
            for n=1:100 mj_step(s.m, s.d) end
            #profilerupdate()
            sensorupdate(s)
        end
    end,
    GLFW.KEY_UP=>(s)->begin  # step back 100
    #    if s.paused
    #       dt = s.m.opt.timestep
    #       s.m.opt.timestep = -dt
    #       #cleartimers(d)
    #       for n=1:100 mj_step(s.m, s.d) end
    #       s.m.opt.timestep = dt
    #       #profilerupdate()
    #       sensorupdate(s)
    #    end
    end,
    GLFW.KEY_ESCAPE=>(s)->begin  # free camera
        s.cam[].type = MJCore.CAMERA_FREE
    end,
    GLFW.KEY_EQUAL=>(s)->begin  # bigger font
        if fontscale < 200
            fontscale += 50
            mjr_makeContext(s.m.m, s.con, fontscale)
        end
    end,
    GLFW.KEY_MINUS=>(s)->begin  # smaller font
        if fontscale > 100
            fontscale -= 50
            mjr_makeContext(s.m.m, s.con, fontscale)
        end
    end,
    GLFW.KEY_LEFT_BRACKET=>(s)->begin  # '[' previous fixed camera or free
        fixedcamtype = s.cam[].type
        if s.m.ncam > 0 && fixedcamtype == MJCore.mjCAMERA_FIXED
            fixedcamid = s.cam[].fixedcamid
            if (fixedcamid  > 0)
                s.cam[].fixedcamid = fixedcamid-1
            elseif fixedcamid == 0
                s.cam[].type = MJCore.mjCAMERA_FREE
                s.cam[].fixedcamid = fixedcamid-1
            end
        end
    end,
    GLFW.KEY_RIGHT_BRACKET=>(s)->begin  # ']' next fixed camera
        if s.m.ncam > 0
            fixedcamtype = s.cam[].type
            fixedcamid = s.cam[].fixedcamid
            if fixedcamid < s.m.ncam - 1
                s.cam[].fixedcamid = fixedcamid+1
                s.cam[].type = MJCore.mjCAMERA_FIXED
            end
        end
    end,
    GLFW.KEY_SEMICOLON=>(s)->begin  # cycle over frame rendering modes
        frame = s.vopt.frame
        s.vopt[].frame = max(0, frame - 1)
    end,
    GLFW.KEY_APOSTROPHE=>(s)->begin  # cycle over frame rendering modes
        frame = s.vopt.frame
        s.vopt[].frame = min(MJCore.mjNFRAME-1, frame+1)
    end,
    GLFW.KEY_PERIOD=>(s)->begin  # cycle over label rendering modes
        label = s.vopt.label
        s.vopt[].label = max(0, label-1)
    end,
    GLFW.KEY_SLASH=>(s)->begin  # cycle over label rendering modes
        label = s.vopt.label
        s.vopt[].label = min(MJCore.mjNLABEL-1, label+1)
    end
)

function quit(s::mjSim)
    s.record !== nothing && finish_recording(s)
    GLFW.SetWindowShouldClose(window, true)
end

function keyboard_visualization_keys(s::mjSim, key::GLFW.Key, mods::Int32)
    #println("mjNVISFLAG: $(MJCore.mjNVISFLAG), mjVISSTRING: $(MJCore.mjVISSTRING)\nmjNRNDFLAG: $(MJCore.mjNRNDFLAG), mjRNDSTRING: $(MJCore.mjRNDSTRING), mjNGROUP: $(MJCore.mjNGROUP)")

    # toggle visualization flag
    # NVISFLAG: 22, VISSTRING: ["Convex Hull" "0" "H"; "Texture" "1" "X"; "Joint" "0" "J"; "Actuator" "0" "U"; "Camera" "0" "Q"; "Light" "0" "Z"; "Tendon" "0" "V"; "Range Finder" "0" "Y"; "Constraint" "0" "N"; "Inertia" "0" "I"; "SCL Inertia" "0" "S"; "Perturb Force" "0" "B"; "Perturb Object" "1" "O"; "Contact Point" "0" "C"; "Contact Force" "0" "F"; "Contact Split" "0" "P"; "Transparent" "0" "T"; "Auto Connect" "0" "A"; "Center of Mass" "0" "M"; "Select Point" "0" "E"; "Static Body" "0" "D"; "Skin" "0" ";"]
    if key != GLFW.KEY_S
        for i=1:MJCore.mjNVISFLAG
            #name = MJCore.mjVISSTRING[1, i]
            #println("Comparing $(Int(key)) to $(MJCore.mjVISSTRING[3, i][1])")
            if Int(key) == Int(MJCore.mjVISSTRING[3, i][1])
                flags = MVector(s.vopt[].flags)
                flags[i] = flags[i] == 0 ? 1 : 0
                s.vopt[].flags = flags
                return
            end
        end
    end

    # toggle rendering flag
    # NRNDFLAG: 9,  RNDSTRING: ["Shadow" "1" "S"; "Wireframe" "0" "W"; "Reflection" "1" "R"; "Additive" "0" "L"; "Skybox" "1" "K"; "Fog" "0" "G"; "Haze" "1" "/"; "Segment" "0" ","; "Id Color" "0" "."], NGROUP: 6
    for i=1:MJCore.mjNRNDFLAG
        if Int(key) == Int(MJCore.mjRNDSTRING[3, i][1])
            flags = MVector(s.scn[].flags)
            flags[i] = flags[i] == 0 ? 1 : 0
            s.scn[].flags = flags
            return
        end
    end

    # toggle geom/site group
    for i=1:MJCore.mjNGROUP
        if Int(key) == i + Int('0')
            if mods & GLFW.MOD_SHIFT == true
                sitegroup = MVector(s.vopt[].sitegroup)
                sitegroup[i] = sitegroup[i] > 0 ? 0 : 1
                s.vopt[].sitegroup[i] = sitegroup
                return
            else
                geomgroup = MVector(s.vopt[].geomgroup)
                geomgroup[i] = geomgroup[i] > 0 ? 0 : 1
                s.vopt[].geomgroup = geomgroup
                return
            end
        end
    end
end

function keyboard(s::mjSim, window::GLFW.Window,
                key::GLFW.Key, scancode::Int32, act::GLFW.Action, mods::Int32)

    # Check for robot command key
    if is_keypad_robotcmd(key)                             # KeyPad keys
        robotcmd = keypadcmd(key, mods)
    elseif is_keyboard_robotcmd(key)
        robotcmd = keyboardcmd(key, mods)
    else
        robotcmd = NO_COMMAND
    end

    # Handle key
    if act == GLFW.RELEASE                              # Key release
        # Stop turning on key release
        is_robotcmd_turn(robotcmd) && end_turn(s.robot)
    elseif act == GLFW.REPEAT                           # Key is held down
        # Repeat keys valid for velocity, height, roll, pitch, and yaw commands
        is_robotcmd_real_valued(robotcmd) && execute_robotcmd(s, robotcmd)
    elseif robotcmd != NO_COMMAND                       # Robot command
        execute_robotcmd(s, robotcmd)
    elseif mods & GLFW.MOD_CONTROL == GLFW.MOD_CONTROL  # Ctrl key commands
        if     key == GLFW.KEY_A alignscale(s)
        #elseif key == GLFW.KEY_L reloadmodel(s)
        elseif key == GLFW.KEY_P print_state(s)
        elseif key == GLFW.KEY_Q quit(s)
        elseif key == GLFW.KEY_V record_video(s)
        elseif key == GLFW.KEY_RIGHT_BRACKET s.showcam = !s.showcam
        end
    elseif mods & GLFW.MOD_ALT == GLFW.MOD_ALT          # Alt key commands
        # Placeholder
    elseif mods & GLFW.MOD_SUPER == GLFW.MOD_SUPER      # Super key commands
        # Placeholder
    else                                                # Other unshifed or shifed key commands
        try
            # Attempt lookup in keycmds dict
            keycmds[key](s) # call anonymous function in keycmds Dict
        catch
            # Key was not in the dictionary or keycmds otherwise threw an exception
            keyboard_visualization_keys(s, key, mods)
        end
    end
end


text:
extract_name(line) = join(split(line)[2:end-1], " ")

start_testset!(dest, line) =
    println(dest, "@testset \"", extract_name(line), "\" begin")
finish_testset!(dest) = println(dest, "end\n")

begins_testset(line) = startswith(line, "testcase")
finishes_testset(line) = line == "}"
begins_decoration_testset(line) = occursin("dec", extract_name(line))
is_comment(line) = startswith(line, "//")
begins_multiline_comment(line) = startswith(line, "/*")

skip_block(src, stopline) = readuntil(src, "\n" * stopline * "\n")
skip_multiline_comment(src) = skip_block(src, "*/")
skip_testset(src) = skip_block(src, "}")

function translate!(dest::IO, src::IO)
    while !eof(src)
        line = rstrip(strip(readline(src)), ';')
        if isempty(line) || is_comment(line)
            nothing
        elseif begins_multiline_comment(line)
            skip_multiline_comment(src)
        elseif begins_testset(line)
            if begins_decoration_testset(line)
                skip_testset(src)
            else
                start_testset!(dest, line)
            end
        elseif finishes_testset(line)
            finish_testset!(dest)
        else
            itl_test = line
            translate!(dest, itl_test)
        end
    end
end

tests_decorated_intervals(case) =
    occursin("]_", case) ||
    occursin("d-numsToInterval", case) ||
    occursin("setDec", case) ||
    occursin("intervalPart", case)
tests_string_input(case) = occursin("textToInterval", case)
tests_logging(case) = occursin("signal", case)

skip_testcase(case) =
    tests_decorated_intervals(case) ||
    tests_string_input(case)

indent(io) = print(io, repeat(' ', 4))

function translate!(dest::IO, itl_test::AbstractString)
    skip_testcase(itl_test) && return
    jl_test = tests_logging(itl_test) ?
        translate_logtest(itl_test) : translate(itl_test)
    indent(dest)
    println(dest, jl_test)
end

function translate_logtest(itl_test)
    itl_test, exc = split(itl_test, " signal ")
    exc != "UndefinedOperation" && return "# " * itl_test
    lhs, rhs = split(itl_test, "=")
    jl_test = rebuild(rebuild_lhs(lhs), rebuild_rhs(rhs))
    return "@test_throws AssertionError " * jl_test
end

isbroken(expr) = !eval(Meta.parse(expr))

"""

This function parses a line into julia code, e.g.

```
add [1, 2] [1, 2] = [2, 4]
```

is parsed into
```
@test +(Interval(1, 2), Interval(1, 2)) === Interval(2, 4)
```
"""
function translate(itl_test)
    lhs, rhs = split(itl_test, "=")
    jl_test = rebuild(rebuild_lhs(lhs), rebuild_rhs(rhs))
    try
        return ifelse(isbroken(jl_test), "@test_broken ", "@test ") * jl_test
    catch
        @warn "caused exception: " * jl_test
        return "# " * itl_test
    end
end

const interval_pattern = r"\[([^\]]+)\]"

function rebuild_lhs(lhs)
    lhs = strip(lhs)
    fname, args = split(lhs, limit = 2)

    # input numbers
    args = replace(args, "infinity" => "Inf")
    args = replace(args, "X" => "x")
    if fname == "b-numsToInterval"
        args = replace(args, ' ' => ',')
        return "Interval($args)"
    end

    # input intervals
    for m in eachmatch(interval_pattern, args)
        args = replace(args, m.match => translate_interval(m[1]))
    end
    args = replace(args, " " => ", ")
    args = replace(args, ",," => ",")
    args = replace(args, "{" => "[")
    args = replace(args, "}" => "]")
    return functions[fname](args)
end

isintstring(x) = !isnothing(tryparse(Int, x))
floatstring(x) = x * ifelse(isintstring(x), ".0", "")

function rebuild_rhs(rhs)
    rhs = strip(rhs)
    rhs = replace(rhs, "infinity" => "Inf")
    rhs = replace(rhs, "X" => "x")
    if '[' ∉ rhs # one or more scalar/boolean values separated by space
        return map(floatstring, split(rhs))
    else # one or more intervals
        intervals = map(
            m -> translate_interval(m[1]),
            eachmatch(interval_pattern, rhs)
        )
        return intervals
    end
end

const special_intervals = Dict(
    "nai" => "nai()",
    "entire" => "entireinterval()",
    "empty" => "emptyinterval()"
)

translate_interval(ival) =
    haskey(special_intervals, ival) ?
    special_intervals[ival] :
    "interval($ival)"

const relaxed_functions = [
        "sinh",
        "cosh",
        "tanh",
        "coth",
        "sech",
        "csch",
        "asinh",
        "acosh",
        "atanh",
        "acoth",
        "atan",
        "cot",
        "acot",
        "sec",
        "csc",
        "exp2",
        "exp10",
        "^",
        "cbrt"
    ]

only_requires_validity(lhs) =
    any(startswith(lhs, f * "(") for f in relaxed_functions)

function rebuild(lhs, rhs::AbstractString)
    rhs == "nai()" && return "isnai($lhs)"
    rhs == "NaN" && return "isnan($lhs)"
    rhs == "true" && return lhs
    rhs == "false" && return lhs[end] == ')' ? "!" * lhs : "!($lhs)"
    only_requires_validity(lhs) && return "$lhs ⊇ $rhs"
    return "$lhs == $rhs"
end

function rebuild(lhs, rhs::Vector)
    length(rhs) == 1 && return rebuild(lhs, rhs[1])
    expr = [rebuild(lhs*"[$i]", r) for (i, r) in enumerate(rhs)]
    return join(expr, " && ")
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

@nospecialize

"""
    call::CallMeta

A simple struct that captures both the return type (`call.rt`)
and any additional information (`call.info`) for a given generic call.
"""
struct CallMeta
    rt::Any
    effects::Effects
    info::Any
end

"""
    info::MethodMatchInfo

Captures the result of a `:jl_matching_methods` lookup for the given call (`info.results`).
This info may then be used by the optimizer to inline the matches, without having
to re-consult the method table. This info is illegal on any statement that is
not a call to a generic function.
"""
struct MethodMatchInfo
    results::MethodLookupResult
end

"""
    info::UnionSplitInfo

If inference decides to partition the method search space by splitting unions,
it will issue a method lookup query for each such partition. This info indicates
that such partitioning happened and wraps the corresponding `MethodMatchInfo` for
each partition (`info.matches::Vector{MethodMatchInfo}`).
This info is illegal on any statement that is not a call to a generic function.
"""
struct UnionSplitInfo
    matches::Vector{MethodMatchInfo}
end

nmatches(info::MethodMatchInfo) = length(info.results)
function nmatches(info::UnionSplitInfo)
    n = 0
    for mminfo in info.matches
        n += nmatches(mminfo)
    end
    return n
end

struct ConstPropResult
    result::InferenceResult
end

struct ConcreteResult
    mi::MethodInstance
    effects::Effects
    result
    ConcreteResult(mi::MethodInstance, effects::Effects) = new(mi, effects)
    ConcreteResult(mi::MethodInstance, effects::Effects, @nospecialize val) = new(mi, effects, val)
end

const ConstResult = Union{ConstPropResult,ConcreteResult}

"""
    info::ConstCallInfo

The precision of this call was improved using constant information.
In addition to the original call information `info.call`, this info also keeps the results
of constant inference `info.results::Vector{Union{Nothing,ConstResult}}`.
"""
struct ConstCallInfo
    call::Union{MethodMatchInfo,UnionSplitInfo}
    results::Vector{Union{Nothing,ConstResult}}
end

"""
    info::MethodResultPure

This struct represents a method result constant was proven to be
effect-free, including being no-throw (typically because the value was computed
by calling an `@pure` function).
"""
struct MethodResultPure
    info::Any
end
let instance = MethodResultPure(false)
    global MethodResultPure
    MethodResultPure() = instance
end

"""
    info::AbstractIterationInfo

Captures all the information for abstract iteration analysis of a single value.
Each (abstract) call to `iterate`, corresponds to one entry in `info.each::Vector{CallMeta}`.
"""
struct AbstractIterationInfo
    each::Vector{CallMeta}
end

const MaybeAbstractIterationInfo = Union{Nothing, AbstractIterationInfo}

"""
    info::ApplyCallInfo

This info applies to any call of `_apply_iterate(...)` and captures both the
info of the actual call being applied and the info for any implicit call
to the `iterate` function. Note that it is possible for the call itself
to be yet another `_apply_iterate`, in which case the `info.call` field will
be another `ApplyCallInfo`. This info is illegal on any statement that is
not an `_apply_iterate` call.
"""
struct ApplyCallInfo
    # The info for the call itself
    call::Any
    # AbstractIterationInfo for each argument, if applicable
    arginfo::Vector{MaybeAbstractIterationInfo}
end

"""
    info::UnionSplitApplyCallInfo

Like `UnionSplitInfo`, but for `ApplyCallInfo` rather than `MethodMatchInfo`.
This info is illegal on any statement that is not an `_apply_iterate` call.
"""
struct UnionSplitApplyCallInfo
    infos::Vector{ApplyCallInfo}
end

"""
    info::InvokeCallInfo

Represents a resolved call to `Core.invoke`, carrying the `info.match::MethodMatch` of
the method that has been processed.
Optionally keeps `info.result::InferenceResult` that keeps constant information.
"""
struct InvokeCallInfo
    match::MethodMatch
    result::Union{Nothing,ConstResult}
end

"""
    info::OpaqueClosureCallInfo

Represents a resolved call of opaque closure, carrying the `info.match::MethodMatch` of
the method that has been processed.
Optionally keeps `info.result::InferenceResult` that keeps constant information.
"""
struct OpaqueClosureCallInfo
    match::MethodMatch
    result::Union{Nothing,ConstResult}
end

"""
    info::OpaqueClosureCreateInfo

This info may be constructed upon opaque closure construction, with `info.unspec::CallMeta`
carrying out inference result of an unreal, partially specialized call (i.e. specialized on
the closure environment, but not on the argument types of the opaque closure) in order to
allow the optimizer to rewrite the return type parameter of the `OpaqueClosure` based on it.
"""
struct OpaqueClosureCreateInfo
    unspec::CallMeta
    function OpaqueClosureCreateInfo(unspec::CallMeta)
        @assert isa(unspec.info, OpaqueClosureCallInfo)
        return new(unspec)
    end
end

# Stmt infos that are used by external consumers, but not by optimization.
# These are not produced by default and must be explicitly opted into by
# the AbstractInterpreter.

"""
    info::ReturnTypeCallInfo

Represents a resolved call of `Core.Compiler.return_type`.
`info.call` wraps the info corresponding to the call that `Core.Compiler.return_type` call
was supposed to analyze.
"""
struct ReturnTypeCallInfo
    info::Any
end

"""
    info::FinalizerInfo

Represents the information of a potential (later) call to the finalizer on the given
object type.
"""
struct FinalizerInfo
    info::Any
    effects::Effects
end

@specialize

max_stars_repo_name: Physolia/julia

max_issues_repo_name: Eurus-Holmes/julia

max_forks_repo_name: Nexuscompute/julia


text:
# # Custom potential

# We solve the 1D Gross-Pitaevskii equation with a custom potential.
# This is similar to [Gross-Pitaevskii equation in one dimension](@ref) and we
# show how to define local potentials attached to atoms, which allows for
# instance to compute forces.
using DFTK
using LinearAlgebra

# First, we define a new element which represents a nucleus generating
# a Gaussian potential.
struct ElementGaussian <: DFTK.Element
    α  # Prefactor
    L  # Width of the Gaussian nucleus
end

# We extend the two methods providing access to the real and Fourier
# representation of the potential to DFTK.
function DFTK.local_potential_real(el::ElementGaussian, r::Real)
    -el.α / (√(2π) * el.L) * exp(- (r / el.L)^2 / 2)
end
function DFTK.local_potential_fourier(el::ElementGaussian, q::Real)
    ## = ∫ V(r) exp(-ix⋅q) dx
    -el.α * exp(- (q * el.L)^2 / 2)
end

# We set up the lattice. For a 1D case we supply two zero lattice vectors
a = 10
lattice = a .* [[1 0 0.]; [0 0 0]; [0 0 0]];

# In this example, we want to generate two Gaussian potentials generated by
# two "nuclei" localized at positions ``x_1`` and ``x_2``, that are expressed in
# ``[0,1)`` in fractional coordinates. ``|x_1 - x_2|`` should be different from
# ``0.5`` to break symmetry and get nonzero forces.
x1 = 0.2
x2 = 0.8
nucleus = ElementGaussian(1.0, 0.5)
atoms = [nucleus => [[x1, 0, 0], [x2, 0, 0]]];

# We setup a Gross-Pitaevskii model
C = 1.0
α = 2;
n_electrons = 1  # Increase this for fun
terms = [Kinetic(),
         AtomicLocal(),
         PowerNonlinearity(C, α),
]
model = Model(lattice; atoms=atoms, n_electrons=n_electrons, terms=terms,
              spin_polarization=:spinless);  # use "spinless electrons"

# We discretize using a moderate Ecut and run a SCF algorithm to compute forces
# afterwards. As there is no ionic charge associated to `nucleus` we have to specify
# a starting density and we choose to start from a zero density.
basis = PlaneWaveBasis(model; Ecut=500, kgrid=(1, 1, 1))
ρ = zeros(eltype(basis), basis.fft_size..., 1)
scfres = self_consistent_field(basis, tol=1e-8, ρ=ρ)
scfres.energies
# Computing the forces can then be done as usual:
hcat(compute_forces(scfres)...)

# Extract the converged total local potential
tot_local_pot = DFTK.total_local_potential(scfres.ham)[:, 1, 1]; # use only dimension 1

# Extract other quantities before plotting them
ρ = scfres.ρ[:, 1, 1, 1]  # converged density, first spin component
ψ_fourier = scfres.ψ[1][:, 1];    # first k-point, all G components, first eigenvector
ψ = G_to_r(basis, basis.kpoints[1], ψ_fourier)[:, 1, 1]
ψ /= (ψ[div(end, 2)] / abs(ψ[div(end, 2)]));

using Plots
x = a * vec(first.(DFTK.r_vectors(basis)))
p = plot(x, real.(ψ), label="real(ψ)")
plot!(p, x, imag.(ψ), label="imag(ψ)")
plot!(p, x, ρ, label="ρ")
plot!(p, x, tot_local_pot, label="tot local pot")


text:
function BatchLSEstimator(arg0::LeastSquaresOptimizer, arg1::Vector{IntegratedPropagatorBuilder})
    return BatchLSEstimator((LeastSquaresOptimizer, Vector{IntegratedPropagatorBuilder}), arg0, arg1)
end

function add_measurement(obj::BatchLSEstimator, arg0::ObservedMeasurement)
    return jcall(obj, "addMeasurement", void, (ObservedMeasurement,), arg0)
end

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

function estimate(obj::BatchLSEstimator)
    return jcall(obj, "estimate", Vector{AbstractIntegratedPropagator}, ())
end

function get_class(obj::Object)
    return jcall(obj, "getClass", Class, ())
end

function get_evaluations_count(obj::BatchLSEstimator)
    return jcall(obj, "getEvaluationsCount", jint, ())
end

function get_iterations_count(obj::BatchLSEstimator)
    return jcall(obj, "getIterationsCount", jint, ())
end

function get_last_estimations(obj::BatchLSEstimator)
    return jcall(obj, "getLastEstimations", Map, ())
end

function get_measurements_parameters_drivers(obj::BatchLSEstimator, arg0::jboolean)
    return jcall(obj, "getMeasurementsParametersDrivers", ParameterDriversList, (jboolean,), arg0)
end

function get_optimum(obj::BatchLSEstimator)
    return jcall(obj, "getOptimum", LeastSquaresOptimizer_Optimum, ())
end

function get_orbital_parameters_drivers(obj::BatchLSEstimator, arg0::jboolean)
    return jcall(obj, "getOrbitalParametersDrivers", ParameterDriversList, (jboolean,), arg0)
end

function get_physical_covariances(obj::BatchLSEstimator, arg0::jdouble)
    return jcall(obj, "getPhysicalCovariances", RealMatrix, (jdouble,), arg0)
end

function get_propagator_parameters_drivers(obj::BatchLSEstimator, arg0::jboolean)
    return jcall(obj, "getPropagatorParametersDrivers", ParameterDriversList, (jboolean,), arg0)
end

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

function notify(obj::Object)
    return jcall(obj, "notify", void, ())
end

function notify_all(obj::Object)
    return jcall(obj, "notifyAll", void, ())
end

function set_convergence_checker(obj::BatchLSEstimator, arg0::ConvergenceChecker)
    return jcall(obj, "setConvergenceChecker", void, (ConvergenceChecker,), arg0)
end

function set_max_evaluations(obj::BatchLSEstimator, arg0::jint)
    return jcall(obj, "setMaxEvaluations", void, (jint,), arg0)
end

function set_max_iterations(obj::BatchLSEstimator, arg0::jint)
    return jcall(obj, "setMaxIterations", void, (jint,), arg0)
end

function set_observer(obj::BatchLSEstimator, arg0::BatchLSObserver)
    return jcall(obj, "setObserver", void, (BatchLSObserver,), arg0)
end

function set_parameters_convergence_threshold(obj::BatchLSEstimator, arg0::jdouble)
    return jcall(obj, "setParametersConvergenceThreshold", void, (jdouble,), arg0)
end

function to_string(obj::Object)
    return jcall(obj, "toString", JString, ())
end

function wait(obj::Object)
    return jcall(obj, "wait", void, ())
end

function wait(obj::Object, arg0::jlong)
    return jcall(obj, "wait", void, (jlong,), arg0)
end

function wait(obj::Object, arg0::jlong, arg1::jint)
    return jcall(obj, "wait", void, (jlong, jint), arg0, arg1)
end