text:
struct abcdJ
    a::Int64
    b::Int64
    c::Int64
    d::Int64
    tJ::Int64
end
function make_int(p_sps,n_sps)
    count = 0
    label_T1 = abcdJ[ ]
    label_T0 = abcdJ[ ]
    tV = [ ]
    for (i,tmp) in enumerate(p_sps)
        n1 = tmp[1]; l1=tmp[2]; j1=tmp[3]
        for (j,tmp2) in enumerate(p_sps)
            if j<i;continue;end
            n2 = tmp2[1]; l2=tmp2[2]; j2=tmp2[3]
            for tJ in div(abs(j1-j2),2):1:div(j1+j2,2)
                if i==j && tJ%2==1;continue;end
                push!(tV,[i,j,tJ])
            end
        end
    end
    for i = 1:length(tV)
        for j = i:length(tV)
            if tV[i][3] != tV[j][3];continue;end
            a,b = tV[i][1:2]
            c,d = tV[j][1:2]
            tJ = tV[j][3]
            push!(label_T1,abcdJ(a,b,c,d,tJ))
            count += 1
        end
    end
    count = 0
    tV = [ ]
    ofst = length(p_sps)
    for (i,tmp) in enumerate(p_sps)
        n1 = tmp[1]; l1=tmp[2]; j1=tmp[3]
        i = i
        for (j,tmp2) in enumerate(n_sps)
            j = j
            if j < i;continue;end
            n2 = tmp2[1]; l2=tmp2[2]; j2=tmp2[3]
            for tJ in div(abs(j1-j2),2):1:div(j1+j2,2)
                if i==j && tJ%2==0;continue;end
                push!(tV,[i,j,tJ])
            end
        end
    end
    for i = 1:length(tV)
        for j = i:length(tV)
            if tV[i][3] != tV[j][3];continue;end
            a,b = tV[i][1:2]
            c,d = tV[j][1:2]
            tJ = tV[j][3]
            count += 1
            push!(label_T0,abcdJ(a,b,c,d,tJ))
        end
    end
    return label_T1,label_T0
end

function int2snt(p_sps,n_sps,label_T1,label_T0,olabels)
    l_TBME = length(olabels)
    ofst = length(p_sps)
    s_from_i = [ Int64[-1,-1] for i = 1:l_TBME]
    facs = [ Float64[0.0,0.0] for i = 1:l_TBME]
    for (j,olabel) in enumerate(olabels)
        a,b,c,d,totJ = @views olabel[1:5]
        fac0 = 0.0; fac1 = 0.0;  rank = "--"

        if (a <= ofst && b <= ofst && c <= ofst && d <= ofst) || (
            a > ofst && b > ofst && c > ofst && d > ofst)
            fac1 = 1.0
            rank = "pp/nn"
            for (i,label) in enumerate(label_T1)
                ta = label.a; tb = label.b; tc = label.c; td = label.d; tJ = label.tJ
                if tJ != totJ; continue;end
                if (a == ta && b == tb && c == tc && d == td ) || (
                    a == ta +ofst && b == tb+ofst && c == tc+ofst && d == td+ofst)
                    s_from_i[j][1] = i
                    facs[j][1] = fac1
                    break
                end
            end
        else
            rank = "pn"
            for (i,label) in enumerate(label_T1)
                ta = label.a; tb = label.b; tc = label.c; td = label.d; tJ = label.tJ
                if totJ != tJ; continue;end
                coef = 0.5 * sqrt( (1+deltaf(ta,tb) * (-1)^tJ) * (1+deltaf(tc,td) * (-1)^tJ))
                if (a == ta && b -ofst == tb && c == tc && d-ofst == td) || (
                    a == tc && b -ofst == td && c == ta && d-ofst == tb) #abcd or cdab
                    fac1 = coef
                elseif (a == tb && b -ofst == ta && c == tc && d-ofst == td) || (
                    a == tc && b -ofst == td && c == tb && d-ofst == ta) #bacd or cdba
                    ja = p_sps[tb][3]; jb = p_sps[ta][3]
                    fac1 = (-1)^(div(ja+jb,2)+tJ+1) .* coef
                elseif (a == ta && b -ofst == tb && c == td && d-ofst == tc) || (
                    a == td && b -ofst == tc && c == ta && d-ofst == tb) #abdc or dcab
                    jc = p_sps[td][3]; jd = p_sps[tc][3]
                    fac1 = (-1)^(div(jc+jd,2)+tJ+1) .* coef
                elseif (a == tb && b -ofst == ta && c == td && d-ofst == tc) || (
                    a == td && b -ofst == tc && c == tb && d-ofst == ta) #badc or dcba
                    ja = p_sps[tb][3]; jb = p_sps[ta][3]
                    jc = p_sps[td][3]; jd = p_sps[tc][3]
                    fac1 = (-1)^(div(ja+jb+jc+jd,2)) .* coef
                else
                    continue
                end
                s_from_i[j][1]=i
                facs[j][1] = fac1
                break
            end
            for (i,label) in enumerate(label_T0)
                ta = label.a; tb = label.b; tc = label.c; td = label.d; tJ = label.tJ
                if totJ != tJ; continue;end
                coef = 0.5 * sqrt( (1-deltaf(ta,tb) * (-1)^tJ) * (1-deltaf(tc,td) * (-1)^tJ))
                if (a == ta && b-ofst == tb && c == tc && d-ofst == td) || (
                    a == tc && b-ofst == td && c == ta && d-ofst == tb)  #abcd or cdab
                    fac0 = coef
                elseif (a == tb && b-ofst == ta && c == tc && d-ofst == td) || (
                    a == tc && b-ofst == td && c == tb && d-ofst == ta) #bacd or cdba
                    ja = p_sps[tb][3]; jb = p_sps[ta][3]
                    fac0 = (-1)^(div(ja+jb,2)+tJ) * coef
                elseif (a == ta && b-ofst == tb && c == td && d-ofst == tc) || (
                    a == td && b-ofst == tc && c == ta && d-ofst == tb) #abdc or dcab
                    jc = p_sps[td][3]; jd = p_sps[tc][3]
                    fac0 = (-1)^(div(jc+jd,2)+tJ) * coef
                elseif (a == tb && b-ofst == ta && c == td && d-ofst == tc) || (
                    a == td && b-ofst == tc && c == tb && d-ofst == ta) #badc or dcba
                    ja = p_sps[tb][3]; jb = p_sps[ta][3]
                    jc = p_sps[td][3]; jd = p_sps[tc][3]
                    fac0 = (-1)^(div(ja+jb+jc+jd,2)) * coef
                else
                    continue
                end
                s_from_i[j][2]=i
                facs[j][2] = fac0
                break
            end
        end
    end
    return s_from_i,facs
end
function readVint(intf,label_T1,label_T0,TF=true) # ad hoc func. to read from **.int
    f = open(intf,"r");tlines = readlines(f);close(f)
    lines = rm_comment(tlines)
    SPEs = Float64[ ]
    ofst = 0
    V0s = zeros(Float64,length(label_T0))
    V1s = zeros(Float64,length(label_T1))

    if intf == ""
        SPEs = zeros(Float64,length(p_sps))
        for i = 1:length(label_T0);V0s[i] = randn(); end
        for i = 1:length(label_T1);V1s[i] = randn(); end
        return SPEs,V1s,V0s
    end
    for line in lines
        if length(split(line)) < 7
            push!(SPEs,parse(Float64,split(line)[2]))
            continue
        end
        if TF;ofst = -length(SPEs);end
        tmp = split(line)
        a = parse(Int64,tmp[1])+ofst
        b = parse(Int64,tmp[2])+ofst
        c = parse(Int64,tmp[3])+ofst
        d = parse(Int64,tmp[4])+ofst
        J = parse(Int64,tmp[5])
        T = parse(Int64,tmp[6])
        V = parse(Float64,tmp[7])
        if T==0
            for (idx,ref) in enumerate(label_T0)
                ta = ref.a; tb = ref.b; tc = ref.c; td = ref.d; tJ = ref.tJ
                if a == ta && b==tb && c==tc && d==td && tJ== J
                    V0s[idx] = V
                    break
                end
            end
        else
            for (idx,ref) in enumerate(label_T1)
                ta = ref.a; tb = ref.b; tc = ref.c; td = ref.d; tJ = ref.tJ
                if a == ta && b==tb && c==tc && d==td && tJ== J
                    V1s[idx] = V
                    break
                end
            end
        end
    end
    return SPEs,V1s,V0s
end
function evalVsnt!(Vsnt,SPEs,V1s,V0s,idx_s_from_i,facs)
    ## adhoc: the same model space for proton/neutron
    #Vsnt = zeros(Float64,length(SPEs)*2 + length(idx_s_from_i))
    ofst = length(SPEs)
    for i = 1:length(SPEs)
        Vsnt[i] = SPEs[i]
        Vsnt[i+ofst] = SPEs[i]
    end
    ofst = 2*length(SPEs)
    fac0 = 0.0; fac1 = 0.0
    for (n,idxs) in enumerate(idx_s_from_i)
        idx_T1,idx_T0 = idxs
        Vsnt[ofst+n] = 0.0
        if idx_T1 != -1
            fac1 = facs[n][1]
            Vsnt[ofst+n] += V1s[idx_T1] .* fac1
        end
        if idx_T0 != -1
            fac0 = facs[n][2]
            Vsnt[ofst+n] += V0s[idx_T0] .* fac0
        end
    end
    return nothing
end

function evalVint!(Vint,SPEs,V1s,V0s)
    ln = length(SPEs)
    @inbounds @simd for i = 1:ln
        Vint[i] = SPEs[i]
    end
    @inbounds @simd for i = 1:length(V1s)
        Vint[ln+i] = V1s[i]
    end
    ln += length(V1s)
    @inbounds @simd for i = 1:length(V0s)
        Vint[ln+i] = V0s[i]
    end
    return nothing
end

function proposal_Vint!(nSPEs,nV1s,nV0s,SPEs,V1s,V0s,varM)
    @inbounds @simd for i = 1:length(SPEs)
        nSPEs[i] = SPEs[i] + varM * randn()
    end
    @inbounds @simd for i = 1:length(V1s)
        nV1s[i] = V1s[i] + varM * randn()
    end
    @inbounds @simd for i = 1:length(V0s)
        nV0s[i] = V0s[i] + varM * randn()
    end
    return nothing
end

# function test()
#     sntf = "../snts/usdb.snt"
#     Anum = 28
#     lp,ln,cp,cn,massop,Aref,pow,p_sps,n_sps,SPEs,olabels,oTBMEs,labels,TBMEs = readsnt(sntf,Anum)

#     label_T1,label_T0 = make_int(p_sps,n_sps)
#     idx_s_from_i,facs = int2snt(p_sps,n_sps,label_T1,label_T0,olabels)

#     intf = "../snts/random_input/ints/sdshl/tmp_1.int"
#     SPEs,V1s,V0s = readVint(intf,label_T1,label_T0)
#     Vsnt = evalVsnt(SPEs,V1s,V0s,idx_s_from_i,facs)
# end
#test()


text:
using Yao, Zygote, YaoFlux, Flux.Optimise

t = ArrayReg(bit"0") + ArrayReg(bit"1")
normalize!(t)

function fid(xs)
    r = zero_state(1)
    U = mat(chain(Rx(xs[1]), Rz(xs[2])))
    return abs(statevec(t)' * U * statevec(r))
end

function train()
    opt = ADAM()
    xs = rand(2)
    for _ in 1:1000
        println(fid(xs))
        Optimise.update!(opt, xs, -fid'(xs))
    end
    return xs
end

train()

ext: jl


text:
primitive type Pixel   <: Signed   64 end
# A constructor to create values of the type MyInt8.
Pixel(n::Int64) = reinterpret(Pixel, n)
_value(p::Pixel) = reinterpret(Int64, p)
valuetype(p::Pixel) = Int64
valuetype(p::Type{Pixel}) = Int64

# Global for multiplying
px = Pixel(1)

# This allows the REPL to show values of type MyInt8.
Base.show(io :: IO, x :: Pixel) = print(io, "$(_value(x))px")

# Pixel × Pixel => Pixel
for (op) in (:+, :-, :rem, :mod, :div, :fld, :cld)
    @eval begin
        import Base.($op)
        function ($op)(x::Pixel, y::Pixel)
            Pixel((valuetype(Pixel))(($op)(_value(x),_value(y))))
        end
    end
end

# Pixel × ... => Pixel
# TODO: How do i feel about float-to-pixel conversion? Currently throws InexactError.
# Relies on float->Pixel conversion
for op in (:*, :^, :/)
    @eval begin
        import Base.($op)
        function ($op)(x::Pixel, y::Number)
            Pixel((valuetype(Pixel))(($op)(_value(x),y)))
        end
    end
end

# ... × Pixel => Pixel
for op in (:*,)
    @eval begin
        import Base.($op)
        function ($op)(x::Number, y::Pixel)
            Pixel((valuetype(Pixel))(($op)(x,_value(y))))
        end
    end
end

# Pixel => Pixel
for op in (:zero, :one,)
    @eval begin
        import Base.($op)
        function ($op)(x::Pixel)
            Pixel((valuetype(Pixel))(($op)(_value(x))))
        end
    end
end

# Pixel × Pixel => ...
for op in (:(==), :!=, :>, :<, :(<=), :(>=), :cmp)
    @eval begin
        import Base.($op)
        function ($op)(x::Pixel, y::Pixel)
            ($op)(_value(x),_value(y))
        end
    end
end

# Pixel => ...
for op in (:sign, :signbit, :isinteger)
    @eval begin
        import Base.($op)
        function ($op)(x::Pixel)
            ($op)(_value(x))
        end
    end
end

import Base: read,write
function read(s::IO, ::Type{Pixel})
    r = read(s,valuetype(Pixel))
    Pixel(r)
end
function write(s::IO, p::Pixel)
    write(s, _value(p))
end

import Base: typemin, typemax
typemin(::Type{Pixel}) = typemin(valuetype(Pixel))
typemax(::Type{Pixel}) = typemax(valuetype(Pixel))

#import Base: convert, promote_rule
## TODO: Does conversion make sense? Probably not, right?
#convert(::Type{Pixel}, x::Pixel) = x
#convert(::Type{Pixel}, x::Integer) = Pixel(convert(valuetype(Pixel),x))
#
##convert(::Type{Integer}, x::Pixel) = (isinteger(x) ? convert(Integer, x.num) : throw(InexactError()))
##convert(::Type{T}, x::Pixel) where {T<:Integer} = (isinteger(x) ? convert(T, x.num) : throw(InexactError()))
#
##convert(::Type{AbstractFloat}, x::Pixel) = float(x.num)/float(x.den)
#convert(::Type{Pixel}, x::Float64) = convert(Pixel, valuetype(Pixel)(x))
#convert(::Type{Pixel}, x::Float32) = convert(Pixel, valuetype(Pixel)(x))
#
## TODO: I don't think promotion makes sense actually..
#promote_rule(::Type{Pixel}, ::Type{S}) where {S<:Integer} = Pixel
#promote_rule(::Type{Pixel}, ::Type{Pixel}) where {S<:Integer} = Pixel
#promote_rule(::Type{Pixel}, ::Type{S}) where {S<:AbstractFloat} = Pixel


text:
import LinearAlgebra
import FileIO
import GeometryBasics
import Observables
import Makie
import GLMakie
import DataFrames
import CSV
import ColorTypes
import FixedPointNumbers
import LinearAlgebra
import DifferentialEquations

using Porta

startframe = 1
frames = 3600
FPS = 60
resolution = (3840, 2160)
segments = 36
speed = 1
scale = 1.0
exportmode = ["gif", "frames", "video"]
exportmode = exportmode[2]
modelname = "gdp_alpha2"

# The scene object that contains other visual objects
#Makie.reasonable_resolution() = (800, 800)
scene = Makie.Scene(backgroundcolor = :white,
                    show_axis = false,
                    resolution = resolution)

fmap(b::S²) = b

"""
λ⁻¹map(p)
Sends a point on the plane back to a point on a unit sphere with the given
point. This is the inverse stereographic projection of a 3-sphere.
"""
function λ⁻¹map(p::S²)
    r3 = ℝ³(Cartesian(p))
    p₁, p₂, p₃ = vec(r3)
    mgnitude² = norm(r3)^2
    x₁ = 2p₁ / (1 + mgnitude²)
    x₂ = 2p₂ / (1 + mgnitude²)
    x₃ = 2p₃ / (1 + mgnitude²)
    x₄ = (-1 + mgnitude²) / (1 + mgnitude²)
    Quaternion(x₁, x₂, x₃, x₄)
end

"""
    getpointonpath(t)

Get a point on a path on S² with the given parameter `t`. t ∈ [0, 1].
"""
getpointonpath(t::Real) = begin
    t = t * 12π
    d = ℯ^cos(t) - 2cos(4t) - sin(t / 12)^5
    x, y = sin(t) * d, cos(t) * d
    ComplexLine(x + im * y)
end

"""
    getbutterflycurve(points)

Get butterfly curve by Temple H. Fay (1989) with the given number of `points`.
"""
function getbutterflycurve(points::Int)
    array = Array{ComplexLine,1}(undef, points)
    # 0 ≤ t ≤ 12π
    for i in 1:points
        t = (i - 1) / points * 12π
        d = ℯ^cos(t) - 2cos(4t) - sin(t / 12)^5
        x, y = sin(t) * d, cos(t) * d
        array[i] = ComplexLine(x + im * y)
    end
    array
end

"""
    Represents a line segment.

fields: p₁ and p₂.
"""
struct Line
    p₁::ℝ³
    p₂::ℝ³
end

function getdistance(point::ℝ³, line::Line)
    p₀ = point
    p₁, p₂ = line.p₁, line.p₂
    norm(cross(p₂ - p₁, p₁ - p₀)) / norm(p₂ - p₁)
end

function decimate(points::Array{Geographic,1}, ϵ::Float64)
    # Find the point with the maximum distance
    dmax = 0
    index = 1
    number = length(points)
    for i in 2:number-1
        line = Line(ℝ³(Cartesian(points[begin])), ℝ³(Cartesian(points[end])))
        point = ℝ³(Cartesian(points[i]))
        d = getdistance(point, line)
        if d > dmax
            index = i
            dmax = d
        end
    end

    array = []

    # If max distance is greater than epsilon, recursively simplify
    if dmax > ϵ
        # Recursive call
        array1 = decimate(points[begin:index], ϵ)
        array2 = decimate(points[index:end], ϵ)

        # Build the result list
        array = [array1[begin:end-1]; array2]
    else
        array = [points[begin]; points[end]]
    end
    # Return the result
    array
end

function getcolor(points::Array{Geographic,1}, color::Array{ColorTypes.RGBA{FixedPointNumbers.Normed{UInt8, 8}}, 2}, α::Float64)
    ϕ = map(x -> x.ϕ, points)
    θ = map(x -> x.θ, points)
    number = length(points)
    ϕ = sum(ϕ) / number
    θ = sum(θ) / number
    height, width = size(color)
    x = Int(floor((ϕ + π) / 2π * width))
    y = height - Int(floor((θ + π / 2) / π * height))
    r, g, b, a = color[y, x].r, color[y, x].g, color[y, x].b, color[y, x].alpha
    Makie.RGBAf(r, g, b, α)
end

# Use QGIS to design a geo map
colorref = FileIO.load("data/gdp/map_colors.png")
color = FileIO.load("data/gdp/map_transparent2.png")
# Made with Natural Earth.
# Free vector and raster map data @ naturalearthdata.com.

# 1. Download cultural data admin 0 from natural earth data
# 2. Install qgis desktop
# 3. Install mmqgis plugin for exporting to CSV

# The path to the dataset
attributespath = "data/gdp/geometry-attributes.csv"
attributes = DataFrames.DataFrame(CSV.File(attributespath))
attributes = DataFrames.sort(attributes, :shapeid, rev = true)
nodespath = "data/gdp/geometry-nodes.csv"
nodes = DataFrames.DataFrame(CSV.File(nodespath))

attributesgroup = DataFrames.groupby(attributes, :NAME)
nodesgroup = DataFrames.groupby(nodes, :shapeid)
number = length(attributesgroup)
ϵ = 5e-3
countries = Dict("shapeid" => [], "name" => [], "gdpmd" => [],
                 "gdpyear" => [], "economy" => [], "partid" => [], "nodes" => [])
for i in 1:number
    shapeid = attributesgroup[i][!, :shapeid][1]
    name = attributesgroup[i][!, :NAME][1]
    gdpmd = attributesgroup[i][!, :GDP_MD][1]
    gdpyear = attributesgroup[i][!, :GDP_YEAR][1]
    economy = attributesgroup[i][!, :ECONOMY][1]
    subdataframe = nodes[nodes.shapeid .== shapeid, :]
    uniquepartid = unique(subdataframe[!, :partid])
    ϵ = name == "Antarctica" ? 1e-3 : 5e-3
    histogram = Dict()
    for id in uniquepartid
        sub = subdataframe[subdataframe.partid .== id, :]
        ϕ = sub.x ./ 180 .* π
        θ = sub.y ./ 180 .* π
        coordinates = map(x -> Geographic(1, x[1], x[2]), eachrow([ϕ θ]))[begin:end-1]
        coordinates = decimate(coordinates, ϵ)
        histogram[id] = length(coordinates)
    end
    partsnumber = max(values(histogram)...)
    index = findfirst(x -> histogram[x] == partsnumber, uniquepartid)
    partid = uniquepartid[index]
    subdataframe = subdataframe[subdataframe.partid .== partid, :]
    ϕ = subdataframe.x ./ 180 .* π
    θ = subdataframe.y ./ 180 .* π
    coordinates = map(x -> Geographic(1, x[1], x[2]), eachrow([ϕ θ]))[begin:end-1]
    println("Length of points: $name : $(length(coordinates))")
    coordinates = decimate(coordinates, ϵ)
    println("Length of points: $name : $(length(coordinates))")
    push!(countries["shapeid"], shapeid)
    push!(countries["name"], name)
    push!(countries["gdpmd"], gdpmd)
    push!(countries["gdpyear"], gdpyear)
    push!(countries["economy"], economy)
    push!(countries["partid"], partid)
    push!(countries["nodes"], coordinates)
end

for i in 1:length(countries["nodes"])
    println(length(countries["nodes"][i]))
end

countries["gdpmd"] = countries["gdpmd"] ./ (max(countries["gdpmd"]...) - min(countries["gdpmd"]...))

s3rotation = Quaternion(1, 0, 0, 0)
config = Biquaternion(ℝ³(0, 0, 0))
solidwhirls = []
ghostwhirls = []
phases = Array{U1,2}(undef, number, 4)
for i in 1:number
    name = countries["name"][i]
    points = countries["nodes"][i]
    gdpmd = countries["gdpmd"][i]
    α₁ = name == "Antarctica" ? 0.1 : 1.0
    α₂ = name == "Antarctica" ? 0.05 : 0.1
    solidcolor = getcolor(points, colorref, α₁)
    ghostcolor = getcolor(points, colorref, α₂)
    solidtop = U1(0)
    solidbottom = U1(gdpmd * 2π)
    ghosttop = solidbottom
    ghostbottom = solidbottom
    phases[i, 1] = solidtop
    phases[i, 2] = solidbottom
    phases[i, 3] = ghosttop
    phases[i, 4] = ghostbottom
    #segments2 = name == "Antarctica" ? 2segments : segments
    solidwhirl = Whirl(scene,
                       points,
                       σmap,
                       fmap,
                       top = solidtop,
                       bottom = solidbottom,
                       s3rotation = s3rotation,
                       config = config,
                       segments = segments,
                       color = solidcolor,
                       transparency = false)
    push!(solidwhirls, solidwhirl)
    ghostwhirl = Whirl(scene,
                       points,
                       σmap,
                       fmap,
                       top = ghosttop,
                       bottom = ghostbottom,
                       s3rotation = s3rotation,
                       config = config,
                       segments = segments,
                       color = ghostcolor,
                       transparency = true)
    push!(ghostwhirls, ghostwhirl)
end

framesprites = []
framesprite1 = Frame(scene,
                     U1(0),
                     σmap,
                     fmap,
                     color,
                     s3rotation = s3rotation,
                     config = config,
                     segments = 2segments,
                     transparency = true)
framesprite2 = Frame(scene,
                     U1(0),
                     σmap,
                     fmap,
                     color,
                     s3rotation = s3rotation,
                     config = config,
                     segments = 2segments,
                     transparency = true)
push!(framesprites, framesprite1)
push!(framesprites, framesprite2)

a, b = ComplexLine(1 + im), ComplexLine(0.5 - 1.5 * im)
c, d = ComplexLine(-0.5 + 1 * im), ComplexLine(-1 + 0.3 * im)
f(z, a, b, c, d) = (a * z + b) / (c * z + d)

"""
    animate(i)

Update the state of observables with the given frame number `i`.
"""
function animate(i::Int)
    factor = 0.666
    step = (i - 1) / frames
    u = getpointonpath(step)
    if isnan(vec(u)[1]) || isnan(vec(u)[2])
        u = Geographic(1, 0, π/2)
    end
    τ = step * speed * 2pi
    ω = 2τ
    u = ℝ³(Cartesian(u))
    q = Quaternion(τ, u)
    a′, b′, c′, d′ = [ComplexLine(Cartesian(rotate(ℝ³(Cartesian(point)), q))).z for point in (a, b, c, d)]
    newmap(x::S²) = begin
        z = ComplexLine(x).z
        z = ComplexLine(f(z, a′, b′, c′, d′))
        if isnan(vec(z)[1]) || isnan(vec(z)[2])
            z = ComplexLine(Geographic(1, 0, -π / 2))
        end
        r, ϕ, θ = vec(Geographic(z))
        if isnan(r)
            z = ComplexLine(Geographic(1, ϕ, θ * 0.99))
        end
        z
    end
    newmap2(x::S²) = begin
        g = Geographic(x)
        r, ϕ, θ = vec(g)
        Geographic(r, ϕ + ω, θ)
    end
    for index in 1:number
        gdpmd = countries["gdpmd"][index] * (factor * 2π)
        solidtop = U1(0)
        solidbottom = U1(min(τ, gdpmd))
        ghosttop = solidbottom
        ghostbottom = τ > gdpmd ? U1(min(τ, 2π * factor)) : ghosttop
        phases[index, 1] = solidtop
        phases[index, 2] = solidbottom
        phases[index, 3] = ghosttop
        phases[index, 4] = ghostbottom
    end
    for (index, item) in enumerate(solidwhirls)
        item.s2tos2map = newmap2
        solidtop = phases[index, 1]
        solidbottom = phases[index, 2]
        points = countries["nodes"][index]
        update(item, points, s3rotation, solidtop, solidbottom)
    end
    for (index, item) in enumerate(ghostwhirls)
        item.s2tos2map = newmap2
        ghosttop = phases[index, 3]
        ghostbottom = phases[index, 4]
        points = countries["nodes"][index]
        update(item, points, s3rotation, ghosttop, ghostbottom)
    end
    update(framesprite1, σmap, newmap2)
    framesprite2.s2tos2map = newmap2
    update(framesprite2, U1(min(factor * 2π, τ)))

    n = ℝ³(0, 0, 1)
    v = ℝ³(-1, -1, sin(τ)^2) * (π / 2)
    v = rotate(v, Quaternion(π / 8 + π / 8 * sin(-τ), ℝ³(0, 0, 1)))
    # v = rotate(v, Quaternion(-π / 8, ℝ³(0, 0, 1)))
    # update eye position
    # scene.camera.eyeposition.val
    upvector = GeometryBasics.Vec3f0(vec(n)...)
    eyeposition = GeometryBasics.Vec3f0(vec(v)...)
    lookat = GeometryBasics.Vec3f0(0, 0, 0)
    #camera = Makie.Camera3D(scene)
    # Makie.update_cam!(scene, camera, eyeposition, lookat, upvector)
    Makie.update_cam!(scene, eyeposition, lookat, upvector)
    scene.center = false # prevent scene from recentering on display
end

function saveframe(filepath, scene, resolution)
    FileIO.save(filepath,
                scene;
                resolution = resolution,
                pt_per_unit = 400.0,
                px_per_unit = 400.0)
end

n = ℝ³(0, 0, 1)
v = ℝ³(-1, -1, 0) * 2.0
v = rotate(v, Quaternion(π / 8, ℝ³(0, 0, 1)))
# v = rotate(v, Quaternion(-π / 8, ℝ³(0, 0, 1)))
# update eye position
# scene.camera.eyeposition.val
upvector = GeometryBasics.Vec3f0(vec(n)...)
eyeposition = GeometryBasics.Vec3f0(vec(v)...)
lookat = GeometryBasics.Vec3f0(0, 0, 0)
#camera = Makie.Camera3D(scene)
# Makie.update_cam!(scene, camera, eyeposition, lookat, upvector)
Makie.update_cam!(scene, eyeposition, lookat, upvector)
scene.center = false # prevent scene from recentering on display
if exportmode ∈ ["gif", "video"]
    outputextension = exportmode == "gif" ? "gif" : "mkv"
    Makie.record(scene, "gallery/$modelname.$outputextension",
                            framerate = FPS) do io
        for i in startframe:frames
            animate(i) # animate the scene
            Makie.recordframe!(io) # record a new frame
            step = (i - 1) / frames
            println("Completed step $(100step).\n")
        end
    end
elseif exportmode == "frames"
    directory = joinpath("gallery", modelname)
    !isdir(directory) && mkdir(directory)
    for i in startframe:frames
        start = time()
        animate(i) # animate the scene
        elapsed = time() - start
        println("Generating frame $i took $elapsed (s).")

        start = time()

        paddingnumber = length(digits(frames))
        imageid = lpad(i, paddingnumber, "0")
        imagename = "$(modelname)_$(resolution[2])p_$(FPS)fps_$imageid.jpeg"
        filepath = joinpath(directory, imagename)
        saveframe(filepath, scene, resolution)
        elapsed = time() - start
        println("Saving file $filepath took $elapsed (s).")

        step = (i - 1) / frames
        println("Completed step $(100step).\n")
        sleep(5)

        stitch() = begin
            part = i ÷ (frames / 5)
            exportdir = joinpath(directory, "export")
            !isdir(exportdir) && mkdir(exportdir)
            WxH = "$(resolution[1])x$(resolution[2])"
            #WxH = "1920x1080"
            commonpart = "$(modelname)_$(resolution[2])p_$(FPS)fps"
            inputname = "$(commonpart)_%0$(paddingnumber)d.jpeg"
            inputpath = joinpath(directory, inputname)
            outputname = "$(commonpart)_$(part).mp4"
            outputpath = joinpath(exportdir, outputname)
            command1 = `ffmpeg -y -f image2 -framerate $FPS -i $inputpath -s $WxH -pix_fmt yuvj420p $outputpath`
            run(command1)
        end

        if isapprox(i % (frames / 5), 0)
            sleep(5)
            stitch()
            sleep(5)
        end
        # break
    end
end

max_stars_repo_name: iamazadi/geometrizer

max_issues_repo_name: iamazadi/Porta

max_forks_repo_name: iamazadi/Porta


text:
export RestMethod

type RestMethod
    description::AbstractString
    f::Function
end

RestMethod(f::Function, d::AbstractString) = RestMethod(d, f)

(r::RestMethod)(args...; kargs...) = r.f(args...; kargs...)


text:
using Documenter, BioServices

makedocs(
    sitename= "BioServices.jl",
    pages = [
        "Home" => "index.md",
        "Manual" => [
            "EUtils" => "man/eutils.md",
            "UMLS" => "man/umls.md",
            "GGGenome" => "man/gggenome.md",
            "Home" => "index.md"
            ]
        ],
    authors = "The BioJulia Organisation",
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true")
)

deploydocs(
    repo 	     = "github.com/BioJulia/BioServices.jl.git",
    push_preview = true,
    devbranch    = "master"
)

max_stars_repo_name: BioJulia/BioServices

max_issues_repo_name: BioJulia/BioServices

max_forks_repo_name: BioJulia/BioServices


text:
# This file is a part of project JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/Materials.jl/blob/master/LICENSE

using FEMMaterials, Pkg, Documenter, Literate, Dates

# automatically generate documentation from examples

"""
    add_datetime(content)
Add page generation time to the end of the content.
"""
function add_datetime(content)
    line =  "\n# Page generated at " * string(DateTime(now())) * "."
    content = content * line
    return content
end

"""
    remove_license(content)
Remove licence strings from source file.
"""
function remove_license(content)
    lines = split(content, '\n')
    function islicense(line)
        occursin("# This file is a part of JuliaFEM.", line) && return false
        occursin("# License is MIT:", line) && return false
        return true
    end
    content = join(filter(islicense, lines), '\n')
    return content
end

function generate_docs(pkg)

    function preprocess(content)
        content = add_datetime(content)
        content = remove_license(content)
    end

    pkg_dir = dirname(dirname(pathof(pkg)))
    exampledir = joinpath(pkg_dir, "examples")
    outdir = joinpath(pkg_dir, "docs", "src", "examples")
    outdir_notebooks = joinpath(pkg_dir, "docs", "src", "notebooks")
    example_pages = []
    for example_file in readdir(exampledir)
        startswith(example_file, "example_") || continue
        Literate.markdown(joinpath(exampledir, example_file), outdir; documenter=true, preprocess=preprocess)
        Literate.notebook(joinpath(exampledir, example_file), outdir_notebooks; documenter=true, preprocess=preprocess, execute=true)
        generated_example_file = joinpath("examples", first(splitext(example_file)) * ".md")
        push!(example_pages, generated_example_file)
    end
    return example_pages

end

function generate_readme(pkg)
    pkg_dir = dirname(dirname(pathof(pkg)))
    function add_header(content)
        lines = read(joinpath(pkg_dir,"docs","readme_header.txt"), String)
        content = lines * content
    end
    function preprocess(content)
        content = add_datetime(content)
        content = remove_license(content)
        content = add_header(content)
    end
    fil = joinpath(pkg_dir, "examples", "example_3dbeam.jl")
    Literate.markdown(fil, pkg_dir; name="README", documenter=false, preprocess=preprocess)
end

example_pages = generate_docs(FEMMaterials)

makedocs(modules=[FEMMaterials],
         format = Documenter.HTML(),
         checkdocs = :all,
         sitename = "FEMMaterials.jl",
         pages = [
                  "index.md",
                  "Examples" => example_pages
                 ]
        )


text:
function main()
    K, S = parse.(Int, split(readline()))

    cnt = 0
    for X in 0:K, Y in 0:K
        if S - K ≤ X + Y ≤ S
            cnt += 1
        end
    end

    println(cnt)
end

main()


text:
mutable struct Impact{T} <: Bound{T}
    Nx::Adjoint{T,SVector{6,T}}
    offset::SVector{6,T}

    function Impact(body::Body{T}, normal::AbstractVector{T};offset::AbstractVector{T} = zeros(3)) where T
        normal = normal / norm(normal)

        # Derived from plane equation a*v1 + b*v2 + distance*v3 = p - offset
        A = Array(svd(skew(normal)).V) # gives two plane vectors
        A[:,3] = normal # to ensure correct sign
        Ainv = inv(A)
        ainv3 = Ainv[3,:]
        Nx = [ainv3;0;0;0]'
        offset = [offset;0;0;0]

        new{T}(Nx, offset), body.id
    end
end

@inline function g(impact::Impact, body::Body, Δt, No)
    impact.Nx[SVector(1, 2, 3)]' * (getx3(body, Δt) - impact.offset[SVector(1, 2, 3)])
end

@inline ∂g∂pos(impact::Impact, No) = impact.Nx
@inline ∂g∂vel(impact::Impact, Δt, No) = impact.Nx * Δt

@inline function schurf(ineqc, impact::Impact, i, body::Body, μ, Δt, No)
    φ = g(impact, body, Δt, No)

    γ1 = ineqc.γ1[i]
    s1 = ineqc.s1[i]

    return impact.Nx' * (γ1 / s1 * φ - μ / s1)
end

@inline function schurD(ineqc, impact::Impact, i, body::Body, Δt)
    Nx = impact.Nx
    Nv = Δt * Nx

    γ1 = ineqc.γ1[i]
    s1 = ineqc.s1[i]

    return Nx' * γ1 / s1 * Nv
end


text:
function getInititeration(simulator::Quantum2dFractalSimulator)
	return simulator.inititeration
end


text:
using PrimitiveOneHot
using Documenter

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

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

deploydocs(;
    repo="github.com/chengchingwen/PrimitiveOneHot.jl",
)


text:
__precompile__()

"""
    Diversity package

The main **Diversity** package provides basic numbers-equivalent
diversity measures (described in
[Hill, 1973](http://www.jstor.org/stable/1934352)),
similarity-sensitive diversity measures (generalised from Hill, and
described in
[Leinster and Cobbold, 2012](http://www.esajournals.org/doi/abs/10.1890/10-2402.1)),
and related alpha, beta and gamma diversity measures at the level of
the metacommunity and its component subcommunities (generalised in
turn from Leinster and Cobbold, and described in
[Reeve et al, 2014](http://arxiv.org/abs/1404.6520)). The diversity
functions exist both with unicode names (e.g. ```ᾱ()```), which are
not automatically exported (as we feel they are too short) and with
matching longer ASCII names (e.g. ```NormalisedAlpha()```), which are.
We also provide functions to calculate appropriate
```subcommunityDiversity()``` and ```metacommunityDiversity()```
values for each measure, a general ```diversity()``` function for
extract any diversity measure at a series of scales.
"""
module Diversity

"""
    Diversity.API submodule

The `Diversity.API` submodule should be `import`ed if you want to create a
new type, partition or metacommunity subtype. Otherwise it can be
ignored.
"""
module API
include("API.jl")
# Base types
export AbstractPartition, AbstractTypes, AbstractMetacommunity

# Base class and functions required for each partition
export _getsubcommunitynames # required
export _countsubcommunities  # optional

# Base class and functions required for each type
export _gettypenames, _calcsimilarity, _getscale      # required
export _counttypes, _calcabundance, _calcordinariness # optional
export _getdiversityname, _addedoutputcols, _getaddedoutput # optional

# Base class and functions required for each metacommunity
export _gettypes, _getpartition, _getabundance  # required
export _getmetaabundance, _getweight            # optional
export _getordinariness!, _getmetaordinariness! # optional
export _hassimilarity

# Function with minimal default implementation for types and partitions
export floattypes
# Functions with standard implementations
export typematch, mcmatch
end

# Add in our user interface and the EcoBase higher level interface
include("Interface.jl")
export gettypes, getpartition
export counttypes, countsubcommunities
export gettypenames, getsubcommunitynames
export getabundance, getmetaabundance, getweight
export getordinariness!, getmetaordinariness!
export hassimilarity, calcsimilarity # Needed because it sometimes doesn't exist unless requested
export getdiversityname, addedoutputcols, getaddedoutput

# Inheritance from EcoBase
include("EcoBase.jl")

include("Iterators.jl")
export TypeIterator, SubcommunityIterator

include("Types.jl")
export UniqueTypes, Species, Taxonomy, GeneralTypes

include("Partition.jl")
export Subcommunities, Onecommunity

include("Metacommunity.jl")
export Metacommunity
export inddiv, subdiv, metadiv

include("EffectiveNumbers.jl")
export qD, qDZ

include("DiversityMeasure.jl")
export DiversityLevel
export individualDiversity, subcommunityDiversity, metacommunityDiversity
export DiversityMeasure, PowerMeanMeasure, RelativeEntropyMeasure
export RawAlpha, NormalisedAlpha
export RawBeta, NormalisedBeta, RawRho, NormalisedRho
export Gamma

"""
    Diversity.ShortNames submodule

We do not directly export `ᾱ`, `α`, `β̄`, `β`, `ρ̄`, `ρ`, `γ` as
they're too short. `γ` actually can't be exported like this -
it'll always just be `Shortnames.γ`, so we export `Γ` instead.

All of these can only be accessed via `Diversity.ShortNames`.
"""
module ShortNames
using Diversity

const α = RawAlpha
const ᾱ = NormalisedAlpha
const β = RawBeta
const β̄ = NormalisedBeta
const ρ = RawRho
const ρ̄ = NormalisedRho
const γ = Gamma

export α, ᾱ, β, β̄, ρ, ρ̄
const Γ = Gamma
export Γ

end

export getName, getASCIIName, getFullName

# Does PhyloBranches need to exist already?
using Requires
function __init__()
    @require Phylo="aea672f4-3940-5932-aa44-993d1c3ff149" begin
        include("Phylogenetics.jl")
        export AbstractPhyloTypes, PhyloBranches, PhyloDistances
        @deprecate(PhyloTypes, PhyloBranches)
    end
end

include("GeneralisedDiversities.jl")
export diversity
export norm_sub_alpha, raw_sub_alpha, norm_sub_beta, raw_sub_beta
export norm_sub_rho, raw_sub_rho, sub_gamma
export norm_meta_alpha, raw_meta_alpha, norm_meta_beta, raw_meta_beta
export norm_meta_rho, raw_meta_rho, meta_gamma

"""
    Diversity.Ecology submodule

The **Diversity.Ecology** module replicates old ecological
diversity measures and generalised versions of them that relate to our
general measures of alpha, beta and gamma diversity at subcommunity
and metacommunity levels. The generalisations of the richness, Shannon
and Simpson are the only standard measures we are aware of whose
subcommunity components sum directly to the corresponding ecosystem
measure (although note that Simpson's index decreases for increased
diversity, so small components are more diverse).
"""
module Ecology

include("Ecology.jl")
export generalisedrichness, richness
export generalisedshannon, shannon
export generalisedsimpson, simpson
export generalisedjaccard, jaccard

end # sub-module Ecology

"""
    Diversity.Jost submodule

Lou Jost's
[diversity](http://dx.doi.org/10.1111/j.2006.0030-1299.14714.x)
[measures](http://www.esajournals.org/doi/abs/10.1890/06-1736.1) are
found in the **Diversity.Jost** module.
"""
module Jost

include("Jost.jl")
export jostbeta, jostalpha
@deprecate(jostα, jostalpha)
@deprecate(jostβ, jostbeta)

end # sub-module Jost

"""
    Diversity.Hill submodule

[Hill numbers](http://www.jstor.org/stable/1934352) are found in the
**Diversity.Hill** package.
"""
module Hill

include("Hill.jl")
export hillnumber

end # sub-module Hill

# Path into package
path(path...; dir::String = "test") = joinpath(@__DIR__, "..", dir, path...)

end # module


text:
"""
Piecewise Legendre polynomial.
"""

import PyCall: PyObject

struct PiecewiseLegendrePoly
    o::PyObject
end

Base.convert(::Type{PiecewiseLegendrePoly}, o::PyObject) = PiecewiseLegendrePoly(o)
(p::PiecewiseLegendrePoly)(x) = p.o(x)

function Base.getindex(p::PiecewiseLegendrePoly, i::Int64)
    1 <= i <= p.o.size || error("Invalid index!")
    PiecewiseLegendrePoly(p.o.__getitem__(i-1))
end

"""
f: Function-like object
   f should takes a real-valued vector and return a real-valued D-dimensional arrays.
   axis option specifies which axis of this array correponds to the axis of the input.

axis:
   By default, the last axis is used for computing overlap.
"""
function overlap(poly::PiecewiseLegendrePoly, f)::Float64
    return poly.o.overlap(f)[1]
end

"""
PiecewiseLegendreFT
"""
struct PiecewiseLegendreFT
    o::PyObject
end

Base.convert(::Type{PiecewiseLegendreFT}, o::PyObject) = PiecewiseLegendreFT(o)
(p::PiecewiseLegendreFT)(x) = p.o(x)

function Base.getindex(p::PiecewiseLegendreFT, i::Int64)
    1 <= i <= p.o.size || error("Invalid index!")
    PiecewiseLegendreFT(p.o.__getitem__(i - 1))
end


text:
using SeqLoggers
using Test
using LoggingExtras

@testset "Utils" begin
@test SeqLoggers.stringify(a="test") == "\"a\":\"test\""
@test SeqLoggers.stringify(a="testa", b="testb") == "\"a\":\"testa\",\"b\":\"testb\""

@test SeqLoggers.to_seq_level(SeqLoggers.Logging.Debug) == "Debug"
@test SeqLoggers.to_seq_level(SeqLoggers.Logging.Info) == "Info"
@test SeqLoggers.to_seq_level(SeqLoggers.Logging.Warn) == "Warning"
@test SeqLoggers.to_seq_level(SeqLoggers.Logging.Error) == "Error"

url1 = SeqLoggers.joinurl("http://localhost:8080", "index.html")
url2 = SeqLoggers.joinurl("http://localhost:8080/", "index.html")
url3 = SeqLoggers.joinurl("http://localhost:8080/", "/index.html")
url4 = SeqLoggers.joinurl("http://localhost:8080", "/index.html")
@test url1 == url2 == url3 == url4 == "http://localhost:8080/index.html"

invalidString = "Test\" \nBadString\\ \r"
validString = "Test' \\nBadString/ \\r"
@test SeqLoggers.replace_invalid_character(invalidString) == validString
end

@testset "SeqLogger" begin

minimalSeqLogger = SeqLogger("http://localhost:8080";)
@test minimalSeqLogger isa SeqLogger
@test minimalSeqLogger.min_level == Logging.Info
minimalSeqLogger.server_url == "http://localhost:5341/api/events/raw"
@test minimalSeqLogger.event_properties[]  == ""
@test minimalSeqLogger.headers == Dict("Content-Type" => "application/vnd.serilog.clef")
@test minimalSeqLogger.batch_size == 10
@test minimalSeqLogger.event_batch == String[]

min_level = Logging.Debug
seqLogger = SeqLogger(
    "http://myhost:1010";
    min_level=min_level,
    api_key="Test",
    App="Trialrun",
    batch_size=3,
    HistoryId=raw"xxxxxxxx-xxxx-xxxx-xxxx-xxxxxxxxxxxx",
    Env="Test"
)

@test seqLogger.min_level == Logging.Debug
@test seqLogger.server_url == "http://myhost:1010/api/events/raw"
@test seqLogger.event_properties[]  ==  "\"App\":\"Trialrun\",\"HistoryId\":\"xxxxxxxx-xxxx-xxxx-xxxx-xxxxxxxxxxxx\",\"Env\":\"Test\""
@test seqLogger.headers == Dict(
    "Content-Type" => "application/vnd.serilog.clef",
    "X-Seq-api_key" => "Test"
)
@test seqLogger.batch_size == 3
@test seqLogger.event_batch == String[]

@test Logging.shouldlog(seqLogger) == true
@test Logging.min_enabled_level(seqLogger) == min_level
@test Logging.catch_exceptions(seqLogger) == false

end

@testset "Add log event properties dynamically" begin
logger = SeqLogger("http://localhost:5341";)
SeqLoggers.event_property!(logger; newProperty="DynamicProperty")
@test logger.event_properties[] == "\"newProperty\":\"DynamicProperty\""
SeqLoggers.event_property!(logger; next="true")
@test logger.event_properties[] == "\"newProperty\":\"DynamicProperty\",\"next\":\"true\""
combinedLogger = TeeLogger(Logging.current_logger(), logger)
SeqLoggers.event_property!(combinedLogger; next3="true3")
@test combinedLogger.loggers[2].event_properties[] ==
    "\"newProperty\":\"DynamicProperty\",\"next\":\"true\",\"next3\":\"true3\""
end

@testset "parse_event_str_from_args" begin

logger = SeqLogger("http://localhost:5341";)
message_args = (
    level=Logging.Info,
    message="Log Message",
    _module=Main,
    group=:file_name,
    file="file_path",
    line=100,
    kwargs=Dict(:a => 1, :b => 2)
)
parsed_msg_str = SeqLoggers.parse_event_str_from_args(logger, message_args)
@test occursin("\"@mt\":\"Log Message\"", parsed_msg_str)
@test occursin("\"@l\":\"Info\"", parsed_msg_str)
@test occursin("\"_file\":\"file_path\"", parsed_msg_str)
@test occursin("\"_line\":\"100\"", parsed_msg_str)
@test occursin("\"a\":\"1\",\"b\":\"2\"", parsed_msg_str)

message_args = (
    level=Logging.Error,
    message="Log Message",
    _module=Main,
    group=:file_name,
    file="file_path",
    line=100,
    kwargs=Dict(:back_trace => "this is a backtrace", :b => 2)
)
parsed_msg_str = SeqLoggers.parse_event_str_from_args(logger, message_args)
@test occursin("\"@mt\":\"Log Message\"", parsed_msg_str)
@test occursin("\"@l\":\"Error\"", parsed_msg_str)
@test occursin("\"@x\":\"this is a backtrace\"", parsed_msg_str)
@test occursin("\"_file\":\"file_path\"", parsed_msg_str)
@test occursin("\"_line\":\"100\"", parsed_msg_str)
@test occursin("\"b\":\"2\"", parsed_msg_str)

end

@testset "Run With logger: ConsoleLogger" begin
logger = ConsoleLogger(stderr, Logging.Info)
run_with_logger(logger) do
    @debug "No Exception"
    @info "No Exception"
    @warn "No Exception"
    @error "No Exception"
end

@test_throws DomainError run_with_logger(logger, -1) do x
    sqrt(x)
end

end


text:
include("Discounts.jl")
include("Cumulants.jl")
include("Policies.jl")

struct GVF{C,D,π}
    cumulant::C
    discount::D
    policy::π
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

# http://ahojukka5.github.io/posts/finite-element-solution-for-one-element-problem/

using JuliaFEM
using JuliaFEM.Testing

@testset "test 2d linear elasticity local matrices" begin
    element = Element(Quad4, [1, 2, 3, 4])
    X = Dict{Int64, Vector{Float64}}(
        1 => [0.0, 0.0],
        2 => [1.0, 0.0],
        3 => [1.0, 1.0],
        4 => [0.0, 1.0])
    u = Dict{Int64, Vector{Float64}}(
        1 => [0.0, 0.0],
        2 => [0.0, 0.0],
        3 => [0.0, 0.0],
        4 => [0.0, 0.0])
    update!(element, "geometry", X)
    update!(element, "displacement", u)
    update!(element, "youngs modulus" =>  288.0, "poissons ratio" => 1/3)
    update!(element, "displacement load", DCTI([4.0, 8.0]))

    problem = Problem(Elasticity, "[0x1] x [0x1] block", 2)
    problem.properties.formulation = :plane_stress
    assemble!(problem, element)
    K = full(problem.assembly.K)
    f = vec(full(problem.assembly.f))

    K_expected = [
        144  54 -90   0 -72 -54  18   0
         54 144   0  18 -54 -72   0 -90
        -90   0 144 -54  18   0 -72  54
          0  18 -54 144   0 -90  54 -72
        -72 -54  18   0 144  54 -90   0
        -54 -72   0 -90  54 144   0  18
         18   0 -72  54 -90   0 144 -54
          0 -90  54 -72   0  18 -54 144]

    f_expected = [1, 2, 1, 2, 1, 2, 1, 2]

    @test isapprox(K, K_expected)
    @test isapprox(f, f_expected)
end


text:
@assert begin
    msh=FinVolMesh{1}()

    nnds=10

    pts=[add_point!(msh, [x]) for x in LinRange(0.0, 1.0, nnds)]

    fcs=[
        add_face!(msh, [pt]) for pt in pts
    ]

    vols=[
        add_volume!(msh, [fcs[i], fcs[i+1]]) for i=1:(length(fcs)-1)
    ]

    process_geometry!(msh)

    u_nodes=collect(LinRange(0.0, 1.0, nnds))
    u_cells=(u_nodes[1:(end-1)].+u_nodes[2:end])./2

    u_faces=Vector{Float64}(undef, nnds)

    face_average!(msh, u_cells, u_faces)

    eqs=isapprox.(u_faces, u_nodes)

    @assert all(eqs[2:(end-1)])

    u_faces[1]=0.0
    u_faces[end]=1.0

    grads=zeros(Float64, nnds-1)
    volume_gradient!(msh, grads, u_faces, 1)

    @assert isapprox(grads, ones(Float64, nnds-1); atol=1e-3)

    flux_ints=zeros(Float64, nnds-1)

    flux_integral!(msh, flux_ints, u_faces)
    @assert isapprox(flux_ints, -1.0*ones(Float64, nnds-1); atol=1e-3)

    true
end


text:
module TestSphere

include("../src/Sphere.jl")
import .Sphere

using Test

const tol = 5e-13

export test_swal_inner_product, test_norm_swal_lap, test_norm_filter

"""
Computes inner product between different Y^s_{lm}

test_swal_inner_product(
   ny::Int64,
   spin::Int64,
   m_ang::Int64
   )::Nothing
"""
function test_swal_inner_product(
      ny::Int64,
      spin::Int64,
      m_ang::Int64
   )::Nothing

   nl = Sphere.num_l(ny)

   vals = Sphere.swal_vals(ny,spin,m_ang)

   println("Test: Y^s_{lm} orthogonality:\tny=$ny\tspin=$spin\tm_ang=$m_ang\tl=[0..$nl]")
   for i=1:nl
      for j=i:nl
         val = Sphere.inner_product(vals[:,i],vals[:,j])
         if i==j
            @test abs(val-1.0) < tol
         else
            @test abs(val) < tol
         end
      end
   end
   return nothing
end

"""
Computes the norm of the difference between the numerical and
exact spin-weighted spherical laplacian for a given
spin-weighted spherical harmonic.

test_norm_swal_lap(
   ny::Int64,
   spin::Int64,
   m_ang::Int64,
   l_ang::Int64
   )::Nothing
"""
function test_norm_swal_lap(
      ny::Int64,
      spin::Int64,
      m_ang::Int64,
      l_ang::Int64
   )::Nothing

   Yv = Sphere.Y_vals(ny)

   swal = [Sphere.swal(spin,m_ang,l_ang,y) for y in Yv]

   swal_lap_v1 = [-(l_ang-spin)*(l_ang+spin+1.0)*v for v in swal]
   swal_lap_v2 = zeros(Float64,ny)

   lap = Sphere.swal_laplacian_matrix(ny,spin,m_ang)
   for j=1:ny
      for k=1:ny
         swal_lap_v2[j] += swal[k]*lap[k,j]
      end
   end

   ## compute integral over interval [-1,1] of difference
   n = Sphere.inner_product(ones(ny),swal_lap_v1 .- swal_lap_v2)

   println("Test: Δ_{S2}Y^s_{lm}:\tny=$ny\tspin=$spin\tm_ang=$m_ang\tl=$l_ang")
   @test abs(n) < tol

   return nothing
end

"""
Test that the swal filter matrix acts as a low pass filter
for angular perturbations.
"""
function test_norm_filter(
      ny::Int64,
      spin::Int64,
      m_ang::Int64,
      l_ang::Int64
   )::Nothing

   Yv = Sphere.Y_vals(ny)

   ## essentially random points multiplied by spherical harmonic;
   ## with the filter the norm should go to zero
   v1 = [0.1*(rand()-0.5)*Sphere.swal(spin,m_ang,l_ang,y) for y in Yv]
   v2 = zeros(Float64,ny)

   filter = Sphere.swal_filter_matrix(ny,spin,m_ang)

   for j=1:ny
      for k=1:ny
         v2[j] += v1[k]*filter[k,j]
      end
   end

   ## compute integral over interval [-1,1] of difference
   i1 = abs(Sphere.inner_product(ones(ny), abs.(v1)))
   i2 = abs(Sphere.inner_product(ones(ny), abs.(v2)))

   println("Test: filter(rand) smoothens:\tny=$ny\tspin=$spin\tm_ang=$m_ang\tl_ang=$l_ang")
   @test i2 < i1

   return nothing
end

end


text:
const K = 2
const N = 10
const α_neal = 10.0

@model function gmm_joint(z, x, ::Type{T}=Float64) where T
    if z === missing
        z = Vector{Int}(undef, N)
    end

    if x === missing
        x = Vector{T}(undef, N)
    end

    # Cluster centers.
    μ ~ filldist(Normal(), K)

    # Cluster association prior.
    w ~ Dirichlet(K, 1.0)

    # Cluster assignments.
    z ~ filldist(DiscreteNonParametric(1:K, w), N)

    # Observations.
    for n = 1:N
        x[n] ~ Normal(μ[z[n]], 1.0)
    end

    return z, x
end

@model function hmm_joint(s, x, ::Type{T}=Float64) where {T<:Real}
    if s === missing
        s = Vector{Int}(undef, N)
    end

    if x === missing
        x = Vector{T}(undef, N)
    end

    # Emission matrix.
    m = Vector{T}(undef, K)

    # Transition matrix.
    T = Vector{Vector{T}}(undef, K)

    # Assign distributions to each element
    # of the transition matrix and the
    # emission matrix.
    for i = 1:K
        T[i] ~ Dirichlet(K, 1.0)
        m[i] ~ Normal(i, 0.5)
    end

    # Observe each point of the input.
    s[1] ~ Categorical(K)
    x[1] ~ Normal(m[s[1]], 0.1)

    for i = 2:N
        s[i] ~ Categorical(T[s[i-1]])
        x[i] ~ Normal(m[s[i]], 0.1)
    end

    return s, x
end

function stickbreak(v)
    K = length(v) + 1
    cumprod_one_minus_v = cumprod(1 .- v)

    return map(1:K) do k
        if k == 1
            v[1]
        elseif k == K
            cumprod_one_minus_v[K - 1]
        else
            v[k] * cumprod_one_minus_v[k - 1]
        end
    end
end

# from https://luiarthur.github.io/TuringBnpBenchmarks/dpsbgmm
@model function imm_stick_joint(z, y, ::Type{T}=Float64) where {T}
    if z === missing
        z = Vector{Int}(undef, N)
    end

    if x === missing
        x = Vector{T}(undef, N)
    end

    crm = DirichletProcess(α_neal)

    v ~ filldist(StickBreakingProcess(crm), K - 1)
    w = stickbreak(v)

    # Cluster assignments
    for n = 1:N
        z[n] ~ Categorical(w)
    end

    # Cluster centers
    L = identity(K)
    μ ~ filldist(Normal(), L)

    # Observations
    for n = 1:N
        y[n] ~ Normal(μ[z[n]], 1.0)
    end
end

# test function
g(θ, x) = cat(θ, x; dims=1)

function geweke_gmm()
    rand_θ_given(x) = rand(StaticConditional(gmm_joint(missing, x), :z))

    test = perform(GewekeTest(2_000), gmm_joint, rand_θ_given, g)
    # plot(test, gmm_joint; size=(300, 300), title="GMM Geweke statistic")
    s = compute_statistic!(test, g)
    @test s.pvalue < 0.05
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "ffts"
version = v"0.9.0"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/anthonix/ffts.git", "fe86885ecafd0d16eb122f3212403d1d5a86e24e"),
    DirectorySource("./bundled")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
for f in ${WORKSPACE}/srcdir/patches/*.patch; do
    atomic_patch -p1 ${f}
done
cd ffts/
mkdir build
cd build
cmake -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} -DCMAKE_BUILD_TYPE=Release -DENABLE_SHARED=ON ..
make
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms(; experimental=true)
# filter failing ARM targets
filter!(p -> arch(p) ≠ "armv7l", platforms)
filter!(p -> arch(p) ≠ "armv6l", platforms)
filter!(p -> !(Sys.isapple(p) && arch(p) == "aarch64"), platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct("libffts", :libffts)
]

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

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6", preferred_gcc_version=v"7")

max_stars_repo_name: pepijndevos/Yggdrasil

max_issues_repo_name: pepijndevos/Yggdrasil

max_forks_repo_name: pepijndevos/Yggdrasil


text:
alpha = rand()+im*rand()
beta = rand()+im*rand()
#beta = 0.0+im*rand()

# Compute covariance matrix:
nt = 5
t = linspace(0,nt-1,nt)

cov = zeros(Complex{Float64},nt,nt)
#cov = zeros(Float64,nt,nt)
for i=1:nt
  for j=i:nt
    cov[i,j]=alpha*exp(-beta*(t[j]-t[i]))
    cov[j,i]=conj(cov[i,j])
  end
  cov[i,i]=real(alpha)
end

ccov=chol(cov)

#covinv=inv(cov)
#
#for i=1:nt
#  for j=1:nt
#    if abs(covinv[i,j]) < eps()
#       covinv[i,j]=0.0+im*0.0
##       covinv[i,j]=0.0
#    end
#  end
#end
#
#covinvchol=chol(covinv)

# Check that cholesky decomposition is correct:
*(ctranspose(ccov),ccov)-cov

# Now try constructing:
a=zeros(Complex64,nt)
b=zeros(Complex64,nt)
v=zeros(Complex64,nt)
matrix = zeros(Complex64,nt,nt)
for i=1:nt
  a[i]=sqrt(alpha)*exp(beta*t[i])
  b[i]=sqrt(alpha)*exp(-beta*t[i])
  v[i]=b[i]/a[i]
end
d = Diagonal(a)
l = zeros(Complex64,nt,nt)
lam=zeros(Complex64,nt,nt)
for i=1:nt
  for j=1:i
    l[i,j]=1.0+im*0.0
  end
  if i < nt
    lam[i,i] = sqrt(v[i]-v[i+1])
  else
    lam[i,i] = sqrt(v[i])
  end
  for j=1:nt
    if j >= i
      matrix[i,j] = a[i]*b[j]
    else
      matrix[i,j] = a[j]*b[i]
    end
  end
end
cholm = *(d,*(l,lam))

max_stars_repo_name: dfm/ess

max_issues_repo_name: dfm/ess

max_forks_repo_name: dfm/ess


text:
cd(@__DIR__)
ks, ctr, a1face, a2face, simTime = KitBase.initialize("config_2d2f.txt")
simTime = KitBase.solve!(ks, ctr, a1face, a2face, simTime)

res = zeros(4)
dt = KitBase.timestep(ks, ctr, 0.0)
KitBase.evolve!(ks, ctr, a1face, a2face, dt; mode = :kcu, bc = :maxwell)
KitBase.update!(ks, ctr, a1face, a2face, dt, zeros(4); coll = :bgk, bc = :extra)
KitBase.update!(ks, ctr, a1face, a2face, dt, zeros(4); coll = :bgk, bc = :period)
KitBase.update!(ks, ctr, a1face, a2face, dt, zeros(4); coll = :bgk, bc = :mirror)

ks, ctr, a1face, a2face, simTime = KitBase.initialize("config_2d1f.txt")
KitBase.reconstruct!(ks, ctr)
KitBase.update!(ks, ctr, a1face, a2face, dt, zeros(4); coll = :bgk, bc = :maxwell)
KitBase.evolve!(ks, ctr, a1face, a2face, dt; mode = :kfvs, bc = :maxwell)
KitBase.evolve!(ks, ctr, a1face, a2face, dt; mode = :kcu, bc = :maxwell)
KitBase.update!(ks, ctr, a1face, a2face, dt, zeros(4); coll = :bgk, bc = :extra)
KitBase.update!(ks, ctr, a1face, a2face, dt, zeros(4); coll = :bgk, bc = :period)
KitBase.update!(ks, ctr, a1face, a2face, dt, zeros(4); coll = :bgk, bc = :mirror)

KitBase.plot_contour(ks, ctr)


text:
module GradientBoosting

using DecisionTree
using StatsBase

function random_subset_selection(X_imp, y_imp, k, index_rm)

    n = size(X_imp,1)
    p = size(X_imp,2)

    # Train arrays
    train_Xs = Array{Float64}(undef, 0, p)
    train_ys = Array{Float64}(undef, 0)

    # Randomize
    shuffleray = hcat(y_imp, X_imp, index_rm)
    shuffleray = shuffleray[shuffle(1:end), :]
    y_shuffled = shuffleray[:,1]
    X_shuffled = shuffleray[:,2:(end-1)]
    index_shuffled = shuffleray[:,end]

    train_ys = y_shuffled[1:k]
    train_Xs = X_shuffled[1:k,:]
    index_selected = index_shuffled[1:k]

    return train_Xs, train_ys, index_selected
end

function fit(y_in, X_in, lr, max_depth, number_of_trees,
    n_subfeatures=0, min_samples_leaf = 5, min_samples_split = 2,
    min_purity_increase = 0.0, pruning_purity = 1.0,
    algo = "General")

    y = deepcopy(y_in)
    X = deepcopy(X_in)
    yi = deepcopy(y)
    Xi = deepcopy(X)

    # Initialize Error
    ri = 0
    n = length(yi)

    # Initialize Predictions with Average
    predf = reshape(zeros(n), :, 1)
    predf .= mean(y)

    # Number of Iterations
    iters = number_of_trees

    # Array for storing all trees in the final model
    finalmodels = Array{Union{Node{Float64,Float64},Leaf{Float64}}}(undef, iters,1)

    # Array for convergence plot
    predictions = zeros(iters,n)

    for i=1:iters

        if algo == "Stochastic"
            #select random subset of the Xs and ys
            Xi, yi, index_sel, check = random_subset_selection(X, y, k, indicies)
        end

        # Fit the Decision Tree
        model = DecisionTree.build_tree(yi, Xi,
                       n_subfeatures,
                       max_depth,
                       min_samples_leaf,
                       min_samples_split,
                       min_purity_increase)

        #add model to final models
        finalmodels[i] = model

        # Predict yi
        predi = apply_tree(model, Xi)

        # Compute the New Prediction and the New Residuals
        # Set the New yi equal to the Residuals
        predf .+= lr*predi
        ri = y_in.-predf
        yi = vec(ri)

        #record predictions for convergence plot
        predictions[i,:] = predf

    end

    return predf, finalmodels, predictions

end

function predict(y_in, X_in, lr, arrayofmodels)
    yi = deepcopy(y_in)
    Xi = deepcopy(X_in)

    n = size(Xi, 1)

    # Initialize Predictions with Average
    predf = reshape(zeros(n), :, 1);
    predf .= mean(yi);

    for i = 1:size(arrayofmodels,1)
        model = arrayofmodels[i]
        predi = apply_tree(model, Xi)
        predf .+= lr*predi
    end
    return predf
end

export fit, predict

end

max_stars_repo_name: ezuccarelli/GradientBoosting

max_issues_repo_name: ezuccarelli/GradientBoosting


text:
function evaluate(obj::LeastSquaresProblem, arg0::RealVector)
    return jcall(obj, "evaluate", LeastSquaresProblem_Evaluation, (RealVector,), arg0)
end

function get_convergence_checker(obj::OptimizationProblem)
    return jcall(obj, "getConvergenceChecker", ConvergenceChecker, ())
end

function get_evaluation_counter(obj::OptimizationProblem)
    return jcall(obj, "getEvaluationCounter", Incrementor, ())
end

function get_iteration_counter(obj::OptimizationProblem)
    return jcall(obj, "getIterationCounter", Incrementor, ())
end

function get_observation_size(obj::LeastSquaresProblem)
    return jcall(obj, "getObservationSize", jint, ())
end

function get_parameter_size(obj::LeastSquaresProblem)
    return jcall(obj, "getParameterSize", jint, ())
end

function get_start(obj::LeastSquaresProblem)
    return jcall(obj, "getStart", RealVector, ())
end


text:
# This file was generated, do not modify it. # hide
feature = ["AA", "BB", "AA", "AA", "BB"]
elscitype(feature)


text:
@testset "Cauchy Endgame" begin
    # Case: 1 regular solution
    @testset "Regular solution" begin
        @polyvar x y
        f = [x^2 - 4, x + 2y - 5]
        tracker, S = coretracker_startsolutions(f; seed = 1234, log_homotopy = true)
        track!(tracker, first(S), 0.0, 20.0)
        eg = HC.CauchyEndgame(tracker.state.x)
        p = copy(tracker.state.x)
        code, m = HC.predict!(p, tracker, eg)
        @test code == HC.CAUCHY_SUCCESS
        @test m == 1
        @test p ≈ [2, 1.5] atol = 1e-8
    end

    # Case 2: Univariate Singular
    @testset "Univariate singular" begin
        @polyvar x
        f = [(x - 3)^5]
        tracker, S = coretracker_startsolutions(f; seed = 1234, log_homotopy = true)
        track!(tracker, first(S), 0.0, 20.0)
        s = copy(tracker.state.x)
        eg = HC.CauchyEndgame(s)
        p = copy(tracker.state.x)
        code, m = HC.predict!(p, tracker, eg)
        @test code == HC.CAUCHY_SUCCESS
        @test m == 5
        @test p ≈ [3] atol = 1e-8
        # test that we can continute tracking from there
        track!(tracker, s, 20.0, 23.0)
        p .= 0
        code, m = HC.predict!(p, tracker, eg)
        @test code == HC.CAUCHY_SUCCESS
        @test m == 5
        @test p ≈ [3] atol = 1e-8

        # use in iterator
        f = [(x - 3)^5]
        tracker, S = coretracker_startsolutions(f; seed = 1234, log_homotopy = true)
        init!(tracker, first(S), 0.0, 25.0)
        # use as an iterator to do cauchy endgame in between,
        # make sure that we can still continue the tracking
        for _ in tracker
            if tracker.state.s > 14
                code, m = HC.predict!(p, tracker, eg)
                @test code == HC.CAUCHY_SUCCESS
                @test m == 5
            end
        end
        @test is_success(tracker.state.status)
        @test current_t(tracker) == 25
    end

    # Case 3
    @testset "Griewank-Osborne" begin
        # affine
        f = equations(griewank_osborne())
        tracker, starts = coretracker_startsolutions(f; seed = 1234, log_homotopy = true)
        S = collect(starts)
        track!(tracker, S[3], 0.0, 20.0)
        s = copy(tracker.state.x)
        eg = HC.CauchyEndgame(s)
        p = copy(tracker.state.x)
        code, m = HC.predict!(p, tracker, eg)
        @test p ≈ [0, 0] atol = 1e-8
        @test m == 3

        # projective
        f = equations(griewank_osborne())
        # random patch
        tracker, starts = coretracker_startsolutions(
            f;
            seed = 1234,
            projective_tracking = true,
            log_homotopy = true,
            patch = RandomPatch(),
        )
        S = collect(starts)
        track!(tracker, S[3], 0.0, 20.0)
        s = copy(tracker.state.x)
        eg = HC.CauchyEndgame(s)
        p = copy(tracker.state.x)
        code, m = HC.predict!(p, tracker, eg)
        @test ProjectiveVectors.affine_chart(p) ≈ [0, 0] atol = 1e-8
        @test m == 3

        # OrthogonalPatch
        tracker, starts = coretracker_startsolutions(
            f;
            seed = 1234,
            projective_tracking = true,
            log_homotopy = true,
            patch = OrthogonalPatch(),
        )
        S = collect(starts)
        track!(tracker, S[3], 0.0, 20.0)
        s = copy(tracker.state.x)
        eg = HC.CauchyEndgame(s)
        p = copy(tracker.state.x)
        code, m = HC.predict!(p, tracker, eg)
        @test ProjectiveVectors.affine_chart(p) ≈ [0, 0] atol = 1e-8
        @test m == 3
    end

    @testset "Limit accuracy error" begin
        # affine
        @polyvar x
        f = [(x - 3)^5]
        tracker, starts = coretracker_startsolutions(
            f;
            seed = 1234, min_step_size = eps()^2, log_homotopy = true,
        )
        S = collect(starts)
        tracker.options.precision = PRECISION_ADAPTIVE
        track(tracker, S[3], 0.0, 30.0)
        tracker.options.precision = PRECISION_FIXED_64
        s = copy(tracker.state.x)
        eg = HC.CauchyEndgame(s)
        p = copy(tracker.state.x)
        code, m = HC.predict!(p, tracker, eg)
        @test code == HC.CAUCHY_TERMINATED_ACCURACY_LIMIT
    end
end


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

module SimilaritySearch
abstract type AbstractSearchContext end

using Parameters

import Base: push!, append!
export AbstractSearchContext, SemiMetric, evaluate, search, searchbatch, getknnresult
include("distances/Distances.jl")

include("db.jl")
include("knnresult.jl")
include("knnresultshift.jl")
@inline Base.length(searchctx::AbstractSearchContext) = length(searchctx.db)
@inline Base.getindex(searchctx::AbstractSearchContext, i::Integer) = searchctx.db[i]
@inline Base.eachindex(searchctx::AbstractSearchContext) = 1:length(searchctx)
@inline Base.eltype(searchctx::AbstractSearchContext) = eltype(searchctx.db)
include("perf.jl")
include("sequential-exhaustive.jl")
include("parallel-exhaustive.jl")
include("opt.jl")
include("searchgraph/SearchGraph.jl")
include("allknn.jl")
include("neardup.jl")

const GlobalKnnResult = [KnnResult(32)]   # see __init__ function at the end of this file

"""
    getknnresult(k::Integer, pools=nothing) -> KnnResult

Generic function to obtain a shared result set for the same thread and avoid memory allocations.
This function should be specialized for indexes and pools that use shared results or threads in some special way.
"""
@inline function getknnresult(k::Integer, pools=nothing)
    res = @inbounds GlobalKnnResult[Threads.threadid()]
    reuse!(res, k)
end

"""
    searchbatch(index, Q, k::Integer=10; parallel=false, pools=GlobalKnnResult) -> indices, distances

Searches a batch of queries in the given index (searches for k neighbors).

- `parallel` specifies if the query should be solved in parallel at object level (each query is sequentially solved but many queries solved in different threads).
- `pool` relevant if `parallel=true`. If it is explicitly given it should be an array of `Threads.nthreads()` preallocated `KnnResult` objects used to reduce memory allocations.
    In most case uses the default is enought, but different pools should be used when indexes use internal indexes to solve queries (e.g., using index's proxies or database objects defined as indexes).

Note: The i-th column in indices and distances correspond to the i-th query in `Q`
Note: The final indices at each column can be `0` if the search process was unable to retrieve `k` neighbors.
"""
function searchbatch(index, Q, k::Integer=10; parallel=false, pools=getpools(index))
    m = length(Q)
    I = zeros(Int32, k, m)
    D = Matrix{Float32}(undef, k, m)
    searchbatch(index, Q, I, D; parallel, pools)
end

"""
    searchbatch(index, Q, I::AbstractMatrix{Int32}, D::AbstractMatrix{Float32}; parallel=false, pools=getpools(index)) -> indices, distances

Searches a batch of queries in the given index and `I` and `D` as output (searches for `k=size(I, 1)`)

"""
function searchbatch(index, Q, I::AbstractMatrix{Int32}, D::AbstractMatrix{Float32}; parallel=false, pools=getpools(index))
    k = size(I, 1)

    if parallel
        Threads.@threads for i in eachindex(Q)
            res, _ = search(index, Q[i], getknnresult(k, pools); pools)
            k_ = length(res)
            I[1:k_, i] .= res.id
            D[1:k_, i] .= res.dist
        end
    else
        @inbounds for i in eachindex(Q)
            res, _ = search(index, Q[i], getknnresult(k, pools); pools)
            k_ = length(res)
            I[1:k_, i] .= res.id
            D[1:k_, i] .= res.dist
        end
    end

    I, D
end

"""
    searchbatch(index, Q, KNN::AbstractVector{KnnResult}; parallel=false, pools=getpools(index)) -> indices, distances

Searches a batch of queries in the given index using an array of KnnResult's; each KnnResult object can specify different `k` values.

"""
function searchbatch(index, Q, KNN::AbstractVector{KnnResult}; parallel=false, pools=getpools(index))
    if parallel
        Threads.@threads for i in eachindex(Q)
            @inbounds search(index, Q[i], KNN[i]; pools)
        end
    else
        @inbounds for i in eachindex(Q)
            search(index, Q[i], KNN[i]; pools)
        end
    end

    KNN
end

function __init__()
    __init__visitedvertices()
    __init__beamsearch()
    __init__neighborhood()

    for _ in 2:Threads.nthreads()
        push!(GlobalKnnResult, KnnResult(32))
    end
end

end

max_stars_repo_name: sadit/NNS

max_issues_repo_name: sadit/NNS

max_forks_repo_name: sadit/NNS


text:
println("loading LowRankModels")
@time @everywhere using LowRankModels

function fit_pca(m,n,k)
	# matrix to encode
	Random.seed!(1)
	A = randn(m,k)*randn(k,n)
	X=randn(k,m)
	Y=randn(k,n)
	losses = fill(QuadLoss(),n)
	r = QuadReg()
	glrm = GLRM(A,losses,r,r,k, X=X, Y=Y)
	glrm = share(glrm)
	p = Params()
	p.max_iter = 10
	X,Y,ch = fit!(glrm)
	println("Convergence history:",ch.objective)
	println("Time/iter:",ch.times[end]/10)
	return A,X,Y,ch
end

@everywhere Random.seed!(1)
fit_pca(100,100,50)


text:
@enum EventKind::UInt8 begin

# the underlying values are also priorites in case the times are equal
# therefore the order of definition implies priority

  QuarantinedEvent
  DetectionEvent
  #DetectedFromQuarantineEvent
  HomeTreatmentSuccessEvent
  QuarantineEndEvent

  GoHospitalEvent

  #DetectedOutsideQuarantineEvent
  #DetectedFromTrackingEvent

  TrackedEvent
  ReleasedEvent

  # the progression events have low priority to let the immediate actions execute
  BecomeInfectiousEvent

  OutsideInfectionEvent
  TransmissionEvent

  MildSymptomsEvent
  HomeTreatmentEvent

  SevereSymptomsEvent
  CriticalSymptomsEvent
  RecoveryEvent
  DeathEvent

  InvalidEvent # should not be executed
end

struct Event
  time::TimePoint                 # 4 bytes
  subject_id::PersonIdx           # 4 bytes
  source_id::PersonIdx            # 4 bytes
  event_kind::EventKind           # 1 byte
  extra::UInt8

  #contact_kind::ContactKind       # 1 byte
  #extension::Bool                 # 1 byte
  #detection_kind::DetectionKind   # 1 byte
                                  # alignment = 4 bytes

  Event() = new(0.0, 0, 0, InvalidEvent, 0)
  Event(::Val{E}, time::Real, subject::Integer) where E = new(time, subject, 0, E, 0)
  Event(::Val{OutsideInfectionEvent}, time::Real, subject::Integer) = new(time, subject, 0, OutsideInfectionEvent, UInt8(OutsideContact))
  Event(::Val{TransmissionEvent}, ::Real, ::Integer) = error("source and contact kind needed for transmission event")
  Event(::Val{TransmissionEvent}, time::Real, subject::Integer, source::Integer, contact_kind::ContactKind) = new(time, subject, source, TransmissionEvent, UInt8(contact_kind))
  Event(::Val{QuarantinedEvent}, time::Real, subject::Integer, extension::Bool) = new(time, subject, 0, QuarantinedEvent, UInt8(extension))
  Event(::Val{TrackedEvent}, ::Real, ::Integer) = error("source and tracking kind must be given for TrackedEvent")
  Event(::Val{TrackedEvent}, time::Real, subject::Integer, source::Integer, tracking_kind::TrackingKind) = new(time, subject, source, TrackedEvent, UInt8(tracking_kind))
  Event(::Val{DetectionEvent}, ::Real, ::Integer) = error("detection kind must be given for detection event")
  Event(::Val{DetectionEvent}, time::Real, subject::Integer, detectionkind::DetectionKind) = new(time, subject, 0, DetectionEvent, UInt8(detectionkind))
end

time(event::Event) = event.time
subject(event::Event) = event.subject_id
source(event::Event) = event.source_id
kind(event::Event) = event.event_kind

contactkind(event::Event) = istransmission(event) ? ContactKind(event.extra) : NoContact
detectionkind(event::Event) = isdetection(event) ? DetectionKind(event.extra) : NoDetection
extension(event::Event) = isquarantine(event) ? Bool(event.extra) : false
trackingkind(event::Event) = istracking(event) ? TrackingKind(event.extra) : NotTracked

function show(io::IO, event::Event)
  print(io, time(event), ":", kind(event), " ", subject(event))
  if TransmissionEvent == kind(event) || OutsideContact == kind(event)
    print(io, " <= ", source(event), " ", contactkind(event))
  elseif QuarantinedEvent == kind(event)
    print(io, " extension=", extension(event))
  elseif TrackedEvent == kind(event)
    print(io, " <= ", source(event))
  elseif InvalidEvent == kind(event)
    print(io, " <= ", source(event), " ", kind(event), " ", event.extra)
  end
end

isdetection(ek::EventKind) = ek == DetectionEvent
isdetection(e::Event) = isdetection(kind(e))
istransmission(ek::EventKind) = ek == TransmissionEvent || ek == OutsideInfectionEvent
istransmission(e::Event) = istransmission(kind(e))
isquarantine(ek::EventKind) = ek == QuarantinedEvent || ek == QuarantineEndEvent
isquarantine(e::Event) = isquarantine(kind(e))
istracking(ek::EventKind) = ek == TrackedEvent
istracking(e::Event) = istracking(kind(e))
ishospitalization(ek::EventKind) = ek == GoHospitalEvent
ishospitalization(e::Event) = ishospitalization(kind(e))
isdeath(ek::EventKind) = ek == DeathEvent
isdeath(e::Event) = isdeath(kind(e))


text:
using DeepQLearning
using POMDPModels
using POMDPSimulators
using POMDPPolicies
using POMDPs
using POMDPModelTools
using Flux
using Random
using StaticArrays
using Test
import CommonRLInterface

RL = CommonRLInterface

Random.seed!(7)
GLOBAL_RNG = MersenneTwister(1) # for test consistency

function evaluate(mdp::Union{MDP,POMDP}, policy, rng, n_ep=100, max_steps=100)
    avg_r = 0.
    sim = RolloutSimulator(rng=rng, max_steps=max_steps)
    for i=1:n_ep
        DeepQLearning.resetstate!(policy)
        avg_r += simulate(sim, mdp, policy)
    end
    return avg_r/=n_ep
end

function evaluate(env::RL.AbstractEnv, policy, rng, n_ep=100, max_steps=100)
    avg_r = 0.
    for i=1:n_ep
        DeepQLearning.resetstate!(policy)
        r = 0.0
        step = 0
        RL.reset!(env)
        while !RL.terminated(env) && step < max_steps
            a = action(policy, RL.observe(env))
            r += RL.act!(env, a)
            step += 1
        end
        avg_r += r
    end
    return avg_r/=n_ep
end

# @testset "TigerPOMDP DDRQN" begin
    pomdp = TigerPOMDP(0.01, -1.0, 0.1, 0.8, 0.95);
    input_dims = reduce(*, size(convert_o(Vector{Float64}, first(observations(pomdp)), pomdp))) + reduce(*, size(convert_a(Vector{Float64}, first(actions(pomdp)), pomdp)))
    @show input_dims
    model = Chain(x->flattenbatch(x), LSTM(input_dims, 4), Dense(4, length(actions(pomdp))))
    max_steps = 10000
    exploration = EpsGreedyPolicy(pomdp, LinearDecaySchedule(start=1.0, stop=0.01, steps=max_steps/2),
                                  rng=GLOBAL_RNG)
    solver = DeepQLearningSolver(qnetwork = model, prioritized_replay=false, max_steps=max_steps,
                             learning_rate=0.0001, exploration_policy = exploration,
                             log_freq=500, target_update_freq = 1000,
                             recurrence=true,trace_length=10, double_q=true, dueling=true, max_episode_length=100)

    policy = solve(solver, pomdp)
    @test size(actionvalues(policy, true)) == (length(actions(pomdp)),)
# end


text:
module VortexDistributions

using Test
using JLD2
using Parameters
using SpecialFunctions
using Interpolations
using LinearAlgebra
using ToeplitzMatrices
using SparseArrays
using FFTW
using FileIO
using ProgressMeter
using LightGraphs
using SimpleWeightedGraphs
import Plots:stroke,scatter!,plot,plot!

const Λ = 0.8249

# topology
export Field, Torus, Sphere,

# vortex groups
VortexGroup, Dipole, Cluster,

# cores
Vortex, CoreShape, Ansatz, Exact, ScalarVortex, PointVortex,

# detection
findvortices, found_near, phase_jumps, phase_jumps!, unwrap, unwrap!, Δ,
find_where, findvortices_jumps, findvortices_grid,
remove_vortices_edge, zoom_interp, zoom_grid,
circ_mask, keep_vortices,

# construction
scalar_ansatz, vortex_array, uniform,
vortex!, dipole_phase, periodic_dipole!,
rand_charge, rand_pointvortex, rand_scalarvortex, rand_vortex, rand_vortexfield,
thetad,

# convenient access
charge, xpos, ypos, pos,

# RCA
distances, periodic_distances, sparse_distances,
spanning_tree, vortex_marker, plot_vortices!,
plot_cluster!, vortex_marker, grow_plus_clusters,
grow_minus_clusers, get_dipoles, seed_clusters

include("types.jl")
include("pointvortex.jl")
include("detection.jl")
include("creation.jl")

# RCA
include("get_dipoles.jl")
include("seed_clusters.jl")
include("grow_minus_clusters.jl")
include("grow_plus_clusters.jl")
# include("get_spanning_trees.jl")
# include("get_negative_spanning_trees.jl")
# include("get_positive_spanning_trees.jl")
# include("find_smallest_spanning_trees.jl")
# include("sort_RCA_structs.jl")
# include("RCA.jl")

# utils
include("utils.jl")

@load joinpath(@__DIR__,"exactcore.jld2") ψi
@load joinpath(@__DIR__,"ansatzcore.jld2") ψa

end


text:
@testset "Top-level interface" begin
    rosenbrock(x) = sum(i -> 100*abs2(x[i+1] - x[i]^2) + abs2(x[i] - 1), Base.OneTo(length(x)-1))

    # Ensure things have compiled before we start testing:
    bboptimize(rosenbrock; SearchRange = (-5.0, 5.0), NumDimensions = 2,
        MaxSteps = 10, TraceMode = :silent)

    @testset "run a simple optimization" begin
        @testset "using bboptimize() with mostly defaults" begin
            res = bboptimize(rosenbrock; SearchRange = (-5.0, 5.0), NumDimensions = 2,
                MaxSteps = 2000, TraceMode = :silent)
            @test best_fitness(res) < 0.1
            xbest = best_candidate(res)
            @test typeof(xbest) == Vector{Float64}
            @test length(xbest) == 2

            # We also mimic some of the Optim.jl api (although it is a bit strange...)
            @test f_minimum(res) < 0.1
            @test minimum(res) == xbest
            @test iteration_converged(res)
        end

        @testset "using bbsetup()/bboptimize() with mostly defaults" begin
            opt = bbsetup(rosenbrock; SearchRange = (-5.0, 5.0), NumDimensions = 2,
                MaxSteps = 2000, TraceMode = :silent)
            @test numruns(opt) == 0
            @test isa(problem(opt), BlackBoxOptim.FunctionBasedProblem)
            res = bboptimize(opt)
            @test numruns(opt) == 1
            @test isa(lastrun(opt), BlackBoxOptim.OptRunController)

            @test best_fitness(res) < 0.1
            xbest = best_candidate(res)
        end

        @testset "using non-population optimizer" begin
            res = bboptimize(rosenbrock; Method=:generating_set_search,
                                            SearchRange = (-5.0, 5.0), NumDimensions = 2,
                                            MaxSteps = 5000, TraceMode = :silent)
            @test best_fitness(res) < 1.0
            xbest = best_candidate(res)
            @test typeof(xbest) == Vector{Float64}
        end

        @testset "using population optimizer" begin
            res = bboptimize(rosenbrock; Method=:adaptive_de_rand_1_bin,
                                            SearchRange = (-5.0, 5.0), NumDimensions = 2,
                                            MaxSteps = 5000, TraceMode = :silent)
            @test best_fitness(res) < 0.1
            xbest = best_candidate(res)
            @test typeof(xbest) == Vector{Float64}
            xpop = population(res)
            @test isa(xpop, BlackBoxOptim.Population)
            @test popsize(xpop) > 0
            @test numdims(xpop) == 2

            # Access a few individual solution vectors in the population...
            @test isa(xpop[1], Array{Float64, 1}) == true # first solution vector
            @test isa(xpop[popsize(xpop)], Array{Float64, 1}) == true # last solution vector
            rand_solution_idx = rand(1:popsize(xpop))
            @test isa(xpop[rand_solution_idx], Array{Float64, 1}) == true # random solution vector

            # and to access their fitness values:
            @test isa(fitness(xpop, 1), Float64)
            @test isa(fitness(xpop, popsize(xpop)), Float64)
            @test isa(fitness(xpop, rand_solution_idx), Float64)

            # Ensure the lowest fitness value is the one returned by best_fitness
            min_fitness_value = minimum(map(i -> fitness(xpop, i), 1:popsize(xpop)))
            @test (min_fitness_value == best_fitness(res))
        end

        @testset "using population optimizer and parallel evaluator" begin
            opt = bbsetup(rosenbrock; Method=:adaptive_de_rand_1_bin,
                                        SearchRange = (-5.0, 5.0), NumDimensions = 2,
                                        MaxSteps = 2000, TraceMode = :silent, Workers=workers())
            res = bboptimize(opt)
            @test isa(BlackBoxOptim.evaluator(lastrun(opt)), BlackBoxOptim.ParallelEvaluator)
        end

        @testset "using population optimizer and multithread evaluator" begin
            if Threads.nthreads() > 1
                opt = bbsetup(rosenbrock; Method=:adaptive_de_rand_1_bin,
                                            SearchRange = (-5.0, 5.0), NumDimensions = 2,
                                            MaxSteps = 2000, TraceMode = :silent, NThreads=Threads.nthreads()-1)
                res = bboptimize(opt)
                @test isa(BlackBoxOptim.evaluator(lastrun(opt)), BlackBoxOptim.MultithreadEvaluator)
            else
                @warn "Top-level MultithreadEvaluator tests require >1 threads, $(Threads.nthreads()) found, use JULIA_NUM_THREADS"
                @test_skip Threads.nthreads() > 1
            end
        end
    end

    @testset "continue running an optimization after it finished" begin
        optctrl = bbsetup(rosenbrock; SearchRange = (-5.0, 5.0), NumDimensions = 100,
            MaxTime = 0.5, TraceMode = :silent)

        res1 = bboptimize(optctrl)
        @test numruns(optctrl) == 1

        res2 = bboptimize(optctrl; MaxTime = 1.0)
        @test numruns(optctrl) == 2

        @test best_fitness(res2) <= best_fitness(res1)

        # parameters should be the same except for MaxTime
        for p in keys(flatten(parameters(res1)))
            if p != :MaxTime
                @test parameters(res1)[p] == parameters(res2)[p]
            end
        end
        @test parameters(res1)[:MaxTime] == 0.5
        @test parameters(res2)[:MaxTime] == 1.0
    end

    @testset "fitness decrease monotonically if optimizing with same optctrl repeatedly" begin
        # There was a bug in 0.3.0 that did not give monotonically decreasing
        # fitness with repeated runs. The bug only happened for very short MaxSteps
        # lengths, such as 10. This test ensures it does not resurface.
        optctrl = bbsetup(rosenbrock; SearchRange = (-5.0, 5.0), NumDimensions = 100,
            MaxSteps = 10, TraceMode = :silent)
        optresults = [bboptimize(optctrl) for _ in 1:10]
        bestfitnesses = map(best_fitness, optresults)
        for i in 2:length(optresults)
            @test bestfitnesses[i-1] >= bestfitnesses[i]
        end
    end

    @testset "return results after interruption" begin
        i = 0
        function rosenbrock_throwing(x)
            i += 1
            if i < 50
                return sum(i -> 100*(x[i+1] - x[i])^2 + (x[i] - 1)^2, 1:(length(x)-1))
            else
                throw(InterruptException())
            end
        end
        @testset ":RecoverResults on" begin
            i = 0
            optctrl = bbsetup(rosenbrock_throwing; SearchRange = (-5.0, 5.0), NumDimensions = 100,
                              MaxSteps=100, TraceMode=:silent, RecoverResults=true)
            res = bboptimize(optctrl)
            @test BlackBoxOptim.isinterrupted(res)
        end

        @testset ":RecoverResults off" begin
            i = 0 # reset the counter, otherwise it will throw in the setup
            optctrl = bbsetup(rosenbrock_throwing; SearchRange = (-5.0, 5.0), NumDimensions = 100,
                              MaxSteps=100, TraceMode=:silent, RecoverResults=false)
            @test_throws InterruptException bboptimize(optctrl)
        end
    end

    @testset "continue running an optimization after serializing to disc" begin
        using Serialization: serialize, deserialize

        optctrl = bbsetup(rosenbrock; SearchRange = (-5.0, 5.0), NumDimensions = 100,
                          MaxTime = 0.5, TraceMode = :silent)
        res1 = bboptimize(optctrl)

        local tempfilename = "./temp" * string(rand(1:10^8)) * ".tmp"

        try # To ensure we delete the temp file afterwards...
            open(tempfilename, "w") do fh
                serialize(fh, optctrl)
            end

            # Try to make sure its not in mem:
            opctrl = nothing; GC.gc()

            local ocloaded
            open(tempfilename, "r") do fh
                ocloaded = deserialize(fh)
            end

            @test numruns(ocloaded) == 1
            res2 = bboptimize(ocloaded; MaxTime = 1.0)
            @test numruns(ocloaded) == 2

            @test best_fitness(res2) <= best_fitness(res1)
        finally
            if isfile(tempfilename)
                rm(tempfilename)
            end
        end
    end

    @testset "TargetFitness option works" begin
        # FIXME use the same (fixed?) random seed to guarantee reproducibility
        result1 = bboptimize(rosenbrock, SearchRange = (-5.0, 5.0), NumDimensions = 5,
                             Method = :de_rand_1_bin, FitnessTolerance = 1e-5,
                             MaxSteps = 1000000, TraceMode = :silent,
                             TargetFitness = 0.0)
        result2 = bboptimize(rosenbrock, SearchRange = (-5.0, 5.0), NumDimensions = 5,
                             Method = :de_rand_1_bin, FitnessTolerance = 1e-5,
                             MaxSteps = 1000000, TraceMode = :silent)
        @test best_fitness(result1) < 1e-5
        @test result1.iterations < result2.iterations
    end

    @testset "callback" begin
        global callbacktimes = Float64[]
        function callbackfn(optctrl)
            global callbacktimes
            push!(callbacktimes, time())
        end
        interval = 0.005
        NumCalls = 30
        res = bboptimize(rosenbrock; SearchRange = (-5.0, 5.0), NumDimensions = 20,
            TraceMode = :silent,
            MaxTime = 1.0 + interval * NumCalls, # Some extra time for startup/compilation etc
            CallbackInterval = interval, CallbackFunction = callbackfn)
        @test length(callbacktimes) >= NumCalls
        for i in 1:(length(callbacktimes) - 1)
            d = (callbacktimes[i+1] - callbacktimes[i])
            @test d >= interval
        end

        # If we call with an interval which is negative there are no callbacks
        prenumcalls = length(callbacktimes)
        res = bboptimize(rosenbrock; SearchRange = (-5.0, 5.0), NumDimensions = 20,
            TraceMode = :silent,
            MaxTime = 1.0 + interval * NumCalls, # Some extra time for startup/compilation etc
            CallbackInterval = -2.0, CallbackFunction = callbackfn)
        @test length(callbacktimes) == prenumcalls
    end

    @testset "callable as objective" begin
        struct CallableObjective1 end
        (::CallableObjective1)(x) = sum(abs2, x)
        res = bboptimize(CallableObjective1(); SearchRange = [(-2.0, 3.0), (-4.0, 5.0)])
        @test best_candidate(res) ≈ zeros(2) atol = 1e-5
    end

    @testset "explicit termination via shutdown method" begin
        cb = function (oc)
            if best_fitness(oc) < 50.0
                BlackBoxOptim.shutdown_optimizer!(oc)
            end
        end
        res = bboptimize(rosenbrock; SearchRange = (-5.0, 5.0), NumDimensions = 2,
                PopulationSize = 10, TraceMode = :silent, MaxFuncEvals = 1000,
                CallbackFunction = cb, CallbackInterval = 0.0)
        reason = BlackBoxOptim.stop_reason(res)
        @test occursin("stopped", reason) && occursin("shutdown method", reason)
    end
end


text:
module TestMaybe

using Test

# TODO: move this to Compat.jl
if try
    >=(1)
    false
catch
    true
end
    @info "Defining `>=(_)`"
    Base.:(>=)(b) = a -> a >= b
end

@testset "$file" for file in sort([
    file for file in readdir(@__DIR__) if match(r"^test_.*\.jl$", file) !== nothing
])
    if file == "test_doctest.jl"
        if lowercase(get(ENV, "JULIA_PKGEVAL", "false")) == "true"
            @info "Skipping doctests on PkgEval."
            continue
        elseif VERSION >= v"1.7-"
            @info "Skipping doctests on Julia $VERSION."
            continue
        end
    end

    include(file)
end

end  # module


text:
using PyCall

println("Running build.jl for the Surprise package.")

# Change that to whatever packages you need.
const PACKAGES = ["numpy", "scikit-surprise"]

# Use eventual proxy info
proxy_arg=String[]
if haskey(ENV, "http_proxy")
    push!(proxy_arg, "--proxy")
    push!(proxy_arg, ENV["http_proxy"])
end

# Import pip
try
    @pyimport pip
catch ee
    typeof(ee) <: PyCall.PyError || rethrow(ee)
    error("""
Python Pip not installed
Please either:
- Install pip
- Rebuild Surprise.jl via `Pkg.build("Surprise")` in the julia REPL
""")
end

try
    @pyimport surprise
catch
    println("Installing required python packages using pip")
    run(`$(PyCall.python) $(proxy_arg) -m pip install --user $(PACKAGES)`)
end


text:
# Define the theoretic model
################### Assemble matrix and vector #####################
# Weak form of the Helmholtz equation :
# Hz: a(p,u,v)=Λ⋅∇v⋅ξ(p)Λ⋅∇u-k²μv⋅u
# Material distribution
#ξ0(x,ϵ1,ϵ2) = x[2]>0 ? 1/ϵ1 : 1/ϵ2
#ξd(p,ϵmin,ϵmax)= 1/(ϵmin + (ϵmax-ϵmin)*p) - 1/ϵmin # in the design region
#a0(u,v,ϵ1,ϵ2,μ,σs,k,LHp,LHn,dpml) = ((x->Λ(x,σs,k,LHp,LHn,dpml))⋅∇(v))⊙((x->ξ0(x,ϵ1,ϵ2))*((x->Λ(x,σs,k,LHp,LHn,dpml))⋅∇(u))) - k^2*μ*(v*u)
#a_d(u,v,ph,k,ϵ1,ϵd) = ∇(v)⊙(((p -> ξd(p, ϵ1, ϵd))∘ph)*∇(u))
# Ez: a(p,u,v)=1/μ*Λ⋅∇v⋅Λ⋅∇u-k²ξ(p)uv
ξd(p,ϵmin,ϵmax)= ϵmin+(ϵmax-ϵmin)*p # in the design region
a0(u,v,ph,ϵ1,ϵ2,μ,σs,k,LHp,LHn,Hd,dpml) = ((x->Λ(x,σs,Hd,k,LHp,LHn,dpml))⋅∇(v))⊙((1/μ)*((x->Λ(x,σs,Hd,k,LHp,LHn,dpml))⋅∇(u))) - k^2*(((p -> ξd(p, ϵ1, ϵ2))∘ph)*(v*u))
#a_d(u,v,ph,k,ϵ1,ϵd) = -k^2*((p -> ξd(p, ϵ1, ϵd))∘ph)*(v*u)
# Source term (Gaussian point source approximation at center)

#b(v,x0,δ) = v*(x->GaussianD(x,x0,δ))

function MatrixA(ph,ϵ1,ϵ2,μ,σs,k,LHp,LHn,dpml,Hd,U,V,dΩ)
    # Assemble the matrix
    return assemble_matrix(U,V) do u,v
        ∫( a0(u,v,ph,ϵ1,ϵ2,μ,σs,k,LHp,LHn,Hd,dpml) )dΩ
    end
end

function MatrixB(x0,δ,Amp,V,dΩ,dΓ)
    l_temp(v) = ∫( Amp*v*x->GaussianY(x,x0,δ) )*dΩ+∫( 0*v )*dΓ
    #op = AffineFEOperator((u,v)->(∫(u*v)*dΩ),l_temp,U,V)
    #op = FESource(l_temp,U,V)
    assemble_vector(l_temp,V)
end

max_stars_repo_name: Tiamat-Tech/hpinn

max_issues_repo_name: louhz/hpinn

max_forks_repo_name: louhz/hpinn


text:
module ExactConversions

export mincommon, maxcommon, exactconv

"""
    mincommon(I, F) :: I
    mincommon(F, I) :: F

Return the smallest number that is representable by both the integral type `I` and the
floating-point type `F`. The type of the result depends on the order of the arguments.
"""
function mincommon end

"""
    maxcommon(I, F) :: I
    maxcommon(F, I) :: F

Return the largest number that is representable by both the integral type `I` and the
floating-point type `F`. The type of the result depends on the order of the arguments.
"""
function maxcommon end

# unsigned integral types

mincommon(U::Type{<:Unsigned}, F::Type{<:AbstractFloat})::U = zero(U)

mincommon(F::Type{<:AbstractFloat}, U::Type{<:Unsigned})::F = F(mincommon(U, F))

function maxcommon(U::Type{<:Unsigned}, F::Type{<:AbstractFloat})::U
    max_exp = min(8sizeof(U)-1, Base.exponent_max(F)-1)
    num_ones = min(8sizeof(U), Base.significand_bits(F)+1)
    ones = (U(1) << num_ones) - U(1)
    ones << (max_exp - num_ones + 1)
end

maxcommon(F::Type{<:AbstractFloat}, U::Type{<:Unsigned})::F = F(maxcommon(U, F))

# signed integral types

function mincommon(S::Type{<:Signed}, F::Type{<:AbstractFloat})::S
    if 8sizeof(S)-1 <= Base.exponent_max(F)-1
        S(1) << (8sizeof(S)-1)
    else
        -maxcommon(S, F)
    end
end

mincommon(F::Type{<:AbstractFloat}, S::Type{<:Signed})::F = F(mincommon(S, F))

function maxcommon(S::Type{<:Signed}, F::Type{<:AbstractFloat})::S
    max_exp = min(8sizeof(S)-2, Base.exponent_max(F)-1)
    num_ones = min(8sizeof(S)-1, Base.significand_bits(F)+1)
    ones = (S(1) << num_ones) - S(1)
    ones << (max_exp - num_ones + 1)
end

maxcommon(F::Type{<:AbstractFloat}, S::Type{<:Signed})::F = F(maxcommon(S, F))

"""
    exactconv(I, x::F) :: Union{I,Nothing}
    exactconv(F, x::I) :: Union{F,Nothing}

Convert between integer and floating-point types. The conversion suceeds only if
it does not change the numeric value, otherwise it returns `nothing`.
"""
function exactconv end

function exactconv(I::Type{<:Union{Signed,Unsigned}}, x::F)::Union{I,Nothing} where F<:AbstractFloat
    if isinteger(x) && mincommon(F, I) <= x <= maxcommon(F, I)
        I(x)
    else
        nothing
    end
end

function exactconv(F::Type{<:AbstractFloat}, x::U)::Union{F,Nothing} where U<:Unsigned
    bw = 8sizeof(U) - leading_zeros(x) - trailing_zeros(x)
    if bw <= Base.significand_bits(F)+1 && mincommon(U, F) <= x <= maxcommon(U, F)
        F(x)
    else
        nothing
    end
end

function exactconv(F::Type{<:AbstractFloat}, x::S)::Union{F,Nothing} where S<:Signed
    if x >= S(0)
        bw = 8sizeof(S) - leading_zeros(x) - trailing_zeros(x)
        if bw <= Base.significand_bits(F)+1 && mincommon(S, F) <= x <= maxcommon(S, F)
            F(x)
        else
            nothing
        end
    elseif x == ~S(0) && x == mincommon(S, F)
        F(x)
    else
        bw = 8sizeof(S) - leading_zeros(-x) - trailing_zeros(-x)
        if bw <= Base.significand_bits(F)+1 && mincommon(S, F) <= x <= maxcommon(S, F)
            F(x)
        else
            nothing
        end
    end
end

end


text:
# SPDX-License-Identifier: MIT

# Helper functions for plotting 2D-Normal pdfs.
# Lucas Ondel, 2021

function normalShape(μ, Σ, σ = 1)
    θ = LinRange(0, 2π, 1000)
    λ, U = eigen(Σ)
    B = U * diagm(σ * sqrt.(λ))
    circle = hcat(sin.(θ), cos.(θ))'
    contour = B * circle .+ μ
    contour[1, :], contour[2, :]
end

function plotnormal!(μ::AbstractVector, Σ::AbstractMatrix; σ = 1, kwargs...)
    plot!(normalShape(μ, Σ, σ), seriestype = [:shape];  kwargs...)
end

max_stars_repo_name: BUTSpeechFIT/PPCA

max_issues_repo_name: BUTSpeechFIT/PPCA

max_forks_repo_name: BUTSpeechFIT/PPCA


text:
function build_op_problem!(op_problem::OperationsProblem{M}; kwargs...) where M<:AbstractOperationsProblem
    sys = get_system(op_problem)
    _build!(op_problem.canonical, op_problem.template, sys; kwargs...)
    return
end

function _build!(canonical::Canonical, ref::OperationsTemplate, sys::PSY.System; kwargs...)

    verbose = get(kwargs, :verbose, true)
    transmission = ref.transmission

    # Build Injection devices
    for device_model in values(ref.devices)
        verbose && @info "Building $(device_model.device_type) with $(device_model.formulation) formulation"
        construct_device!(canonical, sys, device_model, transmission; kwargs...)
    end

    # Build Network
    verbose && @info "Building $(transmission) network formulation"
    construct_network!(canonical, sys, transmission; kwargs...)

    # Build Branches
    for branch_model in values(ref.branches)
        verbose && @info "Building $(branch_model.device_type) with $(branch_model.formulation) formulation"
        construct_device!(canonical, sys, branch_model, transmission; kwargs...)
    end

    #Build Service
    for service_model in values(ref.services)
        #construct_service!(canonical, sys, service_model, transmission; kwargs...)
    end

    # Objective Function
    verbose && @info "Building Objective"
    JuMP.@objective(canonical.JuMPmodel,
                    MOI.MIN_SENSE,
                    canonical.cost_function)

    return

end


text:
module SWqmmm

using JuLIP, QMMM2, NeighbourLists, LinearAlgebra, SHIPs, JuLIPMaterials
using JuLIP.MLIPs: IPSuperBasis

pot() = StillingerWeber()
at_train() = bulk(:Si, cubic=true) * 2

function sketch()
   sw =  pot()
   rcut = cutoff(sw)
   at = at_train()
   return QMMM2.data_djEs_sketch(at, 1e-3, [rcut, rcut, rcut])
end

function training_set()
   sw = pot()
   sw_sketch = sketch()
   QMMM2.eval_dataset!(sw_sketch, sw; key="sw")
end

function get_basis(bo, deg; rinfact = 0.8, wL = 1.5)
   sw = pot()
   r0 = rnn(:Si)
   rcut = cutoff(sw)
   rin = rinfact * r0
   trans = PolyTransform(2, rnn(:Si))
   fcut = PolyCutoff1s(2, rin, rcut)
   spec = SparseSHIP(:Si, bo, deg; wL=wL)
   shipB = SHIPBasis(spec, trans, fcut)
   basis = IPSuperBasis(OneBody(1.0), shipB)
end

function train_ship(bo, deg;
            wL = 1.5,
            rinfact = 0.8,
            weights = Dict("Es" => 100.0, "dEs" => 10.0, "d2Esh" => 1.0),
            rtol = 1e-5
         )
   trainset = training_set()
   basis = get_basis(bo, deg; rinfact=rinfact, wL=wL)
   🚢, fitinfo = QMMM2.lsqfit(basis, trainset, weights; rtol=rtol, key="sw")
   @show fitinfo["rmse"]
   return 🚢, fitinfo
end

function int_config(Rdom, calc = StillingerWeber(), Rbuf = 2 * cutoff(calc))
   # -----------------------------------
   #  construct the configuration
   at = cluster(:Si, Rdom+Rbuf)
   X = positions(at)
   i0 = 1
   x0 = X[i0]
   r = [ norm(x - x0) for x in X]
   Ineig = findall(r .< 3.0)
   Rneig = [ x - x0 for x in X[Ineig] ]
   rup = [1.3575, 1.3575, 1.3575]
   if rup ∈ Rneig
      σ = 1.0
   else
      σ = -1.0
   end
   rleft = σ * [-1.3575, -1.3575, 1.3575]
   dint = 0.5 * (norm(rup) - dot(rleft, rup/norm(rup)))
   xint = x0 - σ * rup * dint
   at = append(at, [xint])
   # hack at.Z => TODO in JuLIP
   empty!(at.Z)
   append!(at.Z, fill(atomic_number(:Si), length(at)) )
   # -----------------------------------
   #  compute the domains
   X = positions(at)
   r = [ norm(x - x0) for x in X]
   set_free!(at, findall(r .<= Rdom))
   set_data!(at, "xcore", x0)
   set_calculator!(at, calc)
   return at
end

sw_eq = JuLIPMaterials.Si.sw_eq()

function edge_config(Rdom, calc = sw_eq, Rbuf = 2 * cutoff(calc))
   # -----------------------------------
   #  construct the configuration
   at, x0 = JuLIPMaterials.Si.edge110(:Si, Rdom + Rbuf;
                                      calc=calc, truncate=true, sym=true)
   set_data!(at, "xcore", x0)
   # -----------------------------------
   #  compute the domains
   X = positions(at)
   r = [ norm(x[1:2] - x0[1:2]) for x in X ]
   Ifree = findall(r .<= Rdom)
   set_free!(at, Ifree)
   set_calculator!(at, calc)
   return at, Ifree
end

end

# # test convergence
# for deg = 2:6
#    🚢, fitinfo = SWqmmm.train_ship(3, deg)
# end