text:
using ITensors,
      Test

include("2d_classical_ising.jl")

function ctmrg(T::ITensor,
               Clu::ITensor,
               Al::ITensor;
               χmax::Int,nsteps::Int)
  Clu = addtags(Clu,"orig","link")
  Al = addtags(Al,"orig","link")
  for i = 1:nsteps

    ## Get the grown corner transfer matrix (CTM)
    Au = replacetags(Al,"down,link","left,link")
    Au = replacetags(Au,"up,link","right,link")
    Au = replacetags(Au,"left,site","up,site")

    Clu⁽¹⁾ = Clu*Al*Au*T

    ## Diagonalize the grown CTM
    ld = firstind(Clu⁽¹⁾,"link,down")
    sd = firstind(Clu⁽¹⁾,"site,down")
    lr = firstind(Clu⁽¹⁾,"link,right")
    sr = firstind(Clu⁽¹⁾,"site,right")

    Cdr,Ur = eigen(Clu⁽¹⁾, (ld, sd), (lr, sr); ishermitian = true,
                                               maxdim = χmax,
                                               lefttags = "link,down,renorm",
                                               righttags = "link,right,renorm")

    ## The renormalized CTM is the diagonal matrix of eigenvalues
    Clu = replacetags(Cdr,"renorm","orig")
    Clu = replacetags(Clu,"down","up")
    Clu = replacetags(Clu,"right","left")
    Clu = Clu/norm(Clu)

    ## Calculate the renormalized half row transfer matrix (HRTM)
    Ud = replacetags(Ur,"right","down")
    Uu = replacetags(Ud,"down","up")

    Al = Al*Uu*T*Ud
    Al = replacetags(Al,"renorm","orig")
    Al = replacetags(Al,"right,site","left,site")
    Al = Al/norm(Al)
  end
  Clu = removetags(Clu,"orig")
  Al = removetags(Al,"orig")
  return Clu,Al
end

@testset "ctmrg" begin
  # Make Ising model MPO
  β = 1.1*βc
  d = 2
  s = Index(d,"site")
  sl = addtags(s,"left")
  sr = addtags(s,"right")
  su = addtags(s,"up")
  sd = addtags(s,"down")

  T = ising_mpo((sl,sr),(su,sd),β)

  χ0 = 1
  l = Index(χ0,"link")
  ll = addtags(l,"left")
  lu = addtags(l,"up")
  ld = addtags(l,"down")

  # Initial CTM
  Clu = ITensor(lu,ll)
  Clu[1,1] = 1.0

  # Initial HRTM
  Al = ITensor(lu,ld,sl)
  Al[lu(1),ld(1),sl(1)] = 1.0
  Al[lu(1),ld(1),sl(2)] = 0.0

  Clu,Al = ctmrg(T,Clu,Al;χmax=30,nsteps=2000)

  # Normalize corner matrix
  trC⁴ = Clu*mapprime(Clu,0,1,"up")*
         mapprime(Clu,0,1)*mapprime(Clu,0,1,"left")
  Clu = Clu/scalar(trC⁴)^(1/4)

  # Normalize MPS tensor
  trA² = Clu*mapprime(Clu,0,1,"up")*Al*
         mapprime(Al,0,1,"link")*
         mapprime(replacetags(mapprime(Clu,0,1,"up"),"up","down"),0,1,"left")*
         replacetags(mapprime(Clu,0,1,"left"),"up","down")
  Al = Al/sqrt(scalar(trA²))

  ## Get environment tensors for a single site measurement
  Ar = mapprime(replacetags(Al,"left","right","site"),0,1,"link")
  Au = replacetags(replacetags(replacetags(Al,"left","up","site"),
                                              "down","left","link"),
                                              "up","right","link")
  Ad  = mapprime(replacetags(Au,"up","down","site"),0,1,"link")
  Cld = mapprime(replacetags(Clu,"up","down"),0,1,"left")
  Cru = mapprime(replacetags(Clu,"left","right"),0,1,"up")
  Crd = replacetags(mapprime(Cru,0,1,"right"),"up","down")

  ## Calculate partition function per site
  κ = scalar(Clu*Al*Cld*Au*T*Ad*Cru*Ar*Crd)
  @test κ≈exp(-β*ising_free_energy(β))

  ## Calculate magnetization
  Tsz = ising_mpo((sl,sr),(su,sd),β;sz=true)
  m = scalar(Clu*Al*Cld*Au*Tsz*Ad*Cru*Ar*Crd)/κ
  @test abs(m)≈ising_magnetization(β)
end

nothing


text:
using RDatasets, DataArrays, Gadfly

plot(dataset("car", "SLID"), x=:Wages, color=:Language,
     Geom.histogram(position=:dodge))

ext: jl


text:
push!(LOAD_PATH, "G:\\Python Programme\\BALL.jl\\src\\" )
using BALL
using Documenter

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

makedocs(;
    modules=[BALL],
    authors="Daniel Klima <dklima@students.uni-mainz.de>",
    repo="https://github.com/pushingPulling/BALL.jl/blob/{commit}{path}#{line}",
    sitename="BALL.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://pushingPulling.github.io/BALL.jl",
        assets=String[],
    ),
    pages=[
        "BALL Documentation" => "index.md",
        "Manual" =>
            [
                "Getting Started" => "start.md",
                "Interfaces" => "interfaces.md",
                "Iteration" => "iteration.md",
                "Advanced Usage" => "advanced_usage.md"
            ],
        "Submodules" =>
            [
                "Concept" => "concept_page.md",
                "Kernel" => "kernel_page.md",
                "Fileformats" => "fileformats_page.md",
                "Molmec" => "molmec_page.md",
                "Structure" => "structure_page.md",
                "QSAR" => "qsar_page.md",
                "Miscellaneous" => "misc.md"
            ],
        "DevDocs" => "devdocs.md"
    ],
)

deploydocs(;
    repo="github.com/pushingPulling/BALL.jl",
)


text:
@testset "E100" begin
    msgs = lintstr("""
    function f(x, y)
        using Base.Meta
        isexpr(x, :call) ? y : 0
    end
    """)
    @test messageset(msgs) == Set([:E100, :E321])
    @test occursin("using expression must be at top level", msgs[1].message)

    msgs = lintstr("""
    function f(x, y)
        import Lint
        isexpr(x, :call) ? y : 0
    end
    """)
    @test messageset(msgs) == Set([:E100, :E321])
    @test occursin("import expression must be at top level", msgs[1].message)

    msgs = lintstr("""
    function f(x, y)
        export f
        isexpr(x, :call) ? y : 0
    end
    """)
    @test messageset(msgs) == Set([:E100, :E321])
    @test occursin("export expression must be at top level", msgs[1].message)
end


text:
using Compat

if VERSION.minor == 7
  @eval using LinearAlgebra, SparseMatrices
end

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

A = Float64[
  16 4 8;
  4 5 -4;
  8 -4 22;
]
a = [16; 18; -22]

F = chol(A)

@show y = F \ a
println()

@show x = F' \ y
println()

@show As = sparse(A)
println()

@show As\a
println()

B = copy(A)
b = copy(a)

B[2,2] += 1.0e10
b[2] = (5 + 1.0e10) * 5.0

F = chol(B)

@show y = F \ b
println()

@show x = F' \ y
println()

@show Bs = sparse(B)
println()

@show Bs\b
println()

end


text:
# frexp,ldexp,significand,exponent
@test frexp(12.8) == (0.8,4)
@test ldexp(0.8,4) == 12.8
@test significand(12.8) == 1.6
@test exponent(12.8) == 3

# degree-based trig functions
for T = (Float32,Float64)
    for x = -400:40:400
        @test_approx_eq_eps sind(convert(T,x))::T convert(T,sin(pi/180*x)) eps(deg2rad(convert(T,x)))
        @test_approx_eq_eps cosd(convert(T,x))::T convert(T,cos(pi/180*x)) eps(deg2rad(convert(T,x)))
    end
    for x = 0.0:180:720
        @test sind(convert(T,x)) === zero(T)
        @test sind(-convert(T,x)) === -zero(T)
    end

    for x = -3:0.3:3
        @test_approx_eq_eps sinpi(convert(T,x))::T convert(T,sin(pi*x)) eps(pi*convert(T,x))
        @test_approx_eq_eps cospi(convert(T,x))::T convert(T,cos(pi*x)) eps(pi*convert(T,x))
    end
    for x = 0.0:1.0:4.0
        @test sinpi(convert(T,x)) === zero(T)
        @test sinpi(-convert(T,x)) === -zero(T)
    end
end

# check type stability
for T = (Float32,Float64,BigFloat)
    for f = (sind,cosd,sinpi,cospi)
        @test Base.return_types(f,(T,)) == [T]
    end
end

# error functions
@test_approx_eq erf(1) 0.84270079294971486934
@test_approx_eq erfc(1) 0.15729920705028513066
@test_approx_eq erfcx(1) 0.42758357615580700442
@test_approx_eq erfi(1) 1.6504257587975428760
@test_approx_eq erfinv(0.84270079294971486934) 1
@test_approx_eq erfcinv(0.15729920705028513066) 1
@test_approx_eq dawson(1) 0.53807950691276841914
# TODO: complex versions only supported on 64-bit for now
@unix_only if WORD_SIZE==64
    @test_approx_eq erf(1+2im) -0.53664356577856503399-5.0491437034470346695im
    @test_approx_eq erfc(1+2im) 1.5366435657785650340+5.0491437034470346695im
    @test_approx_eq erfcx(1+2im) 0.14023958136627794370-0.22221344017989910261im
    @test_approx_eq erfi(1+2im) -0.011259006028815025076+1.0036063427256517509im
    @test_approx_eq dawson(1+2im) -13.388927316482919244-11.828715103889593303im
end
for x in logspace(-200, -0.01)
    @test_approx_eq_eps erf(erfinv(x)) x 1e-12*x
    @test_approx_eq_eps erf(erfinv(-x)) -x 1e-12*x
    @test_approx_eq_eps erfc(erfcinv(2*x)) 2*x 1e-12*x
    if x > 1e-20
        xf = float32(x)
        @test_approx_eq_eps erf(erfinv(xf)) xf 1e-5*xf
        @test_approx_eq_eps erf(erfinv(-xf)) -xf 1e-5*xf
        @test_approx_eq_eps erfc(erfcinv(2xf)) 2xf 1e-5*xf
    end
end

# airy
@test_approx_eq airy(1.8) 0.0470362168668458052247
@test_approx_eq airyprime(1.8) -0.0685247801186109345638
@test_approx_eq airybi(1.8) 2.595869356743906290060
@test_approx_eq airybiprime(1.8) 2.98554005084659907283
@test_throws Base.Math.AmosException airy(200im)
@test_throws Base.Math.AmosException airybi(200)
z = 1.8 + 1.0im
@test_approx_eq airyx(0, z) airy(0, z) * exp(2/3 * z * sqrt(z))
@test_approx_eq airyx(1, z) airy(1, z) * exp(2/3 * z * sqrt(z))
@test_approx_eq airyx(2, z) airy(2, z) * exp(-abs(real(2/3 * z * sqrt(z))))
@test_approx_eq airyx(3, z) airy(3, z) * exp(-abs(real(2/3 * z * sqrt(z))))

# besselh
true_h133 = 0.30906272225525164362 - 0.53854161610503161800im
@test_approx_eq besselh(3,1,3) true_h133
@test_approx_eq besselh(-3,1,3) -true_h133
@test_approx_eq besselh(3,2,3) conj(true_h133)
@test_approx_eq besselh(-3,2,3) -conj(true_h133)
@test_throws Base.Math.AmosException besselh(1,0)

# besseli
true_i33 = 0.95975362949600785698
@test_approx_eq besseli(3,3) true_i33
@test_approx_eq besseli(-3,3) true_i33
@test_approx_eq besseli(3,-3) -true_i33
@test_approx_eq besseli(-3,-3) -true_i33
@test_throws Base.Math.AmosException besseli(1,1000)

# besselj
@test besselj(0,0) == 1
for i = 1:5
    @test besselj(i,0) == 0
    @test besselj(-i,0) == 0
end

j33 = besselj(3,3.)
@test besselj(3,3) == j33
@test besselj(-3,-3) == j33
@test besselj(-3,3) == -j33
@test besselj(3,-3) == -j33

j43 = besselj(4,3.)
@test besselj(4,3) == j43
@test besselj(-4,-3) == j43
@test besselj(-4,3) == j43
@test besselj(4,-3) == j43

@test_approx_eq j33 0.30906272225525164362
@test_approx_eq j43 0.13203418392461221033
@test_throws DomainError    besselj(0.1, -0.4)
@test_approx_eq besselj(0.1, complex(-0.4)) 0.820421842809028916 + 0.266571215948350899im
@test_approx_eq besselj(3.2, 1.3+0.6im) 0.01135309305831220201 + 0.03927719044393515275im
@test_approx_eq besselj(1, 3im) 3.953370217402609396im
@test_throws Base.Math.AmosException besselj(20,1000im)

# besselk
true_k33 = 0.12217037575718356792
@test_approx_eq besselk(3,3) true_k33
@test_approx_eq besselk(-3,3) true_k33
true_k3m3 = -0.1221703757571835679 - 3.0151549516807985776im
@test_throws DomainError besselk(3,-3)
@test_approx_eq besselk(3,complex(-3)) true_k3m3
@test_approx_eq besselk(-3,complex(-3)) true_k3m3
@test_throws Base.Math.AmosException besselk(200,0.01)
# issue #6564
@test besselk(1.0,0.0) == Inf

# bessely
y33 = bessely(3,3.)
@test bessely(3,3) == y33
@test_approx_eq bessely(-3,3) -y33
@test_approx_eq y33 -0.53854161610503161800
@test_throws DomainError bessely(3,-3)
@test_approx_eq bessely(3,complex(-3)) 0.53854161610503161800 - 0.61812544451050328724im
@test_throws Base.Math.AmosException bessely(200.5,0.1)

# issue #6653
for f in (besselj,bessely,besseli,besselk,hankelh1,hankelh2)
    @test_approx_eq f(0,1) f(0,complex128(1))
    @test_approx_eq f(0,1) f(0,complex64(1))
end

# scaled bessel[ijky] and hankelh[12]
for x in (1.0, 0.0, -1.0), y in (1.0, 0.0, -1.0), nu in (1.0, 0.0, -1.0)
    z = complex128(x + y * im)
    z == zero(z) || @test_approx_eq hankelh1x(nu, z) hankelh1(nu, z) * exp(-z * im)
    z == zero(z) || @test_approx_eq hankelh2x(nu, z) hankelh2(nu, z) * exp(z * im)
    (nu < 0 && z == zero(z)) || @test_approx_eq besselix(nu, z) besseli(nu, z) * exp(-abs(real(z)))
    (nu < 0 && z == zero(z)) || @test_approx_eq besseljx(nu, z) besselj(nu, z) * exp(-abs(imag(z)))
    z == zero(z) || @test_approx_eq besselkx(nu, z) besselk(nu, z) * exp(z)
    z == zero(z) || @test_approx_eq besselyx(nu, z) bessely(nu, z) * exp(-abs(imag(z)))
end
@test_throws Base.Math.AmosException hankelh1x(1, 0)
@test_throws Base.Math.AmosException hankelh2x(1, 0)
@test_throws Base.Math.AmosException besselix(-1, 0)
@test_throws Base.Math.AmosException besseljx(-1, 0)
@test besselkx(1, 0) == Inf
@test_throws Base.Math.AmosException besselyx(1, 0)

# beta, lbeta
@test_approx_eq beta(3/2,7/2) 5π/128
@test_approx_eq beta(3,5) 1/105
@test_approx_eq lbeta(5,4) log(beta(5,4))
@test_approx_eq beta(5,4) beta(4,5)
@test_approx_eq beta(-1/2, 3) -16/3
@test_approx_eq lbeta(-1/2, 3) log(16/3)

# gamma, lgamma (complex argument)
@test gamma(1:25) == gamma(Float64[1:25])
@test_approx_eq gamma(1/2) sqrt(π)
@test_approx_eq gamma(-1/2) -2sqrt(π)
@test_approx_eq lgamma(-1/2) log(abs(gamma(-1/2)))
@test_approx_eq lgamma(1.4+3.7im) -3.7094025330996841898 + 2.4568090502768651184im
@test_approx_eq lgamma(1.4+3.7im) log(gamma(1.4+3.7im))

# digamma
for elty in (Float32, Float64)

    @test_approx_eq digamma(convert(elty, 9)) convert(elty, 2.140641477955609996536345)
    @test_approx_eq digamma(convert(elty, 2.5)) convert(elty, 0.7031566406452431872257)
    @test_approx_eq digamma(convert(elty, 0.1)) convert(elty, -10.42375494041107679516822)
    @test_approx_eq digamma(convert(elty, 7e-4)) convert(elty, -1429.147493371120205005198)
    @test_approx_eq digamma(convert(elty, 7e-5)) convert(elty, -14286.29138623969227538398)
    @test_approx_eq digamma(convert(elty, 7e-6)) convert(elty, -142857.7200612932791081972)
    @test_approx_eq digamma(convert(elty, 2e-6)) convert(elty, -500000.5772123750382073831)
    @test_approx_eq digamma(convert(elty, 1e-6)) convert(elty, -1000000.577214019968668068)
    @test_approx_eq digamma(convert(elty, 7e-7)) convert(elty, -1428572.005785942019703646)
    @test_approx_eq digamma(convert(elty, -0.5)) convert(elty, .03648997397857652055902367)
    @test_approx_eq digamma(convert(elty, -1.1)) convert(elty,  10.15416395914385769902271)

    @test_approx_eq digamma(convert(elty, 0.1)) convert(elty, -10.42375494041108)
    @test_approx_eq digamma(convert(elty, 1/2)) convert(elty, -γ - log(4))
    @test_approx_eq digamma(convert(elty, 1)) convert(elty, -γ)
    @test_approx_eq digamma(convert(elty, 2)) convert(elty, 1 - γ)
    @test_approx_eq digamma(convert(elty, 3)) convert(elty, 3/2 - γ)
    @test_approx_eq digamma(convert(elty, 4)) convert(elty, 11/6 - γ)
    @test_approx_eq digamma(convert(elty, 5)) convert(elty, 25/12 - γ)
    @test_approx_eq digamma(convert(elty, 10)) convert(elty, 7129/2520 - γ)
end

# trigamma
for elty in (Float32, Float64)
    @test_approx_eq trigamma(convert(elty, 0.1)) convert(elty, 101.433299150792758817)
    @test_approx_eq trigamma(convert(elty, 1/2)) convert(elty, π^2/2)
    @test_approx_eq trigamma(convert(elty, 1)) convert(elty, π^2/6)
    @test_approx_eq trigamma(convert(elty, 2)) convert(elty, π^2/6 - 1)
    @test_approx_eq trigamma(convert(elty, 3)) convert(elty, π^2/6 - 5/4)
    @test_approx_eq trigamma(convert(elty, 4)) convert(elty, π^2/6 - 49/36)
    @test_approx_eq trigamma(convert(elty, 5)) convert(elty, π^2/6 - 205/144)
    @test_approx_eq trigamma(convert(elty, 10)) convert(elty, π^2/6 - 9778141/6350400)
end

# invdigamma
for elty in (Float32, Float64)
    for val in [0.001, 0.01, 0.1, 1.0, 10.0]
        @test abs(invdigamma(digamma(convert(elty, val))) - convert(elty, val)) < 1e-8
    end
end

@test_approx_eq polygamma(20, 7.) -4.644616027240543262561198814998587152547

# eta, zeta
@test_approx_eq eta(1) log(2)
@test_approx_eq eta(2) pi^2/12
@test_approx_eq zeta(0) -0.5
@test_approx_eq zeta(2) pi^2/6
@test_approx_eq zeta(4) pi^4/90

# quadgk
@test_approx_eq quadgk(cos, 0,0.7,1)[1] sin(1)
@test_approx_eq quadgk(x -> exp(im*x), 0,0.7,1)[1] (exp(1im)-1)/im
@test_approx_eq quadgk(x -> exp(im*x), 0,1im)[1] -1im*expm1(-1)
@test_approx_eq_eps quadgk(cos, 0,BigFloat(1),order=40)[1] sin(BigFloat(1)) 1000*eps(BigFloat)
@test_approx_eq quadgk(x -> exp(-x), 0,0.7,Inf)[1] 1.0
@test_approx_eq quadgk(x -> exp(x), -Inf,0)[1] 1.0
@test_approx_eq quadgk(x -> exp(-x^2), -Inf,Inf)[1] sqrt(pi)
@test_approx_eq quadgk(x -> [exp(-x), exp(-2x)], 0, Inf)[1] [1,0.5]
@test_approx_eq quadgk(cos, 0,0.7,1, norm=abs)[1] sin(1)

# Ensure subnormal flags functions don't segfault
@test any(ccall("jl_zero_subnormals", Uint8, (Uint8,), 1) .== [0x00 0x01])
@test any(ccall("jl_zero_subnormals", Uint8, (Uint8,), 0) .== [0x00 0x01])

# isqrt (issue #4884)
@test isqrt(9223372030926249000) == 3037000498
@test isqrt(typemax(Int128)) == int128("13043817825332782212")
@test isqrt(int128(typemax(Int64))^2-1) == 9223372036854775806
@test isqrt(0) == 0
for i = 1:1000
    n = rand(Uint128)
    s = isqrt(n)
    @test s*s <= n
    @test (s+1)*(s+1) > n
    n = rand(Uint64)
    s = isqrt(n)
    @test s*s <= n
    @test (s+1)*(s+1) > n
end

# useful test functions for relative error
err(z, x) = abs(z - x) / abs(x)
errc(z, x) = max(err(real(z),real(x)), err(imag(z),imag(x)))

for x in -10.2:0.3456:50
    @test 1e-12 > err(digamma(x+0im), digamma(x))
end

# digamma, trigamma, polygamma & zeta test cases (compared to Wolfram Alpha)
@test 1e-13 > err(digamma(7+0im), 1.872784335098467139393487909917597568957840664060076401194232)
@test 1e-13 > errc(digamma(7im), 1.94761433458434866917623737015561385331974500663251349960124 + 1.642224898223468048051567761191050945700191089100087841536im)
@test 1e-13 > errc(digamma(-3.2+0.1im), 4.65022505497781398615943030397508454861261537905047116427511+2.32676364843128349629415011622322040021960602904363963042380im)
@test 1e-13 > err(trigamma(8+0im), 0.133137014694031425134546685920401606452509991909746283540546)
@test 1e-13 > errc(trigamma(8im), -0.0078125000000000000029194973110119898029284994355721719150 - 0.12467345030312762782439017882063360876391046513966063947im)
@test 1e-13 > errc(trigamma(-3.2+0.1im), 15.2073506449733631753218003030676132587307964766963426965699+15.7081038855113567966903832015076316497656334265029416039199im)
@test 1e-13 > err(polygamma(2, 8.1+0im), -0.01723882695611191078960494454602091934457319791968308929600)
@test 1e-13 > errc(polygamma(30, 8.1+2im), -2722.8895150799704384107961215752996280795801958784600407589+6935.8508929338093162407666304759101854270641674671634631058im)
@test 1e-13 > errc(polygamma(3, 2.1+1im), 0.00083328137020421819513475400319288216246978855356531898998-0.27776110819632285785222411186352713789967528250214937861im)
@test 1e-11 > err(polygamma(3, -4.2 + 2im),-0.0037752884324358856340054736472407163991189965406070325067-0.018937868838708874282432870292420046797798431078848805822im)
@test 1e-13 > err(polygamma(13, 5.2 - 2im), 0.08087519202975913804697004241042171828113370070289754772448-0.2300264043021038366901951197725318713469156789541415899307im)
@test 1e-11 > err(polygamma(123, -47.2 + 0im), 5.7111648667225422758966364116222590509254011308116701029e291)
@test 1e-13 > errc(zeta(4.1+0.3im, -3.2+0.1im), -461.95403678374488506025596495576748255121001107881278765917+926.02552636148651929560277856510991293536052745360005500774im)
@test 1e-13 > errc(zeta(4.1+0.3im, 3.2+0.1im), 0.0121197525131633219465301571139288562254218365173899270675-0.00687228692565614267981577154948499247518236888933925740902im)
@test 1e-13 > errc(zeta(4.1, 3.2+0.1im),0.0137637451187986846516125754047084829556100290057521276517-0.00152194599531628234517456529686769063828217532350810111482im)
@test 1e-12 > errc(zeta(1.0001, -4.5e2+3.2im), 9993.89099199843392251301993718413132850540848778561412270571-3.13257480938495907945892330398176989805350557816701044268548im)
@test_throws DomainError zeta(3.1,-4.2)
@test 1e-13 > errc(zeta(3.1,-4.2+0im), -138.06320182025311080661516120845508778572835942189570145952+45.586579397698817209431034568162819207622092308850063038062im)
@test 1e-15 > errc(zeta(3.1+0im,-4.2), zeta(3.1,-4.2+0im))
@test 1e-13 > errc(zeta(3.1,4.2), 0.029938344862645948405021260567725078588893266227472565010234)
@test 1e-13 > err(zeta(27, 3.1), 5.413318813037879056337862215066960774064332961282599376e-14)
@test 1e-13 > err(zeta(27, 2), 7.4507117898354294919810041706041194547190318825658299932e-9)
@test 1e-12 > err(zeta(27, -105.3), -1.311372652244914148556295810515903234635727465138859603e14)
@test polygamma(4, -3.1+Inf*im) == polygamma(4, 3.1+Inf*im) == 0
@test polygamma(4, -0.0) == Inf == -polygamma(4, +0.0)
@test zeta(4, +0.0) == Inf == zeta(4, -0.0)
@test zeta(5, +0.0) == Inf == -zeta(5, -0.0)
@test isa([digamma(x) for x in [1.0]], Vector{Float64})
@test isa([trigamma(x) for x in [1.0]], Vector{Float64})
@test isa([polygamma(3,x) for x in [1.0]], Vector{Float64})
@test 1e-13 > errc(zeta(2 + 1im, -1.1), zeta(2 + 1im, -1.1+0im))
@test 1e-13 > errc(zeta(2 + 1im, -1.1), -1525.8095173321060982383023516086563741006869909580583246557 + 1719.4753293650912305811325486980742946107143330321249869576im)

@test @evalpoly(2,3,4,5,6) == 3+2*(4+2*(5+2*6)) == @evalpoly(2+0im,3,4,5,6)

@test 1e-14 > err(eta(1+1e-9), 0.693147180719814213126976796937244130533478392539154928250926)
@test 1e-14 > err(eta(1+5e-3), 0.693945708117842473436705502427198307157819636785324430166786)
@test 1e-13 > err(eta(1+7.1e-3), 0.694280602623782381522315484518617968911346216413679911124758)
@test 1e-13 > err(eta(1+8.1e-3), 0.694439974969407464789106040237272613286958025383030083792151)
@test 1e-13 > err(eta(1 - 2.1e-3 + 2e-3 * im), 0.69281144248566007063525513903467244218447562492555491581+0.00032001240133205689782368277733081683574922990400416791019im)
@test 1e-13 > err(eta(1 + 5e-3 + 5e-3 * im), 0.69394652468453741050544512825906295778565788963009705146+0.00079771059614865948716292388790427833787298296229354721960im)
@test 1e-12 > errc(zeta(1e-3+1e-3im), -0.5009189365276307665899456585255302329444338284981610162-0.0009209468912269622649423786878087494828441941303691216750im)
@test 1e-13 > errc(zeta(1e-4 + 2e-4im), -0.5000918637469642920007659467492165281457662206388959645-0.0001838278317660822408234942825686513084009527096442173056im)

# Issue #7169: (TODO: better accuracy should be possible?)
@test 1e-9 > errc(zeta(0 + 99.69im), 4.67192766128949471267133846066040655597942700322077493021802+3.89448062985266025394674304029984849370377607524207984092848im)
@test 1e-12 > errc(zeta(3 + 99.69im), 1.09996958148566565003471336713642736202442134876588828500-0.00948220959478852115901654819402390826992494044787958181148im)
@test 1e-9 > errc(zeta(-3 + 99.69im), 10332.6267578711852982128675093428012860119184786399673520976+13212.8641740351391796168658602382583730208014957452167440726im)
@test 1e-13 > errc(zeta(2 + 99.69im, 1.3), 0.41617652544777996034143623540420694985469543821307918291931-0.74199610821536326325073784018327392143031681111201859489991im)

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
#####
##### Terms in the turbulent kinetic energy equation, all at cell centers
#####

@inline ϕ²(i, j, k, grid, ϕ) = ϕ(i, j, k, grid)^2

# Temporary way to get the vertical diffusivity for the TKE equation terms...
# Assumes that the vertical diffusivity is dominated by the CATKE contribution.
@inline shear_production(i, j, k, grid, closure, velocities, diffusivities) = zero(eltype(grid))
@inline buoyancy_flux(i, j, k, grid, closure, velocities, tracers, buoyancy, diffusivities) = zero(eltype(grid))

# Unlike the above, this fallback for dissipation is generically correct (we only want to compute dissipation once)
@inline dissipation(i, j, k, grid, closure, velocities, tracers, buoyancy, clock, tracer_bcs) = zero(eltype(grid))

@inline function shear_production(i, j, k, grid, closure::FlavorOfCATKE, velocities, diffusivities)
    ∂z_u² = ℑxzᶜᵃᶜ(i, j, k, grid, ϕ², ∂zᶠᶜᶠ, velocities.u)
    ∂z_v² = ℑyzᵃᶜᶜ(i, j, k, grid, ϕ², ∂zᶜᶠᶠ, velocities.v)
    νᶻ = ℑzᵃᵃᶜ(i, j, k, grid, diffusivities.Kᵘ)
    return νᶻ * (∂z_u² + ∂z_v²)
end

@inline function buoyancy_flux(i, j, k, grid, closure::FlavorOfCATKE, velocities, tracers, buoyancy, diffusivities)
    κᶻ = ℑzᵃᵃᶜ(i, j, k, grid, diffusivities.Kᶜ)
    N² = ℑzᵃᵃᶜ(i, j, k, grid, ∂z_b, buoyancy, tracers)
    return - κᶻ * N²
end

const VITD = VerticallyImplicitTimeDiscretization

@inline dissipation(i, j, k, grid, closure::FlavorOfCATKE{<:VITD}, args...) = zero(eltype(grid))

@inline function implicit_dissipation_coefficient(i, j, k, grid, closure::FlavorOfCATKE{<:VITD}, velocities, tracers, buoyancy, clock, tracer_bcs)
    e = tracers.e
    FT = eltype(grid)
    @inbounds e⁺ = abs(e[i, j, k])

    ℓ = ℑzᵃᵃᶜ(i, j, k, grid, TKE_mixing_lengthᶜᶜᶠ, closure, velocities, tracers, buoyancy, clock, tracer_bcs)
    Cᴰ = closure.Cᴰ

    # Note:
    #   * because   ∂t e + ⋯ = ⋯ + L e = ⋯ - ϵ,
    #
    #   * then      L e = - ϵ
    #                   = - Cᴰ e³² / ℓ
    #
    #   * and thus    L = - Cᴰ √e / ℓ

    return - Cᴰ * sqrt(e⁺) / ℓ
end

# Fallbacks for explicit time discretization
@inline dissipation(i, j, k, grid, closure::FlavorOfCATKE, velocities, tracers, args...) =
    @inbounds - tracers.e[i, j, k] * implicit_dissipation_coefficient(i, j, k, grid, closure::FlavorOfCATKE, args...)

@inline implicit_dissipation_coefficient(i, j, k, grid, closure::FlavorOfCATKE, args...) = zero(eltype(grid))

#####
##### For closure tuples...
#####

@inline shear_production(i, j, k, grid, closures::Tuple{<:Any}, velocities, diffusivities) =
    shear_production(i, j, k, grid, closures[1], velocities, diffusivities[1])

@inline shear_production(i, j, k, grid, closures::Tuple{<:Any, <:Any}, velocities, diffusivities) =
    shear_production(i, j, k, grid, closures[1], velocities, diffusivities[1]) +
    shear_production(i, j, k, grid, closures[2], velocities, diffusivities[2])

@inline shear_production(i, j, k, grid, closures::Tuple{<:Any, <:Any, <:Any}, velocities, diffusivities) =
    shear_production(i, j, k, grid, closures[1], velocities, diffusivities[1]) +
    shear_production(i, j, k, grid, closures[2], velocities, diffusivities[2]) +
    shear_production(i, j, k, grid, closures[3], velocities, diffusivities[3])

@inline buoyancy_flux(i, j, k, grid, closures::Tuple{<:Any}, velocities, tracers, buoyancy, diffusivities) =
    buoyancy_flux(i, j, k, grid, closures[1], velocities, diffusivities[1])

@inline buoyancy_flux(i, j, k, grid, closures::Tuple{<:Any, <:Any}, velocities, tracers, buoyancy, diffusivities) =
    buoyancy_flux(i, j, k, grid, closures[1], velocities, tracers, buoyancy, diffusivities[1]) +
    buoyancy_flux(i, j, k, grid, closures[2], velocities, tracers, buoyancy, diffusivities[2])

@inline buoyancy_flux(i, j, k, grid, closures::Tuple{<:Any, <:Any, <:Any}, velocities, tracers, buoyancy, diffusivities) =
    buoyancy_flux(i, j, k, grid, closures[1], velocities, tracers, buoyancy, diffusivities[1]) +
    buoyancy_flux(i, j, k, grid, closures[2], velocities, tracers, buoyancy, diffusivities[2]) +
    buoyancy_flux(i, j, k, grid, closures[3], velocities, tracers, buoyancy, diffusivities[3])

@inline dissipation(i, j, k, grid, closures::Tuple{<:Any}, velocities, tracers, buoyancy, diffusivities) =
    dissipation(i, j, k, grid, closures[1], velocities, diffusivities[1])

@inline dissipation(i, j, k, grid, closures::Tuple{<:Any, <:Any}, velocities, tracers, buoyancy, diffusivities) =
    dissipation(i, j, k, grid, closures[1], velocities, tracers, buoyancy, diffusivities[1]) +
    dissipation(i, j, k, grid, closures[2], velocities, tracers, buoyancy, diffusivities[2])

@inline dissipation(i, j, k, grid, closures::Tuple{<:Any, <:Any, <:Any}, velocities, tracers, buoyancy, diffusivities) =
    dissipation(i, j, k, grid, closures[1], velocities, tracers, buoyancy, diffusivities[1]) +
    dissipation(i, j, k, grid, closures[2], velocities, tracers, buoyancy, diffusivities[2]) +
    dissipation(i, j, k, grid, closures[3], velocities, tracers, buoyancy, diffusivities[3])


text:
#################################################
# Rendering
# This controls how failures are displayed
abstract type RenderMode end
struct Diff <: RenderMode end

abstract type BeforeAfter <: RenderMode end
struct BeforeAfterLimited <: BeforeAfter end
struct BeforeAfterFull <: BeforeAfter end
struct BeforeAfterImage <: BeforeAfter end

render_item(::RenderMode, item) = println(item)
function render_item(::BeforeAfterLimited, item)
    show(IOContext(stdout, :limit=>true, :displaysize=>(20,80)), "text/plain", item)
    println()
end
function render_item(::BeforeAfterImage, item)
    str_item = @withcolor ImageInTerminal.encodeimg(ImageInTerminal.SmallBlocks(), ImageInTerminal.TermColor256(), item, 20, 40)[1]
    println("eltype: ", eltype(item))
    println("size: ", map(length, axes(item)))
    println("thumbnail:")
    println.(str_item)
end

## 2 arg form render for comparing
function render(mode::BeforeAfter, reference, actual)
    println("- REFERENCE -------------------")
    render_item(mode, reference)
    println("-------------------------------")
    println("- ACTUAL ----------------------")
    render_item(mode, actual)
    println("-------------------------------")
end
function render(::Diff, reference, actual)
    println("- DIFF ------------------------")
    @withcolor println(deepdiff(reference, actual))
    println("-------------------------------")
end

## 1 arg form render for new content
function render(mode::RenderMode, actual)
    println("- NEW CONTENT -----------------")
    render_item(mode, actual)
    println("-------------------------------")
end

#######################################
# IO
# Right now this basically just extends FileIO to support some things as text files

const TextFile = Union{File{format"TXT"}, File{format"SHA256"}}

function loadfile(T, file::File)
    T(load(file)) # Fallback to FileIO
end

function loadfile(T, file::TextFile)
    read(file.filename,String)
end

function savefile(file::File, content)
    save(file, content) # Fallback to FileIO
end

function savefile(file::TextFile, content)
    write(file.filename, content)
end

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

# Final function
# all other functions should hit one of this eventually
# Which handles the actual testing and user prompting

function _test_reference(rendermode, file::File, actual)
    _test_reference(isequal, rendermode, file, actual)
end

function _test_reference(equiv, rendermode, file::File, actual::T) where T
    path = file.filename
    dir, filename = splitdir(path)
    try
        reference = loadfile(T, file)
        if equiv(reference, actual)
            @test true # to increase test counter if reached
        else # test failed
            println("Test for \"$filename\" failed.")
            render(rendermode, reference, actual)
            if isinteractive()
                if input_bool("Replace reference with actual result (path: $path)?")
                    savefile(file, actual)
                    warn("Please run the tests again for any changes to take effect")
                else
                    @test false
                end
            else
                error("You need to run the tests interactively with 'include(\"test/runtests.jl\")' to update reference images")
            end
        end
    catch ex
        if ex isa SystemError || # File doesn't exist
            (Sys.isapple() && ex isa MethodError) ||  # MethodError is for OSX for some reason
            (ex isa ErrorException && startswith(ex.msg, "unable to open"))

            println("Reference file for \"$filename\" does not exist.")
            render(rendermode, actual)
            if isinteractive()
                if input_bool("Create reference file with above content (path: $path)?")
                    mkpath(dir)
                    savefile(file, actual)
                    warn("Please run the tests again for any changes to take effect")
                else
                    @test false
                end
            else
                error("You need to run the tests interactively with 'include(\"test/runtests.jl\")' to create new reference images")
            end
        else
            rethrow(ex)
        end
    end
end


text:
"""
    write_file!

This function writes a vector to a file.

# Arguments
- `v`: vector to write
- `name`: file name
- `dir`: optional argument to give path to the file
"""
function write_file!(v::Vector{T}, name::String, dir = ".") where T <: Number
    if parallel.is_root()
        open("$dir/$name", "w") do io
            for i in eachindex(v)
                println(io, v[i])
            end
        end
    end
end

function write_file!(v::Vector{Dict{Int,T}}, name::String, dir = ".") where T <: Number
    open("$dir/$(name)_$(parallel.myid()).txt", "w") do io
        ids = keys(v[1])
        j = 1
        for id in ids
            if j > 1
                print(io, ",")
            end
            print(io, id)
            j += 1
        end
        print(io, "\n")

        for i in eachindex(v)
            j = 1
            for id in ids
                if j > 1
                    print(io, ",")
                end
                print(io, v[i][id])
                j += 1
            end
            print(io, "\n")
        end
    end
end


text:
function BubbleSort(x::AbstractVector)

   for i=length(x):-1:1
        for j=2:i
            if x[j-1] > x[j]
                temp = x[j-1]
                x[j-1] = x[j]
                x[j] =temp
            end
        end
    end

    return x
end


text:
import StateSpaceReconstruction:
    DelaunayTriangulation,
    Simplex,
    centroid,
    radius,
    MutableSSimplex

import StaticArrays:
    MMatrix,
    MVector

import Statistics:
    transpose!

import InplaceOps:
    @!

export SuperFast

"""
    simplices_within_circumsphere!(inds::Vector{Int}, idx::Int,
        simplices, image_simplices,
        radii_simplices, radii_image_simplices,
        centroids_simplices, centroids_image_simplices)

Find the indices of the simplices whose centroid lies within the circumsphere of the image
of simplex #`idx`. Store the indices in the pre-allocated vector `inds` (the latter part
of the vector might contain zeros if some simplices do not lie within the circumsphere,
so keep track of this outside the function).
"""
function simplices_within_circumsphere!(inds::Vector{Int}, N, idx::Int,
        simplices, image_simplices,
        radii_simplices, radii_image_simplices,
        centroids_simplices, centroids_image_simplices)

    # Because we're re-using `inds`, re-initialize it with zeros
    inds .= 0

    # We've found zero simplices within the circumsphere to begin with.
    N[] = 0

	n_simplices = length(simplices)

	@inbounds for i = 1:n_simplices
        x = transpose(centroids_image_simplices[idx] .- centroids_simplices[i]) * (centroids_image_simplices[idx] .- centroids_simplices[i])
	    δ = x .- ((radii_image_simplices[idx] + radii_simplices[i])^2)

	    if δ[1] < 0
            N[] += 1
            inds[N[]] = i
	    end
	end

end

"""
Get all simplices in the triangulation.
"""
function get_simplices(pts, DT::DelaunayTriangulation)
    dim = length(DT[1]) - 1

    [Simplex(pts[DT[i]]) for i = 1:length(DT)]
end

function mdet_fm_in_x_view!(M, v, n, startstops)
    @inbounds for i = 0:(n - 1)
        M[startstops[i+1]] = view(v, startstops[i+1])
    end
    @fastmath det(M)
end

"""
Fill the 2nd to last row in the i-th column of `prealloc_matrix` with the
2nd to last row in the `i`-th column of `toprow_simplex` (replacing the `i`th vertex).
"""
function fill_vertex(prealloc_matrix::Array{Float64, 2}, pt,
        startstops::Vector{UnitRange{Int64}}, i::Int)
    @show pt, startstops[i]
    prealloc_matrix[startstops[i]] .= pt
end

"""
Fill the 2nd to last row in the i-th column of `prealloc_matrix` with the
2nd to last row in the `i`-th column of `toprow_simplex` (replacing the `i`th vertex).
"""
function refill_vertex(prealloc_matrix::Array{Float64, 2}, toprow_simplex::Array{Float64, 2},
        startstops::Vector{UnitRange{Int64}}, i::Int)

    prealloc_matrix[startstops[i]] .= toprow_simplex[startstops[i]]
end

"""

## Arguments
- **`pt::Vector{Float64}`**: The point
- **`toprow_simplex::Array{Float64, 2}`**: A matrix representation of the simplex where each
    column represents a vertex, but with a row of ones appended on top of it, so that the
    total dimension is `(dim + 1)-by-(dim + 1)`. This is needed when computing determinants.
- **`dim`**: The dimension of the space.
"""
function iscontained_superoptim2!(pt, toprow_simplex, prealloc_matrix, dim, signs, startstops, found::Array{Float64,0})

    # Replace first vertex with the point
    fill_vertex(prealloc_matrix, pt, startstops, 1)

    signs[1] = sign(mdet_fm_in_x_view!(prealloc_matrix, toprow_simplex, dim + 1, startstops))

    refill_vertex(prealloc_matrix, toprow_simplex, startstops, 1)

    for k = 2:dim
        idxs = ((dim + 1)*(k - 1)+1):((dim + 1)*(k - 1)+ 1 + dim - 1)
        fill_vertex(prealloc_matrix, toprow_simplex, startstops, 1)
    end
    found[] = true
end

function SuperFast(inv_pts; n_sample_pts::Int = 200, sample_randomly::Bool = false)

    # Triangulate all points but the last (we need the last point
    # for the forward projection of the simplices)
    DT = DelaunayTriangulation(inv_pts[1:(end - 1)])

    # Some constants used throughout the function
    n_simplices = length(DT)
    dim = length(DT[1]) - 1

    #=
    # Prepare memory-efficient representations of the simplices
    =#
    simplices = get_simplices(inv_pts[1:(end- 1)], DT)
    image_simplices = get_simplices(inv_pts[2:end], DT)

    radii_simplices = [radius(simplex) for simplex in simplices]
    radii_image_simplices = [radius(simplex) for simplex in image_simplices]
    centroids_simplices = [centroid(simplex) for simplex in simplices]
    centroids_image_simplices = [centroid(simplex) for simplex in image_simplices]

    # Matrix representations of the simplices
    matrix_image_simplices = [hcat(image_simplices[i]...,) for i = 1:n_simplices]

    #=
    Compute the convex expansions coefficients needed to generate points
    within a simplex. Also, update number of points if a shape-preserving
    simplex subdivision algorithm was used (in that case, we get a
    *minimum* of `n_sample_pts` set of coefficients.
    =#
    convex_coeffs = subsample_coeffs(dim, n_sample_pts, sample_randomly)
    n_coeffs::Int = size(convex_coeffs, 2)
    convex_coeffs = [reshape(convex_coeffs[:, i], dim + 1, 1) for i = 1:n_coeffs]

    # The Markov matrix
    M = zeros(Float64, n_simplices, n_simplices)

    inds = zeros(Int, n_simplices)
    n_simplices_within_circumsphere = Array{Int, 0}(undef)

    # Indices to fill when computing whether a point lies inside a simplex or not.
    # The top row will always be a row of ones, so we need to know which other
    # indices to fill.
    starts = [((dim + 1)*i)+2 for i = 0:dim]
    stops = starts .+ (dim - 1)
    startstops = [starts[i]:stops[i] for i = 1:length(starts)]

    # Matrix to re-use when computing whether a point lies inside a simplex.
    prealloc_matrix = zeros(Float64, dim + 1, dim + 1)#zeros(MMatrix{dim + 1, dim + 1})
    toprow_inds = [((dim + 1)*i) + 1 for i = 0:dim]
    [prealloc_matrix[i] = 1.0 for i in toprow_inds]

    # Signs vector when checking convex coefficients
    signs = zeros(Int, dim + 1)

    # Matrix representations of the simplices
    toprow_simplices = [[transpose(ones(dim + 1)); hcat(simplices[i]...,)] for i = 1:n_simplices]

    pt = reshape(zeros(Float64, dim), dim, 1)

    transfermatrix = zeros(Float64, n_simplices, n_simplices)

    found = Array{Float64, 0}(undef)

    for im_idx in 1:n_simplices
        # Find the indices of the simplices who potentially intersect
        # with the image of simplex #i
        simplices_within_circumsphere!(inds, n_simplices_within_circumsphere, im_idx,
               simplices, image_simplices,
               radii_simplices, radii_image_simplices,
               centroids_simplices, centroids_image_simplices)

        # Generate points within the image of simplex #i and check
        # what fraction of those points fall into the simplices.
        # We approximate the intersecting volume between the
        # image simplex and the individual simplices by those
        # fractions.

        # View into the image simplex at index `im_idx`
        im_simplex = view(matrix_image_simplices[im_idx], :, :)

        inds_potentially_intersecting = [inds[i] for i = 1:n_simplices_within_circumsphere[]]

        @show inds_potentially_intersecting
        for k in inds_potentially_intersecting
            # Should be a (dim-by(dim+1) * (dim + 1) multiplication -> (dim, 1) vector
            @! pt = im_simplex * convex_coeffs[k]

            # Keep track of whether the simplex is found
            found[] = false

            # Now, loop over the simplices and see if the point we generated inside the
            # image simplex lies inside any of them, if so, update the counts.
            for idx in n_simplices_within_circumsphere
                iscontained_superoptim2!(pt, toprow_simplices[idx], prealloc_matrix, dim, signs, startstops, found)

                #if found[] == true
                #    break
                #end
            end
        end
    end
    # Need to normalise, because all we have up until now is counts
    # of how many points inside the image simplex falls into
    # the simplices.
    #return transpose(transfermatrix) ./ n_coeffs
end


text:
module HED

# package code goes here
include("detect.jl")

end # module


text:
using CXX

@osx_only begin
    const qtlibdir = "/Users/kfischer/Projects/qt-everywhere-opensource-src-5.3.1/qtbase/~/usr/lib/"
    const QtCore = joinpath(qtlibdir,"QtCore.framework/")
    const QtWidgets = joinpath(qtlibdir,"QtWidgets.framework/")

    addHeaderDir(qtlibdir; isFramework = true, kind = C_System)

    dlopen(joinpath(QtCore,"QtCore_debug"))
    addHeaderDir(joinpath(QtCore,"Headers"), kind = C_System)
    addHeaderDir("/Users/kfischer/Projects/qt-everywhere-opensource-src-5.3.1/qtbase/lib/QtCore.framework/Headers/5.3.1/QtCore")

    cxxinclude("/Users/kfischer/Projects/qt-everywhere-opensource-src-5.3.1/qtbase/lib/QtCore.framework/Headers/5.3.1/QtCore/private/qcoreapplication_p.h")

    dlopen(joinpath(QtWidgets,"QtWidgets"))
    addHeaderDir(joinpath(QtWidgets,"Headers"), kind = C_System)
end

@linux_only begin
    const qtincdir = "/usr/include/qt5"
    const qtlibdir = "/usr/lib/x86_64-linux-gnu/"

    addHeaderDir(qtincdir, kind = C_System)
    addHeaderDir(QtWidgets, kind = C_System)

    dlopen(joinpath(qtlibdir,"libQt5Core.so"), RTLD_GLOBAL)
    dlopen(joinpath(qtlibdir,"libQt5Gui.so"), RTLD_GLOBAL)
    dlopen(joinpath(qtlibdir,"libQt5Widgets.so"), RTLD_GLOBAL)
end

cxxinclude("QApplication", isAngled=true)
cxxinclude("QMessageBox", isAngled=true)

const a = "julia"
x = Ptr{Uint8}[pointer(a),C_NULL]
# This is pretty stupid, but it seems QApplication is capturing the pointer to the reference, so we can't just
# stack allocate it because that won't be valid for exec
ac = [int32(1)]
app = @cxx QApplication(*(pointer(ac)),pointer(x))

@cxx QMessageBox::about(cast(C_NULL,pcpp"QWidget"), pointer("Hello World"), pointer("This is a QMessageBox!"))

@cxx app->exec()


text:
using SafeTestsets

@safetestset "Day01" begin
    using AdventOfCode.TestUtils
    using AdventOfCode.Day01
    using AdventOfCode.Day01: fuel1, fuel2

    @testset "fuel1" begin
        @test fuel1(12) == 2
        @test fuel1(14) == 2
        @test fuel1(1969) == 654
        @test fuel1(100756) == 33583
    end

    @testset "fuel2" begin
        @test fuel2(14) == 2
        @test fuel2(1969) == 966
        @test fuel2(100756) == 50346
    end

    btest("1.1") do
        part1("input01.dat") == 3348909
    end

    btest("1.2") do
        part2("input01.dat") == 5020494
    end
end


text:
using LinearAlgebra

A = [1 2; 3 4]
b = [5,1]
x = A \ b         # Solve Ax = b for x
A*x == b        # Confirm solution is correct

B = [5 7; 1 -3]
X = A \ B          # Solve for two RHS vectors
A*X == B

n = 2000
T = SymTridiagonal(2ones(n), -ones(n))     # n-by-n symmetric tridiagonal

for rep = 1:3 @time T \ randn(n) end       # Very fast since T is a SymTridiagonal
Tfull = Matrix(T)                          # Convert T to a full 2D array
for rep = 1:3 @time Tfull \ randn(n) end   # Now \ is magnitudes slower

x = 0:0.1:10
n = length(x)
y = 3x .- 2 + randn(n)     # Example data: straight line with noise

A = [ones(n) x]            # LHS
ab = A \ y                 # Least-squares solution

using PyPlot
xplot = 0:10;
yplot = @. ab[1] + ab[2] * xplot
plot(x,y,".")
plot(xplot, yplot, "r");


text:
using Documenter, DataFrames

# Build documentation.
# ====================

makedocs(
    # options
    modules = [DataFrames],
    doctest = false,
    clean = false,
    sitename = "DataFrames.jl",
    format = :html,
    pages = Any[
        "Introduction" => "index.md",
        "User Guide" => Any[
            "Getting Started" => "man/getting_started.md",
            "Joins" => "man/joins.md",
            "Split-apply-combine" => "man/split_apply_combine.md",
            "Reshaping" => "man/reshaping_and_pivoting.md",
            "Sorting" => "man/sorting.md",
            "Categorical Data" => "man/categorical.md",
            "Querying frameworks" => "man/querying_frameworks.md"
        ],
        "API" => Any[
            "Types" => "lib/types.md",
            "Functions" => "lib/functions.md"
        ]
    ]
)

# Deploy built documentation from Travis.
# =======================================

deploydocs(
    # options
    repo = "github.com/JuliaData/DataFrames.jl.git",
    target = "build",
    julia = "0.6",
    deps = nothing,
    make = nothing,
)


text:
## Store
@adddevice(Store, OnePort, LinearExpansionForm{LinearDispatchForm},
           :stores, (:S, :T=>:snapshots), joinpath(@__DIR__, "attrs", "stores.csv"))

cost(d::Store) = sum(d[:marginal_cost] .* d[:p]) + sum(d[:capital_cost] .* (d[:e_nom] - getparam(d, :e_nom)))

function addto!(jm::ModelView, m::EnergyModel, d::Store{DF}) where
      {DDF, DF <: LinearExpansionForm{DDF}}

    addto!(jm, m, with_formulation(d, LinearExpansionInvestmentForm))
    addto!(jm, m, with_formulation(d, LinearExpansionDispatchForm{DDF}))
end

function addto!(jm::ModelView, m::EnergyModel, d::Store{DF}) where
      {DF <: LinearExpansionInvestmentForm}

    S = axis(m, d)

    e_nom_min = get(d, :e_nom_min, S)
    e_nom_max = get(d, :e_nom_max, S)

    @variable(jm, e_nom_min[s] <= e_nom[s=S] <= e_nom_max[s])
end

function addto!(jm::ModelView, m::EnergyModel, d::Store{DF}) where
      {DDF, DF <: LinearExpansionDispatchForm{DDF}}

    S = axis(m, d)
    T = axis(m, :snapshots)

    e_nom = get(d, :e_nom, S)
    e_min_pu = get(d, :e_min_pu, S, T)
    e_max_pu = get(d, :e_max_pu, S, T)

    @variable(jm, e[s=S,t=T])

    @constraints jm begin
        e_lower[s=S,t=T], e[s,t] >= e_min_pu[s,t] * e_nom[s]
        e_upper[s=S,t=T], e[s,t] <= e_max_pu[s,t] * e_nom[s]
    end

    addto!_soc(jm, m, d)
end

function addto!(jm::ModelView, m::EnergyModel, d::Store{DF}) where
      {DF <: DispatchForm}

    S = axis(m, d)
    T = axis(m, :snapshots)

    e_nom = get(d, :e_nom, S)
    e_min_pu = get(d, :e_min_pu, S, T)
    e_max_pu = get(d, :e_max_pu, S, T)

    @variable(jm, e_min_pu[s,t] * e_nom[s] <= e[s=S,t=T] <= e_max_pu[s,t] * e_nom[s])

    addto!_soc(jm, m, d)
end

function addto!_soc(jm::ModelView, m::EnergyModel, d::Store)

    S = axis(m, d)
    T = axis(m, :snapshots)

    e = get(d, :e)
    standing_loss = get(d, :standing_loss, S)

    @variable(jm, p[S,T])

    if d[:e_cyclic]
        e_prev = circshift(e, :snapshots=>1)
    else
        e_prev = similar(e, Union{Float64,eltype(e)})
        e_prev[:,1] .= d[:e_initial]
        e_prev[:,2:end] .= e[:,1:end-1]
    end
    @constraint(jm, e_eq[s=S, t=T],
                e[s,t] - (1 - standing_loss[s]) * e_prev[s,t] == p)
end

max_stars_repo_name: coroa/EnergyModels

max_issues_repo_name: coroa/EnergyModels

max_forks_repo_name: coroa/EnergyModels


text:
# Automatically generated using Clang.jl

const CUSTATEVEC_VER_MAJOR = 0
const CUSTATEVEC_VER_MINOR = 1
const CUSTATEVEC_VER_PATCH = 0
const CUSTATEVEC_VERSION = CUSTATEVEC_VER_MAJOR * 1000 + CUSTATEVEC_VER_MINOR * 100 + CUSTATEVEC_VER_PATCH
const custatevecIndex_t = Int64
const custatevecContext = Cvoid
const custatevecHandle_t = Ptr{custatevecContext}

struct custatevecSamplerDescriptor
    data::NTuple{256, UInt8}
end

const custatevecSamplerDescriptor_t = custatevecSamplerDescriptor

struct custatevecAccessorDescriptor
    data::NTuple{1024, UInt8}
end

const custatevecAccessorDescriptor_t = custatevecAccessorDescriptor
const custatevecLoggerCallback_t = Ptr{Cvoid}

@cenum custatevecStatus_t::UInt32 begin
    CUSTATEVEC_STATUS_SUCCESS = 0
    CUSTATEVEC_STATUS_NOT_INITIALIZED = 1
    CUSTATEVEC_STATUS_ALLOC_FAILED = 2
    CUSTATEVEC_STATUS_INVALID_VALUE = 3
    CUSTATEVEC_STATUS_ARCH_MISMATCH = 4
    CUSTATEVEC_STATUS_EXECUTION_FAILED = 5
    CUSTATEVEC_STATUS_INTERNAL_ERROR = 6
    CUSTATEVEC_STATUS_NOT_SUPPORTED = 7
    CUSTATEVEC_STATUS_INSUFFICIENT_WORKSPACE = 8
    CUSTATEVEC_STATUS_SAMPLER_NOT_PREPROCESSED = 9
end

@cenum custatevecPauli_t::UInt32 begin
    CUSTATEVEC_PAULI_I = 0
    CUSTATEVEC_PAULI_X = 1
    CUSTATEVEC_PAULI_Y = 2
    CUSTATEVEC_PAULI_Z = 3
end

@cenum custatevecMatrixLayout_t::UInt32 begin
    CUSTATEVEC_MATRIX_LAYOUT_COL = 0
    CUSTATEVEC_MATRIX_LAYOUT_ROW = 1
end

@cenum custatevecMatrixType_t::UInt32 begin
    CUSTATEVEC_MATRIX_TYPE_GENERAL = 0
    CUSTATEVEC_MATRIX_TYPE_UNITARY = 1
    CUSTATEVEC_MATRIX_TYPE_HERMITIAN = 2
end

@cenum custatevecCollapseOp_t::UInt32 begin
    CUSTATEVEC_COLLAPSE_NONE = 0
    CUSTATEVEC_COLLAPSE_NORMALIZE_AND_ZERO = 1
end

@cenum custatevecComputeType_t::UInt32 begin
    CUSTATEVEC_COMPUTE_DEFAULT = 0
    CUSTATEVEC_COMPUTE_32F = 4
    CUSTATEVEC_COMPUTE_64F = 16
    CUSTATEVEC_COMPUTE_TF32 = 4096
end

@cenum custatevecSamplerOutput_t::UInt32 begin
    CUSTATEVEC_SAMPLER_OUTPUT_RANDNUM_ORDER = 0
    CUSTATEVEC_SAMPLER_OUTPUT_ASCENDING_ORDER = 1
end


text:
import BufferedStreams, Zlib, GZip, Libz

# Source:
# ftp://ftp.ensembl.org/pub/release-81/fasta/homo_sapiens/cds/Homo_sapiens.GRCh38.cds.all.fa.gz
const filename = "Homo_sapiens.GRCh38.cds.all.fa"

# this just tests overhead, but IOStream has an advantage since knows it's a
# nop.
#const output_filename = "/dev/null"
const output_filename = "out.tmp"

const data = readall(open(filename))

const methods = [
    ("IOStream",                      () -> open(output_filename, "w")),
    ("BufferedOutputStream/IOStream", () -> BufferedStreams.BufferedOutputStream(open(output_filename, "w"))),
    ("BufferedOutputStream",          () -> BufferedStreams.BufferedOutputStream()),
    ("IOBuffer",                      () -> IOBuffer()),
    ("GZip",                          () -> GZip.gzopen(output_filename, "w")),
    ("Zlib",                          () -> Zlib.Writer(open(output_filename, "w"))),
    ("Libz",                          () -> Libz.ZlibDeflateOutputStream(open(output_filename, "w"))),
]

function bench(f)
    for (method_name, open_func) in methods
        print("  ", method_name, ":  ")
        f(open_func())
        println(round(f(open_func()), 4))
    end
end

function bench_write(output)
    tic()
    for c in data
        write(output, c)
    end
    close(output)
    return toq()
end

function bench_write_array(output)
    tic()
    write(output, data)
    close(output)
    return toq()
end

println("write")
bench(bench_write)

println("write array")
bench(bench_write_array)


text:
function test_crop_layer(backend::Backend, do_mirror, n_input, T, eps)
  println("-- Testing CropLayer on $(typeof(backend)){$T} $(do_mirror ? "with mirror" : "")...")

  dims = [rand(7:13, 4) for i = 1:n_input]
  input = [rand(T, dims[i]...) for i = 1:n_input]
  crop_size = (7,5)
  input_blob = Blob[make_blob(backend, x) for x in input]
  diff_blob = Blob[NullBlob() for i = 1:n_input]

  println("    > Setup")
  layer = CropLayer(bottoms=Array{Symbol}(n_input), tops=Array{Symbol}(n_input),
      crop_size=crop_size, random_mirror=do_mirror)
  state = setup(backend, layer, input_blob, diff_blob)

  println("    > Forward")
  forward(backend, state, input_blob)
  for i = 1:n_input
    got_output = zeros(T, size(state.blobs[i]))
    copy!(got_output, state.blobs[i])

    w_off = div(size(input[i],1)-crop_size[1],2)
    h_off = div(size(input[i],2)-crop_size[2],2)
    expected_output = input[i][w_off+1:w_off+crop_size[1], h_off+1:h_off+crop_size[2],:,:]
    @test size(expected_output) == size(got_output)
    if do_mirror
      @test all(abs.(got_output - expected_output) .< eps) ||
            all(abs.(got_output - flipdim(expected_output,1)) .< eps)
    else
      @test all(abs.(got_output - expected_output) .< eps)
    end
  end

  shutdown(backend, state)
end
function test_crop_layer_random(backend::Backend, do_mirror, n_input, T, eps)
  println("-- Testing CropLayer{rnd} on $(typeof(backend)){$T} $(do_mirror ? "with mirror" : "")...")
  dims = [rand(9:11, 4) for i = 1:n_input]
  inputs = [rand(T, dims[i]...) for i = 1:n_input]
  crop_size = (9, 9)
  input_blob = Blob[make_blob(backend, x) for x in inputs]
  diff_blob = Blob[NullBlob() for i = 1:n_input]

  println("    > Setup")
  layer = CropLayer(bottoms=Array{Symbol}(n_input), tops=Array{Symbol}(n_input), crop_size=crop_size,
      random_mirror=do_mirror, random_crop=true)
  state = setup(backend, layer, input_blob, diff_blob)

  println("    > Forward")
  forward(backend, state, input_blob)
  for kk = 1:n_input
    got_output = zeros(T, size(state.blobs[kk]))
    copy!(got_output, state.blobs[kk])
    input = inputs[kk]

    matched = false
    for i = 0:size(input,1)-crop_size[1]
      for j = 0:size(input,2)-crop_size[2]
        expected_output = input[i+1:i+crop_size[1], j+1:j+crop_size[2],:,:]
        matched = matched | all(abs.(got_output - expected_output) .< eps)
        if do_mirror
          matched = matched | all(abs.(got_output - flipdim(expected_output,1)) .< eps)
        end
      end
    end
    @test matched
  end

  shutdown(backend, state)
end

function test_crop_layer(backend::Backend, n_input, T, eps)
  test_crop_layer(backend, false, n_input, T, eps)
  test_crop_layer(backend, true, n_input, T, eps)
  test_crop_layer_random(backend, false, n_input, T, eps)
  test_crop_layer_random(backend, true, n_input, T, eps)
end

function test_crop_layer(backend::Backend)
  test_crop_layer(backend, 3, Float64, 1e-10)
  test_crop_layer(backend, 3, Float32, 1e-4)
end

if test_cpu
  test_crop_layer(backend_cpu)
end
if test_gpu
  test_crop_layer(backend_gpu)
end


text:
export UDWT_Filter_Haar, UDWT_Filter_Starck2
export ϕ_filter,ψ_filter,tildeϕ_filter,tildeψ_filter
export udwt, scale, inverse_udwt!, inverse_udwt

import Base: length

"""
    abstract type UDWT_Filter_Biorthogonal{T<:Number} end

Abstract type defining the ϕ, ψ, tildeϕ and tildeψ filters associated
to an undecimated biorthogonal wavelet transform

Subtypes of this structure are:
  * [`UDWT_Filter`](@ref UDWT_Filter)

Associated methods are:
  * [`ϕ_filter`](@ref ϕ_filter)
  * `ψ_filter`
  * `tildeϕ_filter`
  * `tildeψ_filter`
"""
abstract type UDWT_Filter_Biorthogonal{T<:Number} end

"""
    ϕ_filter(c::UDWT_Filter_Biorthogonal)

Returns the ϕ filter
"""
ϕ_filter(c::UDWT_Filter_Biorthogonal)::LinearFilter = c._ϕ_filter
ψ_filter(c::UDWT_Filter_Biorthogonal)::LinearFilter = c._ψ_filter
tildeϕ_filter(c::UDWT_Filter_Biorthogonal)::LinearFilter = c._tildeϕ_filter
tildeψ_filter(c::UDWT_Filter_Biorthogonal)::LinearFilter = c._tildeψ_filter

"""
    abstract type UDWT_Filter{T<:Number} <: UDWT_Filter_Biorthogonal{T} end

A specialization of `UDWT_Filter_Biorthogonal` for **orthogonal** filters.

For orthogonal filters we have:  ϕ=tildeϕ, ψ=tildeψ
"""
abstract type UDWT_Filter{T<:Number} <: UDWT_Filter_Biorthogonal{T}
end

#+UDWT_Filter
tildeϕ_filter(c::UDWT_Filter)::LinearFilter = ϕ_filter(c)

#+UDWT_Filter
tildeψ_filter(c::UDWT_Filter)::LinearFilter = ψ_filter(c)

# Filter examples

#+UDWT_Filter_Available
# Haar filter
struct UDWT_Filter_Haar{T<:AbstractFloat} <: UDWT_Filter{T}
    _ϕ_filter::LinearFilter_Default{T,2}
    _ψ_filter::LinearFilter_Default{T,2}

    #+UDWT_Filter_Available
    # Creates an instance
    UDWT_Filter_Haar{T}() where {T<:Real} = new(LinearFilter_Default{T,2}(SVector{2,T}([+1/2 +1/2]),0),
                                                LinearFilter_Default{T,2}(SVector{2,T}([-1/2 +1/2]),0))
end

#+UDWT_Filter_Available
# Starck2 filter
#
# Defined by Eq. 6 from http://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=4060954
struct UDWT_Filter_Starck2{T<:AbstractFloat} <: UDWT_Filter_Biorthogonal{T}
    _ϕ_filter::LinearFilter_DefaultCentered{T,5}
    _ψ_filter::LinearFilter_DefaultCentered{T,5}
    _tildeϕ_filter::LinearFilter_DefaultCentered{T,1}
    _tildeψ_filter::LinearFilter_DefaultCentered{T,1}

    #+UDWT_Filter_Available
    # Creates an instance
    UDWT_Filter_Starck2{T}() where {T<:Real} = new(LinearFilter_DefaultCentered{T,5}(SVector{5,T}([+1/16 +4/16 +6/16 +4/16 +1/16])),
                                                   LinearFilter_DefaultCentered{T,5}(SVector{5,T}([-1/16 -4/16 +10/16 -4/16 -1/16])),
                                                   LinearFilter_DefaultCentered{T,1}(SVector{1,T}([+1])),
                                                   LinearFilter_DefaultCentered{T,1}(SVector{1,T}([+1])))
end

#+UDWT
# A structure to store 1D UDWT
#
struct UDWT{T<:Number}

    filter::UDWT_Filter_Biorthogonal{T}
    # TODO also store boundary condition

    W::Array{T,2}
    V::Array{T,1}

    #+UDWT
    # Creates an instance
    #
    # *Parameters:*
    # - *filter*: used filter
    # - *scale* : max scale
    # - *n*: signal length
    #
    UDWT{T}(filter::UDWT_Filter_Biorthogonal{T};
            n::Int=0,
            scale::Int=0) where {T<:Number} =
                new(filter,
                    Array{T,2}(undef,n,scale),
                    Array{T,1}(undef,n))
end

#+UDWT
# Returns max scale
scale(udwt::UDWT)::Int = size(udwt.W,2)
#+UDWT
# Returns expected signal length
length(udwt::UDWT)::Int = size(udwt.W,1)

#+UDWT
#
# Performs an 1D undecimated wavelet transform
#
# $$(\mathcal{W}_{j+1}f)[u]=(\bar{g}_j*\mathcal{V}_{j}f)[u]$$
# $$(\mathcal{V}_{j+1}f)[u]=(\bar{h}_j*\mathcal{V}_{j}f)[u]$$
#
function udwt(signal::AbstractArray{T,1},
              filter::UDWT_Filter_Biorthogonal{T};
              scale::Int=3) where {T<:Number}

    @assert scale>=0

    boundary = PeriodicBE
    n = length(signal)
    udwt_domain = UDWT{T}(filter,n=n,scale=scale)
    Ωγ = 1:n

    Vs = Array{T,1}(undef,n)
    Vsp1 = Array{T,1}(undef,n)
    Vs .= signal

    for s in 1:scale
        twoPowScale = 2^(s-1)
        Wsp1 = @view udwt_domain.W[:,s]

        # Computes Vs+1 from Vs
        #
        directConv!(ϕ_filter(filter),
                    twoPowScale,

                    Vs,

                    Vsp1,
                    Ωγ,

                    boundary,
                    boundary)

        # Computes Ws+1 from Ws
        #
        directConv!(ψ_filter(filter),
                    twoPowScale,

                    Vs,

                    Wsp1,
                    Ωγ,

                    boundary,
                    boundary)

        @swap(Vs,Vsp1)

    end

    udwt_domain.V .= Vs

    return udwt_domain
end

# +UDWT
#
# Performs an 1D *inverse* undecimated wavelet transform
#
# *Caveat:* uses a pre-allocated vector =reconstructed_signal=
#
function inverse_udwt!(udwt_domain::UDWT{T},
                       reconstructed_signal::AbstractArray{T,1}) where {T<:Number}

    @assert length(udwt_domain) == length(reconstructed_signal)

    boundary = PeriodicBE
    maxScale = scale(udwt_domain)
    n = length(reconstructed_signal)
    Ωγ = 1:n
    buffer = Array{T,1}(undef,n)

    reconstructed_signal .= udwt_domain.V

    for s in maxScale:-1:1
        twoPowScale = 2^(s-1)

        # Computes Vs from Vs+1
        #
        directConv!(tildeϕ_filter(udwt_domain.filter),
                    -twoPowScale,

                    reconstructed_signal,

                    buffer,
                    Ωγ,

                    boundary,
                    boundary)

        # Computes Ws from Ws+1
        #
        Ws = @view udwt_domain.W[:,s]

        directConv!(tildeψ_filter(udwt_domain.filter),
                    -twoPowScale,

                    Ws,

                    buffer,
                    Ωγ,

                    boundary,
                    boundary,
                    accumulate=true)

        for i in Ωγ
            @inbounds reconstructed_signal[i]=buffer[i]
        end
    end
end

#+UDWT
#
# Performs an 1D *inverse* undecimated wavelet transform
#
# *Returns:* a vector containing the reconstructed signal.
#
function inverse_udwt(udwt_domain::UDWT{T})::Array{T,1} where {T<:Number}
    reconstructed_signal=Array{T,1}(undef,length(udwt_domain))
    inverse_udwt!(udwt_domain,reconstructed_signal)
    return reconstructed_signal
end

max_stars_repo_name: vincent-picaud/DirectConvolution

max_issues_repo_name: vincent-picaud/DirectConvolution

max_forks_repo_name: vincent-picaud/DirectConvolution


text:
module TestIntervals

using Base.Test

using Bio.Intervals,
    Distributions,
    YAML,
    TestFunctions

# Test that an array of intervals is well ordered
function Intervals.isordered{I <: Interval}(intervals::Vector{I})
    for i = 2:length(intervals)
        if !Intervals.isordered(intervals[i-1], intervals[i])
            return false
        end
    end
    return true
end

# Generate an array of n random Interval{Int} object. With sequence names
# samples from seqnames, and intervals drawn to lie in [1, maxpos].
function random_intervals(seqnames, maxpos::Int, n::Int)
    seq_dist = Categorical(length(seqnames))
    strand_dist = Categorical(2)
    length_dist = Normal(1000, 1000)

    intervals = Array(Interval{Int}, n)
    for i in 1:n
        intlen = maxpos
        while intlen >= maxpos || intlen <= 0
            intlen = ceil(Int, rand(length_dist))
        end
        first = rand(1:maxpos-intlen)
        last = first + intlen - 1
        strand = rand(strand_dist) == 1 ? STRAND_POS : STRAND_NEG
        intervals[i] = Interval{Int}(seqnames[rand(seq_dist)],
                                     first, last, strand, i)
    end
    return intervals
end

# A simple interval intersection implementation to test against.
function simple_intersection(intervals_a, intervals_b)
    sort!(intervals_a)
    sort!(intervals_b)

    intersections = Any[]

    i = 1
    j = 1
    while i <= length(intervals_a) && j <= length(intervals_b)
        ai = intervals_a[i]
        bj = intervals_b[j]

        if Intervals.alphanum_isless(ai.seqname, bj.seqname) ||
           (ai.seqname == bj.seqname && ai.last < bj.first)
            i += 1
        elseif Intervals.alphanum_isless(bj.seqname, ai.seqname) ||
               (ai.seqname == bj.seqname && bj.last < ai.first)
            j += 1
        else
            k = j
            while k <= length(intervals_b) && intervals_b[k].first <= ai.last
                if isoverlapping(ai, intervals_b[k])
                    push!(intersections, (ai, intervals_b[k]))
                end
                k += 1
            end
            i += 1
        end
    end

    return intersections
end

function simple_coverage(intervals)
    seqlens = Dict{AbstractString, Int}()
    for interval in intervals
        if get(seqlens, interval.seqname, -1) < interval.last
            seqlens[interval.seqname] = interval.last
        end
    end

    covarrays = Dict{AbstractString, Vector{Int}}()
    for (seqname, seqlen) in seqlens
        covarrays[seqname] = zeros(Int, seqlen)
    end

    for interval in intervals
        arr = covarrays[interval.seqname]
        for i in interval.first:interval.last
            arr[i] += 1
        end
    end

    covintervals = Interval{UInt32}[]
    for (seqname, arr) in covarrays
        i = j = 1
        while i <= length(arr)
            if arr[i] > 0
                j = i + 1
                while j <= length(arr) && arr[j] == arr[i]
                    j += 1
                end
                push!(covintervals,
                      Interval{UInt32}(seqname, i, j - 1, STRAND_BOTH, arr[i]))
                i = j
            else
                i += 1
            end
        end
    end

    return covintervals
end

@testset "Strand" begin
    @testset "Constructor" begin
        @test Strand('?') === STRAND_NA
        @test Strand('+') === STRAND_POS
        @test Strand('-') === STRAND_NEG
        @test Strand('.') === STRAND_BOTH
        @test_throws Exception Strand('x')
    end

    @testset "Conversion" begin
        @test convert(Strand, '?') === STRAND_NA
        @test convert(Strand, '+') === STRAND_POS
        @test convert(Strand, '-') === STRAND_NEG
        @test convert(Strand, '.') === STRAND_BOTH

        @test convert(Char, STRAND_NA) === '?'
        @test convert(Char, STRAND_POS) === '+'
        @test convert(Char, STRAND_NEG) === '-'
        @test convert(Char, STRAND_BOTH) === '.'
    end

    @testset "Order" begin
        @test STRAND_NA < STRAND_POS < STRAND_NEG < STRAND_BOTH
    end

    @testset "Show" begin
        @testset "show" begin
            buf = IOBuffer()
            for s in [STRAND_NA, STRAND_POS, STRAND_NEG, STRAND_BOTH]
                show(buf, s); print(buf, " ")
            end
            @test takebuf_string(buf) == "STRAND_NA STRAND_POS STRAND_NEG STRAND_BOTH "
        end

        @testset "print" begin
            buf = IOBuffer()
            for s in [STRAND_NA, STRAND_POS, STRAND_NEG, STRAND_BOTH]
                print(buf, s)
            end
            @test takebuf_string(buf) == "?+-."
        end
    end
end

@testset "Interval" begin
    @testset "Constructor" begin
        i = Interval("chr1", 10, 20)
        @test seqname(i) == "chr1"
        @test leftposition(i) == 10
        @test rightposition(i) == 20
        @test strand(i) == STRAND_BOTH

        i1 = Interval("chr1", 10, 20, '+')
        i2 = Interval("chr1", 10, 20, STRAND_POS)
        @test i1 == i2

        i1 = Interval("chr2", 5692667, 5701385, '+',        "SOX11")
        i2 = Interval("chr2", 5692667, 5701385, STRAND_POS, "SOX11")
        @test i1 == i2
    end
end

@testset "IntervalCollection" begin

    @testset "Insertion/Iteration" begin
        n = 100000
        intervals = random_intervals(["one", "two", "three"], 1000000, n)
        ic = IntervalCollection{Int}()

        @test isempty(ic)
        @test collect(Interval{Int}, ic) == Interval{Int}[]

        for interval in intervals
            push!(ic, interval)
        end
        @test Intervals.isordered(collect(Interval{Int}, ic))
    end

    @testset "Intersection" begin
        n = 1000
        srand(1234)
        intervals_a = random_intervals(["one", "two", "three"], 1000000, n)
        intervals_b = random_intervals(["one", "three", "four"], 1000000, n)

        # empty versus empty
        ic_a = IntervalCollection{Int}()
        ic_b = IntervalCollection{Int}()
        @test collect(intersect(ic_a, ic_b)) == Any[]

        # empty versus non-empty
        for interval in intervals_a
            push!(ic_a, interval)
        end

        @test collect(intersect(ic_a, ic_b)) == Any[]
        @test collect(intersect(ic_b, ic_a)) == Any[]

        # non-empty versus non-empty
        for interval in intervals_b
            push!(ic_b, interval)
        end

        @test sort(collect(intersect(ic_a, ic_b))) == sort(simple_intersection(intervals_a, intervals_b))
    end

    @testset "Show" begin
        ic = IntervalCollection{Int}()
        show(DevNull, ic)

        push!(ic, Interval{Int}("one", 1, 1000, STRAND_POS, 0))
        show(DevNull, ic)

        intervals = random_intervals(["one", "two", "three"], 1000000, 100)
        for interval in intervals
            push!(ic, interval)
        end
        show(DevNull, ic)

        show(DevNull, STRAND_NA)
        show(DevNull, STRAND_POS)
        show(DevNull, STRAND_NEG)
        show(DevNull, STRAND_BOTH)
    end
end

@testset "Alphanumeric Sorting" begin
    @test sort(["b", "c" ,"a"], lt=Intervals.alphanum_isless) == ["a", "b", "c"]
    @test sort(["a10", "a2" ,"a1"], lt=Intervals.alphanum_isless) == ["a1", "a2", "a10"]
    @test sort(["a10a", "a2c" ,"a3b"], lt=Intervals.alphanum_isless) == ["a2c", "a3b", "a10a"]
    @test sort(["a3c", "a3b" ,"a3a"], lt=Intervals.alphanum_isless) == ["a3a", "a3b", "a3c"]
    @test sort(["a1ac", "a1aa" ,"a1ab"], lt=Intervals.alphanum_isless) == ["a1aa", "a1ab", "a1ac"]

    @test Intervals.alphanum_isless("aa", "aa1")
    @test !Intervals.alphanum_isless("aa1", "aa")

    @test sort(["ac3", "aa", "ab", "aa1", "ac", "ab2"], lt=Intervals.alphanum_isless) ==
        ["aa", "aa1", "ab", "ab2", "ac", "ac3"]
end

@testset "IntervalStream" begin
    @testset "StreamBuffer" begin
        ref = Int[]
        sb = Intervals.StreamBuffer{Int}()
        @test isempty(sb)
        @test length(sb) == 0
        @test_throws Exception shift!(sb)

        ref_shifts = Int[]
        sb_shifts = Int[]

        for i in 1:10000
            if !isempty(ref) && rand() < 0.3
                push!(ref_shifts, shift!(ref))
                push!(sb_shifts, shift!(sb))
            else
                x = rand(Int)
                push!(ref, x)
                push!(sb, x)
            end
        end

        @test length(sb) == length(ref)
        @test [sb[i] for i in 1:length(sb)] == ref
        @test ref_shifts == sb_shifts
        @test_throws BoundsError sb[0]
        @test_throws BoundsError sb[length(sb) + 1]
    end

    @testset "Intersection" begin
        n = 1000
        srand(1234)
        intervals_a = random_intervals(["one", "two", "three"], 1000000, n)
        intervals_b = random_intervals(["one", "three", "four"], 1000000, n)

        ic_a = IntervalCollection{Int}()
        ic_b = IntervalCollection{Int}()

        for interval in intervals_a
            push!(ic_a, interval)
        end

        for interval in intervals_b
            push!(ic_b, interval)
        end

        # non-empty versus non-empty, stream intersection
        it = Intervals.IntervalStreamIntersectIterator{Int, Int,
                IntervalCollection{Int}, IntervalCollection{Int}}(
                ic_a, ic_b, Intervals.alphanum_isless)

        @test sort(collect(it)) == sort(simple_intersection(intervals_a, intervals_b))

        # Interesction edge cases: skipping over whole sequences
        typealias SimpleIntersectIterator
            Intervals.IntervalStreamIntersectIterator{Void, Void,
                Vector{Interval{Void}}, Vector{Interval{Void}}}

        it = SimpleIntersectIterator(
            [Interval("a", 1, 100, STRAND_POS, nothing), Interval("c", 1, 100, STRAND_POS, nothing)],
            [Interval("a", 1, 100, STRAND_POS, nothing), Interval("b", 1, 100, STRAND_POS, nothing)],
            isless)
        @test length(collect(it)) == 1

        it = SimpleIntersectIterator(
            [Interval("c", 1, 100, STRAND_POS, nothing), Interval("d", 1, 100, STRAND_POS, nothing)],
            [Interval("b", 1, 100, STRAND_POS, nothing), Interval("d", 1, 100, STRAND_POS, nothing)],
            isless)
        @test length(collect(it)) == 1

        # unsorted streams are not allowed
        @test_throws Exception begin
            it = SimpleIntersectIterator(
                [Interval("b", 1, 1000, STRAND_POS, nothing),
                 Interval("a", 1, 1000, STRAND_POS, nothing)],
                [Interval("a", 1, 1000, STRAND_POS, nothing),
                 Interval("b", 1, 1000, STRAND_POS, nothing)], isless)
            collect(it)
        end

        @test_throws Exception begin
            it = SimpleIntersectIterator(
                [Interval("a", 1, 1000, STRAND_POS, nothing),
                 Interval("a", 500, 1000, STRAND_POS, nothing),
                 Interval("a", 400, 2000, STRAND_POS, nothing)],
                [Interval("a", 1, 1000, STRAND_POS, nothing),
                 Interval("b", 1, 1000, STRAND_POS, nothing)], isless)
            collect(it)
        end
    end

    @testset "IntervalStream Intersection" begin
        n = 1000
        srand(1234)
        intervals_a = random_intervals(["one", "two", "three"], 1000000, n)
        intervals_b = random_intervals(["one", "two", "three"], 1000000, n)

        ic_a = IntervalCollection{Int}()
        ic_b = IntervalCollection{Int}()

        for interval in intervals_a
            push!(ic_a, interval)
        end

        for interval in intervals_b
            push!(ic_b, interval)
        end

        ItType = Intervals.IntervalStreamIntersectIterator{Int, Int,
            IntervalCollection{Int}, IntervalCollection{Int}}

        @test sort(collect(ItType(ic_a, ic_b, isless))) == sort(simple_intersection(intervals_a, intervals_b))
    end

    @testset "IntervalStream Coverage" begin
        n = 10000
        srand(1234)
        intervals = random_intervals(["one", "two", "three"], 1000000, n)

        ic = IntervalCollection{Int}()
        for interval in intervals
            push!(ic, interval)
        end

        @test sort(simple_coverage(intervals)) == sort(collect(coverage(ic)))
    end
end

@testset "Interval Parsing" begin
    @testset "BED Parsing" begin
        get_bio_fmt_specimens()

        function check_bed_parse(filename)
            # Reading from a stream
            for interval in BEDReader(open(filename))
            end

            # Reading from a regular file
            for interval in open(BEDReader, filename)
            end

            # in-place parsing
            stream = open(BEDReader, filename)
            entry = eltype(stream)()
            while !eof(stream)
                read!(stream, entry)
            end
            close(stream)

            # Check round trip
            output = IOBuffer()
            writer = BEDWriter(output)
            expected_entries = Any[]
            for interval in open(BEDReader, filename)
                write(writer, interval)
                push!(expected_entries, interval)
            end
            flush(writer)

            seekstart(output)
            read_entries = Any[]
            for interval in BEDReader(output)
                push!(read_entries, interval)
            end

            return expected_entries == read_entries
        end

        path = Pkg.dir("Bio", "test", "BioFmtSpecimens", "BED")
        for specimen in YAML.load_file(joinpath(path, "index.yml"))
            valid = get(specimen, "valid", true)
            if valid
                @test check_bed_parse(joinpath(path, specimen["filename"]))
            else
                @test_throws Exception check_bed_parse(joinpath(path, specimen["filename"]))
            end
        end
    end

    @testset "BED Intersection" begin
        # Testing strategy: there are two entirely separate intersection
        # algorithms for IntervalCollection and IntervalStream. Here we test
        # them both by checking that they agree by generating and intersecting
        # random BED files.

        function check_intersection(filename_a, filename_b)
            ic_a = IntervalCollection{BEDMetadata}()
            open(BEDReader, filename_a) do intervals
                for interval in intervals
                    push!(ic_a, interval)
                end
            end

            ic_b = IntervalCollection{BEDMetadata}()
            open(BEDReader, filename_b) do intervals
                for interval in intervals
                    push!(ic_b, interval)
                end
            end

            # This is refactored out to close streams
            fa = open(BEDReader, filename_a)
            fb = open(BEDReader, filename_b)
            xs = sort(collect(intersect(fa, fb)))
            close(fa)
            close(fb)

            ys = sort(collect(intersect(ic_a, ic_b)))

            return xs == ys
        end

        function write_intervals(filename, intervals)
            open(filename, "w") do out
                for interval in sort(intervals)
                    println(out, interval.seqname, "\t", interval.first - 1,
                            "\t", interval.last, "\t", interval.metadata, "\t",
                            1000, "\t", interval.strand)
                end
            end

        end

        n = 10000
        srand(1234)
        intervals_a = random_intervals(["one", "two", "three", "four", "five"],
                                       1000000, n)
        intervals_b = random_intervals(["one", "two", "three", "four", "five"],
                                       1000000, n)

        filename_a = "test_a.bed"
        filename_b = "test_b.bed"
        intempdir() do
            write_intervals(filename_a, intervals_a)
            write_intervals(filename_b, intervals_b)
            @test check_intersection(filename_a, filename_b)
        end

    end
end

@testset "BigBed" begin
    @testset "BED → BigBed → BED round-trip" begin
        path = Pkg.dir("Bio", "test", "BioFmtSpecimens", "BED")
        for specimen in YAML.load_file(joinpath(path, "index.yml"))
            if !get(specimen, "valid", true)
                continue
            end

            # BED → BigBed
            intervals = IntervalCollection(
                open(BEDReader, joinpath(path, specimen["filename"])))
            out = IOBuffer()
            write(BigBedWriter(out), intervals)

            # BigBed → BED
            seekstart(out)
            bb = BigBedReader(out)
            intervals2 = IntervalCollection(bb)

            @test intervals == intervals2
        end
    end

    @testset "BigBed Intersection" begin
        n = 10000
        srand(1234)
        chroms = ["one", "two", "three", "four", "five"]
        intervals = IntervalCollection(
            [Interval{BEDMetadata}(i.seqname, i.first, i.last, STRAND_NA, BEDMetadata())
             for i in random_intervals(chroms, 1000000, n)], true)

        # convert to bigbed in memory
        out = IOBuffer()
        write(BigBedWriter(out), intervals)
        seekstart(out)
        bb = BigBedReader(out)

        # intersection queries
        num_queries = 1000
        queries = random_intervals(chroms, 1000000, num_queries)

        @test all(Bool[IntervalCollection(collect(BEDInterval, intersect(intervals, query))) ==
                       IntervalCollection(collect(BEDInterval, intersect(bb, query)))
                       for query in queries])
    end

    # TODO: test summary information against output from kent's bigBedSummary
end

end # module TestIntervals


text:
dimension(::Type{gt.AbstractSimplex{S, gt.Vec{N, T}}}) where {N,T,S} = Val(N)

dimension(::Type{M}) where {M <: gt.AbstractMesh} = dimension(gt.vertextype(M))

dimension(::Type{gt.Vec{N, T}}) where {N,T} = Val(N)

dimension(::Type{gt.Point{N, T}}) where {N,T} = Val(N)

dimension(::Type{gt.Simplex{M, T}}) where {M,T} = dimension(T)

dimension(::Type{SVector{N, T}}) where {N,T} = Val(N)

dimension(::Type{gt.FlexibleConvexHull{T}}) where {T} = dimension(T)

dimension(::Type{gt.HyperRectangle{N, T}}) where {N,T} = Val(N)


text:
using Test
using LinearAlgebra
using FlightMechanics
using FlightMechanics.Aircrafts
using FlightMechanics.Models

ac = F16()

h = 0.0 * M2FT
psi = 0.0  # rad
gamma = 0.0
turn_rate = 0.0

pos = EarthPosition(0, 0, -h)
env = Environment(pos, atmos = "ISAF16", wind = "NoWind", grav = "const")

#  Stevens, B. L., Lewis, F. L., & Johnson, E. N. (2015). Aircraft control
#  and simulation: dynamics, controls design, and autonomous systems. John Wiley
#  & Sons. (page 193 table 3.6-2)
tas =  500 * FT2M
exp_thtl = 0.137
exp_α    = 2.14
exp_de   = -0.756

# Initial conditions for trimmer
α0 = exp_α * DEG2RAD + 5 * DEG2RAD
β0 = 0.
stick_lon0 = exp_de / 25.0 + 0.5
thtl0 = exp_thtl + 0.3

controls = StickPedalsLeverControls(stick_lon0, 0.5, 0.5, thtl0)

ac_trim, aerostate_trim, state_trim, controls_trim = steady_state_trim(
    ac, controls, env, tas, pos, psi, gamma, turn_rate, α0, 0.0,
    )

fcs_trim = get_fcs(ac)

# TEST: Check results against Stevens (also checked in f16.jl tests)
@test isapprox(ac_trim.pfm.forces, zeros(3), atol = 1e-7)
@test isapprox(ac_trim.pfm.moments, zeros(3), atol = 1e-7)
@test isapprox(get_tas(aerostate_trim), tas)
@test isapprox(get_beta(aerostate_trim), 0.0, atol = 1e-13)

@test isapprox(fcs_trim.thtl.value, exp_thtl, atol = 0.001)
@test isapprox(aerostate_trim.alpha * RAD2DEG, exp_α, atol = 0.01)
@test isapprox(fcs_trim.de.value * RAD2DEG, exp_de, atol = 0.001)

# TEST: calculate aircraft on trimmed aricraft and check that forces and moments
# are the same.
ac_calc = calculate_aircraft(
    ac_trim,
    controls_trim,
    aerostate_trim,
    state_trim,
    get_gravity(env),
    )
@test isapprox(ac_trim.pfm.forces, ac_calc.pfm.forces)
@test isapprox(ac_trim.pfm.moments, ac_calc.pfm.moments)

# TEST: Trim a trimmed aircraft and check if returns the same trim
ac_trim2, aerostate_trim2, state_trim2, controls_trim2 = steady_state_trim(
    ac_trim,
    controls_trim,
    env,
    tas,
    pos,
    psi,
    gamma,
    turn_rate,
    get_alpha(aerostate_trim),
    get_beta(aerostate_trim),
    )

fcs_trim2 = get_fcs(ac)

@test isapprox(get_tas(aerostate_trim2), tas)
@test isapprox(get_beta(aerostate_trim2),  0.0, atol = 1e-13)

@test isapprox(fcs_trim2.thtl.value, exp_thtl, atol = 0.001)
@test isapprox(aerostate_trim2.alpha * RAD2DEG, exp_α, atol = 0.01)
@test isapprox(fcs_trim2.de.value * RAD2DEG, exp_de, atol = 0.001)

@test isapprox(ac_trim2.pfm.forces, ac_trim.pfm.forces)
@test isapprox(ac_trim2.pfm.moments, ac_trim.pfm.moments)


text:
@doc raw"""
    trust_regions(M, F, ∇F, x, H)

evaluate the Riemannian trust-regions solver for optimization on manifolds.
It will attempt to minimize the cost function F on the Manifold M.
If no Hessian H is provided, a standard approximation of the Hessian based on
the gradient ∇F will be computed.
For solving the the inner trust-region subproblem of finding an update-vector,
it uses the Steihaug-Toint truncated conjugate-gradient method.
For a description of the algorithm and more details see

* P.-A. Absil, C.G. Baker, K.A. Gallivan,
    Trust-region methods on Riemannian manifolds, FoCM, 2007.
    doi: [10.1007/s10208-005-0179-9](https://doi.org/10.1007/s10208-005-0179-9)
* A. R. Conn, N. I. M. Gould, P. L. Toint, Trust-region methods, SIAM,
    MPS, 2000. doi: [10.1137/1.9780898719857](https://doi.org/10.1137/1.9780898719857)

# Input
* `M` – a manifold $\mathcal M$
* `F` – a cost function $F \colon \mathcal M \to \mathbb R$ to minimize
* `∇F`- the gradient $\nabla F \colon \mathcal M \to T \mathcal M$ of $F$
* `x` – an initial value $x  ∈  \mathcal M$
* `H` – the hessian $H( \mathcal M, x, \xi)$ of $F$

# Optional
* `retraction` – approximation of the exponential map
* `preconditioner` – a preconditioner (a symmetric, positive definite operator
  that should approximate the inverse of the Hessian)
* `stopping_criterion` – ([`StopWhenAny`](@ref)([`StopAfterIteration`](@ref)`(1000)`,
  [`StopWhenGradientNormLess`](@ref)`(10^(-6))`) a functor inheriting
  from [`StoppingCriterion`](@ref) indicating when to stop.
* `Δ_bar` – the maximum trust-region radius
* `Δ` - the (initial) trust-region radius
* `useRandom` – set to true if the trust-region solve is to be initiated with a
  random tangent vector. If set to true, no preconditioner will be
  used. This option is set to true in some scenarios to escape saddle
  points, but is otherwise seldom activated.
* `ρ_prime` – Accept/reject threshold: if ρ (the performance ratio for the
  iterate) is at least ρ', the outer iteration is accepted.
  Otherwise, it is rejected. In case it is rejected, the trust-region
  radius will have been decreased. To ensure this, ρ' >= 0 must be
  strictly smaller than 1/4. If ρ_prime is negative, the algorithm is not
  guaranteed to produce monotonically decreasing cost values. It is
  strongly recommended to set ρ' > 0, to aid convergence.
* `ρ_regularization` – Close to convergence, evaluating the performance ratio ρ
  is numerically challenging. Meanwhile, close to convergence, the
  quadratic model should be a good fit and the steps should be
  accepted. Regularization lets ρ go to 1 as the model decrease and
  the actual decrease go to zero. Set this option to zero to disable
  regularization (not recommended). When this is not zero, it may happen
  that the iterates produced are not monotonically improving the cost
  when very close to convergence. This is because the corrected cost
  improvement could change sign if it is negative but very small.
* `return_options` – (`false`) – if actiavated, the extended result, i.e. the
  complete [`Options`](@ref) are returned. This can be used to access recorded values.
  If set to false (default) just the optimal value `xOpt` is returned

# Output
* `x` – the last reached point on the manifold

# see also
[`truncated_conjugate_gradient_descent`](@ref)
"""
function trust_regions(
    M::MT,
    F::Function,
    ∇F::Function,
    x,
    H::Union{Function,Missing};
    retraction::Function = exp,
    preconditioner::Function = (M,x,ξ) -> ξ,
    stopping_criterion::StoppingCriterion = StopWhenAny(
        StopAfterIteration(1000),
        StopWhenGradientNormLess(10^(-6))
    ),
    Δ_bar = sqrt(manifold_dimension(M)),
    Δ = Δ_bar/8,
    useRandom::Bool = false,
    ρ_prime::Float64 = 0.1,
    ρ_regularization=1000.,
    return_options=false,
    kwargs... #collect rest
    ) where {MT <: Manifold}
    (ρ_prime >= 0.25) && throw( ErrorException(
        "ρ_prime must be strictly smaller than 0.25 but it is $ρ_prime."
    ))
    (Δ_bar <= 0) && throw( ErrorException("Δ_bar must be positive but it is $Δ_bar.") )
    (Δ <= 0 || Δ > Δ_bar) &&  throw( ErrorException(
        "Δ must be positive and smaller than Δ_bar (=$Δ_bar) but it is $Δ."
    ))
    p = HessianProblem(M,F,∇F,H,preconditioner)
    o = TrustRegionsOptions(
        x,
        stopping_criterion,
        Δ,
        Δ_bar,
        retraction,
        useRandom,
        ρ_prime,
        ρ_regularization
    )
    o = decorate_options(o; kwargs...)
    resultO = solve(p,o)
    if return_options
        return resultO
    else
        return get_solver_result(resultO)
    end
end

initialize_solver!(p::P,o::O) where {P <: HessianProblem, O <: TrustRegionsOptions} = nothing

function step_solver!(p::P,o::O,iter) where {P <: HessianProblem, O <: TrustRegionsOptions}
        # Determine eta0
        if o.useRand==false
                # Pick the zero vector
                eta = zero_tangent_vector(p.M, o.x)
        else
                # Random vector in T_x M (this has to be very small)
                eta = 10.0^(-6)*random_tangent(p.M, o.x)
                while norm(p.M, o.x, eta) > o.Δ
                        # Must be inside trust-region
                        eta = sqrt(sqrt(eps(Float64)))*eta
                end
        end
        # Solve TR subproblem approximately
        opt = truncated_conjugate_gradient_descent(p.M,p.cost,p.gradient,
        o.x,eta,p.hessian,o.Δ;preconditioner=p.precon,useRandom=o.useRand,
        #debug = [:Iteration," ",:Stop],
        return_options=true)
        option = get_options(opt) # remove decorators
        η = get_solver_result(option)
        SR = get_active_stopping_criteria(option.stop)
        Hη = getHessian(p, o.x, η)
        # Initialize the cost function F und the gradient of the cost function
        # ∇F at the point x
        grad = get_gradient(p, o.x)
        fx = get_cost(p, o.x)
        norm_grad = norm(p.M, o.x, grad)
        # If using randomized approach, compare result with the Cauchy point.
        if o.useRand
                # Check the curvature,
                Hgrad = getHessian(p, o.x, grad)
                gradHgrad = inner(p.M, o.x, grad, Hgrad)
                if gradHgrad <= 0
                        tau_c = 1
                else
                        tau_c = min( norm_grad^3 /(o.Δ * gradHgrad), 1)
                end
                # and generate the Cauchy point.
                η_c = (-tau_c * o.Δ / norm_grad) * grad
                Hη_c = (-tau_c * o.Δ / norm_grad) * Hgrad
                # Now that we have computed the Cauchy point in addition to the
                # returned eta, we might as well keep the best of them.
                mdle  = fx + inner(p.M, o.x, grad, η) + .5 * inner(p.M, o.x, Hη, η)
                mdlec = fx + inner(p.M, o.x, grad, η_c) + .5 * inner(p.M, o.x, Hη_c, η_c)
                if mdlec < mdle
                        η = η_c
                        Hη = Hη_c
                end
        end
        # Compute the tentative next iterate (the proposal)
        x_prop  = o.retraction(p.M, o.x, η)
        # Compute the function value of the proposal
        fx_prop = get_cost(p, x_prop)
        # Check the performance of the quadratic model against the actual cost.
        ρnum = fx - fx_prop
        ρden = -inner(p.M, o.x, η, grad) - 0.5*inner(p.M, o.x, η, Hη)
        # Since, at convergence, both ρnum and ρden become extremely small,
        # computing ρ is numerically challenging. The break with ρnum and ρden
        # can thus lead to a large error in rho, making the
        # acceptance / rejection erratic. Meanwhile, close to convergence,
        # steps are usually trustworthy and we should transition to a Newton-
        # like method, with rho=1 consistently. The heuristic thus shifts both
        # rhonum and rhoden by a small amount such that far from convergence,
        # the shift is irrelevant and close to convergence, the ratio rho goes
        # to 1, effectively promoting acceptance of the step.
        ρ_reg = max(1, abs(fx)) * eps(Float64) * o.ρ_regularization
        ρnum = ρnum + ρ_reg
        ρden = ρden + ρ_reg
        model_decreased = (ρden >= 0)
        ρ = ρnum / ρden
        # Choose the new TR radius based on the model performance.
        # If the actual decrease is smaller than 1/4 of the predicted decrease,
        # then reduce the TR radius.
        if ρ < 1/4 || model_decreased == false || isnan(ρ)
                o.Δ = o.Δ/4
        elseif ρ > 3/4 && any([typeof(s) in [StopWhenTrustRegionIsExceeded,StopWhenCurvatureIsNegative] for s in SR] )
                o.Δ = min(2*o.Δ, o.Δ_bar)
        else
                o.Δ = o.Δ
        end
        # Choose to accept or reject the proposed step based on the model
        # performance. Note the strict inequality.
        if model_decreased && ρ > o.ρ_prime
                o.x = x_prop
        end
        return nothing
end
get_solver_result(o::O) where {O <: TrustRegionsOptions} = o.x


text:
@testset "client" begin
    client = SearchClient("idx")
    @test client.index_name == "idx"
    @test typeof(client.client) <: Jedis.Client

    @test typeof(get_search_client()) <: SearchClient

    index_def_hash = IndexDefinition()
    @test typeof(index_def_hash) <: IndexDefinition
    @test index_def_hash.args == ["ON","HASH"]

    index_def_json = IndexDefinition(index_type=IndexType.JSON)
    @test index_def_json.args == ["ON","JSON"]

    @test_throws MethodError IndexDefinition(index_type=3)

    idx = IndexDefinition(
        prefix=["idx:"],
        language_field="language",
        language="english",
        score_field="name",
        score=1.0,
    )

    @test idx.args == ["ON", "HASH", "PREFIX", 1, "idx:", "LANGUAGE_FIELD", "language", "LANGUAGE", "english", "SCORE_FIELD", "name"]

    text_field = TextField("field1")

    @test create_index(
        [text_field];
        definition=idx,
        no_term_offsets=true,
        no_field_flags=true,
        stopwords=["of"],
        max_text_fields=true,
        no_highlight=true,
        no_term_frequencies=true,
        skip_initial_scan=true,
    ) == "OK"

    @test Jedis.execute("FT._LIST") == ["idx"]

    idx = IndexDefinition(
        prefix=["idx:"],
        score=1.0,
    )
    drop_index()
    @test Jedis.execute("FT._LIST") == []
    create_index(
        [text_field];
        definition=idx
    )
    @test Jedis.execute("FT._LIST") == ["idx"]
    drop_index()

    Jedis.execute("FLUSHALL")
end