text:
#===============================================================================

    Implementations of callbacks---pieces of code that are run in-between
    MCMC updates. The following callbacks are implemented here:
        ✔   Saving callback: saves the history of the updated chain as a
            CSV file.
        ✔   REPL feedback callback: for printing progress messages to REPL.

===============================================================================#

"""
Supertype of all callbacks. They specify additional actions that can be
performed before or after each MCMC step or at the very end of sampling.
"""
abstract type Callback end

#NOTE by default nothing to be done, but for instance a plotting callback will
# need to set up the canvas
"""
    init!(::Callback, ws::GlobalWorkspace)

Initialization actions for a callback. By default nothing to be done.
"""
function init!(::Callback, ws::GlobalWorkspace) end

function init!(callbacks::Vector{<:Callback}, ws::GlobalWorkspace)
    for c in callbacks
        init!(c, ws)
    end
end

function update_callbacks!(
        callbacks::Vector{<:Callback},
        global_ws::GlobalWorkspace,
        local_wss::Vector{<:LocalWorkspace},
        step,
        flag,
    )
    for callback in callbacks
        check_if_execute(callback, step, flag) && execute!(
            callback, global_ws, local_wss, step, flag
        )
    end
end

"""
    check_if_execute(callback::Callback, step, flag)

Check if `callback` is supposed to be executed at the mcmc step `step` with
pre- post- update flag `flag`.
"""
check_if_execute(callback::Callback, step, flag) = false

"""
    execute!(
        callback::Callback,
        global_ws::GlobalWorkspace,
        local_wss,
        step,
        flag
    )

Perform actions as specified by the `callback`. `local_wss` is a vector of all
local workspaces, `step` is an iterator from the `MCMCSchedule` and `flag` is
either `::PreMCMCStep` or `::PostMCMCStep`.
"""
function execute!(
        callback::Callback,
        global_ws::GlobalWorkspace,
        local_wss,
        step,
        flag
    )
end

"""
    cleanup!(callback::Callback, ws::GlobalWorkspace, step)

The last call to callback, after all MCMC steps and before exiting the function
`run!`.
"""
function cleanup!(callback::Callback, ws::GlobalWorkspace, step) end

#===============================================================================
                            Saving callback
===============================================================================#
"""
    struct SavingCallback <: Callback
        save_at_the_end::Bool
        save_intermediate::Bool
        save_at_iters::Vector{Int64}
        filename::String
    end

Struct for saving the intermediate or final states of the sampled chain to a
hard drive.

    SavingCallback(;
        save_at_the_end=true,
        save_at_iters=[],
        overwrite_at_save=false,
        filename="mcmc_results",
        add_datestamp=false,
        path=".",
    )

The main constructor for `SavingCallback`.

# Arguments
- `save_at_the_end`: indicates whether to save to a file at the end of mcmc
                     sampling
- `save_at_iters`: specifies at which additional intermediate iterations the
                   chain should be saved to a file
- `overwrite`: set to true to overwrite any already existing file that shares
               the name with the one passed to this function
- `filename`: the main stem of the file's name
- `add_datestamps`: will add the date and time at the time of creating the
                    callback to the filename if set to `true`
- `path`: specifies the directory path to save to
"""
struct SavingCallback <: Callback
    save_at_the_end::Bool
    save_intermediate::Bool
    save_at_iters::Vector{Int64}
    filename::String

    function SavingCallback(;
            save_at_the_end=true,
            save_at_iters=[],
            overwrite_at_save=false,
            filename="mcmc_results",
            add_datestamp=false,
            path=".",
        )
        timestamp = (
            add_datestamp ?
            string(
                "_",
                Dates.format(
                    Dates.now(),
                    "yyyy-mm-dd_HH:MM:SS"
                )
            ) :
            ""
        )
        filename = string(filename, timestamp)
        full_filename = (
            !overwrite_at_save ?
            find_available_name(path, filename, "") :
            joinpath(path, string( filename, ".csv" ))
        )

        new(
            save_at_the_end,
            length(save_at_iters)>0,
            collect(save_at_iters),
            full_filename,
        )
    end
end

"""
    find_available_name(path, filename, disambig_num, extension=".csv")

Check if a chosen filename already exists, if so, then tries appending
consecutive numbers to the end of the file until the first one that is not used
yet is found.
"""
function find_available_name(path, filename, disambig_num, extension=".csv")
    proposal_name = joinpath(path, string( filename, disambig_num, extension ))
    next_num = disambig_num == "" ? "1" : string(parse(Int64, disambig_num) + 1)
    isfile(proposal_name) && (
        return find_available_name(path, filename, next_num, extension)
    )
    proposal_name
end

"""
    init!(callback::SavingCallback)

Create a CSV file for writing into.
"""
function init!(callback::SavingCallback, ws::GlobalWorkspace)
    open(callback.filename, "w") do _ end
end

"""
    check_if_execute(callback::SavingCallback, iter)

Return `true` if the mcmc iteration `iter` is one at which an intermediate save
is to be made.
"""
function check_if_execute(callback::SavingCallback, step, ::PreMCMCStep)
    !callback.save_intermediate && return false
    step.mcmciter in callback.save_at_iters && step.pidx == 1 && return true
    false
end

"""
    cleanup!(callback::SavingCallback, ws::GlobalWorkspace, iter)

Save the entire MCMC chain, history of proposals, acceptance history etc.
"""
function cleanup!(callback::SavingCallback, ws::GlobalWorkspace, local_wss, step)
    callback.save_at_the_end && execute!(callback, ws, local_wss, step, nothing)
end

"""
    execute!(sc::SavingCallback, ws::GlobalWorkspace, iter)

Save the chain of accepted states, proposed states and acceptance history to
the disk.
"""
function execute!(sc::SavingCallback, ws::GlobalWorkspace, local_wss, step, ::Any)
    iter_start = find_starting_idx(sc, step)
    open(sc.filename, "a") do f
        for i in iter_start:(step.mcmciter-1)
            data_to_csv(f, ws, local_wss, i)
        end
    end
end

"""
    find_starting_idx(callback::SavingCallback, iter)

Find the last index for which saving was done.
"""
function find_starting_idx(callback::SavingCallback, step)
    !callback.save_intermediate && return 1
    iter_idx = first(searchsorted(callback.save_at_iters, step.mcmciter))
    iter_idx == 1 && return 1
    callback.save_at_iters[iter_idx-1]
end

"""
    data_to_csv(f, ws::GlobalWorkspace, i)

Write data entries to a CSV file.
"""
function data_to_csv(f, ws::GlobalWorkspace, local_wss, i)
    for j in 1:length(ws.sub_ws.state_history[i])
        idx = "$i, $j, "
        θ = join(["$θₖ, " for θₖ in ws.sub_ws.state_history[i][j]])
        θ° = join(["$θₖ, " for θₖ in ws.sub_ws.state_proposal_history[i][j]])
        ll = join(["$ll, " for ll in local_wss[j].sub_ws.ll_history[i]])
        ll° = join(["$ll, " for ll in local_wss[j].sub_ws°.ll_history[i]])
        a_r = join(["$a_r, " for a_r in local_wss[j].acceptance_history[i]])
        write(f, string(idx, "!, ", θ, "!, ", θ°, "!, ", ll, "!, ", ll°, "!,", a_r, "\n"))
    end
end

#===============================================================================
                        REPL progress printing callback
===============================================================================#
"""
    struct REPLCallback <: Callback
        print_every_k_iter::Int64
        show_all_updates::Bool
        basic_info_only::Bool
    end

Struct with instructions for printing progress messages to REPL.

    REPLCallback(;
        print_every_k_iter=100,
        show_all_upates=true,
        basic_info_only=true,
    )

Base constructor that uses named arguments.
"""
struct REPLCallback <: Callback
    print_every_k_iter::Int64
    show_all_updates::Bool
    basic_info_only::Bool

    function REPLCallback(;
            print_every_k_iter=100,
            show_all_upates=true,
            basic_info_only=true,
        )
        new(print_every_k_iter, show_all_upates, basic_info_only)
    end
end

function init!(callback::REPLCallback, ws::GlobalWorkspace)
    println(join(fill("*", 40)))
    println("Initializing an MCMC chain")
    summary(ws; init=true)
    println("* * *")
end

function check_if_execute(callback::REPLCallback, step, ::PostMCMCStep)
    (step.mcmciter % callback.print_every_k_iter == 0) || return false
    (callback.show_all_updates || step.pidx == 1) && return true
    false
end

function execute!(rc::REPLCallback, ws::GlobalWorkspace, local_wss, step, ::Any)
    lws, M = local_wss[step.pidx], step.mcmciter
    println("- - - - - - - - - - -")
    println(M, ".", step.pidx, " ", name_of_update(lws))
    _λ(x) = (length(x)==1 ? x[1] : x)
    _ll = map(x->round.(_λ(x(lws, M)), sigdigits=4), [ll, ll°, llr])
    a_r = accepted(lws, M) ? "✔" : "✗"
    println("\tll: $(_ll[1]), ll°: $(_ll[2]), llr: $(_ll[3]), a/r: $a_r")
    rc.basic_info_only || begin
        θs = map(x->round.(x(lws); sigdigits=4), [state, state°])
        println("\t\tθ : $(θs[1])")
        println("\t $a_r  θ°: $(θs[2])")
    end
end

function cleanup!(callback::REPLCallback, ws::GlobalWorkspace, local_wss, step)
    println("\n\nMCMC sampling has been successful!")
    println("Doing some clean-up and finishing...")
    println("\n⋆ ⋆ ⋆\n")
end


text:
module PyUtils

function input(prompt="")
    print(prompt)
    return readline()
end

chr(code::Int) = Char(code)

ord(c::Char) = convert(Int, c)

int(x::Union{Char, String}) = parse(Int, x)

int(x::Float64) = (x < 0 ? -1 : 1) * floor(Int, abs(x))

bin(n) = string(n, base=2)

hex(n) = string(n, base=16)

end # module

ext: jl


text:
export layer, layer!, layerfirst!, layer_in_rect!, upward_transform, transform_from_to

Base.Broadcast.broadcastable(l::Layer) = Ref(l)

function Base.:+(l::LayerContent, c::Clip)
    l.clip = c
    l
end

function Base.:+(l::LayerContent, o::Opacity)
    l.opacity = o
    l
end

function Base.:+(l::LayerContent, o::Operator)
    l.operator = o
    l
end

function Base.:+(l::LayerContent, a::T) where T <: Attribute
    insert!(l.attrs, a)
    l
end

function layer(;scale=1, rotation=deg(0), translation=O)
    Layer(
        Transform(scale=scale, rotation=rotation, translation=translation),
        Vector{LayerContent}[], nothing, Attributes(), Clip(nothing), Opacity(1), Operator())
end

function layer(t::Transform)
    Layer(t, Vector{LayerContent}[], nothing, Attributes(), Clip(nothing), Opacity(1), Operator())
end

function layer!(parent::Layer; scale=1, rotation=deg(0), translation=O)
    l = layer(scale=scale, rotation=rotation, translation=translation)
    push!(parent, l)
    l
end

function layer!(parent::Layer, t::Transform)
    l = Layer(t, Vector{LayerContent}[], nothing, Attributes(), Clip(nothing), Opacity(1), Operator())
    push!(parent, l)
    l
end

function layerfirst!(parent::Layer, t::Transform)
    l = Layer(t, Vector{LayerContent}[], nothing, Attributes(), Clip(nothing), Opacity(1), Operator())
    pushfirst!(parent, l)
    l
end

function layer!(f::Function, parent::Layer, deps...)
    l = Layer((f, [deps...]), Vector{LayerContent}[], nothing, Attributes(), Clip(nothing), Opacity(1), Operator())
    push!(parent, l)
    l
end

function layer_in_rect!(l::Layer, sr::Shape{Rect}, wh::Symbol, mode::Symbol=:center; margin=0)
    layer!(l, sr) do r
        d = @match wh begin
            :w => r.width
            :h => r.height
            _ => error("Symbol not allowed")
        end

        @match mode begin
            :center => Transform(translation=r.center, scale=(d - margin)/2, rotation=r.angle)
            :norm => Transform(translation=topleft(r) + P(margin, margin), scale=(d - 2*margin), rotation=r.angle)
        end

    end
end

function Base.push!(l::Layer, lc::LayerContent)
    push!(l.content, lc)
    lc.parent = l
end

function Base.pushfirst!(l::Layer, lc::LayerContent)
    pushfirst!(l.content, lc)
    lc.parent = l
end

function Base.insert!(l::Layer, index::Integer, lc::LayerContent)
    insert!(l.content, index, lc)
    lc.parent = l
end

function upward_transform(l::Layer, cc::Cairo.CairoContext)
    if isnothing(l.parent)
        return gettransform!(l, cc)
    else
        return upward_transform(l.parent, cc) * gettransform!(l, cc)
    end
end

function ancestorchain(l::LayerContent)
    chain = LayerContent[]
    current = l
    while true
        if isnothing(current.parent)
            break
        else
            pushfirst!(chain, current.parent)
            current = current.parent
        end
    end
    chain
end

function lowestcommonancestorchainfromto(s1::LayerContent, s2::LayerContent)
    ancestors1 = ancestorchain(s1)
    ancestors2 = ancestorchain(s2)

    length1 = Base.length(ancestors1)
    length2 = Base.length(ancestors2)
    n = min(length1, length2)
    for i in 1:n
        if ancestors1[i] !== ancestors2[i]
            if i == 1
                error("No common ancestor")
            else
                return ancestors1[i:end], ancestors2[i:end]
            end
        end
    end
    # all ancestors so far were the same, so the last one (i) is the common one
    return ancestors1[n+1:end], ancestors2[n+1:end]
end

function gettransform! end

function transform_from_to(from::LayerContent, to::LayerContent, cc::Cairo.CairoContext)
    from_ancestors, to_ancestors = lowestcommonancestorchainfromto(from, to)

    #start with identity
    t = Transform()

    for a in reverse(from_ancestors)
        t = gettransform!(a, cc) * t
    end

    for a in to_ancestors
        t = inverse(gettransform!(a, cc)) * t
    end

    t
end

function gettransform!(l::Layer, cc::Cairo.CairoContext)
    if typeof(l.transform) <: Transform
        l.transform
    else
        closure = l.transform[1]
        deps = l.transform[2]
        solved_deps = []
        for d in deps
            if typeof(d) <: Shape
                solved = solve!(d, cc)
                t = transform_from_to(d, l, cc)
                transformed = t * solved
                push!(solved_deps, transformed)
            else
                # for other arguments that are not shapes
                # other layers not handled yet
                push!(solved_deps, d)
            end
        end
        l.transform = closure(solved_deps...)
    end
end

# function scale!(l::Layer, s::Real)
#     trans = gettransform!(l)
#     tnew = Transform(trans.scale * s, trans.rotation, trans.translation)
#     l.transform = tnew
# end

function spaces(n)
    s = ""
    for _ in 1:n
        s *= " "
    end
    s
end

function Base.show(io::IO, l::Layer, level::Int=0)
    println(io, "$(spaces(2level))Layer $(l.transform)")
    for c in l.content
        Base.show(io, c, level+1)
        print(io, "\n")
    end
end

function Base.show(io::IO, s::Shape{T}, level::Int=0) where T
    print(io, "$(spaces(2level))Shape{$T}")
end

function Base.show(io::IO, s::Shapes{T}, level::Int=0) where T
    print(io, "$(spaces(2level))Shapes{$T}")
end


text:
function RayTrace(project,itarget,m,c,e)
    # Set (gradient, c0 and x0) of ray
    # m0 = 0
    x0 = project.dat["master"].Easting_HA[itarget]
    c0 = project.dat["master"].Northing_HA[itarget]

    # Find x-crossing of ray and each boundary -> m*x + c = m0*x0 + c0
    x = (c0 .- c)./m

    # Set bounds on x-crossing for earch boundary
    limits = Dict(
        "upper" => maximum([e[1:end-1] e[2:end]],dims=2),
        "lower" => minimum([e[1:end-1] e[2:end]],dims=2),
    )

    # Count number of boundary crossings
    count = sum(intersect(
        findall([x[i].<limits["upper"][i] for i in 1:size(limits["upper"],1)]),
        findall([x[i].>maximum([x0 limits["lower"][i]]) for i in 1:size(limits["lower"],1)]),
    ).>0)

    # Set interior/exterior (odd/even) condition
    interior = isodd(count)

    return interior
end

function GeoRef(project,k)
    # Identify features labelled as "Building" in "descr" column, but also
    # ... exclude "Land" and "Roads Tracks And Paths" in "theme col")
    f = project.dat["gml"][k].summary
    iKeep = intersect(
        findall([length(findall(descr.=="Building"))>0 for descr in f.descr]),
        findall([length(findall(theme.=="Land"))==0 for theme in f.theme]),
        findall([length(findall(descr.=="Roads Tracks And Paths"))==0 for descr in f.theme])
    )
    # Remove all geom with length <4 (adiabatic boundaries - no-rigorous method)
    f = project.dat["gml"][k].summary[iKeep,:]
    iKeep = iKeep[findall(length.([eastings[1] for eastings in f.eastings]).>=5)]

    for j in iKeep
        # Eastings/Northings vectors for current feature
        e = project.dat["gml"][k].summary.eastings[j][1]
        n = project.dat["gml"][k].summary.northings[j][1]

        # Gradients and zero-offsets for each boundary (linear eqn)
        m = diff(n)./diff(e)
        c = n[2:end] - m.*e[2:end]

        # HA Eastings/Northings matches (can be more than one)
        iTarget = intersect(
            findall(project.dat["master"].Easting_HA.>=minimum(e)),
            findall(project.dat["master"].Easting_HA.<=maximum(e)),
            findall(project.dat["master"].Northing_HA.>=minimum(n)),
            findall(project.dat["master"].Northing_HA.<=maximum(n)),
        )

        interior_bool = [RayTrace(project,itarget,m,c,e) for itarget in iTarget]

        interior_pts = iTarget[interior_bool]

        # project.dat["master"].osgb_tile[interior_pts] .= project.dat["gml"][k].tile
        project.dat["master"].iGml[interior_pts] .= j
        project.dat["master"].osgb[interior_pts] .= project.dat["gml"][k].summary.osgb[j]
    end

    return project
end


text:
export randEllipic

"""
-  Look for a `phat`, such that when used in NORTA, it will return correlations approximately `ρ`
- `d` : distribution
- `ρ` : target correlation
- NORTA method for generating Random Vectors with Arbitrary Marginal Distribution
    fand CorrelationMatrices first appears in the following paper
    https://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.48.281&rep=rep1&type=pdf
    The generated vector entries usually have a smaller correlation than what is aimed for.
"""
struct preNORTA
    d::D where D<:S
    ρ::Float64
    ρhat::Float64
    function preNORTA(d,ρ,ρhat=ρ)
        L , H = minmax(sign(ρ),0) # Lower bound, Higher bound
        ϵ = min(0.01,(1-abs(ρ))/10)
        while true
            N = [rand(MvNormal([1 ρhat;ρhat 1])) for _=1:10000]
            X = []
            Y = []

            for i in N
               append!(X,quantile(d,cdf(Normal(),i[1])))
               append!(Y,quantile(d,cdf(Normal(),i[2])))
            end

            temp = cor(X,Y)
            if abs(temp-ρ) < ϵ || H - L < ϵ
                break
            end
            if temp < ρ
                ρhat, L = (ρhat+H)/2 ,  (L+ρhat*4)/5
            else
                ρhat, H = (ρhat+L)/2 , (H+ρhat*4)/5
            end
        end
        new(d,ρ,ρhat)
    end
end
"""
```julia
randEllipic(d, n; r , diag, norm)

randEllipic(n; r, norm)
```
- `d` : default `Normal()`, entry distribution
- `n`  : dimensions
- `r` : default `0.5`, the correlation of ``H_{ij},H_{ji}`` pairs
- `norm` : default `false`, if `norm` set to `true`, then the matrix will be normalized with ``n^{-1/2}``.
- `diag` : default `diag=d`, the distribution for diagonal entries.

# Examples
Generate a random elliptic matrix, with entries from ``\\mathscr{N}(0,1)`` and ``\\rho(H_{ij},H_{ji})=0.5``
```julia
randEllipic(500)

500×500 Matrix{Float64}:
  2.03417    -0.424289    1.28267   …  -0.114754  -1.96059
  0.44479    -1.45563     1.32828       1.00149    0.45786
  1.56525     0.0832211  -0.186738     -1.3914     1.04151
 -0.11633    -0.483301   -1.81348      -1.57536    0.514818
  ⋮                                 ⋱
  1.24274    -0.411623   -1.04984      -0.812778  -1.84479
 -0.0817287  -0.254886    0.674914      0.756269  -0.0296209
 -1.48281     0.51675    -1.58041       0.156923   0.244599
  0.852339    1.04593    -0.119082      1.43634    0.114493
```

Generate a normalized random elliptic matrix, with entries `Poisson(10)` and ``\\rho(H_{ij},H_{ji})=0.1``
```julia
using Distributions
randEllipic(Poisson(10),500, r=0.1 , norm=true)

500×500 Matrix{Float64}:
  0.268328    -0.0413153  -0.0175096   …   0.0190835   0.0201304
 -0.00599949   0.447214   -0.0175878       0.0112805   0.100704
 -0.0258879    0.0219927   0.402492       -0.04749    -0.050853
  0.0219071    0.0119609  -0.00448502      0.0043233   0.0404757
  ⋮                                    ⋱
 -0.0145467    0.0800297   0.00247891      0.0189267   0.071565
  0.016412     0.0334019   0.0663348      -0.0180889   0.023773
 -0.0485914   -0.0575288  -0.0409827       0.491935   -0.0969691
  0.0405447    0.0503843   0.00624668      0.0558304   0.402492

```
"""
function randEllipic(d::D, n::Int; r = 0.5::Float64, diag=d::D, norm = false::Bool) where  D<:S
    if r==1
        return randHermitian(d,n,diag=diag,norm=norm)
    end
    # TODO: take care of the r=-1 case.

    r = preNORTA(d,r).ρhat
    M = zeros(n,n)
    for i in 1:n
        M[i,i] = rand(diag)
        for j in i+1:n
            M[i,j],M[j,i] = rand(MvNormal([1 r;r 1]))
        end
    end

    if norm
        M /= sqrt(n)
    end
    return M
end

function randEllipic(n::Int; r=0.5::Float64, norm = false::Bool)
    return randEllipic(Normal(),n, r=r, norm=norm)
end

"""
```julia
randStochastic(n; type, norm)
```
- `n`: dimension
- `type` : default `type = 3`, `3` for doubly randStochastic, `1` for row, `2` for column stochastic
- `norm` : default `false`, if set to `true`, the matrix will be normalized by ``\\sqrt{n}`` (*not a typo*)
# Examples

Generates a 3 by 3 random doubly stochastic  matrix
```julia
randStochastic(3)

3×3 Matrix{Float64}:
 0.132593  0.216041  0.651367
 0.484097  0.320777  0.195126
 0.261495  0.537825  0.20068
```

Generates a 3 by 3 normalized random row stochastic  matrix
```julia
randStochastic(3,type = 1)

3×3 Matrix{Float64}:
 0.220849   0.146942  0.632209
 0.188052   0.26294   0.549008
 0.0170714  0.524574  0.458355
```

Generates a 3 by 3 normalized random column stochastic  matrix
```julia
randStochastic(3,type=2,norm=true)

3×3 Matrix{Float64}:
 0.583396  0.608739  0.732921
 0.672821  0.078786  0.302657
 0.475834  1.04453   0.696473
```
"""
function randStochastic(n::Int; type = 3 ::Int, norm = false::Bool)
    if type == 3
       M = randUnitary(n)
       M = M.*conj(M)
    else
        M = rand(n,n)
        for i in 1:n
            M[i,:] /= sum(M[i,:])
        end
        if type == 2
            M = M'
        end
    end

    if norm
        M*=sqrt(n)
    end

    return M
end


text:
sbrk(x) = Base.llvmcall(     # Not working
    ("declare i32 @xxx(i32)",
"""%ret.i = call i32 @xxx(i32 %0)
ret i32 %ret.i
"""),
    Int32, Tuple{Int32}, x)

using ExportWebAssembly
#sbrk(x) = @extern(:xxx, Int32, Tuple{Int32}, x)
f(x) = Base.llvmcall(("declare i32 @xxx(i32)", "%ret = call i32 @xxx(i32 %0)\nret i32 %ret\n"), Int32, Tuple{Int32}, x)
irgen(f, Tuple{Int32}, optimize = true)
#write_wasm("sbrk.wasm", sbrk, Tuple{Int32})
#wasm2wast("sbrk.wasm")


text:
###
### Re-naming utilities
###

"""
    create_unique_names(
        model::MOI.ModelLike;
        warn::Bool = false,
        replacements::Vector{Function} = Function[]
    )

Rename variables in `model` to ensure that all variables and constraints have
a unique name. In addition, loop through `replacements` and replace names with
`f(name)`.

If `warn`, print a warning if a variable or constraint is renamed.
"""
function create_unique_names(
    model::MOI.ModelLike;
    warn::Bool = false,
    replacements::Vector{Function} = Function[]
)
    create_unique_variable_names(model, warn, replacements)
    create_unique_constraint_names(model, warn, replacements)
    return
end

function _replace(s::String, replacements::Vector{Function})
    for f in replacements
        s = f(s)
    end
    return s
end

function create_unique_constraint_names(
    model::MOI.ModelLike, warn::Bool, replacements::Vector{Function}
)
    original_names = Set{String}()
    for (F, S) in MOI.get(model, MOI.ListOfConstraints())
        for index in MOI.get(model, MOI.ListOfConstraintIndices{F, S}())
            name = MOI.get(model, MOI.ConstraintName(), index)
            push!(original_names, _replace(name, replacements))
        end
    end
    added_names = Set{String}()
    for (F, S) in MOI.get(model, MOI.ListOfConstraints())
        for index in MOI.get(model, MOI.ListOfConstraintIndices{F, S}())
            original_name = MOI.get(model, MOI.ConstraintName(), index)
            new_name = _replace(
                original_name != "" ? original_name : "c$(index.value)",
                replacements
            )
            if new_name in added_names
                # We found a duplicate name! We could just append a string like
                # "_", but we're going to be clever and loop through the
                # integers to name them appropriately. Thus, if we have three
                # constraints named c, we'll end up with variables named c, c_1,
                # and c_2.
                i = 1
                tmp_name = string(new_name, "_", i)
                while tmp_name in added_names || tmp_name in original_names
                    i += 1
                    tmp_name = string(new_name, "_", i)
                end
                new_name = tmp_name
            end
            push!(added_names, new_name)
            if new_name != original_name
                if warn
                    if original_name == ""
                        @warn("Blank name detected for constraint $(index). " *
                              "Renamed to $(new_name).")
                    else
                        @warn("Duplicate name $(original_name) detected for " *
                              "constraint $(index). Renamed to $(new_name).")
                    end
                end
                MOI.set(model, MOI.ConstraintName(), index, new_name)
            end
        end
    end
end

function create_unique_variable_names(
    model::MOI.ModelLike, warn::Bool, replacements::Vector{Function}
)
    variables = MOI.get(model, MOI.ListOfVariableIndices())
    # This is a list of all of the names currently in the model. We're going to
    # use this to make sure we don't rename a variable to a name that already
    # exists.
    original_names = Set{String}([
        _replace(MOI.get(model, MOI.VariableName(), index), replacements)
        for index in variables
    ])
    # This set of going to store all of the names in the model so that we don't
    # add duplicates.
    added_names = Set{String}()
    for index in variables
        original_name = MOI.get(model, MOI.VariableName(), index)
        new_name = _replace(
            original_name != "" ? original_name : "x$(index.value)",
            replacements
        )
        if new_name in added_names
            # We found a duplicate name! We could just append a string like "_",
            # but we're going to be clever and loop through the integers to name
            # them appropriately. Thus, if we have three variables named x,
            # we'll end up with variables named x, x_1, and x_2.
            i = 1
            tmp_name = string(new_name, "_", i)
            while tmp_name in added_names || tmp_name in original_names
                i += 1
                tmp_name = string(new_name, "_", i)
            end
            new_name = tmp_name
        end
        push!(added_names, new_name)
        if new_name != original_name
            if warn
                if original_name == ""
                    @warn("Blank name detected for variable $(index). Renamed to " *
                          "$(new_name).")
                else
                    @warn("Duplicate name $(original_name) detected for variable " *
                          "$(index). Renamed to $(new_name).")
                end
            end
            MOI.set(model, MOI.VariableName(), index, new_name)
        end
    end
end

###
### Compression utilities
###

function error_mode(mode::String)
    throw(ArgumentError("Compressed mode must be \"r\" or \"w\". Got: $mode."))
end

"""
    abstract type AbstractCompressionScheme end

Base type to implement a new compression scheme for MathOptFormat.

To do so, create a concrete subtype (e.g., named after the compression scheme)
and implement:

    extension(::Val{:your_scheme_extension}) = YourScheme()
    compressed_open(f::Function, filename::String, mode::String, ::YourScheme)
"""
abstract type AbstractCompressionScheme end

struct NoCompression <: AbstractCompressionScheme end

extension(::Val) = NoCompression()

function compressed_open(
    f::Function, filename::String, mode::String, ::NoCompression
)
    return Base.open(f, filename, mode)
end

struct Gzip <: AbstractCompressionScheme end

extension(::Val{:gz}) = Gzip()

function compressed_open(
    f::Function, filename::String, mode::String, ::Gzip
)
    if mode == "w"
        return Base.open(f, CodecZlib.GzipCompressorStream, filename, mode)
    elseif mode == "r"
        return Base.open(f, CodecZlib.GzipDecompressorStream, filename, mode)
    end
    error_mode(mode)
end

struct Bzip2 <: AbstractCompressionScheme end

extension(::Val{:bz2}) = Bzip2()

function compressed_open(
    f::Function, filename::String, mode::String, ::Bzip2
)
    if mode == "w"
        return Base.open(f, CodecBzip2.Bzip2CompressorStream, filename, mode)
    elseif mode == "r"
        return Base.open(f, CodecBzip2.Bzip2DecompressorStream, filename, mode)
    end
    error_mode(mode)
end

struct AutomaticCompression <: AbstractCompressionScheme end

function compressed_open(
    f::Function, filename::String, mode::String, ::AutomaticCompression
)
    ext = Symbol(split(filename, ".")[end])
    return compressed_open(f, filename, mode, extension(Val{ext}()))
end


text:
export ALL_solvers

# Valid combinations

ALL_solvers = Function[]

include("ARCSpectral.jl")
push!(ALL_solvers, ARCSpectral)

include("ARCSpectral_abs.jl")
push!(ALL_solvers, ARCSpectral_abs)

include("TRSpectral.jl")
push!(ALL_solvers, TRSpectral)

include("TRSpectral_abs.jl")
push!(ALL_solvers, TRSpectral_abs)

include("ARCLDLt.jl")
push!(ALL_solvers, ARCLDLt)

include("ARCLDLt_Sham_HO.jl")
push!(ALL_solvers, ARCLDLt_HO_Sham)

include("ARCLDLt_HO_vs_Nwt.jl")
push!(ALL_solvers, ARCLDLt_HO_vs_Nwt)

include("ARCLDLt_abs.jl")
push!(ALL_solvers, ARCLDLt_abs)

include("TRLDLt.jl")
push!(ALL_solvers, TRLDLt)

include("TRLDLt_HO_vs_Nwt_lambda.jl")
push!(ALL_solvers, TRLDLt_HO_vs_Nwt_λ)

include("ARCqKOp1.jl")			   #works in 0.7/1.1
push!(ALL_solvers, ARCqKOp1)

include("TRLDLt_abs.jl")
push!(ALL_solvers, TRLDLt_abs)

include("ARCqKOp.jl")
push!(ALL_solvers, ARCqKOp)

include("ARCqKSparse.jl")
push!(ALL_solvers, ARCqKsparse)

include("ARCqKdense.jl")
push!(ALL_solvers, ARCqKdense)

include("TRKOp.jl")
push!(ALL_solvers, TRKOp)

include("TRKsparse.jl")
push!(ALL_solvers, TRKsparse)

include("TRKdense.jl")
push!(ALL_solvers, TRKdense)

include("ST_TROp.jl")
push!(ALL_solvers, ST_TROp)

include("ST_TRsparse.jl")
push!(ALL_solvers, ST_TRsparse)

include("ST_TRdense.jl")
push!(ALL_solvers, ST_TRdense)

include("ST_ARCOp.jl")
push!(ALL_solvers, ST_ARCOp)

include("ST_ARCsparse.jl")
push!(ALL_solvers, ST_ARCsparse)

include("ST_ARCdense.jl")
push!(ALL_solvers, ST_ARCdense)

# include("ARCMA97.jl")
# push!(ALL_solvers, ARCMA97)
#
# include("ARCMA97_abs.jl")
# push!(ALL_solvers, ARCMA97_abs)
#
# include("TRMA97.jl")
# push!(ALL_solvers, TRMA97)
#
# include("TRMA97_abs.jl")
# push!(ALL_solvers, TRMA97_abs)
#
# include("ARCMA57.jl")
# push!(ALL_solvers, ARCMA57)
#
# include("ARCMA57_Sham_vs_Nwt_lambda.jl")
# push!(ALL_solvers, ARCMA57_Sham_vs_Nwt_λ)
#
# include("ARCMA57_abs.jl")
# push!(ALL_solvers, ARCMA57_abs)
#
# include("TRMA57.jl")
# push!(ALL_solvers, TRMA57)
#
# include("TRMA57_Sham_vs_Nwt_lambda.jl")
# push!(ALL_solvers, TRMA57_Sham_vs_Nwt_λ)
#
# include("TRMA57_abs.jl")
# push!(ALL_solvers, TRMA57_abs)


text:
# https://julialang.github.io/Pkg.jl/v1/creating-packages/#Adding-a-build-step-to-the-package.-1

using SHA

cd(@__DIR__)

# https://github.com/JuliaPackaging/BinaryProvider.jl
# using BinaryProvider

# https://github.com/JuliaPackaging/BinDeps.jl
# using BinDeps
#   ⇒ the future is BinaryBuilder and BinaryProvider

# https://juliapackaging.github.io/BinaryBuilder.jl/latest/
# using BinaryBuilder # TODO: much complicated, many instable
# examples
#   https://github.com/staticfloat/IpoptBuilder/blob/master/build_tarballs.jl
#   BinaryBuilder.run_wizard()

# import Compat.Sys

with_samples   = false
build_samples  = false
force_download = false
keep_download  = false

# for examples, see
#  https://github.com/JuliaInterop/CxxWrap.jl/blob/master/deps/build.jl
#  https://github.com/JuliaMath/GSL.jl/blob/master/deps/build.jl
#  https://github.com/JuliaGraphics/Cairo.jl/blob/master/deps/build.jl
#  https://github.com/jonathanBieler/SimpleDirectMediaLayer.jl/blob/master/deps/build.jl

os_dirname = "unknown"
os_so_extension = ""
Sys.islinux()   &&                        (os_dirname = "linux64"; os_so_extension=".so")
Sys.isapple()   &&                        (os_dirname = "osx32"  ; os_so_extension=".dylib")
Sys.iswindows() && Sys.ARCH == :x86_32 && (os_dirname = "win32"  ; os_so_extension=".dll")
Sys.iswindows() && Sys.ARCH == :x86_64 && (os_dirname = "win64"  ; os_so_extension=".dll")
os_julia_include_dir = joinpath(dirname(Sys.BINDIR),"include","julia")

if build_samples && Sys.iswindows()
  # detect Visual studio https://docs.microsoft.com/de-de/cpp/build/building-on-the-command-line?view=vs-2017#developer_command_file_locations
  found = []
  program_files_x86_dir = ENV["ProgramFiles(x86)"]
  for d in readdir(program_files_x86_dir)
    m = match(r"^Microsoft Visual Studio(.*)$",d)
    if m !== nothing
      s = m.captures[1]
      msvc_dir = joinpath(program_files_x86_dir,m.match)
      if s == "" # MSVC 2017 (and later?)
        for version in readdir(msvc_dir)
          for edition in readdir(joinpath(msvc_dir,version))
            p = joinpath(msvc_dir,version,edition,"VC","Auxiliary","Build")
            ispath(p) && push!(found, (joinpath(version,edition),joinpath(p,"vcvars64.bat")))
          end
        end
      else # MSVC 2015 (and previous?)
        p = joinpath(msvc_dir,"VC")
        ispath(p) && push!(found,(s[2:end],joinpath(p,"vcvars64.bat")))
        p = joinpath(msvc_dir,"VC","bin")
        if ispath(p)
          push!(found,(s[2:end],joinpath(p,"vcvars64.bat")))
          for arch in readdir(p)
            push!(found,(s[2:end],joinpath(p,joinpath(p,arch,"vcvars64.bat"))))
          end
        end
      end
    end
  end
  found = filter(x -> isfile(x[2]),found)
  for (version,path) in found
      println("found Microsoft Visual Studio '$version' in '$path'")
  end
  if isempty(found)
    error("No version of Microsoft Visual Studio found")
  else
    version,path = found[1]
    println("using Microsoft Visual Studio '$version' in '$path'")
    p = run(`$path & cmd`)
  end
end

# prefix = Prefix("./deps") # `/` also on windows?
# @assert ~isempty(prefix.path) ""

openvr_version = "1.9.16"
openvr_dirname = "openvr-$openvr_version"
openvr_dir = joinpath(@__DIR__,openvr_dirname)
openvr_sample_dir = joinpath(openvr_dir,"samples")
openvr_bin_dir = joinpath(openvr_dir,"lib",os_dirname)

ispath(openvr_dir) && force_download && rm(openvr_dir; recursive=true, force=true)

# install(openvr_url, openvr_tar_gz_sha256; prefix=prefix, force=true)
if !ispath(openvr_dir)
  if Sys.islinux() || Sys.isapple()
    openvr_url  = "https://github.com/ValveSoftware/openvr/archive/v$openvr_version.tar.gz"
    openvr_tar_gz_sha256 = "8264e8dc20c926bc1828407a588b0755e15314f48a9b96bd1e550ac030174ecb"
    println("downloading $openvr_url")

    openvr_local_tar_gz="openvr-v$openvr_version.tar.gz"
    (!isfile(openvr_local_tar_gz)) && Base.download(openvr_url, openvr_local_tar_gz)

    println("Checksum...")
    openvr_sha2_256_file = bytes2hex(open(openvr_local_tar_gz) do f;  sha2_256(f); end)
    @assert openvr_tar_gz_sha256 == openvr_sha2_256_file "Checksum mismatch! Expected $openvr_tar_gz_sha256, downloaded $openvr_sha2_256_file"

    println("extracting...")
    run(`tar xf $openvr_local_tar_gz`)
    (!keep_download) && rm(openvr_local_tar_gz)
  elseif Sys.iswindows()
    openvr_url  = "https://github.com/ValveSoftware/openvr/archive/v$openvr_version.zip"
    openvr_zip_sha256 = "e7063fae06220345beb43f8d7a634ebde2be48099e5348f3fc4b506db771770c"
    println("downloading $openvr_url")

    openvr_local_zip="openvr-v$openvr_version.zip"
    (!isfile(openvr_local_zip)) && Base.download(openvr_url, openvr_local_zip)

    println("Checksum...")
    openvr_sha2_256_file = bytes2hex(open(openvr_local_zip) do f;  sha2_256(f); end)
    @assert openvr_zip_sha256 == openvr_sha2_256_file "Checksum mismatch! Expected $openvr_zip_sha256, downloaded $openvr_sha2_256_file"

    println("extracting...")
    # try; using ZipFile; catch e; using Pkg; Pkg.add("ZipFile"); using ZipFile; end
    sevenzip_path = joinpath(Sys.BINDIR, "7z.exe")
    !isfile(sevenzip_path) && (sevenzip_path = joinpath(Sys.BINDIR, "../libexec/7z.exe"))
    @assert isfile(sevenzip_path) "7zip not found in $sevenzip_path."
    run(`$(sevenzip_path) x -y $openvr_local_zip`)

    (!keep_download) && rm(openvr_local_zip)
  end
end

for d in [d for d in readdir(openvr_dir) if d ∉ ["lib", "bin", "headers", "LICENSE", with_samples ? "samples" : "", "README.md"] ]
  println("removing " * joinpath(openvr_dir,d))
  rm(joinpath(openvr_dir,d); recursive=true, force=true)
end

if with_samples
  for d in [d for d in readdir(openvr_sample_dir) if d ∉ ["bin", "hellovr_opengl", "shared", "README.md", "CMakeLists.txt", "thirdparty"] ]
    println("removing " * joinpath(openvr_sample_dir,d))
    rm(joinpath(openvr_sample_dir,d); recursive=true, force=true)
  end

  for d in ["driver_sample", "hellovr_vulkan", "helloworldoverlay", "tracked_camera_openvr_sample"]
    println("re-creating " * joinpath(openvr_sample_dir, d, "CMakeLists.txt"))
    mkpath(joinpath(openvr_sample_dir, d))
    touch(joinpath(openvr_sample_dir, d, "CMakeLists.txt"))
  end

  if build_samples
    mkpath(joinpath(openvr_sample_dir, "build"))
    cd(joinpath(openvr_sample_dir, "build"))
    p = run(`cmake ..`)
    @assert p.exitcode == 0 "Cmake has failed"

    p = run(`make`)
    @assert p.exitcode == 0 "make has failed"
  end

  # run(`$openvr_sample_dir/bin/$os_dirname/hellovr_opengl`)
end

# wrapper_src_dir=joinpath(prefix.path,"wrapper")
# wrapper_build_dir=openvr_bin_dir
# cd(wrapper_src_dir)
# cmd_make="make"
#
# using CxxWrap
#
# oldENV=copy(ENV) # save environment
#   ENV["BUILD_DIR"        ] = wrapper_build_dir
#   ENV["SRC_DIR"          ] = "."
#   ENV["JULIA_INCLUDE_DIR"] = os_julia_include_dir
#   ENV["OPENVR_DIR"       ] = openvr_dir
#   ENV["JLCXX_DIR"        ] = CxxWrap.prefix().path
#   ENV["OS_DIR"           ] = os_dirname
#   ENV["OS_SO_EXTENSION"  ] = os_so_extension
#   ENV["CXX"              ] = "gcc"
#   println("building openvr wrapper with the following settings:")
#   println("  BUILD_DIR         = $(ENV["BUILD_DIR"        ])")
#   println("  SRC_DIR           = $(ENV["SRC_DIR"          ])")
#   println("  JULIA_INCLUDE_DIR = $(ENV["JULIA_INCLUDE_DIR"])")
#   println("  OPENVR_DIR        = $(ENV["OPENVR_DIR"       ])")
#   println("  JLCXX_DIR         = $(ENV["JLCXX_DIR"        ])")
#   println("  OS_DIR            = $(ENV["OS_DIR"           ])")
#   println("  OS_SO_EXTENSION   = $(ENV["OS_SO_EXTENSION"  ])")
#   println("  CXX               = $(ENV["CXX"              ])")
#   p = run(`$cmd_make`)
# for kv in oldENV; ENV[kv[1]] = kv[2]; end  # restore environment
# @assert p.exitcode == 0 "make has failed"


text:
##  Copyright 2022 Myrl Marmarelis
##
##  Licensed under the Apache License, Version 2.0 (the "License");
##  you may not use this file except in compliance with the License.
##  You may obtain a copy of the License at
##
##    http://www.apache.org/licenses/LICENSE-2.0
##
##  Unless required by applicable law or agreed to in writing, software
##  distributed under the License is distributed on an "AS IS" BASIS,
##  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
##  See the License for the specific language governing permissions and
##  limitations under the License.

function scramble_latent_utilities(choice::RankedChoice{R},
    ::Val{N}, spread::T)::RankedChoice{R} where {N,R,T}
  # assumes candidates go from 1:n_candidates
  # throw somewhere in the unit interval? or evenly?
  ceiling = T(1)
  draw_utility = () -> Uniform(T(0), ceiling) |> rand |> T
  utilities = @MVector fill(T(-1), N)
  for rank_index in 1:R
    candidate = choice[rank_index]
    if candidate == 0
      continue
    end
    utility = draw_utility()
    utilities[candidate] = utility
    ceiling = utility
  end
  for leftover in 1:N
    if utilities[leftover] >= 0
      continue
    end
    utility = draw_utility()
    utilities[candidate] = utility
  end
  @show utilities
  # now, finally, we add noise. since candidates further down the line
  # are closer together in the latent space, perturbations are a bigger deal for them.
  for candidate in 1:N
    utilities[candidate] += spread * randn(T)
  end
  ranking = @MVector zeros(Int, N)
  sortperm!(ranking, utilities, order=Base.Reverse)
  ranking[SVector{R}(1:R)] # index by a static vector to keep static
end


text:
s1 = """
123456
abcde
ABCDEFG
"""
d1 = Document(TextDocument(uri"untitled:none", s1, 0), false)
@test get_text(d1) == s1
@test get_offset(d1, 0, 4) == 4
@test get_offset(d1, 1, 2) == 9
@test get_line_offsets(get_text_document(d1)) == [0, 7, 13, 21]
@test get_position_from_offset(d1, 1) == (0, 1)
@test get_position_from_offset(d1, 8) == (1, 1)
@test get_position_from_offset(d1, 15) == (2, 2)
@test get_open_in_editor(d1) == false
set_open_in_editor(d1, true)
@test get_open_in_editor(d1) == true
set_open_in_editor(d1, false)
@test get_open_in_editor(d1) == false
@test is_workspace_file(d1) == false

s2 = """
12μ456
abηde
ABCDEFG"""
d2 = Document(TextDocument(uri"untitled:none", s2, 0), true)
@test get_offset(d2, 0, 4) == 5
@test get_offset(d2, 1, 2) == 10
@test get_line_offsets(get_text_document(d2)) == [0, 8, 15]
@test get_position_from_offset(d2, 1) == (0, 1)
@test get_position_from_offset(d2, 9) == (1, 1)
@test get_position_from_offset(d2, 17) == (2, 2)
@test is_workspace_file(d2) == true

set_text_document!(d2, apply_text_edits(get_text_document(d2), [LanguageServer.TextDocumentContentChangeEvent(Range(1), 0, "12")], 1))
@test get_line_offsets(get_text_document(d2)) == [0, 8, 17]

s4 = "1234\r\nabcd"
d4 = Document(TextDocument(uri"untitled:none", s4, 0), false)
@test_broken get_line_offsets(get_text_document(d4)) == [0, 5]

s5 = "1234\nabcd\n"
d5 = Document(TextDocument(uri"untitled:none", s5, 0), false)
@test get_line_offsets(get_text_document(d5)) == [0, 5, 10]

s6 = "\n"
d6 = Document(TextDocument(uri"untitled:none", s6, 0), false)
@test get_line_offsets(get_text_document(d6)) == [0,1]

@testset "applytextdocumentchanges" begin
    doc = LS.Document(TextDocument(uri"file:///example/path/example.jl", "function foo()", 0), false)
    c1 = LS.TextDocumentContentChangeEvent(LS.Range(LS.Position(0, 14), LS.Position(0, 14)),
                                        0, "\n")
    c2 = LS.TextDocumentContentChangeEvent(LS.Range(LS.Position(1, 0), LS.Position(1, 0)),
                                           0, "    ")
    c3 = LS.TextDocumentContentChangeEvent(missing, missing, "println(\"Hello World\")")

    LS.applytextdocumentchanges(doc, c1)
    @test LS.get_text(doc) == "function foo()\n"
    # Implicitly test for issue #403
    LS.applytextdocumentchanges(doc, c2)
    @test LS.get_text(doc) == "function foo()\n    "
    LS.applytextdocumentchanges(doc, c3)
    @test LS.get_text(doc) == "println(\"Hello World\")"
    # doc currently has only one line, applying change to 2nd line should throw
    @test_throws LanguageServer.LSOffsetError LS.applytextdocumentchanges(doc, c2)
end

@testset "UTF16 handling" begin
    doc = LanguageServer.Document(TextDocument(uri"", "aaa", 0), false)
    @test sizeof(LanguageServer.get_text(doc)) == 3
    @test LanguageServer.get_offset(doc, 0, 0) == 0
    @test LanguageServer.get_position_from_offset(doc, 0) == (0, 0)
    @test LanguageServer.get_offset(doc, 0, 1) == 1
    @test LanguageServer.get_position_from_offset(doc, 1) == (0, 1)
    @test LanguageServer.get_offset(doc, 0, 2) == 2
    @test LanguageServer.get_position_from_offset(doc, 2) == (0, 2)
    @test LanguageServer.get_offset(doc, 0, 3) == 3
    @test LanguageServer.get_position_from_offset(doc, 3) == (0, 3)

    doc = LanguageServer.Document(TextDocument(uri"", "ααα", 0), false)
    @test sizeof(LanguageServer.get_text(doc)) == 6
    @test LanguageServer.get_offset(doc, 0, 0) == 0
    @test LanguageServer.get_position_from_offset(doc, 0) == (0, 0)
    @test LanguageServer.get_offset(doc, 0, 1) == 1
    @test LanguageServer.get_position_from_offset(doc, 1) == (0, 1)
    @test LanguageServer.get_offset(doc, 0, 2) == 3
    @test LanguageServer.get_position_from_offset(doc, 3) == (0, 2)
    @test LanguageServer.get_offset(doc, 0, 3) == 5
    @test LanguageServer.get_position_from_offset(doc, 5) == (0, 3)

    doc = LanguageServer.Document(TextDocument(uri"", "ࠀࠀࠀ", 0), false) # 0x0800
    @test sizeof(LanguageServer.get_text(doc)) == 9
    @test LanguageServer.get_offset(doc, 0, 0) == 0
    @test LanguageServer.get_position_from_offset(doc, 0) == (0, 0)
    @test LanguageServer.get_offset(doc, 0, 1) == 1
    @test LanguageServer.get_position_from_offset(doc, 1) == (0, 1)
    @test LanguageServer.get_offset(doc, 0, 2) == 4
    @test LanguageServer.get_position_from_offset(doc, 4) == (0, 2)
    @test LanguageServer.get_offset(doc, 0, 3) == 7
    @test LanguageServer.get_position_from_offset(doc, 7) == (0, 3)

    doc = LanguageServer.Document(TextDocument(uri"", "𐐀𐐀𐐀", 0), false)
    @test sizeof(LanguageServer.get_text(doc)) == 12
    @test LanguageServer.get_offset(doc, 0, 0) == 0
    @test LanguageServer.get_position_from_offset(doc, 0) == (0, 0)
    @test LanguageServer.get_offset(doc, 0, 2) == 1
    @test LanguageServer.get_position_from_offset(doc, 1) == (0, 2)
    @test LanguageServer.get_offset(doc, 0, 4) == 5
    @test LanguageServer.get_position_from_offset(doc, 5) == (0, 4)
    @test LanguageServer.get_offset(doc, 0, 6) == 9
    @test LanguageServer.get_position_from_offset(doc, 9) == (0, 6)

    doc = LanguageServer.Document(TextDocument(uri"", "𐀀𐀀𐀀", 0), false) # 0x010000
    @test sizeof(LanguageServer.get_text(doc)) == 12
    @test LanguageServer.get_offset(doc, 0, 0) == 0
    @test LanguageServer.get_position_from_offset(doc, 0) == (0, 0)
    @test LanguageServer.get_offset(doc, 0, 2) == 1
    @test LanguageServer.get_position_from_offset(doc, 1) == (0, 2)
    @test LanguageServer.get_offset(doc, 0, 4) == 5
    @test LanguageServer.get_position_from_offset(doc, 5) == (0, 4)
    @test LanguageServer.get_offset(doc, 0, 6) == 9
    @test LanguageServer.get_position_from_offset(doc, 9) == (0, 6)
end

@testset "document link provider" begin
    doc = LS.Document(TextDocument(LS.filepath2uri(@__FILE__), """
    include("test_document.jl")
    include("runtests_does_not_exist.jl")
    """, 0), false)
    links = LS.DocumentLink[]
    LS.find_document_links(LS.getcst(doc), doc, 0, links)
    @test length(links) == 1
    @test links[1].target == LS.filepath2uri(@__FILE__)
end


text:
using ADCME
using PyCall
using LinearAlgebra
using PyPlot
using Random
using FwiFlow
# Random.seed!(233)

function saturation_nn(s0,dporodt,pt,perm,poro,qw,qo,muw,muo,sref,thetaw,configw,thetao,configo,dt,h)
    saturation_nn_ = load_op_and_grad("./build/libSaturationNn","saturation_nn")
    s0,dporodt,pt,perm,poro,qw,qo,muw,muo,sref,thetaw,configw,thetao,configo,dt,h = convert_to_tensor(Any[s0,dporodt,pt,perm,poro,qw,qo,muw,muo,sref,thetaw,configw,thetao,configo,dt,h], [Float64,Float64,Float64,Float64,Float64,Float64,Float64,Float64,Float64,Float64,Float64,Int64,Float64,Int64,Float64,Float64])
    saturation_nn_(s0,dporodt,pt,perm,poro,qw,qo,muw,muo,sref,thetaw,configw,thetao,configo,dt,h)
end

# TODO: specify your input parameters
m = 100
n = 10
h = 0.01
x = zeros(n, m)
y = zeros(n, m)
for i = 1:m
    for j = 1:n
        x[j, i] = h*i
        y[j, i] = h*j
    end
end
t = 3.0
s0 = @. (x^2 + y^2)/(1+x^2+y^2) * exp(-t)
dporodt = exp(t) * ones(n, m)
pt = @. (x^2+y^2)
perm = rand(n, m)
poro = exp(t) * ones(n, m)
qw = ones(n, m)
qo = ones(n, m)
muw = 2.0
muo = 3.0
sref = s0
dt = 0.01
h = 0.1
u3 = sat_op(s0,pt,perm,poro,qw,qo,muw, muo,sref,dt,h)

configw = [1,20,20,20,1]
configo = [1,20,20,20,1]
thetaw = ae_init(configw)
thetao = ae_init(configo)

u = saturation_nn(s0,dporodt,pt,perm,poro,qw,qo,muw,muo,sref,thetaw,configw,thetao,configo,dt,h)
sess = Session(); init(sess)

@show run(sess, u3-u)

# uncomment it for testing gradients
# error()

# TODO: change your test parameter to `m`
#       in the case of `multiple=true`, you also need to specify which component you are testings
# gradient check -- v
function scalar_function(m)
    return sum(saturation_nn(s0,m,pt,perm,poro,qw,qo,muw,muo,sref,thetaw,configw,thetao,configo,dt,h)^3)
    # sum(sat_op(m,pt,perm,poro,qw,qo,muw, muo,sref,dt,h)^2)
end

# TODO: change `m_` and `v_` to appropriate values
m_ = constant(dporodt)
v_ = rand(n, m)
y_ = scalar_function(m_)
dy_ = gradients(y_, m_)
ms_ = Array{Any}(undef, 5)
ys_ = Array{Any}(undef, 5)
s_ = Array{Any}(undef, 5)
w_ = Array{Any}(undef, 5)
gs_ =  @. 10000 / 10^(1:5)

for i = 1:5
    g_ = gs_[i]
    ms_[i] = m_ + g_*v_
    ys_[i] = scalar_function(ms_[i])
    s_[i] = ys_[i] - y_
    w_[i] = s_[i] - g_*sum(v_.*dy_)
end

sess = Session(); init(sess)
sval_ = run(sess, s_)
wval_ = run(sess, w_)
close("all")
loglog(gs_, abs.(sval_), "*-", label="finite difference")
loglog(gs_, abs.(wval_), "+-", label="automatic differentiation")
loglog(gs_, gs_.^2 * 0.5*abs(wval_[1])/gs_[1]^2, "--",label="\$\\mathcal{O}(\\gamma^2)\$")
loglog(gs_, gs_ * 0.5*abs(sval_[1])/gs_[1], "--",label="\$\\mathcal{O}(\\gamma)\$")

plt.gca().invert_xaxis()
legend()
xlabel("\$\\gamma\$")
ylabel("Error")


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule FLANN_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("FLANN")
JLLWrappers.@generate_main_file("FLANN", UUID("48b6455b-4cf5-590d-a543-2d733c79e793"))
end  # module FLANN_jll


text:
# should be constriants not projection, projection is a function on the contrsrinats
abstract type Constraints end
struct BoxConstraint{T} <: Constraints
    lower::T
    upper::T
end
abstract type Scaling end
struct NoScaling <: Scaling end
scale(S::NoScaling, gx) = gx
struct ConstantPDScaling{T} <: Scaling
    S::T
end
scale(S::ConstantPDScaling, gx) = S.S\gx

function update_xλ(prob::OptimizationProblem, S::Scaling, λ, x, gx)
    bounds = bounds(prob)
    min.(max.(bounds.lower, x .- λ.*scale(S, gx)), bounds.upper)
end
#update_xλ(p::BoxConstraint, λ, x, gx) = min.(max.(p.lower, x .- λ.*gx), p.upper)

# 5.13 in Kelley IMO
# what's it in Bertsekas?
function sufficient_decrease(fλ, fx, λ, α, x, xλ, gx)
    lhs = fλ
    # We use (22) from Bersekas (1976) to avoid dividing by λ
    rhs = fx - dot(gx, x.-xλ)*α
    lhs <= rhs
end
function sufficient_decrease(fλ, fx, λ, α, gx)
    lhs = fλ
    rhs = fx - norm(gx, 2)*α*λ
    lhs <= rhs
end

function projected_gradient(prob::OptimizationProblem, f, g, x0, sc::Union{Scaling, NoScaling}=NoScaling(); itermax, ls_itermax, β=0.5, α=1e-6, print_trace=true)

    f0, g0 = f(x0), g(x0)
    if norm(x0.-update_xλ(prob, sc, 1, x0, g0), Inf) <= 1e-6
      return x0, f0, g0
    end

    x = x0

	# 5.12 in Kelley IMO
	# what's it in Bertsekas?
	local fx, gx
    for i = 1:itermax
        fx, gx = f(x), g(x)
        if print_trace
	        println("i = $i    fx: $(fx)    |gx|: $(norm(gx,Inf))")
            println("x1 = $(x[1])")
            println("x2 = $(x[2])")
	    end
        if norm(x.-update_xλ(prob, sc, 1, x, gx)) <= 1e-6
            return x, fx, gx
        end
        for j = 1:ls_itermax
        	λ = β^j
        	xλ = update_xλ(P, sc, λ, x, gx)
            fλ = f(xλ)
            if sufficient_decrease(fλ, fx, λ, α, x, xλ, gx)
            	x = xλ
            	break
            end
        end
    end
    x, fx, gx
end


text:
"""
    struct FrontendOptions <: Any
Hold a pointer to a `clang::FrontendOptions` object.
"""
struct FrontendOptions
    ptr::CXFrontendOptions
end

Base.unsafe_convert(::Type{CXFrontendOptions}, x::FrontendOptions) = x.ptr
Base.cconvert(::Type{CXFrontendOptions}, x::FrontendOptions) = x


text:
using LinearAlgebra

function create_board(nrow = 3, ncol = 3)

    board = Array{String}(undef, nrow, ncol)
    board[:] .= ""

    return board

end

function whos_turn(board)

    if sum(board .!= "") == 0
        turn = "X"
    elseif sum(board .!= "") % 2 == 0
        turn = "X"
    else
        turn = "O"
    end

    return turn

end

function allsame(x)

    return all(y -> y == first(x), x)

end

function checkwinner(board)

    n, m = size(board)
    winner = ""

    if sum(board .== "") == 0
        winner = "Tie"
    end

    # horizontal
    for i in 1:n
        if allsame(board[i, :]) && board[i, 1] != ""
            winner = board[i, 1]
        end
    end
    # vertical
    for j in 1:m
        if allsame(board[:, j]) && board[1, j] != ""
            winner = board[1, j]
        end
    end
    # diagonal
    straightdiag = board[diagind(board)]
    if allsame(straightdiag) && straightdiag[1] != ""
        winner = straightdiag[1]
    end
    revboard = board[:, end:-1:1]
    revdiag = revboard[diagind(revboard)]
    if allsame(revdiag) && revdiag[1] != ""
        winner = revdiag[1]
    end

    return winner
end

function isgameover(board)

    over = false
    if sum(board .== "") == 0
        over = true
    elseif checkwinner(board) != ""
        over = true
    end

    return over

end

function parseinput(input)

    #==
    parse input properly
    ==#

    inputs = split(input, ",")
    inputs_array = Tuple([parse(Int, a) for a in inputs])

    return inputs_array

end

function islegal(boardspot)

    legal = boardspot == ""

    return legal

end

function legal_moves(board)

    legalmoves = Tuple.(findall(i->(i==""), board))

    return legalmoves

end

function play_legal_move!(board)

    movelegal = false
    while movelegal == false

        player = whos_turn(board)
        println("Enter input:")
        input = readline()
        input = parseinput(input)

        if input in legal_moves(board)
            row, col = input
            board[row, col] = player
            movelegal = true
        else
            println("Made illegal move. Please try again")
        end
    end

end

function playgame(nrow = 3, ncol = 3)

    board = create_board(nrow, ncol)
    nmaxturns = nrow * ncol
    gameoutcome = "Tie"

    # start game
    for i in 1:nmaxturns

        println("==========================")
        println("Turn number: $i")
        writedlm(stdout, board,  "\t|\t")
        play_legal_move!(board)
        gameover = isgameover(board)
        println(gameover)
        if gameover
            winner = checkwinner(board)

            if winner == ""
                gameoutcome = "Tie!"
            else
                gameoutcome = "$winner wins!"
            end

            break
        end

    end

    println("==========================")
    writedlm(stdout, board,  "\t|\t")
    return gameoutcome

end

max_stars_repo_name: HariharanJayashankar/tictactoe

max_issues_repo_name: HariharanJayashankar/tictactoe

max_forks_repo_name: HariharanJayashankar/tictactoe


text:
"""
Our home-grown SVG library.

This is very janky. Basically, we have a tree structure (`SVGNode`) which
roughly corresponds to SVG (i.e., each node has attributes and children),
we have a macro (`@svg`) which converts a subset of Julia syntax into that
tree structure, and finally we overload Base.show in order to print out
`SVGNode` into actual SVG XML.

This is used to customize the look of particular semagrams. Currently, all we can customize
is the shape of boxes, but eventually we should be able to customize ports, wires, etc.

Note: Is this necessary? Should we just import SVG generated by other tools?
Are we really ever going to want to programmatically generate SVG?
"""
module SVG
export @svg, SVGNode, write_svg, svgpointstring

using MLStyle

"""
The basic data structure for storing SVG.
"""
struct SVGNode
  tag::Symbol
  attrs::Dict{Symbol,String}
  children::Vector{Union{SVGNode, String}}

  """
  Note that attributes/children can be passed in as non-strings,
  and they get automatically converted to strings.
  """
  function SVGNode(tag::Symbol, attrs::Dict,
                   children::Vector)
    new(tag,
        Dict{Symbol,String}(Symbol(k) => string(v) for (k,v) in attrs),
        to_svgnode.(children))
  end
end

function to_svgnode(x)
  string(x)
end

function to_svgnode(x::SVGNode)
  x
end

"""
This prints an SVG to a display which supports SVG (i.e., a Jupyter notebook)

If the SVG passed in is a fragment (i.e., not a full SVG), then it wraps it
with SVG tag before displaying it.

TODO: this should probably wrap with the standard attributes (i.e. xmlns)
instead of a bare SVG tag.
"""
function Base.show(io::IO, mime::MIME"image/svg+xml", svg::SVGNode)
  if svg.tag == :svg
    write_svg(io, svg)
  else
    write_svg(io, SVGNode(:svg, Dict{Symbol,Any}(), [svg]))
  end
end

function write_svg(io::IO, s::String)
  print(io, s)
end

function write_svg(io::IO, svg::SVGNode)
  print(io,"<$(string(svg.tag))")
  for (key,value) in svg.attrs
    print(io, " $(string(key))=\"$(value)\"")
  end
  print(io, ">")
  for child in svg.children
    write_svg(io, child)
  end
  print(io, "</$(string(svg.tag))>")
end

"""
Converts an SVGNode to the SVG XML syntax
"""
function write_svg(svg::SVGNode)
  sprint(write_svg, svg)
end

"""
Converts Julia syntax to a SVGNode

Example:

```julia
@svg g() do
  circle(r=30)
  circle(r=40)
end
```

For more examples of usage, see Logo.jl
"""
macro svg(block)
  esc(make_svg(block))
end

function strip_lines(expr::Expr; recurse::Bool=false)::Expr
  args = [ x for x in expr.args if !isa(x, LineNumberNode) ]
  if recurse
    args = [ isa(x, Expr) ? strip_lines(x; recurse=true) : x for x in args ]
  end
  Expr(expr.head, args...)
end

function make_svg(block)
  block = strip_lines(block, recurse=true)
  @match block begin
    Expr(:block, lines...) =>
      svg_node(Expr(:do, Expr(:call, :svg), Expr(:(->), (), Expr(:block, lines...))))
    Expr(:do, _args...) => svg_node(block)
    Expr(:call, _args...) => svg_node(block)
    _ => error("Expected a block")
  end
end

function svg_node(stx::Expr)
  (head, body) = @match stx begin
    Expr(:do, head, Expr(:(->), _, Expr(:block, lines...))) => (head, [lines...])
    _ => (stx,Expr[])
  end
  @match head begin
    Expr(:call, tag, args...) =>
      Expr(:call,
           SVGNode,
           Expr(:quote, tag),
           Expr(:call, Dict{Any,Any},
                map((e) ->
                    Expr(:call, :(=>),
                         Expr(:quote, e.args[1]),
                         e.args[2]),
                    args)...),
           Expr(:vect, map(svg_node, body)...))
    Expr(:$, expr) => expr
    _ => error("expected a function call or an escape")
  end
end

svg_node(x) = x

# Utility functions

svgpointstring(X) = join([join(X[i,:],",") for i in 1:(size(X)[1])]," ")

end


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    precompile(CUDArt.CUDArt_gen.checkerror, (UInt32,))
    precompile(CUDArt.init, (Int64,))
    precompile(CUDArt.init, (Array{Int64, 1},))
    precompile(CUDArt.devices, (Array{Any, 1}, Function, Function,))
    precompile(CUDArt.device_reset, (Int64,))
    precompile(CUDArt.devices, (Array{Any, 1}, Function,))
    precompile(CUDArt.launch, (Array{Any, 1}, CUDArt.CuFunction, Int64, Tuple{Int64, Int64}, Tuple{CUDArt.CudaPitchedArray{Float64, 2}, Int64, Int64, Int64, UInt64, Float64},))
    precompile(CUDArt.launch, (Array{Any, 1}, CUDArt.CuFunction, Int64, Int64, Tuple{CUDArt.CudaArray{Float64, 2}, Int64, Float64},))
    precompile(CUDArt.launch, (Array{Any, 1}, CUDArt.CuFunction, Int64, Int64, Tuple{CUDArt.CudaPitchedArray{Int64, 1}, Int64, Int64, Int64, UInt64, Int64},))
    precompile(CUDArt.launch, (Array{Any, 1}, CUDArt.CuFunction, Int64, Int64, Tuple{CUDArt.CudaArray{Int32, 2}, Int64, Int32},))
    precompile(CUDArt.launch, (Array{Any, 1}, CUDArt.CuFunction, Int64, Tuple{Int64, Int64}, Tuple{CUDArt.CudaPitchedArray{Int32, 2}, Int64, Int64, Int64, UInt64, Int32},))
    precompile(CUDArt.launch, (Array{Any, 1}, CUDArt.CuFunction, Int64, Int64, Tuple{CUDArt.CudaArray{Int64, 1}, Int64, Int64},))
    precompile(CUDArt.call, (Type{CUDArt.CudaPitchedArray}, Type{Float64}, Tuple{Int64, Int64},))
    precompile(CUDArt.call, (Type{CUDArt.CudaPitchedArray}, Type{Int64}, Tuple{Int64},))
    precompile(CUDArt.call, (Type{CUDArt.CudaPitchedArray}, Type{UInt16}, Tuple{Int64, Int64},))
    precompile(CUDArt.call, (Type{CUDArt.CudaPitchedArray}, Type{Int32}, Tuple{Int64, Int64},))
    precompile(CUDArt.copy!, (Array{Any, 1}, CUDArt.CudaArray{Int64, 1}, Array{Int64, 1},))
    precompile(CUDArt.cudacopy!, (Array{Any, 1}, CUDArt.CudaPitchedArray{UInt16, 2}, Tuple{Base.UnitRange{Int64}, Base.UnitRange{Int64}}, Array{UInt16, 2}, Tuple{Base.UnitRange{Int64}, Base.UnitRange{Int64}},))
    precompile(CUDArt.cudacopy!, (Array{Any, 1}, CUDArt.CudaPitchedArray{Float64, 2}, Tuple{Base.UnitRange{Int64}, Base.UnitRange{Int64}}, Array{Float64, 2}, Tuple{Base.UnitRange{Int64}, Base.UnitRange{Int64}},))
    precompile(CUDArt.call, (Type{CUDArt.CudaPitchedArray}, Type{Float32}, Tuple{Int64},))
    precompile(CUDArt.copy!, (Array{Any, 1}, CUDArt.CudaArray{Float64, 2}, Array{Float64, 2},))
    precompile(CUDArt.cudacopy!, (Array{Any, 1}, Array{Float64, 2}, Tuple{Base.UnitRange{Int64}, Base.UnitRange{Int64}}, CUDArt.CudaPitchedArray{Float64, 2}, Tuple{Base.UnitRange{Int64}, Base.UnitRange{Int64}},))
    precompile(CUDArt.cudacopy!, (Array{Any, 1}, CUDArt.CudaPitchedArray{Int32, 2}, Tuple{Base.UnitRange{Int64}, Base.UnitRange{Int64}}, Array{Int32, 2}, Tuple{Base.UnitRange{Int64}, Base.UnitRange{Int64}},))
    precompile(CUDArt.cudacopy!, (Array{Any, 1}, CUDArt.CudaPitchedArray{Int64, 1}, Tuple{Base.UnitRange{Int64}}, Array{Int64, 1}, Tuple{Base.UnitRange{Int64}},))
    precompile(CUDArt.cudacopy!, (Array{Any, 1}, Array{UInt16, 2}, Tuple{Base.UnitRange{Int64}, Base.UnitRange{Int64}}, CUDArt.CudaPitchedArray{UInt16, 2}, Tuple{Base.UnitRange{Int64}, Base.UnitRange{Int64}},))
    precompile(CUDArt.cudacopy!, (Array{Any, 1}, Array{Int64, 1}, Tuple{Base.UnitRange{Int64}}, CUDArt.CudaPitchedArray{Int64, 1}, Tuple{Base.UnitRange{Int64}},))
    precompile(CUDArt.cudacopy!, (Array{Any, 1}, CUDArt.CudaArray{UInt16, 2}, Tuple{Base.UnitRange{Int64}, Base.UnitRange{Int64}}, CUDArt.CudaPitchedArray{UInt16, 2}, Tuple{Base.UnitRange{Int64}, Base.UnitRange{Int64}},))
    precompile(CUDArt.cudacopy!, (Array{Any, 1}, Array{Int32, 2}, Tuple{Base.UnitRange{Int64}, Base.UnitRange{Int64}}, CUDArt.CudaPitchedArray{Int32, 2}, Tuple{Base.UnitRange{Int64}, Base.UnitRange{Int64}},))
    precompile(CUDArt.cudacopy!, (Array{Any, 1}, CUDArt.CudaPitchedArray{UInt16, 2}, Tuple{Base.UnitRange{Int64}, Base.UnitRange{Int64}}, CUDArt.CudaArray{UInt16, 2}, Tuple{Base.UnitRange{Int64}, Base.UnitRange{Int64}},))
    precompile(CUDArt.cudacopy!, (Array{Any, 1}, CUDArt.CudaPitchedArray{Float64, 2}, Tuple{Base.UnitRange{Int64}, Base.UnitRange{Int64}}, CUDArt.CudaPitchedArray{Float64, 2}, Tuple{Base.UnitRange{Int64}, Base.UnitRange{Int64}},))
    precompile(CUDArt.copy!, (Array{Any, 1}, CUDArt.CudaArray{UInt16, 2}, Array{UInt16, 2},))
    precompile(CUDArt.copy!, (Array{Any, 1}, CUDArt.CudaArray{Int32, 2}, Array{Int32, 2},))
    precompile(CUDArt.copy!, (Array{Any, 1}, CUDArt.CudaArray{Int32, 2}, CUDArt.HostArray{Int32, 2},))
    precompile(CUDArt.devices, (Function, Array{Int64, 1},))
    precompile(CUDArt.close, (Int64,))
    precompile(CUDArt.close, (Array{Int64, 1},))
    precompile(CUDArt.call, (Type{CUDArt.CudaArray}, Type{Int64}, Tuple{Int64},))
    precompile(CUDArt.call, (Type{CUDArt.CudaArray}, Type{UInt16}, Tuple{Int64, Int64},))
    precompile(CUDArt.call, (Type{CUDArt.CudaArray}, Type{Int32}, Tuple{Int64, Int64},))
    precompile(CUDArt.size, (CUDArt.CudaPitchedArray{Float64, 2}, Int64,))
    precompile(CUDArt.call, (Type{CUDArt.CudaArray}, Type{Float32}, Tuple{Int64},))
    precompile(CUDArt.call, (Type{CUDArt.CudaArray}, Type{Float64}, Tuple{Int64, Int64},))
    precompile(CUDArt.size, (CUDArt.CudaPitchedArray{Int32, 2}, Int64,))
    precompile(CUDArt.capability, (Int64,))
    precompile(CUDArt.checkstrides_pitched, (CUDArt.CudaArray{UInt16, 2},))
    precompile(CUDArt.call, (Type{CUDArt.CudaArray}, Array{Float64, 2},))
    precompile(CUDArt.size, (CUDArt.CudaPitchedArray{UInt16, 2}, Int64,))
    precompile(CUDArt.size, (CUDArt.CudaPitchedArray{Int64, 1}, Int64,))
    precompile(CUDArt.stride, (CUDArt.CudaPitchedArray{Float64, 2}, Int64,))
    precompile(CUDArt.get_dim_x, (Tuple{Int64, Int64},))
    precompile(CUDArt.checkstrides_pitched, (CUDArt.CudaPitchedArray{Int32, 2},))
    precompile(CUDArt.checkstrides_pitched, (CUDArt.CudaPitchedArray{Float64, 2},))
    precompile(CUDArt.stride, (CUDArt.CudaPitchedArray{Int64, 1}, Int64,))
    precompile(CUDArt.checkstrides_pitched, (CUDArt.CudaPitchedArray{Int64, 1},))
    precompile(CUDArt.checkstrides_pitched, (CUDArt.CudaPitchedArray{UInt16, 2},))
    precompile(CUDArt.checkstrides_pitched, (Array{Float64, 2},))
    precompile(CUDArt.stride, (CUDArt.CudaPitchedArray{UInt16, 2}, Int64,))
    precompile(CUDArt.stride, (CUDArt.CudaPitchedArray{Int32, 2}, Int64,))
    precompile(CUDArt.pitchedptr, (CUDArt.CudaArray{UInt16, 2},))
    precompile(CUDArt.size, (CUDArt.CudaArray{UInt16, 2}, Int64,))
    precompile(CUDArt.checkstrides_pitched, (Array{Int32, 2},))
    precompile(CUDArt.checkstrides_pitched, (Array{UInt16, 2},))
    precompile(CUDArt.pitchedptr, (Array{Float64, 2},))
    precompile(CUDArt.pitchedptr, (Array{Int32, 2},))
    precompile(CUDArt.checkstrides_pitched, (Array{Int64, 1},))
    precompile(CUDArt.pitchedptr, (Array{UInt16, 2},))
    precompile(CUDArt.call, (Type{CUDArt.CuModule}, Function, ASCIIString,))
    precompile(CUDArt.pitchedptr, (Array{Int64, 1},))
    precompile(CUDArt.free, (CUDArt.CudaPitchedArray{Int64, 1},))
    precompile(CUDArt.get_dim_y, (Tuple{Int64, Int64},))
    precompile(CUDArt.copy!, (CUDArt.CudaPitchedArray{UInt16, 2}, Tuple{Base.UnitRange{Int64}, Base.UnitRange{Int64}}, Array{UInt16, 2}, Tuple{Base.UnitRange{Int64}, Base.UnitRange{Int64}},))
    precompile(CUDArt.call, (Type{CUDArt.CUDArt_gen.cudaPos}, Type{Int32}, Tuple{Int64, Int64},))
    precompile(CUDArt.copy!, (CUDArt.CudaPitchedArray{Float64, 2}, Tuple{Base.UnitRange{Int64}, Base.UnitRange{Int64}}, Array{Float64, 2}, Tuple{Base.UnitRange{Int64}, Base.UnitRange{Int64}},))
    precompile(CUDArt.copy!, (Array{Int32, 2}, Tuple{Base.UnitRange{Int64}, Base.UnitRange{Int64}}, CUDArt.CudaPitchedArray{Int32, 2}, Tuple{Base.UnitRange{Int64}, Base.UnitRange{Int64}},))
    precompile(CUDArt.call, (Type{CUDArt.CUDArt_gen.cudaExtent}, Type{Int32}, Tuple{Int64, Int64},))
    precompile(CUDArt.copy!, (Array{Float64, 2}, Tuple{Base.UnitRange{Int64}, Base.UnitRange{Int64}}, CUDArt.CudaPitchedArray{Float64, 2}, Tuple{Base.UnitRange{Int64}, Base.UnitRange{Int64}},))
    precompile(CUDArt.call, (Type{CUDArt.CUDArt_gen.cudaPos}, Type{Float64}, Tuple{Int64, Int64},))
    precompile(CUDArt.copy!, (CUDArt.CudaPitchedArray{Int32, 2}, Tuple{Base.UnitRange{Int64}, Base.UnitRange{Int64}}, Array{Int32, 2}, Tuple{Base.UnitRange{Int64}, Base.UnitRange{Int64}},))
    precompile(CUDArt.similar, (CUDArt.CudaPitchedArray{Int64, 1}, Tuple{Int64},))
    precompile(CUDArt.CUDArt_gen.call, (Type{CUDArt.CUDArt_gen.cudaMemcpy3DParms}, Ptr{Void}, CUDArt.CUDArt_gen.cudaPos, CUDArt.CUDArt_gen.cudaPitchedPtr, Ptr{Void}, CUDArt.CUDArt_gen.cudaPos, CUDArt.CUDArt_gen.cudaPitchedPtr, CUDArt.CUDArt_gen.cudaExtent, UInt32,))
    precompile(CUDArt.similar, (CUDArt.CudaArray{Int64, 1}, Tuple{Int64},))
    precompile(CUDArt.copy!, (Array{Int64, 1}, Tuple{Base.UnitRange{Int64}}, CUDArt.CudaPitchedArray{Int64, 1}, Tuple{Base.UnitRange{Int64}},))
    precompile(CUDArt.copy!, (Array{UInt16, 2}, Tuple{Base.UnitRange{Int64}, Base.UnitRange{Int64}}, CUDArt.CudaPitchedArray{UInt16, 2}, Tuple{Base.UnitRange{Int64}, Base.UnitRange{Int64}},))
    precompile(CUDArt.copy!, (CUDArt.CudaPitchedArray{Int64, 1}, Tuple{Base.UnitRange{Int64}}, Array{Int64, 1}, Tuple{Base.UnitRange{Int64}},))
    precompile(CUDArt.call, (Type{CUDArt.CUDArt_gen.cudaExtent}, Type{UInt16}, Tuple{Int64, Int64},))
    precompile(CUDArt.call, (Type{CUDArt.CUDArt_gen.cudaPos}, Type{UInt16}, Tuple{Int64, Int64},))
    precompile(CUDArt.copy!, (CUDArt.CudaPitchedArray{Float64, 2}, Tuple{Base.UnitRange{Int64}, Base.UnitRange{Int64}}, CUDArt.CudaPitchedArray{Float64, 2}, Tuple{Base.UnitRange{Int64}, Base.UnitRange{Int64}},))
    precompile(CUDArt.call, (Type{CUDArt.CUDArt_gen.cudaExtent}, Type{Float64}, Tuple{Int64, Int64},))
    precompile(CUDArt.copy!, (CUDArt.CudaArray{UInt16, 2}, Tuple{Base.UnitRange{Int64}, Base.UnitRange{Int64}}, CUDArt.CudaPitchedArray{UInt16, 2}, Tuple{Base.UnitRange{Int64}, Base.UnitRange{Int64}},))
    precompile(CUDArt.copy!, (CUDArt.CudaPitchedArray{UInt16, 2}, Tuple{Base.UnitRange{Int64}, Base.UnitRange{Int64}}, CUDArt.CudaArray{UInt16, 2}, Tuple{Base.UnitRange{Int64}, Base.UnitRange{Int64}},))
    precompile(CUDArt.free, (CUDArt.CudaPitchedArray{Float32, 1},))
    precompile(CUDArt.call, (Type{CUDArt.CudaArray}, Base.UnitRange{Int64},))
    precompile(CUDArt.call, (Type{CUDArt.CuFunction}, CUDArt.CuModule, ASCIIString,))
    precompile(CUDArt.similar, (CUDArt.CudaPitchedArray{Int64, 1},))
    precompile(CUDArt.similar, (CUDArt.CudaArray{Int64, 1},))
    precompile(CUDArt.free, (CUDArt.CudaPitchedArray{Float64, 2},))
    precompile(CUDArt.call, (Type{CUDArt.CUDArt_gen.cudaExtent}, Type{Int64}, Tuple{Int64},))
    precompile(CUDArt.call, (Type{CUDArt.CUDArt_gen.cudaPos}, Type{Int64}, Tuple{Int64},))
    precompile(CUDArt.free, (CUDArt.CudaPtr{Float32},))
    precompile(CUDArt.getindex, (CUDArt.HostArray{Int32, 2}, Int64, Int64,))
    precompile(CUDArt.free, (CUDArt.CudaPtr{Int32},))
    precompile(CUDArt.free, (CUDArt.CudaPtr{Int64},))
    precompile(CUDArt.free, (CUDArt.CudaPtr{Float64},))
    precompile(CUDArt.to_host, (CUDArt.CudaArray{Float64, 2},))
    precompile(CUDArt.free, (CUDArt.CudaPtr{UInt16},))
    precompile(CUDArt.free, (CUDArt.CudaPtr{UInt8},))
    precompile(CUDArt.call, (Type{CUDArt.CudaPitchedArray}, Base.UnitRange{Int64},))
    precompile(CUDArt.to_host, (CUDArt.CudaPitchedArray{Int32, 2},))
    precompile(CUDArt.call, (Type{CUDArt.CudaPtr}, Ptr{Int64},))
    precompile(CUDArt.to_host, (CUDArt.CudaPitchedArray{UInt16, 2},))
    precompile(CUDArt.pointer, (CUDArt.CudaPitchedArray{Int64, 1},))
    precompile(CUDArt.to_host, (CUDArt.CudaPitchedArray{Float64, 2},))
    precompile(CUDArt.to_host, (CUDArt.CudaPitchedArray{Int64, 1},))
    precompile(CUDArt.call, (Type{CUDArt.CudaArray{Float32, 1}}, CUDArt.CudaPtr{Float32}, Tuple{Int64}, Int64,))
    precompile(CUDArt.copy, (CUDArt.CudaPtr{Float64},))
    precompile(CUDArt.to_host, (CUDArt.CudaArray{Int32, 2},))
    precompile(CUDArt.to_host, (CUDArt.CudaArray{UInt16, 2},))
    precompile(CUDArt.to_host, (CUDArt.CudaArray{Int64, 1},))
    precompile(CUDArt.call, (Type{CUDArt.CudaPtr}, Ptr{UInt8},))
    precompile(CUDArt.device, (Int64,))
    precompile(CUDArt.call, (Type{CUDArt.CudaArray{UInt16, 2}}, CUDArt.CudaPtr{UInt16}, Tuple{Int64, Int64}, Int64,))
    precompile(CUDArt.copy!, (CUDArt.CudaPitchedArray{Float64, 2}, Array{Float64, 2},))
    precompile(CUDArt.call, (Type{CUDArt.CudaPtr}, Ptr{Int32},))
    precompile(CUDArt.call, (Type{CUDArt.CudaPtr}, Ptr{UInt16},))
    precompile(CUDArt.free, (CUDArt.CudaArray{Float64, 2},))
    precompile(CUDArt.call, (Type{CUDArt.CudaArray{Int32, 2}}, CUDArt.CudaPtr{Int32}, Tuple{Int64, Int64}, Int64,))
    precompile(CUDArt.call, (Type{CUDArt.CudaPtr}, Ptr{Float32},))
    precompile(CUDArt.call, (Type{CUDArt.CudaPtr}, Ptr{Float64},))
    precompile(CUDArt.call, (Type{CUDArt.CudaPitchedArray{Int64, 1}}, CUDArt.CUDArt_gen.cudaPitchedPtr, Tuple{Int64}, Int64,))
    precompile(CUDArt.call, (Type{CUDArt.CudaPitchedArray{Float64, 2}}, CUDArt.CUDArt_gen.cudaPitchedPtr, Tuple{Int64, Int64}, Int64,))
    precompile(CUDArt.call, (Type{CUDArt.CudaPitchedArray{Float32, 1}}, CUDArt.CUDArt_gen.cudaPitchedPtr, Tuple{Int64}, Int64,))
    precompile(CUDArt.call, (Type{CUDArt.CudaPitchedArray{UInt16, 2}}, CUDArt.CUDArt_gen.cudaPitchedPtr, Tuple{Int64, Int64}, Int64,))
    precompile(CUDArt.call, (Type{CUDArt.CudaArray{Float64, 2}}, CUDArt.CudaPtr{Float64}, Tuple{Int64, Int64}, Int64,))
    precompile(CUDArt.call, (Type{CUDArt.CudaPitchedArray{Int32, 2}}, CUDArt.CUDArt_gen.cudaPitchedPtr, Tuple{Int64, Int64}, Int64,))
    precompile(CUDArt.call, (Type{CUDArt.CudaArray{Int64, 1}}, CUDArt.CudaPtr{Int64}, Tuple{Int64}, Int64,))
    precompile(CUDArt.copy!, (CUDArt.CudaPitchedArray{Int32, 2}, Array{Int32, 2},))
    precompile(CUDArt.copy!, (Array{Int32, 2}, CUDArt.CudaArray{Int32, 2},))
    precompile(CUDArt.copy!, (CUDArt.CudaPitchedArray{UInt16, 2}, Array{UInt16, 2},))
    precompile(CUDArt.copy!, (CUDArt.CudaArray{UInt16, 2}, CUDArt.CudaPitchedArray{UInt16, 2},))
    precompile(CUDArt.copy!, (CUDArt.CudaPitchedArray{UInt16, 2}, CUDArt.CudaArray{UInt16, 2},))
    precompile(CUDArt.copy!, (CUDArt.CudaPitchedArray{Int64, 1}, Array{Int64, 1},))
    precompile(CUDArt.size, (CUDArt.CudaArray{Int64, 1},))
    precompile(CUDArt.size, (CUDArt.CudaArray{Float64, 2},))
    precompile(CUDArt.size, (CUDArt.CudaPitchedArray{Float64, 2},))
    precompile(CUDArt.copy!, (CUDArt.HostArray{Int32, 2}, CUDArt.CudaArray{Int32, 2},))
    precompile(CUDArt.pointer, (CUDArt.CudaPitchedArray{Float64, 2},))
    precompile(CUDArt.get!, (CUDArt.CudaPitchedArray{Float64, 2}, CUDArt.CudaPitchedArray{Float64, 2}, Tuple{Int64, Base.UnitRange{Int64}}, Float64,))
    precompile(CUDArt.size, (CUDArt.CudaPitchedArray{Int64, 1},))
    precompile(CUDArt.pointer, (CUDArt.CudaArray{Float64, 2},))
    precompile(CUDArt.similar, (CUDArt.CudaPitchedArray{Int64, 1}, Type{Float32},))
    precompile(CUDArt.similar, (CUDArt.CudaArray{Int64, 1}, Type{Float32},))
    precompile(CUDArt.fill!, (CUDArt.CudaArray{Int64, 1}, Int64,))
    precompile(CUDArt.fill!, (CUDArt.CudaArray{Float64, 2}, Int64,))
    precompile(CUDArt.fill!, (CUDArt.CudaPitchedArray{Float64, 2}, Float64,))
    precompile(CUDArt.fill!, (CUDArt.CudaPitchedArray{Int64, 1}, UInt8,))
    precompile(CUDArt.fill!, (CUDArt.CudaArray{Int64, 1}, UInt8,))
    precompile(CUDArt.fill!, (CUDArt.CudaArray{Float64, 2}, Float64,))
    precompile(CUDArt.fill!, (CUDArt.CudaArray{Int32, 2}, Int64,))
    precompile(CUDArt.fill!, (CUDArt.CudaPitchedArray{Int32, 2}, Int64,))
    precompile(CUDArt.eltype, (CUDArt.CudaArray{Int64, 1},))
    precompile(CUDArt.fill!, (CUDArt.CudaPitchedArray{Float64, 2}, Int64,))
    precompile(CUDArt.fill!, (CUDArt.CudaPitchedArray{Int64, 1}, Int64,))
    precompile(CUDArt.eltype, (CUDArt.CudaArray{Float32, 1},))
    precompile(CUDArt.eltype, (CUDArt.CudaPitchedArray{Float32, 1},))
    precompile(CUDArt.eltype, (CUDArt.CudaPitchedArray{Int64, 1},))
    precompile(CUDArt.devices, (Function,))
end


text:
# scal!
scal!(a::Array{T}, p_size) where T = LinearAlgebra.BLAS.scal!(length(a), T(12/p_size), a, 1)
scal!(a::CuArray{T}, p_size) where T = CUBLAS.scal!(length(a), T(1/p_size), a)

# Matvec
dispgemm!(tA::Char, α::Float32, A::AbstractArray{Float32, 2}, b::AbstractArray{Float32, 1}, β::Float32, c::Array{Float32, 1}) = LinearAlgebra.BLAS.gemv!(tA, α, A, b, β, c)
dispgemm!(tA::Char, α::Float32, A::AbstractArray{Float32, 2}, b::AbstractArray{Float32, 1}, β::Float32, c::CuArray{Float32, 1}) = CUBLAS.gemv!(tA, α, A, b, β, c)

# MatMat
dispgemm!(tA::Char, tB::Char, α::Float32, A::AbstractArray{Float32, 2}, b::AbstractArray{Float32, 2}, β::Float32, c::Array{Float32, 2}) = LinearAlgebra.BLAS.gemm!(tA, tB, α, A, b, β, c)
dispgemm!(tA::Char, tB::Char, α::Float32, A::AbstractArray{Float32, 2}, b::AbstractArray{Float32, 2}, β::Float32, c::CuArray{Float32, 2}) = CUBLAS.gemm!(tA, tB, α, A, b, β, c)

# Batched mat-mat
Bgemm!(tA, tB, α, A::AbstractArray, B::AbstractArray, β, C::Array) = NNlib.batched_gemm!(tA, tB, α, A, B, β, C)
Bgemm!(tA, tB, α, A::AbstractArray, B::AbstractArray, β, C::CuArray) = CUBLAS.gemm_strided_batched!(tA, tB, α, A, B, β, C)

max_stars_repo_name: sheevy/XConv

max_issues_repo_name: sheevy/XConv

max_forks_repo_name: sheevy/XConv


text:
# Used by backpropagation contexts during pullback computation.
merge(tp1::Tuple{}, tp2::Tuple{}) = tp1
merge(tp1::Tuple{Nothing}, tp2::Tuple{Nothing}) where T = tp1
merge(tp1::NTuple{N, Float64}, tp2::NTuple{N, Float64}) where N = [tp1[i] + tp2[i] for i in 1 : N]
merge(tp1::Array{Float64}, tp2::NTuple{N, Float64}) where N = [tp1[i] + tp2[i] for i in 1 : N]

# Used during specialization for regenerate/update for caching sites which don't need to be re-visited.
@inline function record_cached!(ctx, addr)
    visit!(ctx, addr)
    sub = get_sub(ctx.prev, addr)
    sc = get_score(sub)
    ctx.score += get_score(sub)
    set_sub!(ctx.tr, addr, sub)
    get_value(sub)
end


text:
using Stan, Turing

@model gdemo(x) = begin
  s ~ InverseGamma(2,3)
  m ~ Normal(0, sqrt(s))
  x[1] ~ Normal(m, sqrt(s))
  x[2] ~ Normal(m, sqrt(s))
  return s, m
end

model_f = gdemo([1.5, 2.0])

chn = sample(model_f, 2000, 1000, false, 1, Stan.Adapt(), Stan.Hmc())


text:
module TestUtils

using FoldsCUDA: ithtype
using Test

function test_ithtype()
    for n in 2:5
        T = Union{(Val{i} for i in 1:n)...}
        @test Union{(ithtype(T, Val(i)) for i in 1:n)...} == T
    end
end

end  # module


text:
using StatsBase, JLD,  Hwloc

initWorkers()
using DynNets , StatsFuns
@everywhere using DynNets, StatsFuns
N_est = 1000
Tvals =  [150 500 5000]
N_ind = length(Tvals)[1]
scoreRescType = "FISHER-EWMA"#""#
useStartVal = true#false#  false
N = 10 ;
GBAdgp =1; GWdgp = 1 ;
GBAest =GBAdgp; GWest = GWdgp;

#Distribuisci i nodi tra i gruppi in parti uguali

groupsIndsDgp = [[DynNets.distributeAinVecN(Array{Int}(1:GWdgp),N)];
                    [DynNets.distributeAinVecN(Array{Int}(1:GBAdgp),GWdgp)]]
groupsIndsEst = [[DynNets.distributeAinVecN(Array{Int}(1:GWest),N)];
                   [DynNets.distributeAinVecN(Array{Int}(1:GBAest),GWest)]]

#the values for the unconidtional means to be used in the DGP are obtained estimating
#an equally spaced degree sequence
est_times = Array{Float64}(zeros(N_est))
est_conv_flag = falses(N_est)
est_flag = falses(N_est)
# Set the value of the parameters common to all replicate of DGP+Estimation

realUm , dgpdegs  =   DynNets.linSpacedPar( DynNets.fooSdErgm1,N; NgroupsW = GWdgp)

sim_B = Array(linspace(0.4,0.9,GBAdgp+1)[2:end]) #
sim_A = Array(linspace(0.02,0.2,GBAdgp+1)[2:end]) #
sim_W = realUm.*(1-sim_B[groupsIndsDgp[2]])
#Initialize the variable to store
sizesParSim = [GWdgp,GBAdgp,GBAdgp]
sizesParEst = [GWest,GBAdgp,GBAest]
for i =1:3
    ex1 = "estPar_$(i)  = SharedArray{Float64}(zeros(N_ind,N_est, $(sizesParEst[i]) ))"
    ex2 = "simPar_$(i) = SharedArray{Float64}(zeros(N_ind,N_est, $(sizesParSim[i])))"
    eval(parse(ex1))
    eval(parse(ex2))
 end
allObs = SharedArray{Int16,4}(N,Tvals[end],N_ind,N_est)
est_times = SharedArray{Float64}(zeros(N_ind,N_est))
est_conv_flag = SharedArray{Bool,2}(N_ind,N_est)
est_flag = SharedArray{Bool,2}(N_ind,N_est)
@progress "T"  for ind = 1:N_ind
#   ind=1
    T = Tvals[ind]
     @sync @parallel   for n = 1:N_est
 #n=1
        simPar_3[ind,n,1:GBAdgp] = sim_A
        simPar_2[ind,n,1:GBAdgp] = sim_B
        simPar_1[ind,n,1:GWdgp] = sim_W
        Model2Sample = DynNets.SdErgm1(ones(T,N),[sim_W,sim_B,sim_A],
                                                    groupsIndsDgp,scoreRescType)
        Model2Sample,Y_T,Fitness_T  =  DynNets.sampl(Model2Sample,T)
        sampleDegsT = squeeze(sum(Y_T,2),2)
        allObs[:,1:T,ind,n] = sampleDegsT'
        Model2Est =  DynNets.SdErgm1(sampleDegsT,[zeros(GWest),zeros(GBAest),zeros(GBAest)],
                                                    groupsIndsEst,scoreRescType)
        #estimate
        useStartVal ?   start_val = [sim_W, sim_B, sim_A] : start_val = zeros(3,3)
        est_times[ind,n] = @elapsed  hat_pars , conv_flag =
                            DynNets.estimate(Model2Est,start_values = start_val )
        estPar_1[ind,n,:] = hat_pars[1]
        estPar_2[ind,n,:] = hat_pars[2]
        estPar_3[ind,n,:] = hat_pars[3]
        est_conv_flag[ind,n] = conv_flag
        est_flag[ind,n] = true

    end

    # save
    save_fold = "./data/estimatesTest/"

    file_name = "EstTest_$(N)_$(GWdgp)_$(GBAdgp)_$(Tvals)_$(N_est)_$(useStartVal)_Miss_$(GWest)_$(GBAest)"* scoreRescType *".jld"
    save_path = save_fold*file_name#
    @save(save_path,groupsIndsEst,groupsIndsDgp,allObs,dgpdegs,simPar_1,simPar_2,simPar_3,estPar_1,estPar_2,estPar_3,est_conv_flag,est_times)
    println([N,GWdgp,GWest,T,sum(est_times[ind,:]),sum(est_conv_flag[ind,:])/sum(est_flag[ind,:])])
end

##

 ##


text:
using Publish
using Artifacts, LazyArtifacts
using BitSAD

# override default theme
Publish.Themes.default() = artifact"flux-theme"

p = Publish.Project(BitSAD)

# needed to prevent error when overwriting
rm("dev", recursive = true, force = true)
rm(p.env["version"], recursive = true, force = true)

# build documentation
deploy(BitSAD; root = "/BitSAD.jl", force = true, label = "dev")


text:
# Meshes a quadraterial area given by four corners with nx * xy * 2 triangles.

# TODO: Add quadraterial option
# TODO: Add test
# TODO: Check if the corners does not represent a convex set.

   # Node number for column i and row j

function meshquad(nx::Int, ny::Int, corners::Matrix{Float64},
                 ele_type::Union(Type{GeoQuad}, Type{GeoTrig}, Type{GeoQTrig}))

    if size(corners) != (2,4)
        error("corner argument needs to be a 2x4 matrix")
    end

    if (nx <= 0 || ny <= 0)
        error("Need at least 1x1 elements")
    end

    mesh = GeoMesh()

    #... you what Lint?
    @lintpragma("Ignore use of undeclared variable :")
    # The four corners
    LL = corners[:, 1]
    UL = corners[:, 2]
    UR = corners[:, 3]
    LR = corners[:, 4]

    if ele_type == GeoQTrig
        np_y = 2*ny
        np_x = 2*nx
    else
        np_y = ny
        np_x = ny
    end

     @inline node_nr(i::Int, j::Int) = (np_x + 1) * (i-1) + j

    # Add the nodes
    for i in 0:np_y
        ratio_bounds = i / np_y

        x0 = LL[1] * (1 - ratio_bounds) + ratio_bounds * UL[1]
        x1 = LR[1] * (1 - ratio_bounds) + ratio_bounds * UR[1]

        y0 = LL[2] * (1 - ratio_bounds) + ratio_bounds * UL[2]
        y1 = LR[2] * (1 - ratio_bounds) + ratio_bounds * UR[2]

        for j in 0:np_x
            ratio = j / np_x
            x = x0 * (1 - ratio) + ratio * x1
            y = y0 * (1 - ratio) + ratio * y1
            push!(mesh, GeoNode2(node_nr(i+1, j+1), [x,y]))
        end
    end

    # Add the elements
    n_elem = 0
       if ele_type == GeoQTrig
        jump = 2
    else
        jump = 1
    end

    for i in 1:jump:jump*ny
        for j in 1:jump:jump*nx
            if ele_type == GeoTrig
                n_elem += 1
                elem1 = GeoTrig(n_elem, Vertex3(node_nr(i, j), node_nr(i+1, j), node_nr(i+1, j+1)))
                push!(mesh, elem1)
                n_elem += 1
                elem2 = GeoTrig(n_elem, Vertex3(node_nr(i, j), node_nr(i+1, j+1), node_nr(i, j+1)))
                push!(mesh, elem2)
            elseif ele_type == GeoQTrig
                n_elem += 1
                elem1 = GeoQTrig(n_elem, Vertex6(node_nr(i, j), node_nr(i+2, j+2), node_nr(i+2, j),
                                                 node_nr(i+1, j+1), node_nr(i+2, j+1), node_nr(i+1, j)))
                push!(mesh, elem1)
                n_elem += 1
                elem1 = GeoQTrig(n_elem, Vertex6(node_nr(i, j), node_nr(i, j+2), node_nr(i+2, j+2),
                                                 node_nr(i, j+1), node_nr(i+1, j+2), node_nr(i+1, j+1)))
                push!(mesh, elem1)
            else
                n_elem += 1
                elem = GeoQuad(n_elem, Vertex4(node_nr(i,j), node_nr(i,j+1), node_nr(i+1,j+1), node_nr(i+1,j)))
                push!(mesh, elem)
            end
        end
    end
    return mesh
end

# Generates what is known as the "Cook membrane"
function gencook(nx, ny, ele_type::Union(Type{GeoQuad}, Type{GeoTrig}, Type{GeoQTrig}), scale = 1.0)
    m = scale * [[0.0; 0.0] [0.0; 44.0] [48.0; 60.0] [48.0; 44.0]]
    meshquad(nx, ny, m, ele_type)
end


text:
if Sys.islinux()
    run(`wget https://www.dropbox.com/s/1dr7g8x270xdrup/whisk-1.1.0d-64bit-Linux.tar.gz`)
    run(`tar -xzf whisk-1.1.0d-64bit-Linux.tar.gz`)
    rm("whisk-1.1.0d-64bit-Linux.tar.gz")
    run(`mv whisk-1.1.0d-Linux whisk`)

    run(`wget https://www.dropbox.com/s/s3ll0taalzqz0c5/quad_hourglass_64.mat`)
elseif Sys.iswindows()
    run(`curl.exe --url "https://www.dropbox.com/s/0dwksb5pihauvz8/whisker.zip?dl=0" -o whisk.zip -L`)
    run(`powershell.exe -nologo -noprofile -command "Expand-Archive -Literalpath ./whisk.zip"`)
    rm("whisk.zip")

    run(`curl.exe --url "https://www.dropbox.com/s/s3ll0taalzqz0c5/quad_hourglass_64.mat?dl=0" -o quad_hourglass_64.mat`)
elseif Sys.isosx()
    println("No Mac support!")
else
    println("Operating system not recognized.")
end


text:
# Computed via wolframalpha.com: SinIntegral[SetPrecision[Table[x,{x, 1,20,1}],20]] and CosIntegral[SetPrecision[Table[x,{x, 1,20,1}],20]]

@testset "trigonometric integrals" begin
    sinintvals = [  0.9460830703671830149, 1.6054129768026934849, 1.848652527999468256,
                    1.75820313894905306, 1.54993124494467414, 1.4246875512805065,
                    1.4545966142480936, 1.5741868217069421, 1.665040075829602,
                    1.658347594218874, 1.578306806945727416, 1.504971241526373371,
                    1.499361722862824564, 1.556211050077665054, 1.618194443708368739,
                    1.631302268270032886, 1.590136415870701122, 1.536608096861185462,
                    1.518630031769363932, 1.548241701043439840]
    cosintvals = [0.3374039229009681347, 0.4229808287748649957, 0.119629786008000328,
                    -0.14098169788693041, -0.19002974965664388, -0.06805724389324713,
                    0.07669527848218452, 0.122433882532010, 0.0553475313331336,
                    -0.045456433004455, -0.08956313549547997948, -0.04978000688411367560,
                    0.02676412556403455504, 0.06939635592758454727, 0.04627867767436043960,
                    -0.01420019012019002240, -0.05524268226081385053, -0.04347510299950100478,
                    0.00515037100842612857, 0.04441982084535331654]
    for x in 1:20
        @test sinint(x)     ≅ sinintvals[x]
        @test sinint(-x)    ≅ -sinintvals[x]
        @test cosint(x)     ≅ cosintvals[x]
    end

    @test sinint(Float16(1)) == Float16(sinint(1))
    @test cosint(Float16(1)) == Float16(cosint(1))
    @test sinint(Float32(1)) == Float32(sinint(1))
    @test cosint(Float32(1)) == Float32(cosint(1))

    @test sinint(1//2) == sinint(0.5)
    @test cosint(1//2) == cosint(0.5)

    @test sinint(1e300)     ≅ π/2
    @test cosint(1e300)     ≅ -8.17881912115908554103E-301
    @test sinint(1e-300)    ≅ 1.0E-300
    @test cosint(1e-300)    ≅ -690.1983122333121

    @test sinint(Inf) == π/2
    @test cosint(Inf) == 0.0
    @test isnan(sinint(NaN))
    @test isnan(cosint(NaN))

    @test_throws ErrorException sinint(big(1.0))
    @test_throws ErrorException cosint(big(1.0))
    @test_throws DomainError cosint(-1.0)
    @test_throws DomainError cosint(Float32(-1.0))
    @test_throws DomainError cosint(Float16(-1.0))
    @test_throws DomainError cosint(-1//2)
    @test_throws MethodError sinint(Complex(1))
    @test_throws MethodError cosint(Complex(1))
end


text:
### A Pluto.jl notebook ###
# v0.16.1

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 3f890c7b-44ad-426a-bf26-32ff17fa8ec4
md"""
# Examples in Julia - Tutorial 1
"""

# ╔═╡ 5bc2c872-fded-4bc8-8044-9eda1822c505
md"""
## Assignment 1

Using the package `Polynomials.jl`, write the function which implements
[Graeffe's method](https://en.wikipedia.org/wiki/Graeffe%27s_method)
(see also [here](http://mathworld.wolfram.com/GraeffesMethod.html))
for computing roots of polynomials with only real roots with simple moduli.

In the function, use Julia's `BigFloat` numbers to overcome the main disadvantage of the method. What is the number of significant decimal digits, and the largest and the smallest number?

Test the function on the [Wilkinson's polynomial](https://en.wikipedia.org/wiki/Wilkinson%27s_polynomial) $\omega(x)$, and the [Chebyshev polynomial](https://en.wikipedia.org/wiki/Chebyshev_polynomials) $T_{50}(x)$
(the latter needs to be transformed in order to apply the method).

Compare your solutions with the exact solutions.
"""

# ╔═╡ 5b542ad3-8471-4adc-9c89-d647000b8ca0
md"""
## Assignment 2

Write the function which computes simple LU factorization (without pivoting) where the matrix is overwritten by the factors.

Make sure that the function also works with block-matrices.

Compare the speed on standard matrices and block-matrices with the built-in LU factorization (which also uses block algorithm AND pivoting). Check the accuracy.
"""

# ╔═╡ dbd0a150-2ad2-49e3-b585-d34bfc2ff3b8
md"""
## Assignment 3

Use the function `eigvals()` to compute the eigenvalues of $k$ random matrices (with uniform and normal distribution of elements) of order $n$.

Plot the results using the macro `@bind` from the package `PlutoUI.jl`.
Use `Plots.jl` for plotting.

Are the eigenvalues random? Can you describe their behaviour? Can random matrices be used to test numerical algorithms?
"""

# ╔═╡ 4cd45baf-67f6-43cc-a9bf-3a89006c3814

# ╔═╡ Cell order:
# ╟─3f890c7b-44ad-426a-bf26-32ff17fa8ec4
# ╟─5bc2c872-fded-4bc8-8044-9eda1822c505
# ╟─5b542ad3-8471-4adc-9c89-d647000b8ca0
# ╟─dbd0a150-2ad2-49e3-b585-d34bfc2ff3b8
# ╠═4cd45baf-67f6-43cc-a9bf-3a89006c3814

max_stars_repo_path: src/Module A - Short Julia Course/T1 Examples in Julia.jl

max_stars_repo_name: ivanslapnicar/GIAN-Applied-NLA-Course

max_issues_repo_name: ivanslapnicar/GIAN-Applied-NLA-Course

max_forks_repo_name: ivanslapnicar/GIAN-Applied-NLA-Course


text:
using DataStructures
using Test
using ChemistryFeaturization
include("../../src/weave_fcns.jl")
using .weave_fcns:
    smiles_atom_features,
    smiles_bond_features,
    chem,
    atom_feat_list,
    bond_feat_bins,
    atom_feat_bins

mol = chem.MolFromSmiles("C1=CC=CC=C1")  # benzene
carbon = get(mol.GetAtoms(), 0)
bond = get(mol.GetBonds(), 0)

@testset "SMILES atom featurization" begin
    @test [k for k in keys(smiles_atom_features(carbon))] == [
        "symbol",
        "degree",
        "implicit_valence",
        "formal_charge",
        "radical_electrons",
        "hybridization",
        "aromaticity",
        "total_H_num",
    ]
    try
        smiles_atom_features(carbon, feature_list = ["invalid_feat"])
    catch e
        @test occursin("is not in supported atomic features list", sprint(showerror, e))
    end
    correct_values = [
        Float32.("C" .== atom_feat_bins["symbol"]),
        Float32.(2 .== atom_feat_bins["degree"]),
        Float32.(1 .== atom_feat_bins["implicit_valence"]),
        Float32(carbon.GetFormalCharge()),
        Float32(carbon.GetNumRadicalElectrons()),
        Float32.(chem.rdchem.HybridizationType.SP2 .== atom_feat_bins["hybridization"]),
        Float32(carbon.GetIsAromatic()),
        Float32.(1 .== atom_feat_bins["total_H_num"]),
    ]
    @test [v for v in values(smiles_atom_features(carbon))] == correct_values
end

@testset "SMILES bond featurization" begin
    @test [k for k in keys(smiles_bond_features(bond))] == ["bond_type", "isConjugated", "isInring"]
    try
        smiles_bond_features(carbon, feature_list = ["invalid_feat"])
    catch e
        @test occursin("is not in supported bond features list", sprint(showerror, e))
    end
    correct_values = [
        Float32.(chem.rdchem.BondType.AROMATIC .== bond_feat_bins["bond_type"]),
        Float32(bond.GetIsConjugated()),
        Float32(bond.IsInRing()),
    ]
    @test [v for v in values(smiles_bond_features(bond))] == correct_values
end


text:
"""
    Tensor{T}

A tensor in a network with index labels of type T.

# Fields
- `adj::Vector{T}`: the labels of the adjacent tensors
- `arr::Array{Float64}`: the contents of the tensor
- `x::Float64`, `y::Float64`: position coordinates of the tensor

The order with which the neighbours are labelled in `adj` may be arbitrary, but must be
compatible with the labelling of the indices of `arr`.
"""
mutable struct Tensor{T}
    adj::Vector{T}
    arr::Array{Float64}
    x::Float64
    y::Float64

    # Native float constructors
    Tensor{T}(adj::Vector{T}, arr::Array{Float64}, x::Float64, y::Float64) where T =
        new{T}(adj, arr, x, y)
    Tensor(adj::Vector{T}, arr::Array{Float64}, x::Float64, y::Float64) where T =
        new{T}(adj, arr, x, y)

    # Promoted constructors
    Tensor{T}(adj::Vector{T}, arr::Array{<:Real}, x::Real, y::Real) where T =
        new{T}(adj, float(arr), float(x), float(y))
    Tensor(adj::Vector{T}, arr::Array{<:Real}, x::Real, y::Real) where T =
        new{T}(adj, float(arr), float(x), float(y))
end

"""
    TensorNetwork

A tensor network with edges labelled by consecutive integers. Alias for
`Vector{Tensor{Int}}`.
"""
const TensorNetwork = Vector{Tensor{Int}}
labels(TN::TensorNetwork) = keys(TN)
tensors(TN::TensorNetwork) = values(TN)

"""
    LabelledTensorNetwork{T}

A tensor network with tensor labels of type `T`. Alias for `Dict{T,Tensor{T}}`. Due to the
potential performacne overheads associated with dealing with arbitrarily labelled tensors
the alternative `TensorNetwork` type should be used in performance sensitive circumstances.
"""
const LabelledTensorNetwork = Dict{T,Tensor{T}} where T
labels(LTN::LabelledTensorNetwork) = keys(LTN)
tensors(LTN::LabelledTensorNetwork) = values(LTN)

# Convert LabelledTensorNetwork into TensorNetwork for easier manipulation
function delabel(LTN::LabelledTensorNetwork{T}) where T
    N = length(LTN)
    # Construct an (arbitrary) mapping from labels to indices
    delabelling = Dict{T,Int}()
    sizehint!(delabelling, N)
    for (i,l) ∈ enumerate(labels(LTN))
        delabelling[l] = i
    end
    # Construct a TensorNetwork with this labelling
    TN = TensorNetwork(undef,N)
    for (l,t) ∈ LTN
        TN[delabelling[l]] = Tensor{Int}(map(λ->delabelling[λ],t.adj), t.arr, t.x, t.y)
    end
    return TN
end

struct InvalidTNError <: Exception
    str::String
end

Base.showerror(io::IO, e::InvalidTNError) = print(io, "Invalid tensor network: ", e.str)

# Check the validity of a tensor network and the tensors therein. If invalid then
# InvalidTNError is thrown, and otherwise nothing is returned
function checkvalid(LTN::LabelledTensorNetwork)
    for l ∈ labels(LTN)
        # Check that the graph degree matches the order of the tensor
        ndims(LTN[l].arr)!=length(LTN[l].adj) && throw(InvalidTNError("graph degree != tensor order"))
        # Check for multi-edges
        !allunique(LTN[l].adj) && throw(InvalidTNError("multi-edge"))
        for (i,k) ∈ enumerate(LTN[l].adj)
            # Check for self-edges
            l==k && throw(InvalidTNError("self-edge"))
            l>k && continue
            # Find edge back
            j = findfirst(isequal(l),LTN[k].adj)
            isnothing(j) && throw(InvalidTNError("directed edge"))
            # Make sure tensor dimensions are consistent along edge
            size(LTN[l].arr,i)!=size(LTN[k].arr,j) && throw(InvalidTNError("dimension mismatch"))
        end
    end
    return nothing
end

function checkvalid(TN::TensorNetwork)
    for u ∈ eachindex(TN)
        ndims(TN[u].arr)!=length(TN[u].adj) && throw(InvalidTNError("graph degree != tensor order"))
        !allunique(TN[u].adj) && throw(InvalidTNError("multi-edge"))
        for (i,v) ∈ enumerate(TN[u].adj)
            u==v && throw(InvalidTNError("self-edge"))
            u>v && continue
            j = findfirst(isequal(u),TN[v].adj)
            isnothing(j) && throw(InvalidTNError("directed edge"))
            size(TN[u].arr,i)!=size(TN[v].arr,j) && throw(InvalidTNError("dimension mismatch"))
        end
    end
    return nothing
end


text:
using Test

@testset "Trixi2Vtk" begin
  include("test_2d.jl")
  include("test_3d.jl")
  include("test_manual.jl")
end


text:
import ForneyLab: ProbabilityDistribution, Univariate, GaussianMeanVariance
import Plots: plot, plot!, xlabel!, ylabel!

function plot_messages(state_prediction::ProbabilityDistribution{Univariate, GaussianMeanVariance},
                       measurement_likelihood::ProbabilityDistribution{Univariate, GaussianMeanVariance},
                       corrected_state_prediction::ProbabilityDistribution{Univariate, GaussianMeanVariance})

    # Extract parameters
    m_p = mean(state_prediction)
    v_p = var(state_prediction)
    m_l = mean(measurement_likelihood)
    v_l = var(measurement_likelihood)
    m_c = mean(corrected_state_prediction)
    v_c = var(corrected_state_prediction)

    # Set plot range
    range_p = [m_p - 8*sqrt(v_p), m_p + 8*sqrt(v_p)]
    range_l = [m_l - 8*sqrt(v_l), m_l + 8*sqrt(v_l)]
    range_c = [m_c - 8*sqrt(v_c), m_c + 8*sqrt(v_c)]
    plot_range = range(minimum([range_p[1], range_l[1], range_c[1]]),
                       maximum([range_p[2], range_l[2], range_c[2]]), length=300)

    # Probabilities over range
    probs_p = [pdf(state_prediction, a) for a in plot_range]
    probs_l = [pdf(measurement_likelihood, a) for a in plot_range]
    probs_c = [pdf(corrected_state_prediction, a) for a in plot_range]

    # Start plotting
    plot(plot_range, probs_p, color="red", label="state prediction", size=(800,500))
    plot!(plot_range, probs_l, color="blue", label="observation likelihood")
    plot!(plot_range, probs_c, color="green", label="state estimate")
    xlabel!("x-coordinate")
    ylabel!("Probability")
end

function pdf(dist::ProbabilityDistribution{Univariate, GaussianMeanVariance}, x::Number)
    return exp(-0.5*(log(2*pi) + log(dist.params[:v]) + (x - dist.params[:m])^2 / dist.params[:v]))
end;

max_stars_repo_name: bertdv/BMLIP