text:
using CharFuncPricing, Test, Printf
using StatsBase

@testset "HestonLong" begin
τ=10.0
forward= 10000.0
params=HestonParams(1.0, 2.0, 0.0025, 0.5, 0.1)
strikes = [100.0001, 101.0, 110.0, 200.0, 1000.0, 10000.0]
cf = DefaultCharFunc(params)
quadratureGL = ModlobQuadrature(1e-13)
pricers = [ALCharFuncPricer(cf, quadratureGL), makeCosCharFuncPricer(cf, τ, 200, 8),
 makeCosCharFuncPricer(cf, τ, 8, tol=1e-8), FlinnCharFuncPricer(cf, τ, tTol = 1e-10, qTol=1e-10),
 FlinnCharFuncPricer(cf, τ, tTol = 1e-40, qTol=1e-10), AdaptiveFlinnCharFuncPricer(cf, τ, qTol=1e-8), ALCharFuncPricer(cf, n=200)]
pricerNames = ["Reference", "Cos (M=200)", "Cos-Adaptive", "Flinn-Truncated (1e-8)", "Flinn-Truncated (1e-40)", "Flinn-Transformed", "Andersen-Lake (n=200)"]
isCall = false
refPrices = map(x-> priceEuropean(pricers[1], isCall, x, forward, τ, 1.0),strikes)
for price in refPrices
    @printf("& %.8e", price)
end
println()
for (pricer,name) in zip(pricers,pricerNames)
    prices = map(x -> priceEuropean(pricer, isCall, x, forward, τ, 1.0), strikes)
    print(name," ")
    for (price,refPrice) in zip(prices,refPrices)
        @printf("& %.2e", price - refPrice)
    end
    println("\\\\")
end
prices = map(x -> priceEuropean(pricers[4], isCall, x, forward, τ, 1.0), strikes)
#@test isapprox(prices[1]-refPrices[1],1.97, atol=1e-2) #bad price for Flinn
@test isapprox(refPrices[1],prices[1], atol=1e-7)
prices = map(x -> priceEuropean(pricers[6], isCall, x, forward, τ, 1.0), strikes)
@test isapprox(prices[1], refPrices[1], atol=1e-7)
prices = map(x -> priceEuropean(pricers[7], isCall, x, forward, τ, 1.0), strikes)
@test isapprox(prices[1], refPrices[1], atol=1e-15)

τ=10.0
forward = 100.0
params = HestonParams(0.0001, 0.1, 0.25, 0.95, 3.0)
strike = 100.0001
cf = DefaultCharFunc(params)
pricers = [ALCharFuncPricer(cf, quadratureGL),  makeCosCharFuncPricer(cf, τ, 8),
makeCosCharFuncPricer(cf, τ, 200, 8), makeCosCharFuncPricer(cf, τ, 1000, 8),
 FlinnCharFuncPricer(cf, τ, tTol = 1e-10, qTol=1e-10),ALCharFuncPricer(cf, n=200)]
 pricerNames = ["Reference",  "Cos-Adaptive","Cos (M=200)","Cos (M=1000)",
  "Flinn-Truncated (1e-10)",  "Andersen-Lake (n=200)"]

isCall = false
refPrices = map(x-> priceEuropean(pricers[1], isCall, x, forward, τ, 1.0),strikes)
for price in refPrices
    @printf("& %.8e", price)
end
println()
for (pricer,name) in zip(pricers,pricerNames)
    prices = map(x -> priceEuropean(pricer, isCall, x, forward, τ, 1.0), strikes)
    print(name," ")
    for (price,refPrice) in zip(prices,refPrices)
        @printf("& %.2e", price - refPrice)
    end
    println("\\\\")
end
prices = map(x -> priceEuropean(pricers[4], isCall, x, forward, τ, 1.0), strikes)
@test isapprox(prices[1]-refPrices[1],0.27, atol=1e-2) #bad price for Cos
prices = map(x -> priceEuropean(pricers[2], isCall, x, forward, τ, 1.0), strikes)
@test isapprox(prices[1], refPrices[1], atol=1e-4) # bad accuracy for adaptive Cos
prices = map(x -> priceEuropean(pricers[5], isCall, x, forward, τ, 1.0), strikes)
@test isapprox(prices[1], refPrices[1], atol=1e-7)

#10.0 0.0001 0.01 1.0 -0.95 3.0 10000.0 9900.0 9882.675792888303
params = HestonParams(0.0001, 0.01, 1.0, -0.95, 3.0)
cf = DefaultCharFunc(params)
pricers=[ALCharFuncPricer(cf, quadratureGL), makeCosCharFuncPricer(cf, τ, 8),
 FlinnCharFuncPricer(cf, τ, tTol = 1e-10, qTol=1e-10), ALCharFuncPricer(cf, n=200)]
 pricerNames=["Reference", "Cos-Adaptive 8", "Flinn Trunc 1e-10", "AndersenLake"]
isCall = false
refPrices = map(x-> priceEuropean(pricers[1], isCall, x, forward, τ, 1.0),strikes)
for price in refPrices
    @printf("& %.8e", price)
end
println()
for (pricer,name) in zip(pricers,pricerNames)
    prices = map(x -> priceEuropean(pricer, isCall, x, forward, τ, 1.0), strikes)
    print(name," ")
    for (price,refPrice) in zip(prices,refPrices)
        @printf("& %.2e", price - refPrice)
    end
    println("\\\\")
end
prices = map(x -> priceEuropean(pricers[2], isCall, x, forward, τ, 1.0), strikes)
@test isless(1e-4, abs(prices[1]-refPrices[1])) #bad price for Cos
prices = map(x -> priceEuropean(pricers[3], isCall, x, forward, τ, 1.0), strikes)
@test isapprox(prices[1], refPrices[1], atol=1e-7)
prices = map(x -> priceEuropean(pricers[4], isCall, x, forward, τ, 1.0), strikes)
@test isapprox(prices[1], refPrices[1], atol=1e-8)

τ=10.0
strike=101.0
params=HestonParams{Float64}(0.0001, 0.01, 0.0001, 0.95, 3.0)
cf = DefaultCharFunc(params)
pricers = [ALCharFuncPricer(cf, quadratureGL), makeCosCharFuncPricer(cf, τ, 200, 8), makeCosCharFuncPricer(cf, τ, 8, tol=1e-8),
FlinnCharFuncPricer(cf, τ, tTol = 1e-10, qTol=1e-10), FlinnCharFuncPricer(cf, τ, tTol = 1e-40, qTol=1e-10),
 AdaptiveFlinnCharFuncPricer(cf, τ, qTol=1e-8), ALCharFuncPricer(cf, n=200)]
pricerNames = ["Reference", "Cos (M=200)", "Cos-Adaptive", "Flinn-Truncated (1e-8)", "Flinn-Truncated (1e-40)", "Flinn-Transformed", "Andersen-Lake (n=200)"]
isCall = false
#bad price for adapt-cos vs cos 8
refPrices = map(x-> priceEuropean(pricers[1], isCall, x, forward, τ, 1.0),strikes)
for price in refPrices
    @printf("& %.8e", price)
end
println()
for (pricer,name) in zip(pricers,pricerNames)
    prices = map(x -> priceEuropean(pricer, isCall, x, forward, τ, 1.0), strikes)
    print(name," ")
    for (price,refPrice) in zip(prices,refPrices)
        @printf("& %.2e", price - refPrice)
    end
    println("\\\\")
    # println(strike," & ",priceC," & ",priceC2," & ",priceF," & ",priceF2, " & ",priceA," \\\\")
end
prices = map(x -> priceEuropean(pricers[2], isCall, x, forward, τ, 1.0), strikes)
@test isless(6.28, abs(prices[1]-refPrices[1])) #bad price for Cos
prices = map(x -> priceEuropean(pricers[3], isCall, x, forward, τ, 1.0), strikes)
@test isapprox(prices[1], refPrices[1], atol=1e-4)
prices = map(x -> priceEuropean(pricers[4], isCall, x, forward, τ, 1.0), strikes)
#@test isless(1.86, abs(prices[1]-refPrices[1])) #bad price for Flinn-truncated
@test isapprox(prices[1], refPrices[1], atol=1e-7)
prices = map(x -> priceEuropean(pricers[6], isCall, x, forward, τ, 1.0), strikes)
@test isapprox(prices[1], refPrices[1], atol=1e-7)
prices = map(x -> priceEuropean(pricers[7], isCall, x, forward, τ, 1.0), strikes)
@test isapprox(prices[1], refPrices[1], atol=1e-8)

τ=10.0
strike=101.0
params= HestonParams{Float64}(1.0, 0.01, 0.25, 0.95, 3.0)
cf = DefaultCharFunc(params)
pricers= [ ALCharFuncPricer(cf, quadratureGL),  makeCosCharFuncPricer(cf, τ, 8),
FlinnCharFuncPricer(HestonCVCharFunc(cf), τ, tTol = 1e-10, qTol=1e-10), ALCharFuncPricer(cf, n=200)]
pricerNames=["Reference","Cos Adapt","Flinn trunc", "AL"]
isCall = false
refPrices = map(x-> priceEuropean(pricers[1], isCall, x, forward, τ, 1.0),strikes)
for price in refPrices
    @printf("& %.8e", price)
end
println()
for (pricer,name) in zip(pricers,pricerNames)
    prices = map(x -> priceEuropean(pricer, isCall, x, forward, τ, 1.0), strikes)
    print(name," ")
    for (price,refPrice) in zip(prices,refPrices)
        @printf("& %.2e", price - refPrice)
    end
    println("\\\\")
    # println(strike," & ",priceC," & ",priceC2," & ",priceF," & ",priceF2, " & ",priceA," \\\\")
end
prices = map(x -> priceEuropean(pricers[3], isCall, x, forward, τ, 1.0), strikes)
@test isapprox(prices[1], refPrices[1], atol=1e-8)
prices = map(x -> priceEuropean(pricers[4], isCall, x, forward, τ, 1.0), strikes)
@test isapprox(prices[1], refPrices[1], atol=1e-8)

#(30.0, 100.0001, HestonParams{Float64}(1.0, 0.1, 1.0, 0.95, 0.1))
τ=30.0
params =  HestonParams{Float64}(1.0, 0.1, 1.0, 0.95, 0.1)
cf = DefaultCharFunc(params)
pricers= [ ALCharFuncPricer(cf, quadratureGL),  makeCosCharFuncPricer(cf, τ, 8),
FlinnCharFuncPricer(HestonCVCharFunc(cf), τ, tTol = 1e-10, qTol=1e-10), ALCharFuncPricer(cf, n=200)]
pricerNames=["Reference","Cos Adapt","Flinn trunc", "AL"]
isCall = false
refPrices = map(x-> priceEuropean(pricers[1], isCall, x, forward, τ, 1.0),strikes)
for price in refPrices
    @printf("& %.8e", price)
end
println()
for (pricer,name) in zip(pricers,pricerNames)
    prices = map(x -> priceEuropean(pricer, isCall, x, forward, τ, 1.0), strikes)
    print(name," ")
    for (price,refPrice) in zip(prices,refPrices)
        @printf("& %.2e", price - refPrice)
    end
    println("\\\\")
    # println(strike," & ",priceC," & ",priceC2," & ",priceF," & ",priceF2, " & ",priceA," \\\\")
end
prices = map(x -> priceEuropean(pricers[3], isCall, x, forward, τ, 1.0), strikes)
@test isapprox(prices[1], refPrices[1], atol=1e-8)
prices = map(x -> priceEuropean(pricers[4], isCall, x, forward, τ, 1.0), strikes)
@test isapprox(prices[1], refPrices[1], atol=1e-8)

τ=0.0025
params = HestonParams{Float64}(0.0001, 0.01, 0.0001, -0.95, 3.0)
cf = DefaultCharFunc(params)
pricers= [ ALCharFuncPricer(cf, quadratureGL),  makeCosCharFuncPricer(cf, τ, 8),
FlinnCharFuncPricer(HestonCVCharFunc(cf), τ, tTol = 1e-10, qTol=1e-10), ALCharFuncPricer(cf, n=10000)]
pricerNames=["Reference","Cos Adapt","Flinn trunc", "AL"]
refPrices = map(x-> priceEuropean(pricers[1], isCall, x, forward, τ, 1.0),strikes)
for price in refPrices
    @printf("& %.8e", price)
end
println()
for (pricer,name) in zip(pricers,pricerNames)
    prices = map(x -> priceEuropean(pricer, isCall, x, forward, τ, 1.0), strikes)
    print(name," ")
    for (price,refPrice) in zip(prices,refPrices)
        @printf("& %.2e", price - refPrice)
    end
    println("\\\\")
    # println(strike," & ",priceC," & ",priceC2," & ",priceF," & ",priceF2, " & ",priceA," \\\\")
end
prices = map(x -> priceEuropean(pricers[3], isCall, x, forward, τ, 1.0), strikes)
@test isapprox(prices[1], refPrices[1], atol=1e-10)
prices = map(x -> priceEuropean(pricers[4], isCall, x, forward, τ, 1.0), strikes)
@test isapprox(prices[1], refPrices[1], atol=1e-15)

#flinn breaks strongly with or without cv. Also due to bad truncatino range (better if we use short mat criteria as well)
τ=0.5
strike = 100.0
params = HestonParams{Float64}(0.0025, 2.0, 0.25, 0.5, 0.0001)
cf = DefaultCharFunc(params)
pricers= [ ALCharFuncPricer(cf, quadratureGL),  makeCosCharFuncPricer(cf, τ, 8),
FlinnCharFuncPricer(HestonCVCharFunc(cf), τ, tTol = 1e-10, qTol=1e-10), ALCharFuncPricer(cf, n=10000)]
pricerNames=["Reference","Cos Adapt","Flinn trunc", "AL"]
refPrices = map(x-> priceEuropean(pricers[1], isCall, x, forward, τ, 1.0),strikes)
for price in refPrices
    @printf("& %.8e", price)
end
println()
for (pricer,name) in zip(pricers,pricerNames)
    prices = map(x -> priceEuropean(pricer, isCall, x, forward, τ, 1.0), strikes)
    print(name," ")
    for (price,refPrice) in zip(prices,refPrices)
        @printf("& %.2e", price - refPrice)
    end
    println("\\\\")
    # println(strike," & ",priceC," & ",priceC2," & ",priceF," & ",priceF2, " & ",priceA," \\\\")
end
prices = map(x -> priceEuropean(pricers[3], isCall, x, forward, τ, 1.0), strikes)
@test isapprox(prices[1], refPrices[1], atol=1e-7)
prices = map(x -> priceEuropean(pricers[4], isCall, x, forward, τ, 1.0), strikes)
@test isapprox(prices[1], refPrices[1], atol=1e-13)

#DE adaptive largest error abs
τ= 0.1
strike= 101.0
forward=100.0
params = HestonParams{Float64}(0.0001, 2.0, 0.04, -0.5, 3.0)
cf = DefaultCharFunc(params)
pricers= [ ALCharFuncPricer(cf, quadratureGL),  makeCosCharFuncPricer(cf, τ, 8),
FlinnCharFuncPricer(HestonCVCharFunc(cf), τ, tTol = 1e-10, qTol=1e-10), ALCharFuncPricer(cf, n=10000)]
pricerNames=["Reference","Cos Adapt","Flinn trunc", "AL"]
refPrices = map(x-> priceEuropean(pricers[1], isCall, x, forward, τ, 1.0),strikes)
for price in refPrices
    @printf("& %.8e", price)
end
println()
for (pricer,name) in zip(pricers,pricerNames)
    prices = map(x -> priceEuropean(pricer, isCall, x, forward, τ, 1.0), strikes)
    print(name," ")
    for (price,refPrice) in zip(prices,refPrices)
        @printf("& %.2e", price - refPrice)
    end
    println("\\\\")
    # println(strike," & ",priceC," & ",priceC2," & ",priceF," & ",priceF2, " & ",priceA," \\\\")
end
prices = map(x -> priceEuropean(pricers[3], isCall, x, forward, τ, 1.0), strikes)
@test isapprox(prices[1], refPrices[1], atol=1e-9)
prices = map(x -> priceEuropean(pricers[4], isCall, x, forward, τ, 1.0), strikes)
@test isapprox(prices[1], refPrices[1], atol=1e-14)

#Interestingly, AL-10000 is likley wrong as AL-0 with tol=1e-20 = TFLinn price, strike = 101 & 200.
τ=30.0
forward = 100.0
params = HestonParams{Float64}(0.0025, 0.1, 0.0001, 0.1, 3.0)
cf = DefaultCharFunc(params)
pricers= [ ALCharFuncPricer(cf, quadratureGL),  makeCosCharFuncPricer(cf, τ, 8),
FlinnCharFuncPricer(HestonCVCharFunc(cf), τ, tTol = 1e-10, qTol=1e-10), ALCharFuncPricer(cf, n=10000)]
pricerNames=["Reference","Cos Adapt","Flinn trunc", "AL"]
refPrices = map(x-> priceEuropean(pricers[1], isCall, x, forward, τ, 1.0),strikes)
for price in refPrices
    @printf("& %.8e", price)
end
println()
for (pricer,name) in zip(pricers,pricerNames)
    prices = map(x -> priceEuropean(pricer, isCall, x, forward, τ, 1.0), strikes)
    print(name," ")
    for (price,refPrice) in zip(prices,refPrices)
        @printf("& %.2e", price - refPrice)
    end
    println("\\\\")
    # println(strike," & ",priceC," & ",priceC2," & ",priceF," & ",priceF2, " & ",priceA," \\\\")
end
prices = map(x -> priceEuropean(pricers[3], isCall, x, forward, τ, 1.0), strikes)
@test isapprox(prices[1], refPrices[1], atol=1e-8)
prices = map(x -> priceEuropean(pricers[4], isCall, x, forward, τ, 1.0), strikes)
@test isapprox(prices[1], refPrices[1], atol=1e-12)

#Largest error on AL-1000 vs 2000
τ=10.0
forward = 100.0
params = HestonParams{Float64}(0.0001, 0.01, 0.0001, 0.95, 1.0)
#(10.0, 100.0001, HestonParams{Float64}(0.0001, 0.01, 0.0001, 0.95, 1.0))
cf = DefaultCharFunc(params)
pricers = [ALCharFuncPricer(cf, quadratureGL), makeCosCharFuncPricer(cf, τ, 200, 8),
makeCosCharFuncPricer(cf, τ, 8, tol=1e-8), FlinnCharFuncPricer(cf, τ, tTol = 1e-10, qTol=1e-10),
 FlinnCharFuncPricer(cf, τ, tTol = 1e-40, qTol=1e-10), AdaptiveFlinnCharFuncPricer(cf, τ, qTol=1e-8), ALCharFuncPricer(cf, n=200)]
pricerNames = ["Reference", "Cos (M=200)", "Cos-Adaptive", "Flinn-Truncated (1e-8)", "Flinn-Truncated (1e-40)", "Flinn-Transformed", "Andersen-Lake (n=200)"]
isCall = false
#bad price for adapt-cos vs cos 8
refPrices = map(x-> priceEuropean(pricers[1], isCall, x, forward, τ, 1.0),strikes)
for price in refPrices
    @printf("& %.8e", price)
end
println()
for (pricer,name) in zip(pricers,pricerNames)
    prices = map(x -> priceEuropean(pricer, isCall, x, forward, τ, 1.0), strikes)
    print(name," ")
    for (price,refPrice) in zip(prices,refPrices)
        @printf("& %.2e", price - refPrice)
    end
    println("\\\\")
    # println(strike," & ",priceC," & ",priceC2," & ",priceF," & ",priceF2, " & ",priceA," \\\\")
end
prices = map(x -> priceEuropean(pricers[2], isCall, x, forward, τ, 1.0), strikes)
@test isless(0.95, abs(prices[1]-refPrices[1])) #bad price for Cos
prices = map(x -> priceEuropean(pricers[3], isCall, x, forward, τ, 1.0), strikes)
@test isless(2e-4, abs(prices[1]-refPrices[1])) #bad price for Cos
prices = map(x -> priceEuropean(pricers[4], isCall, x, forward, τ, 1.0), strikes)
#@test isless(5e-2, abs(prices[1]-refPrices[1])) #bad price for Flinn
@test isapprox(prices[1], refPrices[1], atol=1e-7)
prices = map(x -> priceEuropean(pricers[5], isCall, x, forward, τ, 1.0), strikes)
#@test isless(5e-2, abs(prices[1]-refPrices[1])) #bad price for Flinn
@test isapprox(prices[1], refPrices[1], atol=1e-7)
prices = map(x -> priceEuropean(pricers[6], isCall, x, forward, τ, 1.0), strikes)
@test isapprox(prices[1], refPrices[1], atol=1e-7)
prices = map(x -> priceEuropean(pricers[7], isCall, x, forward, τ, 1.0), strikes)
@test isapprox(prices[1], refPrices[1], atol=1e-8)

#largest error n=200
τ=2.0
strike = 100.0
forward = 200.0
forwards = [100.0,100.0001, 101.0, 110.0, 200.0, 1000.0, 10000.0]
params = HestonParams{Float64}(0.0001, 2.0, 0.0001, 0.95, 1.0)
cf = DefaultCharFunc(params)
pricers = [ALCharFuncPricer(cf, quadratureGL), makeCosCharFuncPricer(cf, τ, 200, 8),
 makeCosCharFuncPricer(cf, τ, 8, tol=1e-8), FlinnCharFuncPricer(cf, τ, tTol = 1e-10, qTol=1e-10),
  FlinnCharFuncPricer(cf, τ, tTol = 1e-40, qTol=1e-10), AdaptiveFlinnCharFuncPricer(cf, τ, qTol=1e-8), ALCharFuncPricer(cf, n=200)]
pricerNames = ["Reference", "Cos (M=200)", "Cos-Adaptive", "Flinn-Truncated (1e-8)", "Flinn-Truncated (1e-40)", "Flinn-Transformed", "Andersen-Lake (n=200)"]
isCall = false
#bad price for adapt-cos vs cos 8
refPrices = map(x-> priceEuropean(pricers[1], isCall, strike, x, τ, 1.0),forwards)
for price in refPrices
    @printf("& %.8e", price)
end
println()
for (pricer,name) in zip(pricers,pricerNames)
    prices = map(x -> priceEuropean(pricer, isCall, strike, x, τ, 1.0), forwards)
    print(name," ")
    for (price,refPrice) in zip(prices,refPrices)
        @printf("& %.2e", price - refPrice)
    end
    println("\\\\")
end
prices = map(x -> priceEuropean(pricers[7], isCall, strike, x, τ, 1.0), forwards)
@test isapprox(prices[5], refPrices[5], atol=1.2e-2)
@test isless(0.01, prices[5]/ refPrices[5])

end

@testset "HestonFlinn" begin
    v0 = 0.0718
    ρ = -0.352
    σ = 0.582
    θ = 0.0762
    κ = 1.542
    params = HestonParams(v0, κ, θ, ρ, σ)

    tol = 1e-8
    τ = 0.1
    forward = 6317.80
    strike = 6317.80
    cf = DefaultCharFunc(params)
    refPricer = makeCosCharFuncPricer(cf, τ, 2048 * 4, 32)

    pricer = FlinnCharFuncPricer(cf, τ, qTol = tol, tTol = tol)
    refPrice = priceEuropean(refPricer, false, strike, forward, τ, 1.0)
    price = priceEuropean(pricer, false, strike, forward, τ, 1.0)
    @test isapprox(refPrice, price, atol = 5e-4)
#    @test isapprox(209.820637, price, atol = 1e-5)

    ccf = HestonCVCharFunc(cf)
    pricer = FlinnCharFuncPricer(ccf, τ, qTol = tol, tTol = tol)
    price = priceEuropean(pricer, false, strike, forward, τ, 1.0)
    @test isapprox(refPrice, price, atol = 1e-6)

    #Make sure BigFloat works
    params = HestonParams(BigFloat("0.0718"), BigFloat("1.542"), BigFloat("0.0762"),BigFloat("-0.352"),BigFloat("0.582"))
    cf = DefaultCharFunc{HestonParams{BigFloat},Complex{BigFloat}}(params)
    τ = BigFloat("0.1")
    strike = BigFloat("6317.80")
    forward = BigFloat("6317.80")
    pricer = AdaptiveFlinnCharFuncPricer(cf, τ, qTol = BigFloat(1e-24))
    priceB = priceEuropean(pricer, false, strike, forward, τ, BigFloat(1.0))
    @test isapprox(refPrice, Float64(priceB), atol = 1e-12)
    pricer = FlinnCharFuncPricer(cf, τ, qTol = BigFloat(1e-24), tTol = BigFloat(1e-24))
    priceC = priceEuropean(pricer, false, strike, forward, τ, BigFloat(1.0))
    @test isapprox(priceB, priceC, atol=BigFloat(1e-18))
end

@testset "HestonLogContinuity" begin
    v0 = 0.08
    ρ = -0.8
    σ = 0.25
    θ = 0.1
    κ = 3.0
    #params from Cui et al for Heston.
    params = HestonParams(v0, κ, θ, ρ, σ)
    #with the original formula from SZ, there is a discontinuity in the log.
    cf = DefaultCharFunc(params)
    sz = CharFuncPricing.evaluateLogCharFuncCui(cf, 2.0 + 0.0im, 15.0)
    lk = CharFuncPricing.evaluateLogCharFunc(cf, 2.0 + 0.0im, 15.0)
    @test isapprox(lk, sz, atol = 1e-13)
    for u = 0.1:0.1:4.0
        lk = CharFuncPricing.evaluateLogCharFunc(cf, u + 0.0im, 15.0)
        sz = CharFuncPricing.evaluateLogCharFuncCui(cf, u + 0.0im, 15.0)
        @test isapprox(lk, sz, atol = 1e-13)
    end
    v0 = 0.010201
    ρ = -0.7
    σ = 0.61
    κ = 6.21
    θ = 0.019
    α = 3.35861
    params = HestonParams{Float64}(sqrt(v0), κ, sqrt(θ), ρ, σ)
    cf = DefaultCharFunc(params)
    sz = CharFuncPricing.evaluateLogCharFunc(cf, 10.0 - (α + 1) * 1im, 10.0)
    lk = CharFuncPricing.evaluateLogCharFuncCui(cf, 10.0 - (α + 1) * 1im, 10.0)
    @test isapprox(lk, sz, atol = 1e-13)
end

@testset "LF2y" begin
    r = 0.0
    q = 0.0
    κ = 1.0
    θ = 0.1
    σ = 1.0
    ρ = -0.9
    v0 = 0.1
    τs = [2.0/52, 2.0/12, 2.0]
    spot = 1.0

    τ = τs[3]
    spot *= exp((r - q) * τ)
    df = exp(-r * τ)
    params = HestonParams{Float64}(v0, κ, θ, ρ, σ)
    cf = DefaultCharFunc(params)
    m = 200
    l = 8
    refPricer = makeCosCharFuncPricer(cf, τ, 2048 * 8, 32)
    pricer = makeCosCharFuncPricer(cf, τ, m, l)
    refPrices = Vector{Float64}(undef, 0)
    prices = Vector{Float64}(undef, 0)
    for strike = 0.4:0.012:1.6
        refPrice = priceEuropean(refPricer, false, strike, spot,τ, df)
        price = priceEuropean(pricer, false, strike, spot, τ, df)
        push!(refPrices, refPrice)
        push!(prices, price)
    end
    rmse = rmsd(prices, refPrices)
    println("RMSE ", rmse)
    @test isless(rmse, 1e-4)
    @test isless(1e-5, rmse)

    for τ in τs
        local refPrices = zeros(length(collect(0.4:0.012:1.6)))
        local prices =  zeros(length(collect(0.4:0.012:1.6)))
        measuredTime = @elapsed begin
        pricer = makeCosCharFuncPricer(cf, τ, 2048 * 8, 32)
        for (i, strike) in enumerate(0.4:0.012:1.6)
        price = priceEuropean(pricer, false, strike, spot,  τ, df)
        refPrices[i] = price
    end
    end
    println("Ref ", measuredTime)

    measuredTime = @elapsed begin
        pricer = makeCosCharFuncPricer(cf, τ, 200,8)
        for (i, strike) in enumerate(0.4:0.012:1.6)
        price = priceEuropean(pricer, false, strike, spot, τ, df)
        prices[i] = price
    end
    end
    rmse = rmsd(prices, refPrices)
    println("Cos ", rmse, " ",measuredTime)

    measuredTime = @elapsed begin
        pricer = FlinnCharFuncPricer(cf, τ)
        for (i, strike) in enumerate(0.4:0.012:1.6)
        price = priceEuropean(pricer, false, strike, spot,  τ, df)
        prices[i] = price
    end
    end
    rmse = rmsd(prices, refPrices)
    println("Truncated Flinn ", rmse, " ",measuredTime)

    measuredTime = @elapsed begin
        pricer = AdaptiveFlinnCharFuncPricer(cf, τ)
        for (i, strike) in enumerate(0.4:0.012:1.6)
        price = priceEuropean(pricer, false, strike, spot, τ, df)
        prices[i] = price
    end
    end
    rmse = rmsd(prices, refPrices)
    println("Transformed Flinn ", rmse, " ",measuredTime)

quadrature = TanhSinhQuadrature(200, eps())
    measuredTime = @elapsed begin
        pricer = ALCharFuncPricer(cf, quadrature)
        for (i, strike) in enumerate(0.4:0.012:1.6)
        price = priceEuropean(pricer, false, strike, spot, τ, df)
        prices[i] = price
        end
    end
    rmse = rmsd(prices, refPrices)
    println("Andersen-Lake ", rmse, " ",measuredTime)

    println(τ)
end
    κ = 1.417
    θ = 0.082
    σ = 0.591
    ρ = -0.384
    v0 = 0.073
end
#
# @testset "LFShortA" begin
#     r = 0.0
#     q = 0.0
#     κ = 0.254
#     θ = 0.320
#     σ = 0.344
#     ρ = -0.557
#     v0 = 0.826
#     τ = 0.0182
#     spot = 1000.0
#
#     spot *= exp((r - q) * τ)
#     df = exp(-r * τ)
#     params = HestonParams(v0, κ, θ, ρ, σ)
#     m = 200
#     l = 8
#     refPricer = makeCosCharFuncPricer(
#         Complex,
#         Float64,
#
#         params,
#         τ,
#         2048*4,
#         32,
#     )
#     pricer =
#         makeCosCharFuncPricer(Complex, Float64,  params, τ, m, l)
#     strike = 1400.0
#         refPrice = priceEuropean(refPricer, true, strike, spot, df)
#         price = priceEuropean(pricer, true, strike, spot, df)
#         println(
#             strike,
#             " ",
#             price,
#             " ",
#             refPrice,
#             " ",
#             price - refPrice,
#             " ",
#             price / refPrice - 1,
#         )
#
#     cinf = (params.v0+params.κ*params.θ*τ)/params.σ *sqrt(1-params.ρ^2)
#     #exp(-cinf*umax) = eps
#     #2*exp(-cinf*mmax*pi/(b-a))/(mmax*pi) = eps => exp(cinf*mmax*pi/(b-a))*mmax*pi/(b-a)*cinf = 2/eps/(b-a)*cinf => cinf*mmax*pi/(b-a) = lambertW( 2/eps/(b-a)*cinf)
#     eps = 1e-8
#     umax = -log(eps)/cinf
#     mmax = ceil(Int,umax*(pricer.b-pricer.a)/pi)+1
#     mmax = ceil(Int, sqrt(   lambertW(2*v0*τ/(eps*eps*(pricer.b-pricer.a)))/(v0*τ*pi)*(pricer.b-pricer.a)) )
#     mmax = max(32, ceil(Int, lambertW(2*cinf/(eps*(pricer.b-pricer.a)))/(cinf*pi)*(pricer.b-pricer.a)))
#     pricer =
#         makeCosCharFuncPricer(Complex, Float64,  params, τ, mmax, l)
#             price = priceEuropean(pricer, true, strike, spot, df)
#             println(
#                 strike,
#                 " ",
#                 price,
#                 " ",
#                 refPrice,
#                 " ",
#                 price - refPrice,
#                 " ",
#                 price / refPrice - 1,
#             )
# end
#
#
# @testset "LF1Y" begin
#     r = 0.0
#     q = 0.0
#     κ = 0.1
#     θ = 0.01
#     σ = 2.0
#     ρ = 0.5
#     v0 = 0.0225
#     τ = 1.0
#     spot = 1000000.0
#
#     spot *= exp((r - q) * τ)
#     df = exp(-r * τ)
#     params = HestonParams(v0, κ, θ, ρ, σ)
#     m = 200
#     l = 10
#     refPricer = makeCosCharFuncPricer(
#         Complex,
#         Float64,
#
#         params,
#         τ,
#         2048*4,
#         32,
#     )
#     pricer =
#         makeCosCharFuncPricer(Complex, Float64,  params, τ, m, l)
#     strike = spot*0.25
#         refPrice = priceEuropean(refPricer, false, strike, spot, df)
#         price = priceEuropean(pricer, false, strike, spot, df)
#         println(
#             strike,
#             " ",
#             price,
#             " ",
#             refPrice,
#             " ",
#             price - refPrice,
#             " ",
#             price / refPrice - 1,
#         )
#
#     cinf = (params.v0+params.κ*params.θ*τ)/params.σ *sqrt(1-params.ρ^2)
#     #exp(-cinf*umax) = eps
#     #2*exp(-cinf*mmax*pi/(b-a))/(mmax*pi) = eps => exp(cinf*mmax*pi/(b-a))*mmax*pi/(b-a)*cinf = 2/eps/(b-a)*cinf => cinf*mmax*pi/(b-a) = lambertW( 2/eps/(b-a)*cinf)
#     eps = 1e-10
#     umax = -log(eps)/cinf
#     mmax = ceil(Int,umax*(pricer.b-pricer.a)/pi)+1
#     mmax = ceil(Int, sqrt(   lambertW(2*v0*τ/(eps*eps*(pricer.b-pricer.a)))/(v0*τ*pi)*(pricer.b-pricer.a)) )
#     mmax = max(32, ceil(Int, lambertW(2*cinf/(eps*(pricer.b-pricer.a)))/(cinf*pi)*(pricer.b-pricer.a)))
#     pricer =
#         makeCosCharFuncPricer(Complex, Float64,  params, τ, mmax, l)
#             price = priceEuropean(pricer, false, strike, spot, df)
#             println(
#                 strike,
#                 " ",
#                 price,
#                 " ",
#                 refPrice,
#                 " ",
#                 price - refPrice,
#                 " ",
#                 price / refPrice - 1,
#             )
# end
@testset "LFShort" begin
    r = 0.0
    q = 0.0
    κ = 1.0
    θ = 0.1
    σ = 1.0
    ρ = -0.9
    v0 = 0.1
    τ = 2.0 / 365
    spot = 1.0

    spot *= exp((r - q) * τ)
    df = exp(-r * τ)
    params = HestonParams{Float64}(v0, κ, θ, ρ, σ)
    cf = DefaultCharFunc(params)
    m = 256
    l = 12
    refPricer = makeCosCharFuncPricer(cf, τ, 1024, 32)
    pricer = makeCosCharFuncPricer(cf, τ, m, l)
    for strike = 1.0:0.025:1.5
        refPrice = priceEuropean(refPricer, true, strike, spot,  τ,df)
        price = priceEuropean(pricer, true, strike, spot, τ, df)
        println(
            strike,
            " ",
            price,
            " ",
            refPrice,
            " ",
            price - refPrice,
            " ",
            price / refPrice - 1,
        )
        @test isapprox(refPrice, price, atol = 1e-15)
    end
end

#Testing https://financepress.com/2019/02/15/heston-model-reference-prices/
@testset "AlanFloat64Set" begin
    r = 0.01
    q = 0.02
    κ = 4.0
    θ = 0.25
    σ = 1.0
    ρ = -0.5
    v0 = 0.04
    τ = 1.0
    spot = 100.0
    strikes = [80.0, 90.0, 100.0, 110.0, 120.0]
    alanPuts = [
        7.958878113256768285213263077598987193482161301733,
        12.017966707346304987709573290236471654992071308187,
        17.055270961270109413522653999411000974895436309183,
        23.017825898442800538908781834822560777763225722188,
        29.811026202682471843340682293165857439167301370697,
    ]
    alanCalls = [
        26.774758743998854221382195325726949201687074848341,
        20.933349000596710388139445766564068085476194042256,
        16.070154917028834278213466703938231827658768230714,
        12.132211516709844867860534767549426052805766831181,
        9.024913483457835636553375454092357136489051667150,
    ]
    spot *= exp((r - q) * τ)
    df = exp(-r * τ)
    params = HestonParams{Float64}(v0, κ, θ, ρ, σ)
    cf = DefaultCharFunc(params)
    l = 32
    m = 1024
    pricer = makeCosCharFuncPricer(cf, τ, m, l)
    for (strike, refCall, refPut) in zip(strikes, alanCalls, alanPuts)
        price = priceEuropean(pricer, false, strike, spot, τ, df)
        println(Float64(strike), " P ", price, " ", price - refPut)
        @test isapprox(Float64(price - refPut), 0, atol = 1e-13)
        price = priceEuropean(pricer, true, strike, spot,  τ,df)
        println(Float64(strike), " C ", price, " ", price - refCall)
        @test isapprox(Float64(price - refCall), 0, atol = 1e-13)
    end
end
@testset "AlanJoshiFloat64Set" begin
    r = 0.01
    q = 0.02
    κ = 4.0
    θ = 0.25
    σ = 1.0
    ρ = -0.5
    v0 = 0.01
    τ = 0.01
    spot = 100.0
    strikes = [90.0, 95.0, 100.0, 105.0, 110.0]
    alanPuts = [
        4.5183603586861772614990106188215872180542e-8,
        0.000461954855653851579672612557018857858641926937,
        0.477781171629504680023239655436072890669645669297,
        5.009501052563650299130635110520904481889436667608,
        10.008998550115123724684210555728039829315964456261,
    ]
    alanCalls = [
        9.989001595065276544935948045293485530832966049263,
        4.989963479738160122154264702582719627807098780529,
        0.467782671512844263098248405184095087949465507760,
        2.527447823194706060519991248106500619490942e-6,
        1.29932760052624920704881258510264466e-13,
    ]
    spot *= exp((r - q) * τ)
    df = exp(-r * τ)
    params = HestonParams{Float64}(v0, κ, θ, ρ, σ)
    cf = DefaultCharFunc(params)
    l = 32
    m = 1024
    pricer = makeCosCharFuncPricer(cf, τ, m, l)
    for (strike, refCall, refPut) in zip(strikes, alanCalls, alanPuts)
        price = priceEuropean(pricer, false, strike, spot, τ, df)
        println(Float64(strike), " P ", price, " ", price - refPut)
        @test isapprox(Float64(price - refPut), 0, atol = 2e-14)
        price = priceEuropean(pricer, true, strike, spot, τ,  df)
        println(Float64(strike), " C ", price, " ", price - refCall)
        @test isapprox(Float64(price - refCall), 0, atol = 2e-14)
    end
end


text:
immutable VonMisesSampler <: Sampleable{Univariate,Continuous}
    μ::Float64
    κ::Float64
    r::Float64
end

# algorithm from
#     DJ Best & NI Fisher (1979). Efficient Simulation of the von Mises
#     Distribution. Journal of the Royal Statistical Society. Series C
#     (Applied Statistics), 28(2), 152-157.
function rand(s::VonMisesSampler)
    f = 0.0
    while true
        t, u = 0.0, 0.0
        while true
            const v, w = rand() - 0.5, rand() - 0.5
            const d, e = v ^ 2, w ^ 2
            if d + e <= 0.25
                t = d / e
                u = 4 * (d + e)
                break
            end
        end
        const z = (1.0 - t) / (1.0 + t)
        f = (1.0 + s.r * z) / (s.r + z)
        const c = s.κ * (s.r - f)
        if c * (2.0 - c) > u || log(c / u) + 1 >= c
            break
        end
    end
    mod(s.μ + (rand() > 0.5 ? acos(f) : -acos(f)), twoπ)
end

immutable VonMisesNormalApproxSampler <: Sampleable{Univariate,Continuous}
    μ::Float64
    σ::Float64
end

rand(s::VonMisesNormalApproxSampler) = mod(s.μ + s.σ * randn(), twoπ)

# normal approximation for large concentrations
VonMisesSampler(μ::Float64, κ::Float64) =
    κ > 700.0 ? VonMisesNormalApproxSampler(μ, sqrt(1.0 / κ)) :
                  begin
                      τ = 1.0 + sqrt(1.0 + 4 * κ ^ 2)
                      ρ = (τ - sqrt(2.0 * τ)) / (2.0 * κ)
                      VonMisesSampler(μ, κ, (1.0 + ρ ^ 2) / (2.0 * ρ))
                  end

ext: jl


text:
using Pkg
Pkg.activate(".")

using GraphPlot, SimpleHypergraphs
import LightGraphs

# In Atom in order to run a line press Ctrl+Enter

# It is also possible to run a block of code
# Simply select the block and press Ctrl+Enter

h = Hypergraph{Float64}(5,4)
h[1:3,1] .= 1.5
h[3,4] = 2.5
h[2,3] = 3.5
h[4,3:4] .= 4.5
h[5,4] = 5.5
h[5,2] = 6.5

display(h)    #you could also type h in the console and press enter

# manipulating the graph structure
add_vertex!(h)
add_hyperedge!(h)
h[5,5] = 1.2
h[6,5] = 1.3

display(h)

#Bipartite representation of a hypergraph h

b = BipartiteView(h)
# note that it is a view and hence
#the data is not actually copied.

# this plot will get shown in Atom as well as in Jupyter
gplot(b,
    nodefillc=vcat(fill("green", nhv(h)), fill("orange",nhe(h)) ),
    nodelabel=1:LightGraphs.nv(b),
    layout=circular_layout)

#several functions of the LightGraphs API
#will work directly on the bipartite
#view b of a hypergraph h

LightGraphs.has_self_loops(b)
LightGraphs.nv(b)
LightGraphs.ne(b)

#However it is also possible to make
#a deep copy of this view.
#In this case the deep copy is just
#a "normal" graph.

b2 = LightGraphs.SimpleGraph(b)

remove_vertex!(h,6)

LightGraphs.nv(b)
LightGraphs.nv(b2)

#Another representation of a hypergraph h
# is a two-section view

t = TwoSectionView(h)

gplot(t, nodelabel=1:LightGraphs.nv(t))

sp = shortest_path(t,1,5)

edge_cols = [ e.src in sp && e.dst in sp ? "green" : "white"  for e in LightGraphs.edges(t)]

gplot(t, nodelabel=1:LightGraphs.nv(t), edgestrokec = edge_cols)


text:
import OpenCL.cl

for (func, typ) in [(:clblasSgemm, cl.CL_float),
                    (:clblasDgemm, cl.CL_double),
                    (:clblasCgemm, CL_float2),
                    (:clblasZgemm, CL_double2)]

    @eval @blasfun $func(order::clblasOrder,
                              transA::clblasTranspose, transB::clblasTranspose,
                              M::Csize_t, N::Csize_t, K::Csize_t,
                              alpha::$typ,
                              A::cl.CL_mem, offA::Csize_t, lda::Csize_t,
                              B::cl.CL_mem, offB::Csize_t, ldb::Csize_t,
                              beta::$typ,
                              C::cl.CL_mem, offC::Csize_t, ldc::Csize_t,
                              n_queues::cl.CL_uint,
                              queues::Ptr{cl.CL_command_queue},
                              n_events_in_wait_list::cl.CL_uint,
                              event_wait_list::Ptr{cl.CL_event},
                              events::Ptr{cl.CL_event})

    @eval @blasfun2 $func(order::clblasOrder,
                               transA::clblasTranspose, transB::clblasTranspose,
                               M::Csize_t, N::Csize_t, K::Csize_t,
                               alpha::$typ,
                               A::cl.CL_mem, offA::Csize_t, lda::Csize_t,
                               B::cl.CL_mem, offB::Csize_t, ldb::Csize_t,
                               beta::$typ,
                               C::cl.CL_mem, offC::Csize_t, ldc::Csize_t)

end


text:
using TestExtra

runtests()


text:
__precompile__()

module RbCall

export _refcnt, _incref, _decref, gc_guard_references

#########################################################################
# prepare

import Libdl

struct Dl_info
    dli_fname::Ptr{UInt8}
    dli_fbase::Ptr{Cvoid}
    dli_sname::Ptr{UInt8}
    dli_saddr::Ptr{Cvoid}
end

hassym(lib, sym) = Libdl.dlsym_e(lib, sym) != C_NULL

proc_handle = unsafe_load(cglobal(:jl_exe_handle, Ptr{Cvoid}))

symbols_present = false
@static if Sys.iswindows()
  # TODO support windows
else
  global symbols_present = hassym(proc_handle, :rb_define_class)
end

if !symbols_present
  # TODO load libruby dynamically
else
  @static if Sys.iswindows()
    # TODO support windows
  else
    libruby_handle = proc_handle
    # Now determine the name of the ruby library that these symbols are from
    some_address_in_libruby = Libdl.dlsym(libruby_handle, :rb_define_class)
    some_address_in_main_exe = Libdl.dlsym_e(proc_handle, Sys.isapple() ? :_mh_execute_header : :main)
    if some_address_in_main_exe == nothing
      # NOTE: we cannot get main symbol from ruby executable in travis-ci
      some_address_in_main_exe = Libdl.dlsym_e(proc_handle, :_start)
    end
    dlinfo1 = Ref{Dl_info}()
    dlinfo2 = Ref{Dl_info}()
    ccall(:dladdr, Cint, (Ptr{Cvoid}, Ptr{Dl_info}), some_address_in_libruby, dlinfo1)
    ccall(:dladdr, Cint, (Ptr{Cvoid}, Ptr{Dl_info}), some_address_in_main_exe, dlinfo2)
    if dlinfo1[].dli_fbase == dlinfo2[].dli_fbase
      const libruby = nothing
    else
      const libruby = unsafe_string(dlinfo1[].dli_fname)
    end
  end
end

#########################################################################
# basic types

@static if sizeof(Clong) == sizeof(Ptr{Cvoid})
  const VALUE = Culong
  const SVALUE = Clong
  const ID = Culong
  const Cintptr_t = Clong
  const Cuintptr_t = Culong
elseif sizeof(Clonglong) == sizeof(Ptr{Cvoid})
  const VALUE = Culonglong
  const SVALUE = Clonglong
  const ID = Culonglong
  const Cintptr_t = Clonglong
  const Cuintptr_t = Culonglong
else
  error("ruby requires sizeof(Ptr{Cvoid}) == sizeof(Clong) or sizeof(Clonglong) to be compiled.")
end

struct RBasic_struct
  flags::VALUE
  klass::VALUE
end

struct RVALUE_struct
  basic::RBasic_struct
  v1::VALUE
  v2::VALUE
  v3::VALUE
end

const RbPtr = Ptr{RVALUE_struct}

#########################################################################
# basic constants

const RUBY_FIXNUM_MAX = typemax(Clong) >> 1
const RUBY_FIXNUM_MIN = typemin(Clong) >> 1

@static if sizeof(VALUE) >= sizeof(Cdouble)
  const RUBY_USE_FLONUM = true
else
  const RUBY_USE_FLONUM = false
  error("currently RbCall.jl only supports Ruby with flonum.")
end

@static if RUBY_USE_FLONUM
  const RUBY_Qfalse = 0x00
  const RUBY_Qtrue  = 0x14
  const RUBY_Qnil   = 0x08
  const RUBY_Qundef = 0x34
  const RUBY_IMMEDIATE_MASK = 0x07
  const RUBY_FIXNUM_FLAG    = 0x01
  const RUBY_FLONUM_MASK    = 0x03
  const RUBY_FLONUM_FLAG    = 0x02
  const RUBY_SYMBOL_FLAG    = 0x0c
else
  const RUBY_Qfalse = 0
  const RUBY_Qtrue  = 2
  const RUBY_Qnil   = 4
  const RUBY_Qundef = 6
  const RUBY_IMMEDIATE_MASK = 0x03
  const RUBY_FIXNUM_FLAG    = 0x01
  const RUBY_FLONUM_MASK    = 0x00
  const RUBY_FLONUM_FLAG    = 0x02
  const RUBY_SYMBOL_FLAG    = 0x0e
end
const RUBY_SPECIAL_SHIFT  = 8

#########################################################################
# basic utilities

RB_FIXNUM_P(v::VALUE) = (x & RUBY_FIXNUM_FLAG) != 0
RB_POSFIXABLE(x::Integer) = x <= RUBY_FIXNUM_MAX
RB_NEGFIXABLE(x::Integer) = x >= RUBY_FIXNUM_MIN
RB_FIXABLE(x::Integer) = RB_POSFIXABLE(x) && RB_NEGFIXABLE(x)

RB_INT2FIX(i::Integer) = (VALUE(i) << 1) | RUBY_FIXNUM_FLAG
RB_LONG2FIX(i::Clong) = RB_INT2FIX(i)

function RB_LONG2NUM(i::Clong)
  if RB_FIXABLE(i)
    return RB_LONG2FIX(i)
  else
    return ccall(:rb_int2big, VALUE, (Cintptr_t,), Cintptr_t(i))
  end
end

function RB_ULONG2NUM(i::Culong)
  if RB_POSFIXABLE(i)
    return RB_LONG2FIX(Clong(i))
  else
    return ccall(:rb_uint2big, VALUE, (Cuintptr_t,), Cuintptr_t(i))
  end
end

RB_LL2NUM(i::Clonglong) = ccall(:rb_ll2inum, VALUE, (Clonglong,), i)
RB_ULL2NUM(i::Culonglong) = ccall(:rb_ull2inum, VALUE, (Culonglong,), i)

DBL2NUM(d::Cdouble) = ccall(:rb_float_new, VALUE, (Cdouble,), d)

HAVE_RB_DBL_COMPLEX_NEW = hassym(libruby_handle, :rb_dbl_complex_new)

#########################################################################
# conversion

include("conversion.jl")

#########################################################################
# gc

const gc_guard_references = Dict{Any,Clong}()

function _refcnt(value::Any)
  return get(gc_guard_references, value , 0)
end

function _incref(value::Any)
  gc_guard_references[value] = _refcnt(value) + 1
  return value
end

function _decref(value::Any)
  if haskey(gc_guard_references, value)
    new_count = _refcnt(value) - 1
    if new_count == 0
      delete!(gc_guard_references, value)
    else
      gc_guard_references[value] = new_count
    end
  end
end

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

end # module RbCall


text:
abstract type AbstractOrbitCoordinate{T} end

#Enable Broadcasting
Base.broadcastable(x::T) where T <: AbstractOrbitCoordinate = (x,)

struct EPRCoordinate{T} <: AbstractOrbitCoordinate{T}
    energy::T # Kinetic Energy
    pitch::T
    r::T
    z::T
    t::T
    m::Float64
    q::Int
end

function EPRCoordinate(T=Float64; amu=H2_amu, q = 1)
    z = zero(T)
    return EPRCoordinate(z, z, z, z, z, amu*mass_u, q)
end

function EPRCoordinate(energy, pitch, r, z; t = zero(z), amu=H2_amu, q = 1)
    return EPRCoordinate(promote(energy, pitch, r, z, t)..., amu*mass_u, q)
end

function EPRCoordinate(M::AbstractEquilibrium, energy, pitch, R ; amu=H2_amu, q=1, dz=0.2)
    rmaxis, zaxis = magnetic_axis(M)
    zmax = zaxis + dz
    zmin = zaxis - dz
    if (M.psi[2]-M.psi[1]) > 0.0
        res = optimize(x->M(R,x), zmin, zmax) # Minimize
    else
        res = optimize(x->-M(R,x), zmin, zmax) # Maximize
    end
    Z = Optim.minimizer(res)
    (Z == zmax || Z == zmin) && error(@sprintf("Unable to find starting Z value with dz = %.2f. Increase dz",dz))
    return EPRCoordinate(promote(energy, pitch, R, Z, 0)..., amu*mass_u, q)
end

function GCParticle(c::EPRCoordinate)
    GCParticle(c.energy,c.pitch,c.r,c.z,c.m,c.q)
end

function EPRParticle(oc::EPRCoordinate; t=oc.t)
    EPRParticle(oc.energy,oc.pitch,oc.r,oc.z,t,oc.m,oc.q)
end

function Base.show(io::IO, c::EPRCoordinate)
    println(io,typeof(c))
    @printf(io," E = %.3f keV\n",c.energy)
    @printf(io," pitch = %.3f\n",c.pitch)
    @printf(io," Rmax = %.3f m",c.r)
end

struct HamiltonianCoordinate{T} <: AbstractOrbitCoordinate{T}
    energy::T #Kinetic + Potential Energy
    mu::T
    p_phi::T
    m::Float64
    q::Int
end

function HamiltonianCoordinate(T=Float64; amu=H2_amu, q=1)
    z = zero(T)
    return HamiltonianCoordinate(z, z, z, amu*mass_u, q)
end

function HamiltonianCoordinate(energy, mu, p_phi; amu=H2_amu, q=1)
    return HamiltonianCoordinate(promote(energy, mu, p_phi)..., amu*mass_u, q)
end

function HamiltonianCoordinate(M::AbstractEquilibrium, KE, pitch, r, z, m, q)
    psi = M(r,z)
    babs = norm(Bfield(M,r,z))
    g = poloidal_current(M,psi)
    pot = electric_potential(M,psi)
    PE = 1e-3*pot

    mc2 = m*c0^2
    KE_j = e0*KE*1e3
    p_rel2 = ((KE_j + mc2)^2 - mc2^2)/(c0*c0)
    p_para = sqrt(p_rel2)*pitch*B0Ip_sign(M)
    p_perp2 = p_rel2*(1-pitch^2)

    mu = p_perp2/(2*m*babs)
    Pphi = p_para*g/babs + q*e0*psi

    return HamiltonianCoordinate(promote(KE+PE,mu,Pphi)...,m,q)
end

function HamiltonianCoordinate(M::AbstractEquilibrium, c::EPRCoordinate)
    HamiltonianCoordinate(M,c.energy,c.pitch,c.r,c.z,c.m,c.q)
end

function HamiltonianCoordinate(M::AbstractEquilibrium, c::HamiltonianCoordinate)
    return c
end

function HamiltonianCoordinate(M::AbstractEquilibrium, KE, p, r, z; amu=H2_amu, q=1)
    HamiltonianCoordinate(M, KE, p, r, z, amu*mass_u, q)
end

function HamiltonianCoordinate(M::AbstractEquilibrium, c::GCParticle)
    HamiltonianCoordinate(M,c.energy,c.pitch,c.r,c.z,c.m,c.q)
end

function Base.show(io::IO, c::HamiltonianCoordinate)
    println(io,typeof(c))
    @printf(io," E = %.3f keV\n",c.energy)
    @printf(io," μ₀ = %.3e\n",c.mu)
    @printf(io," Pᵩ = %.3e",c.p_phi)
end


text:
module SLiSeFilters
export createConfig,createWS,updateWS!,startingPointWS,optimise,residual,Fgradient,hessianApprox,points,descent,levMar,optim,lbfgsb,wise,readFilter,readWS

const R = Float64;
const T = Complex{R};

mutable struct WeightSearch
    intv::Vector{R}
    intv_free
    intv_free_count
    wght::Vector{R}
    wght_free
    wght_free_count
    intv_modifier
end

function createWS( intv,
                   intv_free::BitArray{1},
                   wght,
                   wght_free::BitArray{1},
                   intv_modifier=x->x )
    # Check vector sizes

    if !(length(intv) == length(intv_free) && length(wght) == length(wght_free) && length(intv) == length(wght))
        error("Invalid vector input.\nIt must hold that length(intv) == length(intv_free) && length(wght) == length(wght_free) && length(intv) == length(wght).")
    end

    return SLiSeFilters.WeightSearch(
        vec(intv),
        intv_free,
        length(intv[intv_free]),
        vec(wght),
        wght_free,
        length(wght[wght_free]),
        intv_modifier
    )
end

function startingPointWS( ws::SLiSeFilters.WeightSearch )
    return vcat( ws.intv[ws.intv_free],ws.wght[ws.wght_free] )
end

function updateWS!( ws::SLiSeFilters.WeightSearch, x::Vector )
    ws.intv[ws.intv_free] = x[1:ws.intv_free_count]
    ws.wght[ws.wght_free] = x[ws.intv_free_count+1:end]
    ws.intv_modifier(ws.intv)
    return SLiSeFilters.createConfig( ws.intv, ws.wght );
end

mutable struct Config
    Weights::Array{R,2}    # contains an interval a,b and the corresponding weight w in each column [a;b;w]
    Weights_h::Array{R,2}  # contains the sam as above e for the step-function h --> only in the interval [-1,1]
end

# Some predefined weights
function createConfig( name::String = "gamma" )
    # Gamma-SLISE as per paper
    if (name == "gamma")
        intgl = [ 0.95,  1.05, 1.4, 10];
        wght =  [1, .01, 10, 20];

        # Zeta-SLISE
    elseif (name == "zeta")
        intgl = [   0.95,    0.995,    1.005,     1.05,     1.1,   1.2,    5];
        wght  = [1,       4,        0.5,     4,        1.7,      1,   0.35];

        # Box-SLISE as per paper
    elseif (name == "box")
        intgl = [   0.95,    0.995,    1.005,     1.05,     1.1,   1.3, 1.8,   3];
        wght  = [      1,        4,        2,         4,    0.6,     1, 0.3, 0.1];
    end

    return createConfig(intgl,wght)
end

function createConfig( intgl::Array{R,1}, wght::Array{R,1} )
    intgl_h = intgl;
    wght_h = wght;

    intgl = vcat(-reverse(intgl), intgl);
    wght = vcat(reverse(wght), wght[2:end]);
    intervals = [intgl[1:end-1] intgl[2:end] wght]';

    intgl_h = vcat(intgl_h[ intgl_h .< 1.0], 1);
    wght_h = wght_h[ 1:length( intgl_h ) ];
    intgl_h = vcat(-reverse(intgl_h), intgl_h);
    wght_h = vcat(reverse(wght_h), wght_h[2:end]);
    intervals_h = [intgl_h[1:end-1] intgl_h[2:end] wght_h]';

    return Config(
        intervals,
        intervals_h
    );

end

function updateConfig!( c::SLiSeFilters.Config, intgl::Array{R,1}, wght::Array{R,1} )
    intgl_h = intgl;
    wght_h = wght;

    intgl = vcat(-reverse(intgl), intgl);
    wght = vcat(reverse(wght), wght[2:end]);
    c.Weights = [intgl[1:end-1] intgl[2:end] wght]';

    intgl_h = vcat(intgl_h[ intgl_h .< 1.0], 1);
    wght_h = wght_h[ 1:length( intgl_h ) ];
    intgl_h = vcat(-reverse(intgl_h), intgl_h);
    wght_h = vcat(reverse(wght_h), wght_h[2:end]);
    c.Weights_h = [intgl_h[1:end-1] intgl_h[2:end] wght_h]';

end

function optimise( methodName::String, z, alpha, c, lb = nothing )
    # Cast string to function
    f = getfield(SLiSeFilters, Symbol(methodName))

    if (lb == nothing)
        return f( z,alpha,c )
    else
        try
            return f( z,alpha,c,lb )
        catch e
            println(e)
            return lbfgsb( z,alpha,c,lb )
        end
    end
end

include("./definiteIntegrals.jl")
include("./gradMatrices.jl")
include("./residual.jl")
include("./points.jl")

include("./descent.jl")
include("./levMar.jl")
include("./lbfgsb.jl")
include("./optim.jl")

include("./wcr.jl")
include("./wise.jl")
include("./io.jl")

end # module


text:
using JuMP, Gurobi, OneClassSampling, SVDD, OneClassActiveLearning
include(joinpath(@__DIR__, "config.jl"))

experiment_name="dami_baseline_rand"

solver = with_optimizer(Gurobi.Optimizer; OutputFlag=0, Threads=1)
quality_metrics = Dict(:mcc => matthews_corr, :kappa => cohens_kappa, :f1 => f1_score)

threshold_strategies = [
    OutlierPercentageThresholdStrategy()
]

model_init_strategy = Dict{Symbol, Any}(
                            :svdd => Dict{Symbol, Any}(
                                :ROT_CTax => SVDD.SimpleCombinedStrategy(SVDD.RuleOfThumbScott(), SVDD.BoundedTaxErrorEstimate(0.05, 0.02, 0.98)),
                            ),
                            :sample_qe => Dict{Symbol, Any}(
                                :threshold_strategies => threshold_strategies
                        ))

NUM_RERUNS = 5
KDE_GAMMA = :scott

sampling_strategies = Vector{Dict{Symbol, Any}}()

for r in vcat(0.01:0.01:0.04, 0.05:0.05:1)
    for seed in 1:NUM_RERUNS
        push!(sampling_strategies, Dict{Symbol, Any}(
            :method => RandomRatioSampler(r),
            :seed => seed
        ))
    end
end

max_stars_repo_name: englhardt/ocs-evaluation

max_issues_repo_name: englhardt/ocs-evaluation

max_forks_repo_name: englhardt/ocs-evaluation


text:
using WaterLily
using LinearAlgebra: norm2
include("TwoD_plots.jl")

function circle(radius=8;Re=250,n=10,m=6)
    center, ν = radius*m/2, radius/Re
    body = AutoBody((x,t)->norm2(x .- center) - radius)
    Simulation((n*radius+2,m*radius+2), [1.,0.], radius; ν, body)
end

# sim_gif!(circle(20);duration=10,step=0.25)


text:
using WindMixing: calculate_training_scalings

@testset "Training Scaling" begin
    losses = (u=rand(Float32), v=rand(Float32), T=rand(Float32), ∂u∂z=rand(Float32), ∂v∂z=rand(Float32), ∂T∂z=rand(Float32))
    fractions = (T=rand(Float32), ∂T∂z=rand(Float32), profile=rand(Float32))

    scalings = calculate_training_scalings(losses, fractions)

    u = scalings.u * losses.u
    v = scalings.v * losses.v
    T = scalings.T * losses.T

    ∂u∂z = scalings.∂u∂z * losses.∂u∂z
    ∂v∂z = scalings.∂v∂z * losses.∂v∂z
    ∂T∂z = scalings.∂T∂z * losses.∂T∂z

    @test T / (u + v) ≈ fractions.T / (1 - fractions.T)
    @test ∂T∂z / (∂u∂z + ∂v∂z) ≈ fractions.∂T∂z / (1 - fractions.∂T∂z)
    @test (u + v + T) / (∂u∂z + ∂v∂z + ∂T∂z) ≈ fractions.profile / (1 - fractions.profile)
end


text:
# # Concrete casting heat of hydration
#
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/examples/Heat3D.ipynb)
#

using ReachabilityAnalysis, StructuralDynamicsODESolvers, MAT
using LinearAlgebra, SparseArrays
using LazySets.Arrays
using FastExpm

(@isdefined TARGET_FOLDER) ? nothing : TARGET_FOLDER = ""

# load internal (non-exported) functions
using ReachabilityAnalysis: ReachSolution, discretize, normalize, solve, center, step_size

include("../krylov_concept.jl")

# reduce default tolerance for vertex enumeration
LazySets.set_ztol(Float64, 1e-12)

using Plots
using Plots.PlotMeasures
using LaTeXStrings

using LazySets: center

function load_heat3d_matrices(; path=joinpath(@__DIR__, "Heat3D.mat"))
    vars = matread(path)

    K = vars["K"]
    C = vars["C"]
    timeIncr = vars["timeIncr"]
    fext = vars["fext"]
    qextTamb = vars["qextTamb"][:] # vector
    QhG =  vars["QhG"][:] # vector
    udots=vars["udots"]
    xs=vars["xs"]

    return K, C, qextTamb, QhG, timeIncr
end

function load_heat3d_params(; path=joinpath(@__DIR__, "params.mat"))
    vars       = matread(path)
    rho        = vars["rho"]
    m          = vars["m"]
    QFH        = vars["QFH"]
    Tambmin    = vars["Tambmin"]
    Tambvar    = vars["Tambvar"]
    Tambvarmin = vars["Tambvarmin"]
    Tambvarmax = vars["Tambvarmax"]
    Tini       = vars["Tini"]
    QFHmin     = vars["QFHmin"]
    QFHmax     = vars["QFHmax"]

    return rho, QFHmin, QFHmax, m, Tini, Tambmin, Tambvarmin, Tambvarmax, Tambvar, QFH
end

function load_heat3d_octaveSols(; path=joinpath(@__DIR__, "solsOctave.mat"))
    vars = matread(path)
    solAOctave = vars["Ts3DA"]
    solBOctave = vars["Ts3DB"]
    return solAOctave, solBOctave
end

# compute the maximum temperature at the given nodes
# by default we compute the max temperature in all the nodes,
# but you can specify a vector of indices to search;
# for example pass nodes=1:100 to make the search only on the first 100 nodes
function maxtemp(sol::ReachSolution; nodes=nothing)
    ## initialize search
    Tmax = -Inf; kmax = -1; imax = -1

    if isnothing(nodes)
        nodes = 1:dim(sol)
    end

    ## for each reach-set we compute the max bounds
    ## in every direction
    for (k, R) in enumerate(sol)
        X = set(R)
        q = center(X) + radius_hyperrectangle(X)
        θi = argmax(view(q, nodes))
        θ = q[θi]
        if θ > Tmax
            Tmax = θ
            kmax = k
            imax = θi
        end
    end

    δ = step_size(sol.alg)
    dtmax = tspan(sol[kmax])
    return (Tmax=Tmax, node=imax, step=kmax, dt=dtmax)
end

# compute the minimum temperature on the given nodes
function mintemp(sol::ReachSolution; nodes=nothing)
    ## initialize search
    Tmin = Inf; kmin = -1; imin = -1

    if isnothing(nodes)
        nodes = 1:dim(sol)
    end

    ## for each reach-set we compute the max bounds
    ## in every direction
    for (k, R) in enumerate(sol)
        X = set(R)
        q = center(X) - radius_hyperrectangle(X)
        θi = argmax(view(q, nodes))
        θ = q[θi]
        if θ < Tmin
            Tmin = θ
            kmin = k
            imin = θi
        end
    end

    δ = step_size(sol.alg)
    dtmin = tspan(sol[kmin])
    return (Tmin=Tmin, node=imin, step=kmin, dt=dtmin)
end

## Control nodes
controlNodeA = 11^2*6+1
controlNodeB = 11^2*9+1
nNodes = 11^3

function AhomCalc(C, K, qextTamb, QhG, rho, m)

    invC = inv( Matrix(C) );
    A = - invC * K |> Matrix
    n = size(A, 1)

    b0   =  invC * qextTamb ;
    bext =  invC * QhG * rho * m ;
    bsin =  invC * qextTamb ;

    omega = pi/12.;

    Ahom = [ A                    b0  bext  bsin      zeros(n) ;
             transpose( zeros(n))  0   0    0         0 ;
             transpose( zeros(n))  0  -m    0         0 ;
             transpose( zeros(n))  0   0    0         1 ;
             transpose( zeros(n))  0   0    -omega^2  0 ] ;

    return sparse(Ahom) # cast to sparse matrix
end

function solve_heat3d_rfem(; model=StepIntersect(Forward(inv=false)), X0=nothing, case=1)

    ## load vectors and params
    K, C, qextTamb, QhG, timeIncr = load_heat3d_matrices()
    rho, QFHmin, QFHmax, m, Tini, Tambmin, Tambvarmin, Tambvarmax, Tambvar, QFH = load_heat3d_params()
    usA, usB = load_heat3d_octaveSols()

    δ = timeIncr
    n = size(K, 1)

    NSTEPS = size(usA,2) - 1

    ## compute "homogeneized" system matrix
    Ahom = AhomCalc(C, K, qextTamb, QhG, rho, m)

    QFH     = (QFHmax+QFHmin)*0.5
    Tambvar = (Tambvarmax+Tambvarmin)*0.5

    if case == 1
        ## opcion singleton
        X0 = Singleton( Tini*ones(n)) ×
            Singleton((Tambmin+Tambvar*0.5)*ones( 1 ) )× Singleton( QFH*ones(1) ) ×
            Singleton((-Tambvar*0.5)*ones( 1 ) )× Singleton( zeros(1) )
        X0 = Singleton(LazySets.center(box_approximation(X0)))

    else
        ## opcion distribuida
        X0 = Singleton( Tini*ones(n)) ×
             Interval(Tambvarmin * 0.5 + Tambmin, Tambvarmax * 0.5 + Tambmin ) ×
             Interval( QFHmin, QFHmax ) ×
             Interval(-Tambvarmax * 0.5, -Tambvarmin*0.5) ×
             Singleton( zeros(1) )
        X0 = box_approximation(X0)
    end

    # set propagation method
    method = :sf
    println("Solving with set prop method (Case $case) $method")
    if method == :box
        prob = @ivp(x' = Ahom*x, x(0) ∈ X0)
        alg = BOX(δ=δ, approx_model=model)
        @time sol = solve(prob, alg=alg, NSTEPS=NSTEPS)

    elseif method == :sf
        @time sol = _solve_heat3d_sf(Ahom, X0, δ, NSTEPS, nodeidx=[controlNodeA, controlNodeB])
    end

    if case == 1
        X0 = Singleton(vcat(Tini*ones(n),
                      (Tambmin+Tambvar*0.5)*ones(1) ,
                      QFHmax*ones(1) ,
                      (-Tambvar*0.5)*ones(1),
                      zeros(1)))

        probOrbitMax = @ivp(x' = Ahom*x, x(0) ∈ X0)

        println("Solving with set prop method (Case 1) ORBIT")
        @time solOrbit = solve( probOrbitMax, NSTEPS=NSTEPS, alg=ORBIT(δ=δ))
        return sol, solOrbit

    elseif case == 2

        ## Orbit min
        X0 = Singleton(
              vcat( Tini*ones(n),
                    (Tambmin+Tambvarmin*0.5)*ones(1) ,
                    QFHmin*ones(1) ,
                    (-Tambvarmin*0.5)*ones(1),
                    zeros(1) ) ) ;

        probOrbit   = @ivp(x' = Ahom*x, x(0) ∈ X0)
        println("Solving with set prop method (Case 2) ORBIT min")
        @time solOrbitmin = solve( probOrbit, NSTEPS=NSTEPS, alg=ORBIT(δ=δ));

        ## Orbit min
        X0 = Singleton(
              vcat( Tini*ones(n),
                    (Tambmin+Tambvarmax*0.5)*ones(1) ,
                    QFHmax*ones(1) ,
                    (-Tambvarmax*0.5)*ones(1),
                    zeros(1) ) ) ;

        probOrbit = @ivp(x' = Ahom*x, x(0) ∈ X0)
        println("Solving with set prop method (Case 2) ORBIT max")
        @time solOrbitmax = solve( probOrbit, NSTEPS=NSTEPS, alg=ORBIT(δ=δ))
        return sol, solOrbitmin, solOrbitmax
   end

end

maxtemp_heat3d(sol; nodes=1:dim(sol), kwargs...) = maxtemp(sol, nodes=nodes)
maxtemp_heat3d(; nodes=nothing, kwargs...) = maxtemp(solve_heat3d(kwargs...), nodes=nodes)
mintemp_heat3d(sol; nodes=1:dim(sol), kwargs...) = mintemp(sol, nodes=nodes)
mintemp_heat3d(; nodes=nothing, kwargs...) = mintemp(solve_heat3d(kwargs...), nodes=nodes)

function solve_heat3d_implicit_euler()

    ## load vectors and params
    K, C, qextTamb, QhG, timeIncr = load_heat3d_matrices()
    rho, QFHmin, QFHmax, m, Tini, Tambmin, Tambvarmin, Tambvarmax, Tambvar, QFH = load_heat3d_params()
    usA, usB = load_heat3d_octaveSols()

    δ = timeIncr;  n = size(K, 1)
    NSTEPS = ( size(usA,2)-1 ) ;

    omega = pi/12.;

    M    = zeros(n, n)
    Tambvar = (Tambvarmin + Tambvarmax ) * 0.5 ;
    QFH = (QFHmin+ QFHmax)*0.5;

    R    = [ QhG*QFH*rho*m*exp(-m*δ*(i-1))+qextTamb*(Tambmin+Tambvar*(0.5 + 0.5*sin( omega*δ*(i-1)-pi/2.0 ))) for i in 1:NSTEPS+1 ] ;

    B   = Diagonal(ones(n));

    sys  =  SecondOrderConstrainedLinearControlContinuousSystem(M, Matrix(C), Matrix(K), B, nothing, R)

    U₀ = Tini*ones(n) ;

    prob = InitialValueProblem(sys, (U₀, U₀))
    alg  = BackwardEuler(Δt=δ)

    println("Solving with Backward Euler (Case 1)")
    @time solBackwardEulerSDCaseA = solve(prob, alg, NSTEPS=NSTEPS)

    R    = [ QhG*QFHmax*rho*m*exp(-m*δ*(i-1))+qextTamb*(Tambmin+Tambvarmax*(0.5 + 0.5*sin( omega*δ*(i-1)-pi/2.0 ))) for i in 1:NSTEPS+1 ] ;

    sys  =  SecondOrderConstrainedLinearControlContinuousSystem(M, Matrix(C), Matrix(K), B, nothing, R)
    prob = InitialValueProblem(sys, (U₀, U₀))

    println("Solving with Backward Euler (Case 2)")
    @time solBackwardEulerSDCaseB = solve(prob, alg, NSTEPS=NSTEPS)

    return solBackwardEulerSDCaseA, solBackwardEulerSDCaseB
end

function _solve_heat3d_sf(A, X0, δ, NSTEPS; nodeidx, mk=70, tol=1e-6)
    num_nodes = isa(nodeidx, Int) ? 1 : length(nodeidx)

    Eplus = _build_Eplus_optim(A, X0, δ; m=mk, tol=tol, cutoff=eps(Float64))

    dirs = Vector{LazySets.Arrays.SingleEntryVector{Float64}}()
    k = size(A, 1)
    for j in 1:num_nodes
        node = nodeidx[j]
        push!(dirs, SingleEntryVector(node, k, 1.0))
        push!(dirs, SingleEntryVector(node, k, -1.0))
    end
    cdirs = CustomDirections(dirs)

    propagacion = :expm_fast
    Ω0, Φᵀ = _discretize(A, δ, propagacion, X0, Eplus, mk, tol)

    Mout = Matrix{Float64}(undef, 2*num_nodes, NSTEPS)

    Aᵀδ = copy(transpose(A .* δ))

    for j in 1:num_nodes
        node = nodeidx[j]
        d = zeros(k)
        d[node] = 1.0
        println("Propagando indice $node ...")
        idx = (2(j-1) + 1):2j

        if isnothing(Φᵀ)
            _reach_homog_krylov_LGG09_modif!(view(Mout, idx, :), Ω0, Aᵀδ, d, NSTEPS; m=mk, tol=tol)
        else
            _reach_homog_dir_LGG09_modif!(view(Mout, idx, :), Ω0, Φᵀ, d, NSTEPS)
        end
    end
    sol = [TemplateReachSet(cdirs, view(Mout, :, k), (0.0 .. δ) + (k-1) * δ) for k in 1:NSTEPS] |> Flowpipe
end

function _discretize(A, δ, propagacion, X0::AbstractSingleton, Eplus, mk, tol)
    ΦX0 = _expv(A, X0.element, 1, δ, m=mk, tol=tol) |> Singleton
    Ω0 = CH(X0, minkowski_sum(ΦX0, Eplus))
    return Ω0, nothing
end

function _discretize(A, δ, propagacion, X0::Hyperrectangle, Eplus, mk, tol)
    if propagacion == :expm_base
        Φ = exp(Matrix(transpose(A) .* δ))
        Φ = sparse(Φ)
        droptol!(Φ, eps(Float64))
    elseif propagacion == :expm_fast
        Φ = real.(fastExpm(A .* δ, threshold=eps(Float64), nonzero_tol=eps(Float64)) )
    end
    Ω0 = CH(X0, Φ*X0 ⊕ Eplus)
    return Ω0, transpose(Φ)
end

function _reach_homog_dir_LGG09_modif!(Mout, Ω₀, Φᵀ, d, NSTEPS)
    rᵢ = copy(d)
    rᵢ₊₁ = similar(rᵢ)

    @inbounds for i in 1:NSTEPS
        Mout[1, i] = ρ(rᵢ, Ω₀)
        Mout[2, i] = ρ(-rᵢ, Ω₀)

        # update cache for the next iteration
        mul!(rᵢ₊₁, Φᵀ, rᵢ)
        copy!(rᵢ, rᵢ₊₁)
    end
    return nothing
end

#-----------------------------------------------------------
# Solutions LSD
solSDCase1, solSDCase2 = solve_heat3d_implicit_euler()

# Numerical solutions
tdom                     = times(solSDCase1)
solnodeSDCase1(node)     = displacements(solSDCase1, node)

tdom                     = times(solSDCase2)
solnodeSDCase2(node)     = displacements(solSDCase2, node)

#-----------------------------------------------------------
# Solutions Octave
solAOctave, solBOctave = load_heat3d_octaveSols()

#-----------------------------------------------------------
# Solutions SetBased

## Case 1
solRFEMCase1, solOrbitCase1 = solve_heat3d_rfem(case=1)

## Case 2
solRFEMCase2, solOrbitCase2min, solOrbitCase2max = solve_heat3d_rfem(case=2);

# Check maximum node temp
# @time maxtemp_heat3d(solRFEMCase2, nodes=1:1331)

#-----------------------------------------------------------
# Plots Case 1

timeMark = 200

indTimeMark = 155

tdom[indTimeMark]
solnodeSDCase1(controlNodeA)[indTimeMark]

fig = plot(tdom, solnodeSDCase1(controlNodeA),
           xlab=L"\textrm{Time (h)}",
           lc=:red, ylab=L"\textrm{Temperature } (^\circ\!C)",
           lab="",
           linestyle=:solid,
           legendfontsize=15, legend=:bottomright,
           tickfont=font(18, "Times"), guidefontsize=18,
           xtick=([0, 50, 100, 150, 200], [L"0", L"50", L"100", L"150", L"200"]),
           ytick=([20, 30, 40, 50, 60, 70, 80], [L"20", L"30", L"40", L"50", L"60", L"70", L"80"]),
           xguidefont=font(18, "Times"),
           yguidefont=font(18, "Times"),
           bottom_margin=6mm, left_margin=8mm, right_margin=4mm, top_margin=3mm,
           size=(900, 600))

plot!( tdom, solnodeSDCase1(controlNodeB), lc=:blue, lab="" )

# markers
plot!( [tdom[indTimeMark]], [solnodeSDCase1(controlNodeA)[indTimeMark]], marker=:utriangle, lc=:red, c=:red, lab=L"\textrm{Backward Euler at A}" )
plot!( [tdom[indTimeMark]], [solnodeSDCase1(controlNodeB)[indTimeMark]], marker=:circle, lc=:blue, c=:blue, lab=L"\textrm{Backward Euler at B}" )

plot!( fig, solRFEMCase1, vars=(0, controlNodeA), alpha=.2, lc=:magenta, c=:magenta, lab=L"\textrm{Flowpipe at A}")
plot!( fig, solRFEMCase1, vars=(0, controlNodeB), alpha=.2, lc=:green, c=:green, lab=L"\textrm{Flowpipe at B}")
#plot!(fig, tdom, solAOctave[controlNodeA,1:end], lc=:blue, lab="Backward Euler point A")

# subida
lens!(fig, [9, 12], [39, 43], inset = (1, bbox(0.15, 0.69, 0.28, 0.24)),
           subplot=5,
           tickfont=font(18, "Times"),
           xticks=([9.0, 10.5, 12.0], [ L"9",L"10.5", L"12"]),
           yticks=([39, 41, 43], [ L"39", L"41", L"43"])
)

lens!(fig, [50, 52], [66.5, 67], inset = (1, bbox(0.2, 0.38, 0.18, 0.18)),
           subplot=6,
           tickfont=font(18, "Times"),
           xticks=([50, 51, 52], [L"50", L"51", L"52"]),
           yticks=([66.5, 66.75, 67], [ L"66.5", "", L"67"])
)

lens!(fig, [50, 52], [83.3, 84.5], inset = (1, bbox(0.7, 0.1, 0.2, 0.2)),
           subplot=7,
           tickfont=font(18, "Times"),
           xticks=([50, 51, 52], [L"50", L"51", L"52"]),
           yticks=([83.5, 84.0, 84.5], [ L"83.5", "", L"84.5"])
)

savefig(fig, joinpath(TARGET_FOLDER, "heat3d_Case1.pdf"))

#-----------------------------------------------------------
# Plots Case 2

fig = plot(legendfontsize=16, legend=:topright,
           tickfont=font(18, "Times"), guidefontsize=18,
           xtick=([0, 50, 100, 150, 200], [L"0", L"50", L"100", L"150", L"200"]),
           ytick=([20, 30, 40, 50, 60, 70, 80], [L"20", L"30", L"40", L"50", L"60", L"70", L"80"]),
           xguidefont=font(18, "Times"),
           yguidefont=font(18, "Times"),
           bottom_margin=6mm, left_margin=8mm, right_margin=4mm, top_margin=3mm, size=(900, 600))

plot!(fig, solRFEMCase2, vars=(0, controlNodeA), alpha=.2, lc=:red, c=:red,
    xlab= L"\textrm{Time (h)}",
    ylab = L"\textrm{Temperature } (^\circ\!C)",
    lab=L"\textrm{Flowpipe at A}")#

#plot!(tdom, solnodeSDCase2(controlNodeA), xlab="Time", lc=:red, ylab="Temperature", lab="Backward Euler")
#plot!(fig, tdom, solBOctave[controlNodeA,1:end], lc=:blue, lab="Backward Euler point A")
#plot!( fig, solOrbitCase2max, vars=(0, controlNodeA), alpha=.2, lc=:red, lab="OrbitMax", setiestype=:path)
#plot!( fig, solOrbitCase2min, vars=(0, controlNodeA), alpha=.2, lc=:red, lab="OrbitMin", setiestype=:path)
plot!(fig,  solRFEMCase2, vars=(0, controlNodeB), alpha=.2, lc=:blue, c=:blue, lab=L"\textrm{Flowpipe at B}")#
plot!(fig, tdom, solnodeSDCase1(controlNodeA), lc=:red, lab=L"\textrm{Backward Euler at A}" )
plot!(fig, tdom, solnodeSDCase1(controlNodeB), lc=:blue, lab=L"\textrm{Backward Euler at B}")

savefig(fig, joinpath(TARGET_FOLDER, "heat3d_Case2.pdf"))

#=

julia>     Ahom = AhomCalc(C, K, qextTamb, QhG, rho, m)
1335×1335 Matrix{Float64}:
 -7.28091  0.77986  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.55972  …  0.0  0.0  0.0  0.0  0.0  0.0  4.16147  0.0082213   4.16147    0.0
  ⋮                                ⋮                        ⋮             ⋱  ⋮                        ⋮

  0.0      0.0      0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0         0.0  0.0  0.0  0.0  0.0  0.0  0.0      0.0        -0.0685389  0.0

julia> sparse(Ahom)
1335×1335 SparseMatrixCSC{Float64, Int64} with 10787 stored entries:
⠻⣦⡀⠀⠀⢸
⠀⠈⠻⣦⡀⢸
⠀⠀⠀⠈⠻⣾

julia> 1335 * 1335
1782225

julia> Ahom * Ahom;

julia> Ahom * Ahom;

julia> @time Ahom * Ahom;
  0.100113 seconds (2 allocations: 13.597 MiB)

julia> @time Ahom * Ahom;
  0.109501 seconds (2 allocations: 13.597 MiB)

julia> Ahom_sp = sparse(Ahom);

julia> @time Ahom_sp * Ahom_sp;
  0.133171 seconds (126.53 k allocations: 8.742 MiB, 99.41% compilation time)

julia> @time Ahom_sp * Ahom_sp;
  0.000643 seconds (7 allocations: 1.460 MiB)

=#

max_stars_repo_name: JuliaReach/SetPropagation-FEM-Examples

max_issues_repo_name: JuliaReach/SetPropagation-FEM-Examples

max_forks_repo_name: JuliaReach/SetPropagation-FEM-Examples


text:
module Power

using ..LIKWID:
    LibLikwid, power_initialized, topo_initialized, _powerinfo, powerinfo, TurboBoost, PowerDomain, PowerInfo, init_topology
using Unitful

const POWER_DOMAIN_SUPPORT_STATUS = UInt64(1) << 0
const POWER_DOMAIN_SUPPORT_LIMIT = UInt64(1) << 1
const POWER_DOMAIN_SUPPORT_POLICY = UInt64(1) << 2
const POWER_DOMAIN_SUPPORT_PERF = UInt64(1) << 3
const POWER_DOMAIN_SUPPORT_INFO = UInt64(1) << 4

"""
Initialize power measurements for the given CPU.
Returns the RAPL status, i.e. `false` (no RAPL) or `true` (RAPL working).
"""
function init(cpuid::Integer)
    power_initialized[] && return true
    hasRAPL = LibLikwid.power_init(cpuid)
    if Bool(hasRAPL)
        _powerinfo[] = unsafe_load(LibLikwid.get_powerInfo())
        _build_jl_power()
        power_initialized[] = true
        return true
    end
    return false
end

function init()
    power_initialized[] && return true
    init_topology()
    return init(0)
end

"Finalize power measurements."
function finalize()
    LibLikwid.power_finalize()
    power_initialized[] = false
    _powerinfo[] = nothing
    return nothing
end

function _build_jl_power()
    pi = _powerinfo[]
    tb = TurboBoost(
        pi.turbo.numSteps, unsafe_wrap(Array, pi.turbo.steps, pi.turbo.numSteps)
    )
    pds = ntuple(
        i -> begin
            p = pi.domains[i]
            PowerDomain(
                p.type,
                p.type,
                p.supportFlags,
                p.energyUnit,
                p.tdp,
                p.minPower,
                p.maxPower,
                p.maxTimeWindow,
                !iszero(p.supportFlags & POWER_DOMAIN_SUPPORT_INFO),
                !iszero(p.supportFlags & POWER_DOMAIN_SUPPORT_STATUS),
                !iszero(p.supportFlags & POWER_DOMAIN_SUPPORT_PERF),
                !iszero(p.supportFlags & POWER_DOMAIN_SUPPORT_POLICY),
                !iszero(p.supportFlags & POWER_DOMAIN_SUPPORT_LIMIT),
            )
        end,
        5,
    )
    powerinfo[] = PowerInfo(
        pi.baseFrequency,
        pi.minFrequency,
        tb,
        pi.hasRAPL,
        pi.powerUnit,
        pi.timeUnit,
        pi.uncoreMinFreq,
        pi.uncoreMaxFreq,
        pi.perfBias,
        pds,
    )
    return nothing
end

"""
    get_power_info() -> LIKWID.PowerInfo
Get power / energy information.
"""
function get_power_info()
    if !topo_initialized[]
        init_topology() || error("Couldn't init topology.")
    end
    if !power_initialized[]
        init(0) || error("Couldn't init power.")
    end
    return powerinfo[]
end

"""
Return the start value for a cpu (`cpuid`) for the domain with `domainid`.
"""
function start_power(cpuid::Integer, domainid::Integer)
    pt = LibLikwid.PowerType(domainid)
    pd = LibLikwid.PowerData(pt, 0, 0)
    pd_ref = Ref(pd)
    LibLikwid.power_start(pd_ref, cpuid, pt)
    return Int(pd_ref[].before)
end

"""
Return the stop value for a cpu (`cpuid`) for the domain with `domainid`.
"""
function stop_power(cpuid::Integer, domainid::Integer)
    pt = LibLikwid.PowerType(domainid)
    pd = LibLikwid.PowerData(pt, 0, 0)
    pd_ref = Ref(pd)
    LibLikwid.power_stop(pd_ref, cpuid, pt)
    return Int(pd_ref[].after)
end

"""
    get_power(p_start::Integer, p_stop::Integer, domainid::Integer)
Calculate the μJ from the values retrieved by `start_power()`
and `stop_power()`.
"""
function get_power(p_start::Integer, p_stop::Integer, domainid::Integer)
    pt = LibLikwid.PowerType(domainid)
    pd = LibLikwid.PowerData(pt, p_start, p_stop)
    energy = LibLikwid.power_printEnergy(Ref(pd))
    return energy * u"μJ"
end

"""
    measure(f; cpuid::Integer=0, domainid::Integer)
Measure / calculate the energy for the given `cpuid` and `domainid`
over the execution of the function `f` using [`Power.start_power`](@ref),
[`Power.stop_power`](@ref), etc. under the hood. Automatically
initializes and finalizes the power module.

# Examples
```julia
julia> LIKWID.Power.measure(; cpuid=0, domainid=0) do
           sleep(1)
       end
15.13702392578125 μJ
```
"""
function measure(f; cpuid::Integer=0, domainid::Integer)
    init(cpuid) || error("Couldn't init LIKWIDs power module.")
    try
        p_start = start_power(cpuid, domainid)
        f()
        p_stop = stop_power(cpuid, domainid)

        return get_power(p_start, p_stop, domainid)
    finally
        finalize()
    end
end

end # module


text:
module cg

using BenchmarkTools
using LinearAlgebra

function conjgrad(A, b, x_init; max_iter = Int(ceil(length(b)*5)), residual_tol = 1e-20 )

    n_dim = length(b)

    x0  = copy(x_init)
    r0  = A * x0 - b
    p0  = -r0
    rtr = dot(r0, r0)
    iter = 0

    Ap = similar(x0)
    for outer iter = 1:max_iter
        if rtr <= residual_tol
            break
        end
        # Ap      = A * p0
        mul!(Ap, A, p0)
        alpha   = rtr / dot(p0, Ap)
        @inbounds @simd for ii = 1:n_dim
            x0[ii] += alpha * p0[ii]
            r0[ii] += alpha * Ap[ii]
        end
        # x1      = x0 + alpha * p0
        # r1      = r0 + alpha * Ap
        rtr1    = dot(r0, r0)
        beta    = rtr1 / rtr
        # p1      = - r1 + beta * p0

        @inbounds  @simd  for ii = 1:n_dim
            p0[ii] = p0[ii] * beta - r0[ii]
        end

        # p0  = p1
        rtr = rtr1
        # x0  = x1
        # r0  = r1
    end
    return (x0, iter)
end

function test()

    # n_dim = 3000
    # 0.193567 seconds (6 allocations: 68.711 MiB, 28.44% gc time)
    # 35.281426 seconds (10 allocations: 117.578 KiB)

    n_dim = 5
    A = rand(n_dim, n_dim)
    A = A' * A
    b = rand(n_dim)

    x0 = rand(n_dim)

    # b = @benchmark x_theo = $A \ $b
    # b = @benchmark (x_cg, iter) = conjgrad($A, $b, $x0)

    b1 = @benchmark $A \ $b
    b2 = @benchmark conjgrad($A, $b, $x0)

    return (b1, b2)

    # println("dx: $(x_theo - x_cg)")
    # println("norm dx: $(norm(x_theo - x_cg))")

end

end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
using Images, Colors, ImageView, TestImages

z = ones(10,50);
y = 8; x = 2;
z[y,x] = 0
guidict = imshow(z)
idx = annotate!(guidict, AnnotationText(x, y, "x", color=RGB(0,0,1), fontsize=3))
idx2 = annotate!(guidict, AnnotationPoint(x+10, y, shape='.', size=4, color=RGB(1,0,0)))
idx3 = annotate!(guidict, AnnotationPoint(x+20, y-6, shape='.', size=1, color=RGB(1,0,0), linecolor=RGB(0,0,0), scale=true))
idx4 = annotate!(guidict, AnnotationLine(x+10, y, x+20, y-6, linewidth=2, color=RGB(0,1,0)))
idx5 = annotate!(guidict, AnnotationBox(x+10, y, x+20, y-6, linewidth=2, color=RGB(0,0,1)))
idx6 = annotate!(guidict, AnnotationPoints([(x+25, y-3), (x+24, y-3), (x+23, y-3)],
                                           shape='o', size=3, color=RGB(1,0,1),
                                           linecolor=RGB(0,0,0), scale=true))

sleep(2)
delete!(guidict, idx)
sleep(0.5)

img = testimage("lighthouse")
guidict = imshow(img)
scalebar(guidict, 30; x = 0.1, y = 0.05)


text:
#
# Protein - TMAO (compare new and old implementations)
#

include("../../src/MDDF.jl")
using PDBTools

atoms = PDBTools.readPDB("./structure.pdb")

protein = PDBTools.select(atoms,"protein")
solute = MDDF.Selection(protein,nmols=1 )

water = PDBTools.select(atoms,"water")
solvent = MDDF.Selection(water,natomspermol=3)

options = MDDF.Options()
trajectory = MDDF.PDBTraj("./traj200.pdb",solute,solvent)

pdb = MDDF.mddf(trajectory,options)


text:
using CleanEOSTable
using Base.Test
using HDF5

mktempdir() do path
    println("Downloading...")
    download("https://stellarcollapse.org/EOS/LS220_234r_136t_50y_analmu_20091212_SVNr26.h5.bz2",
             "$path/LS220_234r_136t_50y_analmu_20091212_SVNr26.h5.bz2")
    println("Unpacking...")
    run(`bunzip2 $path/LS220_234r_136t_50y_analmu_20091212_SVNr26.h5.bz2`)

    tablename = joinpath(path, "$path/LS220_234r_136t_50y_analmu_20091212_SVNr26.h5")
    new_tablename = joinpath(path, "$path/LS220_234r_136t_50y_analmu_20091212_SVNr26_corrected.h5")
    cleanEOSTable(tablename, new_tablename)

    println("Comparing tables...")
    table = h5open(tablename, "r")
    table_logenergy = read(table, "logenergy")
    new_table = h5open(new_tablename, "r")
    new_table_logenergy = read(new_table, "logenergy")

    @test count(new_table_logenergy .!= table_logenergy) == 741
end


text:
using BenchmarkTools

const BENCHMARKS_DIR = "../benchmark/"
cp(joinpath(BENCHMARKS_DIR, "benchmark_utils.jl"), joinpath(@__DIR__, "benchmark_utils.jl"), force=true)

benchmark_filepath(benchmark_name, benchmarks_dir=BENCHMARKS_DIR) =
    joinpath(benchmarks_dir, "benchmark_" * benchmark_name * ".jl")

@testset "Performance benchmarks" begin
    @info "Testing performance benchmarks..."

    @testset "Performance benchmark scripts" begin
        @info "  Testing performance benchmark scripts..."

        @testset "Static ocean benchmark" begin
            @info "    Running static ocean benchmarks..."

            replace_strings = [
                ("Ns = [(16, 16, 16), (32, 32, 32), (64, 64, 64), (128, 128, 128), (256, 256, 256)]",
                 "Ns = [(1, 1, 1)]")
            ]

            @test run_script(replace_strings, "static_ocean", benchmark_filepath("static_ocean"))
        end

        @testset "Channel benchmark" begin
            @info "    Running channel benchmarks..."

            replace_strings = [
                ("Ns = [(32, 32, 32), (64, 64, 64), (128, 128, 128), (256, 256, 256)]",
                 "Ns = [(1, 1, 1)]")
            ]

            @test run_script(replace_strings, "channel", benchmark_filepath("channel"))
        end

        @testset "Turbulence closures benchmark" begin
            @info "    Running turbulence closures benchmark..."

            replace_strings = [
                ("Ns = [(32, 32, 32), (256, 256, 128)]",
                 "Ns = [(1, 1, 1)]")
            ]

            @test run_script(replace_strings, "turbulence_closures", benchmark_filepath("turbulence_closures"))
        end

        @testset "Tracers benchmark" begin
            @info "    Running tracers benchmark..."

            replace_strings = [
                ("(32, 32, 32)", "(1, 1, 1)"),
                ("(256, 256, 128)", "(1, 1, 1)"),
                ("test_cases = [(0, 0), (0, 1), (0, 2), (1, 0), (2, 0), (2, 3), (2, 5), (2, 10)]",
                 "test_cases = [(0, 0), (2, 0), (2, 3)]")
            ]

            @test run_script(replace_strings, "tracers", benchmark_filepath("tracers"))
        end
    end

    @testset "Selected performance benchmarks" begin
        @info "  Running selected performance benchmarks..."

        for arch in archs
            sizes = [(16, 16, 16), (32, 32, 32)]
            for sz in sizes
                grid = RegularCartesianGrid(size=sz, extent=(1, 1, 1))
                model = IncompressibleModel(architecture=arch, grid=grid)

                @info "Benchmarking [$arch, Float64, $sz]..."
                b = @benchmark time_step!($model, 1)
                display(b)
                println()

                @test model isa IncompressibleModel
            end
        end
    end
end

rm(joinpath(@__DIR__, "benchmark_utils.jl"))


text:
using SpikingNetworks

parseParameters"""
taupre=20
taupost=20
tmax=50"""

parseParameters"""
N=100
groups=makeGroup([1,2],[N,N])
@group groups 1 tspike[i]=(i-1)*tmax/(N-1)
@group groups 2 tspike[i]=(N-i)*tmax/(N-1)
@connect groups 1 2 post==pre Apre[post,pre]=0.01
@connect groups 1 2 post==pre Apost[post,pre]=-0.01*1.05"""

parseEquations"""
dapre/dt=-apre/taupre
dapost/dt=-apost/taupost
dw/dt=0"""

net=init"""
N=200
apre=zeros(N,N)
apost=zeros(N,N)
w=zeros(N,N)
Apre=Apre
Apost=Apost"""

spike"t≂tspike"

plasticity"""
apre[j,i]+=Apre[j,i]
w[j,i]+=apost[j,i]
apost[i,j]+=Apost[i,j]
w[i,j]+=apre[i,j]"""

net.dt=0.5
@time run!(net,60/net.dt)

using SpikingNetwork: groups,tspike
index=sub2ind(size(net.w),groups[2],groups[1])
using MatlabPlot
mplot(tspike[groups[2]]-tspike[groups[1]],net.w[index])


text:
module ContinuousWavelet
include("transform.jl")
include("synthesis.jl")
include("surrogate.jl")
include("coherence.jl")
include("detrend.jl")
export CauchyWavelet, MorletWavelet, ContinuousWaveletTransform, transform, surrogate, surrogate!,
       coherence, detrend!, detrend, synthesis
end  # module ContinuousWavelet


text:
function test_xF(grid, xF)
    return grid.xF == xF
end

function test_xC(grid, xC)
    return grid.xC == xC
end

function test_yF(grid::Grid2D, yF)
    return grid.yF == yF
end

function test_yC(grid::Grid2D, yC)
    return grid.yC == yC
end

function xdomain_length(grid, Lx)
    return grid.Lx == Lx
end

function ydomain_length(grid::Grid2D, Ly)
    return grid.Ly == Ly
end


text:
struct QuerySynthesisSVMOracle <: Oracle
    classifier::LIBSVM.SVM
end

function QuerySynthesisSVMOracle(init_strategy, data_file::String)
    data, labels = load_data(data_file)
    return QuerySynthesisSVMOracle(init_strategy, data, labels)
end

function QuerySynthesisSVMOracle(init_strategy, data, labels)
    if !isa(init_strategy.gamma_strategy, SVDD.FixedGammaStrategy)
        throw(ArgumentError("Invalid gamma strategy type: $(typeof(init_strategy.gamma_strategy))"))
    end
    if !isa(init_strategy.C_strategy, SVDD.FixedCStrategy)
        throw(ArgumentError("Invalid C strategy type: $(typeof(init_strategy.C_strategy))"))
    end
    oracle = LIBSVM.svmtrain(data, labels; gamma=MLKernels.getvalue(init_strategy.gamma_strategy.kernel.alpha), cost=float(init_strategy.C_strategy.C))
    return QuerySynthesisSVMOracle(oracle)
end

function QuerySynthesisSVMOracle(data_file::String; gamma_search_range_oracle=10.0.^range(-2, stop=2, length=20), C=1, metric=cohens_kappa)
    data, labels = load_data(data_file)
    return QuerySynthesisSVMOracle(data, labels, gamma_search_range_oracle=gamma_search_range_oracle, C=C, metric=metric)
end

function QuerySynthesisSVMOracle(data, labels; gamma_search_range_oracle=10.0.^range(-2, stop=2, length=20), C=1, metric=cohens_kappa)
    best_gamma = 0
    best_score = -Inf
    for gamma in gamma_search_range_oracle
        c = LIBSVM.svmtrain(data, labels; gamma=float(gamma), cost=float(C))
        prediction = LIBSVM.svmpredict(c, data)[1]
        cm = ConfusionMatrix(prediction, labels, pos_class=:inlier, neg_class=:outlier)
        score = metric(cm)
        if score > best_score
            best_gamma = gamma
        end
    end
    return QuerySynthesisSVMOracle(LIBSVM.svmtrain(data, labels, gamma=float(best_gamma), cost=float(C)))
end

function QuerySynthesisSVMOracle(data, labels, params::Dict{Symbol, Any})
    if :init_strategy in keys(params)
        return QuerySynthesisSVMOracle(params[:init_strategy], data, labels)
    else
        additional_params = Dict(k => params[k] for k in [:gamma_search_range_oracle, :C, :metric] if k in keys(params))
        return QuerySynthesisSVMOracle(data, labels; additional_params...)
    end
end

function ask_oracle(oracle::QuerySynthesisSVMOracle, query_objects::Array{T, 2})::Vector{Symbol} where T <: Real
    return LIBSVM.svmpredict(oracle.classifier, query_objects)[1]
end


text:
include("CREST.jl")
include("SDE.jl")
import .CREST
import .SDE
using DataFrames
using SQLite

r = SDE.regionID("The Forge")
#itms = getItemsFromGroup(marketGroupID("Drones"))
itms = SDE.getAllEveItems()
l = get(CREST.latestLog())
o = readtable(l[1])
sort!(o, cols=[:typeID], by=SDE.itemName)
(sells, buys) = groupby(o, :bid)
itms = o[:typeID].data
itms = SDE.getItemsFromGroup(SDE.marketGroupID("Drones"))
q = CREST.getMarketHistory(r, itms)
db = SQLite.DB("testdb.sqlite")
CREST.storeHistoryData(db, q)
item = 28205
df = DataFrame(dbresp)
resp = CREST.getLatestOrders(r, itms)
test = CREST.loadMarketHistory(db, item, now() - Day(5))
# avg = ( buy * bvol + sell * svol )/ Volume
# Volume*avg = buy*bvol+sell*svol
# Volume = bvol + svol
# 2*avg - sell
testr = []
for (n,x) in test
  buy = x[1]
  sell = x[2]
  avg = x[3]
  vol = x[4]
  A = [buy sell; 1 1]
  B = [vol * avg; vol]
  print(A, "\n")
  #print(rank(A), "\n")
  X = A \ B
end
test[1]["LowPrice"].values[1]

CREST.storeMarketData("testdb.sqlite", resp)


text:
module ChemicalIdentifiers

    const DATA_DB = Dict{Symbol,Any}()
    const DATA_INFO = Dict{Symbol,Any}()

    import UUIDs, Unicode, Downloads
    import Arrow
    import Scratch, Preferences
    export search_chemical

    download_cache = ""

    const PKG_UUID = parse(UUIDs.UUID,"fa4ea961-1416-484e-bda2-883ee1634ba5")

"""
    function search_chemical(query,cache=cache=ChemicalIdentifiers.SEARCH_CACHE)

Given a query, performs a search on a database of over 70000 common compounds from CalebBell/Chemicals, returning a Named Tuple with the identifiers of the substance in question.

## Examples
```
julia>using ChemicalIdentifiers
julia>res = search_chemical("water")
(pubchemid = 962, CAS = (7732, 18, 5), formula = "H2O", MW = 18.01528, smiles = "O", InChI = "H2O/h1H2", InChI_key = "XLYOFNOQVPJJNP-UHFFFAOYSA-N", iupac_nam e = "oxidane", common_name = "water")

#worst case scenario, not found on the present databases
julia> @btime search_chemical("dimethylpyruvic acid22",nothing)
  273.700 μs (264 allocations: 15.05 KiB)
missing

#common compound found in the short database
julia> @btime search_chemical("methane",nothing)
  7.075 μs (57 allocations: 2.97 KiB)
(pubchemid = 297, CAS = (74, 82, 8), formula = "CH4", MW = 16.04246, smiles = "C", InChI = "CH4/h1H4", InChI_key = "VNWKTOKETHGBQD-UHFFFAOYSA-N", iupac_name = "methane", common_name = "methane")
```
A query is usually a string and its type is detected automatically when possible. the supported query types are:

### PubChemID:
By using any ``<:Integer` (or a string containing an Integer)
```julia
julia> search_chemical(8003)
(pubchemid = 8003, CAS = (109, 66, 0), formula = "C5H12", MW = 72.14878, smiles =
"CCCCC", InChI = "C5H12/c1-3-5-4-2/h3-5H2,1-2H3", InChI_key = "OFBQJSOFQDEBGM-UHFFFAOYSA-N", iupac_name = "pentane", common_name = "pentane")
```
### CAS registry number:
By using a Tuple of integers or a string with the digits separated by `-` :

```
julia> search_chemical((67,56,1))
(pubchemid = 887, CAS = (67, 56, 1), formula = "CH4O", MW = 32.04186, smiles = "CO", InChI = "CH4O/c1-2/h2H,1H3", InChI_key = "OKKJLVBELUTLKV-UHFFFAOYSA-N", iupac_name = "methanol", common_name = "methanol")

 search_chemical((67,56,1),nothing) == search_chemical("67-56-1",nothing) #true
```
### SMILES:
By using a string starting with `SMILES=` :
```
julia> search_chemical("SMILES=N")
(pubchemid = 222, CAS = (7664, 41, 7), formula = "H3N", MW = 17.03052, smiles = "N", InChI = "H3N/h1H3", InChI_key = "QGZKDVFQNNGYKY-UHFFFAOYSA-N", iupac_name = "azane", common_name = "ammonia")
```
### InChI:
By using a string starting with `InChI=1/` or `InChI=1S/` :

```
julia> search_chemical("InChI=1/C2H4/c1-2/h1-2H2")
(pubchemid = 6325, CAS = (74, 85, 1), formula = "C2H4", MW = 28.05316, smiles = "C=C", InChI = "C2H4/c1-2/h1-2H2", InChI_key = "VGGSQFUCUMXWEO-UHFFFAOYSA-N", iupac_name = "ethene", common_name = "ethene")
```

### InChI key:
By using a string with the pattern `XXXXXXXXXXXXXX-YYYYYYYYFV-P`:
```
julia> search_chemical("IMROMDMJAWUWLK-UHFFFAOYSA-N")
(pubchemid = 11199, CAS = (9002, 89, 5),
formula = "C2H4O", MW = 44.05256, smiles
= "C=CO", InChI = "C2H4O/c1-2-3/h2-3H,1H2", InChI_key = "IMROMDMJAWUWLK-UHFFFAOYSA-N", iupac_name = "ethenol", common_name
= "ethenol")
```
Searches by CAS and PubChemID are a little bit faster thanks to being encoded as native numeric types, other properties are stores as strings.

The package stores each query in `ChemicalIdentifiers.SEARCH_CACHE` as a `Dict{String,Any}`, so subsequent queries on the same (or similar) strings, dont pay the cost of searching in the database.

If you don't want to store the query, you could use `search_chemical(query,nothing)`, or, if you want your own cache to be used, pass your own cache via `search_chemical(query,mycache)`.

"""
  function search_chemical end
  include("data_script.jl")
  include("search_types.jl")
  include("search.jl")

  function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    Base.precompile(search_id_impl,(Int,Symbol))
    Base.precompile(search_id_impl,(String,Symbol))
    Base.precompile(search_id_impl,(Tuple{Int32,Int16,Int16},Symbol))
    Base.precompile(search_chemical_id,(AnyQuery,))
  end

  function clear_download_cache!()
    Preferences.@set_preferences!("CLEAR_CACHE" => "true")
    @info("ChemicalIdentifiers.jl download cache has been marked for deletion; restart your Julia session for this change to take effect.")
  end

  function __init__()
    clear_cache = Preferences.@load_preference("CLEAR_CACHE","false")
    if clear_cache == "true"
      @info "deleting download cache..."
      Scratch.delete_scratch!(PKG_UUID,"databases")
      Preferences.@set_preferences!("CLEAR_CACHE" => "false")
    end
    global download_cache = Scratch.@get_scratch!("databases")

    url_short = "https://github.com/CalebBell/chemicals/raw/master/chemicals/Identifiers/chemical%20identifiers%20pubchem%20small.tsv"
    url_long = "https://github.com/CalebBell/chemicals/raw/master/chemicals/Identifiers/chemical%20identifiers%20pubchem%20large.tsv"
    load_data!(:short,url= url_short)
    load_data!(:long,url = url_long)
    load_db!(:short)
    load_db!(:long)
    return nothing
  end
end


text:
using Base.Threads
function blockranges(nblocks, total_len)
    rem = total_len % nblocks
    main_len = div(total_len, nblocks)

    starts=Int[1]
    ends=Int[]
    for ii in 1:nblocks
        len = main_len
        if rem>0
            len+=1
            rem-=1
        end
        push!(ends, starts[end]+len-1)
        push!(starts, ends[end] + 1)
    end
    @assert ends[end] == total_len
    starts[1:end-1], ends
end

function sumby_multi_rg(by::AbstractVector{T},  val::AbstractVector{S})::Dict{T,S} where {T,S <: Number}
    (starts, ends) = blockranges(nthreads(), length(by))
    res = Dict{T,S}[]
    @threads for (ss, ee) in collect(zip(starts, ends))
        # @inbounds byv = @view by[ss:ee]
        # @inbounds valv = @view val[ss:ee]
        @inbounds push!(res, sumby_radixgroup!(by[ss:ee], val[ss:ee]))
    end

    szero  = zero(S)
    @inbounds res_fnl = res[1]
    for res1 in res[2:end]
        for k in keys(res1)
            res_fnl[k] = get(res_fnl,k,szero) + res1[k]
        end
    end
    return res_fnl
end

function sumby_multi_rs(by::AbstractVector{T},  val::AbstractVector{S})::Dict{T,S} where {T,S <: Number}
    (starts, ends) = blockranges(nthreads(), length(by))
    res = Dict{T,S}[]
    @threads for (ss, ee) in collect(zip(starts, ends))
        # @inbounds byv = @view by[ss:ee]
        # @inbounds valv = @view val[ss:ee]
        @inbounds push!(res, sumby_radixsort!(by[ss:ee], val[ss:ee]))
    end

    szero  = zero(S)
    @inbounds res_fnl = res[1]
    for res1 in res[2:end]
        for k in keys(res1)
            res_fnl[k] = get(res_fnl,k,szero) + res1[k]
        end
    end
    return res_fnl
end

function sumby_multi_van(by::AbstractVector{T},  val::AbstractVector{S})::Dict{T,S} where {T,S <: Number}
    (starts, ends) = blockranges(nthreads(), length(by))
    res = Dict{T,S}[]
    @threads for (ss, ee) in collect(zip(starts, ends))
        # @inbounds byv = @view by[ss:ee]
        # @inbounds valv = @view val[ss:ee]
        @inbounds push!(res, sumby!(by[ss:ee], val[ss:ee]))
    end

    szero  = zero(S)
    @inbounds res_fnl = res[1]
    for res1 in res[2:end]
        for k in keys(res1)
            res_fnl[k] = get(res_fnl,k,szero) + res1[k]
        end
    end
    return res_fnl
end


text:
@deprecate generate(t::Template, pkg::AbstractString) t(pkg)
@deprecate generate(pkg::AbstractString, t::Template) t(pkg)
@deprecate interactive_template() Template(; interactive=true)
@deprecate generate_interactive(pkg::AbstractString) Template(; interactive=true)(pkg)
@deprecate GitHubPages(; kwargs...) Documenter{TravisCI}(; kwargs...)


text:
__precompile__()

module RemarkableFloats

export Markers, MarkableFloat,
       MarkableF64, MarkableF32, MarkableF16,
       Marker0, Marker1, Marker2, Marker3,
       ʊ₀, ʊ₁, ʊ₂, ʊ₃,
       markable, unmarked, marked,
       getmark, setmark, getmarkvalue, setmarkvalue,
       incrmark, decrmark, unmark,
       @incrmark!, @decrmark!, @setmark!, @unmark!,
       has0marks, has1mark, has2marks, has3marks,
       has0or1marks, has1or2marks, has2or3marks,
       has0or2marks, has1or3marks

#=
       @mark!, @unmark!, @addmark!, @submark!, setmark
=#

import Base: reinterpret, hash,
        promote_type, promote_rule, convert,
        string, show,
	isfinite, isinf, isnan, issubnormal,
	signbit, sign, abs, flipsign, copysign,
	(==), (!=), (<), (<=), (>=), (>), isequal, isless,
	(+), (-), (*), (/), (\), (^),
	div, fld, cld, rem, mod, divrem, fldmod,
	ceil, floor, trunc, float

import Base.Math: acos, acosd, acosh, acot, acotd, acoth, acsc, acscd, acsch,
	asec, asecd, asech, asin, asind, asinh, atan, atan2, atand, tanh, beta,
	cbrt, clamp, clamp!, cos, cosc, cosd, cosh, cospi, cot, cotd, coth, csc,
	cscd, csch, deg2rad, exp, exp10, exp2, expm1, exponent, frexp, gamma,
	hypot, lbeta, ldexp, lfact, lgamma, log, log10, log1p, log2, max, min,
	minmax, mod2pi, modf, rad2deg, rem2pi, sec, secd, sech, significand, sin,
	sinc, sincos, sind, sinh, sinpi, sqrt, tan, tand, tanh

# avoid importing Float64, Float32, Float16 (bad things can happen)
# use const qualified equivalences instead

const FLOAT64 = Base.Float64
const FLOAT32 = Base.Float32
const FLOAT16 = Base.Float16

include("type/markers.jl")
include("type/floatbasis.jl")
include("type/ishas.jl")
include("type/getset.jl")
include("type/promote.jl")
include("type/comparison.jl")
include("type/incrdecr.jl")
include("type/macros.jl")
include("type/stringshow.jl")

include("math/predicates.jl")
include("math/prearith.jl")
include("math/arith.jl")
include("math/mathfuncs.jl")

end # RemarkableFloats


text:
const _msise_so_filename =
            string(
                    @__DIR__, "/nrlmsise00",
                    isapple ? "_macos" : islinux ? "_linux" : "_win",
                    ".so",
                )

_check_nrlmsise00() = _check_so(_msise_so_filename)
_check_nrlmsise00()

# export gtd7
gtd7(IYD::Int32, SEC::Float32,
    ALT::Float32, GLAT::Float32, GLONG::Float32,
    STL::Float32, F107A::Float32, F107::Float32,
    AP::Array{Float32,1}, MASS::Int32, Tinf_scl::Float32,
    D::Array{Float32,1}, T::Array{Float32,1},
) = ccall((:gtd7_, _msise_so_filename),
    Nothing,
    (Ref{Int32},   # IYD
        Ref{Float32}, # SEC
        Ref{Float32}, # ALT
        Ref{Float32}, # GLAT
        Ref{Float32}, # GLONG
        Ref{Float32}, # STL
        Ref{Float32}, # F107A
        Ref{Float32}, # F107
        Ptr{Float32}, #  AP
        Ref{Int32},   #  MASS,
        Ref{Float32}, #  Tinf_scl
        Ptr{Float32}, #  D
        Ptr{Float32},),
    IYD, SEC,
    ALT, GLAT, GLONG,
    STL, F107A, F107,
    AP, MASS, Tinf_scl,
    D, T,
)

function test_gtd7()
    aaps_msise  = Array{Float32,1}(undef, 7)
    aaps_msise .= 1.0

    iyd = Int32(19199)
    sec = Float32(12 * 3600.0)
    alt  = Float32(30.0)
    glat = Float32(15.0)
    glon = Float32(36.0)
    hrl_32 = Float32(0.0)
    fbar = Float32(180.0)
    f10p7 = Float32(180.0)
    mmass = Int32(48)
    Tinf_scl = Float32(1.0)
    d    = Array{Float32,1}(undef, 9)
    temp = Array{Float32,1}(undef, 2)

    gtd7(iyd, sec, alt, glat, glon,
                hrl_32, fbar, f10p7, aaps_msise,
                mmass, Tinf_scl, d, temp,
            )
    # println("# ", join(d, ", "))
    # println("# ", join(temp, ", "))
    d, temp

end