text:
function read_data(instance::AbstractString)
  data = Data()
  data.constraint_uncertainty = false
  file = readdlm(string(dirname(@__FILE__),"/../data/SP/",instance))
  println("********** Read instance $(file[2,1]) ****************")
  data.nV = file[4,1]
  data.n = file[6,1]
  coord = file[8:(7+data.nV),1:2]
  data.A = file[(9+data.nV):(8+data.nV+data.n),1:data.nV]
  data.s = file[10+data.nV+data.n,1]
  data.t = file[12+data.nV+data.n,1]
  data.b=zeros(Int64,data.nV)
  for i in 1:data.nV
    data.b[i]=0
    if i==data.s
      data.b[i]=1
    elseif i==data.t
      data.b[i]=-1
    end
  end
  data.from = Vector{Int}(undef,data.n)
  data.to = Vector{Int}(undef,data.n)
  for i in 1:data.n
    data.from[i] = findall(data.A[i,1:data.nV] .== 1)[1]
    data.to[i] = findall(data.A[i,1:data.nV] .== -1)[1]
  end
  data.δ⁺ = Array{Set,1}(undef,data.nV)
  data.δ⁻ = Array{Set,1}(undef,data.nV)
  for i in 1:data.nV
    data.δ⁺[i] = Set()
    data.δ⁻[i] = Set()
  end
  for i in 1:data.n
    push!(data.δ⁺[data.from[i]],i)
    push!(data.δ⁻[data.to[i]],i)
  end
  data.d = ones(data.n)
  data.hc = ones(data.n)
  for i in 1:data.n
    data.hc[i] = norm(coord[data.from[i],1:2]-coord[data.to[i],1:2])
    data.d[i] = data.hc[i]/2
  end
  data.ncons = 0
  return data
end

#-----------------------------------------------------------------------------------

mutable struct Data
  n::Int
  hc::Vector{Float64} # nominal cost
  d::Vector{Float64} # deviation
  from::Vector{Int}
  to::Vector{Int}
  δ⁺::Array{Set,1} #outgoing arcs
  δ⁻::Array{Set,1} #incoming arcs
  nV::Int
  A::Matrix{Int}
  s::Int
  t::Int
  b::Vector{Int} #RHS
  constraint_uncertainty::Bool
  ncons::Int
  heuristic_ub::Float64
  Data() = new()
end

#-----------------------------------------------------------------------------------

function h(x)
  return 1
end

#-----------------------------------------------------------------------------------

function get_all_SP(sp,ub)
  function getAllPathsUtil(arc::Int, u::Int, t::Int, visited::BitArray{1}, path::BitArray{1}, cost::Float64)
    # Mark the current node as visited and store in path
    if cost + sp[u] < ub
      visited[u] = true
      if arc > 0
        path[arc] = true
      end
      # If current vertex is same as destination, then store current path
      if u == t
        push!(X,copy(path))
      else
        for new_arc in data.δ⁺[u]
          if !visited[data.to[new_arc]]
            cost = cost + data.hc[new_arc]
            getAllPathsUtil(new_arc, data.to[new_arc], t, visited, path, cost)
            cost = cost - data.hc[new_arc]
          end
        end
      end
      # Remove current vertex from path[] and mark it as unvisited
      if arc > 0
        path[arc] = false
      end
      visited[u] = false
      return X
    end
  end
  X = Vector{BitArray{1}}()
  path = falses(data.n)
  cost = 0.0
  # Mark all the vertices as not visited
  visited = falses(data.nV)
  # Call the recursive helper function to print all paths
  getAllPathsUtil(0, data.s, data.t, visited, path, cost)
  return X
end

#-----------------------------------------------------------------------------------

function generate_solutions(ub)
  g = DiGraph(data.nV)
  for i in 1:data.n
    add_edge!(g, data.from[i], data.to[i])
  end
  dist = zeros(data.nV,data.nV)
  for i in 1:data.n
    dist[data.from[i],data.to[i]] = data.hc[i]
  end
  sp = Vector(undef,data.nV)
  for v in 1:data.nV
    all_dist = dijkstra_shortest_paths(g, v, dist).dists
    sp[v] = all_dist[data.t]
  end
  print("Computing X ")
  @time X_bool = get_all_SP(sp,ub) # set of all feasible vectors for theSP
  println("|X|=$(length(X_bool))")
  non_dom = trues(length(X_bool))
  @timeit to "dominance" begin
    order = sortperm(data.d, rev=true)
    d = data.d[order]
    arcs_used = Vector{Vector{Int64}}()
    nominal_costs = Vector{Float64}()
    for x in X_bool
      arcs = findall(x)
      push!(nominal_costs,sum(data.hc[arcs]))
      push!(arcs_used,findall(x[order]))
    end
    for i in 1:length(X_bool)
      ci = nominal_costs[i]
      ai = arcs_used[i]
      for j in 1:length(X_bool)
        cj = nominal_costs[j]
        aj = arcs_used[j]
        arcs = setdiff(aj,ai)
        if length(arcs) > 0
          arcs = arcs[1:min(Γ,length(arcs))]
        end
        # We remove solution i if there exists a solution j that is always chaper than solution i. This is exatly what strong dominance means.
        if cj - ci + sum(d[arcs]) < 0
          non_dom[i] = false
          break
        end
      end
    end
  end
  X_bool = X_bool[findall(non_dom)]
  println("After dominance check: |X|=$(length(X_bool))")
  return X_bool
end

#-----------------------------------------------------------------------------------

function add_constraints_X(model::Model, x::Array{VariableRef,1})
  @constraint(model, [n in 1:data.nV], sum(x[i] for i in data.δ⁺[n])-sum(x[i] for i in data.δ⁻[n]) ≥ data.b[n])
end

#-----------------------------------------------------------------------------------

"Create cc with k rows"
function populate_scenarios()
  cc = Array{Float64,2}(undef,1,data.n)
  for l in N cc[1,l] = data.hc[l] end
  for i = 1:k-1
    delta = zeros(data.n)
    for j=1:Γ
      delta[i+j]=1
    end
    new_cost = [data.hc[l] + data.d[l]*delta[l] for l in N]
    cc = [cc; new_cost']
  end
  # Notice the first argument returned is aa, which is empty for this application
  return [], cc
end

#-----------------------------------------------------------------------------------

function add_scenario(δ)
  return [data.hc[l] + data.d[l]*δ[l] for l in N], []
end

#-----------------------------------------------------------------------------------

function build_and_solve_separation(x,v)
  time_remaining = TIME_LIM - (time()-STARTTIME)
  # Compute the worst-case cost:
  # max{η : ξ \in Ξ, η ≤ f(ξ,x^k), ∀ k in [k]}
  separate = create_model(max(0,time_remaining))
  @variable(separate, η)
  @variable(separate, 0 ≤ ξ[N] ≤ 1,container=Array)
  @constraint(separate, sep_cost[k in K], η <= sum((data.hc[l] + data.d[l]*ξ[l])*x[k,l] for l in N))
  @constraint(separate, sum(ξ[l] for l in N) <= Γ)
  @objective(separate, Max, η)
  optimize!(separate)
  if termination_status(separate)==MOI.OPTIMAL
    return value.(ξ),objective_value(separate)
  else
    return [],-Inf
  end
end

#-----------------------------------------------------------------------------------

function heuristic_dualization()
  N = 1:data.n
  heursol = +Inf
  alphas = zeros(k)
  for j in K
    alphas[j] = (2*j)/((k)*(k+1))
  end
  timexstep = 0
  timealphastep = 0
  while true
    xstep = create_model(max(0,TIME_LIM-timexstep-timealphastep))
    @variable(xstep, 0 ≤ w)
    @variable(xstep, 0 ≤ v[i in N])
    @variable(xstep, x[k in K, i in N], Bin, container=Array)
    @objective(xstep,  Min, sum(sum(data.hc[i]*x[k,i]*alphas[k] for i in N) for k in K) + Γ*w + sum(v[i] for i in N))
    @constraint(xstep, dualconst[i in N], w + v[i] ≥ sum(data.d[i]*x[k,i]*alphas[k] for k in K))
    for k in K add_constraints_X(xstep,x[k,:]) end
    timexstep += @elapsed optimize!(xstep)
    if primal_status(xstep) == MOI.FEASIBLE_POINT
      global xsol = JuMP.value.(x)
      if heursol <= JuMP.objective_value(xstep) break end
      heursol = JuMP.objective_value(xstep)
    end
    αstep = create_model(max(0,TIME_LIM-timexstep-timealphastep))
    @variable(αstep, 0 ≤ w)
    @variable(αstep, 0 ≤ α[k in K])
    @variable(αstep, 0 ≤ v[i in N])
    @objective(αstep, Min, sum(sum(data.hc[i]*xsol[k,i]*α[k] for i in N) for k in K) + Γ*w + sum(v[i] for i in N))
    @constraint(αstep, dualconst[i in N], w + v[i] ≥ sum(data.d[i]*xsol[k,i]*α[k] for k in K))
    @constraint(αstep, sum(α[j] for j in K) == 1)
    optimize!(αstep)
    timealphastep += @elapsed optimize!(αstep)
    if primal_status(αstep) == MOI.FEASIBLE_POINT
      alphas=JuMP.value.(α)
    end
    if timexstep+timealphastep > TIME_LIM break end
  end
  return heursol,xsol
end

#-----------------------------------------------------------------------------------

function exact_dualization()
  dualized = create_model(TIME_LIM)
  @variable(dualized, 0 ≤ w)
  @variable(dualized, 0 ≤ v[i in N])
  @variable(dualized, 0 ≤ α[k in K] ≤ 1)
  @variable(dualized, 0 ≤ ρ[k in K,i in N] ≤ 1)
  @variable(dualized, x[k in K, i in N], Bin,container=Array)
  for k in K  add_constraints_X(dualized,x[k,:]) end
  @constraint(dualized, sum(α[k] for k in K) == 1)
  @constraint(dualized,[i in N], w + v[i] ≥ sum(data.d[i]ρ[k,i] for k in K))
  @constraint(dualized,[i in N,k in K], ρ[k,i] ≤ x[k,i])
  @constraint(dualized,[i in N,k in K], ρ[k,i] ≤ α[k])
  @constraint(dualized,[i in N,k in K], ρ[k,i] ≥ x[k,i]+α[k]-1)
  @objective(dualized,Min, sum(sum(data.hc[i]*ρ[k,i] for i in N) for k in K) + Γ*w + sum(v[i] for i in N))
  optimize!(dualized)
  return JuMP.objective_value(dualized), MOI.get(dualized, MOI.RelativeGap()), value.(x)
end

max_stars_repo_name: mjposs/min-max-min

max_issues_repo_name: mjposs/min-max-min

max_forks_repo_name: mjposs/min-max-min


text:
using Test

for path in readdir(isinteractive() ? pwd() : @__FILE__)
    if isdir(path)
        include(joinpath(path, "test", "runtests.jl"))
    end
end

ext: jl

max_stars_repo_name: stellartux/Esolangs

max_issues_repo_name: stellartux/Esolangs

max_forks_repo_name: stellartux/Esolangs


text:
@inline function predictive_stepsize_controller!(integrator, alg)
  # Gustafsson predictive stepsize controller
  if iszero(integrator.EEst)
    q = inv(integrator.opts.qmax)
  else
    gamma = integrator.opts.gamma
    if typeof(alg) <: Union{RKC,IRKC,SERK2}
      fac = gamma
    else
      if alg isa RadauIIA5
        @unpack nl_iters = integrator.cache
        @unpack max_iter = alg
      else
        @unpack nl_iters, max_iter = integrator.cache.nlsolver
      end
      fac = min(gamma,(1+2*max_iter)*gamma/(nl_iters+2*max_iter))
    end
    expo = 1/(get_current_adaptive_order(integrator.alg,integrator.cache)+1)
    qtmp = (integrator.EEst^expo)/fac
    @fastmath q = max(inv(integrator.opts.qmax),min(inv(integrator.opts.qmin),qtmp))
    integrator.qold = q
  end
  q
end

@inline function standard_stepsize_controller!(integrator, alg)
  # Standard stepsize controller
  if iszero(integrator.EEst)
    q = inv(integrator.opts.qmax)
  else
    qtmp = integrator.EEst^(1/(get_current_adaptive_order(integrator.alg,integrator.cache)+1))/integrator.opts.gamma
    @fastmath q = max(inv(integrator.opts.qmax),min(inv(integrator.opts.qmin),qtmp))
    integrator.qold = integrator.dt/q
  end
  q
end

@inline function PI_stepsize_controller!(integrator, alg)
  # PI-controller
  EEst,beta1,q11,qold,beta2 = integrator.EEst, integrator.opts.beta1, integrator.q11,integrator.qold,integrator.opts.beta2
  if iszero(EEst)
    q = inv(integrator.opts.qmax)
  else
    @fastmath q11 = EEst^beta1
    @fastmath q = q11/(qold^beta2)
    integrator.q11 = q11
    @fastmath q = max(inv(integrator.opts.qmax),min(inv(integrator.opts.qmin),q/integrator.opts.gamma))
  end
  q
end

function stepsize_controller!(integrator,alg)
  if ispredictive(alg)
    predictive_stepsize_controller!(integrator, alg)
  elseif isstandard(alg)
    standard_stepsize_controller!(integrator, alg)
  else # Default is PI-controller
    PI_stepsize_controller!(integrator, alg)
  end
end

@inline function predictive_step_accept_controller!(integrator, alg, q)
  if integrator.success_iter > 0
    expo = 1/(get_current_adaptive_order(integrator.alg,integrator.cache)+1)
    qgus=(integrator.dtacc/integrator.dt)*(((integrator.EEst^2)/integrator.erracc)^expo)
    qgus = max(inv(integrator.opts.qmax),min(inv(integrator.opts.qmin),qgus/integrator.opts.gamma))
    qacc=max(q,qgus)
  else
    qacc = q
  end
  if integrator.opts.qsteady_min <= qacc <= integrator.opts.qsteady_max
    qacc = one(qacc)
  end
  integrator.dtacc = integrator.dt
  integrator.erracc = max(1e-2,integrator.EEst)
  return integrator.dt/qacc
end

@inline function standard_step_accept_controller!(integrator, alg, q)
  if q <= integrator.opts.qsteady_max && q >= integrator.opts.qsteady_min
    q = one(q)
  end
  integrator.dt/q
end

@inline function PI_step_accept_controller!(integrator, alg, q)
  if q <= integrator.opts.qsteady_max && q >= integrator.opts.qsteady_min
    q = one(q)
  end
  integrator.qold = max(integrator.EEst,integrator.opts.qoldinit)
  return integrator.dt/q # dtnew
end

function step_accept_controller!(integrator, alg, q)
  if ispredictive(alg)
    predictive_step_accept_controller!(integrator, alg, q)
  elseif isstandard(alg)
    standard_step_accept_controller!(integrator, alg, q)
  else
    PI_step_accept_controller!(integrator, alg, q)
  end
end

function step_reject_controller!(integrator,alg)
  if ispredictive(alg)
    integrator.dt = integrator.success_iter == 0 ? 0.1integrator.dt : integrator.dt/integrator.qold
  elseif isstandard(alg)
    integrator.dt = integrator.qold
  else #PI
    integrator.dt = integrator.dt/min(inv(integrator.opts.qmin),integrator.q11/integrator.opts.gamma)
  end
  return
end

function stepsize_controller!(integrator, alg::JVODE)
  if iszero(integrator.EEst)
    η = integrator.opts.qmax
  else
    η = integrator.cache.η
    integrator.qold = η
  end
  η
end

function step_accept_controller!(integrator,alg::JVODE,η)
  q = inv(η)
  if q <= integrator.opts.qsteady_max && q >= integrator.opts.qsteady_min
    q = one(q)
  end
  return integrator.dt/q  # dtnew
end

function step_reject_controller!(integrator,alg::JVODE)
  integrator.dt *= integrator.qold
end

function stepsize_controller!(integrator, alg::QNDF)
  cnt = integrator.iter
  if cnt <= 3
    q = standard_stepsize_controller!(integrator, alg)
  else
    q = integrator.dt/integrator.cache.h
    integrator.qold = integrator.dt/q
  end
  q
end

@inline function stepsize_controller!(integrator,alg::Union{ExtrapolationMidpointDeuflhard,ImplicitDeuflhardExtrapolation})
  # Dummy function
  # ExtrapolationMidpointDeuflhard's stepsize scaling is stored in the cache;
  # it is computed by  stepsize_controller_internal! (in perfom_step!) resp. stepsize_predictor!
  # (in step_accept_controller! and step_reject_controller!)
  zero(typeof(integrator.opts.qmax))
end

function stepsize_controller_internal!(integrator,alg::Union{ExtrapolationMidpointDeuflhard,ImplicitDeuflhardExtrapolation})
  # Standard stepsize controller
  # Compute and save the stepsize scaling based on the latest error estimate of the current order
  if iszero(integrator.EEst)
    q = inv(integrator.opts.qmax)
  else
    # Update gamma and beta1
    integrator.opts.beta1 = typeof(integrator.opts.beta1)(1 // (2integrator.cache.n_curr + 1))
    integrator.opts.gamma = typeof(integrator.opts.gamma)(1 // 4)^integrator.opts.beta1
    # Compute new stepsize scaling
    qtmp = integrator.EEst^integrator.opts.beta1 / integrator.opts.gamma
    @fastmath q = max(inv(integrator.opts.qmax), min(inv(integrator.opts.qmin), qtmp))
  end
  integrator.cache.Q[integrator.cache.n_curr - alg.n_min + 1] = q
end

function stepsize_predictor!(integrator,alg::Union{ExtrapolationMidpointDeuflhard,ImplicitDeuflhardExtrapolation},n_new::Int64)
  # Compute and save the stepsize scaling for order n_new based on the latest error estimate of the current order.
  if iszero(integrator.EEst)
    q = inv(integrator.opts.qmax)
  else
    # Initialize
    @unpack t,EEst = integrator
    @unpack stage_number = integrator.cache
    tol = integrator.opts.internalnorm(integrator.opts.reltol,t) # Deuflhard's approach relies on EEstD ≈ ||relTol||
    s_curr = stage_number[integrator.cache.n_curr - alg.n_min + 1]
    s_new = stage_number[n_new - alg.n_min + 1]
    # Update gamma and beta1
    integrator.opts.beta1 = typeof(integrator.opts.beta1)(1 // (2integrator.cache.n_curr + 1))
    integrator.opts.gamma = typeof(integrator.opts.gamma)(1 // 4)^integrator.opts.beta1
    # Compute new stepsize scaling
    qtmp = (EEst * tol^(1.0 - s_curr / s_new))^integrator.opts.beta1 / integrator.opts.gamma
    @fastmath q = max(inv(integrator.opts.qmax),min(inv(integrator.opts.qmin),qtmp))
  end
  integrator.cache.Q[n_new - alg.n_min + 1] = q
end

function step_accept_controller!(integrator,alg::Union{ExtrapolationMidpointDeuflhard,ImplicitDeuflhardExtrapolation},q)
  # Compute new order and stepsize, return new stepsize
  @unpack n_min, n_max = alg
  @unpack n_curr, n_old, Q = integrator.cache
  s = integrator.cache.stage_number

  # Compute new order based on available quantities
  tmp = (n_min:n_curr) .- n_min .+ 1 # Index range of quantities computed so far
  dt_new = Vector{eltype(Q)}(undef,length(tmp)+1)
  dt_new[1:end-1] = integrator.dt ./ Q[tmp] # Store for the possible new stepsizes
  dt_new[1:end-1] = max.(abs(integrator.opts.dtmin), min.(abs(integrator.opts.dtmax), abs.(dt_new[1:end-1]))) # Safety scaling

  # n_new is the most efficient order of the last step
  work = s[tmp] ./ dt_new[1:end-1]
  n_new = argmin(work) + n_min - 1

  # Check if n_new may be increased
  if n_new == n_curr < min(n_max, n_old + 1) # cf. win_max in perfom_step! of the last step
    # Predict stepsize scaling for order (n_new + 1)
    stepsize_predictor!(integrator, alg, n_new+1) # Update cache.Q

    # Compute and scale the corresponding stepsize
    dt_new[end] = integrator.dt ./ Q[tmp[end]+1]
    dt_new[end] = max(abs(integrator.opts.dtmin), min(abs(integrator.opts.dtmax), abs.(dt_new[end])))

    # Check if (n_new  + 1) would have been more efficient than n_new
    if work[end] > s[tmp[end]+1] / dt_new[end]
      n_new = n_new + 1
    end
  end

  integrator.cache.n_curr = n_new
  dt_new[n_new - n_min + 1]
end

function step_reject_controller!(integrator, alg::Union{ExtrapolationMidpointDeuflhard,ImplicitDeuflhardExtrapolation})
  # Compute and save reduced stepsize dt_red of order n_old
  # Use the latest error estimate to predict dt_red if an estimate of order n_old is not available
  if integrator.cache.n_curr < integrator.cache.n_old
      stepsize_predictor!(integrator,alg,integrator.cache.n_old) # Update cache.Q
  end
  integrator.cache.n_curr = integrator.cache.n_old # Reset order for redoing the rejected step
  dt_red = integrator.dt / integrator.cache.Q[integrator.cache.n_old - integrator.alg.n_min + 1]
  dt_red = integrator.tdir*max(abs(integrator.opts.dtmin), min(abs(integrator.opts.dtmax), abs(dt_red))) # Safety scaling
  integrator.dt = dt_red
end

@inline function stepsize_controller!(integrator,alg::ExtrapolationMidpointHairerWanner)
  # Dummy function
  # ExtrapolationMidpointHairerWanner's stepsize scaling is stored in the cache;
  # it is computed by  stepsize_controller_internal! (in perfom_step!), step_accept_controller! or step_reject_controller!
  zero(typeof(integrator.opts.qmax))
end

function stepsize_controller_internal!(integrator,alg::ExtrapolationMidpointHairerWanner)
  # Standard stepsize controller
  # Compute and save the stepsize scaling based on the latest error estimate of the current order
  if iszero(integrator.EEst)
    q = inv(integrator.opts.qmax)
  else
    # Update gamma and beta1
    integrator.opts.beta1 = typeof(integrator.opts.beta1)(1 // (2integrator.cache.n_curr + 1))
    integrator.opts.gamma = typeof(integrator.opts.gamma)(65 // 100)^integrator.opts.beta1
    # Compute new stepsize scaling
    qtmp = integrator.EEst^integrator.opts.beta1 / integrator.opts.gamma
    @fastmath q = max(inv(integrator.opts.qmax), min(inv(integrator.opts.qmin), qtmp))
  end
  integrator.cache.Q[integrator.cache.n_curr + 1] = q
end

function step_accept_controller!(integrator,alg::ExtrapolationMidpointHairerWanner,q)
  # Compute new order and stepsize, return new stepsize
  @unpack n_min, n_max = alg
  @unpack n_curr, n_old, Q, sigma = integrator.cache
  s = integrator.cache.stage_number

  # Compute new order based on available quantities
  win_min_old = n_old - 1 # cf. win_min in perfom_step! of the last step
  tmp = win_min_old:(n_curr + 1) # Index range for the new order
  dt_new = fill(zero(eltype(Q)),n_max+1)
  dt_new[tmp] = integrator.dt ./ Q[tmp] # Store for the possible new stepsizes
  dt_new[tmp] = max.(abs(integrator.opts.dtmin), min.(abs(integrator.opts.dtmax), abs.(dt_new[tmp]))) # Safety scaling
  work= Vector{eltype(Q)}(undef,n_max+1) # work[n] is the work for order (n-1)
  work[tmp] = s[tmp] ./ dt_new[tmp]

  # Order selection
  n_new = n_old
  if n_curr == n_min # Enforce n_min + 1 ≦ n_new
    n_new = n_min + 1
  else
    if n_curr <= n_old
      if work[n_curr-1] < sigma * work[n_curr]
        n_new = max(n_curr-1,n_old-1,n_min+1) # Enforce n_min + 1≦ n_new
      elseif work[n_curr] < sigma * work[n_curr-1]
        n_new = min(n_curr+1,n_max-1) # Enforce n_new ≦ n_max - 1
      else
        n_new = n_curr # n_min + 1 ≦ n_curr
      end
    else
      if work[n_old] < sigma *  work[n_old+1]
        n_new = max(n_old-1,n_min+1)  # Enforce n_min + 1 ≦ n_new
      end
      if work[n_curr+1] <  sigma * work[n_new+1]
        n_new = min(n_new+1,n_max-1) # Enforce n_new ≦ n_max - 1
      end
    end
  end
  integrator.cache.n_curr = n_new

  # Stepsize selection
  if n_new == n_curr + 1
    # Compute the new stepsize of order n_new based on the optimal stepsize of order n_curr
    dt_new[n_new+1] = s[n_curr + 2]/s[n_curr + 1 ] * dt_new[n_curr+1]
    dt_new[n_new+1] = max(abs(integrator.opts.dtmin), min(abs(integrator.opts.dtmax), abs(dt_new[n_new+1])))
  end
  dt_new[n_new + 1]
end

function step_reject_controller!(integrator, alg::ExtrapolationMidpointHairerWanner)
  # Compute and save order and stepsize for redoing the current step
  @unpack n_old, n_curr, Q = integrator.cache

  # Order selection
  n_red = n_old
  if n_curr == n_old - 1
    n_red = max(alg.n_min+1,n_old-1) # Enforce n_min + 1 ≦ n_red
  end
  integrator.cache.n_curr = n_red

  # Stepsize selection
  dt_red = integrator.dt / Q[n_red + 1]
  dt_red = integrator.tdir*max(abs(integrator.opts.dtmin), min(abs(integrator.opts.dtmax), abs(dt_red))) # Safety scaling
  integrator.dt = dt_red
end


text:
# This file is adapted from Julia. License is MIT: https://julialang.org/license

function complete_local_path(s, i1, i2)
    cmp = REPL.REPLCompletions.complete_path(s, i2)
    completions = filter!(isdir, [REPL.REPLCompletions.completion_text(p) for p in cmp[1]])
    return completions, cmp[2], !isempty(completions)
end

#=function complete_installed_package(s, i1, i2, project_opt)
    pkgs = project_opt ? API.__installed(PKGMODE_PROJECT) : API.__installed()
    pkgs = sort!(collect(keys(filter((p) -> p[2] != nothing, pkgs))))
    cmp = filter(cmd -> startswith(cmd, s), pkgs)
    return cmp, i1:i2, !isempty(cmp)
end

function complete_remote_package(s, i1, i2)
    cmp = String[]
    julia_version = VERSION
    for reg in Types.registries(;clone_default=false)
        data = Types.read_registry(joinpath(reg, "Registry.toml"))
        for (uuid, pkginfo) in data["packages"]
            name = pkginfo["name"]
            if startswith(name, s)
                compat_data = Operations.load_package_data_raw(
                    VersionSpec, joinpath(reg, pkginfo["path"], "Compat.toml"))
                supported_julia_versions = VersionSpec(VersionRange[])
                for (ver_range, compats) in compat_data
                    for (compat, v) in compats
                        if compat == "julia"
                            union!(supported_julia_versions, VersionSpec(v))
                        end
                    end
                end
                if VERSION in supported_julia_versions
                    push!(cmp, name)
                end
            end
        end
    end
    return cmp, i1:i2, !isempty(cmp)
end=#

function complete_help(options, partial)
    names = String[]
    for cmds in values(super_specs)
         append!(names, [spec.canonical_name for spec in values(cmds)])
    end
    return sort!(unique!(append!(names, collect(keys(super_specs)))))
end

function complete_argument(to_complete, i1, i2, lastcommand, project_opt
                           )::Tuple{Vector{String},UnitRange{Int},Bool}
    if lastcommand == CMD_HELP
        completions = filter(x->startswith(x,to_complete), completion_cache.canonical_names)
        return completions, i1:i2, !isempty(completions)
    #=elseif lastcommand in [CMD_STATUS, CMD_RM, CMD_UP, CMD_TEST, CMD_BUILD, CMD_FREE, CMD_PIN]
        return complete_installed_package(to_complete, i1, i2, project_opt)=#
    elseif lastcommand in [CMD_ADD, CMD_DEVELOP]
        if occursin(Base.Filesystem.path_separator_re, to_complete)
            return complete_local_path(to_complete, i1, i2)
        #=else
            rps = complete_remote_package(to_complete, i1, i2)
            lps = complete_local_path(to_complete, i1, i2)
            return vcat(rps[1], lps[1]), isempty(rps[1]) ? lps[2] : i1:i2, length(rps[1]) + length(lps[1]) > 0=#
        end
    end
    return String[], 0:-1, false
end

function completions(full, index)::Tuple{Vector{String},UnitRange{Int},Bool}
    pre = full[1:index]
    if isempty(pre)
        return completion_cache.commands, 0:-1, false
    end
    x = parse(pre; for_completions=true)
    if x === nothing # failed parse (invalid command name)
        return String[], 0:-1, false
    end
    (key::Symbol, to_complete::String, spec, proj::Bool) = x
    last = split(pre, ' ', keepempty=true)[end]
    offset = isempty(last) ? index+1 : last.offset+1
    if last != to_complete # require a space before completing next field
        return String[], 0:-1, false
    end
    if key == :arg
        return complete_argument(to_complete, offset, index, spec.kind, proj)
    end
    possible::Vector{String} =
        key == :meta ? completion_cache.meta_options :
        key == :cmd ? completion_cache.commands :
        key == :sub ? completion_cache.subcommands[spec] :
        key == :opt ? completion_cache.options[spec.kind] :
        String[]
    completions = filter(x->startswith(x,to_complete), possible)
    return completions, offset:index, !isempty(completions)
end


text:
function abdiv(a::R, b::R) where {R <: AbstractFloat}
    if a == zero(R)
        return convert(R, Inf)
    else
        return b/a
    end
end

function solve_quadratic(a::R, bo2::R, c::R) where {R <: AbstractFloat}
    #bo2 = b over 2
    in_sqrt = bo2^2 - a*c
    if in_sqrt > zero(R)
        return (-bo2 + sqrt(in_sqrt))/a
    else
        return zero(R)
    end
end

function solve_piecewise_quadratic(α1::R, β1::R, κ2x::R) where {R <: AbstractFloat}
    #solve_quadratic(α1^2, -α1*β1, β1^2-κ2x) =
    (sqrt(κ2x)+β1)/α1
end

function solve_piecewise_quadratic(α1::R, α2::R, β1::R, β2::R, κ2x::R) where {R <: AbstractFloat}
    a = α1^2+α2^2
    bo2 = -α1*β1-α2*β2
    c = β1^2+β2^2-κ2x
    τ1t = solve_quadratic(a, bo2, c)
    ∇τ1t1 = α1*τ1t-β1
    ∇τ1t2 = α2*τ1t-β2
    if (∇τ1t1 <= 0) || (∇τ1t2 <= 0)
        tr1 = abdiv(α1, β1)
        tr2 = abdiv(α2, β2)
        if tr1 <= tr2
            τ1t = solve_piecewise_quadratic(α1, β1, κ2x)
        else
            τ1t = solve_piecewise_quadratic(α2, β2, κ2x)
        end
    end
    τ1t
end

function solve_piecewise_quadratic(α1::R, α2::R, α3::R, β1::R, β2::R, β3::R, κ2x::R) where {R <: AbstractFloat}
    a = α1^2+α2^2+α3^2
    bo2 = -α1*β1-α2*β2-α3*β3
    c = β1^2+β2^2+β3^2-κ2x
    τ1t = solve_quadratic(a, bo2, c)
    ∇τ1t1 = α1*τ1t-β1
    ∇τ1t2 = α2*τ1t-β2
    ∇τ1t3 = α3*τ1t-β3
    if (∇τ1t1 <= 0) || (∇τ1t2 <= 0) || (∇τ1t3 <= 0)
        tr1 = abdiv(α1, β1)
        tr2 = abdiv(α2, β2)
        tr3 = abdiv(α3, β3)
        if tr1 <= tr2
            if tr2 <= tr3
                τ1t = solve_piecewise_quadratic(α1, α2, β1, β2, κ2x)
            else
                τ1t = solve_piecewise_quadratic(α1, α3, β1, β3, κ2x)
            end
        else
            if tr1 <= tr3
                τ1t = solve_piecewise_quadratic(α2, α1, β2, β1, κ2x)
            else
                τ1t = solve_piecewise_quadratic(α2, α3, β2, β3, κ2x)
            end
        end
    end
    τ1t
end


text:
__precompile__()

module SQLite_Octo

import Octo

export        QueryBuilders
import .Octo: QueryBuilders

export                 @functions, and, or, not, like, between, asc, desc
import .QueryBuilders: @functions, and, or, not, like, between, asc, desc

export                 Select, Insert, Update, Delete, Join, Over
import .QueryBuilders: Select, Insert, Update, Delete, Join, Over

export                 Create, Alter, Drop, Database, Table, Column, Index, ForeignKey
import .QueryBuilders: Create, Alter, Drop, Database, Table, Column, Index, ForeignKey

export                 Queryable, arrange
import .QueryBuilders: Queryable, arrange

include("sql.jl")

export  SQL, table_from, DefaultArrangement
import .SQL: table_from, DefaultArrangement

end # module


text:
"""
Modifies the implementation of sinpi() and cospi() from IntervalArithmetic
to improve it so that sinpi(1) == 0
"""

using ValidatedNumerics

using ValidatedNumerics.IntervalArithmetic: atomic, SVector, @round

function find_quadrantspi(x::T) where {T}
    temp = IntervalArithmetic.atomic(Interval{T}, x) * 2

    return IntervalArithmetic.SVector(floor(temp.lo), floor(temp.hi))
end

sinpi(a...) = Base.Math.sinpi(a...)

function sinpi(a::Interval{T}) where T
    isempty(a) && return a

    whole_range = Interval{T}(-1, 1)

    diam(a) > 2 && return whole_range

    # The following is equiavlent to doing temp = a / half_pi  and
    # taking floor(a.lo), floor(a.hi)
    lo_quadrant = minimum(find_quadrantspi(a.lo))
    hi_quadrant = maximum(find_quadrantspi(a.hi))

    if hi_quadrant - lo_quadrant > 4  # close to limits
        return whole_range
    end

    lo_quadrant = mod(lo_quadrant, 4)
    hi_quadrant = mod(hi_quadrant, 4)

    # Different cases depending on the two quadrants:
    if lo_quadrant == hi_quadrant
        a.hi - a.lo > 1 && return whole_range  # in same quadrant but separated by almost 2pi
        lo = @round(sinpi(a.lo), sinpi(a.lo)) # Interval(sin(a.lo, RoundDown), sin(a.lo, RoundUp))
        hi = @round(sinpi(a.hi), sinpi(a.hi)) # Interval(sin(a.hi, RoundDown), sin(a.hi, RoundUp))
        return hull(lo, hi)

    elseif lo_quadrant==3 && hi_quadrant==0
        return @round(sinpi(a.lo), sinpi(a.hi)) # Interval(sin(a.lo, RoundDown), sin(a.hi, RoundUp))

    elseif lo_quadrant==1 && hi_quadrant==2
        return @round(sinpi(a.hi), sinpi(a.lo)) # Interval(sin(a.hi, RoundDown), sin(a.lo, RoundUp))

    elseif ( lo_quadrant == 0 || lo_quadrant==3 ) && ( hi_quadrant==1 || hi_quadrant==2 )
        return @round(min(sinpi(a.lo), sinpi(a.hi)), 1)
        # Interval(min(sin(a.lo, RoundDown), sin(a.hi, RoundDown)), one(T))

    elseif ( lo_quadrant == 1 || lo_quadrant==2 ) && ( hi_quadrant==3 || hi_quadrant==0 )
        return @round(-1, max(sinpi(a.lo), sinpi(a.hi)))
        # Interval(-one(T), max(sin(a.lo, RoundUp), sin(a.hi, RoundUp)))

    else #if( lo_quadrant == 0 && hi_quadrant==3 ) || ( lo_quadrant == 2 && hi_quadrant==1 )
        return whole_range
    end
end

cospi(a...) = Base.Math.cospi(a...)

function cospi(a::Interval{T}) where T
    isempty(a) && return a

    whole_range = Interval(-one(T), one(T))

    diam(a) > 2 && return whole_range

    lo_quadrant = minimum(find_quadrantspi(a.lo))
    hi_quadrant = maximum(find_quadrantspi(a.hi))

    if hi_quadrant - lo_quadrant > 4  # close to limits
        return Interval(-one(T), one(T))
    end

    lo_quadrant = mod(lo_quadrant, 4)
    hi_quadrant = mod(hi_quadrant, 4)

    # Different cases depending on the two quadrants:
    if lo_quadrant == hi_quadrant # Interval limits in the same quadrant
        a.hi - a.lo > 1 && return whole_range
        lo = @round(cospi(a.lo), cospi(a.lo))
        hi = @round(cospi(a.hi), cospi(a.hi))
        return hull(lo, hi)

    elseif lo_quadrant == 2 && hi_quadrant==3
        return @round(cospi(a.lo), cospi(a.hi))

    elseif lo_quadrant == 0 && hi_quadrant==1
        return @round(cospi(a.hi), cospi(a.lo))

    elseif ( lo_quadrant == 2 || lo_quadrant==3 ) && ( hi_quadrant==0 || hi_quadrant==1 )
        return @round(min(cospi(a.lo), cospi(a.hi)), 1)

    elseif ( lo_quadrant == 0 || lo_quadrant==1 ) && ( hi_quadrant==2 || hi_quadrant==3 )
        return @round(-1, max(cospi(a.lo), cospi(a.hi)))

    else#if ( lo_quadrant == 3 && hi_quadrant==2 ) || ( lo_quadrant == 1 && hi_quadrant==0 )
        return whole_range
    end
end

# These definitions were missing, too

using TaylorSeries: Taylor1

sinpi(x::Taylor1) = sin(x*π)
cospi(x::Taylor1) = cos(x*π)


text:
using ITensors, Test

import ITensors: nactive

@testset "QN" begin
  @testset "QNVal Basics" begin
    qv = ITensors.QNVal()
    @test !isactive(qv)

    qv = ITensors.QNVal("Sz", 0)
    @test ITensors.name(qv) == ITensors.SmallString("Sz")
    @test val(qv) == 0
    @test modulus(qv) == 1
    @test isactive(qv)

    qv = ITensors.QNVal("A", 1, 2)
    @test ITensors.name(qv) == ITensors.SmallString("A")
    @test val(qv) == 1
    @test modulus(qv) == 2
    @test !isfermionic(qv)

    qv = ITensors.QNVal("Nf", 1, -1)
    @test ITensors.name(qv) == ITensors.SmallString("Nf")
    @test val(qv) == 1
    @test modulus(qv) == -1
    @test isfermionic(qv)
  end

  @testset "QN Basics" begin
    q = QN()
    @test length(sprint(show, q)) > 1

    q = QN(("Sz", 1))
    @test length(sprint(show, q)) > 1
    @test isactive(q[1])
    @test val(q, "Sz") == 1

    q = QN("Sz", 1)
    @test length(sprint(show, q)) > 1
    @test isactive(q[1])
    @test val(q, "Sz") == 1

    q = QN("P", 1, 2)
    @test length(sprint(show, q)) > 1
    @test isactive(q[1])
    @test val(q, "P") == 1
    @test modulus(q, "P") == 2
    @test nactive(q) == 1

    q = QN(("A", 1), ("B", 2))
    @test isactive(q[1])
    @test isactive(q[2])
    @test val(q, "A") == 1
    @test val(q, "B") == 2
    @test modulus(q, "A") == 1
    @test modulus(q, "B") == 1

    q = QN(("B", 2), ("A", 1))
    @test val(q, "A") == 1
    @test val(q, "B") == 2
    @test nactive(q) == 2

    q = QN(("A", 1), ("B", 2), ("C", 3), ("D", 4))
    @test nactive(q) == 4

    @test_throws BoundsError begin
      q = QN(("A", 1), ("B", 2), ("C", 3), ("D", 4), ("E", 5))
    end
  end

  @testset "Comparison" begin
    @test QN() == QN()
    @test QN("A", 1) == QN("A", 1)
    @test QN(("A", 1), ("B", 3)) == QN(("A", 1), ("B", 3))
    @test QN(("A", 1), ("B", 3)) == QN(("B", 3), ("A", 1))

    # Zero value and missing sector treated the same:
    @test QN(("A", 0), ("B", 3)) == QN("B", 3)
    @test QN(("B", 3), ("A", 0)) == QN("B", 3)
  end

  @testset "Arithmetic" begin
    @test QN("Sz", 1) + QN() == QN("Sz", 1)
    @test QN("Sz", 1) + QN("Sz", 2) == QN("Sz", 3)
    @test QN("Sz", 1) + QN("Sz", -2) == QN("Sz", -1)

    @test QN(("A", 1), ("Sz", 0)) + QN(("A", 0), ("Sz", 1)) == QN(("A", 1), ("Sz", 1))

    @test QN("P", 0, 2) + QN("P", 1, 2) == QN("P", 1, 2)
    @test QN("P", 1, 2) + QN("P", 1, 2) == QN("P", 0, 2)

    # Arithmetic involving mixed-label QNs
    @test QN() - QN("Sz", 2) == QN("Sz", -2)
    @test QN("Sz", 2) - QN() == QN("Sz", 2)
    @test QN() - QN(("Sz", 2), ("N", 1)) == QN(("Sz", -2), ("N", -1))
    @test QN("N", 1) - QN("Sz", 2) == QN(("N", 1), ("Sz", -2))
  end

  @testset "Ordering" begin
    z = QN()
    qa = QN(("Sz", 1), ("Nf", 1))
    qb = QN(("Sz", 0), ("Nf", 2))
    qc = QN(("Sz", 1), ("Nf", 2))
    qd = QN(("Sz", 1), ("Nf", 2))
    qe = QN(("Sz", -1), ("Nf", -2))

    @test !(z < z)
    @test !(qa < z)
    @test (z < qa)
    @test (z < qb)
    @test !(qb < z)
    @test (z < qc)
    @test !(qc < z)
    @test (z < qd)
    @test !(qd < z)
    @test !(z < qe)
    @test (qe < z)

    @test !(qa > qb)
    @test qb > qa
    @test !(qb == qa)
    @test (qb < qc)
    @test !(qc < qb)
    @test !(qc == qb)
    @test (qc == qd)
    @test !(qc < qd)
    @test !(qd < qc)
  end

  @testset "Hashing" begin
    @test hash(QN(("Sz", 0))) == hash(QN())
    @test hash(QN("Sz", 0)) == hash(QN("N", 0))
    @test hash(QN(("Sz", 1), ("N", 2))) == hash(QN(("N", 2), ("Sz", 1)))
  end

  @testset "Negative value for mod > 1" begin
    @test QN("T", -1, 3) == QN("T", 2, 3)
    @test QN("T", -2, 3) == QN("T", 1, 3)
    @test QN("T", -3, 3) == QN("T", 0, 3)
  end
end

nothing


text:
## serializing values ##

# dummy types to tell number of bytes used to store length (4 or 1)
abstract LongSymbol
abstract LongTuple
abstract LongExpr
abstract UndefRefTag

const ser_version = 1 # do not make changes without bumping the version #!
const ser_tag = ObjectIdDict()
const deser_tag = ObjectIdDict()
let i = 2
    global ser_tag, deser_tag
    for t = {Symbol, Int8, Uint8, Int16, Uint16, Int32, Uint32,
             Int64, Uint64, Int128, Uint128, Float32, Float64, Char, Ptr,
             DataType, UnionType, Function,
             Tuple, Array, Expr, LongSymbol, LongTuple, LongExpr,
             LineNumberNode, SymbolNode, LabelNode, GotoNode,
             QuoteNode, TopNode, TypeVar, Box, LambdaStaticData,
             Module, UndefRefTag, Task, ASCIIString, UTF8String,
             :reserved6, :reserved7, :reserved8,
             :reserved9, :reserved10, :reserved11, :reserved12,

             (), Bool, Any, :Any, None, Top, Undef, Type,
             :Array, :TypeVar, :Box,
             :lambda, :body, :return, :call, symbol("::"),
             :(=), :null, :gotoifnot, :A, :B, :C, :M, :N, :T, :S, :X, :Y,
             :a, :b, :c, :d, :e, :f, :g, :h, :i, :j, :k, :l, :m, :n, :o,
             :p, :q, :r, :s, :t, :u, :v, :w, :x, :y, :z,
             :add_int, :sub_int, :mul_int, :add_float, :sub_float,
             :mul_float, :unbox, :box,
             :eq_int, :slt_int, :sle_int, :ne_int,
             :arrayset, :arrayref,
             :Core, :Base, :reserved15, :reserved16,
             :reserved17, :reserved18, :reserved19, :reserved20,
             false, true, nothing, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11,
             12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27,
             28, 29, 30, 31, 32}
        ser_tag[t] = int32(i)
        deser_tag[int32(i)] = t
        i += 1
    end
end

# tags >= this just represent themselves, their whole representation is 1 byte
const VALUE_TAGS = ser_tag[()]

writetag(s, x) = write(s, uint8(ser_tag[x]))

function write_as_tag(s, x)
    t = ser_tag[x]
    if t < VALUE_TAGS
        write(s, uint8(0))
    end
    write(s, uint8(t))
end

serialize(s, x::Bool) = write_as_tag(s, x)

serialize(s, ::Ptr) = error("cannot serialize a pointer")

serialize(s, ::()) = write_as_tag(s, ())

function serialize(s, t::Tuple)
    l = length(t)
    if l <= 255
        writetag(s, Tuple)
        write(s, uint8(l))
    else
        writetag(s, LongTuple)
        write(s, int32(l))
    end
    for i = 1:l
        serialize(s, t[i])
    end
end

function serialize(s, x::Symbol)
    if haskey(ser_tag, x)
        return write_as_tag(s, x)
    end
    pname = convert(Ptr{Uint8}, x)
    ln = int(ccall(:strlen, Csize_t, (Ptr{Uint8},), pname))
    if ln <= 255
        writetag(s, Symbol)
        write(s, uint8(ln))
    else
        writetag(s, LongSymbol)
        write(s, int32(ln))
    end
    write(s, pname, ln)
end

function serialize_array_data(s, a)
    elty = eltype(a)
    if elty === Bool && length(a)>0
        last = a[1]
        count = 1
        for i = 2:length(a)
            if a[i] != last || count == 127
                write(s, uint8((uint8(last)<<7) | count))
                last = a[i]
                count = 1
            else
                count += 1
            end
        end
        write(s, uint8((uint8(last)<<7) | count))
    else
        write(s, a)
    end
end

function serialize(s, a::Array)
    writetag(s, Array)
    elty = eltype(a)
    if elty !== Uint8
        serialize(s, elty)
    end
    if ndims(a) != 1
        serialize(s, size(a))
    else
        serialize(s, length(a))
    end
    if isbits(elty)
        serialize_array_data(s, a)
    else
        for i = 1:length(a)
            if isdefined(a, i)
                serialize(s, a[i])
            else
                writetag(s, UndefRefTag)
            end
        end
    end
end

function serialize{T,N,A<:Array}(s, a::SubArray{T,N,A})
    if !isbits(T) || stride(a,1)!=1
        return serialize(s, copy(a))
    end
    writetag(s, Array)
    serialize(s, T)
    serialize(s, size(a))
    serialize_array_data(s, a)
end

function serialize(s, e::Expr)
    l = length(e.args)
    if l <= 255
        writetag(s, Expr)
        write(s, uint8(l))
    else
        writetag(s, LongExpr)
        write(s, int32(l))
    end
    serialize(s, e.head)
    serialize(s, e.typ)
    for a = e.args
        serialize(s, a)
    end
end

function serialize(s, m::Module)
    writetag(s, Module)
    serialize(s, fullname(m))
end

function serialize(s, f::Function)
    writetag(s, Function)
    name = false
    if isgeneric(f)
        name = f.env.name
    elseif isa(f.env,Symbol)
        name = f.env
    end
    if isa(name,Symbol)
        if isdefined(Base,name) && is(f,eval(Base,name))
            write(s, uint8(0))
            serialize(s, name)
            return
        end
        mod = ()
        if isa(f.env,Symbol)
            mod = Core
        elseif !is(f.env.defs, ())
            mod = f.env.defs.func.code.module
        end
        if mod !== ()
            if isdefined(mod,name) && is(f,eval(mod,name))
                # toplevel named func
                write(s, uint8(2))
                serialize(s, mod)
                serialize(s, name)
                return
            end
        end
        write(s, uint8(3))
        serialize(s, f.env)
    else
        linfo = f.code
        @assert isa(linfo,LambdaStaticData)
        write(s, uint8(1))
        serialize(s, f.env)
        serialize(s, linfo)
    end
end

const lambda_numbers = WeakKeyDict()
lnumber_salt = 0
function lambda_number(l::LambdaStaticData)
    global lnumber_salt, lambda_numbers
    if haskey(lambda_numbers, l)
        return lambda_numbers[l]
    end
    # a hash function that always gives the same number to the same
    # object on the same machine, and is unique over all machines.
    ln = hash(lnumber_salt+(uint64(myid())<<44))
    lnumber_salt += 1
    lambda_numbers[l] = ln
    return ln
end

function serialize(s, linfo::LambdaStaticData)
    writetag(s, LambdaStaticData)
    serialize(s, lambda_number(linfo))
    serialize(s, uncompressed_ast(linfo))
    if isdefined(linfo.def, :roots)
        serialize(s, linfo.def.roots)
    else
        serialize(s, {})
    end
    serialize(s, linfo.sparams)
    serialize(s, linfo.inferred)
    serialize(s, linfo.module)
    if isdefined(linfo, :capt)
        serialize(s, linfo.capt)
    else
        serialize(s, nothing)
    end
end

function serialize(s, t::Task)
    if istaskstarted(t) && !istaskdone(t)
        error("cannot serialize a running Task")
    end
    writetag(s, Task)
    serialize(s, t.code)
    serialize(s, t.storage)
    serialize(s, t.state == :queued || t.state == :waiting ? (:runnable) : t.state)
    serialize(s, t.result)
    serialize(s, t.exception)
end

function serialize_type_data(s, t)
    tname = t.name.name
    serialize(s, tname)
    mod = t.name.module
    serialize(s, mod)
    if isdefined(mod,tname) && is(t,eval(mod,tname))
        serialize(s, ())
    else
        serialize(s, t.parameters)
    end
end

function serialize(s, t::DataType)
    if haskey(ser_tag,t)
        write_as_tag(s, t)
    else
        writetag(s, DataType)
        write(s, uint8(0))
        serialize_type_data(s, t)
    end
end

function serialize_type(s, t::DataType)
    if haskey(ser_tag,t)
        writetag(s, t)
    else
        writetag(s, DataType)
        write(s, uint8(1))
        serialize_type_data(s, t)
    end
end

function serialize(s, x)
    if haskey(ser_tag,x)
        return write_as_tag(s, x)
    end
    t = typeof(x)
    serialize_type(s, t)
    if length(t.names)==0 && t.size>0
        write(s, x)
    else
        serialize(s, length(t.names))
        for i in 1:length(t.names)
            if isdefined(x, i)
                serialize(s, getfield(x, i))
            else
                writetag(s, UndefRefTag)
            end
        end
    end
end

## deserializing values ##

function deserialize(s)
    handle_deserialize(s, int32(read(s, Uint8)))
end

function handle_deserialize(s, b)
    if b == 0
        return deser_tag[int32(read(s, Uint8))]
    end
    tag = deser_tag[b]
    if b >= VALUE_TAGS
        return tag
    elseif is(tag,Tuple)
        len = int32(read(s, Uint8))
        return deserialize_tuple(s, len)
    elseif is(tag,LongTuple)
        len = read(s, Int32)
        return deserialize_tuple(s, len)
    end
    return deserialize(s, tag)
end

deserialize_tuple(s, len) = ntuple(len, i->deserialize(s))

deserialize(s, ::Type{Symbol}) = symbol(read(s, Uint8, int32(read(s, Uint8))))
deserialize(s, ::Type{LongSymbol}) = symbol(read(s, Uint8, read(s, Int32)))

function deserialize(s, ::Type{Module})
    path = deserialize(s)
    m = Main
    for mname in path
        if !isdefined(m,mname)
            warn("Module $mname not defined on process $(myid())")  # an error seemingly fails
        end
        m = eval(m,mname)::Module
    end
    m
end

const known_lambda_data = Dict()

function deserialize(s, ::Type{Function})
    b = read(s, Uint8)
    if b==0
        name = deserialize(s)::Symbol
        if !isdefined(Base,name)
            return (args...)->error("function $name not defined on process $(myid())")
        end
        return eval(Base,name)::Function
    elseif b==2
        mod = deserialize(s)::Module
        name = deserialize(s)::Symbol
        if !isdefined(mod,name)
            return (args...)->error("function $name not defined on process $(myid())")
        end
        return eval(mod,name)::Function
    elseif b==3
        env = deserialize(s)
        return ccall(:jl_new_gf_internal, Any, (Any,), env)::Function
    end
    env = deserialize(s)
    linfo = deserialize(s)
    ccall(:jl_new_closure, Any, (Ptr{Void}, Any, Any),
          C_NULL, env, linfo)::Function
end

function deserialize(s, ::Type{LambdaStaticData})
    lnumber = deserialize(s)
    ast = deserialize(s)
    roots = deserialize(s)
    sparams = deserialize(s)
    infr = deserialize(s)
    mod = deserialize(s)
    capt = deserialize(s)
    if haskey(known_lambda_data, lnumber)
        return known_lambda_data[lnumber]
    else
        linfo = ccall(:jl_new_lambda_info, Any, (Any, Any), ast, sparams)
        linfo.inferred = infr
        linfo.module = mod
        linfo.roots = roots
        if !is(capt,nothing)
            linfo.capt = capt
        end
        known_lambda_data[lnumber] = linfo
        return linfo
    end
end

function deserialize(s, ::Type{Array})
    d1 = deserialize(s)
    if isa(d1,Type)
        elty = d1
        d1 = deserialize(s)
    else
        elty = Uint8
    end
    if isa(d1,Integer)
        dims = (int(d1),)
    else
        dims = convert(Dims, d1)::Dims
    end
    if isbits(elty)
        n = prod(dims)::Int
        if elty === Bool && n>0
            A = Array(Bool, dims)
            i = 1
            while i <= n
                b = read(s, Uint8)
                v = bool(b>>7)
                count = b&0x7f
                nxt = i+count
                while i < nxt
                    A[i] = v; i+=1
                end
            end
            return A
        else
            return read(s, elty, dims)
        end
    end
    A = Array(elty, dims)
    for i = 1:length(A)
        tag = int32(read(s, Uint8))
        if tag==0 || !is(deser_tag[tag], UndefRefTag)
            A[i] = handle_deserialize(s, tag)
        end
    end
    return A
end

deserialize(s, ::Type{Expr})     = deserialize_expr(s, int32(read(s, Uint8)))
deserialize(s, ::Type{LongExpr}) = deserialize_expr(s, read(s, Int32))

function deserialize_expr(s, len)
    hd = deserialize(s)::Symbol
    ty = deserialize(s)
    e = Expr(hd)
    e.args = { deserialize(s) for i=1:len }
    e.typ = ty
    e
end

function deserialize(s, ::Type{TypeVar})
    nf_expected = deserialize(s)
    name = deserialize(s)
    lb = deserialize(s)
    ub = deserialize(s)
    TypeVar(name, lb, ub)
end

function deserialize(s, ::Type{UnionType})
    nf_expected = deserialize(s)
    types = deserialize(s)
    Union(types...)
end

function deserialize(s, ::Type{DataType})
    form = read(s, Uint8)
    name = deserialize(s)::Symbol
    mod = deserialize(s)::Module
    params = deserialize(s)
    ty = eval(mod,name)
    if is(params,())
        t = ty
    else
        t = apply_type(ty, params...)
    end
    if form == 0
        return t
    end
    deserialize(s, t)
end

deserialize{T}(s, ::Type{Ptr{T}}) = pointer(T, 0)

function deserialize(s, ::Type{Task})
    t = Task(deserialize(s))
    t.storage = deserialize(s)
    t.state = deserialize(s)
    t.result = deserialize(s)
    t.exception = deserialize(s)
    t
end

# default DataType deserializer
function deserialize(s, t::DataType)
    if length(t.names)==0 && t.size>0
        # bits type
        return read(s, t)
    end
    nf_expected = deserialize(s)
    nf = length(t.names)
    if nf == 0
        return ccall(:jl_new_struct, Any, (Any,Any...), t)
    elseif !t.mutable
        if nf == 1
            return ccall(:jl_new_struct, Any, (Any,Any...), t, deserialize(s))
        elseif nf == 2
            f1 = deserialize(s)
            f2 = deserialize(s)
            return ccall(:jl_new_struct, Any, (Any,Any...), t, f1, f2)
        elseif nf == 3
            f1 = deserialize(s)
            f2 = deserialize(s)
            f3 = deserialize(s)
            return ccall(:jl_new_struct, Any, (Any,Any...), t, f1, f2, f3)
        else
            flds = { deserialize(s) for i = 1:nf }
            return ccall(:jl_new_structv, Any, (Any,Ptr{Void},Uint32),
                         t, flds, nf)
        end
    else
        x = ccall(:jl_new_struct_uninit, Any, (Any,), t)
        for i in 1:length(t.names)
            tag = int32(read(s, Uint8))
            if tag==0 || !is(deser_tag[tag], UndefRefTag)
                setfield!(x, i, handle_deserialize(s, tag))
            end
        end
        return x
    end
end

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
using TypeTraits
using Base.Test

abstract Player

@trait type Entity
  id::Int64
  x::Float64
  y::Float64
end

@mixin type User(Entity) <: Player
  name::String
end

function User(name::String)
  User(0, 0, 0, name)
end

user = User("Foobar")
println(user)


text:
# PREAMBLE

# PKG_SETUP

# ## Setup

using ContactImplicitMPC
using LinearAlgebra

# ## Simulation
s_sim = get_simulation("quadruped", "sine1_2D_lc", "sinusoidal");
s = get_simulation("quadruped", "flat_2D_lc", "flat");
model = s.model
env = s.env
nq = model.nq
s_sim.model.μ_world
s.model.μ_world

# ## Reference Trajectory
ref_traj = deepcopy(get_trajectory(s.model, s.env,
    joinpath(module_dir(), "src/dynamics/quadruped/gaits/gait2.jld2"),
    load_type = :split_traj_alt));

H = ref_traj.H
h = ref_traj.h

# ## MPC setup
N_sample = 5
H_mpc = 10
h_sim = h / N_sample
H_sim = 1000 #4000 #5000
κ_mpc = 1.0e-4

obj = TrackingObjective(model, env, H_mpc,
    q = [Diagonal(1e-2 * [10; 0.02; 0.25; 0.25 * ones(nq-3)]) for t = 1:H_mpc],
    u = [Diagonal(3e-2 * ones(model.nu)) for t = 1:H_mpc],
    γ = [Diagonal(1.0e-100 * ones(model.nc)) for t = 1:H_mpc],
    b = [Diagonal(1.0e-100 * ones(model.nc * friction_dim(env))) for t = 1:H_mpc]);

p = ci_mpc_policy(ref_traj, s, obj,
    H_mpc = H_mpc,
    N_sample = N_sample,
    κ_mpc = κ_mpc,
    n_opts = NewtonOptions(
        r_tol = 3e-4,
        max_iter = 5,
        ),
    mpc_opts = CIMPCOptions(
        altitude_update = true,
        altitude_impact_threshold = 0.05,
        altitude_verbose = true,
        )
    );

# ## Initial conditions
q1_sim, v1_sim = initial_conditions(ref_traj);

# ## Simulator
sim = simulator(s, H_sim, h=h_sim, policy=p);

# ## Simulate
simulate!(sim, q1_sim, v1_sim);
@benchmark simulate!($sim, $q1_sim, $v1_sim)

# ## Visualizer
vis = ContactImplicitMPC.Visualizer()
ContactImplicitMPC.render(vis)

# ## Visualize
plot_surface!(vis, s_sim.env, ylims=[0.3, -0.05])
anim = visualize_meshrobot!(vis, model, sim.traj, sample=5);

# ## Timing result
# Julia is [JIT-ed](https://en.wikipedia.org/wiki/Just-in-time_compilation) so re-run the MPC setup through Simulate for correct timing results.
process!(sim) # Time budget
H_sim * h_sim / sum(sim.stats.dt) # Speed ratio


text:
import ChainRulesTestUtils.test_rrule
import ChainRulesCore: rrule, unthunk
import Yota: isprimitive, CR_CTX

const broadcasted = Broadcast.broadcasted

double_inc(x::Number) = 2x + 1
double_inc(x::AbstractArray) = 2x .+ 1

double_dec(x::Number) = 2x - 1

primitive_test(x; y=1) = x + y
primitive_test(x, y) = x + y
primitive_test2(x; y=1) = x + y
primitive_test2(x, y) = x + y

# rrule(::typeof(primitive_test), x) = primitive_test(x), dy -> (NoTangent(), 1)
rrule(::typeof(primitive_test), x; y=1) = primitive_test(x; y=y), dy -> (NoTangent(), 1)
# rrule(::YotaRuleConfig, ::typeof(primitive_test2), x) = primitive_test2(x), dy -> (NoTangent(), 1)
rrule(::YotaRuleConfig, ::typeof(primitive_test2), x; y=1) = primitive_test2(x; y=y), dy -> (NoTangent(), 1)

@testset "chainrules api" begin
    config = YotaRuleConfig()

    rr = make_rrule(double_inc, 2.0)
    val, pb = rr(config, double_inc, 3.0)
    @test val == 7
    @test pb(1.0) == (ZeroTangent(), 2.0)

    trace(broadcasted, double_dec, [1.0, 2.0])

    rr = make_rrule(broadcasted, double_dec, [1.0, 2.0])
    val, pb = rr(config, broadcasted, double_dec, [3.0, 4.0])
    @test val == [5.0, 7.0]
    @test pb([1, 1]) == (ZeroTangent(), ZeroTangent(), [2.0, 2.0])

    val, pb = rrule_via_ad(config, double_inc, 3.0)
    @test val == 7
    @test pb(1.0) == (ZeroTangent(), 2.0)

    x = rand(3)
    val, pb = rrule_via_ad(config, double_inc, x)
    @test val == double_inc(x)
    dxs = map(unthunk, pb(ones(3)))
    @test dxs == (ZeroTangent(), [2.0, 2.0, 2.0])
    dxs = map(unthunk, pb([1, 2, 3]))
    @test dxs == (ZeroTangent(), [2.0, 4.0, 6.0])

    x = rand(3)
    val, pb = rrule_via_ad(config, broadcasted, double_dec, x)
    @test val == broadcast(double_dec, x)
    dxs = map(unthunk, pb(ones(3)))
    @test dxs == (ZeroTangent(), ZeroTangent(), [2.0, 2.0, 2.0])
    dxs = map(unthunk, pb([1, 2, 3]))
    @test dxs == (ZeroTangent(), ZeroTangent(), [2.0, 4.0, 6.0])

    x, y = rand(2)
    @test isprimitive(CR_CTX, primitive_test, x) == true
    @test isprimitive(CR_CTX, Core.kwfunc(primitive_test), (y=1,), primitive_test, x) == true
    @test isprimitive(CR_CTX, primitive_test, x, y) == false

    @test isprimitive(CR_CTX, primitive_test2, x) == true
    @test isprimitive(CR_CTX, Core.kwfunc(primitive_test2), (y=1,), primitive_test, x,) == true
    @test isprimitive(CR_CTX, primitive_test2, x, y) == false
end


text:
println("Ab) Read in a full multiplet representation from files compatible to Ratip2012 and Grasp2018.")
#
wa = JAC.ManyElectron.Multiplet("from Ratip2012", "../test/files/ratip2012-at-is-a1-csl.inp", "../test/files/ratip2012-at-is-a1-scf.out",
                                                  "../test/files/ratip2012-at-is-a1-relci.mix")

println("Not yet implemented; files do not exist !")
wa = JAC.ManyElectron.Multiplet("from Grasp2018", "../test/files/grasp2018-at-is-a1-csl.inp", "../test/files/grasp2018-at-is-a1-scf.out",
                                                  "../test/files/grasp2018-at-is-a1-rci.mix")


text:
# 6D Vlasov - Poisson Single Species

include("pif_tools.jl")
importall pif

using Sobol

DTYPE=Float64

Np=Int(1e4) #Number of particles
Ns=1 #Number of species
dt=0.05
tmax=20
gamma=1.0 # Gravity <0
q=-1.0 # Electrons
m=1.0

Nx1=1
Nx2=1
Nx3=1

# Default parameters
v01=0;v02=0;eps=0;delta=0;epsi=0; betai=0; betar=0;
sigma1=1;sigma2=1;sigma3=1;

k=0.5; eps=0.1;

L1=DTYPE(2*pi/k)
L2=DTYPE(2*pi/k)
L3=DTYPE(2*pi/k)

# Particles
x1n=zeros(DTYPE,Np)
x2n=zeros(DTYPE,Np)
x3n=zeros(DTYPE,Np)

v1n=zeros(DTYPE,Np)
v2n=zeros(DTYPE,Np)
v3n=zeros(DTYPE,Np)
#qn=q*ones(Float64,Np) # Charge
#mn=ones(Float64,Np) # Mass

fn=zeros(DTYPE,Np)
gn=zeros(DTYPE,Np)

# Initial condition
f0(x1,x2,x3,v1,v2,v3)=((1+ eps*sin.(k*x1+pi/5)).*
    exp.(- (v1.^2./sigma1^2 + v3.^2/sigma3^2)/2 ).*
    (delta*exp.(-(v2-v01).^2/2/sigma2.^2) +
    (1-delta)*exp.( - (v2-v02).^2/2/sigma2^2))
    *(2*pi)^(-1.5))/sigma3/sigma1/sigma2

g0(x1,x2,x3,v1,v2,v3)=(exp.(-v1.^2./sigma1^2/2 -v3.^2./sigma1^2/2).*
    (delta*exp.(-(v2-v01).^2/2/sigma2.^2) +
    (1-delta)*exp.( - (v2-v02).^2/2/sigma2^2))
    *(2*pi)^(-1.5) )/sigma3/sigma1/sigma2/L1/L2/L3;
# Fill with uniform numbers
using Sobol
sob = SobolSeq(6)
skip(sob, 22) # Skip some entries
for i=1:Np
  x1n[i],x2n[i],x3n[i],v1n[i],v2n[i],v3n[i]=next(sob)
end

# Sampling
using SpecialFunctions

norminv(p)=sqrt(2).*erfinv.(2 *p - 1.0)
v1n=norminv(v1n)*sigma1
v2n=norminv(v2n)*sigma2;
v3n=norminv(v3n)*sigma3;
x1n=x1n*L1;x2n=x2n*L2;x3n=x3n*L3;
qn=q*ones(DTYPE,Np)
mn=ones(DTYPE,Np)*m

# Likelihoods
fn=f0(x1n,x2n,x3n,v1n,v2n,v3n)
gn=g0(x1n,x2n,x3n,v1n,v2n,v3n)

#using PyPlot
# figure()
#h = plt[:hist](v1n,v40)

# Coefficients for Runge Kutta
const rksd=[2/3.0, -2/3.0, 1  ]
const rksc=[ 7/24.0, 3/4.0, -1/24.0]

#Particle mesh coupling
pm=pm_pif3d(Nx1,Nx2,Nx3,L1,L2,L3)

print( pm.N1*pm.N2*pm.N3 )
# Fields
rho=Array{Complex{DTYPE}}(pm.N1,pm.N2,pm.N3)
Phi=Array{Complex{DTYPE}}(pm.N1,pm.N2,pm.N3)
E1=Array(Phi);E2=Array(Phi);E3=Array(Phi)
E1n=zeros(Complex{DTYPE},Np)
E2n=zeros(Complex{DTYPE},Np)
E3n=zeros(Complex{DTYPE},Np)
wn=fn./gn

Nt=Int(ceil(tmax/dt))
Epot=zeros(DTYPE,Nt,3)
Ekin=zeros(DTYPE,Nt,3)
Momentum=zeros(DTYPE,Nt,3)
kineticenergy=zeros(DTYPE,Nt)
ttime=collect((0:Nt-1)*dt)
tic()
for tdx=1:Nt

  @inbounds for rkdx=1:length(rksd)
    rho=accum_osde(pm,x1n,x2n,x3n, wn.*qn)./Np

    solve_poisson!(pm,rho*gamma,Phi)

    E1,E2,E3=gradient(pm,-Phi)

    if rkdx==1

    #fieldenergy[tdx]=H1seminorm(pm,Phi)/2.
    #kineticenergy[tdx]=sum( (v1n.^2 .+v2n.^2.+v3n.^2).*wn)/Np/2.
    Epot[tdx,1]=L2norm(pm,E1)/2.
    Epot[tdx,2]=L2norm(pm,E2)/2.
    Epot[tdx,3]=L2norm(pm,E3)/2.
    Ekin[tdx,1]=sum(v1n.^2.*wn)/Np/2.
    Ekin[tdx,2]=sum(v2n.^2.*wn)/Np/2.
    Ekin[tdx,3]=sum(v3n.^2.*wn)/Np/2.
    Momentum[tdx,1]=sum(v1n.*wn)/Np
    Momentum[tdx,2]=sum(v2n.*wn)/Np
    Momentum[tdx,3]=sum(v3n.*wn)/Np
    end

    # eval_vectorfield(pm,x1n,x2n,x3n,E1,E2,E3,E1n,E2n,E3n)
    #eval_scalar(pm,x1n,x2n,x3n,[E1,E2,E3],[E1n,E2n,E3n])
    #  v1n.+=(dt*rksc[rkdx]).*real(E1n).*(qn./mn)
    #  v2n.+=(dt*rksc[rkdx])*real(E2n).*(qn./mn)
    #  v3n.+=(dt*rksc[rkdx])*real(E3n).*(qn./mn)

 integrate_H_E(pm,E1,E2,E3,x1n,x2n,x3n,v1n,v2n,v3n,qn,mn,DTYPE(dt*rksc[rkdx]))
    #  eval_scalar(pm,x1n,x2n,x3n,E1,En)
    # v1n.+=dt*rksc[rkdx]*real(En).*(qn./mn)
    # eval_scalar(pm,x1n,x2n,x3n,E2,En)
    # v2n.+=dt*rksc[rkdx]*real(En).*(qn./mn)
    # eval_scalar(pm,x1n,x2n,x3n,E3,En)
    # v3n.+=dt*rksc[rkdx]*real(En).*(qn./mn)

    # v1n+=dt*rksc[rkdx]*eval_scalar(pm,x1n,x2n,x3n,E1*q/m)
    # v2n+=dt*rksc[rkdx]*eval_scalar(pm,x1n,x2n,x3n,E2*q/m)
    # v3n+=dt*rksc[rkdx]*eval_scalar(pm,x1n,x2n,x3n,E3*q/m)

    x1n.+=dt*rksd[rkdx].*v1n
    x2n.+=dt*rksd[rkdx].*v2n
    x3n.+=dt*rksd[rkdx].*v3n

  end
  x1n=mod(x1n,L1)
  x2n=mod(x2n,L2)
  x3n=mod(x3n,L3)

end
toc()

using PyPlot
# semilogy(ttime,Epot)

energy=gamma*sum(Epot,2)+sum(Ekin,2)
moment_error=sum(abs(Momentum[:,:].-reshape(Momentum[1,:],1,3)),2);

# E1n_=zeros(Float64,Np)
# using BenchmarkTools
#
# @benchmark eval_scalar(pm,x1n,x2n,x3n,E1*q/m,E1n)
# @benchmark eval_scalar2(pm,x1n,x2n,x3n,E1*q/m,E1n_)
# print((E1n_[1]-real(E1n[1]))/real(E1n[1]) )
figure()
semilogy(ttime,abs.((energy-energy[1])./energy[1]))
grid()
ylabel("rel. energy error")
xlabel("time ")
# savefig("$prefix\_energy_error.png")
figure()
semilogy(ttime,moment_error)
xlabel("time "); grid()
ylabel("absolute momentum error")
# savefig("$prefix\_momentum_error.png")

#
figure()
 semilogy(ttime,Epot)

# figure()
# plot(energy)
#
#
# figure()
# semilogy(moment_error)
#
# grid()
#
#
# figure()
# plot(abs(Phi[:]))
# # using PyPlot
# # (L2norm(pm,E1)+L2norm(pm,E2)+L2norm(pm,E3))/2
# #
# # fieldenergy[end]
# Nviz=Np
# figure()
# sc=scatter3D(x1n[1:Nviz],x2n[1:Nviz],x3n[1:Nviz],s=2,c=fn[1:Nviz],depthshade="True",lw=0)
# colorbar(sc)
# tight()
# plot(fieldenergy)
#
# scatter(x1n,x2n,s=3,lw=0)
#
# figure()
# semilogy(ttime,fieldenergy)
# fig=figure()
# ax = gca(projection="3d")
#
# ax[:scatter](x1n,x2n,x3n)
# @time eval_scalar(pm,x1n,x2n,x3n,E1*q/m,En)

# figure()
# plot(x1n,x2n,x3n)

#@time rho=rhs_particle(pm,x1n,x2n,x3n, wn)./Np*q

#Damping rate for linear Landau
#omega=1.415661888604536 - 0.153359466909605im;
#semilogy(ttime[ttime.<10.],0.5*fieldenergy[1]*abs(exp(-im*omega*(ttime[ttime.<10.]-0.4))).^2)

#plot(abs(rho)[:])

#rho

# plot(abs(Phi[:]))

# type vlasov6d
#
#   Np::Int #Number of particles
#   Ns::Int #Number of species
#
#   x1n::Array{Float64,1}
#   x2n::Array{Float64,1}
#   x3n::Array{Float64,1}
#
#   v1n::Array{Float64,1}
#   v2n::Array{Float64,1}
#   v3n::Array{Float64,1}
#   qn::Array{Float64,1} # Charge
#   mn::Array{Float64,1} # Mass
#
#   fn::Array{Float64,1}
#   gn::Array{Float64,1}
#
#   function vlasov6d(Np,Ns)
#     x=zeros(Float64,Int(Np)*Int(Ns))
#     new(Int(Np),Int(Ns),x,x,x, x,x,x, x,x,x,x)
#
#   end
# end
#
#
# function std_Maxwellian(x1,x2,x3,v1,v2,v3)
#   return exp(-(v1.^2+v2.^2+v3.^2)/2 )/(sqrt(2*pi))^3
# end
#
#
#
# function sample_std_Maxwellian(sim::vlasov6d, L1,L2,L3   )
#   sim.x1n=rand(Float64,sim.Np)*L1
#   sim.x2n=rand(Float64,sim.Np)*L2
#   sim.x3n=rand(Float64,sim.Np)*L3
#
#   sim.v1n=randn(Float64,sim.Np)
#   sim.v2n=randn(Float64,sim.Np)
#   sim.v3n=randn(Float64,sim.Np)
#
#   gn=std_Maxwellian(sim.x1n,sim.x2n,sim.x3n,sim.v1n,sim.v2n,sim.v3n)/L1/L2/L3
# end
#
#
# f0=(x1,x2,x3,v1,v2,v3)->(exp(-(v1.^2+v2.^2+v3.^2)/2 )/(sqrt(2*pi))^3
#
#
# #Electron only
# sim=vlasov6d(1e5,1)
# k=0.5;
# epse=0.1
# #Particle mesh coupling
# pm=pif3d(2,2,2,2*pi/k)
#
# sim.qn[:]=1.;sim.mn[:]=1.
# sample_std_Maxwellian(sim,pm.L1,pm.L2,pm.L3)
# sim.fn=(1+epse*cos(k*x)  std_Maxwellian(sim.x1n,sim.x2n,sim.x3n,sim.v1n,sim.v2n,sim.v3n)
#
# rhs=rhs.particle(sim.x1n,sim.x2n,sim.x3n, sim.fn/sim.gn) /Np*Ns
# sim.x1n

max_stars_repo_name: ameresj/VlasovPIF

max_issues_repo_name: ameresj/VlasovPIF

max_forks_repo_name: ameresj/VlasovPIF


text:
"""
hadamarddev(data, rate; [frequency=false], [overlapping=true], [taus=Octave])
Calculates the hadamard deviation

#parameters:
* `<data>`:			The data array to calculate the deviation from either as as phases or frequencies.
* `<rate>`:			The rate of the data given.
* `[frequency]`:		True if `data` contains frequency data otherwise (default) phase data is assumed.
* `[overlapping]`:	True (default) to calculate overlapping deviation, false otherwise.
* `[taus]`:			Taus to calculate the deviation at. This can either be an AllanTauDescriptor type `AllTaus, Decadade, HalfDecade, Octave, HalfOctave, QuarterOctave`, an array of taus to calculate at, a float number to build a custom log-scale on or an integer to build a specific number of log spaced points.

#returns: named tupple (tau, deviation, error, count)
* `tau`:		Taus which where used.
* `deviation`:	Deviations calculated.
* `error`:		Respective errors.
* `count`:		Number of contributing terms for each deviation.
"""
function hadamarddev(
		data::Array{T, 1},
		rate::AbstractFloat;
		frequency::Bool = false,
		overlapping::Bool = true,
		taus::Union{Type{U}, Integer, AbstractFloat, Array{Float64}} = 192) where {T, U <: AllanTauDescriptor}

	#frequency to phase calculation
	if frequency
		data = frequencytophase(data, rate)
	end

	n = length(data)
	if n < 5
		error("Length for `data` in hadamarddev must be at least 5 or greater")
	end

	#tau calculations
	m = taudescription_to_m(taus, rate, n)

	dev = zeros(T, length(m)) #hadamarddeviation
	deverr = zeros(T, length(m)) #hadamarddeviation error
	devcount = zeros(Int, length(m)) #sum term count

	mStride = 1 #overlapping - can be overwritten in loop for consecutive
	@inbounds for (index, τ) in enumerate(m)

		if !overlapping #overwrite stride for consecutive operation
			mStride = τ
		end

		#hadamard deviation: http://www.leapsecond.com/tools/adev_lib.c
		sum = zero(T)
		i = 1
		terms = 0
		while (i + 3 * τ) <= n
			v = data[i + 3 * τ] - 3 * data[i + 2 * τ] + 3 * data[i + τ] - data[i]
			sum += v * v
			i += mStride
			terms += 1
		end
		if terms <= 1 #break the tau loop if no contribution with term-count > 1 is done
			break
		end

		dev[index] = sqrt(sum / (6 * terms)) / τ * rate
		deverr[index] = dev[index] / sqrt(terms)
		devcount[index] = terms
	end

	selector = devcount .> 1 #select only entries, where 2 or more terms contributed to the deviation
	(tau = m[selector] ./ rate, deviation = dev[selector], error = deverr[selector], count = devcount[selector])
end


text:
export vectorizePattern, SamplingOp

"""
 idx contains sampling index (for the first dimension) of a multidimensional Array
 of size "shape". Transform this into idx into the corresponding vector index
"""
function vectorizePattern(idx::Array{Int}, shape::Tuple)
  return [ floor(Int,(i-1)/size(idx,1))*shape[1]+idx[i] for i = 1:length(idx) ]
end

"""
    SamplingOp(pattern::Array{Int}, shape::Tuple)

buildsa `LinearOperator` which only returns the vector elements at positions
indicated by pattern.

# Arguents
* `pattern::Array{Int}` - indices to sample
* `shape::Tuple`        - size of the array to sample
"""
function SamplingOp(pattern::Array{Int}, shape::Tuple)
  ndims(pattern)>1 ?  idx = vectorizePattern(pattern, shape) : idx = pattern

  return opRestriction(idx, prod(shape))
end

function SamplingOp(pattern::Array{Bool})
  return LinearOperator(length(pattern), length(pattern), false, false
                  , x->vec(pattern).*x
                  , nothing
                  , y->vec(pattern).*y )
end

function zeropad(x::Array{T}, pattern::Matrix{Bool}) where T
  y = zeros(T,size(pattern))
  y[pattern] = x[pattern]
  return y
end


text:
# module Jive

onlyonce_evaluated = Dict{String,UInt}()
onlyonce_called = Dict{String,UInt}()

"""
    @onlyonce(block)

used to run the block only once.
"""
macro onlyonce(block)
    h = hash(block)
    if block.head === :block
        node = block.args[1]
        dir = pwd()
        linestr = string(relpath(String(node.file), dir), "#", node.line)
        haskey(onlyonce_evaluated, linestr) && h == onlyonce_evaluated[linestr] && return nothing
        onlyonce_evaluated[linestr] = h
    else
        linestr = string(block)
    end
    quot = quote
        if haskey(Jive.onlyonce_called, $linestr) && $h == Jive.onlyonce_called[$linestr]
            nothing
        else
            Jive.onlyonce_called[$linestr] = $h
            $block
        end
    end
    esc(quot)
end

# module Jive


text:
#####
##### Training.jl
##### Script to run a new training loop on old memory
##### without additional self play and with custom parameters
##### to monitor the performance of loss and related benchmarks
#####

using AlphaZero

include("games.jl")
GAME = get(ENV, "GAME", "double-dummy")
SelectedGame = GAME_MODULE[GAME]
using .SelectedGame: Game, Training

SESSION_DIR = "sessions/$GAME"

include("trials/params.jl")
for trial in trials
    AlphaZero.UserInterface.train_and_monitor(Game, SESSION_DIR, trial.label, 50, Training.benchmark, trial.params, trial.netparams)
end

max_stars_repo_name: btheodorou/Bridge-NN-Research

max_issues_repo_name: btheodorou/Bridge-NN-Research

max_forks_repo_name: btheodorou/Bridge-NN-Research


text:
module Unmarshal

# package code goes here

# Helper function
function prettyPrint(verboseLvl, str)
    tabs = ""
    for cntr=1:verboseLvl
        tabs = tabs * "\t"
    end
    println("$(tabs)$(str)")
end

export unmarshal # returns a reconstructed variable from a JSON parsed string

using JSON
import Missings: Missing, missing
import Nullables: Nullable

function unmarshal(DT :: Type, parsedJson :: String, verbose :: Bool = false, verboseLvl :: Int = 0)
    if (verbose)
        prettyPrint(verboseLvl, "$(DT) (String)")
        verboseLvl+=1
    end
    if DT <: Array
        [parsedJson]
    else
        DT(parsedJson)
    end
end

function unmarshal(::Type{Vector{E}}, parsedJson::Vector, verbose :: Bool = false, verboseLvl :: Int = 0) where E
    if (verbose)
        prettyPrint(verboseLvl, "Vector{$E}")
        verboseLvl+=1
    end

    [(unmarshal(E, field, verbose, verboseLvl) for field in parsedJson)...]
end

unmarshal(::Type{Array{E}}, xs::Vector, verbose :: Bool = false, verboseLvl :: Int = 0) where E = unmarshal(Vector{E}, xs, verbose, verboseLvl)

function unmarshal(::Type{Array{E, N}}, parsedJson::Vector, verbose :: Bool = false, verboseLvl :: Int = 0) where {E, N}
    if (verbose)
        prettyPrint(verboseLvl, "Array{$E, $N}")
        verboseLvl+=1
    end

    cat((unmarshal(Array{E,N-1}, x, verbose, verboseLvl) for x in parsedJson)..., dims=N)
end

"""
    unmarshal(T, dict[, verbose[, verboselvl]])

Reconstructs an object of Type T using the dictionary output of a `JSON.parse.`

Set verbose `true` to get debug information about how the data hierarchy is unmarshalled. This might be useful to track down parsing errors and/or mismatches between the JSON object and the Type definition.

# Example

```jldoctest
julia> using JSON

julia> var = randn(Float64, 5);  # Should work for most other variations of types you can think of

julia> unmarshal(typeof(var), JSON.parse(JSON.json(var)) ) == var
true
```

"""
function unmarshal(DT :: Type, parsedJson :: AbstractDict, verbose :: Bool = false, verboseLvl :: Int = 0)
    if (verbose)
            prettyPrint(verboseLvl, "$(DT) AbstractDict")
        verboseLvl+=1
    end

    if !isconcretetype(DT)
        throw(ArgumentError("Cannot unmarshal a non-leaf type $(DT) without a custom specialization"))
    end

    tup = ()
    for iter in fieldnames(DT)
        DTNext = fieldtype(DT,iter)
        if (verbose)
           prettyPrint(verboseLvl-1, "\\--> $(iter) <: $(DTNext) ")
        end

        if !haskey(parsedJson, string(iter))
            # check whether DTNext is compatible with any scheme for missing values
            val = if DTNext <: Nullable
                DTNext()
            elseif Missing <: DTNext
                missing
            elseif Nothing <: DTNext
                Nothing()
            else
                throw(ArgumentError("Key $(string(iter)) is missing from the structure $DT, and field is neither Nullable nor Missings nor Nothing-compatible"))
            end
        else
            val = unmarshal( DTNext, parsedJson[string(iter)], verbose, verboseLvl)
        end

        tup = (tup..., val)
    end

    DT(tup...)
end

function unmarshal(DT :: Type{T}, parsedJson :: Array{Any,N}, verbose :: Bool = false, verboseLvl :: Int = 0) where {T<:Tuple, N}
    if (verbose)
        prettyPrint(verboseLvl, "$(T) $(N) Dimensions, length $(length(parsedJson))")
        verboseLvl += 1
    end

    ((unmarshal(fieldtype(T,1), field, verbose, verboseLvl) for field in parsedJson)...,)
end

function unmarshal(DT :: Type{T}, parsedJson :: AbstractDict, verbose :: Bool = false, verboseLvl :: Int = 0) where T <: Dict
    if (verbose)
        prettyPrint(verboseLvl, "$(DT) Dict ")
        verboseLvl += 1
    end
    val = DT()
    for iter in keys(parsedJson)
        if (verbose)
           prettyPrint(verboseLvl-1, "\\--> $(iter) ")
        end
        tmp = unmarshal(valtype(DT), parsedJson[iter], verbose, verboseLvl)
        if keytype(DT) <: AbstractString
            val[iter] = tmp
        else
            try
                val[unmarshal(keytype(DT),JSON.parse(iter),verbose, verboseLvl)] = tmp # Use JSON.parse and Unmarshal to cast from type of iter to ketype(DT)
            catch ex
                val[keytype(DT)(iter)] = tmp  # Try direct casting, which will hopefully generate a readable enough exception error if it fails
            end
        end
    end
    val
end

unmarshal(::Type{T}, x::Number, verbose :: Bool = false, verboseLvl :: Int = 0) where T<:Number = T(x)
unmarshal(::Type{Nullable{T}}, x, verbose :: Bool = false, verboseLvl :: Int = 0) where T = Nullable(unmarshal(T, x))
unmarshal(::Type{Nullable{T}}, x::Nothing, verbose :: Bool = false, verboseLvl :: Int = 0) where T = Nullable{T}()
unmarshal(::Type{Union{T,Missing}}, x, verbose :: Bool = false, verboseLvl :: Int = 0) where T = unmarshal(T, x, verbose, verboseLvl)
unmarshal(::Type{Union{T,Nothing}}, x, verbose :: Bool = false, verboseLvl :: Int = 0) where T = unmarshal(T, x, verbose, verboseLvl)

unmarshal(T::Type, x, verbose :: Bool = false, verboseLvl :: Int = 0) =
    throw(ArgumentError("no unmarshal function defined to convert $(typeof(x)) to $(T); consider providing a specialization"))

#The following functions have been moved here to the end, since they appear to be unused by codecov. Consider deleting them unless actual use cases can be identified
function unmarshal(::Type{Vector{E}}, parsedJson::Number, verbose :: Bool = false, verboseLvl :: Int = 0) where E<:Number
    if (verbose)
        prettyPrint(verboseLvl, "Vector{$E}")
        verboseLvl+=1
    end

    [(unmarshal(E, field, verbose, verboseLvl) for field in parsedJson)...]
end

unmarshal(::Type{Array{E}}, xs::Number, verbose :: Bool = false, verboseLvl :: Int = 0) where E<:Number = unmarshal(Vector{E}, xs, verbose, verboseLvl)

function unmarshal(::Type{Array{E, N}}, parsedJson::Number, verbose :: Bool = false, verboseLvl :: Int = 0) where {E<:Number, N}
    if (verbose)
        prettyPrint(verboseLvl, "Array{$E, $N}")
        verboseLvl+=1
    end

    cat(E(parsedJson), dims=N)
end

end # module


text:
#
# This file is part of module
#   ModiaMedia (ModiaMedia/src/ModiaMedia.jl)
#

### Data structures that are additionally available for media <: PartialMixtureMedium -----------------------

saturationPressure(m::AbstractMedium, TSat) = undefinedFunction("saturation pressure", m)


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule xrootdgo_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("xrootdgo")
JLLWrappers.@generate_main_file("xrootdgo", UUID("9d84c17e-11f2-50ef-8cc9-e9701362097f"))
end  # module xrootdgo_jll


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

struct StarsApi <: SwaggerApi
    client::Swagger.Client
end

function _swaggerinternal_starsAdd(_api::StarsApi; channel=nothing, file=nothing, file_comment=nothing, timestamp=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", Dict{String, Any}, "/stars.add", ["slackAuth"])
    Swagger.set_param(_ctx.form, "channel", channel)  # type String
    Swagger.set_param(_ctx.form, "file", file)  # type String
    Swagger.set_param(_ctx.form, "file_comment", file_comment)  # type String
    Swagger.set_param(_ctx.form, "timestamp", timestamp)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded", "application/json"] : [_mediaType])
    return _ctx
end

"""

Adds a star to an item.
Param: channel::String
Param: file::String
Param: file_comment::String
Param: timestamp::String
Return: Dict{String, Any}
"""
function starsAdd(_api::StarsApi; channel=nothing, file=nothing, file_comment=nothing, timestamp=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_starsAdd(_api; channel=channel, file=file, file_comment=file_comment, timestamp=timestamp, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function starsAdd(_api::StarsApi, response_stream::Channel; channel=nothing, file=nothing, file_comment=nothing, timestamp=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_starsAdd(_api; channel=channel, file=file, file_comment=file_comment, timestamp=timestamp, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_starsList(_api::StarsApi; count=nothing, page=nothing, cursor=nothing, limit=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", Dict{String, Any}, "/stars.list", ["slackAuth"])
    Swagger.set_param(_ctx.query, "count", count)  # type String
    Swagger.set_param(_ctx.query, "page", page)  # type String
    Swagger.set_param(_ctx.query, "cursor", cursor)  # type String
    Swagger.set_param(_ctx.query, "limit", limit)  # type Int32
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded"] : [_mediaType])
    return _ctx
end

"""

Lists stars for a user.
Param: count::String
Param: page::String
Param: cursor::String
Param: limit::Int32
Return: Dict{String, Any}
"""
function starsList(_api::StarsApi; count=nothing, page=nothing, cursor=nothing, limit=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_starsList(_api; count=count, page=page, cursor=cursor, limit=limit, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function starsList(_api::StarsApi, response_stream::Channel; count=nothing, page=nothing, cursor=nothing, limit=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_starsList(_api; count=count, page=page, cursor=cursor, limit=limit, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_starsRemove(_api::StarsApi; channel=nothing, file=nothing, file_comment=nothing, timestamp=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", Dict{String, Any}, "/stars.remove", ["slackAuth"])
    Swagger.set_param(_ctx.form, "channel", channel)  # type String
    Swagger.set_param(_ctx.form, "file", file)  # type String
    Swagger.set_param(_ctx.form, "file_comment", file_comment)  # type String
    Swagger.set_param(_ctx.form, "timestamp", timestamp)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded", "application/json"] : [_mediaType])
    return _ctx
end

"""

Removes a star from an item.
Param: channel::String
Param: file::String
Param: file_comment::String
Param: timestamp::String
Return: Dict{String, Any}
"""
function starsRemove(_api::StarsApi; channel=nothing, file=nothing, file_comment=nothing, timestamp=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_starsRemove(_api; channel=channel, file=file, file_comment=file_comment, timestamp=timestamp, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function starsRemove(_api::StarsApi, response_stream::Channel; channel=nothing, file=nothing, file_comment=nothing, timestamp=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_starsRemove(_api; channel=channel, file=file, file_comment=file_comment, timestamp=timestamp, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export starsAdd, starsList, starsRemove


text:
fatalerrors = false
anyerrors = false

my_tests = ["dict.jl",
            "dataframes.jl",
            "grouping.jl",
            "chaining.jl",
            "linqmacro.jl"]

println("Running tests:")

for my_test in my_tests
    try
        include(my_test)
        println("\t\033[1m\033[32mPASSED\033[0m: $(my_test)")
    catch
        global anyerrors = true
        println("\t\033[1m\033[31mFAILED\033[0m: $(my_test)")
        if fatalerrors
            rethrow()
        end
    end
end

if anyerrors
    throw("Tests failed")
end


text:
struct Coordinate{T,F<:Function} <: AbstractArray{T,1}
    u::Vector{T} #Coordinates
    R::F         #Transformation (x,y,z,...) = R(u₁,u₂,u₃,...)
end

# Array Interface for Coordinates
parent(A::T) where {T<:Coordinate} = A.u
size(A::T) where {T<:Coordinate} = size(A.u)
axes(A::T) where {T<:Coordinate} = axes(A.u)
parenttype(::Type{Coordinate{T,F}}) where {T,F} = Vector{T}
IndexStyle(::Type{T}) where {T<:Coordinate} = IndexStyle(parenttype(T))

@propagate_inbounds getindex(A::T, i::Int) where {T<:Coordinate} = A.u[i]
@propagate_inbounds setindex!(A::Coordinate{T}, v::T, i::Int) where {T} = A.u[i] = v

similar(A::Coordinate{T,F}) where {T,F} = Coordinate(Array{T}(undef, size(A)), A.R)
similar(A::Coordinate, ::Type{T}) where {T} = Coordinate(Array{T}(undef, size(A)), A.R)
function similar(A::Coordinate, ::Type{T}, dims::Tuple{Vararg{Int64,N}}) where {T,N}
    return Coordinate(Array{T}(undef,dims), A.R)
end

# Standard Transformations
function Cartesian(u)
    u
end

function CartesianCoordinate(u)
    Coordinate(u,Cartesian)
end

function Cylindrical(u)
    # u -> (r,θ,z) r ∈ [0,∞), θ ∈ [0,2π), z ∈ (-∞,∞)
    r = u[1]
    θ = u[2]
    z = u[3]
    [r*cos(θ), r*sin(θ), z]
end

function CylindricalCoordinate(u)
    Coordinate(u,Cylindrical)
end

function Spherical(u)
    # u -> (r,θ,ϕ) r ∈ [0,∞), θ ∈ [0,2π), ϕ ∈ [0,π]
    r = u[1]
    θ = u[2]
    ϕ = u[3]
    [r*cos(θ)*sin(ϕ), r*sin(θ)*sin(ϕ), r*cos(ϕ)]
end

function SphericalCoordinate(u)
    Coordinate(u,Spherical)
end


text:
##########################################
#         Importing
##########################################

include("FieldArithmetic.jl")

using .FieldArithmetic
import Base: +,-,*,/,^,==,inv,show,zeros,gcd,gcdx,iszero,isone,zero,one,identity
import .FieldArithmetic: order

abstract type EllipticCurve end
abstract type EllipticCurveElem end

##########################################
#
#Elliptic Curve in twisted Hessian form
#
##########################################

struct Affine_EC_TwiHes{T <: CRing, S <: CRingElem} <: EllipticCurve
    parent::T
    a::S
    d::S
end

struct Affine_EC_TwiHesElem{T <: CRingElem, S <: CRing} <: EllipticCurveElem
    x::T
    y::T
    curve::Affine_EC_TwiHes{S, T}
end

# Curve-Constructor
function Affine_EC_TwiHes(R::T, a::S, d::S) where {T <: CRing, S <: CRingElem}
    if isSingular(R,a,d)
        error("Curve is singular!")
    end
    return Affine_EC_TwiHes{T,S}(R, a, d)
end

# Element-Constructor
function (E::Affine_EC_TwiHes{T,S})(x::S, y::S) where {T <: CRing, S <: CRingElem}
    if !validPoint(x, y, E)
        error("Point ($x, $y) is not valid!!!")
    end
    return Affine_EC_TwiHesElem{S, T}(x,y,E)
end

###### Printing ######
function show(io::IO, R::Affine_EC_TwiHes{T,S}) where {T <: CRing, S <: CRingElem}
    print(io, "Elliptic curve [$(R.a)]*x^3 + y^3 + 1 = [$(R.d)]*x*y defined over $(R.parent)")
end

function show(io::IO, P::Affine_EC_TwiHesElem{T,S}) where {T <: CRingElem, S <: CRing}
    print(io, "([$(P.x)], [$(P.y)])")
end

###### Helper-functions ######
function isSingular(R::T, a::S, d::S) where {T <: CRing, S <: CRingElem}
    return (iszero(a*(27*a - d^3)))
end

function validPoint(x::S, y::S, E::Affine_EC_TwiHes{T,S}) where {T <: CRing, S <: CRingElem}
    return E.a*x^3 + y^3 + one(E.parent) == E.d*x*y
end

function identity(E::Affine_EC_TwiHes{T,S}) where {T <: CRing, S <: CRingElem}
    F = E.parent
    return E(zero(F),-one(F))
end

function isIdentity(P::Affine_EC_TwiHesElem{T,S}) where {T <: CRingElem, S <: CRing}
    return P == identity(P.curve)
end

function sameCurve(P::Affine_EC_TwiHesElem{T,S}, Q::Affine_EC_TwiHesElem{T,S}) where {T <: CRingElem, S <: CRing}
    return (P.curve == Q.curve)
end

#Snail speed, but oh well
function order(P::Affine_EC_TwiHesElem{T,S}) where {T <: CRingElem, S <: CRing}
    i = 1
    Pi = P
    while !isIdentity(Pi)
        Pi = Pi + P
        i += 1
    end
    return i
end

###### Basic operators #######

#How to invert (x,0)??
#Answer: A point (x,0) implies that a*x^3 = -1, or x = a^(-1/3), which is only possible if a is a cube.
function -(P::Affine_EC_TwiHesElem{T,S}) where {T <: CRingElem, S <: CRing}
    E = P.curve
    return E(P.x*inv(P.y),inv(P.y))
end

function +(P::Affine_EC_TwiHesElem{T,S}, Q::Affine_EC_TwiHesElem{T,S}) where {T <: CRingElem, S <: CRing}
    if !(sameCurve(P, Q))
        error("Cannot add points from different curves!")
    end
    if P == Q
        return double(P)
    end
    E = P.curve
    a = E.a
    x1, y1, x2, y2 = P.x, P.y, Q.x, Q.y
    z1 = x1 - y1^2 * x2 * y2
    z2 = y1 * y2^2 - a * x1^2 * x2
    z3 = inv(a * x1 * y1 * x2^2 - y2)
    return E(z1*z3, z2*z3)
end

function -(P::Affine_EC_TwiHesElem{T,S}, Q::Affine_EC_TwiHesElem{T,S}) where {T <: CRingElem, S <: CRing}
    return P+(-Q)
end

function double(P::Affine_EC_TwiHesElem{T,S}) where {T <: CRingElem, S <: CRing}
    E = P.curve
    a = E.a
    x, y = P.x, P.y
    z1 = x - y^3 * x
    z2 = y^3 - a * x^3
    z3 = inv(a * y * x^3 - y)
    return E(z1*z3, z2*z3)
end

#Regard E as a Z-module... This is again square and multiply (or strictly speaking double and add)
function *(n::U, P::Affine_EC_TwiHesElem{T,S}) where {T <: CRingElem, S <: CRing, U <: Integer}
    if isIdentity(P) || n == 1
        return P
    end
    E = P.curve
    if n == 0
        return identity(E)
    end
    if n < 0
        P = -P
        n = -n
    end
    bit = U(1) << (ndigits(n, base = 2) - 1)
    Q = P
    bit >>= 1
    while bit != 0
        Q = double(Q)
        if (bit & n) != 0
            Q += P
        end
        bit >>=1
    end
    return Q
end

function *(P::Affine_EC_TwiHesElem{T,S}, n::U) where {T <: CRingElem, S <: CRing, U <: Integer}
    return n*P
end

function ==(P::Affine_EC_TwiHesElem{T,S}, Q::Affine_EC_TwiHesElem{T,S}) where {T <: CRingElem, S <: CRing}
    return sameCurve(P,Q) && (P.x == Q.x) && (P.y == Q.y)
end

############################
#
#   Elliptic Curve Stuff
#
############################

#See notebook for calculation. It is calculated as the j-invariant of the isomorphic Weierstrass-curve
function jInvariant(E::Affine_EC_TwiHes{T,S}) where {T <: CRing, S <: CRingElem}
    a = E.a
    d = E.d
    over = -(216*a*d+d^4)^3
    under = a*(27*a-d^3)^3
    return over*inv(under)
end

function toWeierstrass(E::Affine_EC_TwiHes{T,S}) where {T <: CRing, S <: CRingElem}
    a = E.a
    d = E.d
    w1 = -8*a*d - d^4*inv(27*one(E.parent))
    w2 = 16*a^2 + 40*a*d^3*inv(27*one(E.parent)) - 2*d^6*inv(729*one(E.parent))
    println(E)
    println("is isomorphic (over the algebraic clousure) to the curve")
    println("y^2 = x^3 + [$w1]*x + [$w2]")
end

#FromWeierstrass: Find point of order 3, then transform to triangular form and finally to twisted Hessian

############################
#
#           Other
#
############################
function allPointsNaive(E::Affine_EC_TwiHes{T,S}) where {T <: CRing, S <: CRingElem}
    points = Affine_EC_TwiHesElem{S,T}[]
    for x in allElements(E.parent)
        for y in allElements(E.parent)
            if validPoint(x, y, E)
                push!(points,E(x,y))
            end
        end
    end
    return points
end

#Disregard
function WeierjInvariant(a, b)
    delta = -16*(4*a^3+27*b^2)
    j = -1728*(4*a)^3*inv(delta)
    println(j)
end

############################
#
#        Isogenies
#
############################

struct Affine_isogeny
    domain::EllipticCurve
    codomain::EllipticCurve
    carryP::Function
end

# Constructor
function Affine_isogeny(E::EllipticCurve, P::EllipticCurveElem)
    if isIdentity(2*P)
        return deg2_isogeny(E,P)
    end
    if isIdentity(3*P)
        return deg3_isogeny(E,P)
    end
    error("$P is not of order 2 or 3")
end

function deg3_isogeny(E::EllipticCurve, P::EllipticCurveElem)

end

function deg2_isogeny(E::Affine_EC_TwiHes, P::Affine_EC_TwiHesElem)
    a = E.a
    d = E.d
    F = E.parent
    A = a^2
    s = P.x
    t = P.y
    D = (-d+6*inv(s*t))*inv(s)
    E1 = Affine_EC_TwiHes(F, A, D)
    phi(x,y) = (x*((x-s*t*y^2)*inv(a*s^2*x*y-t)),y*((y*t^2-a*s*x^2)*inv(a*s^2*x*y-t)))
    return Affine_isogeny(E, E1, phi)
end

# Evaluate a point
function (phi::Affine_isogeny)(P::EllipticCurveElem)
    if P.curve != phi.domain
        error("Point $P is not from the right curve")
    end
    (x,y) = phi.carryP(P.x,P.y)
    return phi.codomain(x,y)
end

###### Printing ######
function show(io::IO, phi::Affine_isogeny)
    print(io, "Isogeny φ : E_1 -> E_2 where \nE_1 = $(phi.domain)\nE_2 = $(phi.codomain)")
end


text:
using LibSerialPort

mcu = open( "/dev/ttyUSB0", 9600 )
flush(mcu)

println("Listening for UI to stop VIM...")
function run_listener()
    mcu_message = ""
    while true
       @async mcu_message *= read(mcu, String)
       if occursin("KillVIM!", mcu_message)
           #runs the following command to kill vi!
           #ps aux | grep -ie " vi" | awk '{print $2}' | xargs kill -9
           try
               run(pipeline(`ps aux`, `grep -ie " vi"`, `awk '{print $2}'`, `xargs kill -9`))
           catch
               println("Error thrown. May not be fatal")
           end
           mcu_message = ""
        end
       sleep(0.01)
    end
end

run_listener()

max_stars_repo_path: Julia Code/arduino.jl

max_stars_repo_name: jonathanwcrane/VIMKiller

max_issues_repo_path: Julia Code/arduino.jl

max_issues_repo_name: jonathanwcrane/VIMKiller

max_forks_repo_path: Julia Code/arduino.jl

max_forks_repo_name: jonathanwcrane/VIMKiller


text:
using RecipesBase

function poly_interval(p::AbstractPolynomial)

    # use  restricted domain, if finite
    A,B =  first(domain(p)), last(domain(p))
    if !isinf(A) && !isinf(B)
        if isopen(domain(p))
            Delta = (B-A)/100
            A += Delta
            B -= Delta
        end
        return A:(B-A)/100:B
    end

    # Find points of interest
    zero_pts = roots(p)
    crit_pts = roots(derivative(p, 1))
    infl_pts = roots(derivative(p, 2))
    pts = sort([ real(pt) for pt in [zero_pts; crit_pts; infl_pts] if isreal(pt) ])
    # Choose a range that shows all interesting points with some margin
    min_x, max_x = length(pts) > 0 ? (pts[1], pts[end]) : (-1, 1)
    d = max(max_x - min_x, 1)
    a = min_x - d / 5
    b = max_x + d / 5

    Delta  = b -  a

    return a:Delta/50:b
end

poly_label(p::AbstractPolynomial) = sprint(printpoly, p)

#@recipe function f(p::AbstractPolynomial, x = poly_interval(p))
#    label --> poly_label(p)
#    x, p.(x)
#end

@recipe function f(p::AbstractPolynomial)
    label --> poly_label(p)
    xlims --> extrema(poly_interval(p))
    x -> p(x)
end

@recipe function f(p::AbstractPolynomial, a, b)
    label --> poly_label(p)
    step = (b - a) / 100
    xs = a:step:b
    ys = p.(xs)
    xs, ys
end


text:
"""
The way the Louvain code works, you start with n clusters (singletons), and as nodes
move, many of these become empty. At certain points in the algorithm, it is
helpful to renumber the cluster IDs, so that they go from 1 to M where
M = number of non-empty clusters.

e.g. c = [2 3 5 3 9] ---> c = [1 2 3 2 4]

This version does this with both the Clusters array of arrays AND
the cluster indicator vector c:

c[i] = (integer) cluster ID that node i belongs to
Clusters[j] = (integer array) node IDs for nodes in cluster j
"""
function renumber(c::Vector{Int64},Clusters::Vector{Vector{Int}})

    n = length(c)
    map = sort(unique(c))     # map from new cluster ID to old cluster ID
    old2new = Dict()    # map from old cluster ID to new cluster ID
    for i = 1:length(map)
        old2new[map[i]] = i
    end
    cnew = zeros(Int64,n)

    Clusters = Clusters[map]

    # Rename the clusters
    for i = 1:n
        # newClus = findfirst(x->x == c[i],map)
        newClus = old2new[c[i]]
        cnew[i] = newClus
        push!(Clusters[newClus],i)
    end

    return cnew, Clusters

end

"""
See above function, this function doesn't care about the Clusters array
"""
function renumber(c::Vector{Int64})

    n = length(c)
    map = sort(unique(c))     # map from new cluster ID to old cluster ID
    old2new = Dict()          # map from old cluster ID to new cluster ID
    for i = 1:length(map)
        old2new[map[i]] = i
    end
    cnew = zeros(Int64,n)

    # Rename the clusters
    for i = 1:n
        cnew[i] = old2new[c[i]]
    end

    return cnew

end

# These are some useful functions for converting back and forth between
# different ways to store hypergraphs.
function NeighborList(He2n::SparseArrays.SparseMatrixCSC{Float64,Int64},Hn2e::SparseArrays.SparseMatrixCSC{Float64,Int64})
    """
    NeighborList: given hypergraph (edge x node) indicence and (node x edge)
    incidence, return a list of neighbors for each node.
    """
    m,n = size(He2n)
    rp = He2n.rowval
    ci = He2n.colptr
    Neighbs = Vector{Vector{Int64}}()
    for i = 1:n
        # Get list of edges that node i is in
        Edges = rp[ci[i]:ci[i+1]-1]

        # build a set of neighboring nodes
        NeighbSet = Vector{Int64}()
        for j = 1:length(Edges)
            ej = Edges[j]               # For each node ej adjacent to I...
            nodes = getnodes(Hn2e,ej)   # ...get nodes in edge ej
            append!(NeighbSet,nodes)
        end
        push!(Neighbs, sort(setdiff(unique(NeighbSet),i)))
    end

    return Neighbs
end

function NeighborList(node2edge::Vector{Vector{Int64}},edge2node::Vector{Vector{Int64}})
    """
    NeighborList: given edge to node list and node to edge list, compute
    neighbors
    """
    n = length(node2edge)
    m = length(edge2node)
    Neighbs = Vector{Vector{Int64}}()
    for i = 1:n
        # Get list of edges that node i is in
        Edges = node2edge[i]

        # build a set of neighboring nodes
        NeighbSet = Vector{Int64}()
        for j = 1:length(Edges)
            ej = Edges[j]               # For each node ej adjacent to I...
            nodes = edge2node[ej]       # ...get nodes in edge ej
            append!(NeighbSet,nodes)
        end
        push!(Neighbs, sort(setdiff(unique(NeighbSet),i)))
    end

    return Neighbs
end

function getedges(He2n::SparseArrays.SparseMatrixCSC{Float64,Int64},I::Int64)
    """
    Given the edge-by-node incidence matrix He2n, return the set of hyperedges
    that node I is in.
    """
    first = He2n.colptr[I]
    last = He2n.colptr[I+1]-1
    edges = He2n.rowval[first:last]

    return edges
end

function getnodes(Hn2e::SparseArrays.SparseMatrixCSC{Float64,Int64},J::Int64)
    """
    Given the node-by-edge incidence matrix Hn2e, return the set of nodes
    contained in hyperedge J
    """
    first = Hn2e.colptr[J]
    last = Hn2e.colptr[J+1]-1
    nodes = Hn2e.rowval[first:last]
    mult = Hn2e.nzval[first:last]
    edge = Vector{Int64}()
    # need to adjust for multiplicities
    for t = 1:length(nodes)
        node = nodes[t]
        for k = 1:mult[t]
            push!(edge,node)
        end
    end

    return edge
end

function incidence2elist(H::SparseArrays.SparseMatrixCSC{Float64,Int64},nodelist::Bool=false)
    """
    This converts a hypergraph in incidence matrix form to hyperedge list form.
    Incidence matrix form:  H[e,u] = 1  iff node u is in hyperedge e
    Edgelist form: Hyperedges[j] = array of nodes in hyperedge j

    nodelist == true: means you actually want to get a map from node IDs to hyperedge ids
    """
    if ~nodelist
        # unless you want the node2edge map, transpose first
        H = SparseArrays.sparse(H')
    end
    rp = H.rowval
    ci = H.colptr
    nz = H.nzval
    Hyperedges = Vector{Vector{Int64}}()
    n,m = size(H)

    for i = 1:m
        startedge = ci[i]
        endedge = ci[i+1]-1
        nodes = rp[startedge:endedge]
        mult = nz[startedge:endedge]
        edge = Vector{Int64}()
        # need to adjust for multiplicities
        for t = 1:length(nodes)
            node = nodes[t]
            for k = 1:mult[t]
                push!(edge,node)
            end
        end

        push!(Hyperedges,edge)
    end
    return Hyperedges
end

function elist2incidence(Hyperedges::Vector{Vector{Int64}}, N::Int64)
    """
    Converts a hyperedge list into a binary incidence matrix for the hypergraph.

    This is the exact inverse of incidence2elist
    """
    U = Vector{Int64}()
    E = Vector{Int64}()
    M = length(Hyperedges)
    for enum = 1:length(Hyperedges)
        e = Hyperedges[enum]
        for node in e
            push!(U,node)
            push!(E,enum)
        end
    end

    H = SparseArrays.sparse(E,U,ones(length(U)),M,N)
    return H
end

function hyperedge_formatting(H::hypergraph)
    """
    Convert a hypergraph from dictionary-based format (type 'hypergraph')
    to Hyperedge list with weights vector
    """
    Hyperedges = Vector{Vector{Int64}}()
    weights = Vector{Float64}()
    for key in keys(H.E)
        ed = collect(keys(H.E[key]))
        wt = collect(values(H.E[key]))
        append!(weights,wt)
        append!(Hyperedges,ed)
    end
    return Hyperedges, weights
end

function hypergraph2incidence(H::hypergraph)
    """
    Convert dictionary hypergraph format to binary edge-by-node incidence matrix
    """
    Hyperedges = Vector{Vector{Int64}}()
    weights = Vector{Float64}()
    for key in keys(H.E)
        ed = collect(keys(H.E[key]))
        wt = collect(values(H.E[key]))
        append!(weights,wt)
        append!(Hyperedges,ed)
    end
    N = length(H.D)
    He2n = elist2incidence(Hyperedges,N)
    return He2n, weights
end

function EdgeMap(H::hypergraph)
    """
    Given a hypergraph, extract a map from node indices to hyperedge indices
    """
    n = length(H.D)
    Hyp, weights = hyperedge_formatting(H)
    Hin = elist2incidence(Hyp,n)

    return incidence2elist(Hin,true)
end

max_stars_repo_name: nveldt/TestPkg

max_issues_repo_name: nveldt/TestPkg

max_forks_repo_name: nveldt/TestPkg


text:
module ImplicitArrays

include("compat.jl")

include("block.jl")
export BlockMatrix

include("interaction.jl")
export ConstInteractionFunction, GenericInteractionFunction, InteractionFunction, InteractionMatrix

include("fixedvalue.jl")
export FixedValueArray

include("rowprojection.jl")
export RowProjectionMatrix, RowProjectionVector

end # module


text:
using BinaryBuilder

name = "github_release"
version = v"0.10.0"

# Collection of sources required to build github_release
sources = [
    ArchiveSource("https://github.com/github-release/github-release/archive/refs/tags/v$(version).zip",
                  "369fc1ccdf9e0250b1974cd1089c71f994e12f29d9524205a16d1a2c29364396"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/github-release*
mkdir -p ${bindir}
go mod init github-release
go mod tidy
go mod vendor
go build -o ${bindir}
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    ExecutableProduct("github-release", :github_release),
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; compilers=[:c, :go])

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil