text:
const attributes_setters = Dict(
    :output_style => sass_option_set_output_style,
    :source_comments => sass_option_set_source_comments,
    :source_map_file => sass_option_set_source_map_file,
    :omit_source_map_url => sass_option_set_omit_source_map_url,
    :source_map_embed => sass_option_set_source_map_embed,
    :source_map_contents => sass_option_set_source_map_contents,
    :source_map_file_urls => sass_option_set_source_map_file_urls,
    :source_map_root => sass_option_set_source_map_root,
    :is_indented_syntax_src => sass_option_set_is_indented_syntax_src,
    :include_path => sass_option_set_include_path,
    :plugin_path => sass_option_set_plugin_path,
    :indent => sass_option_set_indent,
    :linefeed => sass_option_set_linefeed,
    :input_path => sass_option_set_input_path,
    :output_path => sass_option_set_output_path,
    :precision => sass_option_set_precision,
)

"""
`compile_file(filename; kwargs...)`

Compile `filename` from sass or scss to a css string. Possible options, given by keyword
arguments, are:

- `output_style`: output style for the generated css code. See `Sass.Style` for options. For example `output_style  = Sass.nested`
- `source_comments`: a boolean to specify whether to insert inline source comments
- `source_map_file`: path to source map file, enables the source map generating used to create sourceMappingUrl
- `omit_source_map_url`: disable sourceMappingUrl in css output
- `source_map_embed`: embed sourceMappingUrl as data uri
- `source_map_contents`: embed include contents in maps
- `source_map_file_urls`: create file urls for sources
- `source_map_root`: pass-through as sourceRoot property
- `is_indented_syntax_src`: treat source_string as sass (as opposed to scss)
- `include_paths` (`AbstractString` or `AbstractArray{<:AbstractString}`)
- `plugin_paths` (`AbstractString` or `AbstractArray{<:AbstractString}`)
- `indent`: string to be used for indentation
- `linefeed`: string to be used to for line feeds
- `input_path`: the input path is used for source map generating. It can be used to define something with string compilation or to overload the input file path. It is set to `stdin` for data contexts and to the input file on file contexts.
- `output_path`: the output path is used for source map generating. LibSass will not write to this file, it is just used to create information in source-maps etc.
- `precision`: precision for outputting fractional numbers

## Examples

```julia
julia> filename = joinpath(Sass.examplefolder, "test.sass");

julia> Sass.compile_file(filename; output_style = Sass.compressed)
"body{font:100% Helvetica,sans-serif;color:#333}\n"
```
"""
function compile_file(filename; input_path = filename, source_map_file = nothing, kwargs...)
    ctx = sass_make_file_context(filename)
    ctx_out = sass_file_context_get_context(ctx)
    options = sass_context_get_options(ctx)
    for (key, val) in kwargs
        setter = get(attributes_setters, key, nothing)
        setter === nothing || setter(options, val)
    end
    sass_option_set_input_path(options, input_path)
    source_map_file === nothing || sass_option_set_source_map_file(options, source_map_file)

    sass_file_context_set_options(ctx, options)

    sass_compile_file_context(ctx)

    status = sass_context_get_error_status(ctx_out)

    status == 0 || error(sass_context_get_error_text(ctx_out))

    css =  sass_context_get_output_string(ctx_out)
    ret = source_map_file === nothing ? css : (css, sass_context_get_source_map_string(ctx_out))

    sass_delete_file_context(ctx)
    ret
end

"""
`compile_file(filename, dest; kwargs...)`

Same as `compile_file(filename; kwargs)` but writes the resulting string in file `dest`.
"""
function compile_file(filename, dest; output_path = dest, source_map_file = nothing, kwargs...)
    if source_map_file === nothing
        css = compile_file(filename; output_path = output_path, kwargs...)
        open(dest, "w") do io
             write(io, css)
         end
    else
        css, src_map = compile_file(filename; output_path = output_path, source_map_file = source_map_file, kwargs...)
        open(dest, "w") do io
            write(io, css)
        end
        open(source_map_file, "w") do io
            write(io, src_map)
        end
    end
 end


text:
using BenchmarkTools
using Piecewise

function piecewise_hardcoded(x)
    if x < -1
        return 0.0
    elseif x < 0
        return 1.0 - x^2
    elseif x < 1
        return x^2 -1.0
    else
        return 0.0
    end
end

polynomials = [
    StaticPolynomial(0.0, 0.0, 0.0),
    p"1.0-x^2",
    p"x^2 - 1.0",
    StaticPolynomial(0.0, 0.0, 0.0),
]

breakpoints = [-1.0, 0.0, 1.0]

mypiecewise = PiecewisePolynomial{3}(polynomials, breakpoints)

println("Hardcoded piecewise polynomial")
@btime piecewise_hardcoded(x) setup=(x=4*rand()-2)
println("Generated piecewise polynomial")
@btime $mypiecewise(x) setup=(x=4*rand()-2)

myderivative = differentiate(mypiecewise)

println("Generated derivative")
@btime $myderivative(x) setup=(x=4*rand()-2)

functions = [
    x -> 0.0,
    x -> 1.0 - x^2,
    x -> x^2 - 1.0,
    x -> 0.0
]

mypiecewisefunction = Piecewise.OrderedPiecewiseFunction{3}(functions, breakpoints)

println("Generated piecewise function")
@btime $mypiecewisefunction(x) setup=(x=4*rand()-2)

op_macro = Piecewise.@ordered_piecewise begin
    -1.0 => x -> 0.0
    0.0 => x -> 1.0 - x^2
    1.0 => x -> x^2- 1.0
    _ => x -> 0.0
end

pp_macro = Piecewise.@piecewise_polynomial begin
    -1.0 => zero(StaticPolynomial{Float64, 3})
    0.0 => p"1.0 - x^2"
    1.0 => p"x^2 - 1.0"
    _ => zero(StaticPolynomial{Float64, 3})
end

pp_deriv = differentiate(pp_macro)
pp_integ = integrate(pp_macro)

println("Macro ordered piecewise")
@btime $op_macro(x) setup=(x=4*rand()-2)

println("Macro piecewise polynomial")
@btime $pp_macro(x) setup=(x=4*rand()-2)

println("Derivative of macro piecewise polynomial")
@btime $pp_deriv(x) setup=(x=4*rand()-2)

println("Integral of macro piecewise polynomial")
@btime $pp_integ(x) setup=(x=4*rand()-2)

ext: jl


text:
#=

Proposal functions for joint
Biogeographic competition model

Ignacio Quintero Mächler

t(-_-t)

May 16 2017

=#

#=
-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
`Y` IID proposal functions
-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
=#

"""
    upnode!(λ1     ::Float64,
            λ0     ::Float64,
            triad  ::Array{Int64,1},
            Y      ::Array{Int64,3},
            stemevs::Array{Array{Float64,1},1},
            bridx_a::Array{Array{UnitRange{Int64},1},1},
            brδt   ::Vector{Vector{Float64}},
            brl    ::Vector{Float64},
            brs    ::Array{Int64,3},
            narea  ::Int64,
            nedge  ::Int64)

Update node and incident branches using discrete
Data Augmentation for all areas using a non-competitive
mutual-independence Markov model.
"""
function upnode!(λ1     ::Float64,
                 λ0     ::Float64,
                 triad  ::Array{Int64,1},
                 Y      ::Array{Int64,3},
                 stemevs::Array{Array{Float64,1},1},
                 bridx_a::Array{Array{UnitRange{Int64},1},1},
                 brδt   ::Vector{Vector{Float64}},
                 brl    ::Vector{Float64},
                 brs    ::Array{Int64,3},
                 narea  ::Int64,
                 nedge  ::Int64)

  @inbounds begin

    # define branch triad
    pr, d1, d2 = triad

    # sample
    samplenode!(λ1, λ0, pr, d1, d2, brs, brl, narea)

    # save extinct
    ntries = 1
    while iszero(sum(view(brs,pr,2,:)))
      samplenode!(λ1, λ0, pr, d1, d2, brs, brl, narea)

      ntries += 1
      if ntries == 500
        return false
      end
    end

    # set new node in Y
    @simd for k in Base.OneTo(narea)
      Y[bridx_a[k][d1][1]] = Y[bridx_a[k][d2][1]] = brs[pr,2,k]
    end

    # sample a consistent history
    createhists!(λ1, λ0, Y, pr, d1, d2, brs, brδt, bridx_a, narea, nedge,
                 stemevs, brl[nedge])

    ntries = 1
    while ifextY(Y, stemevs, triad, brs, brl[nedge], narea, bridx_a, nedge)
      createhists!(λ1, λ0, Y, pr, d1, d2, brs, brδt, bridx_a, narea, nedge,
                   stemevs, brl[nedge])

      ntries += 1
      if ntries == 500
        return false
      end
    end

  end

  return true
end

"""
    samplenode!(λ1   ::Float64,
                λ0   ::Float64,
                pr   ::Int64,
                d1   ::Int64,
                d2   ::Int64,
                brs  ::Array{Int64,3},
                brl  ::Array{Float64,1},
                narea::Int64)

Sample one internal node according to
mutual-independence model transition probabilities.
"""
function samplenode!(λ1   ::Float64,
                     λ0   ::Float64,
                     pr   ::Int64,
                     d1   ::Int64,
                     d2   ::Int64,
                     brs  ::Array{Int64,3},
                     brl  ::Array{Float64,1},
                     narea::Int64)
  @inbounds begin

    # estimate transition probabilities
    pr0_1, pr0_2 = Ptrfast_start(λ1, λ0, brl[pr], Val{0})
    pr1_1, pr1_2 = Ptrfast_start(λ1, λ0, brl[pr], Val{1})
    d10_1, d10_2 = Ptrfast_end(  λ1, λ0, brl[d1], Val{0})
    d11_1, d11_2 = Ptrfast_end(  λ1, λ0, brl[d1], Val{1})
    d20_1, d20_2 = Ptrfast_end(  λ1, λ0, brl[d2], Val{0})
    d21_1, d21_2 = Ptrfast_end(  λ1, λ0, brl[d2], Val{1})

    for k = Base.OneTo(narea)

      if iszero(brs[pr,1,k])
        ppr_1, ppr_2 = pr0_1, pr0_2
      else
        ppr_1, ppr_2 = pr1_1, pr1_2
      end

      if iszero(brs[d1,2,k])
        pd1_1, pd1_2 = d10_1, d10_2
      else
        pd1_1, pd1_2 = d11_1, d11_2
      end

      if iszero(brs[d2,2,k])
        pd2_1, pd2_2 = d20_1, d20_2
      else
        pd2_1, pd2_2 = d21_1, d21_2
      end

      tp = normlize(*(ppr_1, pd1_1, pd2_1),
                    *(ppr_2, pd1_2, pd2_2))::Float64

      # sample the node's character
      brs[pr,2,k] = brs[d1,1,k] = brs[d2,1,k] = coinsamp(tp)::Int64
    end
  end

  return nothing
end

"""
    createhists!(λ::Array{Float64,1}, Y::Array{Int64,3}, pr::Int64, d1::Int64, d2::Int64, brs::Array{Int64,3}, brδt::Array{Array{Float64,1},1}, bridx_a::Array{Array{Array{Int64,1},1},1}, narea::Int64)

Create bit histories for all areas for the branch trio.
"""
function createhists!(λ1     ::Float64,
                      λ0     ::Float64,
                      Y      ::Array{Int64,3},
                      pr     ::Int64,
                      d1     ::Int64,
                      d2     ::Int64,
                      brs    ::Array{Int64,3},
                      brδt   ::Array{Array{Float64,1},1},
                      bridx_a::Array{Array{UnitRange{Int64},1},1},
                      narea  ::Int64,
                      nedge  ::Int64,
                      stemevs::Array{Array{Float64,1},1},
                      stbrl  ::Float64)

  @inbounds begin

    if pr == nedge
      # if stem branch do continuous DA
      mult_rejsam!(stemevs, brs, λ1, λ0, stbrl, narea, nedge)
      for j = Base.OneTo(narea), idx = (d1,d2)
        bit_rejsam!(Y, bridx_a[j][idx], brs[idx,2,j],
                    λ1, λ0, brδt[idx])
      end
    else
      for j = Base.OneTo(narea), idx = (pr,d1,d2)
        bit_rejsam!(Y, bridx_a[j][idx], brs[idx,2,j],
                    λ1, λ0, brδt[idx])
      end
    end
  end

  return nothing
end

"""
    mult_rejsam!(evs  ::Array{Array{Float64,1},1},
                 brs  ::Array{Int64,3},
                 λ1   ::Float64,
                 λ0   ::Float64,
                 t    ::Float64,
                 narea::Int64,
                 nedge::Int64)

  Multi-area branch rejection independent model sampling.
"""
function mult_rejsam!(evs  ::Array{Array{Float64,1},1},
                      brs  ::Array{Int64,3},
                      λ1   ::Float64,
                      λ0   ::Float64,
                      t    ::Float64,
                      narea::Int64,
                      nedge::Int64)

  @simd for k = Base.OneTo(narea)
    rejsam!(evs[k], brs[nedge,1,k], brs[nedge,2,k], λ1, λ0, t)
  end

  return nothing
end

"""
    ifextY(Y      ::Array{Int64,3},
           triad  ::Array{Int64,1},
           narea  ::Int64,
           bridx_a::Array{Array{UnitRange{Int64},1},1})

Return `true` if at some point the species
goes extinct and/or more than one change is
observed after some **δt**, otherwise returns `false`.
"""
function ifextY(Y      ::Array{Int64,3},
                stemevs::Array{Array{Float64,1},1},
                triad  ::Array{Int64,1},
                brs    ::Array{Int64,3},
                stbrl  ::Float64,
                narea  ::Int64,
                bridx_a::Array{Array{UnitRange{Int64},1},1},
                nedge  ::Int64)

  @inbounds begin

    if triad[1] == nedge
      ifext_cont(stemevs, brs, stbrl, narea, nedge) && return true::Bool

      for k in (triad[2],triad[3])
        ifext_disc(Y, k, narea, bridx_a) && return true::Bool
      end
    else

      for k in triad
        ifext_disc(Y, k, narea, bridx_a) && return true::Bool
      end
    end
  end

  return false::Bool
end

"""
    ifext_cont(t_hist::Array{Array{Float64,1},1},
               brs   ::Array{Int64,3},
               t     ::Float64,
               narea ::Int64,
               nedge ::Int64)

Return true if lineage goes extinct.
"""
function ifext_cont(t_hist::Array{Array{Float64,1},1},
                    brs   ::Array{Int64,3},
                    t     ::Float64,
                    narea ::Int64,
                    nedge ::Int64)

  @inbounds begin

    ioc = 0
    # initial occupancy time
    for k in Base.OneTo(narea)
      if brs[nedge,1,k] == 1
        ioc = k
        break
      end
    end

    ioct = t_hist[ioc][1]

    ntries = 0
    while !isapprox(ioct, t, atol = 1.0e-12)

      if ioc == narea
        ioc = 1
      else
        ioc += 1
      end

      tc = 0.0
      cs = brs[nedge,1,ioc]
      for ts in t_hist[ioc]
        tc += ts
        if ioct < tc
          if cs == 1
            ioct  = tc
            ntries = 0
            break
          else
            ntries += 1
            if ntries > narea
              return true
            end
            break
          end
        end
        cs = 1 - cs
      end

    end
  end

  return false::Bool
end

"""
    ifext_disc(Y      ::Array{Int64,3},
               br     ::Int64,
               narea  ::Int64,
               bridx_a::Array{Array{UnitRange{Int64},1},1})

Return `true` if at some point the species
goes extinct and/or more than one change is
observed after some **δt**, otherwise returns `false`.
This specific method is for single branch updates.
"""
function ifext_disc(Y      ::Array{Int64,3},
                    br     ::Int64,
                    narea  ::Int64,
                    bridx_a::Array{Array{UnitRange{Int64},1},1})

  @inbounds begin

    for i = Base.OneTo(length(bridx_a[1][br]::UnitRange{Int64})-1)
      s_e::Int64 = 0            # count current areas
      s_c::Int64 = 0            # count area changes

      for k = Base.OneTo(narea)
        s_e += Y[bridx_a[k][br][i]]::Int64
        if Y[bridx_a[k][br][i]]::Int64 != Y[bridx_a[k][br][i+1]]::Int64
          s_c += 1
        end
      end

      if s_e == 0 || s_c > 1
        return true::Bool
      end
    end

  end

  return false::Bool
end

#=
-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
Y stem node proposal function (continuous DA)
-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
=#

"""
    upstemnode!(λ1   ::Float64,
                λ0   ::Float64,
                nedge::Int64,
                stemevs::Array{Array{Float64,1},1},
                brs  ::Array{Int64,3},
                brl  ::Array{Float64,1},
                narea::Int64)

Update stem node.
"""
function upstemnode!(λ1     ::Float64,
                     λ0     ::Float64,
                     nedge  ::Int64,
                     stemevs::Array{Array{Float64,1},1},
                     brs    ::Array{Int64,3},
                     stbrl  ::Float64,
                     narea  ::Int64)

  @inbounds begin

    # sample
    samplestem!(λ1, λ0, nedge, brs, stbrl, narea)

    # save extinct
    ntries = 1
    while sum(view(brs,nedge,1,:)) < 1
      samplestem!(λ1, λ0, nedge, brs, stbrl, narea)
      ntries += 1
      if ntries == 500
        return false::Bool
      end
    end

    # sample a congruent history
    mult_rejsam!(stemevs, brs, λ1, λ0, stbrl, narea, nedge)

    ntries = 1
    # check if extinct

    while ifext_cont(stemevs, brs, stbrl, narea, nedge)
      mult_rejsam!(stemevs, brs, λ1, λ0, stbrl, narea, nedge)

      ntries += 1
      if ntries == 500
        return false::Bool
      end
    end
  end

  return true::Bool
end

"""
    samplestem!(λ1   ::Float64,
                λ0   ::Float64,
                nedge::Int64,
                brs  ::Array{Int64,3},
                brl  ::Array{Float64,1},
                narea::Int64)

Sample stem node.
"""
function samplestem!(λ1   ::Float64,
                     λ0   ::Float64,
                     nedge::Int64,
                     brs  ::Array{Int64,3},
                     stbrl::Float64,
                     narea::Int64)

 @inbounds begin

    # estimate transition probabilities
    p0 = normlize(Ptrfast_end(λ1, λ0, stbrl, Val{0}))::Float64
    p1 = normlize(Ptrfast_end(λ1, λ0, stbrl, Val{1}))::Float64

    for k = Base.OneTo(narea)
      # sample the node's character
      if iszero(brs[nedge,2,k])
        brs[nedge,1,k] = coinsamp(p0)::Int64
      else
        brs[nedge,1,k] = coinsamp(p1)::Int64
      end
    end
  end

  return nothing
end

#=
-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
Y branch proposal functions
-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
=#

"""
upbranchY!(λ1     ::Float64,
           λ0     ::Float64,
           ω1     ::Float64,
           ω0     ::Float64,
           avg_Δx ::Array{Float64,2},
           br     ::Int64,
           Y      ::Array{Int64,3},
           stemevc::Array{Array{Float64,1},1},
           wareas ::Array{Int64,1},
           bridx_a::Array{Array{UnitRange{Int64},1},1},
           brδt   ::Vector{Vector{Float64}},
           brl    ::Vector{Float64},
           brs    ::Array{Int64,3},
           narea  ::Int64,
           nedge  ::Int64)

Update one branch using discrete Data Augmentation
for all areas with independent
proposals taking into account `Δx` and `ω1` & `ω0`.
"""
function upbranchY!(λ1     ::Float64,
                    λ0     ::Float64,
                    br     ::Int64,
                    Y      ::Array{Int64,3},
                    stemevs::Array{Array{Float64,1},1},
                    bridx_a::Array{Array{UnitRange{Int64},1},1},
                    brδt   ::Vector{Vector{Float64}},
                    stbrl  ::Float64,
                    brs    ::Array{Int64,3},
                    narea  ::Int64,
                    nedge  ::Int64)

  ntries = 1

  # if stem branch
  if br == nedge
    mult_rejsam!(stemevs, brs, λ1, λ0, stbrl, narea, nedge)

    # check if extinct
    while ifext_cont(stemevs, brs, stbrl, narea, nedge)
      mult_rejsam!(stemevs, brs, λ1, λ0, stbrl, narea, nedge)

      ntries += 1
      if ntries == 500
        return false
      end
    end

  else
    createhists!(λ1, λ0, Y, br, brs, brδt, bridx_a, narea)

    # check if extinct
    while ifext_disc(Y, br, narea, bridx_a)
      createhists!(λ1, λ0, Y, br, brs, brδt, bridx_a, narea)

      ntries += 1
      if ntries == 500
        return false
      end
    end
  end

  return true
end

"""
    createhists!(λ1     ::Float64,
                 λ0     ::Float64,
                 ω1     ::Float64,
                 ω0     ::Float64,
                 avg_Δx ::Array{Float64,2},
                 Y      ::Array{Int64,3},
                 br     ::Int64,
                 brs    ::Array{Int64,3},
                 brδt   ::Array{Array{Float64,1},1},
                 bridx_a::Array{Array{UnitRange{Int64},1},1},
                 narea  ::Int64)

Create bit histories for all areas for one single branch
taking into account `Δx` and `ω1` & `ω0` for all areas.
"""
function createhists!(λ1     ::Float64,
                      λ0     ::Float64,
                      Y      ::Array{Int64,3},
                      br     ::Int64,
                      brs    ::Array{Int64,3},
                      brδt   ::Array{Array{Float64,1},1},
                      bridx_a::Array{Array{UnitRange{Int64},1},1},
                      narea  ::Int64)
  @inbounds begin
    for j = Base.OneTo(narea)
      bit_rejsam!(Y, bridx_a[j][br], brs[br,2,j], λ1, λ0, brδt[br])
    end
  end

  return nothing
end

#=
-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
X proposal functions
-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
=#

"""
    uptrioX!(pr   ::Int64,
             d1   ::Int64,
             d2   ::Int64,
             X    ::Array{Float64,2},
             bridx::Array{UnitRange{Int64},1},
             brδt ::Array{Array{Float64,1},1},
             σ²c  ::Float64)

Update the node and adjoining branhces of `trio` using Brownian bridges.
"""
function uptrioX!(pr   ::Int64,
                  d1   ::Int64,
                  d2   ::Int64,
                  X    ::Array{Float64,2},
                  bridx::Array{UnitRange{Int64},1},
                  brδt ::Array{Array{Float64,1},1},
                  brl  ::Array{Float64,1},
                  σ²ϕ  ::Float64,
                  nedge::Int64)
  @inbounds begin

    # if not root
    if pr != nedge

      ipr = bridx[pr]
      id1 = bridx[d1]
      id2 = bridx[d2]

      # update node
      X[id1[1]] =
      X[id2[1]] = trioupd(X[ipr[1]],
                          X[id1[end]],
                          X[id2[end]],
                          brl[pr], brl[d1], brl[d2], σ²ϕ)

      #update branches
      bbX!(X, ipr, brδt[pr], σ²ϕ)
      bbX!(X, id1, brδt[d1], σ²ϕ)
      bbX!(X, id2, brδt[d2], σ²ϕ)

    else
      id1 = bridx[d1]
      id2 = bridx[d2]

      # update node
      X[id1[1]] =
      X[id2[1]] = duoupd(X[id1[end]],
                         X[id2[end]],
                         brl[d1], brl[d2], σ²ϕ)

      # update branches
      bbX!(X, id1, brδt[d1], σ²ϕ)
      bbX!(X, id2, brδt[d2], σ²ϕ)
    end

  end

  return nothing
end

"""
    upbranchX!(j    ::Int64,
               X    ::Array{Float64,2},
               bridx::Array{UnitRange{Int64},1},
               brδt ::Array{Array{Float64,1},1},
               σ²c  ::Float64)

Update a branch j in X using a Brownian bridge.
"""
function upbranchX!(j    ::Int64,
                    X    ::Array{Float64,2},
                    bridx::Array{UnitRange{Int64},1},
                    brδt ::Array{Array{Float64,1},1},
                    σ²ϕ  ::Float64)

  @inbounds bbX!(X, bridx[j], brδt[j], σ²ϕ)

  return nothing
end

"""
    bbX!(X::Array{Float64,2}, idx::UnitRange, t::Array{Float64,1}, σ::Float64)

Brownian bridge simulation function for updating a branch in X in place.
"""
function bbX!(X  ::Array{Float64,2},
              idx::UnitRange,
              t  ::Array{Float64,1},
              σ²ϕ::Float64)

  @inbounds begin

    xf::Float64 = X[idx[end]]

    for i = Base.OneTo(lastindex(t)-1)
      X[idx[i+1]] = (X[idx[i]] + randn()*sqrt((t[i+1] - t[i])*σ²ϕ))::Float64
    end

    invte::Float64 = 1.0/t[end]
    xdif ::Float64 = (X[idx[end]] - xf)

    @simd for i = Base.OneTo(lastindex(t))
      X[idx[i]] = (X[idx[i]] - t[i] * invte * xdif)::Float64
    end
  end

  return nothing
end

max_stars_repo_name: hmorlon/JPANDA

max_issues_repo_name: hmorlon/JPANDA

max_forks_repo_name: hmorlon/JPANDA


text:
using JuLIP, ASE, StaticArrays, NeighbourLists
using Test

X_Ti = vecs([0.0 5.19374 2.59687 3.8953 1.29843 6.49217 7.7906 12.9843 10.3875 11.6859 9.08904 14.2828; 0.0 0.918131 1.83626 -1.11022e-16 0.918131 1.83626 -2.22045e-16 0.918131 1.83626 0.0 0.918131 1.83626; 0.0 0.0 0.0 2.24895 2.24895 2.24895 0.0 0.0 0.0 2.24895 2.24895 2.24895])
C_Ti = (@SMatrix [15.5812 2.47895 0.0; 0.0 2.75439 0.0; 0.0 0.0 4.49791])

# -------------- MatSciPy NeighbourList Patch -------------
using PyCall
import NeighbourLists
matscipy_neighbours = pyimport("matscipy.neighbours")
function asenlist(at::Atoms, rcut)
   pyat = ASEAtoms(at).po
   return matscipy_neighbours[:neighbour_list]("ijdD", pyat, rcut)
end

function matscipy_nlist(at::Atoms{T}, rcut::T; recompute=false, kwargs...) where T <: AbstractFloat
   i, j, r, R = asenlist(at, rcut)
   i = copy(i) .+ 1
   j = copy(j) .+ 1
   r = copy(r)
   R = collect(vecs(copy(R')))
   first = NeighbourLists.get_first(i, length(at))
   NeighbourLists.sort_neigs!(j, r, R, first)
   return NeighbourLists.PairList(positions(at), rcut, i, j, r, R, first)
end
# --------------------------------------------------------

pynlist(at, cutoff) = matscipy_nlist(at, cutoff)
jnlist(at, cutoff) = PairList(positions(at), cutoff, cell(at), pbc(at))

function test_nlist_julip(at, cutoff)
   nlist = jnlist(at, cutoff)
   py_nlist = pynlist(at, cutoff)
   return ( (nlist.i == py_nlist.i) &&
            (nlist.j == py_nlist.j) &&
            (nlist.r ≈ py_nlist.r) &&
            (nlist.R ≈ py_nlist.R) )
end

test_configs = [
   #
   ( "si, cubic, cluster, short",
    set_pbc!(bulk(:Si, cubic=true) * 3, false),
    1.1 * rnn(:Si) ),
   #
   ( "si, cubic, cluster, med",
    set_pbc!(bulk(:Si, cubic=true) * 3, false),
    2.1 * rnn(:Si) ),
   #
   ( "si, non-cubic cell, cluster, med",
    set_pbc!(bulk(:Si) * 5, false),
    2.1 * rnn(:Si) ),
   #
   ( "si, non-cubic cell, pbc",
    set_pbc!(bulk(:Si) * 5, false),
    2.1 * rnn(:Si) ),
   #
   ( "si, non-cubic cell, mixed bc",
    set_pbc!(bulk(:Si) * 5, false),
    2.1 * rnn(:Si) ),
    #
   ("Ti, non-symmetric elongated cell, pbc",
    set_pbc!( set_cell!( Atoms(:Ti, X_Ti), C_Ti ), true ),
    1.45 * rnn(:Ti) ),
   #
   ("Ti (hcp?) canonical cell, pbc",
    set_pbc!( bulk(:Ti) * 4, true ),
    2.3 * rnn(:Ti) ),
   ]

# ----------- A FEW MORE COMPLEX TESTS THAT FAILED AT SOME POINT ---------------

# [1] a left-handed cell orientation
#     the test that failed during Cas' experiments

X = [ 0.00000000e+00  0.00000000e+00  0.00000000e+00
      1.92333044e+00  6.63816518e-17 -1.36000000e+00
      1.92333044e+00  1.92333044e+00 -2.72000000e+00
      3.84666089e+00  1.92333044e+00 -4.08000000e+00 ]'
C = [ 3.84666089   0.           0.
      0.           3.84666089   0.
      0.           0.          -5.44 ]'

# X = [ 0.00000000e+00  0.00000000e+00  0.00000000e+00
#       1.92333044e+00  6.63816518e-17  1.36000000e+00
#       1.92333044e+00  1.92333044e+00  2.72000000e+00
#       3.84666089e+00  1.92333044e+00  4.08000000e+00 ]'
# C = [ 3.84666089   0.           0.
#       0.           3.84666089   0.
#       0.           0.           5.44 ]'
at = Atoms(:Si, collect(vecs(X)))
set_cell!(at, C)
set_pbc!(at, (true,true,true))
atlge = at * (1,1,10)
rcut = 2.3*rnn(:Si)
push!(test_configs, ("Si left-oriented", at, rcut))
push!(test_configs, ("Si left-oriented, large", atlge, rcut))

test_nlist_julip(at, rcut)
test_nlist_julip(atlge, rcut)

# # [2] vacancy in bulk Si
#
# using PyCall
# at = bulk(:Si, cubic=true)
# @pyimport ase.lattice.cubic as cubic
# at2py = cubic.Diamond(symbol = "Si", latticeconstant = 5.43)
# at2py[:get_positions]()
# @test at2py[:get_cell]() ≈ cell(at1)
# @test mat(positions(at1))' ≈ at2py[:get_positions]()[:,[3,2,1]]
#
# at = bulk(:Si, cubic=true)
# at1 = at * 3
# X = mat(positions(at1))
# at1 = deleteat!(at1, 1)
# at2 = deleteat!(set_positions!(at * 3, X[[3,2,1],:]), 1)
# rcut = 2.3 * rnn(:Si)
# test_nlist_julip(at1, rcut)
# test_nlist_julip(at2, rcut)

# [3] Two Ti configuration that seems to be causing problems
C1 = @SMatrix [5.71757 -1.81834e-15 9.74255e-41; -2.85879 4.95156 4.93924e-25; 4.56368e-40 9.05692e-25 9.05629]
X1 = vecs([0.00533847 2.85879 -1.42939 1.42939 0.0 2.85879 -1.42939 1.42939 -1.43e-6 2.85878 -1.42939 1.42939 -1.43e-6 2.85878 -1.42939 1.42939;
          -0.0 -0.0 2.47578 2.47578 0.0 -0.0 2.47578 2.47578 1.65052 1.65052 4.1263 4.1263 1.65052 1.65052 4.1263 4.1263;
          0.00845581 0.0 0.0 0.0 4.52815 4.52815 4.52815 4.52815 2.26407 2.26407 2.26407 2.26407 6.79222 6.79222 6.79222 6.79222]) |> collect
C2 = @SMatrix [5.71757 0.0 0.0; -2.85879 4.95156 0.0; 0.0 0.0 9.05629]
X2 = vecs([0.00534021 2.85879 -1.42939 1.42939 0.0 2.85879 -1.42939 1.42939 0.0 2.85879 -1.4294 1.4294 0.0 2.85879 -1.4294 1.4294;
           0.0 0.0 2.47578 2.47578 0.0 0.0 2.47578 2.47578 1.65052 1.65052 4.1263 4.1263 1.65052 1.65052 4.1263 4.1263;
           0.00845858 0.0 0.0 0.0 4.52815 4.52815 4.52815 4.52815 2.26407 2.26407 2.26407 2.26407 6.79222 6.79222 6.79222 6.79222])  |> collect
at1 = set_cell!(Atoms(:Ti, X1), C1)
at2 = set_cell!(Atoms(:Ti, X2), C2)
rcut = 2.5 * rnn(:Ti)
test_nlist_julip(at1, rcut)
test_nlist_julip(at2, rcut)

# --------------- ACTUALLY RUNNING THE TESTS ------------------

println("JuLIP Configuration tests:")
for (i, (descr, at, cutoff)) in enumerate(test_configs)
  print("TEST $i: $descr => ")
  println(@test test_nlist_julip(at, cutoff))
end


text:
using .VectorizationBase: AbstractSIMD
import .ForwardDiff
import .ChainRulesCore

@generated function SLEEFPirates.tanh_fast(x::ForwardDiff.Dual{T,S,N}) where {T,S,N}
    quote
        $(Expr(:meta,:inline))
        t = tanh_fast(x.value)
        ∂t = vfnmadd_fast(t, t, one(S))
        p = x.partials
        ForwardDiff.Dual(t, ForwardDiff.Partials(Base.Cartesian.@ntuple $N n -> mul_fast(∂t, p[n])))
    end
end
function ChainRulesCore.rrule(::typeof(tanh_fast), x)
    t = tanh_fast(x)
    ∂ = let t = t
        y -> (ChainRulesCore.Zero(), mul_fast(vfnmadd_fast(t, t, one(t)), y))
    end
    t, ∂
end
@generated function SLEEFPirates.sigmoid_fast(x::ForwardDiff.Dual{T,S,N}) where {T,S,N}
    quote
        $(Expr(:meta,:inline))
        s = sigmoid_fast(x.value)
        ∂s = vfnmadd_fast(s,s,s)
        p = x.partials
        ForwardDiff.Dual(s, ForwardDiff.Partials(Base.Cartesian.@ntuple $N n -> mul_fast(∂s, p[n])))
    end
end
function ChainRulesCore.rrule(::typeof(sigmoid_fast), x)
    s = sigmoid_fast(x)
    ∂ = let s = s
        y -> (ChainRulesCore.Zero(), mul_fast(vfnmadd_fast(s, s, s), y))
    end
    s, ∂
end
@generated function VectorizationBase.relu(x::ForwardDiff.Dual{T,S,N}) where {T,S,N}
    quote
        $(Expr(:meta,:inline))
        v = x.value
        z = zero(v)
        cmp = v < z
        r = ifelse(cmp, z, v)
        p = x.partials
        ForwardDiff.Dual(r, ForwardDiff.Partials(Base.Cartesian.@ntuple $N n -> ifelse(cmp, z, p[n])))
    end
end
function ChainRulesCore.rrule(::typeof(relu), v)
    z = zero(v)
    cmp = v < z
    r = ifelse(cmp, z, v)
    ∂ = let cmp = cmp
        y -> (ChainRulesCore.Zero(), ifelse(cmp, zero(y), y))
    end
    r, ∂
end

@generated function init_dual(v::Tuple{Vararg{AbstractSIMD,A}}) where {A}
    res = Expr(:tuple)
    q = Expr(:block, Expr(:meta,:inline))
    for a ∈ 1:A
        v_a = Symbol(:v_,a)
        push!(q.args, Expr(:(=), v_a, Expr(:ref, :v, a)))
        partials = Expr(:tuple)
        for i ∈ 1:A
            push!(partials.args, Expr(:call, i == a ? :one : :zero, v_a))
        end
        push!(res.args, :(ForwardDiff.Dual($v_a, ForwardDiff.Partials($partials))))
    end
    push!(q.args, res)
    q
end

@generated function dual_store!(∂p::Tuple{Vararg{AbstractStridedPointer,A}}, p::AbstractStridedPointer, ∂v, im::Vararg{Any,N}) where {A,N}
    quote
        $(Expr(:meta,:inline))
        v = ∂v.value
        ∂ = ∂v.partials
        Base.Cartesian.@nextract $N im im
        Base.Cartesian.@ncall $N VectorizationBase.vnoaliasstore! p v im  # store
        Base.Cartesian.@nexprs $A a -> begin # for each of `A` partials
            ∂p_a = ∂p[a]
            ∂_a = ∂[a]
            Base.Cartesian.@ncall $N VectorizationBase.vnoaliasstore! ∂p_a ∂_a im # store
        end
        nothing
    end
end

function ∂vmap_singlethread!(
    f::F, ∂y::Tuple{Vararg{DenseArray{T},A}}, y::DenseArray{T},
    args::Vararg{<:DenseArray{<:Base.HWReal},A}
) where {F,T <: Base.HWReal, A}
    N = length(y)
    ptry = VectorizationBase.zero_offsets(stridedpointer(y))
    ptrargs = VectorizationBase.zero_offsets.(stridedpointer.(args))
    ptr∂y = VectorizationBase.zero_offsets.(stridedpointer.(∂y))

    i = 0
    V = VectorizationBase.pick_vector_width(T)
    W = Int(V)
    st = VectorizationBase.static_sizeof(T)
    zero_index = MM{W}(StaticInt(0), st)
    while i < N - ((W << 2) - 1)
        index = VectorizationBase.Unroll{1,W,4,1,W,0x0000000000000000}((i,))
        v = f(init_dual(vload.(ptrargs, index))...)
        dual_store!(ptr∂y, ptry, v, index)
        i = vadd_fast(i, 4W)
    end
    while i < N - (W - 1)
        vᵣ = f(init_dual(vload.(ptrargs, ((MM{W}(i),),)))...)
        dual_store!(ptr∂y, ptry, vᵣ, (MM{W}(i),))
        i = vadd_fast(i, W)
    end
    if i < N
        m = mask(T, N & (W - 1))
        dual_store!(ptr∂y, ptry, f(init_dual(vload.(ptrargs, ((MM{W}(i),),), m))...), (MM{W}(i,),), m)
    end
    nothing
end

struct SIMDMapBack{K,T<:Tuple{Vararg{Any,K}}}
    jacs::T
end
@generated function (b::SIMDMapBack{K,T})(Δ::A) where {K,T,A}
    preloop = Expr(:block, :(jacs = b.jacs))
    loop_body = Expr(:block, :(Δᵢ = Δ[i]))
    ret = Expr(:tuple, ChainRulesCore.Zero(), ChainRulesCore.Zero())
    for k ∈1:K
        jₖ = Symbol(:j_, k)
        push!(preloop.args, :($jₖ = jacs[$k]))
        push!(loop_body.args, :($jₖ[i] *= Δᵢ))
        push!(ret.args, jₖ)
    end
    quote
        $preloop
        @avx for i ∈ eachindex(Δ)
            $loop_body
        end
        $ret
    end
end

function ChainRulesCore.rrule(::typeof(vmap), f::F, args::Vararg{Any,K}) where {F,K}
    out = similar(first(args))
    jacs = map(similar, args)
    ∂vmap_singlethread!(f, jacs, out, args...)
    out, SIMDMapBack(jacs)
end


text:
""" Data structure for undirected wiring diagrams.
"""
module UndirectedWiringDiagrams
export AbstractUndirectedWiringDiagram, UndirectedWiringDiagram,
  outer_box, box, junction, nboxes, njunctions, boxes, junctions, ports,
  ports_with_junction, junction_type, port_type, add_box!, add_junction!,
  add_junctions!, set_junction!, add_wire!, add_wires!, ocompose

using ...CategoricalAlgebra.CSets, ...Present
using ...CategoricalAlgebra.ShapeDiagrams: Span
using ...CategoricalAlgebra.FinSets: FinOrdFunction, pushout
using ...Theories: FreeCategory, dom, codom, compose, ⋅, id

import ..DirectedWiringDiagrams: box, boxes, nboxes, add_box!, add_wire!,
  add_wires!
import ..AlgebraicWiringDiagrams: add_junctions!, ocompose

# Data types
############

@present TheoryUWD(FreeCategory) begin
  Box::Ob
  Port::Ob
  OuterPort::Ob
  Junction::Ob

  box::Hom(Port,Box)
  junction::Hom(Port,Junction)
  outer_junction::Hom(OuterPort,Junction)
end

const AbstractUndirectedWiringDiagram = const AbstractUWD =
  AbstractCSetType(TheoryUWD)
const UntypedUndirectedWiringDiagram = const UntypedUWD =
  CSetType(TheoryUWD, index=[:box, :junction, :outer_junction])

@present TheoryTypedUWD <: TheoryUWD begin
  Type::Ob

  port_type::Hom(Port,Type)
  outer_port_type::Hom(OuterPort,Type)
  junction_type::Hom(Junction,Type)

  compose(junction, junction_type) == port_type
  compose(outer_junction, junction_type) == outer_port_type
end

const TypedUndirectedWiringDiagram = const TypedUWD =
  CSetType(TheoryTypedUWD, data=[:Type],
           index=[:box, :junction, :outer_junction])

# Imperative interface
######################

function UndirectedWiringDiagram(::Type{UWD},
    nports::Int; data_types...) where UWD <: AbstractUWD
  d = UWD(; data_types...)
  add_parts!(d, :OuterPort, nports)
  return d
end
UndirectedWiringDiagram(nports::Int; data_types...) =
  UndirectedWiringDiagram(UntypedUWD, nports; data_types...)

function UndirectedWiringDiagram(::Type{UWD},
    port_types::AbstractVector{T}; data_types...) where {UWD <: AbstractUWD, T}
  d = UWD(; port_type=T, outer_port_type=T, junction_type=T, data_types...)
  nports = length(port_types)
  add_parts!(d, :OuterPort, nports, outer_port_type=port_types)
  return d
end
UndirectedWiringDiagram(port_types::AbstractVector; data_types...) =
  UndirectedWiringDiagram(TypedUWD, port_types; data_types...)

outer_box(::AbstractUWD) = 0
box(d::AbstractUWD, args...) = subpart(d, args..., :box)
junction(d::AbstractUWD, args...; outer::Bool=false) =
  subpart(d, args..., outer ? :outer_junction : :junction)

function junction(d::AbstractUWD, port::Tuple{Int,Int})
  box, nport = port
  box == outer_box(d) ?
    junction(d, nport, outer=true) : junction(d, ports(d, box)[nport])
end

nboxes(d::AbstractUWD) = nparts(d, :Box)
njunctions(d::AbstractUWD) = nparts(d, :Junction)
boxes(d::AbstractUWD) = 1:nboxes(d)
junctions(d::AbstractUWD) = 1:njunctions(d)

ports(d::AbstractUWD; outer::Bool=false) =
  1:nparts(d, outer ? :OuterPort : :Port)
ports(d::AbstractUWD, box) =
  box == outer_box(d) ? (1:nparts(d, :OuterPort)) : incident(d, box, :box)
ports_with_junction(d::AbstractUWD, junction; outer::Bool=false) =
  incident(d, junction, outer ? :outer_junction : :junction)

junction_type(d::AbstractUWD, args...) = subpart(d, args..., :junction_type)
port_type(d::AbstractUWD, args...; outer::Bool=false) =
  subpart(d, args..., outer ? :outer_port_type : :port_type)

function port_type(d::AbstractUWD, port::Tuple{Int,Int})
  box, nport = port
  box == outer_box(d) ?
    port_type(d, nport, outer=true) : port_type(d, ports(d, box)[nport])
end

add_box!(d::AbstractUWD; data...) = add_part!(d, :Box; data...)

function add_box!(d::AbstractUWD, nports::Int; data...)
  box = add_box!(d; data...)
  ports = add_parts!(d, :Port, nports, box=box)
  box
end

function add_box!(d::AbstractUWD, port_types::AbstractVector; data...)
  box = add_box!(d; data...)
  nports = length(port_types)
  ports = add_parts!(d, :Port, nports, box=box, port_type=port_types)
  box
end

add_junction!(d::AbstractUWD; data...) = add_part!(d, :Junction; data...)
add_junction!(d::AbstractUWD, type; data...) =
  add_part!(d, :Junction; junction_type=type, data...)

add_junctions!(d::AbstractUWD, njunctions::Int; data...) =
  add_parts!(d, :Junction, njunctions; data...)
add_junctions!(d::AbstractUWD, types::AbstractVector; data...) =
  add_parts!(d, :Junction, length(types); junction_type=types, data...)

function set_junction!(d::AbstractUWD, port, junction; outer::Bool=false)
  if has_subpart(d, :junction_type)
    ptype, jtype = port_type(d, port, outer=outer), junction_type(d, junction)
    all(ptype .== jtype) || error(
      "Domain error: port type $ptype and junction type $jtype do not match")
  end
  set_subpart!(d, port, outer ? :outer_junction : :junction, junction)
end
set_junction!(d::AbstractUWD, junction; kw...) =
  set_junction!(d, :, junction; kw...)

function set_junction!(d::AbstractUWD, port::Tuple{Int,Int}, junction)
  box, nport = port
  if box == outer_box(d)
    set_junction!(d, nport, junction, outer=true)
  else
    set_junction!(d, ports(d, box)[nport], junction)
  end
end

""" Wire together two ports in an undirected wiring diagram.

A convenience method that creates and sets junctions as needed. Ports are only
allowed to have one junction, so if both ports already have junctions, then the
second port is assigned the junction of the first. The handling of the two
arguments is otherwise symmetric.

FIXME: When both ports already have junctions, the two junctions should be
*merged*. To do this, we must implement `merge_junctions!` and thus also
`rem_part!`.
"""
function add_wire!(d::AbstractUWD, port1::Tuple{Int,Int}, port2::Tuple{Int,Int})
  j1, j2 = junction(d, port1), junction(d, port2)
  if j1 > 0
    set_junction!(d, port2, j1)
  elseif j2 > 0
    set_junction!(d, port1, j2)
  else
    j = has_subpart(d, :junction_type) ?
      add_junction!(d, port_type(d, port1)) : add_junction!(d)
    set_junction!(d, port1, j)
    set_junction!(d, port2, j)
  end
end
add_wire!(d, wire::Pair) = add_wire!(d, first(wire), last(wire))

function add_wires!(d::AbstractUWD, wires)
  for wire in wires
    add_wire!(d, wire)
  end
end

# Operadic interface
####################

function ocompose(f::AbstractUWD, gs::AbstractVector{<:AbstractUWD})
  @assert length(gs) == nboxes(f)
  h = empty(f)
  copy_parts!(h, f, OuterPort=ports(f, outer=true))
  for g in gs
    copy_boxes!(h, g, boxes(g))
  end

  f_junction = FinOrdFunction(
    flat(junction(f, ports(f, i)) for i in boxes(f)), njunctions(f))
  # FIXME: Should use coproduct as monoidal product.
  gs_offset = [0; cumsum(njunctions.(gs))]
  gs_outer = FinOrdFunction(
    flat(junction(g, outer=true) .+ n for (g,n) in zip(gs, gs_offset[1:end-1])),
    gs_offset[end])
  cospan = pushout(Span(f_junction, gs_outer))
  f_inc, g_inc = cospan.left, cospan.right
  junctions = add_junctions!(h, codom(f_inc).n)
  if has_subpart(h, :junction_type)
    set_subpart!(h, [collect(f_inc); collect(g_inc)], :junction_type,
                 [junction_type(f); flat(junction_type(g) for g in gs)])
  end

  f_outer = FinOrdFunction(junction(f, outer=true), njunctions(f))
  # FIXME: Again, should use coproduct.
  gs_junction = FinOrdFunction(
    flat(junction(g) .+ n for (g,n) in zip(gs, gs_offset[1:end-1])),
    gs_offset[end])
  set_junction!(h, collect(f_outer ⋅ f_inc), outer=true)
  set_junction!(h, collect(gs_junction ⋅ g_inc))
  return h
end

function ocompose(f::AbstractUWD, i::Int, g::AbstractUWD)
  @assert 1 <= i <= nboxes(f)
  h = empty(f)
  copy_parts!(h, f, OuterPort=ports(f, outer=true))
  copy_boxes!(h, f, 1:(i-1))
  copy_boxes!(h, g, boxes(g))
  copy_boxes!(h, f, (i+1):nboxes(f))

  f_i = FinOrdFunction(junction(f, ports(f, i)), njunctions(f))
  g_outer = FinOrdFunction(junction(g, outer=true), njunctions(g))
  cospan = pushout(Span(f_i, g_outer))
  f_inc, g_inc = cospan.left, cospan.right
  junctions = add_junctions!(h, codom(f_inc).n)
  if has_subpart(h, :junction_type)
    set_subpart!(h, [collect(f_inc); collect(g_inc)], :junction_type,
                 [junction_type(f); junction_type(g)])
  end

  f_outer = FinOrdFunction(junction(f, outer=true), njunctions(f))
  f_start = FinOrdFunction(junction(f, flat(ports(f, 1:(i-1)))), njunctions(f))
  g_junction = FinOrdFunction(junction(g), njunctions(g))
  f_end = FinOrdFunction(
    junction(f, flat(ports(f, (i+1):nboxes(f)))), njunctions(f))
  set_junction!(h, collect(f_outer ⋅ f_inc), outer=true)
  set_junction!(h, [
    collect(f_start ⋅ f_inc);
    collect(g_junction ⋅ g_inc);
    collect(f_end ⋅ f_inc);
  ])
  return h
end

copy_boxes!(d::AbstractUWD, from::AbstractUWD, boxes) =
  copy_parts!(d, from, Box=boxes, Port=flat(ports(from, boxes)))

flat(vs) = reduce(vcat, vs, init=Int[])

end


text:
function readdata(paths::Array{String, 1}; columns_to_get=[], add_tag::Bool = true, samplestocut::Int = 3)
    formated_data = DataFrame(reshape([], 0, length(columns_to_get)), :auto)
    rename!(formated_data, columns_to_get)

    for file in paths
        formated_data =
            vcat(formated_data, readdata(file, columns_to_get=columns_to_get, add_tag=add_tag, samplestocut=samplestocut))
    end

    return formated_data
end

function readdata(path::String; columns_to_get::Array=[], add_tag::Bool = true, samplestocut::Int = 3)
    df = readdatafromcsv(path)
    meta = getmetadata(df)
    data = getdatawithoutmeta(df)

    lonlat = calculatelonlat(meta)
    data = hcat(data, lonlat)

    if add_tag
        tag_columns = createtagcolumns(path, data)
        data = hcat(data, tag_columns)
    end

    data = cutoutliners(data, samplestocut)
    data[!, :timestep] = 0:nrow(data)-1
    data[!, :timestep] = data[!, :timestep] .* 0.1

    data = data[:, columns_to_get]
end

function readdatafromcsv(path)
    AXIS = ["_x" "_y" "_z"]
    NAMES = [:timestep]
    NAMES = vcat(NAMES, vec(Symbol.(["magnetometr"], AXIS)))
    NAMES = vcat(NAMES, vec(Symbol.(["acceleromter"], AXIS)))
    NAMES = vcat(NAMES, vec(Symbol.(["orientation"], AXIS)))

    df =
        CSV.read(path, DataFrame, header = NAMES, silencewarnings = true, threaded = false)
end

function getmetadata(df)
    METADATA_COLNAMES = [
        "latitude_start",
        "longitude_start",
        "latitude_end",
        "longitude_end",
        "first_sample",
        "last_sample",
    ]

    gapidx = findfirst(occursin.(r"<\d+>", df[!, "timestep"]))
    meta = df[gapidx+1:end, :][:, 1:6]
    rename!(meta, Symbol.(METADATA_COLNAMES))
    meta[!, "latitude_start"] = parse.(Float64, meta[!, "latitude_start"])
    meta[!, "first_sample"] = trunc.(Int, meta[!, "first_sample"])
    meta[!, "last_sample"] = trunc.(Int, meta[!, "last_sample"])
    meta = DataFrame(meta)
end

function getdatawithoutmeta(df)
    gapidx = findfirst(occursin.(r"<\d+>", df[!, "timestep"]))
    data = df[1:gapidx-1, :]
end

function calculatelonlat(meta)
    samples_num = last(meta, 1)[!, "last_sample"][1] + 1
    result = DataFrame(lat = 1:samples_num, lon = 1:samples_num)
    result[!, "lat"] = convert.(Float64, result[!, "lat"])
    result[!, "lon"] = convert.(Float64, result[!, "lon"])

    for corridor in eachrow(meta)
        first_sample_idx = corridor["first_sample"] + 1 # This dataset is 0-index-based
        last_sample_idx = corridor["last_sample"] + 1
        corridor_idxs = first_sample_idx:last_sample_idx

        corridor_samples_num = last_sample_idx - first_sample_idx + 1

        lat_diff = corridor["latitude_end"] - corridor["latitude_start"]
        lon_diff = corridor["longitude_end"] - corridor["longitude_start"]

        lat_step = lat_diff / corridor_samples_num
        lon_step = lon_diff / corridor_samples_num

        result[corridor_idxs, :lon] =
            (result[corridor_idxs, :lon] .- first_sample_idx) * lon_step .+
            corridor["longitude_start"]
        result[corridor_idxs, :lat] =
            (result[corridor_idxs, :lat] .- first_sample_idx) * lat_step .+
            corridor["latitude_start"]
    end

    return result
end

function createtagcolumns(path, data)
    n_rows = nrow(data)
    tag = match(r"(.*[/\\])*(?<id>.*)_(?<sample>.*)\.txt", path)
    if !isnothing(tag)
        return  DataFrame(
            path_id = repeat([tag[:id]], n_rows),
            path_sample = repeat([tag[:sample]], n_rows)
        )
    else
        tag = match(r"(.*[/\\])*(?<pathid>.*)\.txt", path)
        return DataFrame(
            path_id = repeat([tag[:pathid]], n_rows)
        )
    end
end

function cutoutliners(data, samplestocut)
    data = data[setdiff(1:end, 1:samplestocut), :]
    data = data[setdiff(1:end, end-samplestocut:end), :]
end

max_stars_repo_name: SebastianSuchowiak/MagneticLocSuchowiak

max_issues_repo_name: SebastianSuchowiak/MagneticLocSuchowiak

max_forks_repo_name: SebastianSuchowiak/MagneticLocSuchowiak


text:
function check_tolerance(::Type{Sphere2D}, arg0::jdouble)
    return jcall(Sphere2D, "checkTolerance", void, (jdouble,), arg0)
end

function get_dimension(obj::Sphere2D)
    return jcall(obj, "getDimension", jint, ())
end

function get_instance(::Type{Sphere2D})
    return jcall(Sphere2D, "getInstance", Sphere2D, ())
end

function get_sub_space(obj::Sphere2D)
    return jcall(obj, "getSubSpace", Sphere1D, ())
end


text:
@testset "Perfect bipartite matching" begin
  @testset "Uncorrelated" begin
    @testset "Constructor with $i nodes on each side" for i in [2, 5, 10]
      n = i
      ε = (i - 2) / 16
      reward = Distribution[Bernoulli(.5 + ((i == j) ? ε : 0.)) for i in 1:n, j in 1:n]
      instance = UncorrelatedPerfectBipartiteMatching(reward, PerfectBipartiteMatchingNoSolver())

      @test instance.n_arms == n ^ 2
      @test instance.reward == reward

      # Error: non bipartite.
      reward = Distribution[Bernoulli(.5 + ((i == j) ? ε : 0.)) for i in 1:n, j in 1:(n + 2)]
      @test_throws ErrorException UncorrelatedPerfectBipartiteMatching(reward, PerfectBipartiteMatchingNoSolver())

      reward = Distribution[Bernoulli(.5 + ((i == j) ? ε : 0.)) for i in 1:(n + 2), j in 1:n]
      @test_throws ErrorException UncorrelatedPerfectBipartiteMatching(reward, PerfectBipartiteMatchingNoSolver())
    end

    @testset "State with $i nodes on each side" for i in [2, 5, 10]
      n = i
      ε = (i - 2) / 16
      reward = Distribution[Bernoulli(.5 + ((i == j) ? ε : .0)) for i in 1:n, j in 1:n]
      instance = UncorrelatedPerfectBipartiteMatching(reward, PerfectBipartiteMatchingNoSolver())

      state = initial_state(instance)

      @test state.round == 0
      @test state.regret == 0.0
      @test state.reward == 0.0
      @test length(state.arm_counts) == n * n
      @test length(state.arm_reward) == n * n
      @test length(state.arm_average_reward) == n * n

      for i in 1:n
        for j in 1:n
          @test state.arm_counts[(i, j)] == 0
          @test state.arm_reward[(i, j)] == 0.0
          @test state.arm_average_reward[(i, j)] == 0.0
        end
      end
    end

    @testset "Trace with $i nodes on each side" for i in [2, 5, 10]
      n = i
      ε = (i - 2) / 16
      reward = Distribution[Bernoulli(.5 + ((i == j) ? ε : .0)) for i in 1:n, j in 1:n]
      instance = UncorrelatedPerfectBipartiteMatching(reward, PerfectBipartiteMatchingNoSolver())

      trace = initial_trace(instance)

      @test length(trace.states) == 0
      @test length(trace.arms) == 0
      @test length(trace.reward) == 0
      @test length(trace.policy_details) == 0
      @test length(trace.time_choose_action) == 0

      @test eltype(trace.states) == State{Tuple{Int, Int}}
      @test eltype(trace.arms) == Vector{Tuple{Int, Int}}
      @test eltype(trace.reward) == Vector{Float64}
      @test eltype(trace.time_choose_action) == Int
    end

    @testset "Pull with $i nodes on each side" for i in [2, 5, 10]
      n = i
      reward = Distribution[Bernoulli(((i == j) ? 0.0 : 1.0)) for i in 1:n, j in 1:n]
      instance = UncorrelatedPerfectBipartiteMatching(reward, PerfectBipartiteMatchingNoSolver())

      Random.seed!(1)
      @test pull(instance, [(1, 2), (i, i)]) == ([1.0, 0.0], -1) # Reward and regret
    end

    @testset "Check feasibility with 3 nodes on each side" begin
      n = 3
      reward = Distribution[Bernoulli(((i == j) ? 0.0 : 1.0)) for i in 1:n, j in 1:n]
      instance = UncorrelatedPerfectBipartiteMatching(reward, PerfectBipartiteMatchingNoSolver())

      @test is_feasible(instance, [(1, 2), (2, 3), (3, 1)])
      @test ! is_feasible(instance, [(2, 3), (3, 2), (1, 2)])
      @test ! is_feasible(instance, [(2, 3), (3, 2), (2, 1), (1, 2)])
    end

    if ! is_travis
      @testset "LP solver" begin
        @testset "Constructor" for i in [2, 5, 10]
          n = i
          reward = Distribution[Bernoulli(((i == j) ? 0.0 : 1.0)) for i in 1:n, j in 1:n]
          instance = UncorrelatedPerfectBipartiteMatching(reward, PerfectBipartiteMatchingLPSolver(Gurobi.Optimizer))

          @test instance.solver != nothing
          @test instance.solver.model != nothing
          @test size(instance.solver.x, 1) == n * n
        end

        @testset "Solve with $i nodes on each side" for i in [2, 5, 10]
          n = i
          reward = Distribution[Bernoulli(((i == j) ? 0.0 : 1.0)) for i in 1:n, j in 1:n]
          instance = UncorrelatedPerfectBipartiteMatching(reward, PerfectBipartiteMatchingLPSolver(Gurobi.Optimizer))

          Random.seed!(i)
          drawn = Dict((i, j) => rand() for i in 1:n, j in 1:n)
          solution = solve_linear(instance, drawn)
          @test is_feasible(instance, solution)
        end
      end
    end

    @testset "Munkres solver" begin
      @testset "Constructor" for i in [2, 5, 10]
        n = i
        reward = Distribution[Bernoulli(((i == j) ? 0.0 : 1.0)) for i in 1:n, j in 1:n]
        instance = UncorrelatedPerfectBipartiteMatching(reward, PerfectBipartiteMatchingMunkresSolver())

        @test instance.solver != nothing
      end

      @testset "Solve with $i nodes on each side" for i in [2, 5, 10]
        n = i
        reward = Distribution[Bernoulli(((i == j) ? 0.0 : 1.0)) for i in 1:n, j in 1:n]
        instance = UncorrelatedPerfectBipartiteMatching(reward, PerfectBipartiteMatchingMunkresSolver())

        Random.seed!(i)
        drawn = Dict((i, j) => rand() for i in 1:n, j in 1:n)
        solution = solve_linear(instance, drawn)
        @test is_feasible(instance, solution)
      end
    end

    @testset "Hungarian solver" begin
      @testset "Constructor" for i in [2, 5, 10]
        n = i
        reward = Distribution[Bernoulli(((i == j) ? 0.0 : 1.0)) for i in 1:n, j in 1:n]
        instance = UncorrelatedPerfectBipartiteMatching(reward, PerfectBipartiteMatchingHungarianSolver())

        @test instance.solver != nothing
      end

      @testset "Solve with $i nodes on each side" for i in [2, 5, 10]
        n = i
        reward = Distribution[Bernoulli(((i == j) ? 0.0 : 1.0)) for i in 1:n, j in 1:n]
        instance = UncorrelatedPerfectBipartiteMatching(reward, PerfectBipartiteMatchingHungarianSolver())

        Random.seed!(i)
        drawn = Dict((i, j) => rand() for i in 1:n, j in 1:n)
        solution = solve_linear(instance, drawn)
        @test is_feasible(instance, solution)
      end
    end

    @testset "Solver equivalence (size: $i nodes on each side)" for i in [2, 5, 10]
      n = i
      reward = Distribution[Bernoulli(((i == j) ? 0.0 : 1.0)) for i in 1:n, j in 1:n]
      instance_lp = ! is_travis && UncorrelatedPerfectBipartiteMatching(reward, PerfectBipartiteMatchingLPSolver(Gurobi.Optimizer))
      instance_munkres = UncorrelatedPerfectBipartiteMatching(reward, PerfectBipartiteMatchingMunkresSolver())
      instance_hungarian = UncorrelatedPerfectBipartiteMatching(reward, PerfectBipartiteMatchingHungarianSolver())

      Random.seed!(i)
      drawn = Dict((i, j) => rand() for i in 1:n, j in 1:n)

      if ! is_travis
        solution_lp = solve_linear(instance_lp, drawn)
        @test is_feasible(instance_lp, solution_lp)
      end

      solution_munkres = solve_linear(instance_munkres, drawn)
      @test is_feasible(instance_munkres, solution_munkres)

      solution_hungarian = solve_linear(instance_hungarian, drawn)
      @test is_feasible(instance_hungarian, solution_hungarian)

      # All solutions must have the same length, as these are perfect matchings.
      @test length(solution_hungarian) == length(solution_munkres)

      cost_munkres = sum(drawn[o] for o in solution_munkres)
      cost_hungarian = sum(drawn[o] for o in solution_hungarian)

      @test cost_hungarian ≈ cost_munkres
      if ! is_travis
        @test length(solution_lp) == length(solution_hungarian)
        cost_lp = sum(drawn[o] for o in solution_lp)
        @test cost_lp ≈ cost_hungarian
      end
    end
  end

  @testset "Correlated" begin
    @testset "Constructor with $i nodes on each side" for i in [2, 5, 10]
      n = i
      ε = (i - 2) / 16
      μ = vec(Float64[.5 + ((i == j) ? ε : 0.) for i in 1:n, j in 1:n])
      Σ = vec(Float64[((abs(i - j) <= 1) ? sign(i - j) : 0.) for i in 1:n, j in 1:n])
      reward = MvNormal(μ, Σ)
      instance = CorrelatedPerfectBipartiteMatching(reward, PerfectBipartiteMatchingNoSolver())

      @test instance.n_arms == n ^ 2
      @test instance.reward == reward
    end

    @testset "State with $i nodes on each side" for i in [2, 5, 10]
      n = i
      ε = (i - 2) / 16
      μ = vec(Float64[.5 + ((i == j) ? ε : 0.) for i in 1:n, j in 1:n])
      Σ = vec(Float64[((abs(i - j) <= 1) ? sign(i - j) : 0.) for i in 1:n, j in 1:n])
      reward = MvNormal(μ, Σ)
      instance = CorrelatedPerfectBipartiteMatching(reward, PerfectBipartiteMatchingNoSolver())

      state = initial_state(instance)

      @test state.round == 0
      @test state.regret == 0.0
      @test state.reward == 0.0
      @test length(state.arm_counts) == n * n
      @test length(state.arm_reward) == n * n
      @test length(state.arm_average_reward) == n * n

      for i in 1:n
        for j in 1:n
          @test state.arm_counts[(i, j)] == 0
          @test state.arm_reward[(i, j)] == 0.0
          @test state.arm_average_reward[(i, j)] == 0.0
        end
      end
    end

    @testset "Trace with $i nodes on each side" for i in [2, 5, 10]
      n = i
      ε = (i - 2) / 16
      μ = vec(Float64[.5 + ((i == j) ? ε : 0.) for i in 1:n, j in 1:n])
      Σ = vec(Float64[((abs(i - j) <= 1) ? sign(i - j) : 0.) for i in 1:n, j in 1:n])
      reward = MvNormal(μ, Σ)
      instance = CorrelatedPerfectBipartiteMatching(reward, PerfectBipartiteMatchingNoSolver())

      trace = initial_trace(instance)

      @test length(trace.states) == 0
      @test length(trace.arms) == 0
      @test length(trace.reward) == 0
      @test length(trace.policy_details) == 0
      @test length(trace.time_choose_action) == 0

      @test eltype(trace.states) == State{Tuple{Int, Int}}
      @test eltype(trace.arms) == Vector{Tuple{Int, Int}}
      @test eltype(trace.reward) == Vector{Float64}
      @test eltype(trace.time_choose_action) == Int
    end

    @testset "Pull with $i nodes on each side" for i in [2]#, 5, 10]
      n = i
      μ = vec(Float64[.5 + ((i == j) ? 1. : 0.) for i in 1:n, j in 1:n])
      Σ = vec(Float64[((abs(i - j) <= 1) ? sign(i - j) : 0.) for i in 1:n, j in 1:n])
      reward = MvNormal(μ, Σ)
      instance = CorrelatedPerfectBipartiteMatching(reward, PerfectBipartiteMatchingNoSolver())

      Random.seed!(1)
      rewards, regret = pull(instance, [(1, 2), (i, i)])
      @test rewards ≈ [0.5, 1.5] atol=1.e-6
      @test regret ≈ -0.5 atol=1.e-6
    end

    @testset "Check feasibility with 3 nodes on each side" begin
      n = 3
      μ = vec(Float64[.5 + ((i == j) ? 1. : 0.) for i in 1:n, j in 1:n])
      Σ = vec(Float64[((abs(i - j) <= 1) ? sign(i - j) : 0.) for i in 1:n, j in 1:n])
      reward = MvNormal(μ, Σ)
      instance = CorrelatedPerfectBipartiteMatching(reward, PerfectBipartiteMatchingNoSolver())

      @test is_feasible(instance, [(1, 2), (2, 3), (3, 1)])
      @test ! is_feasible(instance, [(2, 3), (3, 2), (1, 2)])
      @test ! is_feasible(instance, [(2, 3), (3, 2), (2, 1), (1, 2)])
    end

    @testset "Check partial acceptability with 3 nodes on each side" begin
      n = 3
      μ = vec(Float64[.5 + ((i == j) ? 1. : 0.) for i in 1:n, j in 1:n])
      Σ = vec(Float64[((abs(i - j) <= 1) ? sign(i - j) : 0.) for i in 1:n, j in 1:n])
      reward = MvNormal(μ, Σ)
      instance = CorrelatedPerfectBipartiteMatching(reward, PerfectBipartiteMatchingNoSolver())

      @test is_partially_acceptable(instance, [(1, 2)])
      @test is_partially_acceptable(instance, [(1, 2), (2, 3)])
      @test is_partially_acceptable(instance, [(1, 2), (2, 3), (3, 1)])
      @test ! is_partially_acceptable(instance, [(2, 3), (3, 2), (1, 2)])
      @test ! is_partially_acceptable(instance, [(2, 3), (3, 2), (2, 1), (1, 2)])
    end

    if ! is_travis
      @testset "LP solver" begin
        @testset "Constructor" for i in [2, 5, 10]
          n = i
          μ = vec(Float64[.5 + ((i == j) ? 1. : 0.) for i in 1:n, j in 1:n])
          Σ = vec(Float64[((abs(i - j) <= 1) ? sign(i - j) : 0.) for i in 1:n, j in 1:n])
          reward = MvNormal(μ, Σ)
          instance = CorrelatedPerfectBipartiteMatching(reward, PerfectBipartiteMatchingLPSolver(Gurobi.Optimizer))

          @test instance.solver != nothing
          @test instance.solver.model != nothing
          @test size(instance.solver.x, 1) == n * n
        end

        @testset "Solve with $i nodes on each side" for i in [2, 5, 10]
          n = i
          μ = vec(Float64[.5 + ((i == j) ? 1. : 0.) for i in 1:n, j in 1:n])
          Σ = vec(Float64[((abs(i - j) <= 1) ? sign(i - j) : 0.) for i in 1:n, j in 1:n])
          reward = MvNormal(μ, Σ)
          instance = CorrelatedPerfectBipartiteMatching(reward, PerfectBipartiteMatchingLPSolver(Gurobi.Optimizer))

          Random.seed!(i)
          drawn = Dict((i, j) => rand() for i in 1:n, j in 1:n)
          solution = solve_linear(instance, drawn)
          @test is_feasible(instance, solution)
        end
      end
    end

    @testset "Munkres solver" begin
      @testset "Constructor" for i in [2, 5, 10]
        n = i
        μ = vec(Float64[.5 + ((i == j) ? 1. : 0.) for i in 1:n, j in 1:n])
        Σ = vec(Float64[((abs(i - j) <= 1) ? sign(i - j) : 0.) for i in 1:n, j in 1:n])
        reward = MvNormal(μ, Σ)
        instance = CorrelatedPerfectBipartiteMatching(reward, PerfectBipartiteMatchingMunkresSolver())

        @test instance.solver != nothing
      end

      @testset "Solve with $i nodes on each side" for i in [2, 5, 10]
        n = i
        μ = vec(Float64[.5 + ((i == j) ? 1. : 0.) for i in 1:n, j in 1:n])
        Σ = vec(Float64[((abs(i - j) <= 1) ? sign(i - j) : 0.) for i in 1:n, j in 1:n])
        reward = MvNormal(μ, Σ)
        instance = CorrelatedPerfectBipartiteMatching(reward, PerfectBipartiteMatchingMunkresSolver())

        Random.seed!(i)
        drawn = Dict((i, j) => rand() for i in 1:n, j in 1:n)
        solution = solve_linear(instance, drawn)
        @test is_feasible(instance, solution)
      end
    end

    @testset "Hungarian solver" begin
      @testset "Constructor" for i in [2, 5, 10]
        n = i
        μ = vec(Float64[.5 + ((i == j) ? 1. : 0.) for i in 1:n, j in 1:n])
        Σ = vec(Float64[((abs(i - j) <= 1) ? sign(i - j) : 0.) for i in 1:n, j in 1:n])
        reward = MvNormal(μ, Σ)
        instance = CorrelatedPerfectBipartiteMatching(reward, PerfectBipartiteMatchingHungarianSolver())

        @test instance.solver != nothing
      end

      @testset "Solve with $i nodes on each side" for i in [2, 5, 10]
        n = i
        μ = vec(Float64[.5 + ((i == j) ? 1. : 0.) for i in 1:n, j in 1:n])
        Σ = vec(Float64[((abs(i - j) <= 1) ? sign(i - j) : 0.) for i in 1:n, j in 1:n])
        reward = MvNormal(μ, Σ)
        instance = CorrelatedPerfectBipartiteMatching(reward, PerfectBipartiteMatchingHungarianSolver())

        Random.seed!(i)
        drawn = Dict((i, j) => rand() for i in 1:n, j in 1:n)
        solution = solve_linear(instance, drawn)
        @test is_feasible(instance, solution)
      end
    end

    @testset "Solver equivalence (size: $i nodes on each side)" for i in [2, 5, 10]
      n = i
      μ = vec(Float64[.5 + ((i == j) ? 1. : 0.) for i in 1:n, j in 1:n])
      Σ = vec(Float64[((abs(i - j) <= 1) ? sign(i - j) : 0.) for i in 1:n, j in 1:n])
      reward = MvNormal(μ, Σ)

      instance_lp = ! is_travis && CorrelatedPerfectBipartiteMatching(reward, PerfectBipartiteMatchingLPSolver(Gurobi.Optimizer))
      instance_munkres = CorrelatedPerfectBipartiteMatching(reward, PerfectBipartiteMatchingMunkresSolver())
      instance_hungarian = CorrelatedPerfectBipartiteMatching(reward, PerfectBipartiteMatchingHungarianSolver())

      Random.seed!(i)
      drawn = Dict((i, j) => rand() for i in 1:n, j in 1:n)

      if ! is_travis
        solution_lp = solve_linear(instance_lp, drawn)
        @test is_feasible(instance_lp, solution_lp)
      end

      solution_munkres = solve_linear(instance_munkres, drawn)
      @test is_feasible(instance_munkres, solution_munkres)

      solution_hungarian = solve_linear(instance_hungarian, drawn)
      @test is_feasible(instance_hungarian, solution_hungarian)

      # All solutions must have the same length, as these are perfect matchings.
      @test length(solution_hungarian) == length(solution_munkres)

      cost_munkres = sum(drawn[o] for o in solution_munkres)
      cost_hungarian = sum(drawn[o] for o in solution_hungarian)

      @test cost_hungarian ≈ cost_munkres

      if ! is_travis
        @test length(solution_lp) == length(solution_hungarian)
        cost_lp = sum(drawn[o] for o in solution_lp)
        @test cost_lp ≈ cost_hungarian
      end
    end
  end
end


text:
function uniformMutation(genotype::IntegerGenotype, rng::Random.AbstractRNG,
                         min::Integer, max::Integer, nGens::Integer=1)

    if min > max
        min, max = max, min
    end

    indRep = genotype._representation
    genLen = length(indRep)
    mutatedIndRep = Array{Integer}(undef, genLen)

    if nGens > genLen
        nGens = genLen
    end

    indexes = randomIndexSelection(genLen, nGens, rng)

    mutatedIndRep = copy(indRep)
    range = max-min

    for i=1:nGens
        mutatedIndRep[indexes[i]] = (rand(rng, UInt64)%range)+1+min
    end

    return IntegerGenotype(mutatedIndRep)
end

max_stars_repo_name: SergioGmezM/GeneticJulia

max_issues_repo_name: SergioGmezM/GeneticJulia

max_forks_repo_name: SergioGmezM/GeneticJulia


text:
md"""
# Model Pruning Example

While NaiveNASflux does not come with any built in search policies, it is still possible to do some cool stuff with it.
Below is a very simple example of parameter pruning.

First we need some boilerplate to create the model and do the training:

"""
@testset "Pruning xor example" begin #src
using NaiveNASflux, Flux, Test
using Flux: train!, mse
import Random
Random.seed!(0)
niters = 50

# To cut down on the verbosity, start by making a  helper function for creating a Dense layer as a graph vertex.
# The keyword argument `layerfun=`[`ActivationContribution`](@ref) will wrap the layer and compute an activity
# based neuron utility metric for it while the model trains.
densevertex(in, outsize, act) = fluxvertex(Dense(nout(in),outsize, act), in, layerfun=ActivationContribution)

# Ok, lets create the model and train it. We overparameterize quite heavily to avoid sporadic test failures :)
invertex = denseinputvertex("input", 2)
layer1 = densevertex(invertex, 32, relu)
layer2 = densevertex(layer1, 1, sigmoid)
original = CompGraph(invertex, layer2)

## Training params, nothing to see here
opt = ADAM(0.1)
loss(g) = (x, y) -> mse(g(x), y)

## Training data: xor truth table: y = xor(x) just so we don't need to download a dataset.
x = Float32[0 0 1 1;
            0 1 0 1]
y = Float32[0 1 1 0]

## Train the model
train!(loss(original), params(original), Iterators.repeated((x,y), niters), opt)
@test loss(original)(x, y) < 0.001

# With that out of the way, lets try three different ways to prune the hidden layer (vertex nr 2 in the graph).
# To make examples easier to compare, lets decide up front that we want to remove half of the hidden layer neurons
# and try out three different ways of how to select which ones to remove.
nprune = 16

# Prune the neurons with lowest utility according to the metric in [`ActivationContribution`](@ref).
# This is the default if no utility function is provided.
pruned_least = deepcopy(original)
Δnout!(pruned_least[2] => -nprune)

# Prune the neurons with higest utility according to the metric in [`ActivationContribution`](@ref).
# This is obviously not a good idea if you want to preserve the accuracy.
pruned_most = deepcopy(original)
Δnout!(pruned_most[2] => -nprune) do v
    vals = NaiveNASlib.defaultutility(v)
    return 2*sum(vals) .- vals # Ensure all values are still > 0, even for last vertex
end

# Prune randomly selected neurons by giving random utility.
pruned_random = deepcopy(original)
Δnout!(v -> rand(nout(v)), pruned_random[2] => -nprune)

# Free lunch anyone?
@test   loss(pruned_most)(x, y)   >
        loss(pruned_random)(x, y) >
        loss(pruned_least)(x, y)  >=
        loss(original)(x, y)

# The metric calculated by [`ActivationContribution`](@ref) is actually quite good in this case.
@test loss(pruned_least)(x, y) ≈ loss(original)(x, y) atol = 1e-5
end #src


text:
# x(t+1) = A*x(t) + B*u(t)
struct Dynamics
    A
    B
    vx_range
    vy_range
    d
end

# ϕ₁*x + ϕ₂*y + ϕ₃*z ≥ -þ
# ϕ₁*x + ϕ₂*y + ϕ₃*z ≥ -þ == -ϕ₁*x - ϕ₂*y - ϕ₃*z ≤ þ
struct Ineq
    ϕ₁::Float64
    ϕ₂::Float64
    ϕ₃::Float64
    þ::Float64
end

struct Region
    name
    r::Vector{Ineq}
    θ::Float64
end

mutable struct Funnel
    name
    params
    pos_prec
    neg_prec
    continuous_prec
    dynamics
    pos_eff
    neg_eff
    end_region
    is_continuous
    cont_ind
    function Funnel(name)
        new(name, [],[],[],[],[],[],[],Dict(),true, 1)
    end
end

mutable struct Graph
    num_levels::Int64
    acts
    μacts
    props
    μprops
    leveled
    initprops
    goalprops
    indexes #[object_index, place_pose_index]
    safe_poses
    has_placement_constraint
    function Graph()
        new(0, Dict(1=>[]), Dict(1=>[]), Dict(1=>Dict()),  Dict(1=>[]), false, Dict(), Dict(), [1,1],[],true)
    end
end


text:
@testset "SVector" begin
    @testset "Inner Constructors" begin
        @test SVector{1,Int}((1,)).data === (1,)
        @test SVector{1,Float64}((1,)).data === (1.0,)
        @test SVector{2,Float64}((1, 1.0)).data === (1.0, 1.0)

        @test_throws Exception SVector{1,Int}()
        @test_throws Exception SVector{2,Int}((1,))
        @test_throws Exception SVector{1,Int}(())
        @test_throws Exception SVector{Int,1}((1,))
    end

    @testset "Outer constructors and macro" begin
        @test SVector{1}((1,)).data === (1,)
        @test SVector{1}((1.0,)).data === (1.0,)
        @test SVector((1,)).data === (1,)
        @test SVector((1.0,)).data === (1.0,)

        @test ((@SVector [1.0])::SVector{1}).data === (1.0,)
        @test ((@SVector [1, 2, 3])::SVector{3}).data === (1, 2, 3)
        @test ((@SVector Float64[1,2,3])::SVector{3}).data === (1.0, 2.0, 3.0)
        @test ((@SVector [i for i = 1:3])::SVector{3}).data === (1, 2, 3)
        @test ((@SVector Float64[i for i = 1:3])::SVector{3}).data === (1.0, 2.0, 3.0)

        @test ((@SVector zeros(2))::SVector{2, Float64}).data === (0.0, 0.0)
        @test ((@SVector ones(2))::SVector{2, Float64}).data === (1.0, 1.0)
        @test isa(@SVector(rand(2)), SVector{2, Float64})
        @test isa(@SVector(randn(2)), SVector{2, Float64})

        @test ((@SVector zeros(Float32, 2))::SVector{2,Float32}).data === (0.0f0, 0.0f0)
        @test ((@SVector ones(Float32, 2))::SVector{2,Float32}).data === (1.0f0, 1.0f0)
        @test isa(@SVector(rand(Float32, 2)), SVector{2, Float32})
        @test isa(@SVector(randn(Float32, 2)), SVector{2, Float32})
    end

    @testset "Methods" begin
        v = @SVector [11, 12, 13]

        @test isimmutable(v) == true

        @test v[1] === 11
        @test v[2] === 12
        @test v[3] === 13

        @test Tuple(v) === (11, 12, 13)

        @test size(v) === (3,)
        @test size(typeof(v)) === (3,)
        @test size(SVector{3}) === (3,)

        @test size(v, 1) === 3
        @test size(v, 2) === 1
        @test size(typeof(v), 1) === 3
        @test size(typeof(v), 2) === 1

        @test length(v) === 3

        @test_throws Exception v[1] = 1
    end
end


text:
using DashBootstrapComponents

form = dbc_row(
    [
        dbc_col(
            [
                dbc_label("Email", html_for = "example-email-grid"),
                dbc_input(
                    type = "email",
                    id = "example-email-grid",
                    placeholder = "Enter email",
                ),
            ],
            width = 6,
        ),
        dbc_col(
            [
                dbc_label("Password", html_for = "example-password-grid"),
                dbc_input(
                    type = "password",
                    id = "example-password-grid",
                    placeholder = "Enter password",
                ),
            ],
            width = 6,
        ),
    ],
    className = "g-3",
)

max_stars_repo_name: glsdown/dash-bootstrap-components

max_issues_repo_name: glsdown/dash-bootstrap-components

max_forks_repo_name: glsdown/dash-bootstrap-components


text:
"""
$(TYPEDEF)

Hierachy of AbstractSoilVC:
- [`BrooksCorey`](@ref)
- [`VanGenuchten`](@ref)
"""
abstract type AbstractSoilVC{FT<:AbstractFloat} end

#######################################################################################################################################################################################################
#
# Changes to this structure
# General
#     2021-Sep-30: define this structure with no default constructor
#     2021-Sep-30: define a constructor to fil the value from VanGenuchten parameters
#
#######################################################################################################################################################################################################
"""
$(TYPEDEF)

Brooks Corey soil parameters

# Fields
$(TYPEDFIELDS)

# Examples
```julia
bc = BrooksCorey{FT}("Test", FT(5), FT(2), FT(0.5), FT(0.1));
bc = BrooksCorey{FT}(VanGenuchten{FT}("Loam"));
```
"""
struct BrooksCorey{FT<:AbstractFloat} <:AbstractSoilVC{FT}
    "Soil type"
    TYPE::String
    "Soil b"
    B::FT
    "Potential at saturation `[MPa]`"
    Ψ_SAT::FT
    "Saturated soil volumetric water content"
    Θ_SAT::FT
    "Residual soil volumetric water content"
    Θ_RES::FT
end

#######################################################################################################################################################################################################
#
# Changes to this structure
# General
#     2021-Sep-30: define this structure with two default constructors from an incomplete parameter set
#
#######################################################################################################################################################################################################
"""
$(TYPEDEF)

van Genuchten soil parameters

# Fields
$(TYPEDFIELDS)

# Examples
```julia
vg = VanGenuchten{FT}("Loam");
vg = VanGenuchten{FT}("Silt");
vg = VanGenuchten{FT}("Test", FT(100), FT(2), FT(0.5), FT(0.1));
```
"""
struct VanGenuchten{FT<:AbstractFloat} <:AbstractSoilVC{FT}
    "Soil type"
    TYPE::String
    "Soil α is related to the inverse of the air entry suction, α > 0"
    Α::FT
    "Soil n is Measure of the pore-size distribution"
    N::FT
    "Soil m = 1 - 1/n"
    M::FT
    "Saturated soil volumetric water content"
    Θ_SAT::FT
    "Residual soil volumetric water content"
    Θ_RES::FT

    # constructors
    VanGenuchten{FT}(name::String, α::FT, n::FT, θ_sat::FT, θ_res::FT) where {FT<:AbstractFloat} = new{FT}(name, α, n, 1-1/n, θ_sat, θ_res)

    VanGenuchten{FT}(name::String) where {FT<:AbstractFloat} = (
        # Parameters from Silt soil
        _p = [ 163.2656, 1.37, 0.46, 0.034];

        # switch name
        if name=="Sand"
            _p = [1479.5945, 2.68, 0.43, 0.045];
        elseif name=="Loamy Sand"
            _p = [1265.3084, 2.28, 0.41, 0.057];
        elseif name=="Sandy Loam"
            _p = [ 765.3075, 1.89, 0.41, 0.065];
        elseif name=="Loam"
            _p = [ 367.3476, 1.56, 0.43, 0.078];
        elseif name=="Sandy Clay Loam"
            _p = [ 602.0419, 1.48, 0.39, 0.100];
        elseif name=="Silt Loam"
            _p = [ 204.0820, 1.41, 0.45, 0.067];
        elseif name=="Silt"
            _p = [ 163.2656, 1.37, 0.46, 0.034];
        elseif name=="Clay Loam"
            _p = [ 193.8779, 1.31, 0.41, 0.095];
        elseif name=="Silty Clay Loam"
            _p = [ 102.0410, 1.23, 0.43, 0.089];
        elseif name== "Sandy Clay"
            _p = [ 275.5107, 1.23, 0.38, 0.100];
        elseif name=="Silty Clay"
            _p = [  51.0205, 1.09, 0.36, 0.070];
        elseif name=="Clay"
            _p = [  81.6328, 1.09, 0.38, 0.068];
        else
            @warn "Soil type $(name) not recognized, use Silt instead.";
            name = "Silt";
        end;

        # return a new struct
        return new{FT}(name, _p[1], _p[2], 1-1/_p[2], _p[3], _p[4])
    );
end

#######################################################################################################################################################################################################
#
# Changes to this structure
# General
#     2021-Sep-30: define this structure with no constructors
#     2021-Oct-19: add soil and air variables
#     2021-Oct-19: sort variable to prognostic and dignostic catergories
#
#######################################################################################################################################################################################################
"""
$(TYPEDEF)

Struct that contains environmental conditions, such as soil moisture and atmospheric vapor pressure. Note that this structure is designed to be containers to interact with other CliMA modules and to
    prescribe values.

# Fields
$(TYPEDFIELDS)

# Examples
```julia
;
```
"""
mutable struct SoilAir{FT<:AbstractFloat}
    # parameters that do not change with time
    "Total area of the soil/air interface `[m²]`"
    AREA::FT
    "Soil color class used for soil albedo calculations"
    COLOR::Int
    "Number of air layers"
    N_AIR::Int
    "Number of soil layers"
    N_SOIL::Int
    "Soil moisture retention curve"
    VC::Union{Vector{BrooksCorey{FT}}, Vector{VanGenuchten{FT}}}
    "Z profile for air `[m]`"
    Z_AIR::Vector{FT}
    "Z profile for soil `[m]`"
    Z_SOIL::Vector{FT}
    "ΔZ profile for air `[m]`"
    ΔZ_AIR::Vector{FT}
    "ΔZ profile for soil `[m]`"
    ΔZ_SOIL::Vector{FT}

    # prognostic variables that change with time
    "CO₂ partial pressure at different air layers `[Pa]`"
    p_CO₂::Vector{FT}
    "H₂O partial pressure at different air layers `[Pa]`"
    p_H₂O::Vector{FT}
    "Temperature at different air layers `[K]`"
    t_air::Vector{FT}
    "Temperature at different soil layers `[K]`"
    t_soil::Vector{FT}
    "Wind speed (total) `[m s⁻¹]`"
    wind::Vector{FT}
    "Wind speed (vertical) `[m s⁻¹]`"
    wind_z::Vector{FT}
    "Soil water content at different soil layers `[m³ m⁻³]`"
    θ::Vector{FT}

    # dignostic variables that change with time
    "Saturated vapor pressure at different air layers `[Pa]`"
    p_H₂O_sat::Vector{FT}
    "Vapor pressure deficit at different air layers `[Pa]`"
    vpd::Vector{FT}
    "Soil water potential at different soil layers `[MPa]`"
    ψ::Vector{FT}

    # caches to speed up calculations

    # constructors
    SoilAir{FT}(z_soil::Vector{FT}, z_air::Vector{FT}, area::FT=FT(1)) where {FT<:AbstractFloat} = (
        _n_air  = length(z_air) - 1;
        _n_soil = length(z_soil) - 1;
        _p_sats = saturation_vapor_pressure(T_25()) * ones(_n_air);
        _vcs    = [VanGenuchten{FT}("Silt") for _i in 1:_n_soil];

        return new{FT}(area, 1, _n_air, _n_soil, _vcs, z_air, z_soil, diff(z_air), diff(z_soil), 41 * ones(_n_air),  0.5 * _p_sats, T_25() * ones(_n_air), T_25() * ones(_n_soil), 2 * ones(_n_air),
                       1 * ones(_n_air), [_vc.Θ_SAT for _vc in _vcs], _p_sats, 0.5 * _p_sats, zeros(_n_soil))
    );
end