text:
"""
    Fermi.Integrals

Module to compute integrals using Lints.jl
"""
module Integrals

using Fermi
using Fermi.Output
using Fermi.Geometry: Molecule
using Fermi.Orbitals: AbstractOrbitals, OrbDict
using Lints
using LinearAlgebra
using TensorOperations

import Base.getindex
import Base.setindex!

"""
    IntegralHelper{T}

Structure to assist with computing and storing integrals.
Accesss like a dictionary e.g.,
    ints["B"]
to get DF ERI integrals. A number of characters can be added to denote orbitals in various
bases, such as MO or NO.
     O               -> occ,alpha,ERI
     o               -> occ,beta,ERI
     V               -> vir,alpha,ERI
     v               -> vir,beta,ERI
     P               -> all,alpha,ERI
     p               -> all,beta,ERI
     B               -> DF-ERI
     μ               -> ao,ERI
     Ω               -> NO,occ,alpha
     ω               -> NO,occ,beta
     U               -> NO,vir,alpha
     u               -> NO,vir,beta
     S               -> AO,overlap
     T               -> AO,kinetic
     V               -> AO,nuclear

# Fields
    cache::Dict{String,Array}                        previously computed integrals
    bname::Dict{String,String}                       basis set names for various purposes
    mol::Molecule                                    attached Molecule object.
    orbs::Dict{String,O} where O < AbstractOrbitals  orbitals of various kinds
    basis::Dict{String,Lints.BasisSetAllocated}      basis sets for various purposes
"""
mutable struct IntegralHelper{T}
    cache::Dict{String,Array}
    bname::Dict{String,String}
    mol::Molecule
    orbs::Fermi.Orbitals.OrbDict
    basis::Dict{String,Lints.BasisSetAllocated}
    type::DataType
    normalize::Bool
end

function IntegralHelper(;mol=Molecule(),orbs=Fermi.Orbitals.OrbDict())
    IntegralHelper{Float64}(;mol=mol,orbs=orbs)
end

function IntegralHelper{T}(;mol=Molecule(),orbs=Fermi.Orbitals.OrbDict()) where T <: AbstractFloat
    cache = Dict{String,Array}()
    bname = Dict{String,String}()
    type = T
    bname["primary"] = Fermi.CurrentOptions["basis"]
    aux = Fermi.CurrentOptions["jkfit"]
    if aux == "auto"
        aux_lookup = Dict{String,String}(
                                         "cc-pvdz" => "cc-pvdz-jkfit",
                                         "cc-pvtz" => "cc-pvtz-jkfit",
                                         "cc-pvqz" => "cc-pvqz-jkfit",
                                         "cc-pv5z" => "cc-pv5z-jkfit"
                                        )
        bname["aux"] = try
            aux_lookup[Fermi.CurrentOptions["basis"]]
        catch KeyError #if we haven't got it programmed, use a large DF basis by default
            "aug-cc-pvqz-rifit"
        end
    else
        bname["aux"] = aux
    end
    #mol = Molecule()
    #orbs = Fermi.Orbitals.OrbDict()
    basis = Dict{String,Lints.BasisSetAllocated}()
    IntegralHelper{T}(cache,bname,mol,orbs,basis,type,false)
end

function normalize!(I::IntegralHelper,normalize::Bool)
    I.normalize=normalize
    for entry in keys(I.cache)
        delete!(I.cache,entry)
    end
end

"""
    aux_ri!(I::IntegralHelper, ri=Fermi.CurrentOptions["rifit"])

Clears the integral cache and switches auxiliary DF integrals to use the
current RI fitting basis set. Used between DF-RHF and DF-post HF.
"""
function aux_ri!(I::IntegralHelper,ri=Fermi.CurrentOptions["rifit"])
    delete!(I.cache,"B") #clear out old aux basis
    GC.gc()
    if ri == "auto"
        aux_lookup = Dict{String,String}(
                                         "cc-pvdz" => "cc-pvdz-rifit",
                                         "cc-pvtz" => "cc-pvtz-rifit",
                                         "cc-pvqz" => "cc-pvqz-rifit",
                                         "cc-pv5z" => "cc-pv5z-rifit"
                                        )
        I.bname["aux"] = try
            aux_lookup[Fermi.CurrentOptions["basis"]]
        catch KeyError
            "aug-cc-pvqz-rifit" # default to large DF basis
        end
    else
        I.bname["aux"] = ri
    end
end
"""
    aux_ri!(I::IntegralHelper, jk=Fermi.CurrentOptions["jkfit"])

Clears the integral cache and switches auxiliary DF integrals to use the
current JK fitting basis set. Used to ensure JK integrals are used in
DF-RHF.
"""
function aux_jk!(I::IntegralHelper,jk=Fermi.CurrentOptions["jkfit"])
    delete!(I.cache,"B") #clear out old aux basis
    GC.gc()
    if jk == "auto"
        aux_lookup = Dict{String,String}(
                                         "cc-pvdz" => "cc-pvdz-jkfit",
                                         "cc-pvtz" => "cc-pvtz-jkfit",
                                         "cc-pvqz" => "cc-pvqz-jkfit",
                                         "cc-pv5z" => "cc-pv5z-jkfit"
                                        )
        I.bname["aux"] = try
            aux_lookup[Fermi.CurrentOptions["basis"]]
        catch KeyError
            "aug-cc-pvqz-rifit" # default to large DF basis
        end
    else
        I.bname["aux"] = jk
    end
end

"""
    getindex(I::IntegralHelper,entry::String)

Called when `helper["foo"]` syntax is used. If the requested entry already
exists, simply return the entry. If not, compute the requested entry.
"""
function getindex(I::IntegralHelper,entry::String)
    try
        I.cache[entry]
    catch KeyError
        compute!(I,entry)
        I.cache[entry]
    end
end

function setindex!(I::IntegralHelper,val,entry)
    I.cache[entry] = val
end

function compute!(I::IntegralHelper,entry::String)
    o = I.orbs
    if entry == "μ" #AO basis eri. 4 index
        I.cache["μ"],_ = aoeri(I.mol,I.bname["primary"],normalize=I.normalize)

    elseif entry == "B" #AO basis eri. 3 index. DF
        I.cache["B"],I.cache["J"] = dfaoeri(I.mol,I.bname["primary"],I.bname["aux"],normalize=I.normalize)

    elseif entry == "S" #AO basis overlap
        I.cache["S"],_ = aooverlap(I.mol,I.bname["primary"],normalize=I.normalize)

    elseif entry == "T" #AO basis kinetic
        I.cache["T"],_ = aokinetic(I.mol,I.bname["primary"],normalize=I.normalize)

    elseif entry == "V" #AO basis nuclear
        I.cache["V"],_ = aonuclear(I.mol,I.bname["primary"],normalize=I.normalize)

    elseif entry == "F"
        D = Fermi.contract(o["[O]"],o["[O]"],"um","vm")
        F = zeros(I.type,size(I["S"]))
        Fermi.HartreeFock.build_fock!(F,
                                                     I["T"] + I["V"],
                                                     D,
                                                     I["μ"],
                                                     o["[O]"])
        I.cache["F"] = F
    elseif entry == "F'"
        D = Fermi.contract(o["[O]"],o["[O]"],"um","vm")
        F = zeros(I.type,size(I["S"]))
        Fermi.HartreeFock.build_fock!(F,
                                                     I["T"] + I["V"],
                                                     D,
                                                     I["B"],
                                                     o["[O]"])
        I.cache["F'"] = F
    elseif 'B' in entry #MO basis eri. 3 index. DF
        aoint = I["B"]
        C1 = try
            o[entry[2:2]]
        catch KeyError
            error("")
        end
        C2 = try
            o[entry[3:3]]
        catch KeyError
            error("")
        end
        C = []
        push!(C,C1)
        push!(C,C2)
        _Pqp = transform_eri(aoint,C...)
        Bqp = Array{Float64}(undef,size(_Pqp))
        Fermi.contract!(Bqp,I.cache["J"],_Pqp,"Qqp","PQ","Pqp")
        I.cache[entry] = Bqp

    elseif 'F' in entry
        df = '\'' in entry
        if df
            F =  I["F'"]
        else
            F = I["F"]
        end
        C1 = try
            o[entry[2:2]]
        catch
            error("")
        end
        C2 = try
            o[entry[3:3]]
        catch
            error("")
        end
        C = []
        push!(C,C1)
        push!(C,C2)
        I[entry] = transform_fock(F,C...)

    else
        aoint  = I["μ"]
        if ('<' in entry || '(' in entry)
            offset = 1
        else
            offset = 0
        end
        if '(' in entry
            notation = "chem"
        else
            notation = "phys"
        end
        C = []
        C1 = try
            o[entry[1+offset:1+offset]]
        catch
            error("")
        end
        C2 = try
            o[entry[2+offset:2+offset]]
        catch
            error("")
        end
        C3 = try
            o[entry[3+offset:3+offset]]
        catch
            error("")
        end
        C4 = try
            o[entry[4+offset:4+offset]]
        catch
            error("")
        end
        push!(C,C1)
        push!(C,C3)
        push!(C,C2)
        push!(C,C4)
        temp = transform_eri(aoint, C...)
        if notation == "phys"
            I.cache[entry] = permutedims(temp,(1,3,2,4))
        else
            I.cache[entry] = temp
        end
    end
end

"""
    aokinetic(molecule::Fermi.Molecule, basis::String)

Computes AO basis kinetic energy ⟨μ|T̂|ν⟩ integrals for the given basis and molecule.
Can be accessed at a higher level by calling

    helper["T"]

where `helper` is bound to the desired molecule and basis set.
"""
function aokinetic(molecule::Molecule, basis::String; normalize=false)#, interconnect::Fermi.Environments.No_IC,
                                                     #   communicator::Fermi.Environments.NoCommunicator,
                                                     #   accelerator::Fermi.Environments.NoAccelerator)

    mol = Fermi.Geometry.to_lints_molecule(molecule)
    bas = Lints.BasisSet(basis, mol)
    @lints begin
        T = Lints.make_T(bas; normalize=normalize)
    end
    T,bas
end
"""
    aooverlap(molecule::Fermi.Molecule, basis::String)

Computes AO basis overlap ⟨p|q⟩ integrals for the given basis and molecule.
Can be accessed at a higher level by calling

    helper["S"]

where `helper` is bound to the desired molecule and basis set.
"""
function aooverlap(molecule::Molecule, basis::String; normalize=false)#, interconnect::Fermi.Environments.No_IC,
                                                     #   communicator::Fermi.Environments.NoCommunicator,
                                                     #   accelerator::Fermi.Environments.NoAccelerator)

    mol = Fermi.Geometry.to_lints_molecule(molecule)
    bas = Lints.BasisSet(basis, mol)
    @lints begin
        S = Lints.make_S(bas; normalize=normalize)
    end
    S,bas
end
#function aodipole(molecule::Molecule, basis::String)#, interconnect::Fermi.Environments.No_IC,
#                                                     #   communicator::Fermi.Environments.NoCommunicator,
#                                                     #   accelerator::Fermi.Environments.NoAccelerator)
#
#    open("/tmp/molfile.xyz","w") do molfile
#        natom = length(molecule.atoms)
#        write(molfile,"$natom\n\n")
#        write(molfile,Fermi.Geometry.get_xyz(molecule))
#    end
#
#    Lints.libint2_init()
#    mol = Lints.Molecule("/tmp/molfile.xyz")
#    bas = Lints.BasisSet(basis, mol)
#    nprim = Lints.max_nprim(bas)
#    l = Lints.max_l(bas)
#    S_engine = Lints.DipoleEngine(nprim,l)
#    sz = Lints.nao(bas)
#    S = zeros(sz,sz)
#    Lints.make_2D(S,S_engine,bas)
#    Lints.libint2_finalize()
#    S_engine = nothing
#    GC.gc()
#
#    S,bas
#end
"""
    aonuclear(molecule::Fermi.Molecule, basis::String)

Computes AO basis nuclear attraction ⟨μ|V̂|ν⟩ integrals for the given basis and molecule.
Can be accessed at a higher level by calling

    helper["V"]

where `helper` is bound to the desired molecule and basis set.
"""
function aonuclear(molecule::Molecule, basis::String; normalize=false)#, interconnect::Fermi.Environments.No_IC,
                                                     #   communicator::Fermi.Environments.NoCommunicator,
                                                     #   accelerator::Fermi.Environments.NoAccelerator)

    mol = Fermi.Geometry.to_lints_molecule(molecule)
    bas = Lints.BasisSet(basis, mol)
    @lints begin
        V = Lints.make_V(bas; normalize=normalize)
    end
    V,bas
end
"""
    aoeri(molecule::Fermi.Molecule, basis::String)

Computes AO basis electron repulsion integrals ⟨μν|Ô₂|ρσ⟩ integrals for the given basis and molecule.
Can be accessed at a higher level by calling

    helper["μ"]

where `helper` is bound to the desired molecule and basis set.
"""
function aoeri(molecule::Molecule, basis::String; normalize=false)#, interconnect::Fermi.Environments.No_IC,
                                                 #       communicator::Fermi.Environments.NoCommunicator,
                                                 #       accelerator::Fermi.Environments.NoAccelerator)

    mol = Fermi.Geometry.to_lints_molecule(molecule)
    bas = Lints.BasisSet(basis,mol)
    @lints begin
        I = Lints.make_ERI4(bas; normalize=normalize)
    end
    I,bas
end

"""
    dfaoeri(molecule::Fermi.Molecule, basis::String)

Computes AO basis density fitted electron repulsion integrals ⟨μν|Ô₂|P⟩J(P,Q)^-1/2 integrals for the given basis and molecule.
Note that the returned integrals DO NOT need to be combined with the Coulomb metric J(P,Q). In common notation, this is B(Q,μ,ν).
Can be accessed at a higher level by calling

    helper["B"]

where `helper` is bound to the desired molecule and basis set.
"""
function dfaoeri(molecule::Molecule, bname::String,dfbname::String; normalize=false)
    mol = Fermi.Geometry.to_lints_molecule(molecule)
    bas = Lints.BasisSet(bname,mol)
    dfbas = Lints.BasisSet(dfbname,mol)
    @lints begin
        Pqp = Lints.make_ERI3(bas,dfbas; normalize=normalize)
        J = Lints.make_ERI2(dfbas; normalize=normalize)
    end
    Jh = Array(Hermitian(J)^(-1/2)) #sometimes Jh becomes complex slightly if J is not ~~exactly~~ hermitian 💔
    Pqp,Jh
end

function transform_fock(F::Array{Float64,2}, O1::O, O2::O) where O <: AbstractOrbitals
    C1 = hcat([orb.C for orb in O1.orbs]...)
    C2 = hcat([orb.C for orb in O2.orbs]...)
    transform_fock(F,C1,C2)
end

function transform_fock(F::Array{Float64,2}, C1::Array{Float64,2}, C2::Array{Float64,2})

    nmo,p = size(C1)
    _,q   = size(C2)

    Q1 = zeros(p,nmo)
    Fermi.contract!(Q1,C1,F,"pv","up","uv")

    Q2 = zeros(p,q)
    Fermi.contract!(Q2,C2,Q1,"pq","vq","pv")

    return Q2
end

function transform_eri(ERI::Array{T,4}, O1::O, O2::O, O3::O, O4::O) where { O <: AbstractOrbitals,
                                                                            T <: AbstractFloat }
    C1 = hcat([orb.C for orb in O1.orbs]...)
    C2 = hcat([orb.C for orb in O2.orbs]...)
    C3 = hcat([orb.C for orb in O3.orbs]...)
    C4 = hcat([orb.C for orb in O4.orbs]...)
    transform_eri(ERI,C1,C2,C3,C4)
end

function transform_eri(ERI::Array{T,3}, O1::O, O2::O) where { O <: AbstractOrbitals,
                                                              T <: AbstractFloat }
    C1 = hcat([orb.C for orb in O1.orbs]...)
    C2 = hcat([orb.C for orb in O2.orbs]...)
    transform_eri(ERI,C1,C2)
end

function transform_eri(ERI::Array{T,4}, C1::Array{Float64,2}, C2::Array{Float64,2}, C3::Array{Float64,2}, C4::Array{Float64,2}) where T <: AbstractFloat

    nmo,i = size(C1)
    _,a   = size(C2)
    _,j   = size(C3)
    _,b   = size(C4)

    C1 = T.(C1)
    C2 = T.(C2)
    C3 = T.(C3)
    C4 = T.(C4)
    Q1 = zeros(T,i,nmo,nmo,nmo)
    Fermi.contract!(Q1,C1,ERI,"ivls","ui","uvls")

    Q2 = zeros(T,i,a,nmo,nmo)
    Fermi.contract!(Q2,C2,Q1,"ials","va","ivls")
    Q1 = nothing

    Q3 = zeros(T,i,a,j,nmo)
    Fermi.contract!(Q3,C3,Q2,"iajs","lj","ials")
    Q2 = nothing

    Q4 = zeros(T,i,a,j,b)
    Fermi.contract!(Q4,C4,Q3,"iajb","sb","iajs")
    Q4
end

function transform_eri(ERI::Array{T,3},C1::Array{Float64,2},C2::Array{Float64,2}) where T <: AbstractFloat
    C1 = T.(C1)
    C2 = T.(C2)
    naux = size(ERI,1)
    nao,i = size(C1)
    _,a = size(C2)
    Bin = zeros(T,naux,i,nao)
    Fermi.contract!(Bin,C1,ERI,"Qiv","ui","Quv")
    Bia = zeros(T,naux,i,a)
    Fermi.contract!(Bia,C2,Bin,"Qia","va","Qiv")
    Bia
end

end #module


text:
using Test
using NeXLParticle
using StaticArrays
using GeometryBasics
using Random
using LinearAlgebra
using Rotations
#using Gadfly

@testset "rough_align" begin
    @testset "rough_align (BADFOOD)" begin
        rnd = MersenneTwister(0xBADF00D)
        npts = 1000
        bounds = Rect2((-5.0, -1.5), (10.0,3.0))
        xy = NeXLParticle.generate_ground_truth(rnd, bounds, npts)
        ps1 = NeXLParticle.measure_particles(rnd, xy, SA[-10.0,5.0], deg2rad(12.0), 0.8, 0.003)
        ps2 = NeXLParticle.measure_particles(rnd, xy, SA[-2.0,-1.0], deg2rad(213.0), 0.78, 0.003)

        (ct1, ct2) = rough_align(ps1, ps2, tol=0.001)
        @test isapprox(ct1.linear, LinearAlgebra.I, atol=1.0e-8)
        @test isapprox(ct1.translation, [9.463509135855892, -5.084006487177425], atol=1.0e-8)
        @test isapprox(ct2.linear, [-0.933518161752927 -0.35853011265085666; 0.35853011265085666 -0.933518161752927], atol=1.0e-8)
        @test isapprox(ct2.translation, [-2.7619203122160743, -0.2997257033588756], atol=1.0e-8 )
        @test isapprox(RotMatrix{2}(deg2rad(12.0-213.0)), ct2.linear, atol = 0.002)

        ts1, ts2 = ct1.(ps1), ct2.(ps2)
        c1, c2 = correspondences(ts1, ts2; tol=0.01)
        @test norm(ts1[c1[1]]-ts2[c2[1]]) < 0.01
        @test norm(ts1[c1[2]]-ts2[c2[2]]) < 0.01
        @test norm(ts1[c1[3]]-ts2[c2[3]]) < 0.01
        @test norm(ts1[c1[end]]-ts2[c2[end]]) < 0.01

        # plot( ( layer(x=map(p->p[1], t), y=map(p->p[2], t), Geom.point, Theme(default_color=c)) for (t,c) in zip((ts1, ts2), ("red", "blue"))) ...)
    end

    @testset "rough_align (0xABCDEF)" begin
        rnd = MersenneTwister(0xABCDEF)
        npts = 2000
        bounds = Rect2((-10.0, -10.0), (20.0,20.0))
        xy = NeXLParticle.generate_ground_truth(rnd, bounds, npts)
        ps1 = NeXLParticle.measure_particles(rnd, xy, SA[2.0,4.0], deg2rad(230.0), 0.99, 0.003)
        ps2 = NeXLParticle.measure_particles(rnd, xy, SA[-2.0,-6.0], deg2rad(23.0), 0.99, 0.003)

        (ct1, ct2) = NeXLParticle.rough_align(ps1, ps2, tol=0.001)

        @test isapprox(ct1.linear, LinearAlgebra.I, atol=1.0e-8)
        @test isapprox(ct1.translation, [-2.439381678523621, -4.83348719078295], atol=1.0e-8)
        @test isapprox(ct2.linear, [-0.8909816833778803 0.4540392492782079; -0.4540392492782079 -0.8909816833778803], atol=1.0e-8)
        @test isapprox(ct2.translation, [0.5030208959078659, -7.088192398044117], atol=1.0e-8 )
        @test isapprox(RotMatrix{2}(deg2rad(230.0-23.0)), ct2.linear, atol = 0.001)

        ts1, ts2 = ct1.(ps1), ct2.(ps2)
        c1, c2 = NeXLParticle.correspondences(ts1, ts2; tol=0.01)
        @test norm(ts1[c1[1]]-ts2[c2[1]]) < 0.01
        @test norm(ts1[c1[2]]-ts2[c2[2]]) < 0.01
        @test norm(ts1[c1[3]]-ts2[c2[3]]) < 0.01
        @test norm(ts1[c1[end]]-ts2[c2[end]]) < 0.01

        # plot( ( layer(x=map(p->p[1], t), y=map(p->p[2], t), Geom.point, Theme(default_color=c)) for (t,c) in zip((ts1, ts2), ("red", "blue"))) ...)
    end

    @testset "rough_align (0xF156F00D)" begin
        rnd = MersenneTwister(0xF156F00D)
        npts = 100
        bounds = Rect2((-10.0, -10.0), (20.0,20.0))
        xy = NeXLParticle.generate_ground_truth(rnd, bounds, npts)
        ps1 = NeXLParticle.measure_particles(rnd, xy, SA[2.0,4.0], deg2rad(230.0), 0.89, 0.003)
        ps2 = NeXLParticle.measure_particles(rnd, xy, SA[-2.0,-6.0], deg2rad(23.0), 0.94, 0.003)

        (ct1, ct2) = NeXLParticle.rough_align(ps1, ps2, tol=0.001)

        @test isapprox(ct1.linear, LinearAlgebra.I, atol=1.0e-8)
        @test isapprox(ct1.translation, [-3.600666188303255, -3.823165241126716], atol=1.0e-8)
        @test isapprox(ct2.linear, [-0.8910316678635093 0.4539411491199855; -0.4539411491199855 -0.8910316678635093], atol=1.0e-8)
        @test isapprox(ct2.translation, [-0.6588973878728955, -6.076703669619125], atol=1.0e-8 )
        @test isapprox(RotMatrix{2}(deg2rad(230.0-23.0)), ct2.linear, atol = 0.001)

        ts1, ts2 = ct1.(ps1), ct2.(ps2)
        c1, c2 = NeXLParticle.correspondences(ts1, ts2; tol=0.01)
        @test norm(ts1[c1[1]]-ts2[c2[1]]) < 0.01
        @test norm(ts1[c1[2]]-ts2[c2[2]]) < 0.01
        @test norm(ts1[c1[3]]-ts2[c2[3]]) < 0.01
        @test norm(ts1[c1[end]]-ts2[c2[end]]) < 0.01

        # plot( ( layer(x=map(p->p[1], t), y=map(p->p[2], t), Geom.point, Theme(default_color=c)) for (t,c) in zip((ts1, ts2), ("red", "blue"))) ...)
    end

    @testset "rough_align (0x7357)" begin
        rnd = MersenneTwister(0x7357)
        npts = 10000
        bounds = Rect2((-5.0, -5.0), (10.0,10.0))
        xy = NeXLParticle.generate_ground_truth(rnd, bounds, npts)
        ps1 = NeXLParticle.measure_particles(rnd, xy, SA[2.0,4.0], deg2rad(113.0), 0.7, 0.001)
        ps2 = NeXLParticle.measure_particles(rnd, xy, SA[-2.0,-6.0], deg2rad(213.0), 0.65, 0.001)

        (ct1, ct2) = NeXLParticle.rough_align(ps1, ps2, tol=0.001)

        @test isapprox(ct1.linear, LinearAlgebra.I, atol=1.0e-8)
        @test isapprox(ct1.translation, [-1.9797624240598497, -5.2545418722493595], atol=1.0e-8)
        @test isapprox(ct2.linear, [-0.17353087223345756 0.9848284299216262; -0.9848284299216262 -0.17353087223345756], atol=1.0e-8)
        @test isapprox(ct2.translation, [5.582172374918754, -4.264693091414748], atol=1.0e-8 )
        @test isapprox(RotMatrix{2}(deg2rad(113.0-213.0)), ct2.linear, atol = 0.001)

        ts1, ts2 = ct1.(ps1), ct2.(ps2)
        c1, c2 = NeXLParticle.correspondences(ts1, ts2; tol=0.01)
        @test length(c1) == 4668 # approx 10000*0.7*0.67
        @test norm(ts1[c1[1]]-ts2[c2[1]]) < 0.01
        @test norm(ts1[c1[2]]-ts2[c2[2]]) < 0.01
        @test norm(ts1[c1[3]]-ts2[c2[3]]) < 0.01
        @test norm(ts1[c1[end]]-ts2[c2[end]]) < 0.01

        #plot( ( layer(x=map(p->p[1], t), y=map(p->p[2], t), Geom.point, Theme(default_color=c)) for (t,c) in zip((ts1, ts2), ("red", "blue"))) ...)
    end
end

ext: jl

max_stars_repo_name: NicholasWMRitchie/NeXLParticle

max_issues_repo_name: NicholasWMRitchie/NeXLParticle

max_forks_repo_name: NicholasWMRitchie/NeXLParticle


text:
struct Point{T <: Integer}
    x::T
    y::T
end

function Base.:(==)(p1::Point, p2::Point)
    return p1.x == p2.x && p1.y == p2.y
end

function Base.show(io::IO, p::Point{T}) where {T}
    print("Point{$T}(x=$(p.x), y=$(p.y))")
end


text:
# Autogenerated wrapper script for GAP_pkg_cvec_jll for x86_64-linux-musl
export cvec

using GAP_jll
using GAP_lib_jll
JLLWrappers.@generate_wrapper_header("GAP_pkg_cvec")
JLLWrappers.@declare_file_product(cvec)
function __init__()
    JLLWrappers.@generate_init_header(GAP_jll, GAP_lib_jll)
    JLLWrappers.@init_file_product(
        cvec,
        "lib/gap/cvec.so",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


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

using Test
using SubModule

@test SubModule.f() == 1

max_stars_repo_name: greimel/julia

max_issues_repo_name: greimel/julia

max_forks_repo_name: greimel/julia


text:
using DifferentialEquations
using Plots; gr()

f(u,p,t) = p[1]*u
u0 = 0.5
#params = [3.0]
p = [3.0]
tspan = (0.0,1.0)
prob = ODEProblem(f,u0,tspan,p)

sol = solve(prob,reltol=1e-7,abstol=1e-7)

plot(sol)

plot(sol,linewidth=5,title="Exponential Growth",
     xaxis="Time",yaxis="Population size",label="Numerical solution")

plot!(sol.t,t->u0*exp(p[1]t),lw=3,ls=:dash,label="True Solution", legend =:bottomright)

############################################################################

f(u,p,t) = p[1]*u*(1-(u/p[2]))
#f(u,p,t) = p[1]*u*(1-u)
u0 = 0.005
#params = [5.0,1.0]
p = [1.2,1.0]
tspan = (0.0, 10.0)
#prob2 = ODEProblem(f, u0, tspan, params)
prob2 = ODEProblem(f, u0, tspan, p)

sol2 = solve(prob2,reltol=1e-8,abstol=1e-8)

plot(sol2)

plot(sol2,linewidth=5,title="Logistic Growth",
     xaxis="Time",yaxis="Population size",label="Numerical Solution", legend=:bottomright)

plot!(sol2.t,t->((p[2])*(u0*exp(p[1]t)))/(p[2]+(u0*(exp(p[1]t)-1))),lw=3,ls=:dash,label="True Solution")


text:
config = MOI.Test.TestConfig()
optimize1!(mock) = MOIU.mock_optimize!(mock, MOI.INFEASIBLE, tuple(), MOI.INFEASIBILITY_CERTIFICATE)
# The test does not check the solution so we just set zeros.
optimize2!(mock) = MOIU.mock_optimize!(mock, zeros(MOI.get(mock, MOI.NumberOfVariables())))
for mock in mocks(optimize1!, optimize2!)
    Tests.motzkin_test(mock, config)
end


text:
##---begin-Plot-functions---------------------------------------------------------------------------
"""
    chromatogram(t, tR, τR)

Calculate the chromatogram as a sum of gaussian peaks over the time `t` for peaks centered at retention times `tR` and with peak width `τR`.
"""
function chromatogram(t::Array{Float64,1}, tR::Array{Float64,1}, τR::Array{Float64,1})
	g(t,tR,τR) = 1/sqrt(2*π*τR^2)*exp(-(t-tR)^2/(2*τR^2))
	chromatograms = Array{Array{Float64,1}}(undef, length(tR))
	for j=1:length(tR)
		chromatograms[j] = g.(t, tR[j], τR[j])
	end
	return sum(chromatograms)
end

"""
    plot_chromatogram(peaklist, tlims; annotation=true, number=true, mirror=false, offset=0.0)

Plot the chromatogram of the peaks listed in `peaklist` over the time tupel `tlims = (t_start, t_end`).

# Arguments
* `peaklist`: DataFrame with the names, retention times and peak widths of the simulated substances.
* `tlims`: Tuple defining the start and end time of the plotted chromatogram.
* `annotation`: Boolean, switching the annotation of the peaks on/off; default = true.
* `number`: Boolean, switching the type of the annotation between _number_ of the substance in the peaklist (`number = true`) or _name_ of the substance (`number = false`); default = true.
* `mirror`: Boolean, if `mirror = true` the chromatogram is multiplied by `-1`; default = false.
* `offset`: Float64, this value is added to the chromatogram; default = 0.0.

# Output
Tupel `(p_chrom, t, chrom)`
* `p_chrom`: the plot of the chromatogram `chrom` over time `t`
* `t`: Array of time of the chromatogram
* `chrom`: Array of the abundance values of the chromatogram
"""
function plot_chromatogram(peaklist, tlims; annotation=true, number=true, mirror=false, offset=0.0)
	t₀ = tlims[1]
	tMax = tlims[2]
	t = t₀:tMax/10000:tMax
	if mirror==true
		sgn = -1.0
	else
		sgn = 1.0
	end
	chrom = sgn.*chromatogram(collect(t), peaklist.tR, peaklist.τR) .+ offset
	chrom_tR = sgn.*chromatogram(peaklist.tR, peaklist.tR, peaklist.τR) .+ offset
	p_chrom = plot(t, chrom, xlabel="time in s", ylabel="abundance", legend=false)
	if annotation==true && number==true
		plot!(p_chrom, annotations = [(peaklist.tR[i], chrom_tR[i], text(i, 10, rotation=0, :center)) for i in 1:length(peaklist.tR)])
	elseif annotation==true && number==false
		plot!(p_chrom, annotations = [(peaklist.tR[i], chrom_tR[i], text(peaklist.Name[i], 10, rotation=90, :center)) for i in 1:length(peaklist.tR)])
	end
	return p_chrom, t, chrom
end

"""
    plot_chromatogram!(p_chrom, peaklist, tlims; annotation=true, number=true, mirror=false, offset=0.0)

Add the chromatogram of the peaks listed in `peaklist` over the time tupel `tlims = (t_start, t_end`) to the plot `p_chrom`.

# Arguments
* `p_chrom`: Plot of an existing chromatogram.
* `peaklist`: DataFrame with the names, retention times and peak widths of the simulated substances.
* `tlims`: Tuple defining the start and end time of the plotted chromatogram.
* `annotation`: Boolean, switching the annotation of the peaks on/off; default = true.
* `number`: Boolean, switching the type of the annotation between _number_ of the substance in the peaklist (`number = true`) or _name_ of the substance (`number = false`); default = true.
* `mirror`: Boolean, if `mirror = true` the chromatogram is multiplied by `-1`; default = false.
* `offset`: Float64, this value is added to the chromatogram; default = 0.0.

# Output
Tupel `(t, chrom)`
* `t`: Array of time of the chromatogram
* `chrom`: Array of the abundance values of the chromatogram
"""
function plot_chromatogram!(p_chrom, peaklist, tlims; t₀=0.0, annotation=true, number=true, mirror=false, offset=0.0)
	t₀ = tlims[1]
	tMax = tlims[2]
	t = t₀:tMax/10000:tMax
	if mirror==true
		sgn = -1.0
	else
		sgn = 1.0
	end
	chrom = sgn.*chromatogram(collect(t), peaklist.tR, peaklist.τR) .+ offset
	chrom_tR = sgn.*chromatogram(peaklist.tR, peaklist.tR, peaklist.τR) .+ offset
	plot!(p_chrom, t, chrom)
	if annotation==true && number==true
		plot!(p_chrom, annotations = [(peaklist.tR[i], chrom_tR[i], text(i, 10, rotation=0, :center)) for i in 1:length(peaklist.tR)])
	elseif annotation==true && number==false
		plot!(p_chrom, annotations = [(peaklist.tR[i], chrom_tR[i], text(peaklist.Name[i], 10, rotation=90, :center)) for i in 1:length(peaklist.tR)])
	end
	return t, chrom
end

"""
	plot_flow(par)

Calculate and plot the flow (in mL/min, normalized) of the carrier gas in a GC Column with a program defined in the parameters `par::GasChromatography.Parameters`.
"""
function plot_flow(par)
	t = 0.0:sum(par.prog.time_steps)/1000.0:sum(par.prog.time_steps)
	F = Array{Float64}(undef, length(t))
	for i=1:length(t)
		F[i] = flow(t[i], par.prog.T_itp, par.prog.Fpin_itp, par.prog.pout_itp, par.col.L, par.col.d, par.col.gas; ng=par.opt.ng, vis=par.opt.vis, control=par.opt.control)
	end
	p_flow = plot(t, F.*60e6, xlabel="time in s", ylabel="column flow in mL/min", legend=false)
	return p_flow
end

"""
	plot_pressure(prog)

Plot the inlet and outlet pressure over time of the program `prog::GasChromatographySimulator.Program`.
"""
function plot_pressure(par)
	t = 0.0:sum(par.prog.time_steps)/1000.0:sum(par.prog.time_steps)
	pin = Array{Float64}(undef, length(t))
	pout = Array{Float64}(undef, length(t))
	for i=1:length(t)
		pin[i] = inlet_pressure(t[i], par.prog.T_itp, par.prog.Fpin_itp, par.prog.pout_itp, par.col.L, par.col.d, par.col.gas; ng=par.opt.ng, vis=par.opt.vis, control=par.opt.control)
		pout[i] = par.prog.pout_itp(t[i])
	end
	p_press = plot(t, pin, xlabel="time in s", ylabel="pressure in Pa", label="inlet", legend=:right)
	plot!(p_press, t, pout, label="outlet")
	return p_press
end

"""
	plot_temperature(par; selector="T(t)")

Plot the temperature program of the GC Column.

# Arguments
* `par::GasChromatographySimulator.Parameters`: parameters of the GC system
* `selector::String`: selection of the plot as:
	* `selector = "T(t)"`: 2D-plot of temperature `T` over time `t` at inlet (`x=0`) and outlet (`x=L`)
	* `selector = "T(x)"`: 2D-plot of temperature `T` over column position `x` at the `time_steps` of the program
	* `selector = "T(x,t)"`: 3D-plot of temperature `T` over column position `x` and `t`
"""
function plot_temperature(par::GasChromatographySimulator.Parameters; selector="T(t)")
	if selector=="T(x)"
		nx = 0.0:par.col.L/1000:par.col.L
		p_temp = plot(xlabel="x in m", ylabel="T in °C", legend=:top)
		for i=1:length(par.prog.time_steps)
			T = par.prog.T_itp.(nx, cumsum(par.prog.time_steps)[i]).-273.15
			plot!(p_temp, nx, T, label="t=$(cumsum(par.prog.time_steps)[i])s")
		end
	elseif selector=="T(t)"
		nt = 0.0:sum(par.prog.time_steps)/1000:sum(par.prog.time_steps)
		T0 = par.prog.T_itp.(0.0, nt).-273.15
		p_temp = plot(nt, T0, xlabel="t in s", ylabel="T in °C", legend=:top, label="inlet", c=:red)
		TL = par.prog.T_itp.(par.col.L, nt).-273.15
		plot!(p_temp, nt, TL, label="outlet", c=:blue)
	elseif selector=="T(x,t)"
		nx = 0.0:par.col.L/1000:par.col.L
		nt = 0.0:sum(par.prog.time_steps)/1000:sum(par.prog.time_steps)
		T = Array{Float64}(undef, length(nx), length(nt))
		for j=1:length(nt)
			for i=1:length(nx)
				T[i,j] = par.prog.T_itp(nx[i], nt[j])-273.15
			end
		end
		p_temp = plot(nx, nt, T', st=:surface, xlabel="x in m", ylabel="t in s", zlabel="T in °C")
	end
	return p_temp
end

# functions from GasChromatographyTools.jl:
"""
	local_plots(xx, yy, sol, par)

Show additional 'local' plots of selected `yy` quantities over selected `xx`
quantities.

# Arguments
* `xx`: Selected quantity shown on the x-axis. Possible values: "z", "t", "T",
  "τ", "σ" and "u".
* `yy`: Selected quantity shown on the y-axis. Possible values: "z", "t", "T",
  "τ", "σ" and "u".
* `sol`: The solution of the simulation.
* `par`: The parameters of the simulated GC-system.
"""
function local_plots(xx, yy, sol, par)
	n = size(sol)[1]

	df_sol = GasChromatographySimulator.sol_extraction(sol, par)
	xvalues = Array{Array{Float64,1}}(undef, n)
	yvalues = Array{Array{Float64,1}}(undef, n)

	p_add = plot(legend=false)
	for i=1:n
		if xx=="z"
			xvalues[i] = df_sol.z[i]
			xlabel = "position z in m"
		elseif xx=="t"
			xvalues[i] = df_sol.t[i]
			xlabel = "time t in s"
		elseif xx=="T"
			xvalues[i] = par.prog.T_itp.(df_sol.z[i], df_sol.t[i]).-273.15
			xlabel = "temperature T in °C"
		elseif xx=="τ"
			xvalues[i] = sqrt.(df_sol.τ²[i])
			xlabel = "peak width τ in s"
		elseif xx=="σ"
			xvalues[i] = velocity(df_sol, i, par).*sqrt.(df_sol.τ²[i])
			xlabel = "band width in m"
		elseif xx=="u"
			xvalues[i] = velocity(df_sol, i, par)
			xlabel = "solute velocity in m/s"
		end
		if yy=="z"
			yvalues[i] = df_sol.z[i]
			ylabel = "position z in m"
		elseif yy=="t"
			yvalues[i] = df_sol.t[i]
			ylabel = "time t in s"
		elseif yy=="T"
			yvalues[i] = par.prog.T_itp.(df_sol.z[i], df_sol.t[i]).-273.15
			ylabel = "temperature T in °C"
		elseif yy=="τ"
			yvalues[i] = sqrt.(df_sol.τ²[i])
			ylabel = "peak width in s"
		elseif yy=="σ"
			yvalues[i] = velocity(df_sol, i, par).*sqrt.(df_sol.τ²[i])
			ylabel = "band width in m"
		elseif yy=="u"
			yvalues[i] = velocity(df_sol, i, par)
			ylabel = "solute velocity in m/s"
		end
		plot!(p_add, xvalues[i], yvalues[i], xlabel=xlabel, ylabel=ylabel, label=par.sub[i].name, m=:o)
	end
	return p_add
end

"""
	velocity(df_sol, i, par)

Calculate the velocity (in m/s) coressponding to solution of the `i-th` sunstance of a
GC-system defined by `par`.
"""
function velocity(df_sol, i, par)
	x = df_sol.z[i]
	t = df_sol.t[i]
	T_itp = par.prog.T_itp
	pin_itp = par.prog.Fpin_itp
	pout_itp = par.prog.pout_itp
	L = par.col.L
	d = par.col.d
	df = par.col.df
	gas = par.col.gas
	ΔCp = par.sub[i].ΔCp
	Tchar = par.sub[i].Tchar
	θchar = par.sub[i].θchar
	φ₀ = par.sub[i].φ₀
	u = Array{Float64}(undef, length(x))
	for j=1:length(x)
		u[j] = 1/GasChromatographySimulator.residency(x[j], t[j], T_itp, pin_itp, pout_itp, L, d, df, gas, ΔCp, Tchar, θchar, φ₀; ng=par.opt.ng, vis=par.opt.vis, control=par.opt.control)
	end
	return u
end


text:
# SHA initial hash values and constants

# Hash constant words K for SHA1
const K1 = UInt32[
    0x5a827999, 0x6ed9eba1, 0x8f1bbcdc, 0xca62c1d6
]

# Initial hash value H for SHA1
const SHA1_initial_hash_value = UInt32[
    0x67452301, 0xefcdab89, 0x98badcfe, 0x10325476, 0xc3d2e1f0
]

# Hash constant words K for SHA-256:
const K256 = UInt32[
    0x428a2f98, 0x71374491, 0xb5c0fbcf, 0xe9b5dba5,
    0x3956c25b, 0x59f111f1, 0x923f82a4, 0xab1c5ed5,
    0xd807aa98, 0x12835b01, 0x243185be, 0x550c7dc3,
    0x72be5d74, 0x80deb1fe, 0x9bdc06a7, 0xc19bf174,
    0xe49b69c1, 0xefbe4786, 0x0fc19dc6, 0x240ca1cc,
    0x2de92c6f, 0x4a7484aa, 0x5cb0a9dc, 0x76f988da,
    0x983e5152, 0xa831c66d, 0xb00327c8, 0xbf597fc7,
    0xc6e00bf3, 0xd5a79147, 0x06ca6351, 0x14292967,
    0x27b70a85, 0x2e1b2138, 0x4d2c6dfc, 0x53380d13,
    0x650a7354, 0x766a0abb, 0x81c2c92e, 0x92722c85,
    0xa2bfe8a1, 0xa81a664b, 0xc24b8b70, 0xc76c51a3,
    0xd192e819, 0xd6990624, 0xf40e3585, 0x106aa070,
    0x19a4c116, 0x1e376c08, 0x2748774c, 0x34b0bcb5,
    0x391c0cb3, 0x4ed8aa4a, 0x5b9cca4f, 0x682e6ff3,
    0x748f82ee, 0x78a5636f, 0x84c87814, 0x8cc70208,
    0x90befffa, 0xa4506ceb, 0xbef9a3f7, 0xc67178f2
]

# Initial hash value H for SHA-224:
const SHA2_224_initial_hash_value = UInt32[
    0xc1059ed8, 0x367cd507, 0x3070dd17, 0xf70e5939,
    0xffc00b31, 0x68581511, 0x64f98fa7, 0xbefa4fa4
]

const SHA2_256_initial_hash_value = UInt32[
    0x6a09e667, 0xbb67ae85, 0x3c6ef372, 0xa54ff53a,
    0x510e527f, 0x9b05688c, 0x1f83d9ab, 0x5be0cd19
]

# Hash constant words K for SHA-384 and SHA-512:
const K512 = UInt64[
    0x428a2f98d728ae22, 0x7137449123ef65cd,
    0xb5c0fbcfec4d3b2f, 0xe9b5dba58189dbbc,
    0x3956c25bf348b538, 0x59f111f1b605d019,
    0x923f82a4af194f9b, 0xab1c5ed5da6d8118,
    0xd807aa98a3030242, 0x12835b0145706fbe,
    0x243185be4ee4b28c, 0x550c7dc3d5ffb4e2,
    0x72be5d74f27b896f, 0x80deb1fe3b1696b1,
    0x9bdc06a725c71235, 0xc19bf174cf692694,
    0xe49b69c19ef14ad2, 0xefbe4786384f25e3,
    0x0fc19dc68b8cd5b5, 0x240ca1cc77ac9c65,
    0x2de92c6f592b0275, 0x4a7484aa6ea6e483,
    0x5cb0a9dcbd41fbd4, 0x76f988da831153b5,
    0x983e5152ee66dfab, 0xa831c66d2db43210,
    0xb00327c898fb213f, 0xbf597fc7beef0ee4,
    0xc6e00bf33da88fc2, 0xd5a79147930aa725,
    0x06ca6351e003826f, 0x142929670a0e6e70,
    0x27b70a8546d22ffc, 0x2e1b21385c26c926,
    0x4d2c6dfc5ac42aed, 0x53380d139d95b3df,
    0x650a73548baf63de, 0x766a0abb3c77b2a8,
    0x81c2c92e47edaee6, 0x92722c851482353b,
    0xa2bfe8a14cf10364, 0xa81a664bbc423001,
    0xc24b8b70d0f89791, 0xc76c51a30654be30,
    0xd192e819d6ef5218, 0xd69906245565a910,
    0xf40e35855771202a, 0x106aa07032bbd1b8,
    0x19a4c116b8d2d0c8, 0x1e376c085141ab53,
    0x2748774cdf8eeb99, 0x34b0bcb5e19b48a8,
    0x391c0cb3c5c95a63, 0x4ed8aa4ae3418acb,
    0x5b9cca4f7763e373, 0x682e6ff3d6b2b8a3,
    0x748f82ee5defb2fc, 0x78a5636f43172f60,
    0x84c87814a1f0ab72, 0x8cc702081a6439ec,
    0x90befffa23631e28, 0xa4506cebde82bde9,
    0xbef9a3f7b2c67915, 0xc67178f2e372532b,
    0xca273eceea26619c, 0xd186b8c721c0c207,
    0xeada7dd6cde0eb1e, 0xf57d4f7fee6ed178,
    0x06f067aa72176fba, 0x0a637dc5a2c898a6,
    0x113f9804bef90dae, 0x1b710b35131c471b,
    0x28db77f523047d84, 0x32caab7b40c72493,
    0x3c9ebe0a15c9bebc, 0x431d67c49c100d4c,
    0x4cc5d4becb3e42b6, 0x597f299cfc657e2a,
    0x5fcb6fab3ad6faec, 0x6c44198c4a475817
]

# Initial hash value H for SHA-384
const SHA2_384_initial_hash_value = UInt64[
    0xcbbb9d5dc1059ed8, 0x629a292a367cd507,
    0x9159015a3070dd17, 0x152fecd8f70e5939,
    0x67332667ffc00b31, 0x8eb44a8768581511,
    0xdb0c2e0d64f98fa7, 0x47b5481dbefa4fa4
]

# Initial hash value H for SHA-512
const SHA2_512_initial_hash_value = UInt64[
    0x6a09e667f3bcc908, 0xbb67ae8584caa73b,
    0x3c6ef372fe94f82b, 0xa54ff53a5f1d36f1,
    0x510e527fade682d1, 0x9b05688c2b3e6c1f,
    0x1f83d9abfb41bd6b, 0x5be0cd19137e2179
]

# Round constants for SHA3 rounds
const SHA3_ROUND_CONSTS = UInt64[
    0x0000000000000001, 0x0000000000008082, 0x800000000000808a,
    0x8000000080008000, 0x000000000000808b, 0x0000000080000001,
    0x8000000080008081, 0x8000000000008009, 0x000000000000008a,
    0x0000000000000088, 0x0000000080008009, 0x000000008000000a,
    0x000000008000808b, 0x800000000000008b, 0x8000000000008089,
    0x8000000000008003, 0x8000000000008002, 0x8000000000000080,
    0x000000000000800a, 0x800000008000000a, 0x8000000080008081,
    0x8000000000008080, 0x0000000080000001, 0x8000000080008008
]

# Rotation constants for SHA3 rounds
const SHA3_ROTC = UInt64[
    1,  3,  6,  10, 15, 21, 28, 36, 45, 55, 2,  14,
    27, 41, 56, 8,  25, 43, 62, 18, 39, 61, 20, 44
]

# Permutation indices for SHA3 rounds (+1'ed so as to work with julia's 1-based indexing)
const SHA3_PILN = Int[
    11, 8,  12, 18, 19, 4, 6,  17, 9,  22, 25, 5,
    16, 24, 20, 14, 13, 3, 21, 15, 23, 10,  7,  2
]

max_stars_repo_name: djsegal/julia-fork

max_issues_repo_name: djsegal/julia-fork

max_forks_repo_name: djsegal/julia-fork


text:
push!(LOAD_PATH,"../src/")

using Documenter
using CureMIDI

makedocs(
    sitename="CureMIDI.jl",
    pages= [
        "index.md",
        "reference.md"
    ]
    )

deploydocs(
    repo = "github.com/yukota/CureMIDI.jl.git",
    )


text:
# @author  Thomas Lehmann
# @file    Queen.jl
#
# Copyright (c) 2016 Thomas Lehmann
#
# Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated
# documentation files (the "Software"), to deal in the Software without restriction, including without limitation
# the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software,
# and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all copies
# or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,
# INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
# IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
# DAMAGES OR OTHER LIABILITY,
# WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
type Queen
    width::Int64
    columns::Array{Int64, 1}
    diagonals1::Array{Bool, 1}
    diagonals2::Array{Bool, 1}
    solutions::Array{Any, 1}

    run::Function

    function Queen(width::Int64)
        this = new(width, Int64[-1 for c = 1: width],
                   [false for d = 1: width*2], [false for d = 1: width*2], [])

        this.run = function(row::Int64)
            for column = 1: this.width
                if this.columns[column] >= 0
                    continue
                end

                ixDiag1 = row + column
                if this.diagonals1[ixDiag1]
                    continue
                end

                ixDiag2 = width - row + column
                if this.diagonals2[ixDiag2]
                    continue
                end

                this.columns[column] = row
                this.diagonals1[ixDiag1] = true
                this.diagonals2[ixDiag2] = true

                if (row == width)
                    append!(this.solutions, [copy(this.columns)])
                else
                    this.run(row + 1)
                end

                this.columns[column] = -1
                this.diagonals1[ixDiag1] = false
                this.diagonals2[ixDiag2] = false
            end
        end

        return this
    end
end

width  = 8
if length(ARGS) == 1
    width = parse(ARGS[1])
end

output = false
queen  = Queen(width)
println("Queen raster ($(queen.width)x$(queen.width))")
start  = now()
queen.run(1)
println("...took $(Int(now() - start)/1000.0) seconds.")
println("...", length(queen.solutions), " solutions found.")

if output
    for solution in queen.solutions
        line = ""
        for (index, value) in enumerate(solution)
            line = string(line," ($index, $value)")
        end
        println(line)
    end
end

max_stars_repo_name: Nachtfeuer/concept-queen

max_issues_repo_name: Nachtfeuer/concept-queen

max_forks_repo_name: Nachtfeuer/concept-queen


text:
function makelogo()

    defaultfont("Helvetica Neue Bold")

    c, l = canvas(92, 92, color = "transparent")

    colors = Colors.JULIA_LOGO_COLORS

    rect!(l, O, 92, 92) + Fill(colors.blue) + Stroke(nothing)

    layers = map(P.([-34, 0, 28], [-5, 0, 11]), -15:15:15, [1.1, 1.0, 0.9]) do p, ang, scale
        layer(translation = p + X(5), rotation = deg(ang), scale = scale)
    end .+ Opacity(0.9)

    for ll in reverse(layers)
        push!(l, ll)
    end

    map(layers) do l
        rect!(l, O, 35, 35)
    end .+ Fill.([colors.red, colors.green, colors.purple])

    textparts = ["Lay", "ere", "d.jl"]
    map(layers, textparts) do l, t
        txt!(l, O, t, 11, :c, :c)
    end .+ Textfill("white")

    circ = circle!(l, O, 45) + Invisible
    l + Clip(circ)

    c
    mkdir("assets")

    Layered.png(c, "assets/logo.png", px_per_pt = 3)
    Layered.svg(c, "assets/logo.svg")
end


text:
using SkewDist, Distributions, Base.Test

println("Testing MvSkewTDist...")

srand(1)
N = 1000000 # Sample size

ξ = [0.0, -0.1]
α = [1.0, 0.5]
Ω = rand_pd(2)
df = 4.0

dist = MvSkewTDist(ξ, Ω, α, df)
rand(dist)

Z = rand(dist, N)

@test_approx_eq_eps mean(Z, 2) mean(dist) 5e-3
@test_approx_eq_eps var(Z, 2) var(dist) 1e-2
@test_approx_eq_eps Base.covm(Z, mean(Z, 2), 2) cov(dist) 1e-2


text:
#push!(LOAD_PATH,"../src/")

using Documenter
using DyadicKDE

makedocs(
    sitename = "DyadicKDE.jl",
    modules = DyadicKDE,
    pages = [
        "Home" => "index.md",
        "Documentation" => "documentation.md",
    ],
)

deploydocs(
    repo = "github.com/WGUNDERWOOD/DyadicKDE.jl.git",
)


text:
using TensorTrains

include("test_tt_tools.jl")
include("test_als.jl")
include("test_mals.jl")
include("test_ordering_schemes.jl")
include("test_models.jl")

max_stars_repo_name: msdupuy/Tensor-Train-Julia

max_issues_repo_name: msdupuy/Tensor-Train-Julia

max_forks_repo_name: msdupuy/Tensor-Train-Julia


text:
# ---
# title: 862. Shortest Subarray with Sum at Least K
# id: problem862
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Binary Search, Queue
# link: <https://leetcode.com/problems/shortest-subarray-with-sum-at-least-k/description/>
# hidden: true
# ---
#
# Return the **length** of the shortest, non-empty, contiguous subarray of `A`
# with sum at least `K`.
#
# If there is no non-empty subarray with sum at least `K`, return `-1`.
#
#
#
# **Example 1:**
#
#
#
#     Input: A = [1], K = 1
#     Output: 1
#
#
# **Example 2:**
#
#
#
#     Input: A = [1,2], K = 4
#     Output: -1
#
#
# **Example 3:**
#
#
#
#     Input: A = [2,-1,2], K = 3
#     Output: 3
#
#
#
#
# **Note:**
#
#   1. `1 <= A.length <= 50000`
#   2. `-10 ^ 5 <= A[i] <= 10 ^ 5`
#   3. `1 <= K <= 10 ^ 9`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using Manifolds
using BenchmarkTools
using StaticArrays
using HybridArrays
using LinearAlgebra
using Random

# Define a parent BenchmarkGroup to contain our suite
SUITE = BenchmarkGroup()

Random.seed!(12334)

function add_manifold(M::Manifold, pts, name;
    test_tangent_vector_broadcasting = true,
    retraction_methods = [],
    inverse_retraction_methods = [],
    point_distributions = [],
    tvector_distributions = [])

    SUITE["manifolds"][name] = BenchmarkGroup()
    tv = log(M, pts[1], pts[2])
    tv1 = log(M, pts[1], pts[2])
    tv2 = log(M, pts[1], pts[3])
    p = similar(pts[1])
    SUITE["manifolds"][name]["similar"] = @benchmarkable similar($(pts[1]))
    SUITE["manifolds"][name]["log"] = @benchmarkable log($M, $(pts[1]), $(pts[2]))
    SUITE["manifolds"][name]["log!"] = @benchmarkable log!($M, $tv, $(pts[1]), $(pts[2]))
    for iretr ∈ inverse_retraction_methods
        SUITE["manifolds"][name]["inverse_retract: "*string(iretr)] = @benchmarkable inverse_retract($M, $(pts[1]), $(pts[2]), $iretr)
        SUITE["manifolds"][name]["inverse_retract!: "*string(iretr)] = @benchmarkable inverse_retract!($M, $tv, $(pts[1]), $(pts[2]), $iretr)
    end
    SUITE["manifolds"][name]["exp"] = @benchmarkable exp($M, $(pts[1]), $tv1)
    SUITE["manifolds"][name]["exp!"] = @benchmarkable exp!($M, $p, $(pts[1]), $tv1)
    for retr ∈ retraction_methods
        SUITE["manifolds"][name]["retract: "*string(retr)] = @benchmarkable retract($M, $(pts[1]), $tv1, $retr)
        SUITE["manifolds"][name]["retract!: "*string(retr)] = @benchmarkable retract!($M, $p, $(pts[1]), $tv1, $retr)
    end
    SUITE["manifolds"][name]["norm"] = @benchmarkable norm($M, $(pts[1]), $tv1)
    SUITE["manifolds"][name]["inner"] = @benchmarkable inner($M, $(pts[1]), $tv1, $tv2)
    SUITE["manifolds"][name]["distance"] = @benchmarkable distance($M, $(pts[1]), $(pts[2]))
    SUITE["manifolds"][name]["isapprox (pt)"] = @benchmarkable isapprox($M, $(pts[1]), $(pts[2]))
    SUITE["manifolds"][name]["isapprox (tv)"] = @benchmarkable isapprox($M, $(pts[1]), $tv1, $tv2)
    SUITE["manifolds"][name]["2 * tv1 + 3 * tv2"] = @benchmarkable 2 * $tv1 + 3 * $tv2
    SUITE["manifolds"][name]["tv = 2 * tv1 + 3 * tv2"] = @benchmarkable $tv = 2 * $tv1 + 3 * $tv2
    if test_tangent_vector_broadcasting
        SUITE["manifolds"][name]["tv = 2 .* tv1 .+ 3 .* tv2"] = @benchmarkable $tv = 2 .* $tv1 .+ 3 .* $tv2
        SUITE["manifolds"][name]["tv .= 2 .* tv1 .+ 3 .* tv2"] = @benchmarkable $tv .= 2 .* $tv1 .+ 3 .* $tv2
    end
    for pd ∈ point_distributions
        distr_name = string(pd)
        distr_name = distr_name[1:min(length(distr_name), 50)]
        SUITE["manifolds"][name]["point distribution "*distr_name] = @benchmarkable rand($pd)
    end
    for tvd ∈ point_distributions
        distr_name = string(tvd)
        distr_name = distr_name[1:min(length(distr_name), 50)]
        SUITE["manifolds"][name]["tangent vector distribution "*distr_name] = @benchmarkable rand($tvd)
    end
end

# General manifold benchmarks
function add_manifold_benchmarks()

    SUITE["manifolds"] = BenchmarkGroup()

    s2 = Manifolds.Sphere(2)
    array_s2 = ArrayManifold(s2)
    r2 = Manifolds.Euclidean(2)

    pts_r2 = [Size(2)([1.0, 1.0]),
              Size(2)([-2.0, 3.0]),
              Size(2)([3.0, -2.0])]

    add_manifold(r2,
                 pts_r2,
                 "Euclidean{2} -- SizedArray")

    add_manifold(r2,
                 [MVector{2,Float64}([1.0, 1.0]),
                  MVector{2,Float64}([-2.0, 3.0]),
                  MVector{2,Float64}([3.0, -2.0])],
                  "Euclidean{2} -- MVector")

    ud_sphere = Manifolds.uniform_distribution(s2, Size(3)([1.0, 0.0, 0.0]))
    gtd_sphere = Manifolds.normal_tvector_distribution(s2, Size(3)([1.0, 0.0, 0.0]), 1.0)

    pts_s2 = [Size(3)([1.0, 0.0, 0.0]),
              Size(3)([0.0, 1.0, 0.0]),
              Size(3)([0.0, 0.0, 1.0])]

    add_manifold(s2,
                 pts_s2,
                 "Sphere{2} -- SizedArray";
                 point_distributions = [ud_sphere],
                 tvector_distributions = [gtd_sphere])

    add_manifold(array_s2,
                 [Size(3)([1.0, 0.0, 0.0]),
                  Size(3)([0.0, 1.0, 0.0]),
                  Size(3)([0.0, 0.0, 1.0])],
                 "ArrayManifold{Sphere{2}} -- SizedArray";
                 test_tangent_vector_broadcasting = false)

    retraction_methods_rot = [Manifolds.PolarRetraction(),
                              Manifolds.QRRetraction()]

    inverse_retraction_methods_rot = [Manifolds.PolarInverseRetraction(),
                                      Manifolds.QRInverseRetraction()]

    so2 = Manifolds.Rotations(2)
    angles = (0.0, π/2, 2π/3)
    add_manifold(so2,
                 [Size(2, 2)([cos(ϕ) -sin(ϕ); sin(ϕ) cos(ϕ)]) for ϕ in angles],
                 "Rotations(2) -- SizedArray",
                 retraction_methods = retraction_methods_rot,
                 inverse_retraction_methods = inverse_retraction_methods_rot)

    so3 = Manifolds.Rotations(3)
    angles = (0.0, π/2, 2π/3)
    add_manifold(so3,
                 [Size(3, 3)([cos(ϕ) -sin(ϕ) 0; sin(ϕ) cos(ϕ) 0; 0 0 1]) for ϕ in angles],
                 "Rotations(3) -- SizedArray",
                 retraction_methods = retraction_methods_rot,
                 inverse_retraction_methods = inverse_retraction_methods_rot)

    m_prod = Manifolds.ProductManifold(s2, r2)
    shape_s2r2_array = Manifolds.ShapeSpecification(Manifolds.ArrayReshaper(), s2, r2)
    shape_s2r2_static = Manifolds.ShapeSpecification(Manifolds.StaticReshaper(), s2, r2)

    pts_prd_base = [[1.0, 0.0, 0.0, 0.0, 0.0],
                    [0.0, 1.0, 0.0, 1.0, 0.0],
                    [0.0, 0.0, 1.0, 0.0, 0.1]]
    pts_prod_static = map(p -> Manifolds.ProductArray(shape_s2r2_static, p), pts_prd_base)
    add_manifold(m_prod, pts_prod_static, "ProductManifold with ProductArray (static)")

    pts_prod_array = map(p -> Manifolds.ProductArray(shape_s2r2_array, p), pts_prd_base)
    add_manifold(m_prod, pts_prod_array, "ProductManifold with ProductArray (array)")

    pts_prod_mpoints = [Manifolds.ProductRepr(p[1], p[2]) for p in zip(pts_s2, pts_r2)]
    add_manifold(m_prod, pts_prod_mpoints, "ProductManifold with MPoint";
        test_tangent_vector_broadcasting = false)

    # vector spaces and bundles
    begin
        T = MVector{3, Float64}
        TB = TangentBundle(s2)

        pts_tb = [ProductRepr(convert(T, [1.0, 0.0, 0.0]), convert(T, [0.0, -1.0, -1.0])),
                  ProductRepr(convert(T, [0.0, 1.0, 0.0]), convert(T, [2.0, 0.0, 1.0])),
                  ProductRepr(convert(T, [1.0, 0.0, 0.0]), convert(T, [0.0, 2.0, -1.0]))]
        add_manifold(TB, pts_tb, "Tangent bundle of S² using MVectors, ProductRepr";
            test_tangent_vector_broadcasting = false)
    end

    # power manifolds
    begin
        Ms = Sphere(2)
        Ms1 = PowerManifold(Ms, 5)
        Ms2 = PowerManifold(Ms, 5, 7)
        Mr = Manifolds.Rotations(3)
        Mr1 = PowerManifold(Mr, 5)
        Mr2 = PowerManifold(Mr, 5, 7)

        types_s1 = [Array{Float64,2},
                    HybridArray{Tuple{3,StaticArrays.Dynamic()}, Float64, 2}]
        types_s2 = [Array{Float64,3},
                    HybridArray{Tuple{3,StaticArrays.Dynamic(),StaticArrays.Dynamic()}, Float64, 3}]

        types_r1 = [Array{Float64,3},
                    HybridArray{Tuple{3,3,StaticArrays.Dynamic()}, Float64, 3}]
        types_r2 = [Array{Float64,4},
                    HybridArray{Tuple{3,3,StaticArrays.Dynamic(),StaticArrays.Dynamic()}, Float64, 4}]

        retraction_methods = [Manifolds.PowerRetraction(Manifolds.ExponentialRetraction())]
        inverse_retraction_methods = [Manifolds.InversePowerRetraction(Manifolds.LogarithmicInverseRetraction())]

        sphere_dist = Manifolds.uniform_distribution(Ms, @SVector [1.0, 0.0, 0.0])
        power_s1_pt_dist = Manifolds.PowerPointDistribution(Ms1, sphere_dist, randn(Float64, 3, 5))
        power_s2_pt_dist = Manifolds.PowerPointDistribution(Ms2, sphere_dist, randn(Float64, 3, 5, 7))
        sphere_tv_dist = Manifolds.normal_tvector_distribution(Ms, (@MVector [1.0, 0.0, 0.0]), 1.0)
        power_s1_tv_dist = Manifolds.PowerFVectorDistribution(TangentBundleFibers(Ms1), rand(power_s1_pt_dist), sphere_tv_dist)
        power_s2_tv_dist = Manifolds.PowerFVectorDistribution(TangentBundleFibers(Ms2), rand(power_s2_pt_dist), sphere_tv_dist)

        id_rot = @SMatrix [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0]
        rotations_dist = Manifolds.normal_rotation_distribution(Mr, id_rot, 1.0)
        power_r1_pt_dist = Manifolds.PowerPointDistribution(Mr1, rotations_dist, randn(Float64, 3, 3, 5))
        power_r2_pt_dist = Manifolds.PowerPointDistribution(Mr2, rotations_dist, randn(Float64, 3, 3, 5, 7))
        rotations_tv_dist = Manifolds.normal_tvector_distribution(Mr, MMatrix(id_rot), 1.0)
        power_r1_tv_dist = Manifolds.PowerFVectorDistribution(TangentBundleFibers(Mr1), rand(power_r1_pt_dist), rotations_tv_dist)
        power_r2_tv_dist = Manifolds.PowerFVectorDistribution(TangentBundleFibers(Mr2), rand(power_r2_pt_dist), rotations_tv_dist)

        trim(s::String) = s[1:min(length(s), 20)]

        for T in types_s1
            pts1 = [convert(T, rand(power_s1_pt_dist)) for _ in 1:3]
            add_manifold(Ms1, pts1, "power manifold S²^(5,), type $(trim(string(T)))";
                test_tangent_vector_broadcasting = true)
        end
        for T in types_s2
            pts2 = [convert(T, rand(power_s2_pt_dist)) for _ in 1:3]
            add_manifold(Ms2, pts2, "power manifold S²^(5,7), type $(trim(string(T)))";
                test_tangent_vector_broadcasting = true)
        end

        for T in types_r1
            pts1 = [convert(T, rand(power_r1_pt_dist)) for _ in 1:3]
            add_manifold(Mr1, pts1, "power manifold SO(3)^(5,), type $(trim(string(T)))";
                test_tangent_vector_broadcasting = true)
        end
        for T in types_r2
            pts2 = [convert(T, rand(power_r2_pt_dist)) for _ in 1:3]
            add_manifold(Mr2, pts2, "power manifold SO(3)^(5,7), type $(trim(string(T)))";
                test_tangent_vector_broadcasting = true)
        end

    end

end

add_manifold_benchmarks()


text:
export PatchedHomotopy, PatchedHomotopyCache

"""
    PatchedHomotopy(H::AbstractHomotopy, patch, v::PVector)

Augment the homotopy `H` with the given patch `v`. This results in the system `[H(x,t); v ⋅ x - 1]`
"""
struct PatchedHomotopy{H<:AbstractHomotopy,PS<:AbstractAffinePatchState} <: AbstractHomotopy
    homotopy::H
    patch::PS
end

function PatchedHomotopy(H::AbstractHomotopy, p::AbstractAffinePatch, x)
    PatchedHomotopy(H, state(p, x))
end

function Base.size(H::PatchedHomotopy)
    m, n = size(H.homotopy)
    (m + nequations(H.patch), n)
end

struct PatchedHomotopyCache{HC,T} <: AbstractHomotopyCache
    cache::HC
    A::Matrix{T} # intermediate storage of the jacobian
    b::Vector{T} # intermediate storage for the evaluation
end

function cache(ph::PatchedHomotopy, x, t)
    H = ph.homotopy
    c = cache(H, x, t)
    PatchedHomotopyCache(c, jacobian(H, x, t, c), evaluate(H, x, t, c))
end

function evaluate!(u, H::PatchedHomotopy, x, t, c::PatchedHomotopyCache)
    M, N = size(H.homotopy)

    # H(x, t)
    evaluate!(c.b, H.homotopy, x, t, c.cache)
    @inbounds for i = 1:M
        u[i] = c.b[i]
    end
    evaluate!(u, H.patch, x)
    u
end
function jacobian!(U, H::PatchedHomotopy, x::PVector, t, c::PatchedHomotopyCache)
    M, N = size(H.homotopy)
    # J_H(x, t)
    jacobian!(c.A, H.homotopy, x, t, c.cache)
    @inbounds for j = 1:N, i = 1:M
        U[i, j] = c.A[i, j]
    end
    jacobian!(U, H.patch, x)
    U
end

function dt!(u, H::PatchedHomotopy, x::PVector, t, c::PatchedHomotopyCache)
    M, N = size(H.homotopy)
    # [H(x,t); v ⋅ x - 1]/∂t = [∂H(x,t)/∂t; 0]
    dt!(c.b, H.homotopy, x, t, c.cache)
    @inbounds for i = 1:M
        u[i] = c.b[i]
    end
    for i = 1:nequations(H.patch)
        u[M+i] = zero(eltype(u))
    end
    u
end

function evaluate_and_jacobian!(
    u,
    U,
    H::PatchedHomotopy,
    x::PVector,
    t,
    c::PatchedHomotopyCache,
)
    M, N = size(H.homotopy)
    evaluate_and_jacobian!(c.b, c.A, H.homotopy, x, t, c.cache)

    @inbounds for j = 1:N, i = 1:M
        U[i, j] = c.A[i, j]
    end
    jacobian!(U, H.patch, x)

    @inbounds for i = 1:M
        u[i] = c.b[i]
    end
    evaluate!(u, H.patch, x)

    nothing
end

function jacobian_and_dt!(U, u, H::PatchedHomotopy, x::PVector, t, c::PatchedHomotopyCache)
    M, N = size(H.homotopy)
    A, b = c.A, c.b
    jacobian_and_dt!(A, b, H.homotopy, x, t, c.cache)
    # jacobian
    @inbounds for j = 1:N, i = 1:M
        U[i, j] = A[i, j]
    end
    jacobian!(U, H.patch, x)

    # dt
    @inbounds for i = 1:M
        u[i] = b[i]
    end
    for i = 1:nequations(H.patch)
        u[M+i] = zero(eltype(u))
    end

    nothing
end

function evaluate(H::PatchedHomotopy, x, t, c::PatchedHomotopyCache)
    M, N = size(H)
    evaluate!(similar(c.b, M), H, x, t, c)
end
function jacobian(H::PatchedHomotopy, x, t, c::PatchedHomotopyCache)
    M, N = size(H)
    jacobian!(similar(c.A, M, N), H, x, t, c)
end
function dt(H::PatchedHomotopy, x, t, c::PatchedHomotopyCache)
    M, N = size(H)
    dt!(similar(c.b, M), H, x, t, c)
end

basehomotopy(H::PatchedHomotopy) = basehomotopy(H.homotopy)