text:
#=
The sum of the primes below 10 is 2 + 3 + 5 + 7 = 17.

Find the sum of all the primes below two million.
=#
using Primes
@time println(sum(primes(2_000_000)))

max_stars_repo_name: daniel-beard/JuliaProjectEuler

max_issues_repo_name: daniel-beard/JuliaProjectEuler

max_forks_repo_name: daniel-beard/JuliaProjectEuler


text:
using forecast
using Test

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

ext: jl


text:
using Test
using MultiScaleTreeGraph
using DataFrames
using MutableNamedTuples
using Graphs

@testset "read_mtg" begin
    include("test-read_mtg.jl")
end

@testset "simple node functions" begin
    include("test-nodes.jl")
end

@testset "descendants()" begin
    include("test-descendants.jl")
end

@testset "ancestors()" begin
    include("test-ancestors.jl")
end

@testset "mutation" begin
    include("test-mutation.jl")
    include("test-insert_node.jl")
    include("test-traverse.jl")
    include("test-transform.jl")
end

@testset "Deletion / Pruning" begin
    include("test-delete-prune.jl")
end

@testset "conversion" begin
    include("test-conversion.jl")
end


text:
# SECTION 3: DECREASE INDI-ACCOUNT'S INCOME
# --------------------------------
# NOTE: we use the expr, macros in SECTION 2

# 3.1 the parameter combinations of this section's reforms
ReformYear = 2020 - env.START_YEAR + 1
idx_plot = (2010:2110) .- env.START_YEAR .+ 1
reform_indiacct_expr = Dict(
   "Benchmark" => :(
      nothing
   ),  # benchmark, no change
   "NoFirmTrans2Indi" => :(
      Pt[:𝕒][ReformYear:end] .= 0.0
   ),  # firm's contrib now only go to pool & retirees's indi acct
   "NoFirmTrans2Retiree" => :(
      # Pt[:𝕒][ReformYear:end] .= 0.0;
      Pt[:𝕓][ReformYear:end] .= 0.0;
   ),  # canceling transfer payment to retirees
   "Everything2Pool" => :(
      # it is equvalent to the parameters let new piM = old piM, where new phi = 0
      tmpNewZeta = (1.0 + Pt[:z][ReformYear-1] * Pt[:η][ReformYear-1]) * (Pt[:ζ][ReformYear-1] + Pt[:ϕ][ReformYear-1]) ;
      tmpNewZeta /= 1.0 + Pt[:z][ReformYear-1] * Pt[:η][ReformYear-1] - Pt[:ϕ][ReformYear-1] ;
      # let new phi = 0, bba = 0, bbb = 0
      Pt[:ϕ][ReformYear:end] .= 0.0;
      Pt[:𝕒][ReformYear:end] .= 0.0;
      Pt[:𝕓][ReformYear:end] .= 0.0;
      # set new zeta to equilize the old piM & new piM
      Pt[:ζ][ReformYear:end] .= tmpNewZeta;
   ),  # hold personal contributions, both personal & firm's contrib all go to pool, i.e. no individual account
   "JustKillIndiAcct" => :(
      Pt[:ϕ][ReformYear:end] .= 0.0;
      Pt[:𝕒][ReformYear:end] .= 0.0;
      Pt[:𝕓][ReformYear:end] .= 0.0;
   ) # just kill individual accounts but change nothing (no contribution, no transfer payments, no benefits)
)

# 3.2 spaces to store results
res_IndiInDesc = Dict(
   "reformNames" => [ "Benchmark", "NoFirmTrans2Indi",       "NoFirmTrans2Retiree",     "Everything2Pool",                      "JustKillIndiAcct" ],
   "plotlegends" => [ "Benchmark", "No transfer to workers", "No transfer to retirees", "Keep contribution but all to pooling", "Canceling individual accounts" ],
   "plotstyles" => [ "-", "-.b", "-r", ":m", "--c" ],
   "gap/exp" => [],
   "gap/in" => [],
   "gap/gdp" => [],
   "gap/taxrev" => [],
   "Dt" => [],
)

# 3.3 evaluate reforms
for tmpreform in res_IndiInDesc["reformNames"]
   # ------ data loading
   include("$(pwd())\\src\\proc_VarsDeclare.jl"); # use pwd() to avoid possible ambiguity of working directory
   include("$(pwd())\\src\\proc_InitPars.jl");
   # ------ reform
   eval( reform_indiacct_expr[tmpreform] )
   # ------ evaluate model & save
   eval(expr_RunModel)
   push!( res_IndiInDesc["gap/exp"], 100 .* (Dt[:LI] ./ Dt[:AggPoolExp])[idx_plot] )
   push!( res_IndiInDesc["gap/in"] , 100 .* (Dt[:LI] ./ Dt[:AggPoolIn])[idx_plot] )
   push!( res_IndiInDesc["gap/gdp"], 100 .* (Dt[:LI] ./ Dt[:Y])[idx_plot] )
   push!( res_IndiInDesc["gap/taxrev"], 100 .* (Dt[:LI] ./ (Dt[:TRc] .+ Dt[:TRw]))[idx_plot] )
   push!( res_IndiInDesc["Dt"], copy(Dt) )
end # end loops

#


text:
import Base: copy, copyto!, ==
import Random: rand, seed!
import RandomNumbers: AbstractRNG, gen_seed, seed_type

"""
```julia
AbstractXorshift64 <: AbstractRNG{UInt64}
```

The base abstract type for `Xorshift64` and `Xorshift64Star`.
"""
abstract type AbstractXorshift64 <: AbstractRNG{UInt64} end

for star in (false, true)
    rng_name = Symbol(string("Xorshift64", star ? "Star" : ""))
    @eval begin
        mutable struct $rng_name <: AbstractXorshift64
            x::UInt64
        end

        function $rng_name(seed::Integer=gen_seed(UInt64))
            r = $rng_name(zero(UInt64))
            seed = init_seed(seed, UInt64)
            seed!(r, seed)
            r
        end

        @inline function xorshift_next(r::$rng_name)
            r.x ⊻= r.x << 18
            r.x ⊻= r.x >> 31
            r.x ⊻= r.x << 11
            $(star ? :(r.x * 2685821657736338717) : :(r.x))
        end
    end
end

@inline seed_type(::Type{T}) where T <: AbstractXorshift64 = UInt64

function copyto!(dest::T, src::T) where T <: AbstractXorshift64
    dest.x = src.x
    dest
end

copy(src::T) where T <: AbstractXorshift64 = copyto!(T(), src)

==(r1::T, r2::T) where T <: AbstractXorshift64 = r1.x == r2.x

function seed!(r::AbstractXorshift64, seed::Integer=gen_seed(UInt64))
    seed == 0 && error("0 cannot be the seed")
    r.x = seed % UInt64
    xorshift_next(r)
    r
end

@inline rand(r::AbstractXorshift64, ::Type{UInt64}) = xorshift_next(r)


text:
# Tested in Julia 1.5.1
using DelimitedFiles
data = readdlm("./data/adventInput08.txt")
function accLoop(data;acc= 0, indx = 1)
    usedIndx = [indx]
    outIndx = false
    for i in 1:size(data)[1]
        op, arg = data[indx,:]
        if op == "acc"
            acc += arg
            indx += 1
            outIndx = indx == size(data)[1]+1
            if outIndx
                break
            end
            indx ∉ usedIndx ? push!(usedIndx, indx) : break
        elseif op == "jmp"
            indx += arg
            outIndx = indx == size(data)[1]+1
            if outIndx
                break
            end
            indx ∉ usedIndx ? push!(usedIndx, indx) : break
        elseif op == "nop"
            indx += 1
            outIndx = indx == size(data)[1]+1
            if outIndx
                break
            end
            indx ∉ usedIndx ? push!(usedIndx, indx) : break
        end
    end
    acc, outIndx
end
function fixingBug(data; acc=0, outIndx = false)
    for (indx,(op, arg)) in enumerate(eachrow(data))
        if op == "jmp"
            data_tmp = copy(data)
            data_tmp[indx,1] = "nop"
            acc, outIndx = accLoop(data_tmp)
            if outIndx
                println(indx)
                break
            end
        elseif op == "nop"
            data_tmp = copy(data)
            data_tmp[indx,1] = "jmp"
            acc, outIndx = accLoop(data_tmp)
            if outIndx
                println(indx)
                break
            end
        end
    end
    acc
end
accLoop(data)
fixingBug(data)


text:
export recreatedata

function recreatedata(file; dir = [])

    if isempty(dir)
        l = load(file)
    else
        l = load(datadir(dir,file))
    end

    startidx = findfirst("/",l[:script])[1]
    endidx = findfirst("#",l[:script])[1]
    scriptname = l[:script][startidx+1:endidx-1]

    commit = l[:gitcommit]

    generate_deffile(; commit = commit, script = scriptname)
end


text:
# ciphertext block {{{1
const ciphertext = filter(isalpha, """
MOMUD EKAPV TQEFM OEVHP AJMII CDCTI FGYAG JSPXY ALUYM NSMYH
VUXJE LEPXJ FXGCM JHKDZ RYICU HYPUS PGIGM OIYHF WHTCQ KMLRD
ITLXZ LJFVQ GHOLW CUHLO MDSOE KTALU VYLNZ RFGBX PHVGA LWQIS
FGRPH JOOFW GUBYI LAPLA LCAFA AMKLG CETDW VOELJ IKGJB XPHVG
ALWQC SNWBU BYHCU HKOCE XJEYK BQKVY KIIEH GRLGH XEOLW AWFOJ
ILOVV RHPKD WIHKN ATUHN VRYAQ DIVHX FHRZV QWMWV LGSHN NLVZS
JLAKI FHXUF XJLXM TBLQV RXXHR FZXGV LRAJI EXPRV OSMNP KEPDT
LPRWM JAZPK LQUZA ALGZX GVLKL GJTUI ITDSU REZXJ ERXZS HMPST
MTEOE PAPJH SMFNB YVQUZ AALGA YDNMP AQOWT UHDBV TSMUE UIMVH
QGVRW AEFSP EMPVE PKXZY WLKJA GWALT VYYOB YIXOK IHPDS EVLEV
RVSGB JOGYW FHKBL GLXYA MVKIS KIEHY IMAPX UOISK PVAGN MZHPW
TTZPV XFCCD TUHJH WLAPF YULTB UXJLN SIJVV YOVDJ SOLXG TGRVO
SFRII CTMKO JFCQF KTINQ BWVHG TENLH HOGCS PSFPV GJOKM SIFPR
ZPAAS ATPTZ FTPPD PORRF TAXZP KALQA WMIUD BWNCT LEFKO ZQDLX
BUXJL ASIMR PNMBF ZCYLV WAPVF QRHZV ZGZEF KBYIO OFXYE VOWGB
BXVCB XBAWG LQKCM ICRRX MACUO IKHQU AJEGL OIJHH XPVZW JEWBA
FWAML ZZRXJ EKAHV FASMU LVVUT TGK
""")
# }}}

# character frequencies {{{1
const letters = Dict{Char, Float32}(
    'E' => 12.702,
    'T' => 9.056,
    'A' => 8.167,
    'O' => 7.507,
    'I' => 6.966,
    'N' => 6.749,
    'S' => 6.327,
    'H' => 6.094,
    'R' => 5.987,
    'D' => 4.253,
    'L' => 4.025,
    'C' => 2.782,
    'U' => 2.758,
    'M' => 2.406,
    'W' => 2.361,
    'F' => 2.228,
    'G' => 2.015,
    'Y' => 1.974,
    'P' => 1.929,
    'B' => 1.492,
    'V' => 0.978,
    'K' => 0.772,
    'J' => 0.153,
    'X' => 0.150,
    'Q' => 0.095,
    'Z' => 0.074)
const digraphs = Dict{AbstractString, Float32}(
    "TH" => 15.2,
    "HE" => 12.8,
    "IN" => 9.4,
    "ER" => 9.4,
    "AN" => 8.2,
    "RE" => 6.8,
    "ND" => 6.3,
    "AT" => 5.9,
    "ON" => 5.7,
    "NT" => 5.6,
    "HA" => 5.6,
    "ES" => 5.6,
    "ST" => 5.5,
    "EN" => 5.5,
    "ED" => 5.3,
    "TO" => 5.2,
    "IT" => 5.0,
    "OU" => 5.0,
    "EA" => 4.7,
    "HI" => 4.6,
    "IS" => 4.6,
    "OR" => 4.3,
    "TI" => 3.4,
    "AS" => 3.3,
    "TE" => 2.7,
    "ET" => 1.9,
    "NG" => 1.8,
    "OF" => 1.6,
    "AL" => 0.9,
    "DE" => 0.9,
    "SE" => 0.8,
    "LE" => 0.8,
    "SA" => 0.6,
    "SI" => 0.5,
    "AR" => 0.4,
    "VE" => 0.4,
    "RA" => 0.4,
    "LD" => 0.2,
    "UR" => 0.2)
const trigraphs = Dict{AbstractString, Float32}(
    "THE" => 18.1,
    "AND" => 7.3,
    "ING" => 7.2,
    "ION" => 4.2,
    "ENT" => 4.2,
    "HER" => 3.6,
    "FOR" => 3.4,
    "THA" => 3.3,
    "NTH" => 3.3,
    "INT" => 3.2,
    "TIO" => 3.1,
    "ERE" => 3.1,
    "TER" => 3.0,
    "EST" => 2.8,
    "ERS" => 2.8,
    "HAT" => 2.6,
    "ATI" => 2.6,
    "ATE" => 2.5,
    "ALL" => 2.5,
    "VER" => 2.4,
    "HIS" => 2.4,
    "HES" => 2.4,
    "ETH" => 2.4,
    "OFT" => 2.2,
    "STH" => 2.1,
    "RES" => 2.1,
    "OTH" => 2.1,
    "ITH" => 2.1,
    "FTH" => 2.1,
    "ONT" => 2.0)
# 1}}}

function decrypt(enc::ASCIIString, key::ASCIIString)
    const enclen = length(enc)
    const keylen = length(key)

    if keylen < enclen
        key = (key^(div(enclen - keylen, keylen) + 2))[1:enclen]
    end

    msg = Array(Char, enclen)

    for i=1:enclen
        msg[i] = Char((Int(enc[i]) - Int(key[i]) + 26) % 26 + 65)
    end

    msg::Array{Char, 1}
end

function cryptanalyze(enc::ASCIIString; maxkeylen::Integer = 20)
    const enclen = length(enc)
    maxkey = ""
    maxdec = ""
    maxscore = 0.0

    for keylen=1:maxkeylen
        key = Array(Char, keylen)
        idx = filter(x -> x % keylen == 0, 1:enclen) - keylen + 1

        for i=1:keylen
            maxsubscore = 0.0

            for j='A':'Z'
                subscore = 0.0

                for k in decrypt(enc[idx], ascii(string(j)))
                    subscore += get(letters, k, 0.0)
                end

                if subscore > maxsubscore
                    maxsubscore = subscore
                    key[i] = j
                end
            end

            idx += 1
        end

        key = join(key)
        const dec = decrypt(enc, key)
        score = 0.0

        for i in dec
            score += get(letters, i, 0.0)
        end

        for i=1:enclen - 2
            const digraph = string(dec[i], dec[i + 1])
            const trigraph = string(dec[i], dec[i + 1], dec[i + 2])

            if haskey(digraphs, digraph)
                score += 2 * get(digraphs, digraph, 0.0)
            end

            if haskey(trigraphs, trigraph)
                score += 3 * get(trigraphs, trigraph, 0.0)
            end
        end

        if score > maxscore
            maxscore = score
            maxkey = key
            maxdec = dec
        end
    end

    (maxkey, join(maxdec))::Tuple{ASCIIString, ASCIIString}
end

key, dec = cryptanalyze(ciphertext)
println("key: ", key, "\n\n", dec)

# post-compilation profiling run
gc()
t = @elapsed cryptanalyze(ciphertext)
println("\nelapsed time: ", t, " seconds")

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
@testset "Mutable interface" begin

    @testset "can_set_first" begin
        @test @inferred(can_set_first(LinRange)) == false
        @test @inferred(can_set_first(LinMRange)) == true
        @test @inferred(can_set_first(StepMRange)) == true
        @test @inferred(can_set_first(UnitMRange)) == true
        @test @inferred(can_set_first(StepMRangeLen)) == true
    end

    @testset "can_set_last" begin
        @test @inferred(can_set_last(LinRange)) == false
        @test @inferred(can_set_last(LinMRange)) == true
        @test @inferred(can_set_last(StepMRange)) == true
        @test @inferred(can_set_last(StepMRangeLen)) == true
        @test @inferred(can_set_last(UnitMRange)) == true
        @test @inferred(can_set_last(OneToMRange)) == true
    end

    @testset "can_set_step" begin
        @test @inferred(can_set_step(LinRange)) == false
        @test @inferred(can_set_step(LinMRange)) == false
        @test @inferred(can_set_step(StepMRange)) == true
        @test @inferred(can_set_step(StepMRangeLen)) == true
    end

    for (r1,b,v,r2) in ((UnitMRange(1,3), true, 2, UnitMRange(2,3)),
                        (StepMRange(1,1,4), true, 2, StepMRange(2,1,4)),
                        (StepSRange(1,1,4), false, nothing, nothing),
                        (LinMRange(1,3,3), true, 2, LinMRange(2,3,3)),
                        (StepMRangeLen(1,1,3), true, 2, StepMRangeLen(2,1,3)),
                        (StepSRangeLen(1,1,3), false, nothing, nothing),
                        ([1,2,3], true, 2, [2,2,3]))
        @testset "set_first-$(r1)" begin
            x = @inferred(can_set_first(r1))
            @test x == b
            if x
                set_first!(r1, v)
                @test r1 == r2
            end
        end
    end

    for (r1,b,v,r2) in ((UnitMRange(1,3), true, 5, UnitMRange(1,5)),
                        (OneToMRange(4), true, 5, OneToMRange(5)),
                        (OneToMRange(4), true, Int32(5), OneToMRange(5)),
                        (StepMRange(1,1,4), true, 5, StepMRange(1,1,5)),
                        (StepSRange(1,1,4), false, nothing, nothing),
                        (LinMRange(1,3,3), true, 4, LinMRange(1,4,3)),
                        (StepMRangeLen(1,1,3), true, 4, StepMRangeLen(1,1,4)),
                        ([1,2,3], true, 2, [1,2,2]))
        @testset "set_last-$(r1)" begin
            x = @inferred(can_set_last(r1))
            @test x == b
            if x
                @test @inferred(set_last!(r1, v)) == r2
            end
        end
    end

    @testset "set_step!" begin
        @testset "has_step" begin
            @test @inferred(has_step(OneToMRange{Int})) == true
            @test @inferred(has_step([])) == false
            @test @inferred(has_step(Vector)) == false
        end
        @test_throws ErrorException set_step!(OneToMRange(10), 2)
        @test_throws ErrorException set_step!(UnitMRange(1, 10), 2)

        for (r1,b,v,r2) in ((UnitMRange(1,3), false, nothing, nothing),
                            (StepMRange(1,1,4), true, 2, StepMRange(1,2,4)),
                            (StepMRange(1,1,4), true, UInt32(2), StepMRange(1,2,4)),
                            (StepMRangeLen(1,1,4), true, 2, StepMRangeLen(1,2,4)),
                           )

            @testset "set_step!-$(r1)" begin
                x = @inferred(can_set_step(r1))
                @test x == b
                if x
                    @test @inferred(set_step!(r1, v)) == r2
                end
            end
        end
    end

    @testset "set_ref!" begin
        @test @inferred(set_ref!(StepMRangeLen(1, 1, 10), 2)) == StepMRangeLen(2, 1, 10)
        @test @inferred(set_ref!(StepMRangeLen(1, 1, 10), UInt32(2))) == StepMRangeLen(2, 1, 10)
    end

    @testset "set_offset!" begin
        @test @inferred(set_offset!(StepMRangeLen(1, 1, 10), 2)) == StepMRangeLen(1, 1, 10, 2)
        @test @inferred(set_offset!(StepMRangeLen(1, 1, 10), UInt32(2))) == StepMRangeLen(1, 1, 10, 2)
    end

    @testset "is_static" begin
        @test @inferred(is_static(Any[])) == false
        @test @inferred(is_static(UnitSRange(1, 10))) == true
    end

    # as_[mutable/immutable/static]
    for (i,m,s) in ((OneTo(4), OneToMRange(4), OneToSRange(4)),
                    (UnitRange(1, 3), UnitMRange(1, 3), UnitSRange(1, 3)),
                    (StepRange(1, 1, 4), StepMRange(1, 1, 4), StepSRange(1, 1, 4)),
                    (StepRangeLen(1, 1, 4), StepMRangeLen(1, 1, 4), StepSRangeLen(1, 1, 4)),
                    (LinRange(1, 4, 4), LinMRange(1, 4, 4), LinSRange(1, 4, 4)),
                   )
        @testset "as_mutable($(typeof(i).name))" begin
            @test ismutable(as_mutable(i)) == true
            @test ismutable(as_mutable(m)) == true
            @test ismutable(as_mutable(s)) == true
        end

        @testset "as_immutable($(typeof(i).name))" begin
            @test isimmutable(as_immutable(i)) == true
            @test isimmutable(as_immutable(m)) == true
            @test isimmutable(as_immutable(s)) == true
        end

        @testset "as_static($(typeof(i).name))" begin
            @test is_static(as_static(i)) == true
            @test is_static(as_static(m)) == true
            @test is_static(as_static(s)) == true
        end
    end
end


text:
# add Dance    # pkg mode

using Dance

start_project("ToDo")


text:
using ReachabilityAnalysis
Dict([
    "opts" => Dict(:T=>20.0, :alg=>GLGM06(δ=0.01)),
    "X0"   => Hyperrectangle(low=[0.002, 0.0, 0.0, 0.0, 0.001, 0.0, 0.0, 0.0],
                            high=[0.0025, 0.0, 0.0, 0.0, 0.0015, 0.0, 0.0, 0.0]),
    "X1"   => Hyperrectangle([0.00225, 0.0, 0.0, 0.0, 0.00125, 0.0, 0.0, 0.0],
                            [0.00025, 0.0, 0.0, 0.0, 0.00025, 0.0, 0.0, 0.0]),
    "info" => Dict(
        "name" => "motor",
        "dim" => 8,
        "linear" => true,
        "hybrid" => false
    )
])

max_stars_repo_name: mforets/ReachabilityLibrary

max_issues_repo_name: mforets/ReachabilityLibrary

max_forks_repo_name: mforets/ReachabilityLibrary


text:
import Distributions.pdf, Distributions.logpdf, Random.rand, Base.length

struct MixedSupport <: ValueSupport; end

"""
    Factored{N} <: Distribution{Multivariate, MixedSupport}

a `Distribution` type that can be used to combine multiple `UnivariateDistribution`'s and sample from them.

Example: it can be used as `prior = Factored(Normal(0,1), Uniform(-1,1))`
"""
struct Factored{N}<:Distribution{Multivariate,MixedSupport}
    p::NTuple{N,UnivariateDistribution}
    Factored(args::UnivariateDistribution...) = new{length(args)}(args)
end
"""
    pdf(d::Factored, x) = begin

Function to evaluate the pdf of a `Factored` distribution object
"""
function pdf(d::Factored{N},x) where N
    s=pdf(d.p[1],x[1])
    for i in 2:N
        s*=pdf(d.p[i],x[i])
    end
    s
end

"""
    logpdf(d::Factored, x) = begin

Function to evaluate the logpdf of a `Factored` distribution object
"""
function logpdf(d::Factored{N},x) where N
    s=logpdf(d.p[1],x[1])
    for i in 2:N
        s+=logpdf(d.p[i],x[i])
    end
    s
end

"""
    rand(rng::AbstractRNG, factoreddist::Factored)

function to sample one element from a `Factored` object
"""
rand(rng::AbstractRNG,factoreddist::Factored{N}) where N = ntuple(i->rand(rng,factoreddist.p[i]),Val(N))

"""
    length(p::Factored) = begin

returns the number of distributions contained in `p`.
"""
length(p::Factored{N}) where N = N

export Factored


text:
using Revise
using ComplexMixtures, PDBTools
const CM = ComplexMixtures

nthreads = Threads.nthreads()
if nthreads == 1
  error(" These tests have to be run with multiple threads as well. ")
end

dir="$(@__DIR__ )/data/NAMD"
atoms = readPDB("$dir/structure.pdb")
protein = Selection(select(atoms,"protein"),nmols=1)
tmao = Selection(select(atoms,"resname TMAO"),natomspermol=14)
water = Selection(select(atoms,"water"),natomspermol=3)

options_single = Options(stride=1,seed=321,nthreads=1,silent=true,lcell=2)
options_multi = Options(stride=1,seed=321,silent=true,lcell=2)

println(" --------------------------------------------------------------")
println(" Compiling - single thread ")
traj = Trajectory("$dir/trajectory.dcd",protein,tmao)
@time R = mddf(traj,options_single);
println(" Compiling - multi-threaded ")
traj = Trajectory("$dir/trajectory.dcd",protein,tmao)
@time R = mddf(traj,options_multi);
println(" Compiling - self single thread ")
traj = Trajectory("$dir/trajectory.dcd",tmao)
@time R = mddf(traj,options_single);
println(" Compiling - self multi-threaded ")
traj = Trajectory("$dir/trajectory.dcd",tmao)
@time R = mddf(traj,options_multi);

println(" --------------------------------------------------------------")
println(" Protein-TMAO - Single threaded ")
traj = Trajectory("$dir/trajectory.dcd",protein,tmao)
@time R = mddf(traj,options_single);

println(" --------------------------------------------------------------")
println(" Protein-TMAO - nthreads = $nthreads ")
traj = Trajectory("$dir/trajectory.dcd",protein,tmao)
@time R = mddf(traj,options_multi);

println(" --------------------------------------------------------------")
println(" Protein-Water - Single threaded ")
traj = Trajectory("$dir/trajectory.dcd",protein,water)
@time R = mddf(traj,options_single);

println(" --------------------------------------------------------------")
println(" Protein-Water - nthreads = $nthreads ")
traj = Trajectory("$dir/trajectory.dcd",protein,water)
@time R = mddf(traj,options_multi);

println(" --------------------------------------------------------------")
println(" Water-TMAO - Single threaded ")
traj = Trajectory("$dir/trajectory.dcd",water,tmao)
@time R = mddf(traj,options_single);

println(" --------------------------------------------------------------")
println(" Water-TMAO - nthreads = $nthreads ")
traj = Trajectory("$dir/trajectory.dcd",water,tmao)
@time R = mddf(traj,options_multi);

println(" --------------------------------------------------------------")
println(" Water-Water - Single threaded ")
traj = Trajectory("$dir/trajectory.dcd",water)
@time R = mddf(traj,options_single);

println(" --------------------------------------------------------------")
println(" Water-Water - nthreads = $nthreads ")
traj = Trajectory("$dir/trajectory.dcd",water)
@time R = mddf(traj,options_multi);
println(" --------------------------------------------------------------")


text:
"""
qqmap(obs::ClimGrid, ref::ClimGrid, fut::ClimGrid; method="Additive", detrend=true, window::Int=15, rankn::Int=50, thresnan::Float64=0.1, keep_original::Bool=false, interp::Function = Linear(), extrap::Function = Flat())

Quantile-Quantile mapping bias correction.

For each julian day of the year (+/- **window** size), a transfer function is estimated through an empirical quantile-quantile mapping. The quantile-quantile transfer function between **ref** and **obs** is etimated on a julian day (and grid-point) basis with a moving window around the julian day. Hence, for a given julian day, the transfer function is then applied to the **fut** dataset for a given julian day.

**Options**

**method::String = "Additive" (default) or "Multiplicative"**. Additive is used for most climate variables. Multiplicative is usually bounded variables such as precipitation and humidity.

**detrend::Bool = true (default)**. A 4th order polynomial is adjusted to the time series and the residuals are corrected with the quantile-quantile mapping.

**window::Int = 15 (default)**. The size of the window used to extract the statistical characteristics around a given julian day.

**rankn::Int = 50 (default)**. The number of bins used for the quantile estimations. The quantiles uses by default 50 bins between 0.01 and 0.99. The bahavior between the bins is controlled by the interp keyword argument. The behaviour of the quantile-quantile estimation outside the 0.01 and 0.99 range is controlled by the extrap keyword argument.

**thresnan::Float64 = 0.1 (default)**. The fraction is missing values authorized for the estimation of the quantile-quantile mapping for a given julian days. If there is more than **treshnan** missing values, the output for this given julian days returns NaNs.

**keep_original::Bool = false (default)**. If **keep_original** is set to true, the values are set to the original values in presence of too many NaNs.

**interp = Interpolations.Linear() (default)**. When the data to be corrected lies between 2 quantile bins, the value of the transfer function is linearly interpolated between the 2 closest quantile estimation. The argument is from Interpolations.jl package.

**extrap = Interpolations.Flat() (default)**. The bahavior of the quantile-quantile transfer function outside the 0.01-0.99 range. Setting it to Flat() ensures that there is no "inflation problem" with the bias correction. The argument is from Interpolation.jl package.
"""
function qqmap(obs::ClimGrid, ref::ClimGrid, fut::ClimGrid; method::String="Additive", detrend::Bool=true, window::Int=15, rankn::Int=50, qmin::Real=0.01, qmax::Real=0.99, thresnan::Float64=0.1, keep_original::Bool=false, interp=Linear(), extrap=Flat())

    # Consistency checks # TODO add more checks for grid definition
    @argcheck size(obs[1], 1) == size(ref[1], 1) == size(fut[1], 1)
    @argcheck size(obs[1], 2) == size(ref[1], 2) == size(fut[1], 2)

    # Modify dates (e.g. 29th feb are dropped/lost by default)
    obs = ClimateTools.dropfeb29(obs)
    ref = ClimateTools.dropfeb29(ref)
    fut = ClimateTools.dropfeb29(fut)

    # Remove trend if specified
    if detrend == true
        # Obs
        obs_polynomials = ClimateTools.polyfit(obs)
        obs = obs - ClimateTools.polyval(obs, obs_polynomials)
        # Ref
        ref_polynomials = ClimateTools.polyfit(ref)
        ref = ref - ClimateTools.polyval(ref, ref_polynomials)
        # Fut
        fut_polynomials = ClimateTools.polyfit(fut)
        poly_values = ClimateTools.polyval(fut, fut_polynomials)
        fut = fut - poly_values
    end

    #Get date vectors
    datevec_obs = get_timevec(obs)
    datevec_ref = get_timevec(ref)
    datevec_fut = get_timevec(fut)

    # Julian days vectors
    obs_jul = Dates.dayofyear.(datevec_obs)
    ref_jul = Dates.dayofyear.(datevec_ref)
    fut_jul = Dates.dayofyear.(datevec_fut)

    # Prepare output array (replicating data type of fut ClimGrid)
    dataout = fill(convert(typeof(fut[1].data[1]), NaN), (size(fut[1], 1), size(fut[1], 2), size(fut[1], 3)))::Array{typeof(fut[1].data[1]),T} where T

    if minimum(ref_jul) == 1 && maximum(ref_jul) == 365
        days = 1:365
    else
        days = minimum(ref_jul) + window:maximum(ref_jul) - window
        start = Dates.monthday(minimum(datevec_obs))
        finish = Dates.monthday(maximum(datevec_obs))
    end

    # Reshape. Allows multi-threading on grid points.
    obsin = reshape(obs[1].data, (size(obs[1].data, 1) * size(obs[1].data, 2), size(obs[1].data, 3)))
    refin = reshape(ref[1].data, (size(ref[1].data, 1) * size(ref[1].data, 2), size(ref[1].data, 3)))
    futin = reshape(fut[1].data, (size(fut[1].data, 1) * size(fut[1].data, 2), size(fut[1].data, 3)))
    dataoutin = reshape(dataout, (size(dataout, 1) * size(dataout, 2), size(dataout, 3)))

    # Looping over grid points using multiple-dispatch calls to qqmap
    Threads.@threads for k = 1:size(obsin, 1)
        # for k = 1:size(obsin, 1)

        obsvec = obsin[k,:]
        refvec = refin[k,:]
        futvec = futin[k,:]

        dataoutin[k, :] = qqmap(obsvec, refvec, futvec, days, obs_jul, ref_jul, fut_jul, method=method, detrend=detrend, window=window, rankn=rankn, qmin=qmin, qmax=qmax, thresnan=thresnan, keep_original=keep_original, interp=interp, extrap=extrap)

    end

    # Apply mask
    dataout = applymask(dataout, obs.msk)

    # Rebuilding ClimGrid
    lonsymbol = Symbol(fut.dimension_dict["lon"])
    latsymbol = Symbol(fut.dimension_dict["lat"])

    dataout2 = AxisArray(dataout, Axis{lonsymbol}(fut[1][Axis{lonsymbol}].val), Axis{latsymbol}(fut[1][Axis{latsymbol}].val), Axis{:time}(datevec_fut))

    timeattrib = fut.timeattrib
    timeattrib["calendar"] = "365_day"

    globalattribs = fut.globalattribs
    if haskey(globalattribs, "history")
        globalattribs["history"] = string(globalattribs["history"], " - ", "Bias-corrected with ClimateTools.jl")
    else
        globalattribs["history"] = string("Bias-corrected with ClimateTools.jl")
    end

    C = ClimGrid(dataout2; longrid = fut.longrid, latgrid = fut.latgrid, msk = fut.msk, grid_mapping = fut.grid_mapping, dimension_dict = fut.dimension_dict, timeattrib = timeattrib, model = fut.model, frequency = fut.frequency, experiment = fut.experiment, run = fut.run, project = fut.project, institute = fut.institute, filename = fut.filename, dataunits = fut.dataunits, latunits = fut.latunits, lonunits = fut.lonunits, variable = fut.variable, typeofvar = fut.typeofvar, typeofcal = "365_day", varattribs = fut.varattribs, globalattribs = globalattribs)

    if detrend == true
        C = C + poly_values
    end

    return C

end

"""
qqmap(obsvec::Array{N,1}, refvec::Array{N,1}, futvec::Array{N,1}, days, obs_jul, ref_jul, fut_jul; method::String="Additive", detrend::Bool=true, window::Int64=15, rankn::Int64=50, thresnan::Float64=0.1, keep_original::Bool=false, interp=Linear(), extrap=Flat())

Quantile-Quantile mapping bias correction for single vector. This is a low level function used by qqmap(A::ClimGrid ..), but can work independently.

"""
function qqmap(obsvec::Array{N,1} where N, refvec::Array{N,1} where N, futvec::Array{N,1} where N, days, obs_jul, ref_jul, fut_jul; method::String="Additive", detrend::Bool=true, window::Int64=15, rankn::Int64=50, qmin::Real=0.01, qmax::Real=0.99, thresnan::Float64=0.1, keep_original::Bool=false, interp=Linear(), extrap=Flat())

    # range over which quantiles are estimated
    P = range(qmin, stop=qmax, length = rankn)
    obsP = similar(obsvec, length(P))
    refP = similar(refvec, length(P))
    sf_refP = similar(refvec, length(P))

    # Prepare output array
    dataout = similar(futvec, (size(futvec)))

    # LOOP OVER ALL DAYS OF THE YEAR
    for ijulian in days

        # idx for values we want to correct
        idxfut = (fut_jul .== ijulian)

        # Find all index of moving window around ijulian day of year
        idxobs = ClimateTools.find_julianday_idx(obs_jul, ijulian, window)
        idxref = ClimateTools.find_julianday_idx(ref_jul, ijulian, window)

        obsval = obsvec[idxobs]
        refval = refvec[idxref]
        futval = futvec[idxfut]

        # Random perturbation is added to avoid division by zero
        if lowercase(method) == "multiplicative"
            # Seed for reproducible results
            Random.seed!(42)

            init     = 0.000001
            interval = 0.00001
            stop     = 0.001
            R = init:interval:stop

            obsval = obsval .+ rand(R, length(obsvec[idxobs]))
            refval = refval .+ rand(R, length(refvec[idxref]))
            futval = futval .+ rand(R, length(futvec[idxfut]))
        end

        obsvalnan = isnan.(obsval)
        refvalnan = isnan.(refval)
        futvalnan = isnan.(futval)

        if (sum(obsvalnan) < (length(obsval) * thresnan)) & (sum(refvalnan) < (length(refval) * thresnan)) & (sum(futvalnan) < (length(futval) * thresnan))

            # Estimate quantiles for obs and ref for ijulian
            quantile!(obsP, obsval[.!obsvalnan], P, sorted = false)
            quantile!(refP, refval[.!refvalnan], P, sorted = false)

            if lowercase(method) == "additive" # used for temperature
                sf_refP .= obsP .- refP
                itp = interpolate((refP,), sf_refP, Gridded(interp))
                itp = extrapolate(itp, extrap) # add extrapolation
                # dataout[idxfut] .= itp(futval) .+ futval
                futval .= itp(futval) .+ futval

                futval[futvalnan] .= 0.0

            elseif lowercase(method) == "multiplicative" # used for precipitation
                sf_refP .= obsP ./ refP
                sf_refP[sf_refP .< 0] .= eps(1.0)
                sf_refP[isnan.(sf_refP)] .= 1.0#eps(1.0)#1.0
                itp = interpolate((refP,), sf_refP, Gridded(interp))
                itp = extrapolate(itp, extrap) # add extrapolation
                futval .= itp(futval) .* futval

                futval[futvalnan] .= 0.0

            else
                error("Wrong method")
            end
            # Replace values with new ones
            dataout[idxfut] .= futval
        else

            if keep_original
                # # Replace values with original ones (i.e. too may NaN values for robust quantile estimation)
                dataout[idxfut] .= futval
            else
                # DO NOTHING (e.g. if there is no reference, we want NaNs and not original values)
                dataout[idxfut] .= NaN
            end
        end
    end

    return dataout

end

"""
    biascorrect_extremes(obs::ClimGrid, ref::ClimGrid, fut::ClimGrid; detrend=false, window::Int=15, rankn::Int=50, thresnan::Float64=0.1, keep_original::Bool=false, interp=Linear(), extrap=Flat(), gev_params::DataFrame, frac=0.25, power=1.0)

Combine the empirial Quantile-Quantile mapping (see [`qqmap`](@ref)) and Generalized Pareto Distribution bias-correction methods.

The tail of the distribution is corrected with a paramatric GPD, using provided parameters μ, σ and ξ contained in gevparams DataFrame. The DataFrame gevparams has column :lat, :lon, :mu, :sigma and :xi, representing the GEV parameters and the spatial location of the parameters. For each grid point, the function extracts the closest GEV parameters.

**Options specific to GPD**

**gevparams::DataFrame** represents the (external) GEV parameters to correct each grid-points.

**frac=0.25** is the fraction where the cutoff happens between QQM and GPD, as defined by **(maximum(x) - minimum(x))*frac** (e.g. For a maximum value of 150mm and a minimum value of 30mm, the linear transition will be between 30mm and 60mm).

**power=1.0** is the shape of the transition.

**Options specific to QQM**

**detrend::Bool = false (default)**. A 4th order polynomial is adjusted to the time series and the residuals are corrected.

**window::Int = 15 (default)**. The size of the window used to extract the statistical characteristics around a given julian day.

**rankn::Int = 50 (default)**. The number of bins used for the quantile estimations. The quantiles uses by default 50 bins between 0.01 and 0.99. The bahavior between the bins is controlled by the interp keyword argument. The behaviour of the quantile-quantile estimation outside the 0.01 and 0.99 range is controlled by the extrap keyword argument.

**thresnan::Float64 = 0.1 (default)**. The fraction is missing values authorized for the estimation of the quantile-quantile mapping for a given julian days. If there is more than **treshnan** missing values, the output for this given julian days returns NaNs.

**keep_original::Bool = false (default)**. If **keep_original** is set to true, the values are set to the original values in presence of too many NaNs.

**interp = Interpolations.Linear() (default)**. When the data to be corrected lies between 2 quantile bins, the value of the transfer function is linearly interpolated between the 2 closest quantile estimation. The argument is from Interpolations.jl package.

**extrap = Interpolations.Flat() (default)**. The bahavior of the quantile-quantile transfer function outside the 0.01-0.99 range. Setting it to Flat() ensures that there is no "inflation problem" with the bias correction. The argument is from Interpolation.jl package.

See also: [`qqmap`](@ref)
"""
function biascorrect_extremes(obs::ClimGrid, ref::ClimGrid, fut::ClimGrid; detrend=false, P::Real=0.95, window::Int=15, rankn::Int=50, qmin::Real=0.01, qmax::Real=0.99, thresnan::Float64=0.1, keep_original::Bool=false, interp=Linear(), extrap=Flat(), gevparams::DataFrame=DataFrame(), frac=0.25, power=1.0)

    # Consistency checks # TODO add more checks for grid definition
    @argcheck size(obs[1], 1) == size(ref[1], 1) == size(fut[1], 1)
    @argcheck size(obs[1], 2) == size(ref[1], 2) == size(fut[1], 2)

    qqmap_base = qqmap(obs, ref, fut, method="multiplicative", detrend=detrend, window=window, rankn=rankn, qmin=qmin, qmax=qmax, thresnan=thresnan, keep_original=keep_original, interp=interp, extrap=extrap)

    # Modify dates (e.g. 29th feb are dropped/lost by default)
    obs = ClimateTools.dropfeb29(obs)
    ref = ClimateTools.dropfeb29(ref)
    fut = ClimateTools.dropfeb29(fut)

    # Remove trend if specified
    if detrend == true
        # Obs
        obs_polynomials = ClimateTools.polyfit(obs)
        obs = obs - ClimateTools.polyval(obs, obs_polynomials)
        # Ref
        ref_polynomials = ClimateTools.polyfit(ref)
        ref = ref - ClimateTools.polyval(ref, ref_polynomials)
        # Fut
        fut_polynomials = ClimateTools.polyfit(fut)
        poly_values = ClimateTools.polyval(fut, fut_polynomials)
        fut = fut - poly_values
    end

    #Get date vectors
    datevec_obs = get_timevec(obs)
    datevec_ref = get_timevec(ref)
    datevec_fut = get_timevec(fut)

    # Condition for a moving window
    movingwindow = maximum(year.(datevec_fut)) - minimum(year.(datevec_fut)) > (maximum(year.(datevec_ref)) - minimum(year.(datevec_ref)))*1.5

    # Prepare output array (replicating data type of fut ClimGrid)
    # dataout = fill(convert(typeof(fut[1].data[1]), NaN), (size(fut[1], 1), size(fut[1], 2), size(fut[1], 3)))::Array{typeof(fut[1].data[1]),T} where T
    dataout = qqmap_base[1].data

    # Reshape. Allows multi-threading on grid points.
    obsin = reshape(obs[1].data, (size(obs[1].data, 1) * size(obs[1].data, 2), size(obs[1].data, 3)))
    refin = reshape(ref[1].data, (size(ref[1].data, 1) * size(ref[1].data, 2), size(ref[1].data, 3)))
    futin = reshape(fut[1].data, (size(fut[1].data, 1) * size(fut[1].data, 2), size(fut[1].data, 3)))
    dataoutin = reshape(dataout, (size(dataout, 1) * size(dataout, 2), size(dataout, 3)))
    latgrid = reshape(obs.latgrid, (size(obs.latgrid, 1) * size(obs.latgrid, 2)))
    longrid = reshape(obs.longrid, (size(obs.longrid, 1) * size(obs.longrid, 2)))
    qqvecin = reshape(qqmap_base[1].data, (size(qqmap_base[1].data, 1) * size(qqmap_base[1].data, 2), size(qqmap_base[1].data, 3)))

    # Looping over grid points using multiple-dispatch calls to qqmap
    # Threads.@threads for k = 1:size(obsin, 1)
    for k = 1:size(obsin, 1)

        obsvec = obsin[k,:]
        refvec = refin[k,:]
        futvec = futin[k,:]

        obsvalnan = isnan.(obsvec)
        refvalnan = isnan.(refvec)
        futvalnan = isnan.(futvec)

        if (sum(obsvalnan) < (length(obsvec) * thresnan)) & (sum(refvalnan) < (length(refvec) * thresnan)) & (sum(futvalnan) < (length(futvec) * thresnan))

            # Estimate threshold
            thres = ClimateTools.get_threshold(obsvec, refvec, thres=P)

            if isempty(gevparams)
                # If empty, we estimate the GPD on the reference dataset
                obsclusters = getcluster(obsvec, thres, 1.0)
                # Estimate GPD parameters of fut clusters
                GPD_obs = gpdfit(obsclusters[!,:Max], threshold = thres)

            else
                GPD_obs = ClimateTools.estimate_gpd(latgrid[k], longrid[k], gevparams, thres)
            end

            if movingwindow

                idxi, idxf, idxi_corr, idxf_corr = build_idx_biascorrect(refvec, futvec)

                for c = 1:length(idxi)

                    futvec_tmp = futvec[idxi[c]:idxf[c]]

                    # Get clusters
                    futclusters = getcluster(futvec_tmp, thres, 1.0)
                    # Estimate GPD parameters of fut clusters
                    GPD_fut = gpdfit(futclusters[!,:Max], threshold = thres)
                    # Calculate quantile values of maximum clusters values
                    fut_cdf = Extremes.cdf.(GPD_fut, futclusters[!,:Max])
                    # Estimation of fut values based on provided GPD parameters
                    newfut = quantile.(GPD_obs, fut_cdf)

                    # Linear transition over a quarter of the distance
                    exIDX = futclusters[!,:Position]

                    # Get the weights based on frac and power values
                    transition = ClimateTools.weight(futvec[exIDX], frac=frac, power=power)

                    # Apply linear transition
                    newfut_trans = (newfut .* transition) .+ (qqvecin[k, idxi[c] .+ exIDX .- 1] .* (1.0 .- transition))

                    # We put the new data in the bias-corrected vector
                    exIDX_corr = findall((idxi[c] .+ futclusters[!,:Position] .-1 .< idxf_corr[c]) .& (idxi[c] .+ futclusters[!,:Position] .- 1 .> idxi_corr[c]))
                    # exIDX_corr = findall(futclusters[!,:Position] .< idxf_corr)
                    # dataoutin[k, idxi+idxi_corr-1:idxi+idxf_corr-1] = qqvecin[k, idxi_corr:idxf_corr]
                    dataoutin[k, idxi[c] .+ futclusters[exIDX_corr,:Position] .- 1] .= newfut_trans[exIDX_corr]

                end
            else
                # refclusters = getcluster(refvec[k,:], thres, 1.0)
                futclusters = getcluster(futvec, thres, 1.0)
                GPD_fut = gpdfit(futclusters[!,:Max], threshold = thres)
                fut_cdf = Extremes.cdf.(GPD_fut, futclusters[!,:Max])
                newfut = quantile.(GPD_obs, fut_cdf)

                # Linear transition over a quarter of the distance
                exIDX = futclusters[!,:Position]

                # Get the weights based on frac and power values
                transition = ClimateTools.weight(futvec[exIDX], frac=frac, power=power)

                # Apply linear transition
                newfut_trans = (newfut .* transition) .+ (qqvecin[k, exIDX] .* (1.0 .- transition))

                # We put the new data in the bias-corrected vector
                # dataoutin[k, :] = qqvecin[k, :]
                dataoutin[k, futclusters[!,:Position]] = newfut_trans

            end

        # else
        #
        #     if keep_original
        #         # # Replace values with original ones (i.e. too may NaN values for robust quantile estimation)
        #         dataoutin[k, idxfut] .= futvec
        #     else
        #         # DO NOTHING (e.g. if there is no reference, we want NaNs and not original values)
        #         dataoutin[k, idxfut] .= NaN
        #     end
        end
    end

    # Apply mask
    dataout = applymask(dataout, obs.msk)

    # Rebuilding ClimGrid
    lonsymbol = Symbol(fut.dimension_dict["lon"])
    latsymbol = Symbol(fut.dimension_dict["lat"])

    dataout2 = AxisArray(dataout, Axis{lonsymbol}(fut[1][Axis{lonsymbol}].val), Axis{latsymbol}(fut[1][Axis{latsymbol}].val), Axis{:time}(datevec_fut))

    timeattrib = fut.timeattrib
    timeattrib["calendar"] = "365_day"

    globalattribs = fut.globalattribs
    if haskey(globalattribs, "history")
        globalattribs["history"] = string(globalattribs["history"], " - ", "Bias-corrected with ClimateTools.jl")
    else
        globalattribs["history"] = string("Bias-corrected for extreme values with ClimateTools.jl")
    end

    C = ClimGrid(dataout2; longrid = fut.longrid, latgrid = fut.latgrid, msk = fut.msk, grid_mapping = fut.grid_mapping, dimension_dict = fut.dimension_dict, timeattrib = timeattrib, model = fut.model, frequency = fut.frequency, experiment = fut.experiment, run = fut.run, project = fut.project, institute = fut.institute, filename = fut.filename, dataunits = fut.dataunits, latunits = fut.latunits, lonunits = fut.lonunits, variable = fut.variable, typeofvar = fut.typeofvar, typeofcal = "365_day", varattribs = fut.varattribs, globalattribs = globalattribs)

    if detrend == true
        C = C + poly_values
    end

    return C

end

"""
    gev2gpd(μ, σ, ξ, thres)

Transform GEV parameters to GPD parameters (Coles, 2001)
"""
function gev2gpd(μ, σ, ξ, thres)

    return σ₂ = σ + ξ * (thres - μ)
end

"""
    weight(x; frac=0.25, power=1.0)

Returns the associated weight used in weighting Quantile-Quantile mapping (QQM) and Generalized Pareto Distribution (GPD) bias-correction methods. Used internally by biascorrect_extremes.

*frac* is the fraction where the cutoff happens between QQM and GPD, as defined by **(maximum(x) - minimum(x))*frac** (e.g. For a maximum value of 150mm and a minimum value of 30mm, the linear transition will be between 30mm and 60mm) and power is the shape of the transition.
"""
function weight(x; frac=0.25, power=1.0)
    ω = ((x .- minimum(x))/((maximum(x) - minimum(x))*frac)).^power
    ω[ω .> 1.0] .= 1.0
    return ω
end

function estimate_gpd(lat, lon, gevparams, thres)

    # Find closest gev parameters
    l1 = (Float64(lat), Float64(lon))
    l2 = (Array(gevparams[!,:lat]), Array(gevparams[!,:lon]))
    dist, idx = findmindist(l1, l2)

    # Extract GEV parameters
    μ = gevparams[idx[1], :mu]
    σ = gevparams[idx[1], :sigma]
    ξ = gevparams[idx[1], :xi]

    # Transform to GPD parameters and build reference GP distribution.
    μ₂ = ClimateTools.gev2gpd(μ, σ, ξ, thres)
    GPD_obs = Extremes.GeneralizedPareto(μ₂, σ, ξ)

    return GPD_obs

end

function build_idx_biascorrect(refvec, futvec)

    # window_length will be roughly the length of REF
    w_length = length(refvec)#(maximum(year.(datevec_ref)) - minimum(year.(datevec_ref)) + 1)*365
    # in case of multiple members, we approximate ref length
    fut_l = length(futvec)

    # The fut series is separated is a few time windows.
    # We also make the windows overlap to replicate a moving window
    R = w_length/2:w_length/3:fut_l
    corrR = range(1, stop=fut_l, length=length(R)-1)

    idxi = Array{Int64}(undef, length(R))
    idxf = Array{Int64}(undef, length(R))
    idxi_corr = Array{Int64}(undef, length(R))
    idxf_corr = Array{Int64}(undef, length(R))

    for c = 1:length(R)

        # idxi = nodes_est[c,1]
        # idxf = nodes_est[c,2]
        # idxi_corr = nodes_corr[c]
        # idxf_corr = nodes_corr[c+1]

        idxi[c] = round(Int, max(R[c]-w_length/2+1,1))
        idxf[c] = round(Int, min(R[c]+w_length/2,fut_l))
        # futvec_tmp = futvec[idxi:idxf];

        # However, we want each moving window to correct a smaller portion of the time series
        idxi_corr[c] = round(Int, idxi[c]+(w_length/2-(w_length/2)/3) - 1)
        # idxi_corr[c] = round(Int, corrR[c])
        if c == 1
            idxi_corr[c] = 1
        end
        # idxf_corr = round(Int, idxi + (w_length/2+(w_length/2)/3) - 1)

        if c == length(R)#maximum(w_length/2:w_length/3:fut_l) || idxi+idxf_corr-1 > fut_l
            idxf_corr[c] = max(round(Int, idxi[c] + w_length-w_length/3),fut_l)
        else
            idxf_corr[c] = round(Int, idxi[c] + (w_length/2+(w_length/2)/3) - 1)
            # idxf_corr[c] = round(Int,corrR[c+1])
        end
        # @show idxi
        # @show idxf
        # @show idxi_corr
        # @show idxf_corr
    end

    idxi_corr[2:end] .-= 1

    return idxi, idxf, idxi_corr, idxf_corr
end

#     w_length = length(refvec)
#     fut_l = length(futvec)
#     nodes_corr = round.(Int,w_length/2:w_length/3:fut_l-w_length/2)
#     N = Array{Int64}(undef, length(nodes_corr)+2)
#     N[1] = 1
#     N[end] = fut_l
#     N[2:end-1] .= nodes_corr

#     N_est = Array{Int64}(undef, length(N), 2)

#     for inode = 1:length(N)

#         if inode == 1
#             N_est[inode, 1] = N[1]
#             N_est[inode, 2] = N[3]
#         elseif inode == length(N)
#             N_est[inode, 1] = N[end-1]
#             N_est[inode, 2] = N[end]
#         # elseif inode == length(N) - 1
#         #     N_est[inode, 1] = N[end-3]
#         #     N_est[inode, 2] = N[end-1]
#         else
#             N_est[inode, 1] = N[inode - 1]
#             N_est[inode, 2] = N[inode + 1]
#         end
#     end

#     return N, N_est
# end

#     # N[1] = 1
#     # N[end] = fut_l

#     for inode = 1:length(nodes_corr)
#         if inode == 1
#             N[inode, 1] = 1
#             N[inode, 2] = nodes_corr[inode]
#         elseif inode == length(nodes_corr)
#             N[inode] = fut_l
#         else
#             N[inode] = 0
#         end

#     end

#     nodes_corr = collect(Base.Iterators.partition(w_length/2:fut_l, round(Int,w_length/3)))
#     nodes_est = collect(Base.Iterators.partition(1:fut_l, round(Int,w_length)))

#     Ni = Vector{Int64}(undef, length(nodes_corr))
#     Nf = Vector{Int64}(undef, length(nodes_corr))

#     for inode = 1:length(nodes_corr)

#         for icorr = 1:length(nodes_est)

#             idxi = findall(nodes_corr[inode][1] .== nodes_est[icorr])
#             idxf = findall(nodes_corr[inode][end] .== nodes_est[icorr])

#             if !isempty(idxi)
#                 Ni[inode] = icorr
#             end
#             if !isempty(idxf)
#                 Nf[inode] = icorr
#             end
#         end
#     end

#     # @assert sum(Ni - Nf) == 0

#     return nodes_corr, nodes_est, Ni, Nf
# end

# for n = 1:30000
#     refvec = rand(10950)
#     futvec = rand(10950 + n)

#     nodes_corr, nodes_est, Ni, Nf = build_idx_biascorrect(refvec, futvec)
# end

# """
#     biascorrect_extremes(obs::AxisArray, ref::AxisArray, fut::AxisArray, μ, σ, ξ)
#
# Correct the tail of the distribution with a paramatric method, using the parameters μ, σ and ξ.
# """
# function biascorrect_extremes(obsvec::AxisArray, refvec::AxisArray, futvec::AxisArray, days, μ::Real, σ::Real, ξ::Real)
#
#     obs_jul = dayofyear.(obsvec[Axis{:time}][:])
#     ref_jul = dayofyear.(refvec[Axis{:time}][:])
#     fut_jul = dayofyear.(futvec[Axis{:time}][:])
#
#     # Prepare output array (replicating data type of fut ClimGrid)
#     dataout = fill(convert(typeof(futvec.data[1]), NaN), length(futvec))::Array{typeof(fut[1].data[1]),T} where T
#
#     if minimum(ref_jul) == 1 && maximum(ref_jul) == 365
#         days = 1:365
#     else
#         days = minimum(ref_jul) + window:maximum(ref_jul) - window
#         start = Dates.monthday(minimum(datevec_obs))
#         finish = Dates.monthday(maximum(datevec_obs))
#     end
#
#     threshold = mean([quantile(obsvec[obsvec .>= 1.0], P) quantile(refvec[refvec .>= 1.0], P)])
#
#     idx_obs = obsvec .> threshold
#     idx_ref = refvec .> threshold
#     idx_fut = futvec .> threshold
#
#     dataout .= qqmap(obsvec.data, refvec.data, futvec.data, days, obs_jul, ref_jul, fut_jul, method = "multiplicative", detrend = false)
#
#
# end


text:
export mnist_image

using Flux.Data.MNIST

function mnist_image(num_train, batch_size, num_test)
	X = float.(hcat(vec.(MNIST.images())...)) .> 0.5
	N = size(X, 2)
	I_train = rand(1:N, num_train); train_X = X[:, I_train]
	I_remain = [i for i in 1:N if !(i∈I_train)]
	I_test = rand(I_remain, num_test); test_X = X[:, I_test]
	batch_data = [train_X[:,i] for i in Iterators.partition(1:num_train, batch_size)]
	batch_data, test_X
end

max_stars_repo_name: HamletWantToCode/Generative

max_issues_repo_name: HamletWantToCode/Generative

max_forks_repo_name: HamletWantToCode/Generative


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

using Materials, Test

@testset "Test Materials.jl" begin
    # @testset "test abstract material" begin
    #     include("test_abstract_material.jl")
    # end
    @testset "test ideal plastic material model" begin
        include("test_idealplastic.jl")
    end
    @testset "test ideal plastic pure shear" begin
        include("test_idealplastic_shear.jl")
    end
    @testset "test uniaxial increment" begin
        include("test_uniaxial_increment.jl")
    end
    @testset "test chaboche pure shear" begin
        include("test_chaboche_shear.jl")
    end
    @testset "test chaboche uniaxial stress" begin
        include("test_chaboche.jl")
    end
    @testset "test biaxial increment" begin
        include("test_biaxial_increment.jl")
    end
    # @testset "test viscoplastic material model" begin
    #     include("test_viscoplastic.jl")
    # end
end


text:
module TestWithOffset

using Base: HasEltype, HasLength, HasShape, IteratorEltype, IteratorSize
using SplittablesBase.Implementations: WithOffset
using Test

@testset "eltype" begin
    @test IteratorEltype(WithOffset(1:2, 0)) isa HasEltype
    @test eltype(WithOffset(1:2, Int128(0))) === Int128
end

@testset "size" begin
    xs = eachindex("abc")
    @test IteratorSize(WithOffset(xs, 1)) isa HasLength
    @test length(WithOffset(xs, 1)) == 3
    @test collect(WithOffset(xs, 1)) == xs .+ 1
end

@testset "size" begin
    xs = reshape(1:6, 2, 3)
    @test IteratorSize(WithOffset(xs, 0)) isa HasShape
    @test size(WithOffset(xs, 1)) == (2, 3)
    @test length(WithOffset(xs, 1)) == 6
    @test collect(WithOffset(xs, 1)) == xs .+ 1
end

end  # module


text:
# CUDAnative run-time library
#
# This module defines method instances that will be compiled into a device-specific image
# and will be available to the CUDAnative compiler to call after Julia has generated code.

module Runtime

using ..CUDAnative
using LLVM

struct MethodInstance
    def::Function
    return_type::Type
    types::Tuple
    name::String
end

instantiate(def, return_type, types, name=String(typeof(def).name.mt.name)) =
    MethodInstance(def, return_type, types, name)

function Base.getproperty(rt::MethodInstance, field::Symbol)
    # overloaded field accessor to get LLVM types (this only works at run-time,
    # ie. after the method instance has been constructed and precompiled)
    if field == :llvm_types
        LLVMType[convert.(LLVMType, typ) for typ in rt.types]
    elseif field == :llvm_return_type
        convert.(LLVMType, rt.return_type)
    else
        return getfield(rt, field)
    end
end

## exception handling

function ptx_report_exception(ex)
    @cuprintf("""
        ERROR: a %s exception occurred during kernel execution.
               Run Julia on debug level 2 for device stack traces.\n""", ex)
    return
end

const report_exception = instantiate(ptx_report_exception, Nothing, (Ptr{Cchar},))

function ptx_report_exception_name(ex)
    @cuprintf("ERROR: a %s exception occurred during kernel execution.\n", ex)
    return
end

const report_exception_name = instantiate(ptx_report_exception_name, Nothing, (Ptr{Cchar},))

function ptx_report_exception_frame(idx, func, file, line)
    @cuprintf(" [%i] %s at %s:%i\n", idx, func, file, line)
    return
end

const report_exception_frame = instantiate(ptx_report_exception_frame, Nothing,
                                           (Cint, Ptr{Cchar}, Ptr{Cchar}, Cint))

end


text:
type Gamma
    _single::Vector{Int}
    _local::Vector{Int}
    _regional::Vector{Int}
    _global::Int
end

type Qrobust
    m::JuMP.Model
    I::UnitRange{Int}
    J::UnitRange{Int}
    d::Vector{JuMP.Variable}
    q::Vector{JuMP.Variable}
    γ::Gamma
end

type RobustDeployment <: DeploymentModel
    m::JuMP.Model
    Q::Qrobust

    I::UnitRange{Int}
    J::UnitRange{Int}

    x::Vector{JuMP.Variable}
    y::Vector{Matrix{JuMP.Variable}}
    z::Vector{Vector{JuMP.Variable}}
    η::JuMP.Variable

    scenarios::Vector{Vector{Int}}
    upperbounds::Vector{Float64}
    lowerbounds::Vector{Float64}
    deployment::Vector{Vector{Int}}

    upptiming::Vector{Float64}
    lowtiming::Vector{Float64}
end
deployment(m::RobustDeployment) = m.deployment[end]

function Gamma(p::DeploymentProblem; α=params.α)
    demand = p.demand[p.train,:]
    γ_single = vec(maximum(demand,1) + 1*(maximum(demand,1) .== 0))
    γ_local = [quantile(Poisson(mean(sum(demand[:,vec(p.adjacency[i,:])], 2))),1-α) for i=1:p.nregions]
    γ_regional = [quantile(Poisson(mean(sum(demand[:,p.coverage[:,i]],2))),1-α) for i in 1:p.nlocations]
    γ_global = quantile(Poisson(mean(sum(demand,2))),1-α)
    Gamma(γ_single,γ_local,γ_regional,γ_global)
end

function Qrobust(problem::DeploymentProblem; α=params.α, verbose=false,
    solver=GurobiSolver(OutputFlag=0, MIPGapAbs=0.9)) #, TimeLimit=30))
    if verbose
        solver=GurobiSolver(OutputFlag=1) #, MIPGapAbs=0.9)
    end
    γ = Gamma(problem, α=α)
    upp_bound = maximum(γ._single)
    I = 1:problem.nlocations
    J = 1:problem.nregions

    m = JuMP.Model(solver=solver)
    JuMP.@variable(m, d[1:problem.nregions]>=0, Int)
    JuMP.@variable(m, p[1:problem.nregions], Bin)
    JuMP.@variable(m, q[1:problem.nlocations], Bin)

    for i in I, j in J
        problem.coverage[j, i] && JuMP.@constraint(m, p[j] <= q[i])
    end

    # Uncertainty
    for j in J
        JuMP.@constraint(m, d[j] <= γ._single[j]*p[j])
        adjacent_regions = filter(k->problem.adjacency[k,j],J)
        JuMP.@constraint(m, sum(d[k] for k in adjacent_regions) <= γ._local[j])
    end
    for i in I
        covered_regions = filter(j->problem.coverage[j,i],J)
        JuMP.@constraint(m, sum(d[j] for j in covered_regions) <= γ._regional[i])
    end
    JuMP.@constraint(m, sum(d[j] for j in J) <= γ._global)

    Qrobust(m, I, J, d, q, γ)
end

function evaluate{T <: Real}(Q::Qrobust, x::Vector{T})
    JuMP.@objective(Q.m, Max, sum(Q.d[j] for j in Q.J) - sum(x[i]*Q.q[i] for i in Q.I))
    status = JuMP.solve(Q.m)
    JuMP.getobjectivevalue(Q.m), Int[round(Int,d) for d in JuMP.getvalue(Q.d)]
end

function evaluate_objvalue{T <: Real}(Q::Qrobust, x::Vector{T})
    JuMP.@objective(Q.m, Max, sum(Q.d[j] for j in Q.J) - sum(x[i]*Q.q[i] for i in Q.I))
    status = JuMP.solve(Q.m)
    JuMP.getobjectivevalue(Q.m)
end

function RobustDeployment(p::DeploymentProblem; α=params.α, eps=params.ε, tol=params.δ,
    solver=GurobiSolver(OutputFlag=0, MIPGapAbs=0.9), verbose=false, master_verbose=false)
    if master_verbose
        solver=GurobiSolver(OutputFlag=1, MIPGapAbs=0.9)
    end
    I = 1:p.nlocations
    J = 1:p.nregions

    warmstart = naive_solution(p)

    m = JuMP.Model(solver=solver)
    JuMP.@variable(m, x[i=1:p.nlocations] >= 0, Int, start=warmstart[i])
    JuMP.@variable(m, η >= 0)
    y = Vector{Matrix{JuMP.Variable}}()
    z = Vector{JuMP.Variable}()

    # Initial Restricted Master Problem
    JuMP.@objective(m, Min, η)
    JuMP.@constraint(m, sum(x[i] for i=I) <= p.nambulances)
    for j in J # coverage over all regions
        JuMP.@constraint(m, sum(x[i] for i in filter(i->p.coverage[j,i], I)) >= 1)
    end

    RobustDeployment(m, Qrobust(p, α=α, verbose=verbose), I, J, x, y, z, η,
                     Vector{Vector{Int}}(), Vector{Float64}(), Vector{Float64}(),
                     Vector{Int}[warmstart], Vector{Float64}(), Vector{Float64}())
end

function add_scenario{T <: Real}(model::RobustDeployment, p::DeploymentProblem, scenario::Vector{T}; tol=params.δ)
    # Create variables yˡ
    push!(model.y, Array(JuMP.Variable, (p.nlocations,p.nregions)))
    l = length(model.y)
    for i in model.I, j in model.J
        model.y[l][i,j] = JuMP.Variable(model.m, 0, p.nambulances, :Int, String("y[$i,$j,$l]"))
    end
    push!(model.z, Array(JuMP.Variable, p.nregions))
    for j in model.J
        model.z[l][j] = JuMP.Variable(model.m, 0, Inf, :Int, String("z[$j,$l]"))
    end

    # (1) η >= 1ᵀ(dˡ + Bᴶyˡ)^+
    JuMP.@constraint(model.m, model.η >= sum(model.z[l][j] for j=model.J) + tol*sum(model.y[l][i,j] for i=model.I, j=model.J))
    for i in model.I # flow constraints at each station
        JuMP.@expression(m, outflow, sum(model.y[l][i,j] for j in filter(j->p.coverage[j,i], model.J)))
        JuMP.@constraint(model.m, model.x[i] >= outflow)
    end
    # (2) yˡ ∈ Y(x)
    for j in model.J # shortfall from satisfying demand/calls
        JuMP.@expression(m, inflow, sum(model.y[l][i,j] for i in filter(i->p.coverage[j,i], model.I)))
        JuMP.@constraint(model.m, model.z[l][j] >= scenario[j] - inflow)
    end
end

function solve(model::RobustDeployment, p::DeploymentProblem; verbose=false, maxiter=params.maxiter, eps=params.ε)
    LB = 0.0
    UB, scenario = evaluate(model.Q, model.deployment[end])
    push!(model.lowerbounds, LB)
    push!(model.upperbounds, UB)
    push!(model.scenarios, scenario)

    for k in 1:maxiter
        verbose && println("iteration $k: LB $LB, UB $UB")
        abs(UB - LB) < eps && break
        verbose && println("  solving Q with $(model.deployment[end])")

        add_scenario(model, p, scenario)
        tic()
        status = JuMP.solve(model.m)
        push!(model.lowtiming, toq())
        @assert status == :Optimal

        LB = JuMP.getobjectivevalue(model.m)

        push!(model.deployment, [round(Int,x) for x in JuMP.getvalue(model.x)])

        tic()
        shortfall, scenario = evaluate(model.Q, model.deployment[end])
        push!(model.upptiming, toq())
        UB = min(UB, shortfall)

        # for tracking convergence later
        push!(model.upperbounds, UB)
        push!(model.scenarios, scenario)
        push!(model.lowerbounds, LB)
    end
end


text:
#==============================================================================#
# CloudWatchEvents.jl
#
# This file is generated from:
# https://github.com/aws/aws-sdk-js/blob/master/apis/events-2015-10-07.normal.json
#==============================================================================#

__precompile__()

module CloudWatchEvents

using AWSCore

"""
    using AWSSDK.CloudWatchEvents.delete_rule
    delete_rule([::AWSConfig], arguments::Dict)
    delete_rule([::AWSConfig]; Name=)

    using AWSCore.Services.events
    events([::AWSConfig], "DeleteRule", arguments::Dict)
    events([::AWSConfig], "DeleteRule", Name=)

# DeleteRule Operation

Deletes the specified rule.

You must remove all targets from a rule using [RemoveTargets](@ref) before you can delete the rule.

When you delete a rule, incoming events might continue to match to the deleted rule. Please allow a short period of time for changes to take effect.

# Arguments

## `Name = ::String` -- *Required*
The name of the rule.

# Exceptions

`ConcurrentModificationException` or `InternalException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/events-2015-10-07/DeleteRule)
"""
@inline delete_rule(aws::AWSConfig=default_aws_config(); args...) = delete_rule(aws, args)

@inline delete_rule(aws::AWSConfig, args) = AWSCore.Services.events(aws, "DeleteRule", args)

@inline delete_rule(args) = delete_rule(default_aws_config(), args)

"""
    using AWSSDK.CloudWatchEvents.describe_event_bus
    describe_event_bus([::AWSConfig], arguments::Dict)
    describe_event_bus([::AWSConfig]; )

    using AWSCore.Services.events
    events([::AWSConfig], "DescribeEventBus", arguments::Dict)
    events([::AWSConfig], "DescribeEventBus", )

# DescribeEventBus Operation

Displays the external AWS accounts that are permitted to write events to your account using your account's event bus, and the associated policy. To enable your account to receive events from other accounts, use [PutPermission](@ref).

# Arguments

# Returns

`DescribeEventBusResponse`

# Exceptions

`ResourceNotFoundException` or `InternalException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/events-2015-10-07/DescribeEventBus)
"""
@inline describe_event_bus(aws::AWSConfig=default_aws_config(); args...) = describe_event_bus(aws, args)

@inline describe_event_bus(aws::AWSConfig, args) = AWSCore.Services.events(aws, "DescribeEventBus", args)

@inline describe_event_bus(args) = describe_event_bus(default_aws_config(), args)

"""
    using AWSSDK.CloudWatchEvents.describe_rule
    describe_rule([::AWSConfig], arguments::Dict)
    describe_rule([::AWSConfig]; Name=)

    using AWSCore.Services.events
    events([::AWSConfig], "DescribeRule", arguments::Dict)
    events([::AWSConfig], "DescribeRule", Name=)

# DescribeRule Operation

Describes the specified rule.

# Arguments

## `Name = ::String` -- *Required*
The name of the rule.

# Returns

`DescribeRuleResponse`

# Exceptions

`ResourceNotFoundException` or `InternalException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/events-2015-10-07/DescribeRule)
"""
@inline describe_rule(aws::AWSConfig=default_aws_config(); args...) = describe_rule(aws, args)

@inline describe_rule(aws::AWSConfig, args) = AWSCore.Services.events(aws, "DescribeRule", args)

@inline describe_rule(args) = describe_rule(default_aws_config(), args)

"""
    using AWSSDK.CloudWatchEvents.disable_rule
    disable_rule([::AWSConfig], arguments::Dict)
    disable_rule([::AWSConfig]; Name=)

    using AWSCore.Services.events
    events([::AWSConfig], "DisableRule", arguments::Dict)
    events([::AWSConfig], "DisableRule", Name=)

# DisableRule Operation

Disables the specified rule. A disabled rule won't match any events, and won't self-trigger if it has a schedule expression.

When you disable a rule, incoming events might continue to match to the disabled rule. Please allow a short period of time for changes to take effect.

# Arguments

## `Name = ::String` -- *Required*
The name of the rule.

# Exceptions

`ResourceNotFoundException`, `ConcurrentModificationException` or `InternalException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/events-2015-10-07/DisableRule)
"""
@inline disable_rule(aws::AWSConfig=default_aws_config(); args...) = disable_rule(aws, args)

@inline disable_rule(aws::AWSConfig, args) = AWSCore.Services.events(aws, "DisableRule", args)

@inline disable_rule(args) = disable_rule(default_aws_config(), args)

"""
    using AWSSDK.CloudWatchEvents.enable_rule
    enable_rule([::AWSConfig], arguments::Dict)
    enable_rule([::AWSConfig]; Name=)

    using AWSCore.Services.events
    events([::AWSConfig], "EnableRule", arguments::Dict)
    events([::AWSConfig], "EnableRule", Name=)

# EnableRule Operation

Enables the specified rule. If the rule does not exist, the operation fails.

When you enable a rule, incoming events might not immediately start matching to a newly enabled rule. Please allow a short period of time for changes to take effect.

# Arguments

## `Name = ::String` -- *Required*
The name of the rule.

# Exceptions

`ResourceNotFoundException`, `ConcurrentModificationException` or `InternalException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/events-2015-10-07/EnableRule)
"""
@inline enable_rule(aws::AWSConfig=default_aws_config(); args...) = enable_rule(aws, args)

@inline enable_rule(aws::AWSConfig, args) = AWSCore.Services.events(aws, "EnableRule", args)

@inline enable_rule(args) = enable_rule(default_aws_config(), args)

"""
    using AWSSDK.CloudWatchEvents.list_rule_names_by_target
    list_rule_names_by_target([::AWSConfig], arguments::Dict)
    list_rule_names_by_target([::AWSConfig]; TargetArn=, <keyword arguments>)

    using AWSCore.Services.events
    events([::AWSConfig], "ListRuleNamesByTarget", arguments::Dict)
    events([::AWSConfig], "ListRuleNamesByTarget", TargetArn=, <keyword arguments>)

# ListRuleNamesByTarget Operation

Lists the rules for the specified target. You can see which of the rules in Amazon CloudWatch Events can invoke a specific target in your account.

# Arguments

## `TargetArn = ::String` -- *Required*
The Amazon Resource Name (ARN) of the target resource.

## `NextToken = ::String`
The token returned by a previous call to retrieve the next set of results.

## `Limit = ::Int`
The maximum number of results to return.

# Returns

`ListRuleNamesByTargetResponse`

# Exceptions

`InternalException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/events-2015-10-07/ListRuleNamesByTarget)
"""
@inline list_rule_names_by_target(aws::AWSConfig=default_aws_config(); args...) = list_rule_names_by_target(aws, args)

@inline list_rule_names_by_target(aws::AWSConfig, args) = AWSCore.Services.events(aws, "ListRuleNamesByTarget", args)

@inline list_rule_names_by_target(args) = list_rule_names_by_target(default_aws_config(), args)

"""
    using AWSSDK.CloudWatchEvents.list_rules
    list_rules([::AWSConfig], arguments::Dict)
    list_rules([::AWSConfig]; <keyword arguments>)

    using AWSCore.Services.events
    events([::AWSConfig], "ListRules", arguments::Dict)
    events([::AWSConfig], "ListRules", <keyword arguments>)

# ListRules Operation

Lists your Amazon CloudWatch Events rules. You can either list all the rules or you can provide a prefix to match to the rule names.

# Arguments

## `NamePrefix = ::String`
The prefix matching the rule name.

## `NextToken = ::String`
The token returned by a previous call to retrieve the next set of results.

## `Limit = ::Int`
The maximum number of results to return.

# Returns

`ListRulesResponse`

# Exceptions

`InternalException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/events-2015-10-07/ListRules)
"""
@inline list_rules(aws::AWSConfig=default_aws_config(); args...) = list_rules(aws, args)

@inline list_rules(aws::AWSConfig, args) = AWSCore.Services.events(aws, "ListRules", args)

@inline list_rules(args) = list_rules(default_aws_config(), args)

"""
    using AWSSDK.CloudWatchEvents.list_targets_by_rule
    list_targets_by_rule([::AWSConfig], arguments::Dict)
    list_targets_by_rule([::AWSConfig]; Rule=, <keyword arguments>)

    using AWSCore.Services.events
    events([::AWSConfig], "ListTargetsByRule", arguments::Dict)
    events([::AWSConfig], "ListTargetsByRule", Rule=, <keyword arguments>)

# ListTargetsByRule Operation

Lists the targets assigned to the specified rule.

# Arguments

## `Rule = ::String` -- *Required*
The name of the rule.

## `NextToken = ::String`
The token returned by a previous call to retrieve the next set of results.

## `Limit = ::Int`
The maximum number of results to return.

# Returns

`ListTargetsByRuleResponse`

# Exceptions

`ResourceNotFoundException` or `InternalException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/events-2015-10-07/ListTargetsByRule)
"""
@inline list_targets_by_rule(aws::AWSConfig=default_aws_config(); args...) = list_targets_by_rule(aws, args)

@inline list_targets_by_rule(aws::AWSConfig, args) = AWSCore.Services.events(aws, "ListTargetsByRule", args)

@inline list_targets_by_rule(args) = list_targets_by_rule(default_aws_config(), args)

"""
    using AWSSDK.CloudWatchEvents.put_events
    put_events([::AWSConfig], arguments::Dict)
    put_events([::AWSConfig]; Entries=)

    using AWSCore.Services.events
    events([::AWSConfig], "PutEvents", arguments::Dict)
    events([::AWSConfig], "PutEvents", Entries=)

# PutEvents Operation

Sends custom events to Amazon CloudWatch Events so that they can be matched to rules.

# Arguments

## `Entries = [[ ... ], ...]` -- *Required*
The entry that defines an event in your system. You can specify several parameters for the entry such as the source and type of the event, resources associated with the event, and so on.
```
 Entries = [[
        "Time" =>  timestamp,
        "Source" =>  ::String,
        "Resources" =>  [::String, ...],
        "DetailType" =>  ::String,
        "Detail" =>  ::String
    ], ...]
```

# Returns

`PutEventsResponse`

# Exceptions

`InternalException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/events-2015-10-07/PutEvents)
"""
@inline put_events(aws::AWSConfig=default_aws_config(); args...) = put_events(aws, args)

@inline put_events(aws::AWSConfig, args) = AWSCore.Services.events(aws, "PutEvents", args)

@inline put_events(args) = put_events(default_aws_config(), args)

"""
    using AWSSDK.CloudWatchEvents.put_permission
    put_permission([::AWSConfig], arguments::Dict)
    put_permission([::AWSConfig]; Action=, Principal=, StatementId=)

    using AWSCore.Services.events
    events([::AWSConfig], "PutPermission", arguments::Dict)
    events([::AWSConfig], "PutPermission", Action=, Principal=, StatementId=)

# PutPermission Operation

Running `PutPermission` permits the specified AWS account to put events to your account's default *event bus*. CloudWatch Events rules in your account are triggered by these events arriving to your default event bus.

For another account to send events to your account, that external account must have a CloudWatch Events rule with your account's default event bus as a target.

To enable multiple AWS accounts to put events to your default event bus, run `PutPermission` once for each of these accounts.

The permission policy on the default event bus cannot exceed 10KB in size.

# Arguments

## `Action = ::String` -- *Required*
The action that you are enabling the other account to perform. Currently, this must be `events:PutEvents`.

## `Principal = ::String` -- *Required*
The 12-digit AWS account ID that you are permitting to put events to your default event bus. Specify "*" to permit any account to put events to your default event bus.

If you specify "*", avoid creating rules that may match undesirable events. To create more secure rules, make sure that the event pattern for each rule contains an `account` field with a specific account ID from which to receive events. Rules with an account field do not match any events sent from other accounts.

## `StatementId = ::String` -- *Required*
An identifier string for the external account that you are granting permissions to. If you later want to revoke the permission for this external account, specify this `StatementId` when you run [RemovePermission](@ref).

# Exceptions

`ResourceNotFoundException`, `PolicyLengthExceededException`, `InternalException` or `ConcurrentModificationException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/events-2015-10-07/PutPermission)
"""
@inline put_permission(aws::AWSConfig=default_aws_config(); args...) = put_permission(aws, args)

@inline put_permission(aws::AWSConfig, args) = AWSCore.Services.events(aws, "PutPermission", args)

@inline put_permission(args) = put_permission(default_aws_config(), args)

"""
    using AWSSDK.CloudWatchEvents.put_rule
    put_rule([::AWSConfig], arguments::Dict)
    put_rule([::AWSConfig]; Name=, <keyword arguments>)

    using AWSCore.Services.events
    events([::AWSConfig], "PutRule", arguments::Dict)
    events([::AWSConfig], "PutRule", Name=, <keyword arguments>)

# PutRule Operation

Creates or updates the specified rule. Rules are enabled by default, or based on value of the state. You can disable a rule using [DisableRule](@ref).

If you are updating an existing rule, the rule is completely replaced with what you specify in this `PutRule` command. If you omit arguments in `PutRule`, the old values for those arguments are not kept. Instead, they are replaced with null values.

When you create or update a rule, incoming events might not immediately start matching to new or updated rules. Please allow a short period of time for changes to take effect.

A rule must contain at least an EventPattern or ScheduleExpression. Rules with EventPatterns are triggered when a matching event is observed. Rules with ScheduleExpressions self-trigger based on the given schedule. A rule can have both an EventPattern and a ScheduleExpression, in which case the rule triggers on matching events as well as on a schedule.

Most services in AWS treat : or / as the same character in Amazon Resource Names (ARNs). However, CloudWatch Events uses an exact match in event patterns and rules. Be sure to use the correct ARN characters when creating event patterns so that they match the ARN syntax in the event you want to match.

# Arguments

## `Name = ::String` -- *Required*
The name of the rule that you are creating or updating.

## `ScheduleExpression = ::String`
The scheduling expression. For example, "cron(0 20 * * ? *)" or "rate(5 minutes)".

## `EventPattern = ::String`
The event pattern. For more information, see [Events and Event Patterns](http://docs.aws.amazon.com/AmazonCloudWatch/latest/events/CloudWatchEventsandEventPatterns.html) in the *Amazon CloudWatch Events User Guide*.

## `State = "ENABLED" or "DISABLED"`
Indicates whether the rule is enabled or disabled.

## `Description = ::String`
A description of the rule.

## `RoleArn = ::String`
The Amazon Resource Name (ARN) of the IAM role associated with the rule.

# Returns

`PutRuleResponse`

# Exceptions

`InvalidEventPatternException`, `LimitExceededException`, `ConcurrentModificationException` or `InternalException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/events-2015-10-07/PutRule)
"""
@inline put_rule(aws::AWSConfig=default_aws_config(); args...) = put_rule(aws, args)

@inline put_rule(aws::AWSConfig, args) = AWSCore.Services.events(aws, "PutRule", args)

@inline put_rule(args) = put_rule(default_aws_config(), args)

"""
    using AWSSDK.CloudWatchEvents.put_targets
    put_targets([::AWSConfig], arguments::Dict)
    put_targets([::AWSConfig]; Rule=, Targets=)

    using AWSCore.Services.events
    events([::AWSConfig], "PutTargets", arguments::Dict)
    events([::AWSConfig], "PutTargets", Rule=, Targets=)

# PutTargets Operation

Adds the specified targets to the specified rule, or updates the targets if they are already associated with the rule.

Targets are the resources that are invoked when a rule is triggered.

You can configure the following as targets for CloudWatch Events:

*   EC2 instances

*   AWS Lambda functions

*   Streams in Amazon Kinesis Streams

*   Delivery streams in Amazon Kinesis Firehose

*   Amazon ECS tasks

*   AWS Step Functions state machines

*   AWS Batch jobs

*   Pipelines in Amazon Code Pipeline

*   Amazon Inspector assessment templates

*   Amazon SNS topics

*   Amazon SQS queues, including FIFO queues

*   The default event bus of another AWS account

Note that creating rules with built-in targets is supported only in the AWS Management Console.

For some target types, `PutTargets` provides target-specific parameters. If the target is an Amazon Kinesis stream, you can optionally specify which shard the event goes to by using the `KinesisParameters` argument. To invoke a command on multiple EC2 instances with one rule, you can use the `RunCommandParameters` field.

To be able to make API calls against the resources that you own, Amazon CloudWatch Events needs the appropriate permissions. For AWS Lambda and Amazon SNS resources, CloudWatch Events relies on resource-based policies. For EC2 instances, Amazon Kinesis streams, and AWS Step Functions state machines, CloudWatch Events relies on IAM roles that you specify in the `RoleARN` argument in `PutTargets`. For more information, see [Authentication and Access Control](http://docs.aws.amazon.com/AmazonCloudWatch/latest/events/auth-and-access-control-cwe.html) in the *Amazon CloudWatch Events User Guide*.

If another AWS account is in the same region and has granted you permission (using `PutPermission`), you can send events to that account by setting that account's event bus as a target of the rules in your account. To send the matched events to the other account, specify that account's event bus as the `Arn` when you run `PutTargets`. If your account sends events to another account, your account is charged for each sent event. Each event sent to antoher account is charged as a custom event. The account receiving the event is not charged. For more information on pricing, see [Amazon CloudWatch Pricing](https://aws.amazon.com/cloudwatch/pricing/).

For more information about enabling cross-account events, see [PutPermission](@ref).

**Input**, **InputPath** and **InputTransformer** are mutually exclusive and optional parameters of a target. When a rule is triggered due to a matched event:

*   If none of the following arguments are specified for a target, then the entire event is passed to the target in JSON form (unless the target is Amazon EC2 Run Command or Amazon ECS task, in which case nothing from the event is passed to the target).

*   If **Input** is specified in the form of valid JSON, then the matched event is overridden with this constant.

*   If **InputPath** is specified in the form of JSONPath (for example, `\$.detail`), then only the part of the event specified in the path is passed to the target (for example, only the detail part of the event is passed).

*   If **InputTransformer** is specified, then one or more specified JSONPaths are extracted from the event and used as values in a template that you specify as the input to the target.

When you specify `InputPath` or `InputTransformer`, you must use JSON dot notation, not bracket notation.

When you add targets to a rule and the associated rule triggers soon after, new or updated targets might not be immediately invoked. Please allow a short period of time for changes to take effect.

This action can partially fail if too many requests are made at the same time. If that happens, `FailedEntryCount` is non-zero in the response and each entry in `FailedEntries` provides the ID of the failed target and the error code.

# Arguments

## `Rule = ::String` -- *Required*
The name of the rule.

## `Targets = [[ ... ], ...]` -- *Required*
The targets to update or add to the rule.
```
 Targets = [[
        "Id" => <required> ::String,
        "Arn" => <required> ::String,
        "RoleArn" =>  ::String,
        "Input" =>  ::String,
        "InputPath" =>  ::String,
        "InputTransformer" =>  [
            "InputPathsMap" =>  ::Dict{String,String},
            "InputTemplate" => <required> ::String
        ],
        "KinesisParameters" =>  ["PartitionKeyPath" => <required> ::String],
        "RunCommandParameters" =>  ["RunCommandTargets" => <required> [[
                "Key" => <required> ::String,
                "Values" => <required> [::String, ...]
            ], ...]],
        "EcsParameters" =>  [
            "TaskDefinitionArn" => <required> ::String,
            "TaskCount" =>  ::Int
        ],
        "BatchParameters" =>  [
            "JobDefinition" => <required> ::String,
            "JobName" => <required> ::String,
            "ArrayProperties" =>  ["Size" =>  ::Int],
            "RetryStrategy" =>  ["Attempts" =>  ::Int]
        ],
        "SqsParameters" =>  ["MessageGroupId" =>  ::String]
    ], ...]
```

# Returns

`PutTargetsResponse`

# Exceptions

`ResourceNotFoundException`, `ConcurrentModificationException`, `LimitExceededException` or `InternalException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/events-2015-10-07/PutTargets)
"""
@inline put_targets(aws::AWSConfig=default_aws_config(); args...) = put_targets(aws, args)

@inline put_targets(aws::AWSConfig, args) = AWSCore.Services.events(aws, "PutTargets", args)

@inline put_targets(args) = put_targets(default_aws_config(), args)

"""
    using AWSSDK.CloudWatchEvents.remove_permission
    remove_permission([::AWSConfig], arguments::Dict)
    remove_permission([::AWSConfig]; StatementId=)

    using AWSCore.Services.events
    events([::AWSConfig], "RemovePermission", arguments::Dict)
    events([::AWSConfig], "RemovePermission", StatementId=)

# RemovePermission Operation

Revokes the permission of another AWS account to be able to put events to your default event bus. Specify the account to revoke by the `StatementId` value that you associated with the account when you granted it permission with `PutPermission`. You can find the `StatementId` by using [DescribeEventBus](@ref).

# Arguments

## `StatementId = ::String` -- *Required*
The statement ID corresponding to the account that is no longer allowed to put events to the default event bus.

# Exceptions

`ResourceNotFoundException`, `InternalException` or `ConcurrentModificationException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/events-2015-10-07/RemovePermission)
"""
@inline remove_permission(aws::AWSConfig=default_aws_config(); args...) = remove_permission(aws, args)

@inline remove_permission(aws::AWSConfig, args) = AWSCore.Services.events(aws, "RemovePermission", args)

@inline remove_permission(args) = remove_permission(default_aws_config(), args)

"""
    using AWSSDK.CloudWatchEvents.remove_targets
    remove_targets([::AWSConfig], arguments::Dict)
    remove_targets([::AWSConfig]; Rule=, Ids=)

    using AWSCore.Services.events
    events([::AWSConfig], "RemoveTargets", arguments::Dict)
    events([::AWSConfig], "RemoveTargets", Rule=, Ids=)

# RemoveTargets Operation

Removes the specified targets from the specified rule. When the rule is triggered, those targets are no longer be invoked.

When you remove a target, when the associated rule triggers, removed targets might continue to be invoked. Please allow a short period of time for changes to take effect.

This action can partially fail if too many requests are made at the same time. If that happens, `FailedEntryCount` is non-zero in the response and each entry in `FailedEntries` provides the ID of the failed target and the error code.

# Arguments

## `Rule = ::String` -- *Required*
The name of the rule.

## `Ids = [::String, ...]` -- *Required*
The IDs of the targets to remove from the rule.

# Returns

`RemoveTargetsResponse`

# Exceptions

`ResourceNotFoundException`, `ConcurrentModificationException` or `InternalException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/events-2015-10-07/RemoveTargets)
"""
@inline remove_targets(aws::AWSConfig=default_aws_config(); args...) = remove_targets(aws, args)

@inline remove_targets(aws::AWSConfig, args) = AWSCore.Services.events(aws, "RemoveTargets", args)

@inline remove_targets(args) = remove_targets(default_aws_config(), args)

"""
    using AWSSDK.CloudWatchEvents.test_event_pattern
    test_event_pattern([::AWSConfig], arguments::Dict)
    test_event_pattern([::AWSConfig]; EventPattern=, Event=)

    using AWSCore.Services.events
    events([::AWSConfig], "TestEventPattern", arguments::Dict)
    events([::AWSConfig], "TestEventPattern", EventPattern=, Event=)

# TestEventPattern Operation

Tests whether the specified event pattern matches the provided event.

Most services in AWS treat : or / as the same character in Amazon Resource Names (ARNs). However, CloudWatch Events uses an exact match in event patterns and rules. Be sure to use the correct ARN characters when creating event patterns so that they match the ARN syntax in the event you want to match.

# Arguments

## `EventPattern = ::String` -- *Required*
The event pattern. For more information, see [Events and Event Patterns](http://docs.aws.amazon.com/AmazonCloudWatch/latest/events/CloudWatchEventsandEventPatterns.html) in the *Amazon CloudWatch Events User Guide*.

## `Event = ::String` -- *Required*
The event, in JSON format, to test against the event pattern.

# Returns

`TestEventPatternResponse`

# Exceptions

`InvalidEventPatternException` or `InternalException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/events-2015-10-07/TestEventPattern)
"""
@inline test_event_pattern(aws::AWSConfig=default_aws_config(); args...) = test_event_pattern(aws, args)

@inline test_event_pattern(aws::AWSConfig, args) = AWSCore.Services.events(aws, "TestEventPattern", args)

@inline test_event_pattern(args) = test_event_pattern(default_aws_config(), args)

end # module CloudWatchEvents

#==============================================================================#
# End of file
#==============================================================================#


text:
module Test

export @test, @test_fails, @test_throws, @test_approx_eq, @test_approx_eq_eps, @inferred

abstract Result
type Success <: Result
    expr
end
type Failure <: Result
    expr
end
type Error <: Result
    expr
    err
    backtrace
end

default_handler(r::Success) = nothing
default_handler(r::Failure) = error("test failed: $(r.expr)")
default_handler(r::Error)   = rethrow(r)

handler() = get(task_local_storage(), :TEST_HANDLER, default_handler)

with_handler(f::Function, handler) =
    task_local_storage(f, :TEST_HANDLER, handler)

import Base.showerror

showerror(io::IO, r::Error) = showerror(io, r, {})
function showerror(io::IO, r::Error, bt)
    println(io, "test error during $(r.expr)")
    showerror(io, r.err, r.backtrace)
end

function do_test(body,qex)
    handler()(try
        body() ? Success(qex) : Failure(qex)
    catch err
        Error(qex,err,catch_backtrace())
    end)
end

function do_test_throws(body, qex, bt, extype)
    handler()(try
        body()
        Failure(qex)
    catch err
        if extype == nothing
            Base.warn("@test_throws without an exception type is deprecated\n\t Use @test_throws $(typeof(err)) $(qex)", bt = bt)
            Success(qex)
        else
            if isa(err, extype)
                Success(qex)
            else
                Failure(qex)
            end
        end
    end)
end

macro test(ex)
    :(do_test(()->($(esc(ex))),$(Expr(:quote,ex))))
end

macro test_throws(args...)
    ex = nothing
    extype = nothing
    # Users should pass (ExceptionType, Expression) but we give a warning to users that only pass (Expression)
    if length(args) == 1
        ex = args[1]
    elseif length(args) == 2
        ex = args[2]
        extype = args[1]
    end
    :(do_test_throws(()->($(esc(ex))),$(Expr(:quote,ex)),backtrace(),$(esc(extype))))
end

macro test_fails(ex)
    Base.warn_once("@test_fails is deprecated, use @test_throws instead.")
    :(@test_throws $ex Exception)
end

approx_full(x::AbstractArray) = x
approx_full(x::Number) = x
approx_full(x) = full(x)

function test_approx_eq(va, vb, Eps, astr, bstr)
    va = approx_full(va)
    vb = approx_full(vb)
    if length(va) != length(vb)
        error("lengths of ", astr, " and ", bstr, " do not match: ",
              "\n  ", astr, " (length $(length(va))) = ", va,
              "\n  ", bstr, " (length $(length(vb))) = ", vb)
    end
    diff = real(zero(eltype(va)))
    for i = 1:length(va)
        xa = va[i]; xb = vb[i]
        if isfinite(xa) && isfinite(xb)
            diff = max(diff, abs(xa-xb))
        elseif !isequal(xa,xb)
            error("mismatch of non-finite elements: ",
                  "\n  ", astr, " = ", va,
                  "\n  ", bstr, " = ", vb)
        end
    end

    if !isnan(Eps) && !(diff <= Eps)
        sdiff = string("|", astr, " - ", bstr, "| <= ", Eps)
        error("assertion failed: ", sdiff,
              "\n  ", astr, " = ", va,
              "\n  ", bstr, " = ", vb,
              "\n  difference = ", diff, " > ", Eps)
    end
end

array_eps{T}(a::AbstractArray{Complex{T}}) = eps(float(maximum(x->(isfinite(x) ? abs(x) : nan(T)), a)))
array_eps(a) = eps(float(maximum(x->(isfinite(x) ? abs(x) : nan(x)), a)))

test_approx_eq(va, vb, astr, bstr) =
    test_approx_eq(va, vb, 1E4*length(va)*max(array_eps(va), array_eps(vb)), astr, bstr)

macro test_approx_eq_eps(a, b, c)
    :(test_approx_eq($(esc(a)), $(esc(b)), $(esc(c)), $(string(a)), $(string(b))))
end

macro test_approx_eq(a, b)
    :(test_approx_eq($(esc(a)), $(esc(b)), $(string(a)), $(string(b))))
end

macro inferred(ex)
    ex.head == :call || error("@inferred requires a call expression")
    quote
        vals = ($([esc(ex.args[i]) for i = 2:length(ex.args)]...),)
        inftypes = Base.return_types($(ex.args[1]), ($([:(typeof(vals[$i])) for i = 1:length(ex.args)-1]...),))
        @assert length(inftypes) == 1
        result = $(esc(ex.args[1]))(vals...)
        rettype = typeof(result)
        is(rettype, inftypes[1]) || error("return type $rettype does not match inferred return type $(inftypes[1])")
        result
    end
end

end # module

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
"""Calculate the R integral from Popel divided by x, (x = q*L0, L0 = zero-temp London depth).

Parameters
----------
a : float

b : float

x : float

Returns
-------
r : float
    The R integral from Popel, divided by x.

Note
----
See R. Pöpel (1989), doi: 10.1063/1.343622 for more details."""
function intR(a, b, x)

    z2 = a^2+b^2

    #for small x
    if x < 0.01*sqrt(z2)
        r = b/(3.0*z2) #This is really r/x

    #for large x
    elseif x > 100*sqrt(z2)
        r = (π*(1+(b^2-a^2)/x^2)/4-b/x)/x #This is really r/x

    #in between x
    else
        #calculate all the terms of r
        r = (1/x^2)*(-0.5*b*x+0.25*a*b*log((z2+x^2+2*a*x)/(z2+x^2-2*a*x))+
        0.25*(x^2+b^2-a^2)*atan2(2*b*x,(z2-x^2)))/x #This is really r/x
    end

    return r
end

"""Calculate the R integral from Popel divided by x, (x = q*L0, L0 = zero-temp London depth).

Parameters
----------
a : float

b : float

x : float

Returns
-------
s : float
    The S integral from Popel, divided by x.

Note
----
See R. Pöpel (1989), doi: 10.1063/1.343622 for more details."""
function intS(a, b, x)

    z2 = a^2+b^2

    #for small x
    if x < 0.01*sqrt(z2)
        s = a/(3.0*z2) #This is really s/x

    #for large x
    elseif x > 100*sqrt(z2)
        s = (a/x - a*b*π/(2*x^2))/x #This is really s/x

    #in between x
    else
        #calculate all the terms of s
        s = (1/x^2)*(0.5*(a*x)+
            0.125*(x^2+b^2-a^2)*log((z2+x^2+2*a*x)/(z2+x^2-2*a*x))-
            0.5*b*a*atan2(2*b*x,(z2-x^2)))/x #This is really s/x
    end

    return s
end

"""Similar to intS, but divided by a. This is specifically useful for when
a=0, as lim(a->0) intS/a should not blow up!"""
function intS_over_a(b, x)

    z2 = b^2

    #for small x
    if x < 0.01*sqrt(z2)
        s = 1.0/(3.0*z2) #This is really s/(x*a)

    #for large x
    elseif x > 100*sqrt(z2)
        s = (1.0/x - b*π/(2*x^2))/x #This is really s/(x*a)

    #in between x
    else
        #calculate all the terms of s
        s = (1/x^2)*(x-0.5*b*atan2(2*b*x,(z2-x^2)))/x #This is really s/(x*a)
    end

    return s
end


text:
"""
    LdagLSparseSuperopProblem <: AbstractProblem

Problem or finding the steady state of a ℒdagℒ matrix
"""
struct LdagLSparseSuperopProblem{B, SM} <: HermitianMatrixProblem where {B<:Basis,
                                                 SM<:SparseMatrixCSC}
    HilbSpace::B            # 0
    LdagL::SM
    ρss
end

"""
    LdagLSparseSuperopProblem([T=STD_REAL_PREC], lindbladian)

Creates a problem for minimizing the cost function 𝒞 = ∑|ρ(σ)|²|⟨⟨σ|ℒ'ℒ |ρ⟩⟩|².
Computes |⟨⟨σ|ℒ'ℒ |ρ⟩⟩| by building the sparse superoperator, which can be done
for sizes up to dimℋ < 500. For more than 9 spins it is reccomended to use the
command LRhoSparseOpProblem

`lindbladian` can either be the lindbladian on a graph, a QuantumOptics superoperator
or the Hamiltonian and a vector of collapse operators.

`T=STD_REAL_PREC` by default is the numerical precision used. It should match that of
the network.
"""
LdagLSparseSuperopProblem(args...) = LdagLSparseSuperopProblem(STD_REAL_PREC, args...)
LdagLSparseSuperopProblem(T::Type{<:Number}, gl::GraphLindbladian) =
    LdagLSparseSuperopProblem(T, liouvillian(gl))
LdagLSparseSuperopProblem(T::Type{<:Number}, Ham::DataOperator, cops::Vector) =
    LdagLSparseSuperopProblem(T, liouvillian(Ham, cops))
LdagLSparseSuperopProblem(T::Type{<:Number}, Liouv::SparseSuperOperator) =
    LdagLSparseSuperopProblem(first(Liouv.basis_l), Liouv.data'*Liouv.data, 0.0)

basis(prob::LdagLSparseSuperopProblem) = prob.HilbSpace

function compute_Cloc(prob::LdagLSparseSuperopProblem, net::MatrixNet, σ, lnψ=net(σ), σp=deepcopy(σ))
    ℒdagℒ = prob.LdagL
    set_index!(σp, index(σ))

    #### Now compute E(S) = Σₛ⟨s|Hψ⟩/⟨s|ψ⟩
    C_loc = zero(Complex{real(out_type(net))})
    # Iterate through all elements in row i_σ of the matrix computing
    # ⟨i_σ|ℒdagℒψ⟩ = Σ_{i_σp} ⟨i_σ|ℒdagℒ|i_σp⟩⟨i_σp|ψ⟩
    # NOTE: ℒdagℒ is CSC, but I would like a CSR matrix. Since it is hermitian I
    # can simply take the conjugate of the elements in the columns
    i_σ = index(σ)
    for row_id = ℒdagℒ.colptr[i_σ]:(ℒdagℒ.colptr[i_σ+1]-1)
      # Find nonzero elements s by doing <i_sp|ℒdagℒ|i_σ>
      i_σp = ℒdagℒ.rowval[row_id]
      # BackConvert to int
      set_index!(σp, i_σp)
      # Compute the log(ψ(σ)/ψ(σ')), by only computing differences.
      log_ratio = logψ(net, σp) - lnψ
      # Conj because I am taking the transpose... see the note above.
      C_loc += conj(ℒdagℒ.nzval[row_id]) * exp(log_ratio)
    end
    C_loc
end

Base.show(io::IO, p::LdagLSparseSuperopProblem) = print(io,
    "LdagLSparseSuperopProblem on space : $(basis(p)), computing the energy of LdagL with the sparse liouvillian")


text:
models = (
    SVMLinearClassifier,
    SVMClassifier,
    SVMNuClassifier,
)

fparams = (
    SVMLinearClassifier=(:coef, :intercept, :classes),
    SVMClassifier=(:support, :support_vectors, :n_support, :dual_coef, :coef,
                   :intercept, :fit_status, :classes),
    SVMNuClassifier=(:support, :support_vectors, :n_support, :dual_coef, :coef,
                     :intercept, :fit_status, :classes),
)

@testset "Fit/Predict" begin
    X, y = simple_binaryclf()
    for mod in models
        m = mod()
        acc, f = test_clf(m, X, y)
        @test acc >= 0.75
        @test keys(MB.fitted_params(m, f)) == getproperty(fparams, Symbol(mod))
    end
end

@testset "Docstrings" begin
    for mod in models
        m = mod()
        @test !isempty(MB.docstring(m))
    end
end

models = (
    SVMLinearRegressor,
    SVMRegressor,
    SVMNuRegressor,
)

fparams = (
    SVMLinearRegressor=(:coef, :intercept),
    SVMRegressor=(:support, :support_vectors, :dual_coef, :coef,
                  :intercept, :fit_status),
    SVMNuRegressor=(:support, :support_vectors, :dual_coef, :coef,
                    :intercept),
)

@testset "Fit/Predict" begin
    X, y, ls = simple_regression()
    for mod in models
        m = mod()
        r, f = test_regression(m, X, y, ls)
        @test r < 5
        @test keys(MB.fitted_params(m, f)) == getproperty(fparams, Symbol(mod))
    end
end