text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

using Base.Iterators.Enumerate

"""
    asyncmap(f, c...; ntasks=0, batch_size=nothing)

Uses multiple concurrent tasks to map `f` over a collection (or multiple
equal length collections). For multiple collection arguments, `f` is
applied elementwise.

`ntasks` specifies the number of tasks to run concurrently.
Depending on the length of the collections, if `ntasks` is unspecified,
up to 100 tasks will be used for concurrent mapping.

`ntasks` can also be specified as a zero-arg function. In this case, the
number of tasks to run in parallel is checked before processing every element and a new
task started if the value of `ntasks_func` is less than the current number
of tasks.

If `batch_size` is specified, the collection is processed in batch mode. `f` must
then be a function that must accept a `Vector` of argument tuples and must
return a vector of results. The input vector will have a length of `batch_size` or less.

The following examples highlight execution in different tasks by returning
the `object_id` of the tasks in which the mapping function is executed.

First, with `ntasks` undefined, each element is processed in a different task.
```
julia> tskoid() = object_id(current_task());

julia> asyncmap(x->tskoid(), 1:5)
5-element Array{UInt64,1}:
 0x6e15e66c75c75853
 0x440f8819a1baa682
 0x9fb3eeadd0c83985
 0xebd3e35fe90d4050
 0x29efc93edce2b961

julia> length(unique(asyncmap(x->tskoid(), 1:5)))
5
```

With `ntasks=2` all elements are processed in 2 tasks.
```
julia> asyncmap(x->tskoid(), 1:5; ntasks=2)
5-element Array{UInt64,1}:
 0x027ab1680df7ae94
 0xa23d2f80cd7cf157
 0x027ab1680df7ae94
 0xa23d2f80cd7cf157
 0x027ab1680df7ae94

julia> length(unique(asyncmap(x->tskoid(), 1:5; ntasks=2)))
2
```

With `batch_size` defined, the mapping function needs to be changed to accept an array
of argument tuples and return an array of results. `map` is used in the modified mapping
function to achieve this.
```
julia> batch_func(input) = map(x->string("args_tuple: ", x, ", element_val: ", x[1], ", task: ", tskoid()), input)
batch_func (generic function with 1 method)

julia> asyncmap(batch_func, 1:5; ntasks=2, batch_size=2)
5-element Array{String,1}:
 "args_tuple: (1,), element_val: 1, task: 9118321258196414413"
 "args_tuple: (2,), element_val: 2, task: 4904288162898683522"
 "args_tuple: (3,), element_val: 3, task: 9118321258196414413"
 "args_tuple: (4,), element_val: 4, task: 4904288162898683522"
 "args_tuple: (5,), element_val: 5, task: 9118321258196414413"
```

!!! note
    Currently, all tasks in Julia are executed in a single OS thread co-operatively. Consequently,
    `ayncmap` is beneficial only when the mapping function involves any I/O - disk, network, remote
    worker invocation, etc.

"""
function asyncmap(f, c...; ntasks=0, batch_size=nothing)
    return async_usemap(f, c...; ntasks=ntasks, batch_size=batch_size)
end

function async_usemap(f, c...; ntasks=0, batch_size=nothing)
    ntasks = verify_ntasks(c[1], ntasks)
    batch_size = verify_batch_size(batch_size)

    if batch_size !== nothing
        exec_func = batch -> begin
            # extract the Refs from the input tuple
            batch_refs = map(x->x[1], batch)

            # and the args tuple....
            batched_args = map(x->x[2], batch)

            results = f(batched_args)
            foreach(x -> (batch_refs[x[1]].x = x[2]), enumerate(results))
        end
    else
        exec_func = (r,args) -> (r.x = f(args...))
    end
    chnl, worker_tasks = setup_chnl_and_tasks(exec_func, ntasks, batch_size)
    return wrap_n_exec_twice(chnl, worker_tasks, ntasks, exec_func, c...)
end

batch_size_err_str(batch_size) = string("batch_size must be specified as a positive integer. batch_size=", batch_size)
function verify_batch_size(batch_size)
    if batch_size === nothing
        return batch_size
    elseif isa(batch_size, Number)
        batch_size = Int(batch_size)
        batch_size < 1 && throw(ArgumentError(batch_size_err_str(batch_size)))
        return batch_size
    else
        throw(ArgumentError(batch_size_err_str(batch_size)))
    end
end

function verify_ntasks(iterable, ntasks)
    if !((isa(ntasks, Number) && (ntasks >= 0)) || isa(ntasks, Function))
        err = string("ntasks must be specified as a positive integer or a 0-arg function. ntasks=", ntasks)
        throw(ArgumentError(err))
    end

    if ntasks == 0
        chklen = iteratorsize(iterable)
        if (chklen == HasLength()) || (chklen == HasShape())
            ntasks = max(1,min(100, length(iterable)))
        else
            ntasks = 100
        end
    end
    return ntasks
end

function wrap_n_exec_twice(chnl, worker_tasks, ntasks, exec_func, c...)
    # The driver task, creates a Ref object and writes it and the args tuple to
    # the communication channel for processing by a free worker task.
    push_arg_to_channel = (x...) -> (r=Ref{Any}(nothing); put!(chnl,(r,x));r)

    if isa(ntasks, Function)
        map_f = (x...) -> begin
            # check number of tasks every time, and start one if required.
            # number_tasks > optimal_number is fine, the other way around is inefficient.
            if length(worker_tasks) < ntasks()
                start_worker_task!(worker_tasks, exec_func, chnl)
            end
            push_arg_to_channel(x...)
        end
    else
        map_f = push_arg_to_channel
    end
    maptwice(map_f, chnl, worker_tasks, c...)
end

function maptwice(wrapped_f, chnl, worker_tasks, c...)
    # first run, returns a collection of Refs
    asyncrun_excp = nothing
    local asyncrun
    try
        asyncrun = map(wrapped_f, c...)
    catch ex
        if isa(ex,InvalidStateException)
            # channel could be closed due to exceptions in the async tasks,
            # we propagate those errors, if any, over the `put!` failing
            # in asyncrun due to a closed channel.
            asyncrun_excp = ex
        else
            rethrow(ex)
        end
    end

    # close channel and wait for all worker tasks to finish
    close(chnl)

    # check and throw any exceptions from the worker tasks
    foreach(x->(v=wait(x); isa(v, Exception) && throw(v)), worker_tasks)

    # check if there was a genuine problem with asyncrun
    (asyncrun_excp !== nothing) && throw(asyncrun_excp)

    if isa(asyncrun, Ref)
        # scalar case
        return asyncrun.x
    else
        # second run, extract values from the Refs and return
        return map(ref->ref.x, asyncrun)
    end
end

function setup_chnl_and_tasks(exec_func, ntasks, batch_size=nothing)
    if isa(ntasks, Function)
        nt = ntasks()
        # start at least one worker task.
        if nt == 0
            nt = 1
        end
    else
        nt = ntasks
    end

    # Use an unbuffered channel for communicating with the worker tasks. In the event
    # of an error in any of the worker tasks, the channel is closed. This
    # results in the `put!` in the driver task failing immediately.
    chnl = Channel(0)
    worker_tasks = []
    foreach(_ -> start_worker_task!(worker_tasks, exec_func, chnl, batch_size), 1:nt)
    yield()
    return (chnl, worker_tasks)
end

function start_worker_task!(worker_tasks, exec_func, chnl, batch_size=nothing)
    t = @schedule begin
        retval = nothing

        try
            if isa(batch_size, Number)
                while isopen(chnl)
                    # The mapping function expects an array of input args, as it processes
                    # elements in a batch.
                    batch_collection=Any[]
                    n = 0
                    for exec_data in chnl
                        push!(batch_collection, exec_data)
                        n += 1
                        (n == batch_size) && break
                    end
                    if n > 0
                        exec_func(batch_collection)
                    end
                end
            else
                for exec_data in chnl
                    exec_func(exec_data...)
                end
            end
        catch e
            close(chnl)
            retval = e
        end
        retval
    end
    push!(worker_tasks, t)
end

# Special handling for some types.
function asyncmap(f, s::AbstractString; kwargs...)
    s2 = Array{Char,1}(length(s))
    asyncmap!(f, s2, s; kwargs...)
    return convert(String, s2)
end

# map on a single BitArray returns a BitArray if the mapping function is boolean.
function asyncmap(f, b::BitArray; kwargs...)
    b2 = async_usemap(f, b; kwargs...)
    if eltype(b2) == Bool
        return BitArray(b2)
    end
    return b2
end

# TODO: Optimize for sparse arrays
# For now process as regular arrays and convert back
function asyncmap(f, s::AbstractSparseArray...; kwargs...)
    sa = map(Array, s)
    return sparse(asyncmap(f, sa...; kwargs...))
end

mutable struct AsyncCollector
    f
    results
    enumerator::Enumerate
    ntasks
    batch_size
    nt_check::Bool     # check number of tasks on every iteration

    AsyncCollector(f, r, en::Enumerate, ntasks, batch_size) = new(f, r, en, ntasks, batch_size, isa(ntasks, Function))
end

"""
    AsyncCollector(f, results, c...; ntasks=0, batch_size=nothing) -> iterator

Returns an iterator which applies `f` to each element of `c` asynchronously
and collects output into `results`.

Keyword args `ntasks` and `batch_size` have the same behavior as in
[`asyncmap`](@ref). If `batch_size` is specified, `f` must
be a function which operates on an array of argument tuples.

!!! note
    `next(::AsyncCollector, state) -> (nothing, state)`. A successful return
    from `next` indicates that the next element from the input collection is
    being processed asynchronously. It blocks until a free worker task becomes
    available.

!!! note
    `for _ in AsyncCollector(f, results, c...; ntasks=1) end` is equivalent to
    `map!(f, results, c...)`.
"""
function AsyncCollector(f, results, c...; ntasks=0, batch_size=nothing)
    AsyncCollector(f, results, enumerate(zip(c...)), ntasks, batch_size)
end

mutable struct AsyncCollectorState
    chnl::Channel
    worker_tasks::Array{Task,1}
    enum_state      # enumerator state
end

function start(itr::AsyncCollector)
    itr.ntasks = verify_ntasks(itr.enumerator, itr.ntasks)
    itr.batch_size = verify_batch_size(itr.batch_size)
    if itr.batch_size !== nothing
        exec_func = batch -> begin
            # extract indexes from the input tuple
            batch_idxs = map(x->x[1], batch)

            # and the args tuple....
            batched_args = map(x->x[2], batch)

            results = f(batched_args)
            foreach(x -> (itr.results[batch_idxs[x[1]]] = x[2]), enumerate(results))
        end
    else
        exec_func = (i,args) -> (itr.results[i]=itr.f(args...))
    end
    chnl, worker_tasks = setup_chnl_and_tasks((i,args) -> (itr.results[i]=itr.f(args...)), itr.ntasks, itr.batch_size)
    return AsyncCollectorState(chnl, worker_tasks, start(itr.enumerator))
end

function done(itr::AsyncCollector, state::AsyncCollectorState)
    if !isopen(state.chnl) || done(itr.enumerator, state.enum_state)
        close(state.chnl)

        # wait for all tasks to finish
        foreach(x->(v=wait(x); isa(v, Exception) && throw(v)), state.worker_tasks)
        empty!(state.worker_tasks)
        return true
    else
        return false
    end
end

function next(itr::AsyncCollector, state::AsyncCollectorState)
    if itr.nt_check && (length(state.worker_tasks) < itr.ntasks())
        start_worker_task!(state.worker_tasks, itr.f, state.chnl)
    end

    # Get index and mapped function arguments from enumeration iterator.
    (i, args), state.enum_state = next(itr.enumerator, state.enum_state)
    put!(state.chnl, (i, args))

    return (nothing, state)
end

"""
    AsyncGenerator(f, c...; ntasks=0, batch_size=nothing) -> iterator

Apply `f` to each element of `c` using at most `ntasks` asynchronous tasks.

Keyword args `ntasks` and `batch_size` have the same behavior as in
[`asyncmap`](@ref). If `batch_size` is specified, `f` must
be a function which operates on an array of argument tuples.

!!! note
    `collect(AsyncGenerator(f, c...; ntasks=1))` is equivalent to
    `map(f, c...)`.
"""
mutable struct AsyncGenerator
    collector::AsyncCollector
end

function AsyncGenerator(f, c...; ntasks=0)
    AsyncGenerator(AsyncCollector(f, Dict{Int,Any}(), c...; ntasks=ntasks))
end

mutable struct AsyncGeneratorState
    i::Int
    collector_state::AsyncCollectorState
end

start(itr::AsyncGenerator) = AsyncGeneratorState(0, start(itr.collector))

# Done when source async collector is done and all results have been consumed.
function done(itr::AsyncGenerator, state::AsyncGeneratorState)
    done(itr.collector, state.collector_state) && isempty(itr.collector.results)
end

function next(itr::AsyncGenerator, state::AsyncGeneratorState)
    state.i += 1

    results_dict = itr.collector.results
    while !haskey(results_dict, state.i)
        if done(itr.collector, state.collector_state)
            # `done` waits for async tasks to finish. if we do not have the index
            # we are looking for, it is an error.
            !haskey(results_dict, state.i) && error("Error processing index ", i)
            break;
        end
        _, state.collector_state = next(itr.collector, state.collector_state)
    end
    r = results_dict[state.i]
    delete!(results_dict, state.i)

    return (r, state)
end

# pass-through iterator traits to the iterable
# on which the mapping function is being applied
iteratorsize(itr::AsyncGenerator) = iteratorsize(itr.collector.enumerator)
size(itr::AsyncGenerator) = size(itr.collector.enumerator)
length(itr::AsyncGenerator) = length(itr.collector.enumerator)

"""
    asyncmap!(f, results, c...; ntasks=0, batch_size=nothing)

Like [`asyncmap`](@ref), but stores output in `results` rather than
returning a collection.
"""
function asyncmap!(f, r, c1, c...; ntasks=0, batch_size=nothing)
    foreach(identity, AsyncCollector(f, r, c1, c...; ntasks=ntasks, batch_size=batch_size))
    c
end

max_stars_repo_name: visr/julia

max_issues_repo_name: visr/julia

max_forks_repo_name: visr/julia


text:
# Testing part 2

include("../scripts-part2/elastic_wave_3D_testing.jl")

# Reference test using ReferenceTests.jl
"Compare all dict entries"
comp(d1, d2) = keys(d1) == keys(d2) && all([ isapprox(v1, v2; atol = 1e-5) for (v1,v2) in zip(values(d1), values(d2))])
inds = Int.(ceil.(LinRange(1, length(xc), 12)))
d = Dict(:X=> xc[inds], :H=>P[inds, inds, 15])

@testset "Reference Test - Elastic Wave 3D" begin
    @test_reference "reftest-files/test_2.bson" d by=comp
end

ext: jl


text:
struct HMC{TV, TGV, TF<:AbstractFloat, TM} <: SecondOrderMetropolisSampler
    V::TV
    ∇V!::TGV
    β::TF
    M::TM
    Δt::TF
    nΔt::Int
    sqrtM::TM
    sqrtβ::TF
end

"""
    HMC(V, ∇V!, β, M, Δt, nΔt)

Set up the HMC sampler for Boltzmann.

### Fields

* V     - Potential
* ∇V!   - In place gradient of the potential
* β     - Inverse temperature
* Δt    - Time step
* nΔt   - Number of time steps to use in each Verlet run
"""
function HMC(V::TV, ∇V!::TGV, β::TF, M::TM, Δt::TF, nΔt::Int) where {TV, TGV, TF<:AbstractFloat,TM}
    sqrtM = sqrt.(M);
    sqrtβ = sqrt(β);
    return HMC(V, ∇V!, β, M, Δt, nΔt, sqrtM,sqrtβ);
end

mutable struct HMCState{Tq, TF<:AbstractFloat} <:SecondOrderMetropolisSamplerState
    x::Tq
    x_proposal::Tq
    V::TF
    V_proposal::TF
    E::TF
    E_proposal::TF
    ∇V::Tq
    p::Tq
    p_mid::Tq
    invMp::Tq
    accept::Int
end

function InitState!(x₀, sampler::HMC)
    V = sampler.V(x₀);
    return HMCState(x₀, copy(x₀), V, V, V, V, similar(x₀), similar(x₀), similar(x₀),similar(x₀), Int(0));
end

function InitState(x₀, sampler::HMC)
    V = sampler.V(x₀);
    return HMCState(copy(x₀), copy(x₀), V, V, V, V, similar(x₀), similar(x₀), similar(x₀),similar(x₀), Int(0));
end

#
# function Verlet!(q, p, ∇V!, M, Δt, nΔt, p_mid, ∇V)
function Verlet!(state::HMCState, sampler::HMC)

    sampler.∇V!(state.∇V, state.x_proposal);
    for _ in 1:sampler.nΔt
        @. state.p_mid = state.p - 0.5 * sampler.Δt * state.∇V;
        @. state.x_proposal = state.x_proposal + sampler.Δt * state.p_mid/sampler.M;
        sampler.∇V!(state.∇V, state.x_proposal);
        @. state.p = state.p_mid - 0.5 * sampler.Δt * state.∇V;
    end
    state
end

function UpdateState!(state::HMCState, sampler::HMC)

    # prepare initial state
    @. state.x_proposal = state.x;
    @. state.p = (sampler.sqrtM/sampler.sqrtβ)  * randn();
    @. state.invMp = state.p / sampler.M;
    state.E = state.V + 0.5 * (state.p ⋅ state.invMp);
    # run Verlet integrator
    Verlet!(state, sampler);
    # compute energy
    state.V_proposal = sampler.V(state.x_proposal);
    @. state.invMp = state.p / sampler.M;
    state.E_proposal = state.V_proposal + 0.5 * (state.p ⋅ state.invMp);
    # accept/reject
    a = min(1, exp( sampler.β *(state.E - state.E_proposal)));
    if rand()<a
        @. state.x = state.x_proposal;
        state.V = state.V_proposal;
        state.accept = 1;
    else
        state.accept = 0;
    end
    state
end

max_stars_repo_name: cortner/JuBasicMD

max_issues_repo_name: cortner/JuBasicMD

max_forks_repo_name: cortner/JuBasicMD


text:
mw(model::EoSModel) = paramvals(model.params.Mw)

function group_molecular_weight(groups::GroupParam,mw,z = @SVector [1.])
    n = groups.n_flattenedgroups
    res = zero(first(z))
    Σz = sum(z)
    @inbounds for i in 1:length(groups.components)
        ni = n[i]
        gi = groups.i_groups[i]
        mwi = zero(res)
        for idx in 1:length(gi)
            j = gi[idx]
            mwi += mw[j]*ni[j]
        end
        res +=z[i]*mwi
    end
    return 0.001*res/Σz
end

comp_molecular_weight(mw,z = @SVector [1.]) = 0.001*mapreduce(*,+,mw,z)

include("initial_guess.jl")
include("differentials.jl")
include("VT.jl")
include("property_solvers/property_solvers.jl")
include("pT.jl")
include("unitful_base.jl")
include("unitful_methods.jl")


text:
"""
    nested_step!(e<:GMC_NS_Ensemble)

Step and update ensemble e by sampling the e.contour-bounded prior mass via Galilean Monte Carlo.

"""

function nested_step!(e::GMC_NS_Ensemble, tuners::Dict{Int64,τ_PID}, cache::Union{GMC_NS_Model,Nothing})
    N = length(e.models); N′=N #number of sample models/particles on the posterior surface
    i = length(e.log_Li) #iterate number, index for last values
    j = i+1 #index for newly pushed values

    e.contour, least_likely_idx = findmin([model.log_Li for model in e.models])
    Li_model = e.models[least_likely_idx]
    m = deserialize(Li_model.path)
    t=tuners[Li_model.trajectory]

    #SELECT NEW MODEL, SAVE TO ENSEMBLE DIRECTORY, CREATE RECORD AND PUSH TO ENSEMBLE

    cache===nothing ? (model_selected=false; cache_insert=false) : (model_selected=true; cache_insert=true)
    samples=1; sample_limit=Int64(floor(length(e.models)*e.GMC_exhaust_σ))

    while !model_selected && samples<=sample_limit
        if t.τ > e.GMC_τ_death && m.i <= e.GMC_chain_κ #while the particle's timestep is greater than the minimum to keep a chain alive, and the chain is below the maximum allowed length, sample by advancing it along a galilean trajectory
            model_selected,cache = galilean_search!(m,e,t,least_likely_idx,samples,sample_limit,cache) #ensemble is modified by galilean_search! directly
        elseif (N-1)<e.GMC_Nmin #if removing the particle would put the sample under the minimum size, resample by diffusion from an existing particle
            model_selected=diffusion_search!(e,tuners,least_likely_idx,N,e.GMC_τ_death) #ensemble modified by diffusion_search!
            samples+=1
        else #if not, push to the posterior samples and let the trajectory chain die
            model_selected=true
            e.sample_posterior && push!(e.posterior_samples, e.models[least_likely_idx])#if sampling posterior, push the model record to the ensemble's posterior samples vector
            deleteat!(e.models, least_likely_idx) #remove least likely model record from the active particles in the ensemble
            N′-=1 #update for trapezoidal log_wi calc
        end
    end

    if cache_insert #in this case a cached precalculated reflection model can be used
        insert_reflection!(e, cache) #insert the model into the ensemble
        process_report!(t, true)
        cache=nothing #reset the cache
    end

    samples==sample_limit && (return 1, cache) #failed to find a new sample after sampling up to the limit, return an error code

    #UPDATE ENSEMBLE QUANTITIES
    push!(e.log_Li, e.contour) #log likelihood of the least likely model - the current ensemble ll contour at Xi
    push!(e.log_Xi, -i/N) #log Xi - crude estimate of the iterate's enclosed prior mass
    push!(e.log_wi, logsubexp(e.log_Xi[i], -j/N′) - log(2)) #log width of prior mass spanned by the last step-trapezoidal approx
    push!(e.log_Liwi, lps(e.log_Li[j],e.log_wi[j])) #log likelihood + log width = increment of evidence spanned by iterate
    push!(e.log_Zi, logaddexp(e.log_Zi[i],e.log_Liwi[j]))    #log evidence
    #information- dimensionless quantity. cf. MultiNest @ https://github.com/farhanferoz/MultiNest/blob/master/MultiNest_v3.12/nested.F90- MultiNest now gives info only at the final step
    Hj=lps(
        (exp(lps(e.log_Liwi[j],-e.log_Zi[j])) * e.log_Li[j]), #information contribution of this step
        (exp(lps(e.log_Zi[i],-e.log_Zi[j])) * lps(e.Hi[i],e.log_Zi[i])), #rescale last information by evidence
        -e.log_Zi[j])
    Hj === -Inf ? push!(e.Hi,0.) : push!(e.Hi, Hj) #prevent problems from early strings of models with -Inf log likelihoods

    return 0, cache
end


text:
import Base: download
import Dates: unix2datetime, now

"""
    download(rf::RemoteFile;
             quiet::Bool=false,
             verbose::Bool=false,
             force::Bool=false,
             retries::Int=0)

Download `rf`.

- `quiet`: Do not print messages.
- `verbose`: Print all messages.
- `force`: Force download and overwrite existing files.
- `force_update`: Overwrite existing files even if they are equal.
- `retries`: Override the number of retries in `rf` if `retries != 0`
"""
function download(rf::RemoteFile; kwargs...)
    for (i, backend) in enumerate(rf.backends)
        try
            if i == 1
                download(backend, rf; kwargs...)
            else
                download(backend, rf; retries=1, kwargs...)
            end
            return
        catch err
            err isa DownloadError && rf.try_backends && continue
            rethrow(err)
        end
    end
    throw(DownloadError("All available backends failed."))
end

function download(backend, rf::RemoteFile; verbose::Bool=false, quiet::Bool=false,
    force::Bool=false, force_update::Bool=false, retries::Int=0)
    file = joinpath(rf.dir, rf.file)
    retries = retries != 0 ? retries : rf.retries

    if isfile(file) && !force
        if !isoutdated(rf)
            verbose && @info "File '$(rf.file)' is up-to-date."
            return
        end
    end

    tempfile = tempname()

    if !isdir(rf.dir)
        mkpath(rf.dir)
    end

    !quiet && @info "Downloading file '$(rf.file)' from '$(rf.uri)' with $(nameof(backend))."
    tries = 0
    success = false
    while tries < retries
        tries += 1
        try
            download(backend, string(rf.uri), tempfile, verbose=verbose)
            success = true
        catch err
            if err isa DownloadError
                if !quiet
                    @warn "Download failed with error: $(err.msg)"
                    @info "Retrying in $(rf.wait) seconds."
                end
                sleep(rf.wait)
            else
                rethrow(err)
            end
        end
    end

    if success
        update = true
        if isfile(file) && !force
            if !force_update && samecontent(tempfile, file)
                update = false
                verbose && @info "File '$(rf.file)' has not changed. Update skipped."
            else
                verbose && @info "File '$(rf.file)' was successfully updated."
            end
        else
            verbose && @info "File '$(rf.file)' was successfully downloaded."
        end
        update && mv(tempfile, file, force=true)
    else
        if rf.failed == :warn && isfile(file)
            @warn "Download of '$(rf.uri)' failed after $(retries) retries. Local file was not updated."
        elseif rf.failed == :error || (rf.failed == :warn && !isfile(file))
            throw(DownloadError("Download of '$(rf.uri)' failed after $(retries) retries."))
        end
    end
    rm(tempfile, force=true)
end

function samecontent(file1, file2)
    h1 = hash(open(read, file1, "r"))
    h2 = hash(open(read, file2, "r"))
    h1 == h2
end


text:
using Telegraph
using Test

@testset "Telegraph.jl" begin
    λ₁ = .1
    λ₂ = .3
    Δt = 1e-2
    λ₁Δt = λ₁*Δt
    λ₂Δt = λ₂*Δt
    s = rand([0,1])
    # Single test
    T = update(s,λ₁Δt,λ₂Δt)
    @test (T == 1) | (T == 0) # test output is either a 1 or a 0

    # Several test
    t = 100
    t⃗ = range(0,t,step=Δt)
    st⃗ = size(t⃗,1)
    st = rand([0,1])
    state = []
    push!(state,st)
    for i in 2:st⃗
            st = update(st,λ₁Δt,λ₂Δt)
            push!(state,st)
        end

    @test state isa Vector
    @test minimum(state) == 0
    @test maximum(state) == 1
    @test size(state,1) == st⃗

    # test all indices
    gi = get_all_indices(state)
    @test gi isa Vector
    @test gi[end] == st⃗
    @test gi[1] == 1

    # test indices
    gibm = get_indices(state, t⃗)
    @test gibm isa telegraph_points

    tele_time_dist = get_times(gibm,t⃗,st)
    @test tele_time_dist isa telegraph_time_dist

end


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

# UNIT WEIGHTS

struct UnitWeights{T<:Real} <: AbstractWeights{Int, T, V where V<:Vector{T}}
    len::Int
end

Base.eltype(wv::UnitWeights{T}) where T = T
Base.sum(wv::UnitWeights{T}) where T = convert(T, length(wv))
Base.isempty(wv::UnitWeights) = iszero(wv.len)
Base.length(wv::UnitWeights) = wv.len
Base.size(wv::UnitWeights) = Tuple(length(wv))

Base.convert(::Type{Vector}, wv::UnitWeights{T}) where {T} = ones(T, length(wv))

@propagate_inbounds function Base.getindex(wv::UnitWeights{T}, i::Integer) where T
    @boundscheck checkbounds(wv, i)
    one(T)
end

@propagate_inbounds function Base.getindex(wv::UnitWeights{T}, i::AbstractArray{<:Int}) where T
    @boundscheck checkbounds(wv, i)
    UnitWeights{T}(length(i))
end

Base.getindex(wv::UnitWeights{T}, ::Colon) where {T} = UnitWeights{T}(wv.len)

@inline function varcorrection(w::UnitWeights, corrected::Bool=false)
    corrected ? (1 / (w.len - 1)) : (1 / w.len)
end

function Base.sum(A::AbstractArray, w::UnitWeights; dims::Union{Colon, Int}=:)
    a = (dims === :) ? length(A) : size(A, dims)
    a != length(w) && throw(DimensionMismatch("Inconsistent array dimension."))
    return sum(A, dims=dims)
end

function mean(A::AbstractArray, w::UnitWeights; dims::Union{Colon, Int}=:)
    a = (dims === :) ? length(A) : size(A, dims)
    a != length(w) && throw(DimensionMismatch("Inconsistent array dimension."))
    return mean(A, dims=dims)
end

Base.isequal(x::UnitWeights, y::UnitWeights) = isequal(x.len, y.len)
Base.:(==)(x::UnitWeights, y::UnitWeights)   = (x.len == y.len)

max_stars_repo_name: lbittarello/Microeconometrics

max_issues_repo_name: lbittarello/Microeconometrics

max_forks_repo_name: lbittarello/Microeconometrics


text:
"""
    PetersFlap <: NoStateModel

Linear, steady-state, two-dimensional control surface model with parameters
``c_{n,\\delta}``, ``c_{a,\\delta}``, and ``c_{m,\\delta}``.
"""
struct PetersFlap <: NoStateModel end

"""
    PetersFlap()

Initialize an object of type [`PetersFlap`](@ref)
"""
PetersFlap()

# --- Traits --- #

number_of_parameters(::Type{PetersFlap}) = 3

# --- Convenience Methods --- #

function set_parameters!(p, model::PetersFlap; cnd, cad, cmd)

    d
    ε
    δ

    p[1] = cnd
    p[2] = cad
    p[3] = cmd

    return p
end

function separate_parameters(model::PetersFlap, p)

    return (cnd=p[1], cad=p[2], cmd=p[3])
end

# --- Internal Methods --- #

function petersflap_loads(b, U, ρ, cnδ, caδ, cmδ, δ)

    T = petersflap_angular_deflection_coefficients(θ, n)
    h = T*δ
    Σ = (1:N)'*h
    L = 2*pi*ρ*U^2*Σ
    D =
    M = 2*pi*ρ*U^2*h[1]

    n = 2
    hn = petersflap_displacement_field
    2*sin(n*θ)/(n*pi) * ε + b/pi*((sin(n-1)*θ)/(n-1) - (1+d)*sin(n*θ)/n + sin((n+1)*θ)/(n+1))*δ

    L = ρ*U^2*b*cnδ*δ
    M = ρ*U^2*b^2*cmδ*δ
    A/(2*pi*ρ*b) = -(v + d*ω - λ0 - u*α0)^2 - (v + d*ω - λ0 - u*α0)*u/b*Σ +
        ω1*U2*h[1] -1/2*vdot*h[1] + b/4*ωdot*h[2]

    A/(2*pi*ρ*b) = -(v + d*ω - λ0 - u*α0)^2 + (Σ/b*u*v + (Σ/b*d + h[1])*u*ω - Σ/b*u*λ0 - u^2*α0*Σ/b - 1/2*vdot*h[1]) +
         + b/4*ωdot*h[2]

    return L, D, M
end

function petersflap_displacement_field(b, θ, ε, δ, n)
    if n == 0
        return θ/pi*ε + b/pi*(sin(θ) - (1+d)/2*θ)*δ
    elseif n == 1
        return 2*sin(θ)/π*ε + b/pi*(θ - (1+d)*sin(θ) + sin(2*θ)/2)*δ
    else
        return 2*sin(n*θ)/(n*pi) * ε + b/pi*((sin(n-1)*θ)/(n-1) -
            (1+d)*sin(n*θ)/n + sin((n+1)*θ)/(n+1))*δ
    end
end

function petersflap_linear_deflection_coefficients(θ, n)
    return SVector(ntuple(n -> petersflap_linear_deflection_coefficient(θ, n), n)...)
end

function petersflap_linear_deflection_coefficient(θ, n)
    if iszero(n)
        return θ/pi
    else
        return 2*sin(θ)/(n*pi)
    end
end

function petersflap_angular_deflection_coefficients(θ, n)
    return SVector(ntuple(n -> petersflap_angular_deflection_coefficient(θ, n), n)...)
end

function petersflap_angular_deflection_coefficient(θ, n)
    if iszero(n)
        return sin(θ) - θ*cos(θ)
    else
        return b/pi*sin((n+1)*θ)/(n+1) - 2/n*cos(θ)*sin(n*θ) + sin((n-1)*θ)/(n-1)
    end
end


text:
@testset "Solvers" begin
  @testset "CookieCutter" begin
    problem = SimulationProblem(RegularGrid(100,100), (:facies => Int, :property => Float64), 3)
    solver = CookieCutter(DummySimSolver(:facies=>NamedTuple()),
                          Dict(0=>DummySimSolver(), 1=>DummySimSolver()))

    @test sprint(show, solver) == "CookieCutter"
    @test sprint(show, MIME"text/plain"(), solver) == "CookieCutter\n  └─facies ⇨ DummySimSolver\n    └─0 ⇨ DummySimSolver\n    └─1 ⇨ DummySimSolver\n"

    Random.seed!(1234)
    solution = solve(problem, solver)

    if visualtests
      @plottest plot(solution,size=(800,600)) joinpath(datadir,"cookiecutter.png") !isCI
    end
  end

  @testset "SeqSim" begin
    Random.seed!(1234)
    sdata = georef((z=rand(100),), 100*rand(2,100))
    sgrid = RegularGrid(100,100)

    prob1 = SimulationProblem(sgrid, :z => Float64, 3)
    prob2 = SimulationProblem(sdata, sgrid, :z, 3)

    solver = SeqSim(:z => (estimator=DummyEstimator(),
                           neighborhood=BallNeighborhood(10.),
                           minneighbors=1, maxneighbors=10,
                           marginal=Normal(), path=LinearPath(),
                           mapping=NearestMapping()))

    Random.seed!(1234)
    usol = solve(prob1, solver)
    csol = solve(prob2, solver)

    if visualtests
      @plottest plot(usol,size=(900,300)) joinpath(datadir,"seqsim.png") !isCI
    end
  end

  @testset "PointwiseLearn" begin
    # synthetic data
    Random.seed!(1234)
    f(x,y) = sin(4*(abs(x)+abs(y))) < 0 ? 1 : 0
    X = [sin(i/10) for i in 1:100, j in 1:100]
    Y = [sin(j/10) for i in 1:100, j in 1:100]
    Z = categorical(f.(X,Y))
    ϵ₁ = 0.1randn(Float64, size(X))
    ϵ₂ = 0.1randn(Float64, size(Y))

    # source and target data
    S = georef((X=X,Y=Y,Z=Z))
    T = georef((X=X+ϵ₁,Y=Y+ϵ₂))

    # view versions
    inds = shuffle(1:nelms(S))
    Sv = view(S, inds)
    Tv = view(T, inds)

    # classification task
    𝓉 = ClassificationTask((:X,:Y), :Z)

    # learning problems
    𝒫₁ = LearningProblem(S, T, 𝓉)
    𝒫₂ = LearningProblem(Sv, Tv, 𝓉)

    # pointwise solver
    ℒ = PointwiseLearn(dtree)

    R₁ = solve(𝒫₁, ℒ)
    R₂ = solve(𝒫₂, ℒ)

    # error is small
    @test mean(S[:Z] .!= R₁[:Z]) < 0.15
    @test mean(Sv[:Z] .!= R₂[:Z]) < 0.15

    if visualtests
      for (i,s) in enumerate([(S,R₁), (Sv,R₂)])
        @plottest begin
          p1 = plot(s[1], (:Z,))
          p2 = plot(s[2])
          plot(p1, p2, size=(800,400))
        end joinpath(datadir,"pointlearn$i.png") !isCI
      end
    end
  end
end


text:
using FromFile
@from "ProgramConstants.jl" import CONST_TYPE

# Define a serialization format for the symbolic equations:
mutable struct Node
    #Holds operators, variables, constants in a tree
    degree::Int #0 for constant/variable, 1 for cos/sin, 2 for +/* etc.
    constant::Bool #false if variable
    val::CONST_TYPE
    # ------------------- (possibly undefined below)
    feature::Int #Either const value, or enumerates variable.
    op::Int #enumerates operator (separately for degree=1,2)
    l::Node
    r::Node

    Node(d::Int, c::Bool, v::CONST_TYPE) = new(d, c, v)
    Node(d::Int, c::Bool, v::CONST_TYPE, f::Int) = new(d, c, v, f)
    Node(d::Int, c::Bool, v::CONST_TYPE, f::Int, o::Int, l::Node) = new(d, c, v, f, o, l)
    Node(d::Int, c::Bool, v::CONST_TYPE, f::Int, o::Int, l::Node, r::Node) = new(d, c, v, f, o, l, r)
end

Node(val::CONST_TYPE) =                                                     Node(0, true,                       val                                     ) #Leave other values undefined
"""
    Node(feature::Int)

Create a variable node using feature `feature::Int`
"""
Node(feature::Int) =                                                        Node(0, false, convert(CONST_TYPE, 0f0), feature                            )
"""
    Node(op::Int, l::Node)

Apply unary operator `op` (enumerating over the order given) to `Node` `l`
"""
Node(op::Int, l::Node) =                                                    Node(1, false, convert(CONST_TYPE, 0f0),       0,      op,        l         )
"""
    Node(op::Int, l::Union{AbstractFloat, Int})

Short-form for creating a scalar/variable node, and applying a unary operator
"""
Node(op::Int, l::Union{AbstractFloat, Int}) =                               Node(1, false, convert(CONST_TYPE, 0f0),       0,      op,  Node(l)         )
"""
    Node(op::Int, l::Node, r::Node)

Apply binary operator `op` (enumerating over the order given) to `Node`s `l` and `r`
"""
Node(op::Int, l::Node, r::Node) =                                           Node(2, false, convert(CONST_TYPE, 0f0),       0,      op,        l,       r)
"""
    Node(op::Int, l::Union{AbstractFloat, Int}, r::Node)

Short-form to create a scalar/variable node, and apply a binary operator
"""
Node(op::Int, l::Union{AbstractFloat, Int}, r::Node) =                      Node(2, false, convert(CONST_TYPE, 0f0),       0,      op,  Node(l),       r)
"""
    Node(op::Int, l::Node, r::Union{AbstractFloat, Int})

Short-form to create a scalar/variable node, and apply a binary operator
"""
Node(op::Int, l::Node, r::Union{AbstractFloat, Int}) =                      Node(2, false, convert(CONST_TYPE, 0f0),       0,      op,        l, Node(r))
"""
    Node(op::Int, l::Union{AbstractFloat, Int}, r::Union{AbstractFloat, Int})

Short-form for creating two scalar/variable node, and applying a binary operator
"""
Node(op::Int, l::Union{AbstractFloat, Int}, r::Union{AbstractFloat, Int}) = Node(2, false, convert(CONST_TYPE, 0f0),       0,      op,  Node(l), Node(r))
"""
    Node(val::AbstractFloat)

Create a scalar constant node
"""
Node(val::AbstractFloat) =                                                  Node(convert(CONST_TYPE, val))
"""
    Node(var_string::String)

Create a variable node, using the format `"x1"` to mean feature 1
"""
Node(var_string::String) =                                                  Node(parse(Int, var_string[2:end]))
"""
    Node(var_string::String, varMap::Array{String, 1})

Create a variable node, using a user-passed format
"""
Node(var_string::String, varMap::Array{String, 1}) =                        Node([i for (i, _variable) in enumerate(varMap) if _variable==var_string][1]::Int)

# Copy an equation (faster than deepcopy)
function copyNode(tree::Node)::Node
   if tree.degree == 0
       if tree.constant
           return Node(tree.val)
        else
           return Node(tree.feature)
        end
   elseif tree.degree == 1
       return Node(tree.op, copyNode(tree.l))
    else
        return Node(tree.op, copyNode(tree.l), copyNode(tree.r))
   end
end


text:
using GeoPhyInv
using BenchmarkTools
using LinearAlgebra
using Test
GInterp=GeoPhyInv.Interpolation

function dptest(nx,nz,mx,mz,Battrib)
	pa=GInterp.Kernel([nx, nz], [mx, mz],Battrib)
	ny=randn(length(nz), length(nx));
	my=zeros(length(mz), length(mx));
	# interp
	@time GInterp.interp_spray!(ny, my, pa, :interp)
	myp = randn(size(my));
	# spray
	nyp=zeros(length(nz), length(nx));
	@time GInterp.interp_spray!(nyp, myp, pa, :spray)
	# dot product test
	@test LinearAlgebra.dot(my, myp) ≈ LinearAlgebra.dot(ny, nyp)
end

@testset "GInterp to a larger grid using  nx==1 or nz==1" begin
	nx=range(1.,step=20.,length=1);
	mx=range(1.,stop=40.,length=10);
	nz=range(1.,stop=20.,length=30);
	mz=range(4.,stop=16.,length=30);

	dptest(nx,nz,mx,mz,:B1)

	nx=range(4.,stop=16.,length=20);
	mx=range(1.,stop=30.,length=30);
	nz=range(1.,step=5.,length=1);
	mz=range(1.,stop=10.,length=66);

	dptest(nx,nz,mx,mz,:B1)
end

@testset "GInterp on same grid for nx==1 or nz==1" begin
	nx=range(1.,stop=1.,length=1);
	mx=range(1.,stop=1.,length=1);
	nz=range(1.,stop=20.,length=30);
	mz=range(4.,stop=16.,length=30);

	dptest(nx,nz,mx,mz,:B1)

	nx=range(4.,stop=16.,length=20);
	mx=range(1.,stop=30.,length=30);
	nz=range(1.,stop=1.,length=1);
	mz=range(1.,stop=1.,length=1);

	dptest(nx,nz,mx,mz,:B1)
end

@testset "no extrapolations" begin
	nx=range(21,stop=100,length=80);
	mx=range(1,stop=100,length=100);
	ny=ones(length(nx));
	c=randn()
	my=c.*ones(length(mx))

	pa=GInterp.P_core([nx], [mx], :B1)

	GInterp.interp_spray!(ny, my, pa, :interp)
	@test all(my[1:20] .== c)

	nx=range(3,stop=7,length=64);
	nz=range(3,stop=7,length=95);
	mx=range(1,stop=10,length=10);
	mz=range(1,stop=10,length=10);
	ny=ones(length(nz), length(nx));
	c=randn()
	my=c.*ones(length(mz), length(mx))

	pa=GInterp.P_core([nx, nz], [mx, mz], :B1)

	GInterp.interp_spray!(ny, my, pa, :interp)

	@test all(my[my .≠1] .== c)

end

# also testing behaviour when out of bounds!!
nx=range(1,stop=20,length=100);
mx=range(3,stop=25,length=200);
nz=range(1,stop=20,length=30);
mz=range(5,stop=27,length=40);

println("=====================================")
for Battrib in [:B1, :B2]
	global nx, mx

	pa=GInterp.Kernel([nx], [mx], Battrib)
	println("=====================================")
	## 1D
	ny=randn(length(nx));
	# interpolate (my=f(ny))
	my=zeros(length(mx));
	@time GInterp.interp_spray!(ny, my, pa, :interp)

	myp = randn(length(my));

	# spray (nyp=fˣ(myp))
	nyp=zeros(length(nx));
	@time GInterp.interp_spray!(nyp, myp, pa, :spray)

	# dot product test
	@test LinearAlgebra.dot(my, myp) ≈ LinearAlgebra.dot(ny, nyp)
end

@testset "GInterp on the same grid should not change for B1" begin
	Battrib=:B1
	pa=GInterp.Kernel([nx, nz], [nx, nz], Battrib)
	println("=====================================")
	ny=randn(length(nz), length(nx));
	my=zeros(length(nz), length(nx));
	@time GInterp.interp_spray!(ny, my, pa, :interp)
	@test ny≈my

	nyvec=vcat(vec(ny),vec(ny))
	myvec=zeros(length(my)*2)
	@time GInterp.interp_spray!(nyvec, myvec, pa, :interp, 2)
	@test nyvec≈myvec

	myp = randn(size(my));
	nyp=zeros(length(nz), length(nx));
	@time GInterp.interp_spray!(nyp, myp, pa, :spray)
	@test LinearAlgebra.dot(my, myp) ≈ LinearAlgebra.dot(ny, nyp)
end

for Battrib in [:B1, :B2]
	global nx, nz, mx, mz
	pa=GInterp.Kernel([nx, nz], [mx, mz], Battrib)
	## 2D
	ny=randn(length(nz), length(nx));
	# interpolate
	my=zeros(length(mz), length(mx));
	@time GInterp.interp_spray!(ny, my, pa, :interp)

	myp = randn(size(my));

	# spray
	nyp=zeros(length(nz), length(nx));
	@time GInterp.interp_spray!(nyp, myp, pa, :spray)

	# dot product test
	@test LinearAlgebra.dot(my, myp) ≈ LinearAlgebra.dot(ny, nyp)

	nmod=3
	ny=randn(length(nz)*length(nx)*nmod);
	# interpolate
	my=zeros(length(mz)*length(mx)*nmod);
	@time GInterp.interp_spray!(ny, my, pa, :interp,nmod)

	myp = randn(size(my));

	# spray
	nyp=zeros(length(nz)*length(nx)*nmod);
	@time GInterp.interp_spray!(nyp, myp, pa, :spray,nmod)

	# dot product test
	@test LinearAlgebra.dot(my, myp) ≈ LinearAlgebra.dot(ny, nyp)

end


text:
function fct(a,b,c)

    @show a,b,c

    function fct1(v)
        a,b,c  = (a,b,c) .+ v
    end

    fct1(4)

    @show a, b, c
end

fct(1,2,3)


text:
# using SafeTestsets
using MinimaxAI
using Test
# MinimaxAI.playBones(ai = 2,N=20)

@testset "Seeing if it works" begin
    @test playBones(N = 2, ai =2) == nothing
    @test playBones(N = 25, ai = 2) == nothing
    @test playTicTacToe(ai = 2) == nothing
    @test playConnect4(ai = 2, method = "minimaxAB", abDepth = 3) == nothing
    @test playBones(N = 25, ai = 2, method = "minimaxAB", abDepth = 3) == nothing
    @test playTicTacToe(ai = 2 , method = "minimaxAB", abDepth = 5) == nothing
end

# @safetestset "Generic tests" begin include("minimax_tests.jl") end


text:
function test_Species()
    @test_throws TypeError Species(QE, ME, z->12, electroncollisionfrequency)
    @test_throws TypeError Species(QE, ME, z -> waitprofile(z, 75, 0.32), z->12)
end

function test_BField()
    bfield = BField(50e-6, π/2, 3π/2)
    @test isbits(bfield)
    @test (@test_logs (:warn, "BField magnitude of exactly 0 is not supported."*
                       " Setting B = 1e-16.") BField(0, π/2, 3π/2).B) == 1e-16

    #==
    dip and az
    ==#
    B = 50000e-9
    true_dip = deg2rad(-78)
    true_az = deg2rad(240)
    bfield = BField(B, true_dip, true_az)

    @test LMP.dip(bfield) ≈ true_dip
    @test mod2pi(LMP.azimuth(bfield)) ≈ true_az

    true_dip = deg2rad(62.4)
    true_az = deg2rad(18.3)
    bfield = BField(B, true_dip, true_az)

    @test LMP.dip(bfield) ≈ true_dip
    @test mod2pi(LMP.azimuth(bfield)) ≈ true_az

    #==
    direction cosines
    ==#
    # dip angle 90° up from horizontal
    bfield = BField(50e-6, π/2, 3π/2)
    @test bfield.dcn == -1

    # dip angle 90° down from horizontal
    bfield = BField(50e-6, -π/2, 3π/2)
    @test bfield.dcn == 1
end

function test_Ground()
    @test isbits(Ground(15, 0.001))

    @test LMP.GROUND[3] == Ground(10, 1e-4)
end

function test_waitprofile()
    zs = 50e3:2e3:86e3
    hbs = [(65, 0.2), (65, 0.25), (68, 0.25), (70, 0.25), (72, 0.3), (72, 0.35), (75, 0.35),
           (78, 0.35), (78, 0.4), (82, 0.5), (82, 0.55), (85, 0.6), (85, 0.65), (88, 0.8)]

    N(z, h′, β) = 1.43e13*exp(-0.15*h′)*exp((β - 0.15)*(z/1000 - h′))

    for hb in hbs, z in zs
        @test waitprofile(z, hb[1], hb[2]) ≈ N(z, hb[1], hb[2])
    end

    # Check if zero below cutoff height
    @test waitprofile(30e3, 60, 0.35; cutoff_low=40e3) == 0
    @test waitprofile(30e3, 60, 0.35; cutoff_low=40e3) isa Float64
    @test waitprofile(30_000, 60, 0.35; cutoff_low=40e3) isa Float64

    # Check default threshold
    @test waitprofile(110e3, 70, 0.5) == 1e12
end

@testset "Geophysics.jl" begin
    @info "Testing Geophysics"

    test_Species()
    test_BField()
    test_Ground()
    test_waitprofile()
    @test electroncollisionfrequency(110e3) isa Float64
    @test ioncollisionfrequency(110e3) isa Float64
end


text:
const _basis{N} = SVector{N,SVector{N,Float64}}

"""
    _basis{N}()

Returns a set of `N` `N`-dimensional zero vectors, which can serve to indicate that the basis
vectors of a dataset arecons missing, or that the structure is aperiodic.
"""
function _basis{N}() where N
    return zeros(SVector{N,SVector{N,Float64}})
end

"""
    _allsame(itr)

Returns `true` if all the elements of an iterator are identical.
"""
_allsame(itr) = all(x -> x == itr[1], itr)

"""
    _tovectors(M::AbstractMatrix)

Converts a matrix of basis vectors to a vector of vectors. This is the preferred representation of
basis vectors throughout Crystalline.jl.
"""
_tovectors(M::AbstractMatrix) = [M[:,n] for n in 1:size(M)[2]]

"""
    _tomatrix(V::AbstractVector{<:AbstractVector})

Converts a vector of basis vectors to a matrix. The preferred represenation of basis vectors in
Crystalline.jl is a vector of vectors, so this should only be used when matrix-specific operations
are needed.
"""
_tomatrix(V::AbstractVector{<:AbstractVector}) = hcat(V...)

"""
    _is_linearly_independent(vecs::AbstractMatrix{<:Number}) -> Bool

Returns `true` if a matrix consists of linearly independent vectors, `false` otherwise.

This test is performed by checking that the rank of the matrix is equal to the smallest dimension
of the matrix.
"""
function _is_linearly_independent(vecs::AbstractMatrix{<:Number})
    return LinearAlgebra.rank(vecs) == minimum(size(vecs))
end

"""
    _is_linearly_independent(vecs::NTuple{N,NTuple{M,<:Number}}) -> Bool
    _is_linearly_independent(vecs::NTuple{M,<:Number}...) -> Bool
    _is_linearly_independent(vecs::AbstractVector{<:AbstractVector{<:Number}}) -> Bool

Returns `true` if a set of vectors is linearly independent, false otherwise.
"""
function _is_linearly_independent(vecs::NTuple{N,NTuple{M,<:Number}}) where {M,N}
    # if M is greater than N, they cannot be linearly independent
    M < N && return false
    mat = [vecs[n][m] for m = 1:M, n = 1:N]
    return _is_linearly_independent(mat)
end

_is_linearly_independent(vecs::NTuple{M,<:Number}...) where M = _is_linearly_independent(vecs)

function _is_linearly_independent(vecs::AbstractVector{<:AbstractVector{<:Number}})
    # Get number of vectors
    N = length(vecs)
    # Get length of each vector in the vector
    M = let l = [length(v) for v in vecs]
        _allsame(l) ? l[1] : throw(ErrorException("Vectors have inconsistent dimensions"))
    end
    # if M is greater than N, they cannot be linearly independent
    M < N && return false
    # Convert to a matrix for rank checking
    mat = [vecs[m][n] for n = 1:N, m = 1:M]
    return _is_linearly_independent(mat)
end

"""
    _reducecoords(B::AbstractMatrix{<:Number}, v::AbstractVector{<:Number})

Converts a vector `v` from Cartesian coordinates to reduced coordinates relative to the basis
vectors `B` of a unit cell, expressed as a collection of vectors.
"""
function _reducecoords(b::AbstractVector{AbstractVector{<:Number}}, v::AbstractVector{<:Number})
    B = hcat(b...)
    return inv(B) * v
end

"""
    printbasis(io::IO, basis::AbstractVector{<:AbstractVector{<:Real}})

Prints the basis vectors in the following format:

```
Basis vectors:
a:    2.847216    0.000000    0.000000
b:   -1.423608    2.465762   -0.000000
c:    0.000000    0.000000    4.538054
```
"""
function printbasis(io::IO, basis::AbstractVector{<:AbstractVector{<:Real}})
    tostr(x::Number) = lpad(@sprintf("%f", x), 12)
    # Check if basis vectors are empty or zero
    if isempty(basis) || any(all.(iszero, v for v in basis))
        println(io, "Aperiodic (no basis vectors specified)")
    else
        println(io, "Basis vectors:")
        for (n,v) in enumerate(basis)
            # '`' is 0x60, and all letters come after this in Unicode
            # e.g. 'a' is 0x61, 'b' is 0x62, 'c' is 0x63...
            println(io, ('`' + n) * ':' * join(tostr.(v)))
        end
    end
    return nothing
end


text:
#=
# 7. Lid-driven cavity flow
In this notebook we will simulate the flow in a square cavity with a top moving wall.
=#
using ViscousFlow
#-
using Plots

#=
## Problem specification
Take $Re=100$ for example:
=#
Re = 100

#=
## Discretization
Note that the rectangle function used for making the cavity shape requires a specified half length. The
immersed boundary projection method for internal flow requires the size of the domain to
be at least a step size greater at the boundaries (i.e. `halflength + Δx`).
=#
Δt,Δx = setstepsizes(Re,gridRe=1.0)
halflength=0.5
domain_lim=halflength+1.01*Δx
xlim, ylim = (-domain_lim,domain_lim),(-domain_lim,domain_lim)

#=
## Cavity Geometry
A square cavity can be created using the `Rectangle()` function with the half length defined above.
The `shifted=true` argument ensures that points are not placed at the corners, where
they have ill-defined normal vectors.
=#
body = Rectangle(halflength,halflength,1.5*Δx,shifted=true)
plot(body,fillrange=nothing)

#=
## Boundary Condition at the moving wall
Assign velocity to the top boundary. The `LidDrivenCavity()` function can be used
to specify the velocity value at the top wall. Note : Non-dimensional velocity = 1
=#
m = ViscousFlow.LidDrivenCavity(1.0)

#=
## Construct the system structure
The last two input `flow_side` and `static_points` must specified so the default setting in the
 `NavierStokes()` function can be overwritten.

`static_points` is set to true because the cavity wall points are not actually moving.
=#
sys = NavierStokes(Re,Δx,xlim,ylim,Δt,body,m,flow_side = InternalFlow,static_points = true)

#=
Initialize
=#
u0 = newstate(sys)

#=
Set up integrator
=#
tspan = (0.0,10.0)
integrator = init(u0,tspan,sys)

#=
## Solve
=#
step!(integrator,10)

#=
## Examine
Plot the vorticity and streamlines
=#

plot(
plot(vorticity(integrator),sys,title="Vorticity",clim=(-10,10),color=:turbo,linewidth=1.5,ylim=ylim,levels=-6:0.5:5),
plot(streamfunction(integrator),sys,title="Streamfunction",color=:black,levels=vcat(0.009:0.01:0.11,0.1145,0.11468,0.11477),ylim=ylim)
   )

#=
Make a movie:
=#
sol = integrator.sol;
@gif for (u,t) in zip(sol.u,sol.t)
    plot(vorticity(u,sys,t),sys,clim=(-10,10),levels=range(-10,10,length=30),color=:turbo)
end every 5


text:
@testset "Geospatial" begin
  println(testname("Geospatial tests"))
  print("  [  ] Create a geospatial index ")
  result = geo_index(db, "geodd", "geoidx",
    "function(doc){if(doc.geometry&&doc.geometry.coordinates){st_index(doc.geometry);}}"
  )
  @test result["ok"] == true
  println("\r  [OK] Create a geospatial index")

  print("  [  ] Get geospatial index info ")
  # result = geo_indexinfo(db, "geodd", "geoidx")
  # @test haskey(result, "geo_index") == true
  # println("\r  [OK] Get geospatial index info")

  if geo_database != ""
    geodb = connectdb(cl, geo_database)
    print("  [  ] Radius geospatial query ")
    result = geo_query(geodb, "geodd", "geoidx";
      lat    = 42.357963,
      lon    = -71.063991,
      radius = 10000.0,
      limit  = 200)

    @test length(result["rows"]) == 200
    println("\r  [OK] Radius geospatial query")

    print("  [  ] Elliptic geospatial query ")
    result = geo_query(geodb, "geodd", "geoidx";
      lat    = 42.357963,
      lon    = -71.063991,
      rangex = 100.0,
      rangey = 500.0,
      limit  = 1)

    @test length(result["rows"]) == 1
    println("\r  [OK] Elliptic geospatial query")

    print("  [  ] Polygon geospatial query ")
    result = geo_query(geodb, "geodd", "geoidx";
      g="POLYGON ((-71.0537124 42.3681995 0,-71.054399 42.3675178 0,-71.0522962 42.3667409 0,-71.051631 42.3659324 0,-71.051631 42.3621431 0,-71.0502148 42.3618577 0,-71.0505152 42.3660275 0,-71.0511589 42.3670263 0,-71.0537124 42.3681995 0))")
    @test length(result["rows"]) == 2
    println("\r  [OK] Polygon geospatial query")

    print("  [  ] Radius geospatial query with skip, legacy format and stale=true ")
    result = geo_query(geodb, "geodd", "geoidx";
      skip   = 100,
      format = "legacy",
      nearest = true,
      stale  = true,
      lat    = 42.357963,
      lon    = -71.063991,
      radius = 10000.0,
      limit  = 200)
    @test result["type"] == "FeatureCollection"
    println("\r  [OK] Radius geospatial query with skip, legacy format and stale=true")

    print("  [  ] Radius geospatial query with bookmark ")
    result = geo_query(geodb, "geodd", "geoidx";
      skip   = 100,
      format = "legacy",
      nearest = true,
      stale  = true,
      lat    = 42.357963,
      lon    = -71.063991,
      radius = 10000.0,
      limit  = 200,
      bookmark = result["bookmark"])
    @test result["type"] == "FeatureCollection"
    println("\r  [OK] Radius geospatial query with bookmark")

    print("  [  ] Geospatial query with unknown format should fail ")
    @test_throws ErrorException geo_query(geodb, "geodd", "geoidx";
      format = "never-heard-of-it",
      lat    = 42.357963,
      lon    = -71.063991,
      radius = 10000.0)
    println("\r  [OK] Geospatial query with unknown format should fail")

    print("[  ] Geospatial query with unknown relation should fail ")
    @test_throws ErrorException geo_query(geodb, "geodd", "geoidx";
      relation = "never-heard-of-it",
      lat    = 42.357963,
      lon    = -71.063991,
      radius = 10000.0)
    println("\r  [OK] Geospatial query with unknown relation should fail")

  else
    println("  ** Skipping geospatial query tests")
    println("  ** Replicate https://education.cloudant.com/crimes and set the variable COUCH_GEO_DATABASE")
  end
end

max_stars_repo_name: xpqz/couchzilla

max_issues_repo_name: xpqz/couchzilla

max_forks_repo_name: xpqz/couchzilla


text:
push!(LOAD_PATH, "/Users/megorov/Desktop/projects/stanford/sniper/egorov/sniperRepo/new-version")

using Collada_
using UrbanMap_
using SniperPOMDP_
using SniperMDP_
using ValueIterations
using Helpers
using RayCaster_
using Simulator
using SimServers
using RolloutPolicies

map = UrbanMap(1, [0.5], [0.5], [0.42], [0.22], 10, 10)
#t   = [(2.0,8.0), (2.0,1.0), (5.0,2.0), (7.0,3.0)]
t   = [(2.0,8.0), (2.0,1.0)]
#t   = [(7.0,7.0), (2.0,14.0), (8.0,4.0), (11.0,13.0)]
mdp = SniperMDP(map, nSnipers = 1, nMonitors = 2, targets=t)

policy = RolloutPolicy(mdp, d=10, n=3000, c=0.4)

# create the simulation type
#init  = 8_888_888_803 # sniper at (3,1), all monitors at (9,9)
init = 470_872
sim  = MDPSimulation(mdp, init, 100)

# set up the server
socket = 4444
server = SniperServer(socket)

# begin simulation and send results to client
#simAndSend(server, sim, mdp, policy)

# batch sim
nSims = 10
obs = zeros(Int64, nSims)
kills = zeros(Int64, nSims)
rewards = zeros(n)

for i = 1:nSims
    println("On sim: $i")
    res = simulate(policy, mdp, 100, initState=init)
    rewards[i] = sum(res.rewards)
    obs[i] = res.nObs
    kills[i] = res.nKills
end


text:
abstract type AbstractTree{NDIMS} <: AbstractContainer end

# Type traits to obtain dimension
@inline Base.ndims(::AbstractTree{NDIMS}) where NDIMS = NDIMS

# Auxiliary methods to allow semantic queries on the tree
# Check whether cell has parent cell
has_parent(t::AbstractTree, cell_id::Int) = t.parent_ids[cell_id] > 0

# Count number of children for a given cell
n_children(t::AbstractTree, cell_id::Int) = count(x -> (x > 0), @view t.child_ids[:, cell_id])

# Check whether cell has any child cell
has_children(t::AbstractTree, cell_id::Int) = n_children(t, cell_id) > 0

# Check whether cell is leaf cell
is_leaf(t::AbstractTree, cell_id::Int) = !has_children(t, cell_id)

# Check whether cell has specific child cell
has_child(t::AbstractTree, cell_id::Int, child::Int) = t.child_ids[child, cell_id] > 0

# Check if cell has a neighbor at the same refinement level in the given direction
has_neighbor(t::AbstractTree, cell_id::Int, direction::Int) = t.neighbor_ids[direction, cell_id] > 0

# Check if cell has a coarse neighbor, i.e., with one refinement level lower
function has_coarse_neighbor(t::AbstractTree, cell_id::Int, direction::Int)
  return has_parent(t, cell_id) && has_neighbor(t, t.parent_ids[cell_id], direction)
end

# Check if cell has any neighbor (same-level or lower-level)
function has_any_neighbor(t::AbstractTree, cell_id::Int, direction::Int)
  return has_neighbor(t, cell_id, direction) || has_coarse_neighbor(t, cell_id, direction)
end

# Check if cell is own cell, i.e., belongs to this MPI rank
is_own_cell(t::AbstractTree, cell_id) = true

# Return cell length for a given level
length_at_level(t::AbstractTree, level::Int) = t.length_level_0 / 2^level

# Return cell length for a given cell
length_at_cell(t::AbstractTree, cell_id::Int) = length_at_level(t, t.levels[cell_id])

# Return minimum level of any leaf cell
minimum_level(t::AbstractTree) = minimum(t.levels[leaf_cells(t)])

# Return maximum level of any leaf cell
maximum_level(t::AbstractTree) = maximum(t.levels[leaf_cells(t)])

# Check if tree is periodic
isperiodic(t::AbstractTree) = all(t.periodicity)
isperiodic(t::AbstractTree, dimension) = t.periodicity[dimension]

# Auxiliary methods for often-required calculations
# Number of potential child cells
n_children_per_cell(::AbstractTree{NDIMS}) where NDIMS = 2^NDIMS
n_children_per_cell(dims::Integer) = 2^dims

# Number of directions
#
# Directions are indicated by numbers from 1 to 2*ndims:
# 1 -> -x
# 2 -> +x
# 3 -> -y
# 4 -> +y
# 5 -> -z
# 6 -> +z
@inline n_directions(::AbstractTree{NDIMS}) where NDIMS = 2 * NDIMS
# TODO: Taal performance, 1:n_directions(tree) vs. Base.OneTo(n_directions(tree)) vs. SOneTo(n_directions(tree))
@inline eachdirection(tree::AbstractTree) = Base.OneTo(n_directions(tree))

# For a given direction, return its opposite direction
#
# dir -> opp
#  1  ->  2
#  2  ->  1
#  3  ->  4
#  4  ->  3
#  5  ->  6
#  6  ->  5
opposite_direction(direction::Int) = direction + 1 - 2 * ((direction + 1) % 2)

# For a given child position (from 1 to 8) and dimension (from 1 to 3),
# calculate a child cell's position relative to its parent cell.
#
# Essentially calculates the following
#         dim=1 dim=2 dim=3
# child     x     y     z
#   1       -     -     -
#   2       +     -     -
#   3       -     +     -
#   4       +     +     -
#   5       -     -     +
#   6       +     -     +
#   7       -     +     +
#   8       +     +     +
child_sign(child::Int, dim::Int) = 1 - 2 * (div(child + 2^(dim - 1) - 1, 2^(dim-1)) % 2)

# For each child position (1 to 8) and a given direction (from 1 to 6), return
# neighboring child position.
adjacent_child(child::Int, direction::Int) = [2 2 3 3 5 5;
                                              1 1 4 4 6 6;
                                              4 4 1 1 7 7;
                                              3 3 2 2 8 8;
                                              6 6 7 7 1 1;
                                              5 5 8 8 2 2;
                                              8 8 5 5 3 3;
                                              7 7 6 6 4 4][child, direction]

# For each child position (1 to 8) and a given direction (from 1 to 6), return
# if neighbor is a sibling
function has_sibling(child::Int, direction::Int)
  return (child_sign(child, div(direction + 1, 2)) * (-1)^(direction - 1)) > 0
end

# Obtain leaf cells that fulfill a given criterion.
#
# The function `f` is passed the cell id of each leaf cell
# as an argument.
function filter_leaf_cells(f, t::AbstractTree)
  filtered = Vector{Int}(undef, length(t))
  count = 0
  for cell_id in 1:length(t)
    if is_leaf(t, cell_id) && f(cell_id)
      count += 1
      filtered[count] = cell_id
    end
  end

  return filtered[1:count]
end

# Return an array with the ids of all leaf cells
leaf_cells(t::AbstractTree) = filter_leaf_cells((cell_id)->true, t)

# Return an array with the ids of all leaf cells for a given rank
leaf_cells_by_rank(t::AbstractTree, rank) = leaf_cells(t)

# Return an array with the ids of all local leaf cells
local_leaf_cells(t::AbstractTree) = leaf_cells(t)

# Count the number of leaf cells.
count_leaf_cells(t::AbstractTree) = length(leaf_cells(t))

# Store cell id in each cell to use for post-AMR analysis
function reset_original_cell_ids!(t::AbstractTree)
  t.original_cell_ids[1:length(t)] .= 1:length(t)
end

# Refine entire tree by one level
refine!(t::AbstractTree) = refine!(t, leaf_cells(t))

# Refine given cells and rebalance tree.
#
# Note 1: Rebalancing is iterative, i.e., neighboring cells are refined if
#         otherwise the 2:1 rule would be violated, which can cause more
#         refinements.
# Note 2: Rebalancing currently only considers *Cartesian* neighbors, not diagonal neighbors!
function refine!(t::AbstractTree, cell_ids)
  # Reset original cell ids such that each cell knows its current id
  reset_original_cell_ids!(t)

  # Refine all requested cells
  refined = refine_unbalanced!(t, cell_ids)
  refinement_count = length(refined)

  # Iteratively rebalance the tree until it does not change anymore
  while length(refined) > 0
    refined = rebalance!(t, refined)
    refinement_count += length(refined)
  end

  # Determine list of *original* cell ids that were refined
  # Note: original_cell_ids contains the cell_id *before* refinement. At
  # refinement, the refined cell's original_cell_ids value has its sign flipped
  # to easily find it now.
  @views refined_original_cells = (
      -t.original_cell_ids[1:length(t)][t.original_cell_ids[1:length(t)] .< 0])

  # Check if count of refinement cells matches information in original_cell_ids
  @assert refinement_count == length(refined_original_cells) (
      "Mismatch in number of refined cells")

  return refined_original_cells
end

# Refine all leaf cells with coordinates in a given rectangular box
function refine_box!(t::AbstractTree{NDIMS}, coordinates_min, coordinates_max) where NDIMS
  for dim in 1:NDIMS
    @assert coordinates_min[dim] < coordinates_max[dim] "Minimum coordinates are not minimum."
  end

  # Find all leaf cells within box
  cells = filter_leaf_cells(t) do cell_id
    return (all(coordinates_min .< t.coordinates[:, cell_id]) &&
            all(coordinates_max .> t.coordinates[:, cell_id]))
  end

  # Refine cells
  refine!(t, cells)
end

# Convenience method for 1D
function refine_box!(t::AbstractTree{1}, coordinates_min::Real, coordinates_max::Real)
  return refine_box!(t, [convert(Float64, coordinates_min)], [convert(Float64, coordinates_max)])
end

# For the given cell ids, check if neighbors need to be refined to restore a rebalanced tree.
#
# Note 1: Rebalancing currently only considers *Cartesian* neighbors, not diagonal neighbors!
# Note 2: The current algorithm assumes that a previous refinement step has
#         created level differences of at most 2. That is, before the previous
#         refinement step, the tree was balanced.
function rebalance!(t::AbstractTree, refined_cell_ids)
  # Create buffer for newly refined cells
  to_refine = zeros(Int, n_directions(t) * length(refined_cell_ids))
  count = 0

  # Iterate over cell ids that have previously been refined
  for cell_id in refined_cell_ids
    # Go over all potential neighbors of child cell
    for direction in eachdirection(t)
      # Continue if refined cell has a neighbor in that direction
      if has_neighbor(t, cell_id, direction)
        continue
      end

      # Continue if refined cell has no coarse neighbor, since that would
      # mean it there is no neighbor in that direction at all (domain
      # boundary)
      if !has_coarse_neighbor(t, cell_id, direction)
        continue
      end

      # Otherwise, the coarse neighbor exists and is not refined, thus it must
      # be marked for refinement
      coarse_neighbor_id = t.neighbor_ids[direction, t.parent_ids[cell_id]]
      count += 1
      to_refine[count] = coarse_neighbor_id
    end
  end

  # Finally, refine all marked cells...
  refined = refine_unbalanced!(t, unique(to_refine[1:count]))

  # ...and return list of refined cells
  return refined
end

# Refine given cells without rebalancing tree.
#
# Note: After a call to this method the tree may be unbalanced!
function refine_unbalanced!(t::AbstractTree, cell_ids) end

# Wrap single-cell refinements such that `sort(...)` does not complain
refine_unbalanced!(t::AbstractTree, cell_id::Int) = refine_unbalanced!(t, [cell_id])

# Coarsen entire tree by one level
function coarsen!(t::AbstractTree)
  # Special case: if there is only one cell (root), there is nothing to do
  if length(t) == 1
    return Int[]
  end

  # Get list of unique parent ids for all leaf cells
  parent_ids = unique(t.parent_ids[leaf_cells(t)])
  coarsen!(t, parent_ids)
end

# Coarsen given *parent* cells (= these cells must have children who are all
# leaf cells) while retaining a balanced tree.
#
# A cell to be coarsened might cause an unbalanced tree if the neighboring cell
# was already refined. Since it is generally not desired that cells are
# coarsened without specifically asking for it, these cells will then *not* be
# coarsened.
function coarsen!(t::AbstractTree, cell_ids::AbstractArray{Int})
  # Return early if array is empty
  if length(cell_ids) == 0
    return Int[]
  end

  # Reset original cell ids such that each cell knows its current id
  reset_original_cell_ids!(t)

  # To maximize the number of cells that may be coarsened, start with the cells at the highest level
  sorted_by_level = sort(cell_ids, by = i -> t.levels[i])

  # Keep track of number of cells that were actually coarsened
  n_coarsened = 0

  # Local function to adjust cell ids after some cells have been removed
  function adjust_cell_ids!(cell_ids, coarsened_cell_id, count)
    for (id, cell_id) in enumerate(cell_ids)
      if cell_id > coarsened_cell_id
        cell_ids[id] = cell_id - count
      end
    end
  end

  # Iterate backwards over cells to coarsen
  while true
    # Retrieve next cell or quit
    if length(sorted_by_level) > 0
      coarse_cell_id = pop!(sorted_by_level)
    else
      break
    end

    # Ensure that cell has children (violation is an error)
    if !has_children(t, coarse_cell_id)
      error("cell is leaf and cannot be coarsened to: $coarse_cell_id")
    end

    # Ensure that all child cells are leaf cells (violation is an error)
    for child in 1:n_children_per_cell(t)
      if has_child(t, coarse_cell_id, child)
        if !is_leaf(t, t.child_ids[child, coarse_cell_id])
          error("cell $coarse_cell_id has child cell at position $child that is not a leaf cell")
        end
      end
    end

    # Check if coarse cell has refined neighbors that would prevent coarsening
    skip = false
    # Iterate over all children (which are to be removed)
    for child in 1:n_children_per_cell(t)
      # Continue if child does not exist
      if !has_child(t, coarse_cell_id, child)
        continue
      end
      child_id = t.child_ids[child, coarse_cell_id]

      # Go over all neighbors of child cell. If it has a neighbor that is *not*
      # a sibling and that is not a leaf cell, we cannot coarsen its parent
      # without creating an unbalanced tree.
      for direction in eachdirection(t)
        # Continue if neighbor would be a sibling
        if has_sibling(child, direction)
          continue
        end

        # Continue if child cell has no neighbor in that direction
        if !has_neighbor(t, child_id, direction)
          continue
        end
        neighbor_id = t.neighbor_ids[direction, child_id]

        if !has_children(t, neighbor_id)
          continue
        end

        # If neighbor is not a sibling, is existing, and has children, do not coarsen
        skip = true
        break
      end
    end
    # Skip if a neighboring cell prevents coarsening
    if skip
      continue
    end

    # Flip sign of cell to be coarsened to such that we can easily find it
    t.original_cell_ids[coarse_cell_id] = -t.original_cell_ids[coarse_cell_id]

    # If a coarse cell has children that are all leaf cells, they must follow
    # immediately due to depth-first ordering of the tree
    count = n_children(t, coarse_cell_id)
    @assert count == n_children_per_cell(t) "cell $coarse_cell_id does not have all child cells"
    remove_shift!(t, coarse_cell_id + 1, coarse_cell_id + count)

    # Take into account shifts in tree that alters cell ids
    adjust_cell_ids!(sorted_by_level, coarse_cell_id, count)

    # Keep track of number of coarsened cells
    n_coarsened += 1
  end

  # Determine list of *original* cell ids that were coarsened to
  # Note: original_cell_ids contains the cell_id *before* coarsening. At
  # coarsening, the coarsened parent cell's original_cell_ids value has its sign flipped
  # to easily find it now.
  @views coarsened_original_cells = (
      -t.original_cell_ids[1:length(t)][t.original_cell_ids[1:length(t)] .< 0])

  # Check if count of coarsened cells matches information in original_cell_ids
  @assert n_coarsened == length(coarsened_original_cells) (
      "Mismatch in number of coarsened cells")

  return coarsened_original_cells
end

# Wrap single-cell coarsening such that `sort(...)` does not complain
coarsen!(t::AbstractTree, cell_id::Int) = coarsen!(t::AbstractTree, [cell_id])

# Coarsen all viable parent cells with coordinates in a given rectangular box
function coarsen_box!(t::AbstractTree{NDIMS}, coordinates_min::AbstractArray{Float64},
                     coordinates_max::AbstractArray{Float64}) where NDIMS
  for dim in 1:NDIMS
    @assert coordinates_min[dim] < coordinates_max[dim] "Minimum coordinates are not minimum."
  end

  # Find all leaf cells within box
  leaves = filter_leaf_cells(t) do cell_id
    return (all(coordinates_min .< t.coordinates[:, cell_id]) &&
            all(coordinates_max .> t.coordinates[:, cell_id]))
  end

  # Get list of unique parent ids for all leaf cells
  parent_ids = unique(t.parent_ids[leaves])

  # Filter parent ids to be within box
  parents = filter(parent_ids) do cell_id
    return (all(coordinates_min .< t.coordinates[:, cell_id]) &&
            all(coordinates_max .> t.coordinates[:, cell_id]))
  end

  # Coarsen cells
  coarsen!(t, parents)
end

# Convenience method for 1D
function coarsen_box!(t::AbstractTree{1}, coordinates_min::Real, coordinates_max::Real)
  return coarsen_box!(t, [convert(Float64, coordinates_min)], [convert(Float64, coordinates_max)])
end

# Return coordinates of a child cell based on its relative position to the parent.
function child_coordinates(::AbstractTree{NDIMS}, parent_coordinates, parent_length::Number, child::Int) where NDIMS
  # Calculate length of child cells and set up data structure
  child_length = parent_length / 2
  coordinates = MVector{NDIMS, Float64}(undef)

  # For each dimension, calculate coordinate as parent coordinate + relative position x length/2
  for d in 1:NDIMS
    coordinates[d] = parent_coordinates[d] + child_sign(child, d) * child_length / 2
  end

  return coordinates
end

# Reset range of cells to values that are prone to cause errors as soon as they are used.
#
# Rationale: If an invalid cell is accidentally used, we want to know it as soon as possible.
function invalidate!(t::AbstractTree, first::Int, last::Int) end
invalidate!(t::AbstractTree, id::Int) = invalidate!(t, id, id)
invalidate!(t::AbstractTree) = invalidate!(t, 1, length(t))

# Delete connectivity with parents/children/neighbors before cells are erased
function delete_connectivity!(t::AbstractTree, first::Int, last::Int)
  @assert first > 0
  @assert first <= last
  @assert last <= t.capacity + 1

  # Iterate over all cells
  for cell_id in first:last
    # Delete connectivity from parent cell
    if has_parent(t, cell_id)
      parent_id = t.parent_ids[cell_id]
      for child in 1:n_children_per_cell(t)
        if t.child_ids[child, parent_id] == cell_id
          t.child_ids[child, parent_id] = 0
          break
        end
      end
    end

    # Delete connectivity from child cells
    for child in 1:n_children_per_cell(t)
      if has_child(t, cell_id, child)
        t.parent_ids[t._child_ids[child, cell_id]] = 0
      end
    end

    # Delete connectivity from neighboring cells
    for direction in eachdirection(t)
      if has_neighbor(t, cell_id, direction)
        t.neighbor_ids[opposite_direction(direction), t.neighbor_ids[direction, cell_id]] = 0
      end
    end
  end
end

# Move connectivity with parents/children/neighbors after cells have been moved
function move_connectivity!(t::AbstractTree, first::Int, last::Int, destination::Int)
  @assert first > 0
  @assert first <= last
  @assert last <= t.capacity + 1
  @assert destination > 0
  @assert destination <= t.capacity + 1

  # Strategy
  # 1) Loop over moved cells (at target location)
  # 2) Check if parent/children/neighbors connections are to a cell that was moved
  #    a) if cell was moved: apply offset to current cell
  #    b) if cell was not moved: go to connected cell and update connectivity there

  offset = destination - first
  has_moved(n) = (first <= n <= last)

  for source in first:last
    target = source + offset

    # Update parent
    if has_parent(t, target)
      # Get parent cell
      parent_id = t.parent_ids[target]
      if has_moved(parent_id)
        # If parent itself was moved, just update parent id accordingly
        t.parent_ids[target] += offset
      else
        # If parent was not moved, update its corresponding child id
        for child in 1:n_children_per_cell(t)
          if t.child_ids[child, parent_id] == source
            t.child_ids[child, parent_id] = target
          end
        end
      end
    end

    # Update children
    for child in 1:n_children_per_cell(t)
      if has_child(t, target, child)
        # Get child cell
        child_id = t.child_ids[child, target]
        if has_moved(child_id)
          # If child itself was moved, just update child id accordingly
          t.child_ids[child, target] += offset
        else
          # If child was not moved, update its parent id
          t.parent_ids[child_id] = target
        end
      end
    end

    # Update neighbors
    for direction in eachdirection(t)
      if has_neighbor(t, target, direction)
        # Get neighbor cell
        neighbor_id = t.neighbor_ids[direction, target]
        if has_moved(neighbor_id)
          # If neighbor itself was moved, just update neighbor id accordingly
          t.neighbor_ids[direction, target] += offset
        else
          # If neighbor was not moved, update its opposing neighbor id
          t.neighbor_ids[opposite_direction(direction), neighbor_id] = target
        end
      end
    end
  end
end

# Raw copy operation for ranges of cells.
#
# This method is used by the higher-level copy operations for AbstractContainer
function raw_copy!(target::AbstractTree, source::AbstractTree, first::Int, last::Int, destination::Int) end

# Reset data structures by recreating all internal storage containers and invalidating all elements
function reset_data_structures!(t::AbstractTree{NDIMS}) where NDIMS end


text:
"""
    metricreturntype(metric, x, y=x)

Get the expected return type of `metric(x, y)`. This is inferred statically if
possible, otherwise the function is actually called with the supplied arguments.
Always returns a concrete type.
"""
function metricreturntype(metric, x, y=x)
    # The following should return a concrete type if metric is type stable:
    T = Core.Compiler.return_type(metric, Tuple{typeof(x), typeof(y)})
    isconcretetype(T) && return T
    return typeof(metric(x, y))
end


text:
# we use warmstart to hint at a solution that has better objective value in pipe
# dimensioning than the dual bound of junction location?!

using MathOptInterface
using JuMP
using SCIP

using PipeLayout
using PipeLayout.JuncLoc
using PipeLayout.PipeDim
import PipeLayout: ploss_coeff_nice

const MOI = MathOptInterface

# instance data
pipelength = 100.0
width = height = pipelength / sqrt(2)
terminals = [Node(0,0), Node(width,height), Node(2width,0), Node(3width,height)]
demand = [400, -600, 600, -400]
bounds = fill(Bounds(40, 80), 4)
diams = map(Diameter, [1, 2], [1, 2])
inst = Instance(terminals, demand, bounds, diams, ploss_coeff_nice)

_scip = JuMP.optimizer_with_attributes(SCIP.Optimizer, "limits/gap" => 0.001)

function pipedim()
    println("### Pipe Dimensioning LP model, path topology ###")

    # topology for path
    path = Topology(terminals, [Arc(2,1), Arc(2,3), Arc(4,3)])

    # instead of calling `optimize` we do the steps manually
    solver = PipeDim.LP(_scip)
    model, π, l = PipeDim.make_model(inst, path, solver.lpsolver)
    JuMP.optimize!(model)

    println(model)
    @show JuMP.value.(l)
    @show JuMP.value.(π)
end

function juncloc()
    println("### Junction Location NLP model, FST with fixed locations ###")

    # topology data for FST
    junctions = [Node(0,-1), Node(0,-2)]
    nodes = vcat(terminals, junctions)
    fstw = Topology(nodes, [Arc(5,6), Arc(5,1), Arc(2,5), Arc(3,6), Arc(6,4)])

    # instead of calling `optimize` we do the steps manually
    solver = JuncLoc.NLP(_scip)
    model, x, y, L, l, π = JuncLoc.make_nlp(inst, fstw, solver)

    # fix values for the Steiner node positions
    # 1st steiner is at 2nd term
    @constraint(model, x[5] == inst.nodes[2].x)
    @constraint(model, y[5] == inst.nodes[2].y)
    # 2nd steiner is at 3rd term
    @constraint(model, x[6] == inst.nodes[3].x)
    @constraint(model, y[6] == inst.nodes[3].y)

    JuMP.optimize!(model)

    println(model)
    @show JuMP.value.(L)
    @show JuMP.value.(l)
    @show JuMP.value.(π)
end

pipedim()
println()
juncloc()

# demands = [400, -600, 600, -400]
#
# LP:                flow   flow^2
# 1: 1 --> 2  (a)    400     40000
# 2: 2 --> 3  (b)    200    160000
# 3: 3 --> 4  (c)    400     40000
#
#      2   4
#     / \ /
#    1   3
#
#
# NLP:                flow   flow squared
# 1: 5 --> 6 (b)      200       40000
# 2: 5 --> 1 (a)      400      160000
# 3: 2 --> 5 (-)      600      360000
# 4: 3 --> 6 (-)     -600      360000
# 5: 6 --> 4 (c)     -400      160000
#
#       2     _4
#     _ 5 -- 6
#    1       3
#


text:
@testset "NoiseApproximation" begin

using DiffEqNoiseProcess, DiffEqBase, StochasticDiffEq
using Test

using DiffEqProblemLibrary.SDEProblemLibrary: importsdeproblems; importsdeproblems()
import DiffEqProblemLibrary.SDEProblemLibrary: prob_sde_linear, prob_sde_2Dlinear

prob = prob_sde_linear
integrator = init(prob,EM(),dt=0.01)

W = NoiseApproximation(integrator)

dt = 0.1
calculate_step!(W,dt,nothing,nothing)
dWold = W.dW
@test W.curW == W[1]
@test W.curt == 0.0
accept_step!(W,dt,nothing,nothing)
@test W.curW == 0.5 + dWold
@test W.curt == dt
@test W.curW + W.dW == W[end]
@test W.curW == W[11]

W = NoiseApproximation(integrator)
for i in 1:10
  accept_step!(W,dt,nothing,nothing)
end
W.t[end] == 1.0

W = NoiseApproximation(integrator)
prob = NoiseProblem(W,(0.0,1.0))
sol = solve(prob;dt=0.1)

prob = prob_sde_2Dlinear
integrator = init(prob,EM(),dt=0.01)
W = NoiseApproximation(integrator)
prob = NoiseProblem(W,(0.0,1.0))
sol = solve(prob;dt=0.1)

end


text:
## Show the wikipedia categories graph
using SparseArrays
using LinearAlgebra
using DelimitedFiles
ei,ej = copy.(eachcol(Int.(readdlm("wiki-cats.edges"))))
xy = readdlm("wiki-cats.xy")
n = size(xy,1)
A = sparse(ei,ej,1,n,n) # we only store half the edges here
A = max.(A,A')
dropzeros!(A)
##
issymmetric(A)
##
ei,ej = findnz(triu(A,1))
using GraphRecipes
using Plots
graphplot(ei, ej, x =xy[:,1], y=xy[:,2],
  markercolor=:black, markerstrokecolor=:white,
  size=(1200,1200),dpi=400,
  markersize=2, linecolor=1, linealpha=0.2, linewidth=0.5,
  markeralpha=0.3,
  axis_buffer=0.02, background=nothing)
#savefig()
##
using MatrixNetworks
pr = pagerank(A,0.85)
##
p=sortperm(pr)
scatter(xy[p,1],xy[p,2],marker_z=-abs.(log10.(pr[p])),markerstrokewidth=0)
## Now show the entire graph
p=sortperm(pr)
ei, ej = findnz(triu(A,1))
graphplot(ei, ej, x =xy[:,1], y=xy[:,2],
  markercolor=:black, markerstrokecolor=:white,
  size=(1200,1200),dpi=300,
  markersize=0, linecolor=1, linealpha=0.2, linewidth=0.5,
  colorbar=false,
  axis_buffer=0.02, background=nothing)
scatter!(xy[p,1],xy[p,2],marker_z=-abs.(log10.(pr[p])),markerstrokewidth=0,
  markersize=-abs.(log10.(pr[p])).+6)
savefig("wiki-cats-pagerank.png")
##


text:
include("smote_exs.jl")
include("ub_smote.jl")
include("Boosting.jl")
include("simtools.jl")

# running single simulation
@time runsim(500, 20, 2.0, 0.1, seed = 111, ξ = 0.5)

# running multiple simulations
simresults = run_simulations(50, 500, 10, -2.5, 0.05, ξ = 0.6)
summarise_sims(simresults, 5)

n = 100
p = 10
# X = hcat(ones(n), randn(n, p-1))

Σ = [1.0 0.8 0.4 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
     0.8 1.0 0.6 0.2 0.0 0.0 0.0 0.0 0.0 0.0;
     0.4 0.6 1.0 0.6 0.4 0.0 0.0 0.0 0.0 0.0;
     0.0 0.2 0.6 1.0 0.8 0.0 0.0 0.0 0.0 0.0;
     0.0 0.0 0.4 0.8 1.0 0.0 0.0 0.0 0.0 0.0;
     0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0;
     0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0;
     0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0;
     0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0;
     0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0]

mvn = MvNormal(ones(p), Σ)

srand(round(Int, time()))

X = rand(mvn, n)'
cor(X)

ϵ = rand(Normal(2, 0.5), n)                                  # gives 10% minority class,
# ϵ = rand(Normal(-2.5, 0.5), n)                               # gives 5% minority class,
# ϵ = rand(Normal(4, 0.5), n)                                  # gives 15% minority class,

β = [-3, -4, 1, -5, -4, 0.0, 0.0, 0.0, 0.0, 0.0]
η = X*β + ϵ                                                    # linear predictor w/ error
pr = 1.0 ./ (1.0 + exp(-η))                                    # inv-logit

# simulate outcome variable
y = map(π -> rand(Binomial(1, π)), pr)
y[y .== 0] = -1

# y2, X2 = get_subset(y, X, 450, 0.10)
srand(111)
model, coeffs = build_adaboost_stumps(y, X, 50, 1.0);          # 50 boosting rounds, "sub-sample" 100%

# apply learned model
y_hat = apply_adaboost_stumps(model, coeffs, X)
mean(y .== y_hat)

mean(y)

# get the probability of each label
# apply_adaboost_stumps_proba(model, coeffs, [5.9, 3.0, 5.1, 1.9, 1.0, 2.1], [0, 1])

srand(111)
model, coeffs2 = build_adaboost_stumps(y, X, 50, 1.0, 0.2);

y_hat = apply_adaboost_stumps(model, coeffs2, X)
mean(y .== y_hat)
1 - mean(y)

# apply_adaboost_stumps_proba(model, coeffs, [5.9, 3.0, 5.1, 1.9, 1.0, 2.1], [0, 1])

find(y)

# run n-fold cross validation for boosted stumps, using 7 iterations and 3 folds
# accuracy = nfoldCV_stumps(y, X, 7, 3)

using DataFrames

ds = readtable("example_data.csv");
y = convert(Vector, ds[:, 1]);
X = convert(Array, ds[:, 2:end]);

using DecisionTree
using Distributions

n = 500
p = 4
# X = hcat(ones(n), randn(n, p-1))

Σ = [1.0 0.6 0.2 0.0;
          0.6 1.0 0.6 0.2;
          0.2 0.6 1.0 0.6;
          0.0 0.2 0.6 1.0]

mvn = MvNormal([1.0, -1.0, -3.0, 3.0], Σ)

X = rand(mvn, n)'
cor(X)

β = vcat([1.5, 1.1, -10.5, -100.5])
η = X*β                                                      # linear predictor
pr = 1.0 ./ (1.0 + exp(-η))                                  # inv-logit

# simulate outcome variable
y = map(π -> rand(Binomial(1, π)), pr)
mean(y)

srand(111)
model, coeffs = DecisionTree.build_adaboost_stumps(y, X, 50);

# apply learned model
y_hat = DecisionTree.apply_adaboost_stumps(model, coeffs, X)
mean(y .== y_hat)
mean(y)

# get the probability of each label
apply_adaboost_stumps_proba(model, coeffs, [5.9, 3.0, 5.1, 1.9], [0, 1])

max_stars_repo_name: paulstey/boosting

max_issues_repo_name: paulstey/boosting

max_forks_repo_name: paulstey/boosting


text:
"""
`JuLIP.FIO` : provides some basis file IO. All sub-modules and derived
packages should use these interface functions so that the file formats
can be changed later. This submodule provides

- load_dict
- save_dict
- read_dict
"""
module FIO

using JSON, ZipFile
using SparseArrays: SparseMatrixCSC

export read_dict, write_dict,
       zip_dict, unzip_dict,
       load_json, save_json,
       load_dict, save_dict

#######################################################################
#                     Conversions to and from Dict
#######################################################################

# eventually we want to be able to serialise all JuLIP types to Dict
# and back and those Dicts may only contain elementary data types
# this will then allow us to load them back via read_dict without
# having to know the type in the code

"""
`read_dict(D::Dict) -> ?`

Looks for a key `__id__` in `D` whose value is used to dynamically dispatch
the decoding to
```julia
read_dict(Val(Symbol(D["__id__"])), D)
```
That is, a user defined type must implement this `convert(::Val, ::Dict)`
utility function. The typical code would be
```julia
module MyModule
    struct MyStructA
        a::Float64
    end
    Dict(A::MyStructA) = Dict( "__id__" -> "MyModule_MyStructA",
                               "a" -> a )
    MyStructA(D::Dict) = A(D["a"])
    Base.convert(::Val{:MyModule_MyStructA})
end
```
The user is responsible for choosing an id that is sufficiently unique.

The purpose of this function is to enable the loading of more complex JSON files
and automatically converting the parsed Dict into the appropriate composite
types. It is intentional that a significant burden is put on the user code
here. This will maximise flexibiliy, e.g., by introducing version numbers,
and being able to read old version in the future.
"""
function read_dict(D::Dict)
    if !haskey(D, "__id__")
        error("JuLIP.FIO.read_dict: `D` has no key `__id__`")
    end
    if haskey(D, "__v__")
      return read_dict(Val(Symbol(D["__id__"])),
                       Val(Symbol(D["__v__"])),
                       D)
    end
    return read_dict(Val(Symbol(D["__id__"])), D)
end

read_dict(::Val{sym}, D::Dict) where {sym} =
    error("ACEbase.FIO.read_dict no implemented for symbol $(sym)")

#######################################################################
#                     JSON
#######################################################################

function load_json(fname::AbstractString)
    return JSON.parsefile(fname)
end

function save_json(fname::AbstractString, D::Dict; indent=0)
    f = open(fname, "w")
    JSON.print(f, D, indent)
    close(f)
    return nothing
end

load_dict = load_json
save_dict = save_json

function zip_dict(fname::AbstractString, D::Dict; indent=0)
   zipdir = ZipFile.Writer(fname)
   fptr = ZipFile.addfile(zipdir, "dict.json"; method=ZipFile.Deflate)
   write(fptr, JSON.json(D, indent))
   close(fptr)
   close(zipdir)
end

function unzip_dict(fname::AbstractString; verbose=false)
   # using global here is a weird workaround for a bug as suggest in
   # https://github.com/fhs/ZipFile.jl/issues/14
   global _zipdir = ZipFile.Reader(fname)
   if length(_zipdir.files) != 1
      error("Expected exactly one file in `$(fname)`")
   end
   fptr = _zipdir.files[1]
   verbose && @show fptr.name
   return JSON.parse(fptr)
   close(fptr)
   close(_zipdir)
end

## Some useful utility functions

# Datatype
write_dict(T::Type) = Dict("__id__" => "Type", "T" => string(T))
read_dict(::Val{:Type}, D::Dict) = Meta.eval(Meta.parse(D["T"]))

# Complex Vector

function write_dict(A::Vector{T}) where {T <: Number}
   D = Dict("__id__" => "ACE_VectorOfNumber",
                 "T" => write_dict(T),
              "real" => real.(A))
   if T <: Complex
      D["imag"] = imag.(A)
   end
   return D
end

function read_dict(::Val{:ACE_VectorOfNumber}, D::Dict)
   T = read_dict(D["T"])
   A = T.(D["real"])
   if T <: Complex
      A[:] .+= im .* D["imag"]
   end
   return A
end

write_dict(A::Vector{T}) where {T} =
         Dict("__id__" => "ACE_Vector",
                "vals" => write_dict.(A))

read_dict(::Val{:ACE_Vector}, D::Dict) =
         read_dict.(D["vals"])

# Matrix

write_dict(A::Matrix{T}) where {T <: Number} =
    Dict("__id__" => "JuLIP_Matrix",
         "T"      => write_dict(T),
         "nrows"  => size(A, 1),
         "ncols"  => size(A, 2),
         "data"   => A[:])

function read_dict(::Val{:JuLIP_Matrix}, D::Dict)
   T = read_dict(Val(:Type), D["T"])
   return reshape(T.(D["data"]), D["nrows"], D["ncols"])
end

# SparseMatrixCSC

write_dict(A::SparseMatrixCSC{TF, TI}) where {TF, TI} =
   Dict("__id__" => "SparseMatrixCSC",
        "TI" => write_dict(TI),
        "colptr" => A.colptr,
        "rowval" => A.rowval,
        "nzval" => write_dict(A.nzval),
        "m" => A.m,
        "n" => A.n )

function read_dict(::Val{:SparseMatrixCSC}, D::Dict)
   TI = read_dict(D["TI"])
   return SparseMatrixCSC(Int(D["m"]), Int(D["n"]),
                           TI.(D["colptr"]), TI.(D["rowval"]),
                           read_dict(D["nzval"]))
end

##
end


text:
using Test, DiscreteExteriorCalculus
const DEC = DiscreteExteriorCalculus
using LinearAlgebra: I, norm
using SparseArrays: sparse, spzeros

@testset "circumcenter_hodge and exterior_derivative" begin
    # setup
    begin
        m = Metric(2)
        n = 10
        _, tcomp = DEC.triangulated_lattice(n * [1,0], n * [.5, .5 * sqrt(3)], n, n)

        comp = tcomp.complex
        @test simplicial(comp)
        @test one_sided(m, comp)
        @test pairwise_delaunay(m, comp)
        @test well_centered(m, comp)
        mesh = Mesh(tcomp, circumcenter(m))
        N, K = 2, 3
        @test typeof(mesh) <: Mesh{N, K}
    end
    # test hodge star
    begin
        ★s = [DEC.circumcenter_hodge(m, mesh, k, true)
            for k in 1:length(mesh.primal.complex.cells)]
        @test map(t -> size(t, 1), ★s) == map(length, mesh.primal.complex.cells)
        @test map(t -> size(t, 2), ★s) == map(length, mesh.primal.complex.cells)
        dual_★s = [DEC.circumcenter_hodge(m, mesh, k, false)
            for k in 1:length(mesh.dual.complex.cells)]
        @test map(t -> size(t, 1), dual_★s) == map(length, mesh.dual.complex.cells)
        @test map(t -> size(t, 2), dual_★s) == map(length, mesh.dual.complex.cells)
        for (m1, m2, s) in zip(★s, reverse(dual_★s), [1,-1,1])
            @test count(!iszero, m1 * m2 - I*s) == 0 # ★★ ∝ I
        end
    end
    # test exterior derivative
    begin
        ds = [DEC.exterior_derivative(mesh.primal.complex, k)
            for k in 1:length(mesh.primal.complex.cells)]
        @test map(t -> size(t, 2), ds) == map(length, mesh.primal.complex.cells)
        @test map(t -> size(t, 1), ds) == [map(length, mesh.primal.complex.cells)[2:end]..., 0]
        @test all([count(!iszero, ds[i+1] * ds[i]) == 0 for i in 1:(length(ds)-1)]) # d² = 0
        dual_ds = [DEC.exterior_derivative(mesh.dual.complex, k)
            for k in 1:length(mesh.dual.complex.cells)]
        @test map(t -> size(t, 2), dual_ds) == map(length, mesh.dual.complex.cells)
        @test map(t -> size(t, 1), dual_ds) == [map(length, mesh.dual.complex.cells)[2:end]..., 0]
        @test all([count(!iszero, dual_ds[i+1] * dual_ds[i]) == 0 for i in 1:(length(dual_ds)-1)]) # d² = 0
        # use children to write an alternative definition of exterior_derivative
        function _exterior_derivative(comp::CellComplex{N, K}, k::Int) where {N, K}
            @assert 1 <= k <= K
            if k == K
                return spzeros(Int, 0, length(comp.cells[k]))
            else
                row_inds, col_inds, vals = Int[], Int[], Int[]
                for (row_ind, cell) in enumerate(comp.cells[k+1])
                    for c in cell.children
                        o = c.parents[cell]
                        col_ind = findfirst(isequal(c), comp.cells[k])
                        push!(row_inds, row_ind); push!(col_inds, col_ind); push!(vals, 2 * o - 1)
                    end
                end
                num_rows = length(comp.cells[k+1])
                num_cols = length(comp.cells[k])
                return sparse(row_inds, col_inds, vals, num_rows, num_cols)
            end
        end
        @test [_exterior_derivative(mesh.primal.complex, k)
            for k in 1:length(mesh.primal.complex.cells)] == ds
        @test [_exterior_derivative(mesh.dual.complex, k)
            for k in 1:length(mesh.dual.complex.cells)] == dual_ds
    end
    # test differential_operator
    begin
        differential_operator_sequence(m, mesh, "★d★d", 1, true) ==
            [dual_★s[3], dual_ds[2],★s[2],ds[1]]
        @test differential_operator(m, mesh, "★d★d", 1, true) ==
            dual_★s[3] * dual_ds[2] * ★s[2] * ds[1]
        v = ones(length(mesh.primal.complex.cells[1]))
        @test norm(differential_operator(m, mesh, "★d★d", 1, true, v) -
            dual_★s[3] * dual_ds[2] * ★s[2] * ds[1] * v) < 1e-14
        for primal in [true, false]
            for k in 1:2
                @test count(!iszero, differential_operator(m, mesh, "dd", k, primal)) == 0 # d² = 0
            end
            for (k, s) in zip(1:3, [1,-1,1])
                ★★ = differential_operator(m, mesh, "★★", k, primal)
                @test count(!iszero, ★★ - I*s) == 0 # ★★ ∝ I
                ★★ = differential_operator(m, mesh, "★", K-k+1, !primal) *
                    differential_operator(m, mesh, "★", k, primal)
                @test count(!iszero, ★★ - I*s) == 0 # ★★ ∝ I
            end
        end
    end
end


text:
# This file is part of Kpax3. License is MIT.

function test_basic_functions_shuffleunits()
  x = [1; 2; 3; 4; 5; 6]
  y = [1; 2; 3; 4; 5; 6]

  S = 4

  N = 1000000

  p1 = 1 / factorial(length(x))
  p2 = 1 / factorial(S)

  v1 = 0
  v2 = 0

  for i in 1:N
    Kpax3.shuffleunits!(x)

    if x == [3; 1; 2; 4; 5; 6]
      v1 += 1
    end

    Kpax3.shuffleunits!(y, S)

    if y == [2; 1; 3; 4; 5; 6]
      v2 += 1
    end
  end

  v1 /= N
  v2 /= N

  @test isapprox(v1, p1, atol=0.005)
  @test isapprox(v2, p2, atol=0.005)

  nothing
end

test_basic_functions_shuffleunits()


text:
import Pkg
using Base64

const default_binder_url = "https://mybinder.org/v2/gh/fonsp/pluto-on-binder/v$(string(PLUTO_VERSION))"

const cdn_version_override = nothing
# const cdn_version_override = "2a48ae2"

if cdn_version_override !== nothing
    @warn "Reminder to fonsi: Using a development version of Pluto for CDN assets. The binder button might not work. You should not see this on a released version of Pluto." cdn_version_override
end

"""
See [PlutoSliderServer.jl](https://github.com/JuliaPluto/PlutoSliderServer.jl) if you are interested in exporting notebooks programatically.
"""
function generate_html(;
        version=nothing, pluto_cdn_root=nothing,
        notebookfile_js="undefined", statefile_js="undefined",
        slider_server_url_js="undefined", binder_url_js=repr(default_binder_url),
        disable_ui=true
    )::String

    original = read(project_relative_path("frontend", "editor.html"), String)

    cdn_root = if pluto_cdn_root === nothing
        if version === nothing
            version = PLUTO_VERSION
        end
        "https://cdn.jsdelivr.net/gh/fonsp/Pluto.jl@$(something(cdn_version_override, string(PLUTO_VERSION)))/frontend/"
    else
        pluto_cdn_root
    end

    @debug "Using CDN for Pluto assets:" cdn_root

    cdnified = replace(
    replace(original,
        "href=\"./" => "href=\"$(cdn_root)"),
        "src=\"./" => "src=\"$(cdn_root)")

    result = replace(cdnified,
        "<!-- [automatically generated launch parameters can be inserted here] -->" =>
        """
        <script data-pluto-file="launch-parameters">
        window.pluto_notebookfile = $(notebookfile_js)
        window.pluto_disable_ui = $(disable_ui ? "true" : "false")
        window.pluto_slider_server_url = $(slider_server_url_js)
        window.pluto_binder_url = $(binder_url_js)
        window.pluto_statefile = $(statefile_js)
        </script>
        <!-- [automatically generated launch parameters can be inserted here] -->
        """
    )

    return result
end

function generate_html(notebook; kwargs...)::String
    state = notebook_to_js(notebook)

    notebookfile_js = let
        notebookfile64 = base64encode() do io
            save_notebook(io, notebook)
        end

        "\"data:text/julia;charset=utf-8;base64,$(notebookfile64)\""
    end

    statefile_js = let
        statefile64 = base64encode() do io
            pack(io, state)
        end

        "\"data:;base64,$(statefile64)\""
    end

    generate_html(; statefile_js=statefile_js, notebookfile_js=notebookfile_js, kwargs...)
end


text:
#---------------------------------------------------------------------#
#This function computes the Grid and Elements in 1D.
#Written by F.X. Giraldo on April 19, 2019
#           Department of Applied Mathematics
#           Naval Postgraduate School
#           Monterey; CA 93943-5216
#---------------------------------------------------------------------#
function compute_norms(qn,qe,npoin,DFloat)

    norms=zero(DFloat,3)

    l1_top, l1_bot = 0, 0
    l2_top, l2_bot = 0, 0
    l8_top, l8_bot =-1000, -1000
    for i=1:npoin
        l1_top=l1_top + abs(qn[i]-qe[i])
        l1_bot=l1_bot + abs(qe[i])
        l2_top=l2_top + (qn[i]-qe[i])^2
        l2_bot=l2_bot + qe[i]^2
        l8_top=max(l8_top, abs(qn[i]-qe[i]))
        l8_bot=max(l8_bot, abs(qe[i]))
    end
    norms[1]=( l1_top/l1_bot )
    norms[2]=max( sqrt( l2_top/l2_bot ), machine_zero)
    norms[3]=( l8_top/l8_bot )

    return (norms)
end


text:
abstract type FrequencyEstimator <: InverterComponent end

@def freq_estimation_ports begin
    state_input = [:vd_cap,:vq_cap, :δθ_vsm]
    #TODO: Move PLL to PCC, i.e. move v_cap (D'Arco v_o), to inner inputs
    inner_input = [Vdo_var,Vqo_var,δdqRI_var,ω_freq_estimator_var]
end

"""
Parameters of a Phase-Locked Loop (PLL) for VSM
# Conmutable structor
```julia
PLL(ω_lp, kp_pll, ki_pll)
```

# Arguments
*  `ω_lp`::Float64 : PLL low-pass filter frequency (rad/sec)
*  `kp_pll`::Float64 : PLL proportional gain
*  `ki_pll`::Float64 : PLL integral gain
"""
mutable struct PLL <: FrequencyEstimator
    ω_lp::Float64
    kp_pll::Float64
    ki_pll::Float64
    n_states::Int64
    states::Vector{Symbol}
    ports::Ports

        function PLL(ω_lp::Float64,
                     kp_pll::Float64,
                     ki_pll::Float64)

            @freq_estimation_ports

            new(ω_lp,
                kp_pll,
                ki_pll,
                4,
                [:vpll_d, :vpll_q, :ε_pll, :δθ_pll],
                Ports(state_input, inner_input))

        end

end


text:
using ConvDiffMIPDECO
using Test
using jInv.Mesh

domain = [0. 3. 0 1. 0 2.]
n      = 3*[7 9 12] .- 1
M      = getRegularMesh(domain,n)
Mass, Mass_const, SM = getFEMMatrices3D(M)

e = ones(prod(M.n.+1))
@test abs(prod((domain[2:2:end]-domain[1:2:end])) - dot(e,Mass*e))/dot(e,Mass*e) < 1e-2

f = getFEMsource3D(M)
v = Mass_const*f
u = SM\v

max_stars_repo_name: JuliaInv/ConvDiffMIPDECO

max_issues_repo_name: JuliaInv/ConvDiffMIPDECO

max_forks_repo_name: JuliaInv/ConvDiffMIPDECO