text:
"""
Stores the results from running the benchmarks on a package.

The following (unexported) methods are defined on a `BenchmarkResults` (written below as `results`):

* `name(results)::String` - The commit of the package benchmarked
* `commit(results)::String` - The commit of the package benchmarked. If the package repository was dirty, the string `"dirty"` is returned.
* `juliacommit(results)::String` - The commit of the Julia executable that ran the benchmarks
* `benchmarkgroup(results)::BenchmarkGroup` - a [`BenchmarkGroup`](https://github.com/JuliaCI/BenchmarkTools.jl/blob/master/doc/manual.md#the-benchmarkgroup-type)
   contaning the results of the benchmark.
* `date(results)::DateTime` - Tthe time when the benchmarks were executed
* `benchmarkconfig(results)::BenchmarkConfig` - The [`BenchmarkConfig`](@ref) used for the benchmarks.

`BenchmarkResults` can be exported to markdown using the function [`export_markdown`](@ref).
"""
struct BenchmarkResults
    name::String
    commit::String
    benchmarkgroup::BenchmarkGroup
    date::DateTime
    julia_commit::String
    vinfo::String
    benchmarkconfig::BenchmarkConfig
end

name(results::BenchmarkResults) = results.name
commit(results::BenchmarkResults) = results.commit
juliacommit(results::BenchmarkResults) = results.julia_commit
benchmarkgroup(results::BenchmarkResults) = results.benchmarkgroup
date(results::BenchmarkResults) = results.date
benchmarkconfig(results::BenchmarkResults) = results.benchmarkconfig
Base.versioninfo(results::BenchmarkResults) = results.vinfo

function Base.show(io::IO, results::BenchmarkResults)
    print(io, "Benchmarkresults:\n")
    println(io, "    Package: ", results.name)
    println(io, "    Date: ", Base.Dates.format(results.date, "d u Y - HH:MM"))
    println(io, "    Package commit: ", results.commit[1:min(length(results.commit), 6)])
    println(io, "    Julia commit: ", results.julia_commit[1:6])
    iob = IOBuffer()
    ioc = IOContext(iob)
    show(ioc, MIME("text/plain"), results.benchmarkgroup)
    println(io,   "    BenchmarkGroup:")
    print(join("        " .* split(String(take!(iob)), "\n"), "\n"))
end

"""
    writeresults(file::String, results::BenchmarkResults)

Writes the [`BenchmarkResults`](@ref) to `file`.
"""
function writeresults(file::String, results::BenchmarkResults)
    open(file, "w") do io
        JSON.print(io,
            Dict(
                "name" => results.name,
                "commit" => results.commit,
                "benchmarkgroup" => sprint(BenchmarkTools.save, results.benchmarkgroup),
                "date" => results.date,
                "julia_commit" => results.julia_commit,
                "vinfo" => results.vinfo,
                "benchmarkconfig" => results.benchmarkconfig
            )
        )
    end
end

"""
    readresults(file::String)

Reads the [`BenchmarkResults`](@ref) stored in `file` (given as a path).
"""
function readresults(file::String)
    d = JSON.parsefile(file)
    BenchmarkResults(
        d["name"],
        d["commit"],
        BenchmarkTools.load(IOBuffer(d["benchmarkgroup"]))[1],
        DateTime(d["date"]),
        d["julia_commit"],
        d["vinfo"],
        BenchmarkConfig(d["benchmarkconfig"]),
    )
end

"""
    export_markdown(file::String, results::Union{BenchmarkResults, BenchmarkJudgement})
    export_markdown(io::IO,       results::Union{BenchmarkResults, BenchmarkJudgement})

Writes the `results` to `file` or `io` in markdown format.

See also: [`BenchmarkResults`](@ref), [`BenchmarkJudgement`](@ref)
"""
function export_markdown(file::String, results::BenchmarkResults)
    open(file, "w") do f
        export_markdown(f, results)
    end
end

function export_markdown(io::IO, results::BenchmarkResults)
    env_str = if isempty(benchmarkconfig(results).env)
        "None"
    else
        join(String[string("`", k, " => ", v, "`") for (k, v) in benchmarkconfig(results).env], " ")
    end

    jlcmd = benchmarkconfig(results).juliacmd
    flags = length(jlcmd) <= 1 ? [] : jlcmd[2:end]
    julia_command_flags = if isempty(flags)
        "None"
    else
        """`$(join(flags, ","))`"""
    end

    println(io, """
                # Benchmark Report for *$(name(results))*

                ## Job Properties
                * Time of benchmark: $(Base.Dates.format(date(results), "d u Y - H:M"))
                * Package commit: $(commit(results)[1:min(6, length(commit(results)))])
                * Julia commit: $(juliacommit(results)[1:min(6, length(juliacommit(results)))])
                * Julia command flags: $julia_command_flags
                * Environment variables: $env_str
                """)

      println(io, """
                ## Results
                Below is a table of this job's results, obtained by running the benchmarks.
                The values listed in the `ID` column have the structure `[parent_group, child_group, ..., key]`, and can be used to
                index into the BaseBenchmarks suite to retrieve the corresponding benchmarks.
                The percentages accompanying time and memory values in the below table are noise tolerances. The "true"
                time/memory value for a given benchmark is expected to fall within this percentage of the reported value.
                An empty cell means that the value was zero.
                """)

    entries = BenchmarkTools.leaves(benchmarkgroup(results))
    entries = entries[sortperm(map(x -> string(first(x)), entries))]

    cw = [2, 4, 7, 6, 11]
    for (ids, t) in entries
        _update_col_widths!(cw, ids, t)
    end

    print(io, """
                | ID$(" "^(cw[1]-2)) | time$(" "^(cw[2]-4)) | GC time$(" "^(cw[3]-7)) | memory$(" "^(cw[4]-6)) | allocations$(" "^(cw[5]-11)) |
                |---$("-"^(cw[1]-2))-|-----$("-"^(cw[2]-4)):|--------$("-"^(cw[3]-7)):|-------$("-"^(cw[4]-6)):|------------$("-"^(cw[5]-11)):|
                """)

    for (ids, t) in entries
        println(io, _resultrow(ids, t, cw))
    end
    println(io)
    println(io, """
                ## Benchmark Group List
                Here's a list of all the benchmark groups executed by this job:
                """)

    for id in unique(map(pair -> pair[1][1:end-1], entries))
        println(io, "- `", _idrepr(id), "`")
    end

    println(io)
    println(io, "## Julia versioninfo")
    print(io, "```\n", versioninfo(results), "```")

    return nothing
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct EncryptionServices <: SwaggerModel
    blob::Any # spec type: Union{ Nothing, EncryptionService } # spec name: blob
    file::Any # spec type: Union{ Nothing, EncryptionService } # spec name: file
    table::Any # spec type: Union{ Nothing, EncryptionService } # spec name: table
    queue::Any # spec type: Union{ Nothing, EncryptionService } # spec name: queue

    function EncryptionServices(;blob=nothing, file=nothing, table=nothing, queue=nothing)
        o = new()
        validate_property(EncryptionServices, Symbol("blob"), blob)
        setfield!(o, Symbol("blob"), blob)
        validate_property(EncryptionServices, Symbol("file"), file)
        setfield!(o, Symbol("file"), file)
        validate_property(EncryptionServices, Symbol("table"), table)
        setfield!(o, Symbol("table"), table)
        validate_property(EncryptionServices, Symbol("queue"), queue)
        setfield!(o, Symbol("queue"), queue)
        o
    end
end # type EncryptionServices

const _property_map_EncryptionServices = Dict{Symbol,Symbol}(Symbol("blob")=>Symbol("blob"), Symbol("file")=>Symbol("file"), Symbol("table")=>Symbol("table"), Symbol("queue")=>Symbol("queue"))
const _property_types_EncryptionServices = Dict{Symbol,String}(Symbol("blob")=>"EncryptionService", Symbol("file")=>"EncryptionService", Symbol("table")=>"EncryptionService", Symbol("queue")=>"EncryptionService")
Base.propertynames(::Type{ EncryptionServices }) = collect(keys(_property_map_EncryptionServices))
Swagger.property_type(::Type{ EncryptionServices }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_EncryptionServices[name]))}
Swagger.field_name(::Type{ EncryptionServices }, property_name::Symbol) =  _property_map_EncryptionServices[property_name]

function check_required(o::EncryptionServices)
    true
end

function validate_property(::Type{ EncryptionServices }, name::Symbol, val)
end

ext: jl


text:
cd(@__DIR__)

using Pkg

pkg"activate --temp"

Pkg.add("PackageCompiler")

using PackageCompiler

platform = Int==Int64 ? "x64" : "x86"

output_path = joinpath(@__DIR__, "..", "build", "juliaup", platform)

mkpath(output_path)

create_app(".", output_path, filter_stdlibs=true, force=true, precompile_execution_file="precompile_app.jl")

max_stars_repo_name: piever/juliaup

max_issues_repo_name: piever/juliaup

max_forks_repo_name: piever/juliaup


text:
struct Best end # best e for all combinations of a_k

# best e for all combinations of subset of a_k's at positions of ones in c
struct SubSet{T<:Integer}
    a::Vector{T}      # subset to optimize
end

struct SlidingWindow
    width :: Int
end

struct MabeyBest end # optimized version of Best

function search(::Best,s::Integer,n::Integer,::Vector=[])
    a = Array{Int}(1:s-1)
    a_to_n = Array{BigInt}(BigInt.(a).^n)
    lhs = BigInt(s)^n
    bc,best_error = best_combination(a_to_n, lhs)
    collect(OnePositions(bc)),best_error
end

function search(ss::SubSet{T}, s::T, n::Integer, a::Vector{T}) where {T<:Integer}
    @assert length(a)<s
    lhs =  BigInt(s)^n
    a_not_in_subset = setdiff(a,ss.a)
    for a_k in a_not_in_subset
        lhs-=BigInt(a_k)^n
    end
    bc_subset,best_error = best_combination([BigInt(x)^n for x in ss.a], lhs)
    best_a = a_not_in_subset
    sizehint!(best_a,length(best_a)+count_ones(bc_subset))
    for k in OnePositions(bc_subset)
        push!(best_a,ss.a[k])
    end
    best_a, best_error
end

function best_combination(a_to_n,
                          lhs::BigInt,
                          best_error::BigInt = BigInt(lhs)
                          )
    bc = 0
    for combination in 1:2^length(a_to_n)-1
        error = BigInt(lhs)
        for k in OnePositions(combination)
            error-=a_to_n[k]
        end
        if abs(error) < abs(best_error)
            best_error = error
            bc = combination
        end
    end
    bc,best_error
end

function search(sw::SlidingWindow, s::Int, n::Int, a::Vector{Int})
    window_max = s-1
    window_min = s-1-sw.width
    ss = SubSet([x for x in window_max:-1:window_min])
    e = BigInt(0)
    for i in 1:window_min
        a,e = search(ss,s,n,a)
        ss.a.-=1
    end
    a,e
end

function search(::MabeyBest, s::Int, n::Int, ::Vector{Int}=Vector{Int}([]))
    max_k = s-1
    all_a_to_n = [BigInt(x)^n for x in 1:max_k]
    lhs = BigInt(s)^n
    select_a = if max_k<65
        UInt64(0)
    elseif max_k<129
        UInt128(0)
    else
        BigInt(0)
    end
    best = Tracker(select_a,lhs)
    _mabey_best!(select_a, best, lhs, max_k, all_a_to_n)
    best_a, best_error = best()
    Int.(collect(OnePositions(best_a))), best_error
end
function _mabey_best!(select_a::T, best::Tracker, lhs::BigInt, max_k::Int, all_a_to_n) where T<:Integer
    # select_a and best are modified
    is_error_zero(best) && return nothing
    if max_k == 0
        best(select_a, lhs)
        return nothing
    end
    split,e = find_split(lhs, max_k, all_a_to_n)
    if split-1==max_k
        all_ones_to_split = select_a | T(2)^(split-1)-1
        best(all_ones_to_split, e)
        return nothing
    end
    upper_view_a_to_n = view(all_a_to_n, split:max_k)
    lower_view_a_to_n = view(all_a_to_n, 1:split-1)
    upper_combinations_range = 1:T(2)^(max_k-split)
    for upper_combinations in upper_combinations_range
        new_select_a = select_a ⊻ upper_combinations<<(split-1)
        new_lhs = lhs-sum(upper_view_a_to_n[OnePositions(upper_combinations)])
        if new_lhs>0
            _mabey_best!(new_select_a, best, new_lhs, split-1, lower_view_a_to_n)
        else
            best(new_select_a, new_lhs)
        end
    end
    nothing
end
function find_split(lhs::BigInt, max_k::Int, all_a_to_n)
    e = lhs
    for k in 1:max_k
        previous_error = e
        e-=all_a_to_n[k]
        e<=0 && return k,previous_error
    end
    return max_k,e
end
function true_positions(x::BitArray{1})
    positions = Vector{Int}()
    for i in eachindex(x)
        i[x] && push!(positions,i)
    end
    positions
end

function max_s_for_all_a(n)
    s = 1
    old_s_to_n = BigInt(1)
    e = old_s_to_n
    while true
        s+=1
        new_s_to_n = BigInt(s)^n
        e = e + new_s_to_n - 2*old_s_to_n
        e<0 && return s-1
        old_s_to_n = new_s_to_n
    end
end


text:
using GeometricIntegrators
using GeometricIntegrators.Solutions
using Test

dt = Float64
nt = 10

@testset "$(rpad("Dataseries 1d",80))" begin
    nd = 1
    ni = 1

    ds = SDataSeries(dt, nt)
    @test typeof(ds) != typeof(SDataSeries(dt, nd, nt))
    @test typeof(ds) != typeof(SDataSeries(dt, nd, nt, ni))
    @test typeof(ds) <: AbstractArray{dt,1}
    @test firstindex(ds)   == 0
    @test firstindex(ds,1) == firstindex(ds.d,1) - 1
    @test firstindex(ds,2) == 1
    @test lastindex(ds)    == nt
    @test lastindex(ds,1)  == lastindex(ds.d,1) - 1
    @test lastindex(ds,2)  == 1
    @test strides(ds)      == (1,)
    @test stride(ds,1)     == 1
    @test axes(ds)   == (0:nt,)
    @test axes(ds,1) == 0:nt
    @test axes(ds,2) == 1:1
    @test size(ds.d) == (nt+1,)
    @test size(ds.d) == size(ds)
    @test size(ds.d,1) == nt+1
    @test ndims(ds)  == 1
    @test eltype(ds) == dt
    @test parent(ds) == ds.d

    @test eachindex(ds) == 0:nt
    @test eachindex(IndexLinear(), ds) == 0:nt
    @test eachindex(IndexCartesian(), ds) == CartesianIndices((0:nt,))

    for i in 0:nt
        ds[i] = i
    end

    @test ds.d[1] == ds[0]
    @test ds.d[end] == ds[nt]
    @test ds.d[end] == ds[end]
    @test ds.d == collect(0:nt)

    reset!(ds)
    @test ds[0] == ds[end]

    d = rand(nt+1)
    ds = SDataSeries(d)
    @test ds.d == d

    @test similar(ds).d == zero(d)

    for i in 0:nt
        set_data!(ds, dt(i), i)
    end

   @test ds.d == collect(0:nt)
end

@testset "$(rpad("Dataseries 2d",80))" begin
    nd = 2
    ni = 1

    ds = SDataSeries(dt, nd, nt, ni)
    @test typeof(ds) == typeof(SDataSeries(dt, nd, nt))
    @test typeof(ds) <: AbstractArray{dt,2}
    @test firstindex(ds)   == firstindex(ds.d)
    @test firstindex(ds,1) == firstindex(ds.d,1)
    @test firstindex(ds,2) == firstindex(ds.d,2) - 1
    @test firstindex(ds,3) == 1
    @test lastindex(ds)    == lastindex(ds.d)
    @test lastindex(ds,1)  == lastindex(ds.d,1)
    @test lastindex(ds,2)  == lastindex(ds.d,2) - 1
    @test lastindex(ds,3)  == 1
    @test strides(ds)      == (1,nd)
    @test stride(ds,1)     == 1
    @test axes(ds)   == (1:nd, 0:nt)
    @test axes(ds,1) == 1:nd
    @test axes(ds,2) == 0:nt
    @test axes(ds,3) == 1:1
    @test size(ds.d) == (nd, nt+1)
    @test size(ds.d) == size(ds)
    @test size(ds.d,1) == nd
    @test size(ds.d,2) == nt+1
    @test ndims(ds)  == 2
    @test eltype(ds) == dt
    @test parent(ds) == ds.d

    for i in 0:nt
        ds[1,i] = i
        ds[2,i] = 2^i
    end

    tx = zeros(nd)
    for i in 0:nt
        get_data!(ds, tx, i)
        @test tx == Array{eltype(tx)}([i, 2^i])
    end

    @test ds.d[1,1] == ds[1,0]
    @test ds.d[1:ds.nd,1] == ds[1:ds.nd,0]

    @test ds[1,end] == ds.d[1,end]
    @test ds[1,end] == ds[1,nt]
    @test ds[1:ds.nd,nt] == ds[1:ds.nd,end]

    @test ds.d[:,1] == ds[:,0]

    reset!(ds)
    @test ds[1,0] == ds[1,end]

    tx = rand(nd)
    ds[:,0] .= tx
    @test ds.d[:,1] == tx

    d = rand(nd,nt+1)
    ds = SDataSeries(d)
    @test ds.d == d

    @test similar(ds).d == zero(d)
end

@testset "$(rpad("Dataseries 3d",80))" begin
    nd = 2
    ni = 2

    ds = SDataSeries(dt, nd, nt, ni)
    @test typeof(ds) <: AbstractArray{dt,3}
    @test firstindex(ds)   == firstindex(ds.d)
    @test firstindex(ds,1) == firstindex(ds.d,1)
    @test firstindex(ds,2) == firstindex(ds.d,2) - 1
    @test firstindex(ds,3) == firstindex(ds.d,3)
    @test firstindex(ds,4) == 1
    @test lastindex(ds)    == lastindex(ds.d)
    @test lastindex(ds,1)  == lastindex(ds.d,1)
    @test lastindex(ds,2)  == lastindex(ds.d,2) - 1
    @test lastindex(ds,3)  == lastindex(ds.d,3)
    @test lastindex(ds,4)  == 1
    @test strides(ds)      == (1,nd,nd*(nt+1))
    @test stride(ds,1)     == 1
    @test stride(ds,2)     == nd
    @test axes(ds)   == (1:nd, 0:nt, 1:ni)
    @test axes(ds,1) == 1:nd
    @test axes(ds,2) == 0:nt
    @test axes(ds,3) == 1:ni
    @test axes(ds,4) == 1:1
    @test size(ds.d) == (nd, nt+1, ni)
    @test size(ds.d) == size(ds)
    @test size(ds.d,1) == nd
    @test size(ds.d,2) == nt+1
    @test size(ds.d,3) == ni
    @test ndims(ds)  == 3
    @test eltype(ds) == dt
    @test parent(ds) == ds.d

    for j in 1:ni
        for i in 1:nt
            ds[1,i,j] = j*i
            ds[2,i,j] = j*i^2
        end
    end

    tx = zeros(nd)
    ty = zeros(nd,ni)
    tz = zeros(nd,ni)
    for i in 1:nt
        for j in 1:ni
            get_data!(ds, tx, i, j)
            @test tx == Array{eltype(tx)}([j*i, j*i^2])
            tz[1,j] = j*i
            tz[2,j] = j*i^2
        end
        get_data!(ds, ty, i)
        @test ty == tz
    end

    @test ds.d[1,1,1] == ds[1,0,1]
    @test ds.d[1:ds.nd,1,1] == ds[1:ds.nd,0,1]

    @test ds[1,end,1] == ds.d[1,end,1]
    @test ds[1,nt,1] == ds[1,end,1]
    @test ds[1:ds.nd,nt,1] == ds[1:ds.nd,end,1]

    @test ds.d[:,1,1] == ds[:,0,1]
    @test ds.d[:,1,:] == ds[:,0,:]

    reset!(ds)
    @test ds[1,0,1] == ds[1,end,1]

    tx = rand(nd)
    ds[:,0,1] .= tx
    @test ds.d[:,1,1] == tx

    d = rand(nd,nt+1,ni)
    ds = SDataSeries(d)
    @test ds.d == d

    @test similar(ds).d == zero(d)
end


text:
function assemble(iop,qdict;compress=Matrix)
    out        = compress(iop)
    correction = singular_weights(iop,qdict)
    axpy!(1,correction,out) # out <-- out + correction
end

function assemble(iop;compress=Matrix,k=2)
    Y = target_surface(iop)
    dict = Dict{DataType,SingularQuadratureRule}()
    for (E,q) in etypes(Y)
        # use a finer quadrature to compute singular integrals. Oversemple by a
        # factor of k
        q     = Integration.refine(q,k)
        # combine the oversample quadrature with a singularity handler to obtain
        # an appropriate SingularQuadratureRule
        shand = Kress(;order=3)
        qsing = SingularQuadratureRule(q,shand)
        # associate to each element type E a singular quadrature rule
        dict[E] = qsing
    end
    return assemble(iop,dict;compress)
end

function singular_weights(iop::IntegralOperator,qdict)
    X,Y = iop.X, iop.Y
    T = eltype(iop)
    K = iop.kernel
    dict_near = near_interaction_list(qnodes(X),Y;atol=0.1)
    Is    = Int[]
    Js    = Int[]
    Vs    = T[]
    for (E,list_near) in dict_near # for each element type
        qstd      = Y.etype2qrule[E]
        q         = qdict[E] # singular quadrature rule
        _singular_weights!(Is,Js,Vs,iop,E,qstd,q,list_near)
    end
    return Sp = sparse(Is,Js,Vs,size(iop)...)
end

function _singular_weights!(Is,Js,Vs,iop,E,qstd,q,list_near)
    X,Y = iop.X, iop.Y
    K   = iop.kernel
    ui        = qnodes(qstd)
    elt2dof   = Y.elt2dof[E]
    for (n,el) in enumerate(elements(Y)[E]) # for each element of type E in the mesh Y
        j_glob  = elt2dof[:,n]
        for (i,jloc) in list_near[n] # for each near observation point
            x   = X.qnodes[i]   # observation point
            vs  = ui[jloc]    # regularization center in parameter space
            if kernel_type(K) == SingleLayer()
                k    = (v) -> K(x,el(v))*measure(el,v)
            elseif kernel_type(K) == DoubleLayer()
                k    = (v) -> K(x,el(v),normal(el,v)) * measure(el,v)
            elseif kernel_type(K) == AdjointDoubleLayer()
                νx   = X.qnormals[i] # point x must also have a normal
                k    = (v) -> K(x,el(v),νx) * measure(el,v)
            elseif kernel_type(K) == HyperSingular()
                νx  = X.qnormals[i] # point x must also have a normal
                k   = (v) -> K(x,el(v),νx,normal(el,v)) * measure(el,v)
            else
                error("unkown kernel type")
            end
            w  = singular_weights(k,qstd,q,vs)
            append!(Is,fill(i,length(j_glob)))
            append!(Js,j_glob)
            w  = w - iop[i,j_glob]
            append!(Vs,w)
        end
    end
end

max_stars_repo_name: WaveProp/WaveProp

max_issues_repo_name: WaveProp/WaveProp

max_forks_repo_name: WaveProp/WaveProp


text:
# The library of functions for the simple SV model

using Statistics, Random

# method that generates the sample
function auxstat(θ, reps)
    stats = zeros(reps,11)
    if InSupport(θ)
        n = 500
        burnin = 100
        for rep = 1:reps
            y = SVmodel(θ, n, burnin)
            stats[rep,:] = auxstat(y)
        end
    end
    stats
end

# method that uses pre-generated shocks
function auxstat(θ, ϵ, u)
    ϕ = θ[1]
    ρ = θ[2]
    σ = θ[3]
    reps = size(ϵ,2)
    stats = zeros(reps,11)
    if InSupport(θ)
        n = 500
        burnin = 100
        for rep = 1:reps
            hlag = 0.0
            ys = zeros(n)
            @inbounds for t = 1:burnin+n
                h = ρ*hlag + σ*u[t, rep]
                y = ϕ*exp(h/2.0)*ϵ[t, rep]
                if t > burnin
                    ys[t-burnin] = y
                end
                hlag = h
            end
            stats[rep,:] = auxstat(ys)
        end
    end
    stats
end

# method for a given sample
function auxstat(y)
	s = std(y)
	y = abs.(y)
	m = mean(y)
	s2 = std(y)
	y = y ./ s2
	k = std((y).^2.0)
	c = cor(y[1:end-1],y[2:end])
	# ratios of quantiles of moving averages to detect clustering
	q = try
	    q = quantile((ma(y,3)[3:end]), [0.25, 0.75])
	catch
	    q = [1.0, 1.0]
	end
	c1 = log(q[2]/q[1])
	stats = sqrt(size(y,1)) .* vcat(m, s, s2, k, c, c1, HAR(y))
end

function SVmodel(θ, n, burnin)
    ϕ = θ[1]
    ρ = θ[2]
    σ = θ[3]
    hlag = 0.0
    ys = zeros(n)
    if InSupport(θ)
        @inbounds for t = 1:burnin+n
            h = ρ*hlag + σ*randn()
            y = ϕ*exp(h/2.0)*randn()
            if t > burnin
                ys[t-burnin] = y
            end
            hlag = h
        end
    end
    ys
end

function TrueParameters()
    [exp(-0.736/2.0), 0.9, 0.363]
end

function PriorSupport()
    lb = [0.05, 0.0, 0.05]
    ub = [2.0, 0.999, 1.0]
    lb,ub
end

# prior checks that we're in the bounds, and that the unconditional std. dev. of log vol is not too high
# returns 1 if this is true, zero otherwise. Value is not important, as it's constant
function Prior(θ)
    InSupport(θ) ? 1.0 : 0.0
end

# check if parameter is in support. In this case, we require
# the bounds, and that the unconditional variance of the volatility
# shock be limited
function InSupport(θ)
    lb,ub = PriorSupport()
    if all(θ .>= lb) & all(θ .<= ub)
	return  (θ[3]/sqrt(1.0 - θ[2]^2.0) < 5.0)
    else
    return false
    end
end

function PriorDraw()
    lb, ub = PriorSupport()
    ok = false
    θ = 0.0
    while !ok
        θ = (ub-lb).*rand(size(lb,1)) + lb
        ok = InSupport(θ)
    end
    return θ
end

# taken from https://github.com/mcreel/Econometrics
# # returns the variable (or matrix), lagged p times,
# with the first p rows filled with ones (to avoid divide errors)
# remember to drop those rows before doing analysis
function lag(x::Array{Float64,2},p::Int64)
	n,k = size(x)
	lagged_x = [ones(p,k); x[1:n-p,:]]
end

function lag(x::Array{Float64,1},p::Int64)
	n = size(x,1)
	lagged_x = [ones(p); x[1:n-p]]
end

# taken from https://github.com/mcreel/Econometrics
# returns the variable (or matrix), lagged from 1 to p times,
# with the first p rows filled with ones (to avoid divide errors)
# remember to drop those rows before doing analysis
function  lags(x::Array{Float64,2},p)
	n, k = size(x)
	lagged_x = zeros(eltype(x),n,p*k)
	for i = 1:p
		lagged_x[:,i*k-k+1:i*k] = lag(x,i)
	end
    return lagged_x
end

function  lags(x::Array{Float64,1},p)
	n = size(x,1)
	lagged_x = zeros(eltype(x), n,p)
	for i = 1:p
		lagged_x[:,i] = lag(x,i)
	end
    return lagged_x
end

# taken from https://github.com/mcreel/Econometrics
# compute moving average using p most recent values, including current value
function ma(x, p)
    m = zeros(size(x))
    for i = p:size(x,1)
        m[i] = mean(x[i-p+1:i])
    end
    return m
end

# auxiliary model: HAR-RV
# Corsi, Fulvio. "A simple approximate long-memory model
# of realized volatility." Journal of Financial Econometrics 7,
# no. 2 (2009): 174-196.
function HAR(y)
    ylags = lags(y,10)
    X = [ones(size(y,1)) ylags[:,1]  mean(ylags[:,1:4],dims=2) mean(ylags[:,1:10],dims=2)]
    # drop missings
    y = y[11:end]
    X = X[11:end,:]
    βhat = X\y
    σhat = std(y-X*βhat)
    vcat(βhat,σhat)
end

max_stars_repo_name: Hiroakiyusheng/Econometrics

max_issues_repo_name: Hiroakiyusheng/Econometrics

max_forks_repo_name: Hiroakiyusheng/Econometrics


text:
using ConvertBNG
using Printf

const NUM_COORDS = 1_000_000
const REPEAT_RUNS = 50
const LONS_LATS = Ref{Array{Float64,2}}()
const ENS = Ref{Array{Float64,2}}()

function initial_setup(num_coords)
    # London bounding box
    N = 51.691874117
    E =  0.334015564
    S = 51.286760163
    W = -0.510375069

    # Create random testset of lon, lats
    lons = rand(W:0.01:E, num_coords, 1)
    lats = rand(S:0.01:N, num_coords, 1)

    # Create random testset of eastings, northings
    e_max, n_max = convert_bng(E, N)
    e_min, n_min = convert_bng(W, S)
    eastings = rand(e_min:0.1:e_max, num_coords, 1)
    northings = rand(n_min:0.1:n_max, num_coords, 1)

    return hcat(lons, lats), hcat(eastings, northings)
end

function benchmark(num_coords::Integer, repeat_runs::Integer)
    # Run benchmark for convert_bng
    println("Running convert_bng benchmark for $repeat_runs repeat runs of $num_coords points...")
    starttime = time()
    for x in 1:repeat_runs
        convert_bng(LONS_LATS[])
    end
    totaltime = @sprintf("%.2f", time() - starttime)
    println("Time taken: $totaltime s")

    # Run benchmark for convert_lonlat
    println("Running convert_lonlat benchmark for $repeat_runs repeat runs of $num_coords points...")
    starttime = time()
    for x in 1:repeat_runs
        convert_lonlat(ENS[])
    end
    totaltime = @sprintf("%.2f", time() - starttime)
    println("Time taken: $totaltime s")
    print("Benchmarks complete!")
end

# Setup test coordinate arrays
LONS_LATS[], ENS[] = initial_setup(NUM_COORDS)

# Run benchmarks
benchmark(NUM_COORDS, REPEAT_RUNS)


text:
export tforeach;

"""
	tforeach(f::Function, c::AbstractArray)::Nothing

Multi-threaded version of [foreach(f, c)](https://docs.julialang.org/en/v1.2/base/collections/#Base.foreach).
"""
function tforeach(f::Function, c...)::Nothing
	ensureThreaded();
	Threads.@threads for i in minimum(eachindex.(c))
		f(getindex.(c, i)...);
	end
end


text:
struct SAFTVRQMieParam <: EoSParam
    segment::SingleParam{Float64}
    sigma::PairParam{Float64}
    lambda_a::PairParam{Float64}
    lambda_r::PairParam{Float64}
    epsilon::PairParam{Float64}
    Mw::SingleParam{Float64}
end

abstract type SAFTVRQMieModel <: SAFTVRMieModel end
@newmodel SAFTVRQMie SAFTVRQMieModel SAFTVRQMieParam

export SAFTVRQMie
function SAFTVRQMie(components::Array{<:Any,1}; idealmodel::Type=BasicIdeal, userlocations::Array{String,1}=String[], verbose=false)
    params = getparams(components, ["SAFT/SAFTVRQMie"]; userlocations=userlocations, verbose=verbose)

    params["Mw"].values .*= 1E-3
    Mw = params["Mw"]
    segment = params["m"]
    params["sigma"].values .*= 1E-10
    sigma = sigma_LorentzBerthelot(params["sigma"])
    epsilon = epsilon_HudsenMcCoubrey(params["epsilon"], sigma)
    lambda_a = lambda_LorentzBerthelot(params["lambda_a"])
    lambda_r = lambda_LorentzBerthelot(params["lambda_r"])

    sites = SiteParam(Dict("e" => params["n_e"], "H" => params["n_H"]))

    packagedparams = SAFTVRQMieParam(segment, sigma, lambda_a, lambda_r, epsilon, Mw)
    references = ["todo"]

    model = SAFTVRQMie(packagedparams, sites, idealmodel; references=references, verbose=verbose)
    return model
end

function a_res(model::SAFTVRQMieModel, V, T, z)
    return @f(a_mono)
end

function d(model::SAFTVRQMieModel, V, T, z, i)
    ϵ = model.params.epsilon.diagvalues
    σ = model.params.sigma.diagvalues
    λr = model.params.lambda_r.diagvalues
    λa = model.params.lambda_a.diagvalues
    Mw = model.params.Mw.values

    Di = ħ^2/(12*k_B*T*Mw[i]/N_A*σ[i]^2)
    βi = ϵ[i]/T
    u(x) = @f(C,i,i)*((x^-λr[i]-x^-λa[i])+
              Di*(@f(Q1,λr[i])*x^-(λr[i]+2)-@f(Q1,λa[i])*x^-(λa[i]+2))+
              Di^2*(@f(Q2,λr[i])*x^-(λr[i]+4)-@f(Q2,λa[i])*x^-(λa[i]+4)))
    du(x) = -@f(C,i,i)*((λr[i]*x^-(λr[i]+1)-λa[i]*x^-(λa[i]+1))+
                Di*(@f(Q1,λr[i])*(λr[i]+2)*x^-(λr[i]+3)-@f(Q1,λa[i])*(λa[i]+2)*x^-(λa[i]+3))+
                Di^2*(@f(Q2,λr[i])*(λr[i]+4)*x^-(λr[i]+5)-@f(Q2,λa[i])*(λa[i]+4)*x^-(λa[i]+5)))
    d2u(x) = @f(C,i,i)*((λr[i]*(λr[i]+1)*x^-(λr[i]+2)-λa[i]*(λa[i]+1)*x^-(λa[i]+2))+
                Di*(@f(Q1,λr[i])*(λr[i]+2)*(λr[i]+3)*x^-(λr[i]+4)-@f(Q1,λa[i])*(λa[i]+2)*(λa[i]+3)*x^-(λa[i]+4))+
                Di^2*(@f(Q2,λr[i])*(λr[i]+4)*(λr[i]+5)*x^-(λr[i]+6)-@f(Q2,λa[i])*(λa[i]+4)*(λa[i]+5)*x^-(λa[i]+6)))

    f(x) = exp(-βi*u(x))
    g(x) = -βi*exp(-βi*u(x))*du(x)
    h(x) = βi*exp(-βi*u(x))*(βi*du(x)^2-d2u(x))
    x_min = @f(Halley,f,g,h,1.)
    σ_effi = @f(σeff,i,i)/σ[i]
    x = SAFTVRQMieconsts.x
    w = SAFTVRQMieconsts.w
    return σ[i]*(σ_effi-(σ_effi-x_min)/2*∑(w[i]*f((σ_effi-x_min)/2*x[i]+(σ_effi+x_min)/2) for i in 1:21))
end

function σeff(model::SAFTVRQMieModel, V, T, z, i, j)
    σ = model.params.sigma.values
    λr = model.params.lambda_r.values
    λa = model.params.lambda_a.values
    Mwij = (model.params.Mw.values[i] + model.params.Mw.values[j])/2 # check
    Dij = ħ^2/(12*k_B*T*Mwij/N_A*σ[i,j]^2)

    f(x) = Dij^2*@f(Q2,λr[i,j])+Dij*@f(Q1,λr[i,j])*x^2+x^4-x^(λr[i,j]-λa[i,j])*(Dij^2*@f(Q2,λa[i,j])+Dij*@f(Q1,λa[i,j])*x^2+x^4)
    g(x) = 2*Dij*@f(Q1,λr[i,j])*x+4x^3-x^(λr[i,j]-λa[i,j]-1)*((λr[i,j]-λa[i,j])*Dij^2*@f(Q2,λa[i,j])+(λr[i,j]-λa[i,j]+2)*Dij*@f(Q1,λa[i,j])*x^2+(λr[i,j]-λa[i,j]+4)*x^4)
    h(x) = 2*Dij*@f(Q1,λr[i,j])+12x^2-x^(λr[i,j]-λa[i,j]-2)*((λr[i,j]-λa[i,j]-1)*(λr[i,j]-λa[i,j])*Dij^2*@f(Q2,λa[i,j])+((λr[i,j]-λa[i,j]-1)*(λr[i,j]-λa[i,j]+2)+2*(λr[i,j]-λa[i,j]+2))*Dij*@f(Q1,λa[i,j])*x^2+((λr[i,j]-λa[i,j]-1)*(λr[i,j]-λa[i,j]+4)+4*(λr[i,j]-λa[i,j]+4))*x^4)

    return σ[i,j]*@f(Halley,f,g,h,1.)
end

function ϵeff(model::SAFTVRQMieModel, V, T, z, i, j)
    ϵ = model.params.epsilon.values
    σ = model.params.sigma.values
    λr = model.params.lambda_r.values
    λa = model.params.lambda_a.values
    Mwij = (model.params.Mw.values[i] + model.params.Mw.values[j])/2 # check
    Dij = ħ^2/(12*k_B*T*Mwij/N_A*σ[i,j]^2)

    u(x) = @f(C,i,j)*((x^-λr[i,j]-x^-λa[i,j])+
            Dij*(@f(Q1,λr[i,j])*x^-(λr[i,j]+2)-@f(Q1,λa[i,j])*x^-(λa[i,j]+2))+
            Dij^2*(@f(Q2,λr[i,j])*x^-(λr[i,j]+4)-@f(Q2,λa[i,j])*x^-(λa[i,j]+4)))
    f(x) = (λr[i,j]+4)*Dij^2*@f(Q2,λr[i,j])+(λr[i,j]+2)*Dij*@f(Q1,λr[i,j])*x^2+λr[i,j]*x^4-x^(λr[i,j]-λa[i,j])*((λa[i,j]+4)*Dij^2*@f(Q2,λa[i,j])+(λa[i,j]+2)*Dij*@f(Q1,λa[i,j])*x^2+λa[i,j]*x^4)
    g(x) = 2*(λr[i,j]+2)*Dij*@f(Q1,λr[i,j])*x+4λr[i,j]*x^3-x^(λr[i,j]-λa[i,j]-1)*((λa[i,j]+4)*(λr[i,j]-λa[i,j])*Dij^2*@f(Q2,λa[i,j])+(λa[i,j]+2)*(λr[i,j]-λa[i,j]+2)*Dij*@f(Q1,λa[i,j])*x^2+λa[i,j]*(λr[i,j]-λa[i,j]+4)*x^4)
    h(x) = 2*(λr[i,j]+2)*Dij*@f(Q1,λr[i,j])+12λr[i,j]*x^2-x^(λr[i,j]-λa[i,j]-2)*((λa[i,j]+4)*(λr[i,j]-λa[i,j]-1)*(λr[i,j]-λa[i,j])*Dij^2*@f(Q2,λa[i,j])+(λa[i,j]+2)*((λr[i,j]-λa[i,j]-1)*(λr[i,j]-λa[i,j]+2)+2*(λr[i,j]-λa[i,j]+2))*Dij*@f(Q1,λa[i,j])*x^2+λa[i,j]*((λr[i,j]-λa[i,j]-1)*(λr[i,j]-λa[i,j]+4)+4*(λr[i,j]-λa[i,j]+4))*x^4)

    σ_min = @f(Halley,f,g,h,(λr[i,j]/λa[i,j])^(1/(λr[i,j]-λa[i,j])))
    return -ϵ[i,j]*u(σ_min)
end

function Halley(model::SAFTVRQMieModel, V, T, z, f_, g_, h_, x0)
    tolx = 1.
    tolf = 1.
    f0 = f_(x0)
    if f0 < 1e-16
        return x0
    else
        while tolf > 1e-16 && tolx > 1e-15
            d = f_(x0)/g_(x0)*(1-f_(x0)*h_(x0)/(2*g_(x0)^2))^-1
            x0=x0-d
            f0 = f_(x0)
            tolf = abs(f0)
            tolx = abs(d)
        end
        return x0
    end
end

function ζ_x(model::SAFTVRQMieModel, V, T, z)
    return π/6*@f(ρ_S)*∑(@f(x_S,i)*@f(x_S,j)*(@f(d,i,j))^3 for i in @comps for j in @comps)
end

function a_1(model::SAFTVRQMieModel, V, T, z, i, j)
    ϵ = model.params.epsilon.values
    σ = model.params.sigma.values
    λr = model.params.lambda_r.values
    λa = model.params.lambda_a.values
    Mwij = (model.params.Mw.values[i] + model.params.Mw.values[j])/2 # check
    Dij = ħ^2/(12*k_B*T*Mwij/N_A*σ[i,j]^2)
    x_0ij = @f(x_0,i,j)
    x_0effij = @f(x_0eff,i,j)
    return 2*π*ϵ[i,j]*@f(d,i,j)^3*@f(C,i,j)*@f(ρ_S)*
    ( (x_0ij^λa[i,j]*(@f(aS_1,λa[i,j])+@f(B,λa[i,j],x_0effij))-
        x_0ij^λr[i,j]*(@f(aS_1,λr[i,j])+@f(B,λr[i,j],x_0effij)))+
        (x_0ij^(λa[i,j]+2)*@f(Q1,λa[i,j])*(@f(aS_1,λa[i,j]+2)+@f(B,λa[i,j]+2,x_0effij))-
         x_0ij^(λr[i,j]+2)*@f(Q1,λr[i,j])*(@f(aS_1,λr[i,j]+2)+@f(B,λr[i,j]+2,x_0effij)))*Dij+
        (x_0ij^(λa[i,j]+4)*@f(Q2,λa[i,j])*(@f(aS_1,λa[i,j]+4)+@f(B,λa[i,j]+4,x_0effij))-
         x_0ij^(λr[i,j]+4)*@f(Q2,λr[i,j])*(@f(aS_1,λr[i,j]+4)+@f(B,λr[i,j]+4,x_0effij)))*Dij^2 )
end

function Q1(model::SAFTVRQMieModel, V, T, z, λ)
    return λ*(λ-1)
end

function Q2(model::SAFTVRQMieModel, V, T, z, λ)
    return 1/2*(λ+2)*(λ+1)*λ*(λ-1)
end

function x_0eff(model::SAFTVRQMieModel, V, T, z, i, j)
    return @f(σeff,i,j)/@f(d,i,j)
end

function a_2(model::SAFTVRQMieModel, V, T, z, i, j)
    ϵ = model.params.epsilon.values
    σ = model.params.sigma.values
    λr = model.params.lambda_r.values
    λa = model.params.lambda_a.values
    Mwij = (model.params.Mw.values[i] + model.params.Mw.values[j])/2 # check
    Dij = ħ^2/(12*k_B*T*Mwij/N_A*σ[i,j]^2)

    x_0ij = @f(x_0,i,j)
    x_0effij = @f(x_0eff,i,j)
    return π*@f(KHS)*(1+@f(χ,i,j))*@f(ρ_S)*ϵ[i,j]^2*@f(d,i,j)^3*@f(C,i,j)^2*(x_0ij^(2*λa[i,j])*(@f(aS_1,2*λa[i,j])+@f(B,2*λa[i,j],x_0effij))-
        x_0ij^(λa[i,j]+λr[i,j]+8)*(2*@f(Q2,λa[i,j])*@f(Q2,λr[i,j]))*(@f(aS_1,λa[i,j]+λr[i,j]+8)+@f(B,λa[i,j]+λr[i,j]+8,x_0effij))*Dij^4)
end

function χ(model::SAFTVRQMieModel, V, T, z,i,j)
    λr = model.params.lambda_r.values
    λa = model.params.lambda_a.values
    σ = model.params.sigma.values
    ϵ = model.params.epsilon.values
    σeffij = @f(σeff,i,j)
    ϵeffij = @f(ϵeff,i,j)
    Mwij = (model.params.Mw.values[i] + model.params.Mw.values[j])/2 # check
    Dij = ħ^2/(12*k_B*T*Mwij/N_A*σ[i,j]^2)

    ζst_ = ζst(model, V, T, z)
    α = @f(C,i,j)*ϵ[i,j]/ϵeffij*
    (((σ[i,j]/σeffij)^λa[i,j]/(λa[i,j]-3)-(σ[i,j]/σeffij)^λr[i,j]/(λr[i,j]-3))+
        Dij*((σ[i,j]/σeffij)^(2+λa[i,j])*@f(Q1,λa[i,j])/(λa[i,j]-1)-
             (σ[i,j]/σeffij)^(2+λr[i,j])*@f(Q1,λr[i,j])/(λr[i,j]-1))+
        Dij^2*((σ[i,j]/σeffij)^(4+λa[i,j])*@f(Q2,λa[i,j])/(λa[i,j]+1)-
               (σ[i,j]/σeffij)^(4+λr[i,j])*@f(Q2,λr[i,j])/(λr[i,j]+1)))
    return @f(f,α,1)*ζst_+@f(f,α,2)*ζst_^5+@f(f,α,3)*ζst_^8
end

function ζst(model::SAFTVRQMieModel, V, T, z)
    return @f(ρ_S)*π/6*∑(@f(x_S,i)*@f(x_S,j)*(@f(σeff,i,j))^3 for i in @comps for j in @comps)
end

function a_3(model::SAFTVRQMieModel, V, T, z, i, j)
    λr = model.params.lambda_r.values
    λa = model.params.lambda_a.values
    σ = model.params.sigma.values
    ϵ = model.params.epsilon.values
    σeffij = @f(σeff,i,j)
    ϵeffij = @f(ϵeff,i,j)
    Mwij = (model.params.Mw.values[i] + model.params.Mw.values[j])/2 # check
    Dij = ħ^2/(12*k_B*T*Mwij/N_A*σ[i,j]^2)

    ζst_ = ζst(model, V, T, z)
    α = @f(C,i,j)*ϵ[i,j]/ϵeffij*
    (((σ[i,j]/σeffij)^λa[i,j]/(λa[i,j]-3)-(σ[i,j]/σeffij)^λr[i,j]/(λr[i,j]-3))+
        Dij*((σ[i,j]/σeffij)^(2+λa[i,j])*@f(Q1,λa[i,j])/(λa[i,j]-1)-
             (σ[i,j]/σeffij)^(2+λr[i,j])*@f(Q1,λr[i,j])/(λr[i,j]-1))+
        Dij^2*((σ[i,j]/σeffij)^(4+λa[i,j])*@f(Q2,λa[i,j])/(λa[i,j]+1)-
               (σ[i,j]/σeffij)^(4+λr[i,j])*@f(Q2,λr[i,j])/(λr[i,j]+1)))
    return -ϵeffij^3*@f(f,α,4)*ζst_*exp(@f(f,α,5)*ζst_+@f(f,α,6)*ζst_^2)
end

const SAFTVRQMieconsts = (
    x = [-0.9956571630258080807355,-0.973906528517171720078,-0.9301574913557082260012,-0.8650633666889845107321,-0.7808177265864168970637,-0.6794095682990244062343,-0.562757134668604683339,
        -0.4333953941292471907993,-0.294392862701460198131,-0.1488743389816312108848,0,0.1488743389816312108848,0.2943928627014601981311,0.4333953941292471907993,
        0.562757134668604683339,0.6794095682990244062343,0.7808177265864168970637,0.865063366688984510732,0.9301574913557082260012,0.973906528517171720078,0.9956571630258080807355],
    w = [0.0116946388673718742781,0.0325581623079647274788,0.0547558965743519960314,0.075039674810919952767,0.093125454583697605535,0.1093871588022976418992,0.123491976262065851078,
        0.134709217311473325928,0.142775938577060080797,0.1477391049013384913748,0.149445554002916905665,0.1477391049013384913748,0.1427759385770600807971,0.134709217311473325928,
        0.123491976262065851078,0.109387158802297641899,0.093125454583697605535,0.075039674810919952767,0.05475589657435199603138,0.032558162307964727479,0.0116946388673718742781],
)


text:
using Literate
using Weave
path = @__DIR__
cd(joinpath(path,".."))
filename = "DesitionTree"
inputfile = joinpath(path,filename*".jl")
outputdir = joinpath(path,"..","Doc")
Literate.markdown(inputfile, outputdir)
mdfile = joinpath(path,"..","Doc",filename*".md")
weave(mdfile, out_path=outputdir)


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct VirtualNetworkGatewayConnectionPropertiesFormat <: SwaggerModel
    authorizationKey::Any # spec type: Union{ Nothing, String } # spec name: authorizationKey
    virtualNetworkGateway1::Any # spec type: Union{ Nothing, VirtualNetworkGateway } # spec name: virtualNetworkGateway1
    virtualNetworkGateway2::Any # spec type: Union{ Nothing, VirtualNetworkGateway } # spec name: virtualNetworkGateway2
    localNetworkGateway2::Any # spec type: Union{ Nothing, LocalNetworkGateway } # spec name: localNetworkGateway2
    connectionType::Any # spec type: Union{ Nothing, String } # spec name: connectionType
    routingWeight::Any # spec type: Union{ Nothing, Int32 } # spec name: routingWeight
    sharedKey::Any # spec type: Union{ Nothing, String } # spec name: sharedKey
    connectionStatus::Any # spec type: Union{ Nothing, String } # spec name: connectionStatus
    tunnelConnectionStatus::Any # spec type: Union{ Nothing, Vector{TunnelConnectionHealth} } # spec name: tunnelConnectionStatus
    egressBytesTransferred::Any # spec type: Union{ Nothing, Int64 } # spec name: egressBytesTransferred
    ingressBytesTransferred::Any # spec type: Union{ Nothing, Int64 } # spec name: ingressBytesTransferred
    peer::Any # spec type: Union{ Nothing, SubResource } # spec name: peer
    enableBgp::Any # spec type: Union{ Nothing, Bool } # spec name: enableBgp
    usePolicyBasedTrafficSelectors::Any # spec type: Union{ Nothing, Bool } # spec name: usePolicyBasedTrafficSelectors
    ipsecPolicies::Any # spec type: Union{ Nothing, Vector{IpsecPolicy} } # spec name: ipsecPolicies
    resourceGuid::Any # spec type: Union{ Nothing, String } # spec name: resourceGuid
    provisioningState::Any # spec type: Union{ Nothing, String } # spec name: provisioningState

    function VirtualNetworkGatewayConnectionPropertiesFormat(;authorizationKey=nothing, virtualNetworkGateway1=nothing, virtualNetworkGateway2=nothing, localNetworkGateway2=nothing, connectionType=nothing, routingWeight=nothing, sharedKey=nothing, connectionStatus=nothing, tunnelConnectionStatus=nothing, egressBytesTransferred=nothing, ingressBytesTransferred=nothing, peer=nothing, enableBgp=nothing, usePolicyBasedTrafficSelectors=nothing, ipsecPolicies=nothing, resourceGuid=nothing, provisioningState=nothing)
        o = new()
        validate_property(VirtualNetworkGatewayConnectionPropertiesFormat, Symbol("authorizationKey"), authorizationKey)
        setfield!(o, Symbol("authorizationKey"), authorizationKey)
        validate_property(VirtualNetworkGatewayConnectionPropertiesFormat, Symbol("virtualNetworkGateway1"), virtualNetworkGateway1)
        setfield!(o, Symbol("virtualNetworkGateway1"), virtualNetworkGateway1)
        validate_property(VirtualNetworkGatewayConnectionPropertiesFormat, Symbol("virtualNetworkGateway2"), virtualNetworkGateway2)
        setfield!(o, Symbol("virtualNetworkGateway2"), virtualNetworkGateway2)
        validate_property(VirtualNetworkGatewayConnectionPropertiesFormat, Symbol("localNetworkGateway2"), localNetworkGateway2)
        setfield!(o, Symbol("localNetworkGateway2"), localNetworkGateway2)
        validate_property(VirtualNetworkGatewayConnectionPropertiesFormat, Symbol("connectionType"), connectionType)
        setfield!(o, Symbol("connectionType"), connectionType)
        validate_property(VirtualNetworkGatewayConnectionPropertiesFormat, Symbol("routingWeight"), routingWeight)
        setfield!(o, Symbol("routingWeight"), routingWeight)
        validate_property(VirtualNetworkGatewayConnectionPropertiesFormat, Symbol("sharedKey"), sharedKey)
        setfield!(o, Symbol("sharedKey"), sharedKey)
        validate_property(VirtualNetworkGatewayConnectionPropertiesFormat, Symbol("connectionStatus"), connectionStatus)
        setfield!(o, Symbol("connectionStatus"), connectionStatus)
        validate_property(VirtualNetworkGatewayConnectionPropertiesFormat, Symbol("tunnelConnectionStatus"), tunnelConnectionStatus)
        setfield!(o, Symbol("tunnelConnectionStatus"), tunnelConnectionStatus)
        validate_property(VirtualNetworkGatewayConnectionPropertiesFormat, Symbol("egressBytesTransferred"), egressBytesTransferred)
        setfield!(o, Symbol("egressBytesTransferred"), egressBytesTransferred)
        validate_property(VirtualNetworkGatewayConnectionPropertiesFormat, Symbol("ingressBytesTransferred"), ingressBytesTransferred)
        setfield!(o, Symbol("ingressBytesTransferred"), ingressBytesTransferred)
        validate_property(VirtualNetworkGatewayConnectionPropertiesFormat, Symbol("peer"), peer)
        setfield!(o, Symbol("peer"), peer)
        validate_property(VirtualNetworkGatewayConnectionPropertiesFormat, Symbol("enableBgp"), enableBgp)
        setfield!(o, Symbol("enableBgp"), enableBgp)
        validate_property(VirtualNetworkGatewayConnectionPropertiesFormat, Symbol("usePolicyBasedTrafficSelectors"), usePolicyBasedTrafficSelectors)
        setfield!(o, Symbol("usePolicyBasedTrafficSelectors"), usePolicyBasedTrafficSelectors)
        validate_property(VirtualNetworkGatewayConnectionPropertiesFormat, Symbol("ipsecPolicies"), ipsecPolicies)
        setfield!(o, Symbol("ipsecPolicies"), ipsecPolicies)
        validate_property(VirtualNetworkGatewayConnectionPropertiesFormat, Symbol("resourceGuid"), resourceGuid)
        setfield!(o, Symbol("resourceGuid"), resourceGuid)
        validate_property(VirtualNetworkGatewayConnectionPropertiesFormat, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        o
    end
end # type VirtualNetworkGatewayConnectionPropertiesFormat

const _property_map_VirtualNetworkGatewayConnectionPropertiesFormat = Dict{Symbol,Symbol}(Symbol("authorizationKey")=>Symbol("authorizationKey"), Symbol("virtualNetworkGateway1")=>Symbol("virtualNetworkGateway1"), Symbol("virtualNetworkGateway2")=>Symbol("virtualNetworkGateway2"), Symbol("localNetworkGateway2")=>Symbol("localNetworkGateway2"), Symbol("connectionType")=>Symbol("connectionType"), Symbol("routingWeight")=>Symbol("routingWeight"), Symbol("sharedKey")=>Symbol("sharedKey"), Symbol("connectionStatus")=>Symbol("connectionStatus"), Symbol("tunnelConnectionStatus")=>Symbol("tunnelConnectionStatus"), Symbol("egressBytesTransferred")=>Symbol("egressBytesTransferred"), Symbol("ingressBytesTransferred")=>Symbol("ingressBytesTransferred"), Symbol("peer")=>Symbol("peer"), Symbol("enableBgp")=>Symbol("enableBgp"), Symbol("usePolicyBasedTrafficSelectors")=>Symbol("usePolicyBasedTrafficSelectors"), Symbol("ipsecPolicies")=>Symbol("ipsecPolicies"), Symbol("resourceGuid")=>Symbol("resourceGuid"), Symbol("provisioningState")=>Symbol("provisioningState"))
const _property_types_VirtualNetworkGatewayConnectionPropertiesFormat = Dict{Symbol,String}(Symbol("authorizationKey")=>"String", Symbol("virtualNetworkGateway1")=>"VirtualNetworkGateway", Symbol("virtualNetworkGateway2")=>"VirtualNetworkGateway", Symbol("localNetworkGateway2")=>"LocalNetworkGateway", Symbol("connectionType")=>"String", Symbol("routingWeight")=>"Int32", Symbol("sharedKey")=>"String", Symbol("connectionStatus")=>"String", Symbol("tunnelConnectionStatus")=>"Vector{TunnelConnectionHealth}", Symbol("egressBytesTransferred")=>"Int64", Symbol("ingressBytesTransferred")=>"Int64", Symbol("peer")=>"SubResource", Symbol("enableBgp")=>"Bool", Symbol("usePolicyBasedTrafficSelectors")=>"Bool", Symbol("ipsecPolicies")=>"Vector{IpsecPolicy}", Symbol("resourceGuid")=>"String", Symbol("provisioningState")=>"String")
Base.propertynames(::Type{ VirtualNetworkGatewayConnectionPropertiesFormat }) = collect(keys(_property_map_VirtualNetworkGatewayConnectionPropertiesFormat))
Swagger.property_type(::Type{ VirtualNetworkGatewayConnectionPropertiesFormat }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_VirtualNetworkGatewayConnectionPropertiesFormat[name]))}
Swagger.field_name(::Type{ VirtualNetworkGatewayConnectionPropertiesFormat }, property_name::Symbol) =  _property_map_VirtualNetworkGatewayConnectionPropertiesFormat[property_name]

const _allowed_VirtualNetworkGatewayConnectionPropertiesFormat_connectionType = ["IPsec", "Vnet2Vnet", "ExpressRoute", "VPNClient"]

const _allowed_VirtualNetworkGatewayConnectionPropertiesFormat_connectionStatus = ["Unknown", "Connecting", "Connected", "NotConnected"]

function check_required(o::VirtualNetworkGatewayConnectionPropertiesFormat)
    (getproperty(o, Symbol("virtualNetworkGateway1")) === nothing) && (return false)
    true
end

function validate_property(::Type{ VirtualNetworkGatewayConnectionPropertiesFormat }, name::Symbol, val)
    if name === Symbol("connectionType")
        Swagger.validate_param(name, "VirtualNetworkGatewayConnectionPropertiesFormat", :enum, val, _allowed_VirtualNetworkGatewayConnectionPropertiesFormat_connectionType)
    end
    if name === Symbol("connectionStatus")
        Swagger.validate_param(name, "VirtualNetworkGatewayConnectionPropertiesFormat", :enum, val, _allowed_VirtualNetworkGatewayConnectionPropertiesFormat_connectionStatus)
    end
end


text:
module Readinputfile_module
    include("data_module.jl")
    using .Data_module

    mutable struct Readinputfile_variables
        data::Data_module.Data_val
        lines::Array{String, 1}
        lineindex::UInt64
    end
end


text:
using Random, Distributions, Plots; pyplot()
Random.seed!(0)

lambda  = 1/4.5
expDist = Exponential(1/lambda)
n, N    = 10, 10^6

means     = Array{Float64}(undef, N)
variances = Array{Float64}(undef, N)

for i in 1:N
    data = rand(expDist,n)
    means[i] = mean(data)
    variances[i] = var(data)
end

println("Actual mean: ",mean(expDist),
		"\nMean of sample means: ",mean(means))
println("Actual variance: ",var(expDist),
	"\nMean of sample variances: ",mean(variances))

stephist(means, bins=200, c=:blue, normed=true,
	label="Histogram of Sample Means")
stephist!(variances, bins=600, c=:red, normed=true,
	label="Histogram of Sample Variances", xlims=(0,40), ylims=(0,0.4),
    	xlabel = "Statistic value", ylabel = "Density")

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
abstract type AbstractPeriodicOrbitProblem end

# Periodic orbit computations by finite differences
abstract type AbstractPOFDProblem <: AbstractPeriodicOrbitProblem end
# Periodic orbit computations by shooting
abstract type AbstractShootingProblem <: AbstractPeriodicOrbitProblem end

# get the number of time slices
@inline getMeshSize(pb::AbstractPeriodicOrbitProblem) = pb.M

"""
$(SIGNATURES)

Compute the period of the periodic orbit associated to `x`.
"""
@inline getPeriod(sh::AbstractPeriodicOrbitProblem, x, par = nothing) = extractPeriod(x)
@inline extractPeriod(x::AbstractVector) = x[end]
@inline extractPeriod(x::BorderedArray)  = x.p

# update a problem with arguments
function updateForBS(prob::AbstractPeriodicOrbitProblem, F, dF, hopfpt, ζr, M, orbitguess_a, period) end

# update a section with a problem, we do nothing by default
updateSection!(prob::AbstractPeriodicOrbitProblem, x, par) = @warn "Not yet implemented!"

Base.size(pb::AbstractPOFDProblem) = (pb.M, pb.N)
onGpu(pb::AbstractPOFDProblem) = pb.ongpu
hasHessian(pb::AbstractPOFDProblem) = pb.d2F == nothing
isInplace(pb::AbstractPOFDProblem) = pb.isinplace

function applyF(pb::AbstractPOFDProblem, dest, x, p)
	if isInplace(pb)
		pb.F(dest, x, p)
	else
		dest .= pb.F(x, p)
	end
	dest
end

function applyJ(pb::AbstractPOFDProblem, dest, x, p, dx)
	if isInplace(pb)
		pb.J(dest, x, p, dx)
	else
		dest .= apply(pb.J(x, p), dx)
	end
	dest
end

# function to extract trajectories from branch
getPeriodicOrbit(br::AbstractBranchResult, ind::Int) = getPeriodicOrbit(br.functional, br.sol[ind].x, setParam(br, br.sol[ind].p))

"""
$(SIGNATURES)

This function generates an initial guess for the solution of the problem `pb` based on the orbit `t -> orbit(t)` for t ∈ [0,1] and the period `period`.
"""
function generateSolution(pb::AbstractPeriodicOrbitProblem, orbit, period)
	M = getMeshSize(pb)
	orbitguess_a = [orbit(t) for t in LinRange(0, 2pi, M + 1)[1:M]]
	# append period at the end of the initial guess
	orbitguess_v = reduce(vcat, orbitguess_a)
	if pb  isa PoincareShootingProblem
		return vec(orbitguess_v)
	else
		return vcat(vec(orbitguess_v), period) |> vec
	end
end
####################################################################################################
# this struct allows to have a unified interface with Shooting methods in term of plotting
@with_kw_noshow struct SolPeriodicOrbit{Ts, Tu}
	t::Ts
	u::Tu
end
Base.getindex(sol::SolPeriodicOrbit, i...) = getindex(sol.u, i...)
####################################################################################################
"""
$(TYPEDEF)

Define a structure to interface the Jacobian of the periodic orbits Problem with the Floquet computation methods. If we use the same code as for `newton` (see below) but in `continuation`, it is difficult to tell the eigensolver not to use the jacobian but the monodromy matrix instead.

$(TYPEDFIELDS)
"""
struct FloquetWrapper{Tpb, Tjacpb, Torbitguess, Tp}
	pb::Tpb
	jacpb::Tjacpb
	x::Torbitguess
	par::Tp
end
FloquetWrapper(pb, x, par) = FloquetWrapper(pb, dx -> pb(x, par, dx), x, par)

# jacobian evaluation
(shjac::FloquetWrapper)(dx) = apply(shjac.jacpb, dx)

# this is for use with BorderingBLS with checkPrecision = true
apply(shjac::FloquetWrapper, dx) = apply(shjac.jacpb, dx)

# specific linear solver to dispatch
struct FloquetWrapperLS{T} <: AbstractLinearSolver
	linsolver::T
end
# this constructor prevents from having FloquetWrapperLS(FloquetWrapperLS(ls))
FloquetWrapperLS(ls::FloquetWrapperLS) = ls
(ls::FloquetWrapperLS)(J, rhs; kwargs...) = ls.linsolver(J, rhs; kwargs...)
(ls::FloquetWrapperLS)(J::FloquetWrapper, rhs; kwargs...) = ls.linsolver(J.jacpb, rhs; kwargs...)
(ls::FloquetWrapperLS)(J::FloquetWrapper, rhs1, rhs2) = ls.linsolver(J.jacpb, rhs1, rhs2)

# this is for the use of MatrixBLS
LinearAlgebra.hcat(shjac::FloquetWrapper, dR) = hcat(shjac.jacpb, dR)

####################################################################################################
const DocStrjacobianPOSh = """
- `jacobianPO` Specify the choice of the linear algorithm, which must belong to `[:autodiffMF, :MatrixFree, :autodiffDense, :autodiffDenseAnalytical, :FiniteDifferences]`. This is used to select a way of inverting the jacobian dG
    - For `:MatrixFree`, we use an iterative solver (e.g. GMRES) to solve the linear system. The jacobian was specified by the user in `prob`.
    - For `:autodiffMF`, we use iterative solver (e.g. GMRES) to solve the linear system. We use Automatic Differentiation to compute the (matrix-free) derivative of `x -> prob(x, p)`.
    - For `:autodiffDense`. Same as for `:autodiffMF` but the jacobian is formed as a dense Matrix. You can use a direct solver or an iterative one using `options`.
    - For `:FiniteDifferencesDense`, same as for `:autodiffDense` but we use Finite Differences to compute the jacobian of `x -> prob(x, p)` using the `δ = 1e-8` which can be passed as an argument.
    - For `:autodiffDenseAnalytical`. Same as for `:autodiffDense` but the jacobian is using a mix of AD and analytical formula.
    - For `:FiniteDifferences`, use Finite Differences to compute the jacobian of `x -> prob(x, p)` using the `δ = 1e-8` which can be passed as an argument.
"""
##########################
# newton wrapper
function buildJacobian(prob::AbstractShootingProblem, orbitguess, par, jacobianPO = :MatrixFree)
	if jacobianPO == :autodiffDenseAnalytical
		_J = prob(Val(:JacobianMatrix), orbitguess, par)
		jac = (x, p) -> prob(Val(:JacobianMatrixInplace), _J, x, p)
	elseif jacobianPO == :autodiffDense
		jac = (x, p) -> ForwardDiff.jacobian(z -> prob(z, p), x)
	elseif jacobianPO == :autodiffMF
		jac = (x, p) -> (dx -> ForwardDiff.derivative(z -> prob((@. x + z * dx), p), 0))
	elseif jacobianPO == :FiniteDifferencesDense
		jac = (x, p) -> finiteDifferences(z -> prob(z, p), x; δ = 1e-8)
	else
		jac = (x, p) -> (dx -> prob(x, p, dx))
	end
end

"""
$(SIGNATURES)

This is the Newton-Krylov Solver for computing a periodic orbit using (Standard / Poincaré) Shooting method.
Note that the linear solver has to be apropriately set up in `options`.

# Arguments

Similar to [`newton`](@ref) except that `prob` is either a [`ShootingProblem`](@ref) or a [`PoincareShootingProblem`](@ref). These two problems have specific options to be tuned, we refer to their link for more information and to the tutorials.

- `prob` a problem of type `<: AbstractShootingProblem` encoding the shooting functional G.
- `orbitguess` a guess for the periodic orbit. See [`ShootingProblem`](@ref) and See [`PoincareShootingProblem`](@ref) for information regarding the shape of `orbitguess`.
- `par` parameters to be passed to the functional
- `options` same as for the regular [`newton`](@ref) method.

# Optional argument
$DocStrjacobianPOSh

# Output:
- solution
- history of residuals
- flag of convergence
- number of iterations
"""
function newton(prob::AbstractShootingProblem, orbitguess, par, options::NewtonPar;
		jacobianPO = :MatrixFree, kwargs...)
	jac = buildJacobian(prob, orbitguess, par, jacobianPO)
	return newton(prob, jac, orbitguess, par, options; kwargs...)
end

"""
$(SIGNATURES)

This is the deflated Newton-Krylov Solver for computing a periodic orbit using a (Standard / Poincaré) Shooting method.

# Arguments

Similar to [`newton`](@ref) except that `prob` is either a [`ShootingProblem`](@ref) or a [`PoincareShootingProblem`](@ref).

# Optional argument
$DocStrjacobianPOSh

# Output:
- solution
- history of residuals
- flag of convergence
- number of iterations
"""
function newton(
	prob::AbstractShootingProblem,
	orbitguess::vectype, par,
	options::NewtonPar{T,S,E},
	defOp::DeflationOperator{Tp,Tdot,T,vectype};
	jacobianPO = :MatrixFree,
	kwargs...,
) where {T,Tp,Tdot,vectype,S,E}
	jac = buildJacobian(prob, orbitguess, par, jacobianPO)
	return newton(prob, jac, orbitguess, par, options, defOp; kwargs...)
end

####################################################################################################
# Continuation for shooting problems

"""
$(SIGNATURES)

This is the continuation method for computing a periodic orbit using a (Standard / Poincaré) Shooting method.

# Arguments

Similar to [`continuation`](@ref) except that `prob` is either a [`ShootingProblem`](@ref) or a [`PoincareShootingProblem`](@ref). By default, it prints the period of the periodic orbit.

# Optional argument
- `δ = 1e-8` used for finite differences
$DocStrjacobianPOSh
"""
function continuation(
	prob::AbstractShootingProblem,
	orbitguess, par, lens::Lens, contParams::ContinuationPar, linearAlgo::AbstractBorderedLinearSolver;
	jacobianPO = :MatrixFree,
	updateSectionEveryStep = 0,
	δ = 1e-8,
	kwargs...,
)
	@assert jacobianPO in
			(:autodiffMF, :MatrixFree, :autodiffDense, :autodiffDenseAnalytical, :FiniteDifferencesDense, :FiniteDifferences) "This jacobian is not defined. Please chose another one."

	if computeEigenElements(contParams)
		contParams = @set contParams.newtonOptions.eigsolver = FloquetQaD(contParams.newtonOptions.eigsolver)
	end

	options = contParams.newtonOptions

	# change the user provided finalise function by passing prob in its parameters
	_finsol = modifyPOFinalise(prob, kwargs, updateSectionEveryStep)
	_recordsol = modifyPORecord(prob, kwargs, par, lens)
	_plotsol = modifyPOPlot(prob, kwargs)

	if jacobianPO == :autodiffDenseAnalytical
		_J = prob(Val(:JacobianMatrix), orbitguess, par)
		jac = (x, p) -> (prob(Val(:JacobianMatrixInplace), _J, x, p); FloquetWrapper(prob, _J, x, p));
	elseif jacobianPO == :autodiffDense
		jac = (x, p) -> FloquetWrapper(prob, ForwardDiff.jacobian(z -> prob(z, p), x), x, p)
	elseif jacobianPO == :FiniteDifferencesDense
		jac = (x, p) -> FloquetWrapper(prob, finiteDifferences(z -> prob(z, p), x), x, p)
	elseif jacobianPO == :autodiffMF
		jac = (x, p) -> FloquetWrapper(prob, (dx -> ForwardDiff.derivative(z -> prob(x .+ z .* dx, p), 0)), x, p)
	elseif jacobianPO == :FiniteDifferences
		jac = (x, p) -> FloquetWrapper(prob, dx -> (prob(x .+ δ .* dx, p) .- prob(x, p)) ./ δ, x, p)
	else
		jac = (x, p) -> FloquetWrapper(prob, x, p)
	end

	# we have to change the Bordered linearsolver to cope with our type FloquetWrapper
	linearAlgo = @set linearAlgo.solver = FloquetWrapperLS(linearAlgo.solver)

	branch, u, τ = continuation(
		prob, jac,
		orbitguess, par, lens,
		(@set contParams.newtonOptions.linsolver = FloquetWrapperLS(options.linsolver)), linearAlgo;
		kwargs...,
		recordFromSolution = _recordsol,
		finaliseSolution = _finsol,
		plotSolution = _plotsol
		)
	return setproperties(branch; type = :PeriodicOrbit, functional = prob), u, τ
end

"""
$(SIGNATURES)

This is the continuation routine for computing a periodic orbit using a (Standard / Poincaré) Shooting method.

# Arguments

Similar to [`continuation`](@ref) except that `prob` is either a [`ShootingProblem`](@ref) or a [`PoincareShootingProblem`](@ref). By default, it prints the period of the periodic orbit.

# Optional argument
- `δ = 1e-8` used for finite differences
$DocStrjacobianPOSh

"""
function continuation(prob::AbstractPeriodicOrbitProblem, orbitguess, par, lens::Lens, _contParams::ContinuationPar; linearAlgo = nothing, kwargs...)
	_linearAlgo = isnothing(linearAlgo) ?  MatrixBLS() : linearAlgo
	return continuation(prob, orbitguess, par, lens, _contParams, _linearAlgo; kwargs...)
end

####################################################################################################
"""
$(SIGNATURES)

Perform automatic branch switching from a Hopf bifurcation point labelled `ind_bif` in the list of the bifurcated points of a previously computed branch `br::ContResult`. It first computes a Hopf normal form.

# Arguments

- `F, dF, d2F, d3F`: function `(x, p) -> F(x, p)` and its differentials `(x, p, dx) -> d1F(x, p, dx)`, `(x, p, dx1, dx2) -> d2F(x, p, dx1, dx2)`... These are used to compute the Hopf normal form.
- `br` branch result from a call to `continuation`
- `ind_hopf` index of the bifurcation point in `br`
- `contParams` parameters for the call to `continuation`
- `prob` problem used to specify the way the periodic orbit is computed. It can be [`PeriodicOrbitTrapProblem`](@ref), [`ShootingProblem`](@ref) or [`PoincareShootingProblem`](@ref) .

# Optional arguments

- `Jᵗ` is the jacobian adjoint, used for computation of the eigen-elements of the jacobian adjoint, needed to compute the spectral projector for the Hopf normal form. For matrix-free methods, `transpose` is not readily available and the user must provide a dedicated method. In the case of Matrix / Sparse based jacobian, `Jᵗ` should not be passed as it is computed internally more efficiently, i.e. it avoid recomputing the jacobian as it would be if you pass `Jᵗ = (x, p) -> transpose(dF(x, p))`
- `δ = 1e-8` used for finite differences
- `δp` used to specify a particular guess for the parameter on the bifurcated branch which is otherwise determined by `contParams.ds`. This allows to use a step larger than `contParams.dsmax`.
- `ampfactor = 1` factor which alter the amplitude of the bifurcated solution. Useful to magnify the bifurcated solution when the bifurcated branch is very steep.
- `usedeflation = true` whether to use nonlinear deflation (see [Deflated problems](@ref)) to help finding the guess on the bifurcated branch
- `updateSectionEveryStep = 0` updates the section every `updateSectionEveryStep` step during continuation
- `jacobianPO` specify the way the jacobian is computed. See corresponding options (online docs) in the regular call `continuation` for the problem `prob` used to define the periodic orbit.
- all `kwargs` from [`continuation`](@ref)

A modified version of `prob` is passed to `plotSolution` and `finaliseSolution`.

!!! note "Linear solver"
	You have to be careful about the options `contParams.newtonOptions.linsolver`. In the case of Matrix-Free solver, you have to pass the right number of unknowns `N * M + 1`. Note that the options for the preconditioner are not accessible yet.

!!! tip "Jacobian transpose"
	The adjoint of the jacobian `J` is computed internally when `Jᵗ = nothing` by using `transpose(J)` which works fine when `J` is an `AbstractArray`. In this case, do not pass the jacobian adjoint like `Jᵗ = (x, p) -> transpose(d_xF(x, p))` otherwise the jacobian would be computed twice!

!!! warning "Hessian"
	The hessian of `F`, when `d2F` is not `nothing`, is computed with Finite differences.
"""
function continuation(F, dF, d2F, d3F, br::AbstractBranchResult, ind_bif::Int, _contParams::ContinuationPar, prob::AbstractPeriodicOrbitProblem ; Jᵗ = nothing, δ = 1e-8, δp = nothing, ampfactor = 1, usedeflation = false, nev = _contParams.nev, updateSectionEveryStep = 0, kwargs...)
	# compute the normal form of the branch point
	verbose = get(kwargs, :verbosity, 0) > 1 ? true : false
	verbose && (println("--> Considering bifurcation point:"); _show(stdout, br.specialpoint[ind_bif], ind_bif))

	cb = get(kwargs, :callbackN, cbDefault)

	hopfpt = hopfNormalForm(F, dF, d2F, d3F, br, ind_bif; Jᵗ = Jᵗ, δ = δ, nev = nev, verbose = verbose)

	# compute predictor for point on new branch
	ds = isnothing(δp) ? _contParams.ds : δp
	Ty = typeof(ds)
	pred = predictor(hopfpt, ds; verbose = verbose, ampfactor = Ty(ampfactor))

	verbose && printstyled(color = :green, "#"^61*
			"\n--> Start branching from Hopf bif. point to periodic orbits.
			 \n--> Method = \n", prob,
			"\n--> Bifurcation type = ", hopfpt.type,
			"\n----> Hopf param     = ", br.specialpoint[ind_bif].param,
			"\n----> newp           = ", pred.p, ", δp = ", pred.p - br.specialpoint[ind_bif].param,
			"\n----> amplitude      = ", pred.amp,
			"\n----> period         = ", pred.period, "\n")

	# we compute a phase so that the constraint equation
	# < u(0) − u_hopf, ψ > is satisfied, i.e. equal to zero.
	ζr = real.(hopfpt.ζ)
	ζi = imag.(hopfpt.ζ)
	# this phase is for POTrap problem constraint to be satisfied
	ϕ = atan(dot(ζr, ζr), dot(ζi, ζr))
	verbose && printstyled(color = :green, "----> phase ϕ        = ", ϕ / pi, "⋅π\n")

	M = getMeshSize(prob)
	orbitguess_a = [pred.orbit(t - ϕ) for t in LinRange(0, 2pi, M + 1)[1:M]]

	# build the variable to hold the functional for computing PO based on finite differences
	probPO, orbitguess = reMake(prob, F, dF, br.params, hopfpt, ζr, orbitguess_a, abs(2pi/pred.ω); orbit = pred.orbit)

	if _contParams.newtonOptions.linsolver isa GMRESIterativeSolvers
		_contParams = @set _contParams.newtonOptions.linsolver.N = length(orbitguess)
	end

	if usedeflation
		verbose &&
			println("\n--> Attempt branch switching\n--> Compute point on the current branch...")
		prob isa PoincareShootingProblem &&
			@warn "Poincaré Shooting does not work very well with stationary states."
		optn = _contParams.newtonOptions

		# we start with the case of zero amplitude
		orbitzeroamp_a = [hopfpt.x0 for _ = 1:M]
		# this factor prevent shooting jacobian from being singular at fixed points
		# Tfactor = (prob isa AbstractPOFDProblem) || (prob isa PoincareShootingProblem) ? 0 : 0.001
		if probPO isa PoincareShootingProblem
			Tfactor = 0
		elseif prob isa AbstractPOFDProblem
			Tfactor = 100 / abs(2pi / pred.ω)
		else
			Tfactor = 0.001
		end

		_, orbitzeroamp = reMake(prob, F, dF, br.params, hopfpt, ζr, orbitzeroamp_a,	Tfactor * abs(2pi / pred.ω))
		sol0, _, flag, _ = newton(probPO, orbitzeroamp, setParam(br, pred.p), optn; callback = cb, kwargs...)

		# find the bifurcated branch using deflation
		if ~(probPO isa PoincareShootingProblem)
			deflationOp = DeflationOperator(2, (x, y) -> dot(x[1:end-1], y[1:end-1]), 1.0, [sol0]; autodiff = true)
		else
			deflationOp = DeflationOperator(2, (x, y) -> dot(x, y) / M, 1.0, [sol0]; autodiff = true)
		end

		verbose && println("\n--> Compute point on bifurcated branch...")
		# @assert 1==9 "Mauvais dispatch ICI. Va en L199 Deflation.jl"
		solbif, _, flag, _ = newton(
			probPO, orbitguess, setParam(br, pred.p), (@set optn.maxIter = 10 * optn.maxIter), deflationOp; callback = cb, kwargs...)
		@assert flag "Deflated newton did not converge"
		orbitguess .= solbif

		branch, u, τ = continuation(
			probPO, orbitguess,
			setParam(br, pred.p), br.lens,
			_contParams;
			kwargs...,
			updateSectionEveryStep = updateSectionEveryStep,
		)

		return Branch(branch, hopfpt), u, τ
	end

	# perform continuation
	branch, u, τ = continuation(
		probPO, orbitguess, setParam(br, pred.p),
		br.lens, _contParams;
		kwargs..., # put this first to be overwritten just below!
		updateSectionEveryStep = updateSectionEveryStep,
	)
	return Branch(branch, hopfpt), u, τ
end

####################################################################################################
# Branch switching from Bifs of PO
"""
$(SIGNATURES)

Branch switching at a Bifurcation point of a branch of periodic orbits (PO) specified by a `br::AbstractBranchResult`. The functional used to compute the PO is `br.functional`. A deflated Newton-Krylov solver can be used to improve the branch switching capabilities.

# Arguments
- `br` branch of periodic orbits computed with a [`PeriodicOrbitTrapProblem`](@ref)
- `ind_bif` index of the branch point
- `_contParams` parameters to be used by a regular [`continuation`](@ref)

# Optional arguments
- `δp = 0.1` used to specify a particular guess for the parameter in the branch which is otherwise determined by `contParams.ds`. This allows to use a step larger than `contParams.dsmax`.
- `ampfactor = 1` factor which alter the amplitude of the bifurcated solution. Useful to magnify the bifurcated solution when the bifurcated branch is very steep.
- `usedeflation = true` whether to use nonlinear deflation (see [Deflated problems](@ref)) to help finding the guess on the bifurcated branch
- `jacobianPO = :BorderedLU` linear solver used for the Newton-Krylov solver when applied to [`PeriodicOrbitTrapProblem`](@ref).
- `recordFromSolution = (u, p) -> u[end]`, print method used in the bifurcation diagram, by default this prints the period of the periodic orbit.
- `linearAlgo = BorderingBLS()`, same as for [`continuation`](@ref)
- `kwargs` keywords arguments used for a call to the regular [`continuation`](@ref) and the ones specific to periodic orbits (POs).
"""
function continuation(br::AbstractBranchResult, ind_bif::Int, _contParams::ContinuationPar;
			δp = 0.1, ampfactor = 1,
			usedeflation = false,
			linearAlgo = nothing,
			kwargs...)

	bppt = br.specialpoint[ind_bif]
	bptype = bppt.type
	@assert bptype in (:pd, :bp) "Branching from $(bptype) not possible yet."

	# @assert br.functional isa AbstractShootingProblem
	@assert abs(bppt.δ[1]) == 1 "Only simple bifurcation points are handled"

	verbose = get(kwargs, :verbosity, 0) > 0

	verbose && printstyled(color = :green, "#"^61*
			"\n--> Start branching from $(bptype) point to periodic orbits.
			 \n--> Bifurcation type = ", bppt.type,
			"\n----> bif. param     = ", bppt.param,
			"\n----> period at bif. = ", getPeriod(br.functional, bppt.x, setParam(br, bppt.param)),
			"\n----> newp           = ", bppt.param + δp, ", δp = ", δp,
			"\n----> amplitude      = ", ampfactor,
			"\n")

	_linearAlgo = isnothing(linearAlgo) ? BorderingBLS(_contParams.newtonOptions.linsolver) : linearAlgo

	bifpt = br.specialpoint[ind_bif]

	# we copy the problem for not mutating the one passed by the user
	pb = deepcopy(br.functional)

	# let us compute the kernel
	λ = (br.eig[bifpt.idx].eigenvals[bifpt.ind_ev])
	verbose && print("--> computing nullspace of Periodic orbit problem...")
	ζ = geteigenvector(br.contparams.newtonOptions.eigsolver, br.eig[bifpt.idx].eigenvec, bifpt.ind_ev)
	# we normalize it by the sup norm because it could be too small/big in L2 norm
	# TODO: user defined scaleζ
	ζ ./= norm(ζ, Inf)
	verbose && println("Done!")

	# compute the full eigenvector
	ζ_a = MonodromyQaD(Val(:ExtractEigenVector), pb, bifpt.x, setParam(br, bifpt.param), real.(ζ))
	ζs = reduce(vcat, ζ_a)

	## predictor
	pbnew, orbitguess = predictor(pb, bifpt, ampfactor, real.(ζs), bifpt.type)
	newp = bifpt.param + δp

	pbnew(orbitguess, setParam(br, newp))[end] |> abs > 1 && @warn "PO Trap constraint not satisfied"

	# a priori, the following do not overwrite the options in br
	# hence the results / parameters in br are kept intact
	if pb isa AbstractShootingProblem
		if _contParams.newtonOptions.linsolver isa GMRESIterativeSolvers
			@set! _contParams.newtonOptions.linsolver.N = length(orbitguess)
		elseif _contParams.newtonOptions.linsolver isa FloquetWrapperLS
			if _contParams.newtonOptions.linsolver.linsolver isa GMRESIterativeSolvers
				@set! _contParams.newtonOptions.linsolver.linsolver.N = length(orbitguess)
			end
		end
	end

	if usedeflation
		verbose && println("\n--> Attempt branch switching\n--> Compute point on the current branch...")
		optn = _contParams.newtonOptions
		# find point on the first branch
		sol0, _, flag, _ = newton(pbnew, bifpt.x, setParam(br, newp), optn; kwargs...)

		# find the bifurcated branch using deflation
		@assert pbnew isa AbstractPOFDProblem || pbnew isa ShootingProblem "Deflated newton is not available for your problem. Try Trapezoid / collocation method of ShootingProblem"
		deflationOp = DeflationOperator(2, (x, y) -> dot(x[1:end-1], y[1:end-1]), 1.0, [sol0]; autodiff = true)
		verbose && println("\n--> Compute point on bifurcated branch...")
		solbif, _, flag, _ = newton(pbnew, orbitguess, setParam(br, newp),
			(@set optn.maxIter = 10 * optn.maxIter),
			deflationOp; kwargs...,)
		@assert flag "Deflated newton did not converge"
		orbitguess .= solbif
	end

	# perform continuation
	branch, u, τ = continuation( pbnew, orbitguess, setParam(br, newp), br.lens, _contParams;
		kwargs..., # put this first to be overwritten just below!
		linearAlgo = _linearAlgo,
	)

	# create a branch
	bppo = Pitchfork(bifpt.x, bifpt.param, setParam(br, bifpt.param), br.lens, ζ, ζ, nothing, :nothing)

	return Branch(setproperties(branch; type = :PeriodicOrbit, functional = br.functional), bppo), u, τ
end


text:
print("Please enter a number to factorialize: ")

n = parse(Int64, readline())

counter = n
total = 1

while counter > 1
	global total *= counter
	global counter -= 1
end

println("$n ! = $total")


text:
# All said and done, this was not a very sensible way to impute data

using HTTP
using URIs
using JSON

# Some basic constants to do URI
const GEOCODE_SCHEME = "https"
const GEOCODE_HOST = "maps.googleapis.com"
const GEOCODE_PATH = "/maps/api/geocode/json"
const GEOCODE_ADDRESS_KEY = "address"
const GEOCODE_API_KEY = "key"
const GEOCODE_ADDRESS_SUFFIX = "Victoria Australia"

"""
    generate_geocode_uri(address::String, suburb::String, api_key::String)::URI

Construct `URI` based on our constants and given variables
"""
generate_geocode_uri(address::String, suburb::String, api_key::String)::URI =
    URI(scheme=GEOCODE_SCHEME, host=GEOCODE_HOST, path=GEOCODE_PATH, query=Dict([(GEOCODE_ADDRESS_KEY, string(address, " ", suburb, " ", GEOCODE_ADDRESS_SUFFIX)), (GEOCODE_API_KEY, api_key)]))

"""
    get_geocode(address, suburb, api_key)

Send a GET request to Google's Geocoding API based on the given vars and our constants
"""
get_geocode(address, suburb, api_key) = HTTP.get(generate_geocode_uri(address, suburb, api_key))

"""
    get_coordinates(address::String, suburb::String, api_key::String)::Tuple{Float64, Float64}

Takes a street address and a suburb, and returns a `Tuple` (Latitude, Longitude) of that house's coordinates - as per Google's Geocoding API. It assumes the house is in Victoria, Australia.

Note that this is (I believe) synchronous.

# Examples
```jldoctest
julia>get_coordinates("85 studley st", "abbotsford", "ImaginaryAPIKey")
(-37.8015848, 144.99679)
```
"""
function get_coordinates(address::String, suburb::String, api_key::String)::Tuple{Float64, Float64}
    location = get_geocode(address, suburb, api_key).body |> String |> JSON.parse |> x->x["results"][1]["geometry"]["location"]
    return location["lat"], location["lng"]
end


text:
lines(s) = split(s, "\n")

function memoize(f)
  mem = d()
  (args...) -> Base.@get!(mem, args, f(args...))
end

function memoize_debounce(f, expiry = 1)
  mem = d()
  function (args...)
    if haskey(mem, args)
      result, t = mem[args]
      if t + expiry > time()
        mem[args] = result, time()
        return result
      end
    end
    result = f(args...)
    mem[args] = result, time()
    return result
  end
end

function withpath(f, path)
  tls = task_local_storage()
  hassource = haskey(tls, :SOURCE_PATH)
  hassource && (path′ = tls[:SOURCE_PATH])
  tls[:SOURCE_PATH] = path
  try
    return f()
  finally
    hassource ?
      (tls[:SOURCE_PATH] = path′) :
      delete!(tls, :SOURCE_PATH)
  end
end


text:
immutable CachedSwapKey
  index::SwapIndex
  fixing::Date
  tenor::Dates.Period
end

Base.isequal(A::CachedSwapKey, B::CachedSwapKey) = A.index == B.index && A.fixing == B.fixing && A.tenor == B.tenor
# Base.hash(A::CachedSwapKey) = hash(A.index + A.fixing + A.tenor)

type GSR{TermStructureConsistentModelType, P <: StochasticProcess1D, PARAM1 <: Parameter, PARAM2 <: Parameter, T <: YieldTermStructure} <: Gaussian1DModel{TermStructureConsistentModelType}
  lazyMixin::LazyMixin
  modT::TermStructureConsistentModelType
  stateProcess::P
  evaluationDate::Date
  enforcesTodaysHistoricFixings::Bool
  reversion::PARAM1
  sigma::PARAM2
  volatilities::Vector{Quote}
  reversions::Vector{Quote}
  volstepdates::Vector{Date}
  volsteptimes::Vector{Float64}
  volsteptimesArray::Vector{Float64}
  ts::T
  swapCache::Dict{CachedSwapKey, VanillaSwap}
  common::ModelCommon
end

function GSR(ts::YieldTermStructure, volstepdates::Vector{Date}, volatilities::Vector{Float64}, reversion::Float64, T::Float64 = 60.0)
  volatilityQuotes = Quote[Quote(v) for v in volatilities]
  reversions = [Quote(reversion)]

  volsteptimes, volsteptimesArray = update_times(ts, volstepdates)

  sigma = PiecewiseConstantParameter(volsteptimes, NoConstraint())
  if length(reversions) == 1
    reversion = ConstantParameter([reversions[1].value], NoConstraint())
  else
    reversion = PiecewiseConstantParameter(volsteptimes, NoConstraint())
  end

  for i in eachindex(sigma.times)
    set_params!(sigma, i, volatilities[i])
  end

  for i in eachindex(get_data(reversion))
    set_params!(reversion, i, reversions[i].value)
  end

  stateProcess = GsrProcess(volsteptimesArray, sigma.times, get_data(reversion), T)

  swapCache = Dict{CachedSwapKey, VanillaSwap}()

  # TOOD registration with quotes in reversions and volatilities

  return GSR(LazyMixin(), TermStructureConsistentModelType(), stateProcess, Date(), true, reversion, sigma, volatilityQuotes, reversions, volstepdates,
        volsteptimes, volsteptimesArray, ts, swapCache, ModelCommon())
end

function update_times(ts::TermStructure, volstepdates::Vector{Date})
  volsteptimes = zeros(length(volstepdates))
  volsteptimesArray = zeros(length(volsteptimes))
  for i in eachindex(volstepdates)
    volsteptimes[i] = time_from_reference(ts, volstepdates[i])
    volsteptimesArray[i] = volsteptimes[i]
    if i == 1
      volsteptimes[1] > 0.0 || error("volsteptimes must be positive")
    else
      volsteptimes[i] > volsteptimes[i - 1] || error("volsteptimes must be strictly increasing")
    end
  end

  # TODO flush cache for process

  return volsteptimes, volsteptimesArray
end

function update_times!(model::GSR)
  model.volsteptimes, model.volsteptimesArray = update_times(model.ts, model.volstepdates)
  flush_cache!(model.stateProcess)

  return model
end

function update_state!(model::GSR)
  for i in eachindex(get_data(model.sigma))
    set_params!(model.sigma, i, model.volatilities[i].value)
  end
  for i in eachindex(get_data(model.reversion))
    set_params!(model.reversion, i, model.reversions[i].value)
  end

  flush_cache!(model.stateProcess)
end

function generate_arguments!(model::GSR)
  flush_cache!(model.stateProcess)
  notify_observers!(model)

  return model
end

function set_params!(model::GSR, params::Vector{Float64})
  if check_params_equal(model, params)
    return model
  end
  n = length(get_data(model.reversion))
  # first set reversion
  for i = 1:n
    set_params!(model.reversion, i, params[i])
  end
  # now set sigma
  for i = n + 1:length(get_data(model.sigma))+n
    set_params!(model.sigma, i -n, params[i])
  end

  generate_arguments!(model)

  return model
end

# additional methods #
function perform_calculations!(model::GSR)
  if model.evaluationDate == Date() # not calculated
   model.evaluationDate = settings.evaluation_date
  end

  update_times!(model)
  update_state!(model)

  return model
 end

 function get_params(model::GSR)
   # res = Vector{Float64}(length(get_data(model.reversion)) + length(get_data(model.sigma)))
   res = Vector{Float64}()
   append!(res, get_data(model.reversion))
   append!(res, get_data(model.sigma))

   return res
 end

 get_volatilities(model::GSR) = get_data(model.sigma)

 function move_volatility(model::GSR, i::Int)
   res = trues(length(model.reversions) + length(model.volatilities))
   res[length(model.reversions) + i] = false

   return res
 end

 function gaussian_polynomial_integral(model::Gaussian1DModel, a::Float64, b::Float64, c::Float64, d::Float64, e::Float64, y0::Float64, y1::Float64)
   aa = 4.0 * a
   ba = 2.0 * M_SQRT2 * b
   ca = 2.0 * c
   da = M_SQRT2 * d

   x0 = y0 * M_SQRT1_2
   x1 = y1 * M_SQRT1_2

   return (0.125 * (3.0 * aa + 2.0 * ca + 4.0 * e) * erf(x1) - 1.0 / (4.0 * M_SQRTPI) * exp(-x1 * x1) * (2.0 * aa * x1 * x1 * x1 + 3.0 * aa * x1 +
                 2.0 * ba * (x1 * x1 + 1.0) + 2.0 * ca * x1 + 2.0 * da)) - (0.125 * (3.0 * aa + 2.0 * ca + 4.0 * e) * erf(x0) - 1.0 / (4.0 * M_SQRTPI) * exp(-x0 * x0) *
                (2.0 * aa * x0 * x0 * x0 + 3.0 * aa * x0 + 2.0 * ba * (x0 * x0 + 1.0) + 2.0 * ca * x0 + 2.0 * da))
end

function gaussian_shifted_polynomial_integral(model::Gaussian1DModel, a::Float64, b::Float64, c::Float64, d::Float64, e::Float64, h::Float64, x0::Float64, x1::Float64)
   return gaussian_polynomial_integral(model, a, -4.0 * a * h + b, 6.0 * a * h * h - 3.0 * b * h + c, -4 * a * h * h * h + 3.0 * b * h * h - 2.0 * c * h + d,
        a * h * h * h * h - b * h * h * h + c * h * h - d * h + e, x0, x1)
end

function numeraire_impl(model::GSR, t::Float64, y::Float64, yts::YieldTermStructure)
  calculate!(model)

  p = model.stateProcess

  if t == 0.0
    return isa(yts, NullTermStructure) ? discount(model.ts, p.T) : discount(yts, p.T)
  end

  return zerobond(model, p.T, t, y, yts)
end

numeraire(model::Gaussian1DModel, t::Float64, y::Float64, yts::YieldTermStructure) = numeraire_impl(model, t, y, yts)
numeraire(model::Gaussian1DModel, referenceDate::Date, y::Float64, yts::YieldTermStructure) = numeraire(model, time_from_reference(model.ts, referenceDate), y, yts)

function zerobond_impl(model::GSR, T::Float64, t::Float64, y::Float64, yts::YieldTermStructure)
  calculate!(model)

  if t == 0.0
   return isa(yts, NullTermStructure) ? discount(model.ts, T) : discount(yts, T)
  end

  p = model.stateProcess

  # teststDv = std_deviation(p, 0.0, 0.0, t)
  # testEpxt = expectation(p, 0.0, 0.0, t)
  #
  # println(teststDv)
  # println(testEpxt)
  # error("DIE")

  x = y * std_deviation(p, 0.0, 0.0, t) + expectation(p, 0.0, 0.0, t)
  gtT = G!(p, t, T, x)

  d = isa(yts, NullTermStructure) ? discount(model.ts, T) / discount(model.ts, t) : discount(yts, T) / discount(yts, t)

  # println("x: ", x)
  # println("gtT: ", gtT)
  # println("d: ", d)

  return d * exp(-x * gtT - 0.5 * y!(p, t) * gtT * gtT)
end

zerobond(model::Gaussian1DModel, T::Float64, t::Float64, y::Float64, yts::YieldTermStructure) = zerobond_impl(model, T, t, y, yts)
zerobond(model::Gaussian1DModel, maturity::Date, referenceDate::Date, y::Float64, yts::YieldTermStructure) =
        zerobond(model, time_from_reference(yts, maturity), referenceDate != Date() ? time_from_reference(yts, referenceDate) : 0.0, y, yts)

function forward_rate(model::Gaussian1DModel, fixing::Date, referenceDate::Date, y::Float64, iborIdx::IborIndex)
  calculate!(model)

  if fixing <= (model.evaluationDate - Dates.Day(1))
   return fixing(iborIdx, fixing)
  end

  yts = iborIdx.ts

  valueDate = value_date(iborIdx, fixing)
  endDate = advance(iborIdx.tenor.period, iborIdx.fixingCalendar, valueDate, iborIdx.convention)

  dcf = year_fraction(iborIdx.dc, valueDate, endDate)

  return (zerobond(model, valueDate, referenceDate, y, yts) - zerobond(model, endDate, referenceDate, y, yts)) / (dcf * zerobond(model, endDate, referenceDate, y, yts))
end

function y_grid(model::Gaussian1DModel, stdDevs::Float64, gridPoints::Int, T::Float64 = 1.0, t::Float64 = 0.0, y::Float64 = 0.0)
  result = zeros(2 * gridPoints + 1)

  stdDev_0_T = std_deviation(model.stateProcess, 0.0, 0.0, T)
  e_0_T = expectation(model.stateProcess, 0.0, 0.0, T)

  if t < eps()
   stdDev_t_T = stdDev_0_T
   e_t_T = e_0_T
  else
   stdDev_0_t = std_deviation(model.stateProcess, 0.0, 0.0, t)
   stdDev_t_T = std_deviation(model.stateProcess, t, 0.0, T - t)
   e_0_t = expectation(model.stateProcess, 0.0, 0.0, t)
   x_t = y * stdDev_0_t + e_0_t
   e_t_T = expectation(model.stateProcess, t, x_t, T - t)
  end

  h = stdDevs / gridPoints

  for j = -gridPoints:gridPoints
   result[j + gridPoints + 1] = (e_t_T + stdDev_t_T * j * h - e_0_T) / stdDev_0_T
  end

  return result
end

function calibrate_volatilities_iterative!{H <: CalibrationHelper}(model::GSR, helpers::Vector{H}, method::OptimizationMethod, endCriteria::EndCriteria,
                                          constraint::Constraint = PositiveConstraint(), weights::Vector{Float64} = Vector{Float64}())
  for i in eachindex(helpers)
    h = H[helpers[i]]
    calibrate!(model, h, method, endCriteria, constraint, weights, move_volatility(model, i))
  end

  return model
end


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

name = "CPUInfo"
version = v"0.0.20200612"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/pytorch/cpuinfo.git", "63b254577ed77a8004a9be6ac707f3dccc4e1fd9"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd cpuinfo
mkdir build
cd build
cmake \
    -DCMAKE_INSTALL_PREFIX=$prefix \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DCMAKE_BUILD_TYPE=Release \
    -DCMAKE_INSTALL_LIBDIR=$libdir \
    -DCPUINFO_BUILD_UNIT_TESTS=OFF \
    -DCPUINFO_BUILD_MOCK_TESTS=OFF \
    -DCPUINFO_BUILD_BENCHMARKS=OFF \
    -DCPUINFO_LIBRARY_TYPE=shared \
    ..
cmake --build . -- -j $nproc
make install
if [[ $target == *-w64-mingw32 ]]; then
    install -Dvm 755 libcpuinfo.dll "${libdir}/libcpuinfo.dll"
fi
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()
filter!(p -> !(Sys.isapple(p) && arch(p) == "aarch64"), platforms) # aarch64-macos unsupported

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

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

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


text:
const DIR = [SVector( 0, 1), # north
             SVector( 1, 0), # east
             SVector( 0,-1), # south
             SVector(-1, 0)] # west

const PING = 5
const ENGAGE = 6

const ACT_LETTER = ['N', 'E', 'S', 'W', 'P', 'G']

const END_KILL = SubState([-1,-1], [-1,-1], -1, false)


text:
using DataFrames

function expand_grid(iters...)
    vals = Base.Iterators.product(iters...)
    # l = Base.Iterators.take(vals, 1) |> collect
    reduce(vcat, DataFrame([v]) for v in vals)
end

max_stars_repo_name: konkam/GibbsTypePriors

max_issues_repo_name: konkam/GibbsTypePriors

max_forks_repo_name: konkam/GibbsTypePriors


text:
# This file contains types and structures for DEA model
"""
    AbstractDEAModel
An abstract type representing a DEA model.
"""
abstract type AbstractDEAModel end

"""
    nbos(model::Abstract DEAModel)
Return number of observations of a DEA model.
# Examples
```jldoctest
julia> X = [5 13; 16 12; 16 26; 17 15; 18 14; 23 6; 25 10; 27 22; 37 14; 42 25; 5 17];

julia> Y = [12; 14; 25; 26; 8; 9; 27; 30; 31; 26; 12];

julia> deaio = dea(X, Y);

julia> nobs(deaio)
11
```
"""
nobs(model::AbstractDEAModel) = model.n

"""
    ninputs(model::AbstractDEAModel)
Return number of inputs of a DEA model.
# Examples
```jldoctest
julia> X = [5 13; 16 12; 16 26; 17 15; 18 14; 23 6; 25 10; 27 22; 37 14; 42 25; 5 17];

julia> Y = [12; 14; 25; 26; 8; 9; 27; 30; 31; 26; 12];

julia> deaio = dea(X, Y);

julia> ninputs(deaio)
2
```
"""
ninputs(model::AbstractDEAModel) = model.m

"""
    noutputs(model::AbstractDEAModel)
Return number of outputs of a DEA model.
# Examples
```jldoctest
julia> X = [5 13; 16 12; 16 26; 17 15; 18 14; 23 6; 25 10; 27 22; 37 14; 42 25; 5 17];

julia> Y = [12; 14; 25; 26; 8; 9; 27; 30; 31; 26; 12];

julia> deaio = dea(X, Y);

julia> noutputs(deaio)
1
```
"""
noutputs(model::AbstractDEAModel) = model.s

"""
    names(model::AbstractDEAModel)
Return the names of the decision making units (DMUs)
# Examples
```jldoctest
julia> X = [5 13; 16 12; 16 26; 17 15; 18 14; 23 6; 25 10; 27 22; 37 14; 42 25; 5 17];

julia> Y = [12; 14; 25; 26; 8; 9; 27; 30; 31; 26; 12];

julia> dmunames = ["A", "B", "C", "D", "E", "F", "G", "H", "I", "J", "K"]

julia> deaio = dea(X, Y, names = dmunames);

julia> names(deaio)
11-element Vector{String}:
 "A"
 "B"
 "C"
 "D"
 ⋮
 "H"
 "I"
 "J"
 "K"
```
"""
function Base.names(model::AbstractDEAModel)

    xnobs = nobs(model)

    if isdefined(model, :dmunames)

        if model.dmunames === nothing
            # If model have no names, return numeric sequence
            retnames = ["$i" for i in 1:xnobs]
        else
            xnameslength = length(model.dmunames)

            if xnameslength == xnobs
                retnames = model.dmunames
            elseif xnameslength < xnobs
                # If length of names is lower than number of observations, append numbers to match
                @warn("Length of names lower than number of observations")
                retnames = [model.dmunames; ["$i" for i in (xnameslength + 1):xnobs]]
            elseif xnameslength > xnobs
                # If length of names is greater than number of observations, split
                @warn("Length of names greater than number of observations")
                retnames = model.dmunames[1:xnobs]
            end
        end
    else
        # Return numeric sequence
        retnames = ["$i" for i in 1:xnobs]
    end

    return retnames
end