text:
using Test
using Metrics
using util

function createRecommendationList(correct_pos; range=100)
    extra_id = range * 2
    dup_id = extra_id + 2
    query = Report(extra_id + 1, dup_id, 1000, 10, [[0]])

    recommendation_list= Vector{Tuple{Report, Float64}}(undef, range)
    recommendation_list[1] =  (Report(1, 1, 1000, 10, [[0]]), 100.0)
    recommendation_list[2] =  (Report(extra_id+3, 1, 1000, 10, [[0]]), 100.0)

    recommendation_list[3] =  (Report(2, 2, 1000, 10, [[0]]), 100.0)
    recommendation_list[4] =  (Report(extra_id+4, 2, 1000, 10, [[0]]), 100.0)
    recommendation_list[5] =  (Report(extra_id+5, 2, 1000, 10, [[0]]), 100.0)

    recommendation_list[6] =  (Report(3, 3, 1000, 10, [[0]]), 100.0)

    recommendation_list[7] =  (Report(4, 1, 1000, 10, [[0]]), 100.0)
    recommendation_list[8] =  (Report(extra_id+6, 1, 1000, 10, [[0]]), 100.0)
    recommendation_list[9] =  (Report(extra_id+7, 1, 1000, 10, [[0]]), 100.0)

    for i in 9:range
        recommendation_list[i] = (Report(i, i, 1000, 10, [[0]]), 100.0)
    end

    if correct_pos > 0
        recommendation_list[correct_pos] = (Report(extra_id+8, dup_id, 1000, 10, [[0]]), 100.0)
    end

    return (query, recommendation_list)
end

map_rr = createMAP_RecallRate(;group_by_master=false)

correct_pos_lis = [1,1,1,2,2,2,3,3,3,4,4,4,5,5,6,7,19,20,20,21,50,60,44,43,80,99, -1]

for correct_pos in correct_pos_lis
    # println("Correct: $(correct_pos)")

    query1, recommendation_list = createRecommendationList(correct_pos)
    # println("\tdup_id: $(recommendation_list[correct_pos][1].dup_id)")
    # println("\tdup_id: $([recommendation_list[i][1].dup_id for i in 1:length(recommendation_list)])")
    update_metric!(map_rr, query1, recommendation_list)
end

result = compute(map_rr)

rr_values = Vector{Float64}(undef, 20)

for (k,v) in result["rr"]
    rr_values[k] = v
end

@test rr_values == [3.0,6.0,9.0,12.0,14.0,15.0,16.0,16.0,16.0,16.0,16.0,16.0,16.0,16.0,16.0,16.0,16.0,16.0,17.0, 19.0]./length(correct_pos_lis)
@test result["map"] == 0.2690750073903208

map_rr = createMAP_RecallRate(;group_by_master=true)

for correct_pos in correct_pos_lis
    query, recommendation_list = createRecommendationList(correct_pos)
    update_metric!(map_rr, query, recommendation_list)
end
result = compute(map_rr)

rr_values = Vector{Float64}(undef, 20)

for (k,v) in result["rr"]
    rr_values[k] = v
end

@test rr_values == [3.0,9.0,15.0,16.0,16.0,16.0,16.0,16.0,16.0,16.0,16.0,16.0,16.0,17.0,19.0,20.0,20.0,20.0,20.0, 20.0]./length(correct_pos_lis)
@test result["map"] == 0.3197627552402343

max_stars_repo_name: irving-muller/FaST

max_issues_repo_name: irving-muller/FaST

max_forks_repo_name: irving-muller/FaST


text:
#=
Utilities for combining and manipulating TS objects using their indexes
=#

import Base: hcat, vcat, merge

function overlaps(x::AbstractArray, y::AbstractArray, n::Int=1)
    if n == 1
        xx = falses(x)
        @inbounds for i = 1:size(x,1), j = 1:size(y,1)
            if x[i] == y[j]
                xx[i] = true
            end
        end
        return xx
    elseif n == 2
        yy = falses(y)
        @inbounds for i = 1:size(x,1), j = 1:size(y,1)
            if x[i] == y[j]
                yy[i] = true
            end
        end
        return yy
    else
        error("Argument `n` must be either 1 (x) or 2 (y).")
    end
end

function partner(x::TS, y::TS)
    yy = !overlaps(x.index, y.index) .* NaN
    yy[!isnan(yy),:] = y.values
    return ts(yy, x.index, y.fields)
end

function ojoin(x::TS, y::TS)
    idx = union(x.index, y.index)
    xna = setdiff(idx, x.index)
    yna = setdiff(idx, y.index)
    xi = sortperm([x.index; xna])
    yi = sortperm([y.index; yna])
    xvals = [x.values; fill(NaN, (length(xna), size(x,2)))][xi,:]
    yvals = [y.values; fill(NaN, (length(yna), size(y,2)))][yi,:]
    return ts([xvals yvals], idx, [x.fields; y.fields])
end

function ijoin(x::TS, y::TS)
    idx = intersect(x.index, y.index)
    return ts([x[idx].values y[idx].values], idx, [x.fields; y.fields])
end

ljoin(x::TS, y::TS) = ts([x.values partner(x,y).values], x.index, [x.fields; y.fields])
rjoin(x::TS, y::TS) = ts([partner(y,x).values y.values], y.index, [x.fields; y.fields])
hcat(x::TS, y::TS) = ojoin(x, y)
hcat(x::TS) = x
function hcat(series::TS...)
    out = series[1]
    @inbounds for j = 2:length(series)
        out = [out series[j]]
    end
    return out
end
function vcat(x::TS, y::TS)
    @assert size(x,2) == size(y,2) "Dimension mismatch: Number of columns must be equal."
    return TS([x.values;y.values], [x.index;y.index], x.fields)
end
function vcat(series::TS...)
    out = series[1]
    @inbounds for j = 2:length(series)
        out = vcat(out, series[j])
    end
    return out
end

function merge(x::TS, y::TS; join::Char='o')
    @assert join == 'o' || join == 'i' || join == 'l' || join == 'r' "`join` must be 'o', 'i', 'l', or 'r'."
    if join == 'o'
        return ojoin(x, y)
    elseif join == 'i'
        return ijoin(x, y)
    elseif join == 'l'
        return ljoin(x, y)
    elseif join == 'r'
        return rjoin(x, y)
    end
end

#===============================================================================
                COMBINING/MERGING WITH OTHER TYPES
===============================================================================#
hcat(x::TS, y::AbstractArray) = ojoin(x, ts(y, x.index))
hcat(y::AbstractArray, x::TS) = ojoin(ts(y, x.index), x)
hcat(x::TS, y::Number) = ojoin(x, ts(fill(y,size(x,1)), x.index))
hcat(y::Number, x::TS) = ojoin(ts(fill(y,size(x,1)), x.index), x)
function hcat{V}(series::TS, arrs::AbstractArray{V}...)
    n = size(series,1)
    k = length(arrs)
    out = zeros(V, (n,k))
    @inbounds for j = 1:k
        out[:,j] = arrs[j]
    end
    return [series out]
end
function hcat{V<:Number}(series::TS, nums::V...)
    n = size(series,1)
    k = length(nums)
    out = zeros(V, (n,k))
    @inbounds for j = 1:k
        out[:,j] = fill(n, nums[j])
    end
    return [series out]
end

ext: jl


text:
#states constructors
struct State{T<:AbstractVector}
    st::T
    ope::Op
    State(st,ope) = log(2,length(st)) != dim(ope) ? error("lenght of vector does not match dimension") : new{typeof(st)}(st/sqrt(st'*st),ope)
end

st(s::State) = s.st

ope(s::State) = s.ope

dim(s::State) = dim(s.ope)

basis(s::State) = basis(s.ope)

typ(s::State) = eltype(s.st)

function rhosp(sta::State, prec=15)
    n = dim(ope(sta))
    rhospv = spzeros(typ(sta),n,n)
    estate = st(sta)
    o = ope(sta)
    for i in 1:n
        for j in 1:n
            rhospv[i,j] = round(estate'*ada(o, i, j)*estate, digits = prec)
        end
    end
    return rhospv
end

function kqsp(sta::State, prec=15)
    n = dim(ope(sta))
    k = spzeros(typ(sta),n,n)
    estate = st(sta)
    o = ope(sta)
    for i in 1:n
        for j in 1:n
            k[i,j] = round(estate'*aa(o, j, i)*estate, digits = prec)
        end
    end
    return k
end

rhoqsp(s::State, prec=15) = [rhosp(s, prec) kqsp(s, prec); kqsp(s, prec)' I-rhosp(s, prec)']

function non_zero(c, prec=15)
    a = similar(c, Int)
    cof = similar(c, Float64)
    count = 1
    @inbounds for i in eachindex(c)
        a[count] = i
        cof[count] = round(c[i], digits = prec)
        count += (round(c[i], digits = prec) != zero(eltype(c)))
    end
    return resize!(a, count-1), resize!(cof, count-1)
end

#states constructors for both states_fixed and states
#disadvantage is that working with unmutables is less efficient
#also I believe it would be better to have both structs
# mutable struct State{T<:AbstractVector}
#     st::T
#     ope::Op
#     nume::Union{Int64,Missing}
# end
#
# function State(st,ope;nume=missing)
#     return State(st,ope,nume)
# end

# non_zero inputs a vector c and a precision (number of digits)
# and outputs two vectors: the first one indicating the indexes
# of non zero elements (up to precision prec) and the second
# the coefficient in each of these non zero indexes


text:
using Flux, CUDA, Test
using Flux: pullback

@testset for R in [RNN, GRU, LSTM]
  m = R(10, 5) |> gpu
  x = gpu(rand(10))
  (m̄,) = gradient(m -> sum(m(x)), m)
  Flux.reset!(m)
  θ = gradient(() -> sum(m(x)), params(m))
  @test x isa CuArray
  @test θ[m.cell.Wi] isa CuArray
  @test_broken collect(m̄[].cell[].Wi) == collect(θ[m.cell.Wi])
end

@testset "RNN" begin
  @testset for R in [RNN, GRU, LSTM], batch_size in (1, 5)
    rnn = R(10, 5)
    curnn = fmap(gpu, rnn)

    Flux.reset!(rnn)
    Flux.reset!(curnn)
    x = batch_size == 1 ?
      rand(Float32, 10) :
      rand(Float32, 10, batch_size)
    cux = gpu(x)

    y, back = pullback((r, x) -> r(x), rnn, x)
    cuy, cuback = pullback((r, x) -> r(x), curnn, cux)

    @test y ≈ collect(cuy)

    ȳ = randn(size(y))
    m̄, x̄ = back(ȳ)
    cum̄, cux̄ = cuback(gpu(ȳ))

    @test x̄ ≈ collect(cux̄)
    @test_broken m̄[].cell[].Wi ≈ collect(cum̄[].cell[].Wi)
    @test_broken m̄[].cell[].Wh ≈ collect(cum̄[].cell[].Wh)
    @test_broken m̄[].cell[].b ≈ collect(cum̄[].cell[].b)
    if m̄[].state isa Tuple
      for (x, cx) in zip(m̄[].state, cum̄[].state)
        @test x ≈ collect(cx)
      end
    else
      @test m̄[].state ≈ collect(cum̄[].state)
    end

    Flux.reset!(rnn)
    Flux.reset!(curnn)
    ohx = batch_size == 1 ?
      Flux.onehot(rand(1:10), 1:10) :
      Flux.onehotbatch(rand(1:10, batch_size), 1:10)
    cuohx = gpu(ohx)
    y = (rnn(ohx); rnn(ohx))

    cuy = (curnn(cuohx); curnn(cuohx))
    @test y ≈ collect(cuy)
  end
end


text:
#!/usr/bin/env julia

function puaHexStr(x)
 s = lpad( string(x, base=16), 2 , "0")
 "&#xF00$s;"
end

println.( puaHexStr.(0:0xff) )

max_stars_repo_name: sven-oly/LanguageTools

max_issues_repo_name: sven-oly/LanguageTools

max_forks_repo_name: sven-oly/LanguageTools


text:
using FilesLH
using Test

# ------------  Helpers

include("test_setup.jl")

# --------  Tests

include("file_transfer_test.jl")
include("computer_test.jl")
include("sbatch_test.jl")

@testset "FilesLH" begin
    file_transfer_test()
    computer_test()
    sbatch_test()
    include("directories_test.jl")
    include("file_names_test.jl")
end

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


text:
using SimpleAiboWebAPI
using Test

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


text:
"""
      SatelliteInertia

      Orientation represented with
      Modified Rodrigues Parameters.
"""

struct SatelliteInertia{I, T} <: Model{I, T}
      n::Int
      m::Int
      d::Int
end

function f(model::SatelliteInertia, z, u, w)
      # states
      r = view(z, 1:3)
      ω = view(z, 4:6)

      # controls
      τ = view(u, 1:3)
	  J = Diagonal(view(u, 4:6))

      SVector{6}([0.25 * ((1.0 - r' * r) * ω - 2.0 * cross(ω, r) + 2.0 * (ω' * r) * r);
                  J \ (τ - cross(ω, J * ω))])
end

function kinematics(model::SatelliteInertia, q)
	p = @SVector [0.1, 0.0, 0.0]
	k = MRP(view(q, 1:3)...) * p
	return k
end

# model
n, m, d = 6, 6, 0

model = SatelliteInertia{Midpoint, FixedTime}(n, m, d)

# visuals
function visualize!(vis, p::SatelliteInertia, q; Δt = 0.1)
	setvisible!(vis["/Background"], true)
	setprop!(vis["/Background"], "top_color", RGBA(0.0, 0.0, 0.0, 1.0))
	setprop!(vis["/Background"], "bottom_color", RGBA(0.0, 0.0, 0.0, 1.0))
	setvisible!(vis["/Axes"], false)
	setvisible!(vis["/Grid"], false)

    setobject!(vis[:satellite],
    	Rect(Vec(-0.25, -0.25, -0.25),Vec(0.5, 0.5, 0.5)),
    	MeshPhongMaterial(color = RGBA(1.0, 1.0, 1.0, 1.0)))

    arrow_x = ArrowVisualizer(vis[:satellite][:arrow_x])
    mat = MeshPhongMaterial(color=RGBA(1.0, 0.0, 0.0, 1.0))
    setobject!(arrow_x, mat)
    settransform!(arrow_x,
    	Point(0.0, 0.0, 0.0),
    	Vec(0.75, 0.0, 0.0),
    	shaft_radius=0.05,
    	max_head_radius=0.1)

    arrow_y = ArrowVisualizer(vis[:satellite][:arrow_y])
    mat = MeshPhongMaterial(color=RGBA(0.0, 1.0, 0.0, 1.0))
    setobject!(arrow_y, mat)
    settransform!(arrow_y,
    	Point(0.0, 0.0, 0.0),
    	Vec(0.0, 0.75, 0.0),
    	shaft_radius=0.05,
    	max_head_radius=0.1)

    arrow_z = ArrowVisualizer(vis[:satellite][:arrow_z])
    mat = MeshPhongMaterial(color=RGBA(0.0, 0.0, 1.0, 1.0))
    setobject!(arrow_z, mat)
    settransform!(arrow_z,
    	Point(0.0, 0.0, 0.0),
    	Vec(0.0, 0.0, 0.75),
    	shaft_radius=0.05,
    	max_head_radius=0.1)

    anim = MeshCat.Animation(convert(Int, floor(1.0 / Δt)))

     for t = 1:length(q)
    	 MeshCat.atframe(anim, t) do
    		 settransform!(vis["satellite"],
    			   compose(Translation(0.0 * [-0.25; -0.25; -0.25]...),
    					 LinearMap(MRP(q[t][1:3]...))))
    	 end
     end

     MeshCat.setanimation!(vis, anim)
end

max_issues_repo_name: thowell/DirectMotionPlanning


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Desription
# ==============================================================================
#
#   Tests related to conversion from Euler angles to quaternion.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# File: ./src/conversions/angle_to_quat.jl
# ========================================

# Functions: angle_to_quat
# ------------------------

@testset "Angle about an axis => Quaternion" begin
    for T in (Float32, Float64)
        # Sample an angle to use in all rotations.
        ang = _rand_ang(T)

        # Rotations
        # =========

        # X Axis
        # ------

        # Create a quaternion that rotates about X axis.
        q = angle_to_quat(ang, :X)
        @test eltype(q) === T

        # Create a vector that does not have X component.
        v = SVector(0, randn(T), randn(T))

        # Rotate the reference using the quaternion.
        v_r = vect(q \ v * q)
        @test eltype(v_r) === T

        # Get the sine of the angle between the representations.
        sin_ang = ((v_r × v) / norm(v)^2)[1]

        # Test the angle between the two representations.
        @test abs(sin(ang) - sin_ang) ≈ 0 atol = √(eps(T))

        # Y Axis
        # ------

        # Create a quaternion that rotates about Y axis.
        q = angle_to_quat(ang, :Y)
        @test eltype(q) === T

        # Create a vector that does not have Y component.
        v = SVector(randn(T), 0, randn(T))

        # Rotate the reference using the quaternion.
        v_r = vect(q \ v * q)
        @test eltype(v_r) === T

        # Get the sine of the angle between the representations.
        sin_ang = ((v_r × v) / norm(v)^2)[2]

        # Test the angle between the two representations.
        @test abs(sin(ang) - sin_ang) ≈ 0 atol = √(eps(T))

        # Z Axis
        # ------

        # Create a quaternion that rotates about Z axis.
        q = angle_to_quat(ang, :Z)
        @test eltype(q) === T

        # Create a vector that does not have Z component.
        v = SVector(randn(T), randn(T), 0)

        # Rotate the reference using the quaternion.
        v_r = vect(q \ v * q)
        @test eltype(v_r) === T

        # Get the sine of the angle between the representations.
        sin_ang = ((v_r × v) / norm(v)^2)[3]

        # Test the angle between the two representations.
        @test abs(sin(ang) - sin_ang) ≈ 0 atol = √(eps(T))
    end
end

@testset "Euler angles => Quaternion" begin
    for T in (Float32, Float64)
        # Two rotations
        # ======================================================================

        for rot_seq in valid_rot_seqs_2angles
            # Sample Euler angles.
            θ₁ = _rand_ang(T)
            θ₂ = _rand_ang(T)

            # Convert to DCM.
            q = angle_to_quat(θ₁, θ₂, rot_seq)
            @test eltype(q) === T

            # Create the DCM by composing the rotations.
            rot_seq_str = string(rot_seq)
            q₁ = angle_to_quat(θ₁, Symbol(rot_seq_str[1]))
            q₂ = angle_to_quat(θ₂, Symbol(rot_seq_str[2]))
            qe = q₁ * q₂

            # Compare.
            @test q ≈ qe
        end

        # Three rotations
        # ======================================================================

        for rot_seq in valid_rot_seqs
            # Sample Euler angles.
            ea = EulerAngles(_rand_ang(T), _rand_ang(T), _rand_ang(T), rot_seq)

            # Convert to quaternion.
            q = angle_to_quat(ea)
            @test eltype(q) === T

            # Create the quaternion by composing the rotations.
            rot_seq_str = string(rot_seq)

            q₁ = angle_to_quat(ea.a1, Symbol(rot_seq_str[1]))
            q₂ = angle_to_quat(ea.a2, Symbol(rot_seq_str[2]))
            q₃ = angle_to_quat(ea.a3, Symbol(rot_seq_str[3]))
            qe = q₁ * q₂ * q₃

            # Make sure q0 is positive.
            (qe.q0 < 0) && (qe = -qe)

            # Compare.
            @test q ≈ qe
        end
    end
end

@testset "Euler angles => Quaternion (Promotion)" begin
    # Check if promotion is working as intended.
    quat = angle_to_quat(Int64(1), 0.0f0, Float64(0))
    @test eltype(quat) === Float64

    quat = angle_to_quat(Int64(1), 0.0f0, 0.0f0)
    @test eltype(quat) === Float32

    quat = angle_to_quat(Int64(1), Int32(0), Float16(0))
    @test eltype(quat) === Float16
end

@testset "Euler angles => Quaternion (Errors)" begin
    @test_throws ArgumentError angle_to_quat(1, :V)
    @test_throws ArgumentError angle_to_quat(1, 2, :ZZ)
    @test_throws ArgumentError angle_to_quat(1, 2, 3, :ZZY)
end

# Functions: smallangle_to_quat
# -----------------------------

@testset "Small Euler angles => Quaternion" begin
    for T in (Float32, Float64)
        # Sample three small angles.
        θx = _rand_ang(T) * T(0.001)
        θy = _rand_ang(T) * T(0.001)
        θz = _rand_ang(T) * T(0.001)

        # Create the quaternion.
        q = smallangle_to_quat(θx, θy, θz)

        # Expected result.
        qe = Quaternion(1, θx / 2, θy / 2, θz / 2)
        qe = qe / norm(qe)

        @test qe ≈ q
    end
end

@testset "Small Euler angles => Quaternion (Promotion)" begin
    # Check if promotion is working as intended.
    quat = smallangle_to_quat(Int64(1), 0.0f0, Float64(0))
    @test eltype(quat) === Float64

    quat = smallangle_to_quat(Int64(1), 0.0f0, 0.0f0)
    @test eltype(quat) === Float32

    quat = smallangle_to_quat(Int64(1), Int32(0), Float16(0))
    @test eltype(quat) === Float16
end


text:
using DemoPackageMCF
using Documenter

DocMeta.setdocmeta!(DemoPackageMCF, :DocTestSetup, :(using DemoPackageMCF); recursive=true)

makedocs(;
    modules=[DemoPackageMCF],
    authors="mike freund <m.freund@wustl.edu> and contributors",
    repo="https://github.com/mcfreund/DemoPackageMCF.jl/blob/{commit}{path}#{line}",
    sitename="DemoPackageMCF.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://mcfreund.github.io/DemoPackageMCF.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/mcfreund/DemoPackageMCF.jl",
    devbranch="main",
)


text:
"""
Simple string-based CI code
"""
module StringCI
using Arpack
using TimerOutputs
using TensorOperations

include("Helpers.jl")
include("DeterminantStrings.jl")
include("FCI.jl")
include("TDMs.jl")

N = 30
binom_coeff = Array{Int,2}(undef,N+1,N+1)
for i in 0:N
    for j in i:N
        binom_coeff[j+1,i+1] = calc_nchk(j,i)
    end
end

end

max_stars_repo_name: nmayhall-vt/FermiCG

max_issues_repo_name: nmayhall-vt/FermiCG

max_forks_repo_name: nmayhall-vt/FermiCG


text:
using Sundials, Plots
using Test

pyplot(size=(700,700))

ProjDir = dirname(@__FILE__)
cd(ProjDir) do

  function sysfn(y, fy)
      fy[1] = y[1]^2 + y[2]^2 - 1.0
      fy[2] = y[2] - y[1]^2
  end

  res = Sundials.kinsol(sysfn, ones(2))
  res |> display
  println()

  @test round(res, 4) == round([0.786153, 0.618035], 4)

end


text:
using DFTK
using Unitful
using UnitfulAtomic
using JLD2
using PyCall
using MPI

if mpi_nprocs() == 1
    setup_threading()
else
    disable_threading()
end

# This tries to replicate the computational setup used for generating
# the materials project data as closely as possible. For details on
# the parameter choices see
#
#    https://docs.materialsproject.org/methodology/total-energies/
#
#
# All functions are MPI-ready, i.e. can be used through an MPI-parallelised julia session.

#
# If convergence issues on a system, try one of them (in this order and separately)
#    (a) is_metal=false
#    (b) use_experimental=false
#    (c) use_experimental=false, lower the damping (to 0.6, 0.4 and below if needed.)
#

function compute_formation_energy_per_atom(file::String; kwargs_bulk=(; ),
                                           kwargs_compound=(; ), kwargs...)
    scfres = run_if_needed(file; kwargs..., kwargs_compound...)

    bulk_energy_per_atom = Dict(map(scfres.basis.model.atoms) do (element, _)
        scfres_bulk = run_if_needed(element.symbol; kwargs..., kwargs_bulk...)
        n_atoms_bulk = sum(el_pos -> length(el_pos[2]), scfres_bulk.basis.model.atoms)
        element => scfres_bulk.energies.total / n_atoms_bulk
    end)

    atoms   = scfres.basis.model.atoms
    n_atoms = sum(el_pos -> length(el_pos[2]), atoms)
    energy_bulk = sum(el_pos -> bulk_energy_per_atom[el_pos[1]] * length(el_pos[2]), atoms)
    formation_energy_hartree = (scfres.energies.total - energy_bulk) / n_atoms
    auconvert(u"eV", formation_energy_hartree)
end

function run_if_needed(file::String; kwargs...)
    outfile = joinpath(@__DIR__, file * ".jld2")
    if !isfile(outfile)
        mpi_master() && println("#\n# -- $file\n#")
        scfres = run_from_file(file; kwargs...)
        store_results(scfres, outfile)
        scfres
    else
        load_results(outfile)
    end
end
function run_if_needed(element::Symbol; kwargs...)
    mkpath(joinpath(@__DIR__, "reference"))
    outfile = joinpath(@__DIR__, "reference", "$element.jld2")
    if !isfile(outfile)
        mpi_master() && println("#\n# -- $element\n#")
        scfres = run_bulk(element; kwargs...)
        store_results(scfres, outfile)
        scfres
    else
        load_results(outfile)
    end
end

function run_bulk(element::Symbol; kwargs...)
    bulkstructure = pyimport("ase.build").bulk(string(element))
    lattice = load_lattice(bulkstructure)
    atoms   = load_atoms(bulkstructure)
    run_dftk(lattice, atoms; kwargs...)
end

function run_from_file(file::AbstractString; kwargs...)
    lattice = load_lattice(file)
    atoms   = load_atoms(file)
    run_dftk(lattice, atoms; kwargs...)
end

# Store the results from `run_from_file`
function store_results(scfres, file)
    basis_master = DFTK.gather_kpts(scfres.basis)
    if mpi_master()
        JLD2.jldopen(file, "w") do jld
            jld["basis"]    = basis_master
            jld["ρ"]        = scfres.ρ
            jld["energies"] = scfres.energies
        end
    end
end

function load_results(file)
    JLD2.jldopen(file, "r") do jld
        (energies=jld["energies"], basis=jld["basis"], ρ=jld["ρ"])
    end
end

function run_dftk(lattice, atoms; use_experimental=true, is_metal=true, damping=0.8, debug=false, kwargs...)
    n_atoms = sum(el_pos -> length(el_pos[2]), atoms)

    # Check there are no atoms in the atoms list that require GGA+U
    has_ggau_metal = any(atoms) do (element, positions)
        element.symbol in (:V, :Cr, :Mn, :Fe, :Co, :Ni, :W, :Mo)
    end
    has_oxygen   = any(el_pos -> el_pos[1].symbol == :O, atoms)
    has_flouride = any(el_pos -> el_pos[1].symbol == :F, atoms)
    if (has_oxygen || has_flouride) && has_ggau_metal
        error("This compound requires GGA+U which is not implemented in DFTK.")
    end

    # Pseudos: Select the ones with largest number of electrons
    atoms = map(atoms) do (element, positions)
        symbol = element.symbol
        pspid = sort(list_psp(symbol, functional="pbe"), by=t->t.n_elec_valence)[end]
        ElementPsp(symbol, psp=load_psp(pspid.identifier)) => positions
    end

    magnetic_moments = map(atoms) do (element, positions)
        element => fill(default_magnetic_moment(element.symbol), length(positions))
    end
    # https://docs.materialsproject.org/methodology/total-energies/#calculation-details

    # They used the Tetrahedron method + Blöchl corrections ... which we don't have in DFTK.
    smearing = Smearing.Gaussian()
    temperature = 0.2u"eV"
    model = model_PBE(lattice, atoms; temperature, magnetic_moments, smearing)

    if debug
        kgrid = kgrid_from_minimal_n_kpoints(lattice, 6)
        Ecut = 10
    else
        kgrid = kgrid_from_minimal_n_kpoints(lattice, min(500, ceil(Int, 1000 / n_atoms)))
        Ecut  = 520u"eV"
    end
    basis = PlaneWaveBasis(model; Ecut, kgrid)

    ρ       = guess_density(basis, magnetic_moments)
    tol     = austrip(n_atoms * 5e-5u"eV")
    n_bands = DFTK.default_n_bands(basis.model)

    if use_experimental
        mixing    = is_metal ? KerkerMixing() : SimpleMixing()
        damping   = :adaptive
        algorithm = ("adaptive damping", DFTK.scf_potential_mixing_adaptive)
        extra     = (; )
    else
        mixing  = LdosMixing()
        extra = (damping=damping, )
        algorithm = ("fixed damping", self_consistent_field)
    end

    if mpi_master()
        magmom_compact = [element.symbol => values for (element, values) in magnetic_moments]
        display(basis)
        println()
        println()
        println("Solver parameters:")
        println("    magmom     : $magmom_compact")
        println("    tol        : $tol")
        println("    n_bands    : $n_bands")
        println("    algorithm  : $(algorithm[1])")
        println("    mixing     : $mixing")
        println("    damping    : $damping")
        println("    kwargs     : $(Dict(kwargs...))")
        println()
        flush(stdout)
    end

    DFTK.reset_timer!(DFTK.timer)
    run_scf = algorithm[2]
    scfres = run_scf(basis; ρ, tol, mixing, n_bands, maxiter=200, extra..., kwargs...)
    if mpi_master()
        println()
        println(DFTK.timer)
        println()
    end
    scfres
end

function default_magnetic_moment(element::Symbol)
    # TODO Consider adding this to PeriodicTable
    #      ... or DFTK ... or both
    get(Dict(
        :H  => 1.0,
        :He => 0.0,
        :Li => 1.0,
        :Be => 0.0,
        :B  => 1.0,
        :C  => 2.0,
        :N  => 3.0,
        :O  => 2.0,
        :F  => 1.0,
        :Ne => 0.0,
        :Na => 1.0,
        :Mg => 0.0,
        :Al => 1.0,
        :Si => 2.0,
        :P  => 3.0,
        :S  => 2.0,
        :Cl => 1.0,
        :Ar => 0.0,
        :K  => 1.0,
        :Ca => 0.0,
        :Sc => 1.0,
        :Ti => 2.0,
        :V  => 3.0,
        :Cr => 6.0,
        :Mn => 5.0,
        :Fe => 4.0,
        :Co => 3.0,
        :Ni => 2.0,
        :Cu => 1.0,
        :Zn => 0.0,
        :Ga => 1.0,
        :Ge => 2.0,
        :As => 3.0,
        :Se => 2.0,
        :Br => 1.0,
        :Kr => 0.0,
        :Rb => 1.0,
        :Sr => 0.0,
        :Y  => 1.0,
        :Zr => 2.0,
        :Nb => 5.0,
        :Mo => 6.0,
        :Tc => 5.0,
        :Ru => 4.0,
        :Rh => 3.0,
        :Pd => 0.0,
        :Ag => 1.0,
        :Cd => 0.0,
        :In => 1.0,
        :Sn => 2.0,
        :Sb => 3.0,
        :Te => 2.0,
        :I  => 1.0,
        :Xe => 0.0,
        :Cs => 1.0,
        :Ba => 0.0,
        :La => 1.0,
        :Ce => 1.0,
        :Pr => 3.0,
        :Nd => 4.0,
        :Pm => 5.0,
        :Sm => 6.0,
        :Eu => 7.0,
        :Gd => 8.0,
        :Tb => 5.0,
        :Dy => 4.0,
        :Ho => 3.0,
        :Er => 2.0,
        :Tm => 1.0,
        :Yb => 0.0,
        :Lu => 1.0,
        :Hf => 2.0,
        :Ta => 3.0,
        :W  => 4.0,
        :Re => 5.0,
        :Os => 4.0,
        :Ir => 3.0,
        :Pt => 2.0,
        :Au => 1.0,
        :Hg => 0.0,
        :Tl => 1.0,
        :Pb => 2.0,
        :Bi => 3.0,
        :Po => 2.0,
        :At => 1.0,
        :Rn => 0.0,
        :Fr => 1.0,
        :Ra => 0.0,
        :Ac => 1.0,
        :Th => 2.0,
        :Pa => 3.0,
        :U  => 4.0,
        :Np => 5.0,
        :Pu => 6.0,
        :Am => 7.0,
        :Cm => 8.0,
        :Bk => 5.0,
        :Cf => 4.0,
        :Es => 4.0,
        :Fm => 2.0,
        :Md => 1.0,
        :No => 0.0,
    ), element, 0)
end

max_stars_repo_name: mfherbst/DFTK-materialsproject

max_issues_repo_name: mfherbst/DFTK-materialsproject

max_forks_repo_name: mfherbst/DFTK-materialsproject


text:
# wrappers of the low-level dense CUSOLVER functionality

using ..CUBLAS: cublasfill, cublasop, cublasside

using LinearAlgebra
using LinearAlgebra: BlasInt, checksquare
using LinearAlgebra.LAPACK: chkargsok

# convert Char {N,V} to cusolverEigMode_t
function cusolverjob(jobz::Char)
    if jobz == 'N'
        return CUSOLVER_EIG_MODE_NOVECTOR
    end
    if jobz == 'V'
        return CUSOLVER_EIG_MODE_VECTOR
    end
    throw(ArgumentError("unknown cusolver eigmode $jobz."))
end

# po
## potrf
for (bname, fname,elty) in ((:cusolverDnSpotrf_bufferSize, :cusolverDnSpotrf, :Float32),
                            (:cusolverDnDpotrf_bufferSize, :cusolverDnDpotrf, :Float64),
                            (:cusolverDnCpotrf_bufferSize, :cusolverDnCpotrf, :ComplexF32),
                            (:cusolverDnZpotrf_bufferSize, :cusolverDnZpotrf, :ComplexF64))
    @eval begin
        function LinearAlgebra.LAPACK.potrf!(uplo::Char,
                        A::CuMatrix{$elty})
            cuuplo  = cublasfill(uplo)
            n       = checksquare(A)
            lda     = max(1, stride(A, 2))
            bufSize = Ref{Cint}(0)
            $bname(dense_handle(), cuuplo, n, A, lda, bufSize)

            buffer = CuArray{$elty}(undef, bufSize[])
            devinfo = CuArray{Cint}(undef, 1)
            $fname(dense_handle(), cuuplo, n, A, lda, buffer, bufSize[], devinfo)
            unsafe_free!(buffer)

            info = BlasInt(_getindex(devinfo, 1))
            unsafe_free!(devinfo)
            chkargsok(info)

            A, info
        end
    end
end

## potrs
for (fname,elty) in ((:cusolverDnSpotrs, :Float32),
                     (:cusolverDnDpotrs, :Float64),
                     (:cusolverDnCpotrs, :ComplexF32),
                     (:cusolverDnZpotrs, :ComplexF64))
    @eval begin
        function LinearAlgebra.LAPACK.potrs!(uplo::Char,
                        A::CuMatrix{$elty},
                        B::CuVecOrMat{$elty})
            cuuplo = cublasfill(uplo)
            n = checksquare(A)
            if size(B, 1) != n
                throw(DimensionMismatch("first dimension of B, $(size(B,1)), must match second dimension of A, $n"))
            end
            nrhs = size(B,2)
            lda  = max(1, stride(A, 2))
            ldb  = max(1, stride(B, 2))

            devinfo = CuArray{Cint}(undef, 1)
            $fname(dense_handle(), cuuplo, n, nrhs, A, lda, B, ldb, devinfo)

            info = _getindex(devinfo, 1)
            unsafe_free!(devinfo)
            chkargsok(BlasInt(info))

            B
        end
    end
end

#getrf
for (bname, fname,elty) in ((:cusolverDnSgetrf_bufferSize, :cusolverDnSgetrf, :Float32),
                            (:cusolverDnDgetrf_bufferSize, :cusolverDnDgetrf, :Float64),
                            (:cusolverDnCgetrf_bufferSize, :cusolverDnCgetrf, :ComplexF32),
                            (:cusolverDnZgetrf_bufferSize, :cusolverDnZgetrf, :ComplexF64))
    @eval begin
        function getrf!(A::CuMatrix{$elty})
            m,n     = size(A)
            lda     = max(1, stride(A, 2))
            bufSize = Ref{Cint}(0)
            $bname(dense_handle(), m, n, A, lda, bufSize)

            buffer  = CuArray{$elty}(undef, bufSize[])
            devipiv = CuArray{Cint}(undef, min(m,n))
            devinfo = CuArray{Cint}(undef, 1)
            $fname(dense_handle(), m, n, A, lda, buffer, devipiv, devinfo)
            unsafe_free!(buffer)

            info = _getindex(devinfo, 1)
            unsafe_free!(devinfo)
            if info < 0
                throw(ArgumentError("The $(info)th parameter is wrong"))
            elseif info > 0
                throw(LinearAlgebra.SingularException(info))
            end

            A, devipiv
        end
    end
end

#geqrf
for (bname, fname,elty) in ((:cusolverDnSgeqrf_bufferSize, :cusolverDnSgeqrf, :Float32),
                            (:cusolverDnDgeqrf_bufferSize, :cusolverDnDgeqrf, :Float64),
                            (:cusolverDnCgeqrf_bufferSize, :cusolverDnCgeqrf, :ComplexF32),
                            (:cusolverDnZgeqrf_bufferSize, :cusolverDnZgeqrf, :ComplexF64))
    @eval begin
        function geqrf!(A::CuMatrix{$elty})
            m, n    = size(A)
            lda     = max(1, stride(A, 2))
            bufSize = Ref{Cint}(0)
            $bname(dense_handle(), m, n, A, lda, bufSize)

            buffer  = CuArray{$elty}(undef, bufSize[])
            tau  = CuArray{$elty}(undef, min(m, n))
            devinfo = CuArray{Cint}(undef, 1)
            $fname(dense_handle(), m, n, A, lda, tau, buffer, bufSize[], devinfo)
            unsafe_free!(buffer)

            info = _getindex(devinfo, 1)
            unsafe_free!(devinfo)
            if info < 0
                throw(ArgumentError("The $(info)th parameter is wrong"))
            end

            A, tau
        end
    end
end

#sytrf
for (bname, fname,elty) in ((:cusolverDnSsytrf_bufferSize, :cusolverDnSsytrf, :Float32),
                            (:cusolverDnDsytrf_bufferSize, :cusolverDnDsytrf, :Float64),
                            (:cusolverDnCsytrf_bufferSize, :cusolverDnCsytrf, :ComplexF32),
                            (:cusolverDnZsytrf_bufferSize, :cusolverDnZsytrf, :ComplexF64))
    @eval begin
        function sytrf!(uplo::Char,
                        A::CuMatrix{$elty})
            cuuplo = cublasfill(uplo)
            n      = checksquare(A)
            lda = max(1, stride(A, 2))
            bufSize = Ref{Cint}(0)
            $bname(dense_handle(), n, A, lda, bufSize)

            buffer  = CuArray{$elty}(undef, bufSize[])
            devipiv = CuArray{Cint}(undef, n)
            devinfo = CuArray{Cint}(undef, 1)
            $fname(dense_handle(), cuuplo, n, A, lda, devipiv, buffer, bufSize[], devinfo)
            unsafe_free!(buffer)

            info = _getindex(devinfo, 1)
            unsafe_free!(devinfo)
            if info < 0
                throw(ArgumentError("The $(info)th parameter is wrong"))
            elseif info > 0
                throw(LinearAlgebra.SingularException(info))
            end

            A, devipiv
        end
    end
end

#getrs
for (fname,elty) in ((:cusolverDnSgetrs, :Float32),
                     (:cusolverDnDgetrs, :Float64),
                     (:cusolverDnCgetrs, :ComplexF32),
                     (:cusolverDnZgetrs, :ComplexF64))
    @eval begin
        function getrs!(trans::Char,
                        A::CuMatrix{$elty},
                        ipiv::CuVector{Cint},
                        B::CuVecOrMat{$elty})
            cutrans = cublasop(trans)
            n = size(A, 1)
            if size(A, 2) != n
                throw(DimensionMismatch("LU factored matrix A must be square!"))
            end
            if size(B, 1) != n
                throw(DimensionMismatch("first dimension of B, $(size(B,1)), must match second dimension of A, $n"))
            end
            nrhs = size(B, 2)
            lda  = max(1, stride(A, 2))
            ldb  = max(1, stride(B, 2))

            devinfo = CuArray{Cint}(undef, 1)
            $fname(dense_handle(), cutrans, n, nrhs, A, lda, ipiv, B, ldb, devinfo)

            info = _getindex(devinfo, 1)
            unsafe_free!(devinfo)
            if info < 0
                throw(ArgumentError("The $(info)th parameter is wrong"))
            end
            B
        end
    end
end

#ormqr
for (bname, fname, elty) in ((:cusolverDnSormqr_bufferSize, :cusolverDnSormqr, :Float32),
                             (:cusolverDnDormqr_bufferSize, :cusolverDnDormqr, :Float64),
                             (:cusolverDnCunmqr_bufferSize, :cusolverDnCunmqr, :ComplexF32),
                             (:cusolverDnZunmqr_bufferSize, :cusolverDnZunmqr, :ComplexF64))
    @eval begin
        function ormqr!(side::Char,
                        trans::Char,
                        A::CuMatrix{$elty},
                        tau::CuVector{$elty},
                        C::CuVecOrMat{$elty})
            cutrans = cublasop(trans)
            cuside  = cublasside(side)
            if side == 'L'
                m   = size(A, 1)
                ldc = size(C, 1)
                n   = size(C, 2)
                if m > ldc
                    Ctemp = CuArrays.zeros($elty, m - ldc, n)
                    C = [C; Ctemp]
                    ldc = m
                end
                lda = m
            else
                m   = size(C, 1)
                n   = size(C, 2)
                ldc = m
                lda = n
            end
            k       = length(tau)

            bufSize = Ref{Cint}(0)
            $bname(dense_handle(), cuside, cutrans, m, n, k, A, lda, tau, C, ldc, bufSize)

            buffer  = CuArray{$elty}(undef, bufSize[])
            devinfo = CuArray{Cint}(undef, 1)
            $fname(dense_handle(), cuside, cutrans, m, n, k, A, lda, tau, C, ldc, buffer,
                   bufSize[], devinfo)
            unsafe_free!(buffer)

            info = _getindex(devinfo, 1)
            unsafe_free!(devinfo)
            if info < 0
                throw(ArgumentError("The $(info)th parameter is wrong"))
            end

            side == 'L' ? C : C[:, 1:minimum(size(A))]
        end
    end
end

#orgqr
for (bname, fname, elty) in ((:cusolverDnSorgqr_bufferSize, :cusolverDnSorgqr, :Float32),
                             (:cusolverDnDorgqr_bufferSize, :cusolverDnDorgqr, :Float64),
                             (:cusolverDnCungqr_bufferSize, :cusolverDnCungqr, :ComplexF32),
                             (:cusolverDnZungqr_bufferSize, :cusolverDnZungqr, :ComplexF64))
    @eval begin
        function orgqr!(A::CuMatrix{$elty}, tau::CuVector{$elty})
            m = size(A , 1)
            n = min(m, size(A, 2))
            lda = max(1, stride(A, 2))
            k = length(tau)

            bufSize = Ref{Cint}(0)
            $bname(dense_handle(), m, n, k, A, lda, tau, bufSize)

            buffer  = CuArray{$elty}(undef, bufSize[])
            devinfo = CuArray{Cint}(undef, 1)
            $fname(dense_handle(), m, n, k, A, lda, tau, buffer, bufSize[], devinfo)
            unsafe_free!(buffer)

            info = _getindex(devinfo, 1)
            unsafe_free!(devinfo)
            if info < 0
                throw(ArgumentError("The $(info)th parameter is wrong"))
            end

            if n < size(A, 2)
                A[:, 1:n]
            else
                A
            end
        end
    end
end

#gebrd
for (bname, fname, elty, relty) in ((:cusolverDnSgebrd_bufferSize, :cusolverDnSgebrd, :Float32, :Float32),
                                    (:cusolverDnDgebrd_bufferSize, :cusolverDnDgebrd, :Float64, :Float64),
                                    (:cusolverDnCgebrd_bufferSize, :cusolverDnCgebrd, :ComplexF32, :Float32),
                                    (:cusolverDnZgebrd_bufferSize, :cusolverDnZgebrd, :ComplexF64, :Float64))
    @eval begin
        function gebrd!(A::CuMatrix{$elty})
            m, n    = size(A)
            lda     = max(1, stride(A, 2))

            bufSize = Ref{Cint}(0)
            $bname(dense_handle(), m, n, bufSize)

            buffer  = CuArray{$elty}(undef, bufSize[])
            devinfo = CuArray{Cint}(undef, 1)
            k       = min(m, n)
            D       = CuArray{$relty}(undef, k)
            E       = CuArrays.zeros($relty, k)
            TAUQ    = CuArray{$elty}(undef, k)
            TAUP    = CuArray{$elty}(undef, k)
            $fname(dense_handle(), m, n, A, lda, D, E, TAUQ, TAUP, buffer, bufSize[], devinfo)
            unsafe_free!(buffer)

            info = _getindex(devinfo, 1)
            unsafe_free!(devinfo)
            if info < 0
                throw(ArgumentError("The $(info)th parameter is wrong"))
            end

            A, D, E, TAUQ, TAUP
        end
    end
end

for (bname, fname, elty, relty) in ((:cusolverDnSgesvd_bufferSize, :cusolverDnSgesvd, :Float32, :Float32),
                                    (:cusolverDnDgesvd_bufferSize, :cusolverDnDgesvd, :Float64, :Float64),
                                    (:cusolverDnCgesvd_bufferSize, :cusolverDnCgesvd, :ComplexF32, :Float32),
                                    (:cusolverDnZgesvd_bufferSize, :cusolverDnZgesvd, :ComplexF64, :Float64))
    @eval begin
        function gesvd!(jobu::Char,
                        jobvt::Char,
                        A::CuMatrix{$elty})
            m, n    = size(A)
            if m < n
                throw(ArgumentError("CUSOLVER's gesvd currently requires m >= n"))
            end
            lda     = max(1, stride(A, 2))

            lwork   = Ref{Cint}(0)
            $bname(dense_handle(), m, n, lwork)

            if jobu === 'S' && m > n
                U   = CuArray{$elty}(undef, m, n)
            elseif jobu === 'N'
                U   = CuArray{$elty}(undef, m, 0)
            else
                U   = CuArray{$elty}(undef, m, m)
            end
            ldu     = max(1, stride(U, 2))

            S       = CuArray{$relty}(undef, min(m, n))

            if jobvt === 'S' && m < n
                Vt  = CuArray{$elty}(undef, m, n)
            elseif jobvt === 'N'
                Vt  = CuArray{$elty}(undef, 0, n)
            else
                Vt  = CuArray{$elty}(undef, n, n)
            end
            ldvt    = max(1, stride(Vt, 2))

            work    = CuArray{$elty}(undef, lwork[])
            rwork   = CuArray{$relty}(undef, min(m, n) - 1)
            devinfo = CuArray{Cint}(undef, 1)
            $fname(dense_handle(), jobu, jobvt, m, n, A, lda, S, U, ldu, Vt, ldvt,
                   work, lwork[], rwork, devinfo)
            unsafe_free!(work)
            unsafe_free!(rwork)

            info = _getindex(devinfo, 1)
            unsafe_free!(devinfo)
            if info < 0
                throw(ArgumentError("The $(info)th parameter is wrong"))
            end

            return U, S, Vt
        end
    end
end

for (bname, fname, elty, relty) in ((:cusolverDnSgesvdj_bufferSize, :cusolverDnSgesvdj, :Float32, :Float32),
                                    (:cusolverDnDgesvdj_bufferSize, :cusolverDnDgesvdj, :Float64, :Float64),
                                    (:cusolverDnCgesvdj_bufferSize, :cusolverDnCgesvdj, :ComplexF32, :Float32),
                                    (:cusolverDnZgesvdj_bufferSize, :cusolverDnZgesvdj, :ComplexF64, :Float64))
    @eval begin
        function gesvdj!(jobz::Char,
                         econ::Int,
                         A::CuMatrix{$elty};
                         tol::$relty=eps($relty),
                         max_sweeps::Int=100)
            cujobz  = cusolverjob(jobz)
            m,n     = size(A)
            lda     = max(1, stride(A, 2))

            # Warning! For some reason, the solver needs to access U and V even
            # when only the values are requested
            if jobz === 'V' && econ == 1 && m > n
                U   = CuArray{$elty}(undef, m, n)
            else
                U   = CuArray{$elty}(undef, m, m)
            end
            ldu     = max(1, stride(U, 2))

            S       = CuArray{$relty}(undef, min(m, n))

            if jobz === 'V' && econ == 1 && m < n
                V   = CuArray{$elty}(undef, n, m)
            else
                V   = CuArray{$elty}(undef, n, n)
            end
            ldv     = max(1, stride(V, 2))

            params  = Ref{gesvdjInfo_t}(C_NULL)
            cusolverDnCreateGesvdjInfo(params)
            cusolverDnXgesvdjSetTolerance(params[], tol)
            cusolverDnXgesvdjSetMaxSweeps(params[], max_sweeps)

            lwork   = Ref{Cint}(0)
            $bname(dense_handle(), cujobz, Cint(econ), m, n, A, lda, S, U, ldu, V, ldv,
                   lwork, params[])

            work  = CuArray{$elty}(undef, lwork[])
            devinfo = CuArray{Cint}(undef, 1)
            $fname(dense_handle(), cujobz, econ, m, n, A, lda, S, U, ldu, V, ldv,
                   work, lwork[], devinfo, params[])
            unsafe_free!(work)

            info = _getindex(devinfo, 1)
            unsafe_free!(devinfo)
            if info < 0
                throw(ArgumentError("The $(info)th parameter is wrong"))
            end

            cusolverDnDestroyGesvdjInfo(params[])

            U, S, V
        end
    end
end

for (jname, bname, fname, elty, relty) in ((:syevd!, :cusolverDnSsyevd_bufferSize, :cusolverDnSsyevd, :Float32, :Float32),
                                           (:syevd!, :cusolverDnDsyevd_bufferSize, :cusolverDnDsyevd, :Float64, :Float64),
                                           (:heevd!, :cusolverDnCheevd_bufferSize, :cusolverDnCheevd, :ComplexF32, :Float32),
                                           (:heevd!, :cusolverDnZheevd_bufferSize, :cusolverDnZheevd, :ComplexF64, :Float64))
    @eval begin
        function $jname(jobz::Char,
                        uplo::Char,
                        A::CuMatrix{$elty})
            cuuplo  = cublasfill(uplo)
            cujobz  = cusolverjob(jobz)
            n       = checksquare(A)
            lda     = max(1, stride(A, 2))
            W       = CuArray{$relty}(undef, n)

            bufSize = Ref{Cint}(0)
            $bname(dense_handle(), cujobz, cuuplo, n, A, lda, W, bufSize)

            buffer  = CuArray{$elty}(undef, bufSize[])
            devinfo = CuArray{Cint}(undef, 1)
            $fname(dense_handle(), cujobz, cuuplo, n, A, lda, W,
                   buffer, bufSize[], devinfo)
            unsafe_free!(buffer)

            info = _getindex(devinfo, 1)
            unsafe_free!(devinfo)
            if info < 0
                throw(ArgumentError("The $(info)th parameter is wrong"))
            end

            if jobz == 'N'
                return W
            elseif jobz == 'V'
                return W, A
            end
        end
    end
end

for (jname, bname, fname, elty, relty) in ((:sygvd!, :cusolverDnSsygvd_bufferSize, :cusolverDnSsygvd, :Float32, :Float32),
                                           (:sygvd!, :cusolverDnDsygvd_bufferSize, :cusolverDnDsygvd, :Float64, :Float64),
                                           (:hegvd!, :cusolverDnChegvd_bufferSize, :cusolverDnChegvd, :ComplexF32, :Float32),
                                           (:hegvd!, :cusolverDnZhegvd_bufferSize, :cusolverDnZhegvd, :ComplexF64, :Float64))
    @eval begin
        function $jname(itype::Int,
                        jobz::Char,
                        uplo::Char,
                        A::CuMatrix{$elty},
                        B::CuMatrix{$elty})
            cuuplo  = cublasfill(uplo)
            cujobz  = cusolverjob(jobz)
            nA, nB  = checksquare(A, B)
            if nB != nA
                throw(DimensionMismatch("Dimensions of A ($nA, $nA) and B ($nB, $nB) must match!"))
            end
            n       = nA
            lda     = max(1, stride(A, 2))
            ldb     = max(1, stride(B, 2))
            W       = CuArray{$relty}(undef, n)

            bufSize = Ref{Cint}(0)
            cuitype = cusolverEigType_t(itype)
            $bname(dense_handle(), cuitype, cujobz, cuuplo, n, A, lda, B, ldb, W, bufSize)

            buffer  = CuArray{$elty}(undef, bufSize[])
            devinfo = CuArray{Cint}(undef, 1)
            $fname(dense_handle(), cuitype, cujobz, cuuplo, n, A, lda, B, ldb, W,
                   buffer, bufSize[], devinfo)
            unsafe_free!(buffer)

            info = _getindex(devinfo, 1)
            unsafe_free!(devinfo)
            if info < 0
                throw(ArgumentError("The $(info)th parameter is wrong"))
            end

            if jobz == 'N'
                return W
            elseif jobz == 'V'
                return W, A, B
            end
        end
    end
end

for (jname, bname, fname, elty, relty) in ((:sygvj!, :cusolverDnSsygvj_bufferSize, :cusolverDnSsygvj, :Float32, :Float32),
                                           (:sygvj!, :cusolverDnDsygvj_bufferSize, :cusolverDnDsygvj, :Float64, :Float64),
                                           (:hegvj!, :cusolverDnChegvj_bufferSize, :cusolverDnChegvj, :ComplexF32, :Float32),
                                           (:hegvj!, :cusolverDnZhegvj_bufferSize, :cusolverDnZhegvj, :ComplexF64, :Float64))
    @eval begin
        function $jname(itype::Int,
                        jobz::Char,
                        uplo::Char,
                        A::CuMatrix{$elty},
                        B::CuMatrix{$elty};
                        tol::$relty=eps($relty),
                        max_sweeps::Int=100)
            cujobz  = cusolverjob(jobz)
            cuuplo  = cublasfill(uplo)
            nA, nB  = checksquare(A, B)
            if nB != nA
                throw(DimensionMismatch("Dimensions of A ($nA, $nA) and B ($nB, $nB) must match!"))
            end
            n       = nA
            lda     = max(1, stride(A, 2))
            ldb     = max(1, stride(B, 2))
            W       = CuArray{$relty}(undef, n)
            params  = Ref{syevjInfo_t}(C_NULL)
            cusolverDnCreateSyevjInfo(params)
            cusolverDnXsyevjSetTolerance(params[], tol)
            cusolverDnXsyevjSetMaxSweeps(params[], max_sweeps)

            bufSize = Ref{Cint}(0)
            cuitype = cusolverEigType_t(itype)
            $bname(dense_handle(), cuitype, cujobz, cuuplo, n, A, lda, B, ldb, W,
                   bufSize, params[])

            buffer  = CuArray{$elty}(undef, bufSize[])
            devinfo = CuArray{Cint}(undef, 1)
            $fname(dense_handle(), cuitype, cujobz, cuuplo, n, A, lda, B, ldb, W,
                   buffer, bufSize[], devinfo, params[])
            unsafe_free!(buffer)

            info = _getindex(devinfo, 1)
            unsafe_free!(devinfo)
            if info < 0
                throw(ArgumentError("The $(info)th parameter is wrong"))
            end

            cusolverDnDestroySyevjInfo(params[])

            if jobz == 'N'
                return W
            elseif jobz == 'V'
                return W, A, B
            end
        end
    end
end


text:
######## Hybrid kernel polynomial method ########
# Functions:

# fraunhofer
#     methods fullKPM, hybridKPM, exact

# Here we compute: <Ic> = <Ic>KPM - <Ic>KPM_proj_ABS + Ic_exact_ABS
#                  <Ic> = Ic_KPM_continuum + Ic_exact_ABS

const σ0τ0 = @SMatrix[1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1]
const σ0τz = @SMatrix[1 0 0 0; 0 1 0 0; 0 0 -1 0; 0 0 0 -1]
const ħoe = ustrip(u"eV/nA",ħ/e)
const a0 = 0.246

fraunhofer_kpm(Bperprange::StepRangeLen, p; Δϕpoints = 3, kw...) =
    fraunhofer_kpm(collect(Bperprange), p, collect(-1.:2/Δϕpoints:1); kw...)

"""
"the adaptive maximum value search for CPRs is enabled with the kw `adaptive`.
Disabled by default
"""
fraunhofer_kpm(Bperp::Array{T,1}, p; kw...) where {T} =
    fraunhofer_kpm(Bperp, p, missing; kw...)
function fraunhofer_kpm(Bperp::Array{T,1}, p, Δϕ::Union{Array{T,1}, Missing};
    adaptive = false, abssubspace = missing, kw...) where {T}
    println("computing current matrix...")
    Bperp, ic = icϕ_kpm(Bperp, p, Δϕ, abssubspace; adaptive = adaptive, kw...)
    if adaptive == false
        icmax = similar(Bperp)
        for i in 1:length(Bperp)
            println("it: ", i/length(Bperp))
            icmax[i] = maximum(abs.(ic[i,:]))
        end
        return Bperp, icmax, ic
    else
        return Bperp, ic, abs.(ic)
    end
end

"computes the fraunhofer pattern only considering the states inside the parent gap using
exact diagonalization"
fraunhofer_abs_exact(Bperp::Array{T,1}, p; kw...) where {T} =
    fraunhofer_abs_exact(Bperp, p, missing; kw...)
function fraunhofer_abs_exact(Bperp::Array{T,1}, p,
        Δϕ::Union{Array{T,1}, Missing}; kw...) where {T}
    println("computing current matrix...")
    icmax = SharedArray(similar(Bperp))
    ic = icϕ_exactdiag(Bperp, p, Δϕ; kw...)
    for i in 1:length(Bperp)
        println("it: ", i/length(Bperp))
        icmax[i] = maximum(abs.(ic[:, i]))
    end
    return Bperp, icmax, ic
end

"""
    `icϕ_exactdiag(Bperplist::Array{T,1}, p, Δϕ::Union{Array{T,1}, Missing}; kw...)`
supercurrent sweep with both Bperp and Δϕ. See: `supercurrent_exactdiag()`
"""
function icϕ_exactdiag(Bperplist::Array{T,1}, p,
    Δϕ::Union{Array{T,1}, Missing}; kw...) where {T}
    Δϕlist = ifelse(isa(Δϕ, Missing) == true, -0.5:0.25:π+0.5, Δϕ)
    I = zeros(Float64, length(Δϕlist), length(Bperplist))
    [I[:, i] = supercurrent_exactdiag(collect(Δϕlist),
        reconstruct(p, B_perp = Bperplist[i]); kw...)
            for i in 1:length(Bperplist)]
    return I
end

"""
    `supercurrent_exactdiag(Δϕlist, p = Params() ; nev = 10)`
Computes the supercurrent for a fixed Bperp value specified in p. If `nev::Missing` it
iteratively finds all ABS (i.e. with energies inside a lengthscale which we set us the gap).
It returns the supercurrent contribution of these subgap subspace. Note that this subspace
is chosen to always capture all the ABS contribution although it may contain also continuum
states, because the parent gap will be larger that the induced gap. However, once the
subspace is specified, we can always substract this contribution for the full KPM
contribution.
    Two different methods implemented
        `method = :free_energy`
        `method = :current_cut` (`sum(j_ik)` along a cut)
"""
function supercurrent_exactdiag(Δϕlist, p = Params(); nev = 10,
    method = :free_energy, kw...)
    println("we are using the following method: ", method)
    @unpack Δ = p
    hpar, = paramhams(p; flat = false)
    if method == :free_energy
        f = SharedArray(zeros(Float64, length(Δϕlist)))
        @sync @distributed for i in 1:length(Δϕlist)
            println(i/length(Δϕlist))
            f[i] = f_e(hpar, Δϕlist[i], nev, Δ; kw...)
        end
        ip = interpolate((Δϕlist,), f, Gridded(Linear()));
        deriv =[Interpolations.gradient(ip, i)[1] for i in Δϕlist]
        return @. 2/ħoe .* deriv
    else
        I = SharedArray(zeros(Float64, length(Δϕlist)))
        @sync @distributed for i in 1:length(Δϕlist)
            println(i/length(Δϕlist))
             st = cut_st(hpar, Δϕlist[i], nev, Δ)
            I[i] = supercurrent_cut(p, st, hpar())
        end
        return I
    end
end

"computes the supercurrent matrix (Bperp, Δϕ) if adaptive = false or the critical
supercurrent vs magnetic field if adaptive = true, using a KPM routine
see: adaptive_max_finder documentation. If adaptive is true and a ketmodel is given
for the ABSsubspace it will use the hybridKPM approach, otherwise it will
run a  KPM calculation over the whole spectrum"
function icϕ_kpm(Bperp::Array{T,1}, p, Δϕ::Union{Array{T,1}, Missing},
        abssubspace::Union{Quantica.KetModel, Missing};
        adaptive = false, rk = 10, kw...) where {T}
    # method = (isa(abssubspace, Missing) ? :fullKPM, :hybridKPM)
    if adaptive == false
        if isa(Δϕ, Missing)
            @warn "provide a valid Δϕ range"
        else
            # ic = SharedArray(Array{Float64, 2}(undef, length(Bperp), length(Δϕ)))
            # @sync @distributed for j in 1:length(Δϕ)
            ic = Array{Float64, 2}(undef, length(Bperp), length(Δϕ))
            for j in 1:length(Δϕ)
                println("phase: ", Δϕ[j])
                [ic[i, j] = supercurrent_kpm(rk, abssubspace,
                    reconstruct(p, B_perp = Bperp[i]), Δϕ = Δϕ[j]; kw...)[2]
                 for i in 1:length(Bperp)]
            end
        end
        return Bperp, ic
    else
        icmax = SharedArray(similar(Bperp))
        # local definition of generating functions for the adaptive calculation path
        gen_model(x::Number; it = 1, kw...) =
            supercurrent_kpm(rk, abssubspace,
                reconstruct(p, B_perp = Bperp[it]), Δϕ = x; kw...)[2]
            # /rk comes from normalization of the stochastic trace
        gen_model(x::Array; kw...) =
            [gen_model(x[i]; kw ...) for i in 1:length(x)]
        @sync @distributed for i in 1:length(Bperp)
            @unpack bandrange = kw
            brange = _bandrange(reconstruct(p, B_perp = Bperp[i]), bandrange)
            kw1 = [key=> ifelse(key === :bandrange, brange, value) for (key, value) in kw]
            # update kw with bandrange kwarg so there is no need to compute it for each ϕ,
            # it only varies with Bperp
            icmax[i] = adaptive_max_finder(gen_model, fourier_model, it = i; kw1...)[1]
        end
        # I could have used pmap for parallelization but since each task requires
        # the same effort
        return Bperp, icmax
    end
end

"""
it should be called icϕ_continuum_kpm
"""
function icϕ_hybrid(Bperp::Array{T,1}, p, Δϕ::Union{Array{T,1}, Missing};
    adaptive = false, nev = 10, rk = 10, bandrange = missing, kw...) where {T}
    method =(isa(nev, Missing) ? :adaptive_nev : :nev)
    @unpack Δ = p
    if adaptive == false
        if isa(Δϕ, Missing)
            @warn "provide a valid Δϕ range"
        else
            ic_continuum = Array{Float64, 2}(undef, length(Bperp),
            length(Δϕ))
            @sync @distributed for i in 1:length(Bperp)
                p = reconstruct(p, B_perp = Bperp[i])
                hpar, = paramhams(p; flat = false)
                brange = _bandrange(p, bandrange)
                for j in 1:length(Δϕ)
                    println("phase: ", Δϕ[j])
                    st = eigen_st(hpar, Δϕ[j], nev, Δ,
                        method = method)[2]
                    abssubspace = subspace_kets(st, hpar())
                    ic_continuum[i, j] =
                        supercurrent_kpm(rk, abssubspace, p, Δϕ = Δϕ[j],
                        bandrange = brange; kw...)[2]
                end
            end
        end
    else

        icmax = SharedArray(similar(Bperp))
        # local definition of generating functions for the adaptive
        # calculation path
        gen_model(x::Number; it = 1, kw...) =
            supercurrent_kpm(rk, abssubspace,
                reconstruct(p, B_perp = Bperp[it]), Δϕ = x; kw...)[2]
            # /rk comes from normalization of the stochastic trace
        gen_model(x::Array; kw...) =
            [gen_model(x[i]; kw ...) for i in 1:length(x)]
        @sync @distributed for i in 1:length(Bperp)
            @unpack bandrange = kw
            brange = _bandrange(reconstruct(p, B_perp = Bperp[i]),
                bandrange)
            kw1 = [key=> ifelse(key === :bandrange, brange, value)
                for (key, value) in kw]
            # update kw with bandrange kwarg so there is no need
            # to compute it for each ϕ,
            # it only varies with Bperp
            icmax[i] = adaptive_max_finder(gen_model, fourier_model,
            it = i; kw1...)[1]
        end
        # I could have used pmap for parallelization but since each
        # task requires the same effort
        # it wouldn't
        return Bperp, icmax
    end

    return Bperp, ic_continuum
end

"""
Computes the supercurrent for a hamiltonian built with parameters `p::Params` using `rk`
randomvectors with non-zero values only in the right lead (this results in a better
stochastic trace estimation with less resources). If `rk` given by the user do not ensure
the convergence it will be overwritten see: stochastictracecriteria. It uses distributed
with a loop over the phases instead of intel mklsparse
"""
function supercurrent_kpm(rk, abssubspace, p = Params(); Δϕ = 1,
    bandrange = missing, kw...)
    @unpack Δ, Ls, scale = p
    # bandranges = _bandrange(p, bandrange)
    hfunc, dhfunc = paramhams(p; flat = false)
    # rk = criteria(p, rk, dhfunc)
    if isa(abssubspace, Missing)
        return [Δϕ, av(π*Δϕ, Δ, bandrange, hfunc, dhfunc; kw...)]
    else
        return [Δϕ,
        av(π*Δϕ, Δ, bandrange, hfunc, dhfunc, rk = rk; kw...) -
        av(π*Δϕ, Δ, bandrange, hfunc, dhfunc, abssubspace; kw...)]
    end
end

"""
computes the thermodynamic average of the current operator using KPM. If
abssubspace is not an argument it computes the whole KPM, else it
substracts the contribution from ABS states
"""
av(ϕ, Δ, bandranges, hfunc, dhfunc; rk = 10, kw...) = 2/ħoe *
    real(averageKPM(hfunc(phi = ϕ),dhfunc(phi = ϕ); flat = Val(false),
        ket = ketsmodel(rk), bandrange = bandranges, kw...))/rk
av(ϕ, Δ, bandranges, hfunc, dhfunc, abssubspace::Vector; kw...) =
    sum(x -> av(ϕ, Δ, bandranges, hfunc, dhfunc, x; kw...), abssubspace)
av(ϕ, Δ, bandranges, hfunc, dhfunc, abssubspace::Quantica.Ket; kw...) =
    2/ħoe * real(averageKPM(hfunc(phi = ϕ),dhfunc(phi = ϕ);
    flat = Val(false), ket = abssubspace, bandrange = bandranges, kw...))

"""
random ket model for the stochastic evaluation of the two KPM traces
"""
ketsmodel(numkets) = randomkets(numkets,  r -> cis(2π * rand()) * σ0τ0)
"""
    subspace_kets(sp::Quantica.Spectrum, h::Quantica.Hamiltonian)
builds a ket array of the states contained in spectrum that should
be filtered out in a hybrid KPM calculation.
"""
subspace_kets(sp::Quantica.Spectrum, h::Quantica.Hamiltonian) =
    subspace_kets(sp.states, h)
function subspace_kets(states::Matrix{ComplexF64},
    h::Quantica.Hamiltonian)
    st = [reshape(states[:,1], 4, :) for i in 1:size(states,2)]
    ket_0 = ket(h)
    basis_ket = [ket_0 for i in 1:size(st)[1]]
    [basis_ket[i].amplitudes[j]= st[i][:,j] for i in 1:length(basis_ket)
    for j in 1:length(st[1][1,:])]
    return basis_ket
end

"""
returns the free energy computed using the exact eigenvalues for the exact evaluation of
the supercurrent over a finite set of eigenvectors
"""
function f_e(hpar::Quantica.ParametricHamiltonian, Δϕ,
    nev::Union{Int64, Missing}, Δ, dacp = false; kw...)
    method =(isa(nev, Missing) ? :adaptive_nev : :nev)
    if dacp == false
        λ = negative_eigen(hpar, Δϕ, nev, Δ, method = method; kw...)
    else
        λ = negative_eigen_dacp(hpar, Δϕ, gap; maxdeg = 2, kw...)
    end
    f = -sum(λ) #this minus comes from the tanh(1/β) when T->0
    return f
end

"""
returns the free energy computed using the exact eigenvalues and the andreev eigenstates
associated
"""
function f_e_st(hpar::Quantica.ParametricHamiltonian, Δϕ,
    nev::Union{Int64, Missing}, Δ)
    method =(isa(nev, Missing) ? :adaptive_nev : :nev)
    λ, st = negative_eigen(hpar, Δϕ, nev, Δ, method = method)
    f = -sum(λ) #this minus comes from the tanh(1/β) when T->0
    return f, st
end

function cut_st(hpar::Quantica.ParametricHamiltonian, Δϕ,
    nev::Union{Int64, Missing}, Δ)
    method =(isa(nev, Missing) ? :adaptive_nev : :nev)
    λ, st = eigen_st(hpar, Δϕ, nev, Δ/2, method = method,
    negative = true)
    # return convert(Matrix{ComplexF64}, qr(st).Q)
    # QR factorization to avoid degeneracies
    return st
end
"""
computes the eigenvalues and eigenvectors withing the gap using an adaptive procedure for
setting the number of nevs calculated with the shift invert method
"""
function eigen_st(hpar, Δϕ, nev, gap; method = :nev, nev0 = 4,
        negative = false)
    sp = spectrum(hpar(phi = Δϕ), method =
         ArpackPackage(nev = ifelse(method == :nev, nev, nev0),
        sigma = -0.001im))
    λ = real(sp.energies)
    st = sp.states
    if method == :adaptive_nev
        if maximum(abs.(λ)) < gap
            nev0 *= 2
            println(nev0)
            if nev0 < 401
                λ, st = eigen_st(hpar, Δϕ, missing, nev0 = nev0,
                    gap, method = :adaptive_nev)
            else @warn "we are missing subgap states!" end
        else nothing end
    else nothing end
    return ifelse(method == :nev, λ, λ[λ .< gap]),
        ifelse(negative == false, st, st[:,findall(x -> x < 0, λ)])
end

"""
returns the NEGATIVE eigenvalues withing the gap using for the exact supercurrent
calculation an adaptive procedure for setting the number of nevs calculated with the shift
invert method
"""
function negative_eigen(hpar, Δϕ, nev, gap; method = :nev, exact = false, nev0 = 5)
    if exact == true
        sp = spectrum(hpar(phi = Δϕ))
    else
        sp = spectrum(hpar(phi = Δϕ), method =
            ArpackPackage(nev = ifelse(method == :nev, nev, nev0),
            sigma = -0.001im))
    end
    λ = real(sp.energies)
    λneg = λ[λ.<=0];
    if method == :adaptive_nev
        if maximum(abs.(λneg)) < gap
            nev0 *= 2
            println(nev0)
            if nev0 < 401
                λneg = negative_eigen(hpar, Δϕ, missing,
                    nev0 = nev0, gap, method =:adaptive_nev)
            else @warn "we are missing subgap states!" end
        else nothing end
    else nothing end
    return ifelse(method == :nev, λneg, λneg[λneg .< gap])
end

"""
    `negative_eigen_dacp(hpar, Δϕ, gap; maxdeg = 2, kw...)`
returns the negative eigenvalues inside the interval (-a, a) using the DACP method,
which is suited to deal with large matrices
"""
function negative_eigen_dacp(hpar, Δϕ, gap; maxdeg = 2, kw...)
    a = gap * 2 # we select both gap and continuum spectrum
    eigs_dacp = DACPdiagonaliser(hpar(phi = Δϕ), a, numkets = maxdeg; kw...)
    λ = real(eigs_dacp)
    return λ[λ.<=0];
end

"""
current_cut. I do not set very precisely the boundaries but it doesn't matter since non
interface sites will have zero matrix elements H_ab
"""
edge_positions(h, myregion) =
    [pos for pos in sitepositions(h, region = myregion)]
edge_indices_l(h, scale, a0) =
    [ind for ind in siteindices(h, region = (r -> - 3 * scale * a0/2 < r[1] < 0))]
edge_indices_r(h, scale, a0) =
    [ind for ind in siteindices(h, region = (r -> 0 <= r[1] < 9*scale*a0/√3))]

"""
computes the supercurrent between two sites i and j
"""
function supercurrent_cut(p, states, h)
    eigenstates = subspace_kets(states, h)
    dn = h.harmonics[1].dn
    lat_positions = Quantica.allsitepositions(h.lattice)
    I = 0
    for i in 1:length(eigenstates)
        I += supercurrent_e(h, eigenstates[i], dn, lat_positions, p.scale, p.a0)
    end
    return I
end

"""
computes the supercurrent for a given energy state, psi, along a vertical (x = 0) cut
(this could have been made more efficiently selecting only those sites conected by hoppings
with a hopping selector and a region - future optimization)
"""
function supercurrent_e(h, psi, dn, lat_positions, scale, a0)
    I_e = 0
    for i in edge_indices_l(h, scale, a0)
        for j in edge_indices_r(h, scale, a0)
           I_e += supercurrent_ij(h, psi, i, j, dn, lat_positions)
        end
    end
    return I_e
end

"computes the supercurrent due to a given state psi across the i - j link"
supercurrent_ij(h, psi, i, j, dn, lat_positions) =
    im*(psi[j]' * h[dn][i, j]' * σ0τz  * psi[i] -
        psi[i]' * σ0τz * h[dn][i, j] * psi[j])


text:
function IMPLEMENT_ONE_ARG_FUNC(meth, symnm; lib=:basic_)
     @eval begin
        function ($meth)(b::SymbolicType)
            a = Basic()
            ccall(($(string(lib,symnm)), libsymengine), Void, (Ptr{Basic}, Ptr{Basic}), &a, &b)
            return a
        end
    end
end

function IMPLEMENT_TWO_ARG_FUNC(meth, symnm; lib=:basic_)
     @eval begin
        function ($meth)(b1::SymbolicType, b2::Number)
            a = Basic()
            b1, b2 = promote(b1, b2)
            ccall(($(string(lib,symnm)), libsymengine), Void, (Ptr{Basic}, Ptr{Basic}, Ptr{Basic}), &a, &b1, &b2)
            return a
        end
    end
end

## import from base one argument functions
## these are from cwrapper.cpp, one arg func
## Where are exp? log?, sqrt?
for (meth, libnm) in [
                      (:abs,:abs),
                      (:sin,:sin),
                      (:cos,:cos),
                      (:tan,:tan),
                      (:csc,:csc),
                      (:sec,:sec),
                      (:cot,:cot),
                      (:asin,:asin),
                      (:acos,:acos),
                      (:asec,:asec),
                      (:acsc,:acsc),
                      (:atan,:atan),
                      (:acot,:acot),
                      (:sinh,:sinh),
                      (:cosh,:cosh),
                      (:tanh,:tanh),
                      (:csch,:csch),
                      (:sech,:sech),
                      (:coth,:coth),
                      (:asinh,:asinh),
                      (:acosh,:acosh),
                      (:asech,:asech),
                      (:acsch,:acsch),
                      (:atanh,:atanh),
                      (:acoth,:acoth),
                      (:zeta,:zeta),
                      (:gamma,:gamma),
                      (:eta,:dirichlet_eta),
                      ]
    eval(Expr(:import, :Base, meth))
    IMPLEMENT_ONE_ARG_FUNC(meth, libnm)
end

# export not import
for  (meth, libnm) in [
                       (:lambertw,:lambertw)   # in add-on packages, not base
                       ]
    IMPLEMENT_ONE_ARG_FUNC(meth, libnm)
    eval(Expr(:export, meth))
end

## add these in until they are wrapped
Base.exp(a::SymbolicType) = E^a
function Base.log(a::SymbolicType) # super hacky
    u = symbols(gensym())
    v = a^u
    diff(v, u) / v
end
Base.sqrt(a::SymbolicType) = a^(1//2)
Base.cbrt(a::SymbolicType) = a^(1//3)

for (meth, fn) in [(:sind, :sin), (:cosd, :cos), (:tand, :tan), (:secd, :sec), (:cscd, :csc), (:cotd, :cot)]
    eval(Expr(:import, :Base, meth))
    @eval begin
        $(meth)(a::SymbolicType) = $(fn)(a*PI/180)
    end
end

## Number theory module from cppwrapper
for (meth, libnm) in [(:gcd, :gcd),
                      (:lcm, :lcm),
                      (:mod, :mod),
                      (:div, :quotient),
                      ]
    eval(Expr(:import, :Base, meth))
    IMPLEMENT_TWO_ARG_FUNC(meth, libnm, lib=:ntheory_)
end

Base.binomial(n::Basic, k::Number) = binomial(N(n), N(k))  #ntheory_binomial seems wrong
Base.rem(a::SymbolicType, b::SymbolicType) = a - (a ÷ b) * b
Base.factorial(n::SymbolicType, k) = factorial(N(n), N(k))

## but not (:fibonacci,:fibonacci), (:lucas, :lucas) (Basic type is not the signature)
for (meth, libnm) in [(:nextprime,:nextprime)
                      ]
    IMPLEMENT_ONE_ARG_FUNC(meth, libnm, lib=:ntheory_)
    eval(Expr(:export, meth))
end


text:
using BinaryBuilder

name = "SuiteSparse"
version = v"5.4.0"

# Collection of sources required to build SuiteSparse
sources = [
    ArchiveSource("https://github.com/DrTimothyAldenDavis/SuiteSparse/archive/v$(version).tar.gz",
                  "d9d62d539410d66550d0b795503a556830831f50087723cb191a030525eda770"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/SuiteSparse-*

# Apply Jameson's shlib patch
atomic_patch -p1 ${WORKSPACE}/srcdir/patches/SuiteSparse-shlib.patch

# Disable OpenMP as it will probably interfere with blas threads and Julia threads
FLAGS+=(INSTALL="${prefix}" INSTALL_LIB="${libdir}" INSTALL_INCLUDE="${prefix}/include" CFOPENMP=)

if [[ ${target} == *mingw32* ]]; then
    FLAGS+=(UNAME=Windows)
    FLAGS+=(LDFLAGS="${LDFLAGS} -L${libdir} -shared")
else
    FLAGS+=(UNAME="$(uname)")
    FLAGS+=(LDFLAGS="${LDFLAGS} -L${libdir}")
fi

if [[ ${nbits} == 64 ]] && [[ ${target} != aarch64* ]]; then
    SUN="-DSUN64 -DLONGBLAS='long long'"

    FLAGS+=(BLAS="-lopenblas64_" LAPACK="-lopenblas64_")
else
    FLAGS+=(BLAS="-lopenblas" LAPACK="-lopenblas")
fi

# Disable METIS in CHOLMOD by passing -DNPARTITION and avoiding linking metis
#FLAGS+=(MY_METIS_LIB="-lmetis" MY_METIS_INC="${prefix}/include")
FLAGS+=(UMFPACK_CONFIG="$SUN" CHOLMOD_CONFIG+="$SUN -DNPARTITION" SPQR_CONFIG="$SUN")

make -j${nproc} -C SuiteSparse_config "${FLAGS[@]}" library config

for proj in SuiteSparse_config AMD BTF CAMD CCOLAMD COLAMD CHOLMOD LDL KLU UMFPACK RBio SPQR; do
    make -j${nproc} -C $proj "${FLAGS[@]}" library CFOPENMP="$CFOPENMP"
    make -j${nproc} -C $proj "${FLAGS[@]}" install CFOPENMP="$CFOPENMP"
done

# For now, we'll have to adjust the name of the OpenBLAS library on macOS and FreeBSD.
# Eventually, this should be fixed upstream
if [[ ${target} == *-apple-* ]] || [[ ${target} == *freebsd* ]]; then
    echo "-- Modifying library name for OpenBLAS"

    for nm in libcholmod libspqr libumfpack; do
        # Figure out what version it probably latched on to:
        if [[ ${target} == *-apple-* ]]; then
            OPENBLAS_LINK=$(otool -L ${libdir}/${nm}.dylib | grep libopenblas64_ | awk '{ print $1 }')
            install_name_tool -change ${OPENBLAS_LINK} @rpath/libopenblas64_.dylib ${libdir}/${nm}.dylib
        elif [[ ${target} == *freebsd* ]]; then
            OPENBLAS_LINK=$(readelf -d ${libdir}/${nm}.so | grep libopenblas64_ | sed -e 's/.*\[\(.*\)\].*/\1/')
            patchelf --replace-needed ${OPENBLAS_LINK} libopenblas64_.so ${libdir}/${nm}.so
        fi
    done
fi

# Delete the extra soversion libraries built. https://github.com/JuliaPackaging/Yggdrasil/issues/7
if [[ "${target}" == *-mingw* ]]; then
    rm -f ${libdir}/lib*.*.${dlext}
    rm -f ${libdir}/lib*.*.*.${dlext}
fi

install_license LICENSE.txt

# Compile SuiteSparse_wrapper shim
cd $WORKSPACE/srcdir/SuiteSparse_wrapper
"${CC}" -O2 -shared -fPIC -I${prefix}/include SuiteSparse_wrapper.c -o ${libdir}/libsuitesparse_wrapper.${dlext} -L${libdir} -lcholmod
"""

# 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 = [
    LibraryProduct("libsuitesparseconfig",   :libsuitesparseconfig),
    LibraryProduct("libamd",                 :libamd),
    LibraryProduct("libbtf",                 :libbtf),
    LibraryProduct("libcamd",                :libcamd),
    LibraryProduct("libccolamd",             :libccolamd),
    LibraryProduct("libcolamd",              :libcolamd),
    LibraryProduct("libcholmod",             :libcholmod),
    LibraryProduct("libldl",                 :libldl),
    LibraryProduct("libklu",                 :libklu),
    LibraryProduct("libumfpack",             :libumfpack),
    LibraryProduct("librbio",                :librbio),
    LibraryProduct("libspqr",                :libspqr),
    LibraryProduct("libsuitesparse_wrapper", :libsuitesparse_wrapper),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("OpenBLAS_jll"),
#    Dependency("METIS_jll"),
]

# Build the tarballs
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
#push!(LOAD_PATH,"Simulation")
#import Simulation
include("Simulation/Simulation.jl")

using Random
using Profile
using FileIO
using JLD2

using Distributions
using DataFrames

let
  individuals_df = load("../../data/simulations/wroclaw-population-orig.jld2")["individuals_df"]

  if(!(:social_competence in names(individuals_df)))
    #Produce some test data for friendship kernel
    individuals_df.social_competence = [Float32(x) for x in rand(Beta(3.0, 3.0), nrow(individuals_df))]
  end

  #loads and pre-processes data that is constant during the simulation time

  mild_detection_prob = 0.1
  constant_kernel_param = 1.0
  household_kernel_param = 0.3
  tracking_prob = 0.3
  tracking_delay = 5

  param_rng = MersenneTwister(1)

  params = Simulation.load_params(
    param_rng,
    population = individuals_df,

    mild_detection_prob = mild_detection_prob,

    constant_kernel_param = constant_kernel_param,
    household_kernel_param = household_kernel_param,

    backward_tracking_prob = tracking_prob,
    backward_detection_delay = tracking_delay/2,

    forward_tracking_prob = tracking_prob,
    forward_detection_delay = tracking_delay/2,

    testing_time = tracking_delay/2
  )

  state = Simulation.SimState(Simulation.num_individuals(params), seed=123)
  sizehint!(state.queue, 10^6)

  push!(state.queue, Simulation.Event(Val(Simulation.OutsideInfectionEvent),0.0,1))

  function run_some(num_iters::Integer)
    for iter in 1:num_iters
        if isempty(state.queue)
            @info "queue empty before real testing starts, after $iter events"
            break
        end
        event = pop!(state.queue)
        result = Simulation.execute!(state, params, event)
    end
  end

  #@code_warntype run_some(10000)

  run_some(10000)

  Profile.clear_malloc_data()

  run_some(10000)
end


text:
"""
    Restriction(base_freq::Vector{Float64},
                SubstitutionRates::Vector{Float64})::Tuple{Array{Float64,2}, Array{Float64,1}, Array{Float64,2}}

Calculate the eigenvalue decomposition of the Q matrix of the restriction site model.
The `SubstitutionRates` are ignored, and just for call stability.

The function returns the Eigenvectors, Eigenvalues, inverse of eigenvectors and
    the scale factor for expected number changes per site
"""
function Restriction(base_freq::Vector{Float64}, SubstitutionRates::Vector{Float64})::Tuple{Array{Float64,2}, Array{Float64,1}, Array{Float64,2}, Float64}
    Nbases = length(base_freq)
    Q::Array{Float64,2} = ones(Nbases,Nbases)
    Q[diagind(Nbases,Nbases)] .= -1
    Q .*= reverse(base_freq)
    D, U = eigen(Q)
    Uinv = inv(U)
    mu::Float64 =  1.0 / (2.0 * prod(base_freq))
    return U, D, Uinv, mu
end

"""
    JC(base_freq::Vector{Float64},
        SubstitutionRates::Vector{Float64})::Tuple{Array{Float64,2}, Array{Float64,1}, Array{Float64,2}}

Calculate the eigenvalue decomposition of the Q matrix of the Jukes-Cantor model.
The `SubstitutionRates` are ignored, and just for call stability.

The function returns the Eigenvectors, Eigenvalues and inverse of eigenvectors.
"""
function JC(base_freq::Vector{Float64}, SubstitutionRates::Vector{Float64})::Tuple{Array{Float64,2}, Array{Float64,1}, Array{Float64,2}, Float64}
    Nbases = length(base_freq)
    Q::Array{Float64,2} = ones(Nbases,Nbases)
    μ = SubstitutionRates[1]
    off_diag = μ/Nbases
    diag = off_diag * (Nbases-1)
    Q .= off_diag
    Q[diagind(Nbases,Nbases)] .= -diag
    D, U = eigen(Q)
    Uinv = inv(U)
    mu = 1.0
    return U, D, Uinv, mu
end

"""
    GTR(base_freq::Vector{Float64},
        SubstitutionRates::Vector{Float64})::Tuple{Array{Float64,2}, Array{Float64,1}, Array{Float64,2}}

Calculate the eigenvalue decomposition of the Q matrix of the General Time Reversible model.

The function returns the Eigenvectors, Eigenvalues and inverse of eigenvectors.
"""
function GTR(base_freq::Vector{Float64}, SubstitutionRates::Vector{Float64})::Tuple{Array{Float64,2}, Array{Float64,1}, Array{Float64,2}, Float64}
    Nbases = length(base_freq)
    Q::Array{Float64,2} = setmatrix(SubstitutionRates)
    Q = mapslices(x-> x.*base_freq, Q, dims=2)
    dia = sum(Q, dims=2)
    Q[diagind(Nbases,Nbases)] = -dia
    D, U = eigen(Q)
    Uinv = inv(U)
    return U, D, Uinv, 1/sum(dia)
end

"""
    function freeK(base_freq::Vector{Float64},SubstitutionRates::AbstractArray)::Tuple{Array{N,2}, Array{N,1}, Array{N,2}, M} where {N <: Number, M <: Number}

FreeK model of substitution.

"""
function freeK(
      base_freq::Vector{Float64},
      SubstitutionRates::AbstractArray,
      )::Tuple{Array, Array, Array, Float64}# where {N <: Number, M <: Number}
      Nrates = length(SubstitutionRates)
      Nbases = Int(ceil(sqrt(Nrates)))
      Q = zeros(Nbases, Nbases)
      counter = 1
      for i in 1:Nbases, j in 1:Nbases
            if i!=j
                  Q[j,i] = SubstitutionRates[counter]
                  counter += 1
            end
      end
      dia = sum(Q, dims = 2)
      Q[diagind(Nbases, Nbases)] = -dia
      D, U = eigen(Q)
      Uinv = inv(U)
      return U, D, Uinv, 1 / sum(dia)
end

## Helper Functions ##
function setmatrix(vec_vals::Array{T,1})::Array{T,2} where T
    n = length(vec_vals)
    s = round((sqrt(8n+1)+1)/2)
    s*(s-1)/2 == n || error("setmatrix: length of vector is not triangular")
    k = 0
    Q = [ i<j ? (k+=1; vec_vals[k]) : 0 for i=1:s, j=1:s ]
    Q.+transpose(Q)
end


text:
module StabilityCriteria

using LinearAlgebra: diag

using Crystallography:
    CrystalSystem, Cubic, Hexagonal, Tetragonal, Trigonal, Orthorhombic, Monoclinic

using LinearElasticity: StiffnessMatrix, ComplianceMatrix

export isstable

ispositive(x) = x > zero(x)

function stability_criteria(::Cubic, c::StiffnessMatrix)
    c₁₁, c₁₂, c₄₄ = c[1, 1], c[1, 2], c[4, 4]
    return (  # Must satisfy all criteria!
        c₁₁ > abs(c₁₂),
        c₁₁ > -2c₁₂,
        ispositive(c₄₄),
    )
end
function stability_criteria(::Hexagonal, c::StiffnessMatrix)
    c₁₁, c₁₂, c₁₃, c₃₃, c₄₄, c₆₆ = c[1, 1], c[1, 2], c[1, 3], c[3, 3], c[4, 4], c[6, 6]
    return (  # Must satisfy all criteria!
        2c₆₆ == c₁₁ - c₁₂,
        c₁₁ > abs(c₁₂),
        2c₁₃^2 < c₃₃ * (c₁₁ + c₁₂),
        ispositive(c₄₄),
        ispositive(c₆₆),
    )
end
function stability_criteria(::Tetragonal, c::StiffnessMatrix)
    c₁₁, c₁₂, c₁₃, c₁₆, c₃₃, c₄₄, c₆₆ =
        c[1, 1], c[1, 2], c[1, 3], c[1, 6], c[3, 3], c[4, 4], c[6, 6]
    if iszero(c₁₆)  # Tetragonal (I) class
        return (
            c₁₁ > abs(c₁₂),
            2c₁₃^2 < c₃₃ * (c₁₁ + c₁₂),
            ispositive(c₄₄),
            ispositive(c₆₆),
        )
    else
        return (  # Tetragonal (II) class
            c₁₁ > abs(c₁₂),
            2c₁₃^2 < c₃₃ * (c₁₁ + c₁₂),
            ispositive(c₄₄),
            2c₁₆^2 < c₆₆ * (c₁₁ - c₁₂),
        )
    end
end
function stability_criteria(::Trigonal, c::StiffnessMatrix)
    c₁₁, c₁₂, c₁₃, c₁₄, c₁₅, c₃₃, c₄₄, c₆₆ =
        c[1, 1], c[1, 2], c[1, 3], c[1, 4], c[1, 5], c[3, 3], c[4, 4], c[6, 6]
    return (
        c₁₁ > abs(c₁₂),
        ispositive(c₄₄),
        2c₁₃^2 < c₃₃ * (c₁₁ + c₁₂),
        c₄₄ * (c₁₁ - c₁₂) == 2c₄₄ * c₆₆,
        if iszero(c₁₅)  # Rhombohedral (I) class
            2c₁₄^2 < c₄₄ * (c₁₁ - c₁₂)
        else  # Rhombohedral (II) class
            2(c₁₄^2 + c₁₅^2) < c₄₄ * (c₁₁ - c₁₂)
        end,
    )
end
function stability_criteria(::Orthorhombic, c::StiffnessMatrix)
    c₁₁, c₂₂, c₃₃, c₄₄, c₅₅, c₆₆ = diag(c)
    c₁₂, c₁₃, c₂₃ = c[1, 2], c[1, 3], c[2, 3]
    return (
        map(ispositive, (c₁₁, c₄₄, c₅₅, c₆₆))...,
        c₁₁ * c₂₂ > c₁₂^2,
        c₁₁ * c₂₂ * c₃₃ + 2c₁₂ * c₁₃ * c₂₃ > c₁₁ * c₂₃^2 + c₂₂ * c₁₃^2 + c₃₃ * c₁₂^2,
    )
end
function stability_criteria(::Monoclinic, c::StiffnessMatrix)
    c₁₁, c₂₂, c₃₃, c₄₄, c₅₅, c₆₆ = diag(c)
    c₁₂, c₁₃, c₁₅, c₂₃, c₂₅, c₃₅, c₄₆ =
        c[1, 2], c[1, 3], c[1, 5], c[2, 3], c[2, 5], c[3, 5], c[4, 6]
    g =
        c₁₁ * c₂₂ * c₃₃ - c₁₁ * c₂₃ * c₂₃ - c₂₂ * c₁₃ * c₁₃ - c₃₃ * c₁₂ * c₁₂ +
        2c₁₂ * c₁₃ * c₂₃
    return (
        map(ispositive, (c₁₁, c₂₂, c₃₃, c₄₄, c₅₅, c₆₆))...,
        c₁₁ + c₂₂ + c₃₃ > -2(c₁₂ + c₁₃ + c₂₃),
        c₃₃ * c₅₅ > c₃₅^2,
        c₄₄ * c₆₆ > c₄₆^2,
        c₂₂ + c₃₃ > 2c₂₃,
        c₂₂ * (c₃₃ * c₅₅ - c₃₅^2) + 2c₂₃ * c₂₅ * c₃₅ - c₂₃^2 * c₅₅ > c₂₅^2 * c₃₃,
        2(
            c₁₅ * c₂₅ * (c₃₃ * c₁₂ - c₁₃ * c₂₃) +
            c₁₅ * c₃₅ * (c₂₂ * c₁₃ - c₁₂ * c₂₃) +
            c₂₅ * c₃₅ * (c₁₁ * c₂₃ - c₁₂ * c₁₃)
        ) - (
            c₁₅ * c₁₅ * (c₂₂ * c₃₃ - c₂₃^2) +
            c₂₅ * c₂₅ * (c₁₁ * c₃₃ - c₁₃^2) +
            c₃₅ * c₃₅ * (c₁₁ * c₂₂ - c₁₂^2)
        ) > -c₅₅ * g,
    )
end
stability_criteria(C::CrystalSystem, s::ComplianceMatrix) = stability_criteria(C, inv(s))

isstable(C::CrystalSystem, x::Union{StiffnessMatrix,ComplianceMatrix}) =
    all(stability_criteria(C, x))

end


text:
module Curl

export
    with_handle,
    Easy,
        set_url,
        set_method,
        set_verbose,
        set_debug,
        set_body,
        set_upload_size,
        set_seeker,
        set_ca_roots_path,
        set_timeout,
        add_headers,
        enable_upload,
        enable_progress,
        upload_data,
        get_protocol,
        get_effective_url,
        get_response_status,
        get_response_info,
        get_curl_errstr,
        status_ok,
    Multi,
        add_handle,
        remove_handle

using LibCURL
using LibCURL: curl_off_t, libcurl
# not exported: https://github.com/JuliaWeb/LibCURL.jl/issues/87

# constants that LibCURL should have but doesn't
const CURLE_PEER_FAILED_VERIFICATION = 60
const CURLSSLOPT_REVOKE_BEST_EFFORT = 1 << 3

# these are incorrectly defined on Windows by LibCURL:
if Sys.iswindows()
    const curl_socket_t = Base.OS_HANDLE
    const CURL_SOCKET_TIMEOUT = Base.INVALID_OS_HANDLE
else
    const curl_socket_t = Cint
    const CURL_SOCKET_TIMEOUT = -1
end

# definitions affected by incorrect curl_socket_t (copied verbatim):
function curl_multi_socket_action(multi_handle, s, ev_bitmask, running_handles)
    ccall((:curl_multi_socket_action, libcurl), CURLMcode, (Ptr{CURLM}, curl_socket_t, Cint, Ptr{Cint}), multi_handle, s, ev_bitmask, running_handles)
end
function curl_multi_assign(multi_handle, sockfd, sockp)
    ccall((:curl_multi_assign, libcurl), CURLMcode, (Ptr{CURLM}, curl_socket_t, Ptr{Cvoid}), multi_handle, sockfd, sockp)
end

# additional curl_multi_socket_action method
function curl_multi_socket_action(multi_handle, s, ev_bitmask)
    curl_multi_socket_action(multi_handle, s, ev_bitmask, Ref{Cint}())
end

using FileWatching
using NetworkOptions
using Base: OS_HANDLE, preserve_handle, unpreserve_handle

include("utils.jl")

function __init__()
    @check curl_global_init(CURL_GLOBAL_ALL)
end

const CURL_VERSION_INFO = unsafe_load(curl_version_info(CURLVERSION_NOW))
const SSL_VERSION = unsafe_string(CURL_VERSION_INFO.ssl_version)
const SYSTEM_SSL =
    Sys.isapple() && startswith(SSL_VERSION, "SecureTransport") ||
    Sys.iswindows() && startswith(SSL_VERSION, "Schannel")

const CURL_VERSION_STR = unsafe_string(curl_version())
let m = match(r"^libcurl/(\d+\.\d+\.\d+)\b", CURL_VERSION_STR)
    m !== nothing || error("unexpected CURL_VERSION_STR value")
    curl = m.captures[1]
    julia = "$(VERSION.major).$(VERSION.minor)"
    global const CURL_VERSION = VersionNumber(curl)
    global const USER_AGENT = "curl/$curl julia/$julia"
end

include("Easy.jl")
include("Multi.jl")

function with_handle(f, handle::Union{Multi, Easy})
    try f(handle)
    finally
        Curl.done!(handle)
    end
end

setopt(easy::Easy, option::Integer, value) =
    @check curl_easy_setopt(easy.handle, option, value)
setopt(multi::Multi, option::Integer, value) =
    @check curl_multi_setopt(multi.handle, option, value)

end # module

max_stars_repo_name: StefanKarpinski/Downloader

max_issues_repo_name: StefanKarpinski/Downloader


text:
# Convert # of Rata Die days to proleptic Gregorian calendar y,m,d,w
# Reference: http://mysite.verizon.net/aesir_research/date/date0.htm

# Accessor functions
value(sp::SimpleDate) = sp.instant.periods.value
value(p::Period) = p.value
function days(q::QuarterDate)
    y1, q1 = divrem(value(q), 4)
    m1 = 3 * q1
    return value(DayDate(y1, m1, daysinmonth(m1)))
end
function days(m::MonthDate)
    y1, m1 = divrem(value(m), 12)
    return value(DayDate(y1, m1, daysinmonth(m1)))
end
days(w::WeekDate) = 7 * value(w)
days(d::DayDate) = value(d)
weeks(w::WeekDate) = value(w)
months(m::MonthDate) = value(m)
quarters(q::QuarterDate) = value(q)
semesters(s::SemesterDate) = value(s)
year(y::YearDate) = value(y) + 1
year(s::SemesterDate) = div(value(s) - 1, 2) + 1
year(q::QuarterDate) = div(value(q) - 1, 4) + 1
year(m::MonthDate) = div(value(m) - 1, 12) + 1
year(w::WeekDate) = Dates.year((value(w) * 7))
year(d::DayDate) = Dates.year(value(d))
semester(s::SemesterDate) = rem(abs(value(s)) - 1, 2) + 1
semester(q::QuarterDate) = rem(abs(value(q)) - 1, 2) + 1
semester(m::MonthDate) = div(rem(abs(value(m)) - 1, 12), 2) + 1
semester(d::DayDate) = semester(value(d))
quarter(q::QuarterDate) = rem(abs(value(q)) - 1, 4) + 1
quarter(m::MonthDate) = div(rem(abs(value(m)) - 1, 12), 4) + 1
quarter(d::DayDate) = quarter(value(d))
month(m::MonthDate) = rem(abs(value(m)) - 1, 12) + 1
month(d::DayDate) = month(value(d))
week(w::WeekDate) = Dates.week(value(w) * 7)
week(d::DayDate) = Dates.week(value(d))
day(d::DayDate) = Dates.day(value(d))
undated(u::UndatedDate) = value(u)

# accessor functions for number of days
year(d::Int64) = Dates.year(d)

function semester(d::Int64)
    y, m = Dates.yearmonth(d)
    return (m < 7) ? 1 : 2
end

function quarter(d::Int64)
    y, m = Dates.yearmonth(d)
    return div(m - 1, 3) + 1
end

function month(d::Int64)
    y, m = Dates.yearmonth(d)
    return m
end

week(d) = Dates.week(d)
day(d) = Dates.day(d)

dayofmonth(dt::DayDate) = day(dt)

yearmonth(dt::DayDate) = Dates.yearmonth(days(dt))
monthday(dt::DayDate) = Dates.monthday(days(dt))
yearmonthday(dt::DayDate) = Dates.yearmonthday(days(dt))

# Documentation for exported accessors
for func in (:year, :month, :quarter)
    name = string(func)
    @eval begin
        @doc """
            $($name)(dt::Dates.TimeType) -> Int64

        The $($name) of a `Date` or `DateTime` as an [`Int64`](@ref).
        """ $func(dt::Dates.TimeType)
    end
end

"""
    week(dt::Dates.TimeType) -> Int64

Return the [ISO week date](https://en.wikipedia.org/wiki/ISO_week_date) of a `Date` or
`DateTime` as an [`Int64`](@ref). Note that the first week of a year is the week that
contains the first Thursday of the year, which can result in dates prior to January 4th
being in the last week of the previous year. For example, `week(Date(2005, 1, 1))` is the 53rd
week of 2004.

# Examples
```jldoctest
julia> Dates.week(Date(1989, 6, 22))
25

julia> Dates.week(Date(2005, 1, 1))
53

julia> Dates.week(Date(2004, 12, 31))
53
```
"""
week(dt::Dates.TimeType)

for func in (:day, :dayofmonth)
    name = string(func)
    @eval begin
        @doc """
            $($name)(dt::Dates.TimeType) -> Int64

        The day of month of a `Date` or `DateTime` as an [`Int64`](@ref).
        """ $func(dt::Dates.TimeType)
    end
end


text:
struct UnindexedDimsArray{T, N, IA<:IndexedDimsArray{T, N}} <: AbstractArray{T, N}
    data::IA
end

Base.parent(arr::UnindexedDimsArray) = arr.data
Base.size(arr::UnindexedDimsArray) = size(parent(arr))

function Base.getindex(arr::UnindexedDimsArray, inds::Int...)
    return getindex(parent(parent(arr)), inds...)
end

function Base.setindex!(arr::UnindexedDimsArray, val, inds::Int...)
    return setindex!(parent(parent(arr)), val, inds...)
end

# function Base.similar(arr::UnindexedDimsArray, )


text:
# # Foundations of Bilevel Programming: Example 2.1
# This example is from the book _Foundations of Bilevel Programming_ by Stephan
# Dempe, Chapter 3.2, Page 25. [url](https://www.springer.com/gp/book/9781402006319)
# Moving the bound on x to lower level

# Model of the problem

# First level
# ```
# Min 3x + y
# s.t.
# 0 <= y <=8
# ```
# Second level
# ```
# Min -x
# s.t.
# x + y <= 8
# 4x + y >= 8
# 2x + y <= 13
# 2x - y <= 0
# x <= 5
# ```

using BilevelJuMP
using Ipopt
using JuMP
using Test

model = BilevelModel(Ipopt.Optimizer, mode = BilevelJuMP.ProductMode(1e-9))

# First we need to create all of the variables in the upper and lower problems:

@variable(Lower(model), x, start = 3.5 * 8 / 15)

#-

@variable(Upper(model), y, start = 8 / 15)

# Then we can add the objective and constraints of the upper problem:

@objective(Upper(model), Min, 3x + y)
#-

@constraint(Upper(model), y <= 8)
@constraint(Upper(model), y >= 0)

# Followed by the objective and constraints of the lower problem:

@objective(Lower(model), Min, -x)

#-
@constraint(Lower(model), x +  y <= 8)
@constraint(Lower(model), 4x +  y >= 8)
@constraint(Lower(model), 2x +  y <= 13)
@constraint(Lower(model), 2x - 7y <= 0)
@constraint(Lower(model), x <= 5)

# Now we can solve the problem and verify the solution again that reported by
# Dempe.

optimize!(model)

@test objective_value(model) ≈ 3 * (3.5 * 8 / 15) + (8 / 15) atol=1e-3
@test value(x) ≈ 3.5 * 8 / 15 atol=1e-6
@test value(y) ≈ 8 / 15 atol=1e-6

# TODO: why are these commented out?    #src
# @test dual(l2) #≈ [0] atol=atol       #src
# @test dual(l4) #≈ [0] atol=atol       #src


text:
module DesignDriveCoilFiltGUI
using FFTW,Gtk,PyPlot,Optim
using LinearAlgebra
using PyPlot
using CSV
using Reexport

@reexport using DesignDriveCoilFilt

# using Interpolations

include("Filter_Designer_GUI.jl")
export useFilterDesignGUI

end


text:
using DataDeps

register(DataDep(
  "poker-hand",
  "http://archive.ics.uci.edu/ml/datasets/Poker+Hand",
  "http://archive.ics.uci.edu/ml/machine-learning-databases/poker/poker-hand-training-true.data",
  "37becdf87d5f8cbf2b91d6471e965a25b86cb4a6d878c0f94a4025969fca464f",
  post_fetch_method=(path -> begin
    UCIData.process_dataset(path,
        target_index=11,
        feature_indices=1:10,
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
import Base: map, AbstractVector, AbstractMatrix
import Distances: pairwise
export BlockMean, BlockKernel, BlockCrossKernel

"""
    BlockMean <: MeanFunction

(Vertically) concatentated mean functions.
"""
struct BlockMean <: MeanFunction
    μ::Vector
end
BlockMean(μs::Vararg{<:MeanFunction}) = BlockMean([μs...])
length(μ::BlockMean) = sum(length.(μ.μ))
==(μ::BlockMean, μ′::BlockMean) = μ.μ == μ′.μ
@noinline eachindex(μ::BlockMean) = BlockData(eachindex.(μ.μ))
map(μ::BlockMean, X::BlockData) = BlockVector(map.(μ.μ, blocks(X)))

# Define the zero element.
zero(μ::BlockMean) = BlockMean(zero.(μ.μ))

"""
    BlockCrossKernel <: CrossKernel

A cross kernel comprising lots of other kernels.
"""
struct BlockCrossKernel <: CrossKernel
    ks::Matrix
end
BlockCrossKernel(ks::AbstractVector) = BlockCrossKernel(reshape(ks, length(ks), 1))
function BlockCrossKernel(ks::Adjoint{T, AbstractVector{T}} where T)
    return BlockCrossKernel(reshape(ks, 1, length(ks)))
end
size(k::BlockCrossKernel, N::Int) = N == 1 ?
    sum(size.(k.ks[:, 1], Ref(1))) :
    N == 2 ? sum(size.(k.ks[1, :], Ref(2))) : 1
==(k::BlockCrossKernel, k′::BlockCrossKernel) = k.ks == k′.ks
function eachindex(k::BlockCrossKernel, N::Int)
    if N == 1
        return BlockData(eachindex.(k.ks[:, 1], 1))
    elseif N == 2
        return BlockData(eachindex.(k.ks[1, :], 2))
    else
        throw(error("N ∉ {1, 2}"))
    end
end
map(k::BlockCrossKernel, X::BlockData) = BlockVector(map.(diag(k.ks), blocks(X)))
function map(k::BlockCrossKernel, X::BlockData, X′::BlockData)
    return BlockVector(map.(diag(k.ks), blocks(X), blocks(X′)))
end
function pairwise(k::BlockCrossKernel, X::BlockData, X′::BlockData)
    return BlockMatrix(broadcast(pairwise, k.ks, blocks(X), reshape(blocks(X′), 1, :)))
end
zero(k::BlockCrossKernel) = BlockCrossKernel(zero.(k.ks))

"""
    BlockKernel <: Kernel

A kernel comprising lots of other kernels. This is represented as a matrix whose diagonal
elements are `Kernels`, and whose off-diagonal elements are `CrossKernel`s. In the absence
of determining at either either compile- or construction-time whether or not this actually
constitutes a valid Mercer kernel, we take the construction of this type to be a promise on
the part of the caller that the thing they are constructing does indeed constitute a valid
Mercer kernel.

`ks_diag` represents the kernels on the diagonal of this matrix-valued kernel, and `ks_off`
represents the elements in the rest of the matrix. Only the upper-triangle will actually
be used.
"""
struct BlockKernel <: Kernel
    ks_diag::Vector{<:Kernel}
    ks_off::Matrix{<:CrossKernel}
end
length(k::BlockKernel) = sum(size.(k.ks_diag, 1))
eachindex(k::BlockKernel) = BlockData(eachindex.(k.ks_diag))
==(k::BlockKernel, k′::BlockKernel) = k.ks_diag == k′.ks_diag && k.ks_off == k′.ks_off

map(k::BlockKernel, X::BlockData) = BlockVector(map.(k.ks_diag, blocks(X)))
function map(k::BlockKernel, X::BlockData, X′::BlockData)
    return BlockVector(map.(k.ks_diag, blocks(X), blocks(X′)))
end

Base.adjoint(z::Zeros{T}) where {T} = Zeros{T}(reverse(size(z)))

function pairwise(k::BlockKernel, X::BlockData)
    bX = blocks(X)
    Σ = BlockArray(undef_blocks, AbstractMatrix{Float64}, length.(bX), length.(bX))
    for q in eachindex(k.ks_diag)
        setblock!(Σ, unbox(pairwise(k.ks_diag[q], bX[q])), q, q)
        for p in 1:q-1
            setblock!(Σ, pairwise(k.ks_off[p, q], bX[p], bX[q]), p, q)
            setblock!(Σ, getblock(Σ, p, q)', q, p)
        end
    end
    return LazyPDMat(Symmetric(Σ))
end
function pairwise(k::BlockKernel, X::BlockData, X′::BlockData)
    bX, bX′ = blocks(X), blocks(X′)
    Ω = BlockArray(undef_blocks, AbstractMatrix{Float64}, length.(bX), length.(bX′))
    for q in eachindex(k.ks_diag), p in eachindex(k.ks_diag)
        if p == q
            setblock!(Ω, pairwise(k.ks_diag[p], bX[p], bX′[p]), p, p)
        elseif p < q
            setblock!(Ω, pairwise(k.ks_off[p, q], bX[p], bX′[q]), p, q)
        else
            setblock!(Ω, pairwise(k.ks_off[q, p], bX[p], bX′[q]), p, q)
        end
    end
    return Ω
end
@noinline zero(k::BlockKernel) = BlockKernel(zero.(k.ks_diag), zero.(k.ks_off))


text:
function format_double(::Type{CompositeFormat}, arg0::jdouble, arg1::NumberFormat, arg2::StringBuffer, arg3::FieldPosition)
    return jcall(CompositeFormat, "formatDouble", StringBuffer, (jdouble, NumberFormat, StringBuffer, FieldPosition), arg0, arg1, arg2, arg3)
end

function get_default_number_format(::Type{CompositeFormat})
    return jcall(CompositeFormat, "getDefaultNumberFormat", NumberFormat, ())
end

function get_default_number_format(::Type{CompositeFormat}, arg0::Locale)
    return jcall(CompositeFormat, "getDefaultNumberFormat", NumberFormat, (Locale,), arg0)
end

function parse_and_ignore_whitespace(::Type{CompositeFormat}, arg0::JString, arg1::ParsePosition)
    return jcall(CompositeFormat, "parseAndIgnoreWhitespace", void, (JString, ParsePosition), arg0, arg1)
end

function parse_fixedstring(::Type{CompositeFormat}, arg0::JString, arg1::JString, arg2::ParsePosition)
    return jcall(CompositeFormat, "parseFixedstring", jboolean, (JString, JString, ParsePosition), arg0, arg1, arg2)
end

function parse_next_character(::Type{CompositeFormat}, arg0::JString, arg1::ParsePosition)
    return jcall(CompositeFormat, "parseNextCharacter", jchar, (JString, ParsePosition), arg0, arg1)
end

function parse_number(::Type{CompositeFormat}, arg0::JString, arg1::NumberFormat, arg2::ParsePosition)
    return jcall(CompositeFormat, "parseNumber", Number, (JString, NumberFormat, ParsePosition), arg0, arg1, arg2)
end