text:
include(joinpath(@__DIR__, "config.jl"))
using JuMP, Gurobi, MLKernels

solver = with_optimizer(Gurobi.Optimizer; OutputFlag=0, Threads=1)

### data settings ###
data_sets = ["Annthyroid", "Arrhythmia", "Arrhythmia", "Cardiotocography", "Glass",
             "HeartDisease", "Hepatitis", "InternetAds", "Ionosphere", "PageBlocks",
             "Parkinson", "Pima", "SpamBase", "Stamps", "WPBC", "Wilt"]
data_dirs = [joinpath("dami", x) for x in data_sets]
data_info = "dami_outlier_detection"

### learning scenario ###
initial_pool_strategy = [("Pnin", Dict(:n => 25))]
split_strategy = [("Sl", Dict())]
initial_sample_strat = :nn_biased_sample

#### models ####
models = [Dict(:type => :SVDDneg, :param => Dict{Symbol, Any}())]
init_strategies = [WangCombinedInitializationStrategy(solver, 2.0.^range(-4, stop=4, step=1.0), FixedCStrategy(1.0))]
classify_precision = SVDD.OPT_PRECISION

#### oracle ####
oracle_param = [
    Dict{Symbol, Any}(
        :type => PoolOracle,
        :param => Dict{Symbol, Any}()
    )
]

#### query strategies ####
query_strategies = [
    Dict(:type => :RandomPQs, :param => Dict{Symbol, Any}()),
    Dict(:type => :RandomOutlierPQs, :param => Dict{Symbol, Any}()),
    Dict(:type => :HighConfidencePQs, :param => Dict{Symbol, Any}()),
    Dict(:type => :DecisionBoundaryPQs, :param => Dict{Symbol, Any}())
]

num_al_iterations = 100
num_resamples_initial_pool = 10
exp_name = "eval_part_2_pqs"

max_stars_repo_name: englhardt/des-evaluation

max_issues_repo_name: englhardt/des-evaluation

max_forks_repo_name: englhardt/des-evaluation


text:
using FFTW

mutable struct Point2D{T<:Real}
    x::T
    y::T
end

mutable struct Point3D{T<:Real}
    x::T
    y::T
    z::T
end

function gaussian(data, mean, sigma)
    return exp(-(sum((data - mean).^2)) / (2 * (sigma^2)))
end

"""
[h1, h2],[w1, w2]の総和
"""
function getRangeSum(
        accumArray::Array{Float64, 2},
        h1::Int,
        w1::Int,
        h2::Int,
        w2::Int)::Float64
    ret = accumArray[h2, w2]
    if (h1 > 1) ret -= accumArray[h1-1, w2] end
    if (w1 > 1) ret -= accumArray[h2, w1-1] end
    if (h1 > 1 && w1 > 1) ret += accumArray[h1-1, w1-1] end

    return ret
end

function bayesEstimationWithFFTConvolution(
        observed::Array{Float64, 2},
        calculated::Array{Float64, 2},
        sigma::Float64)::Float64

    imgH, imgW = size(observed)
    N_pix = imgH * imgW
    plausibleSum = Float64(0)

    C_oo = sum(observed .^ 2)
    C_ccAcc = [calculated calculated; calculated calculated] .^ 2

    for h in 1:2imgH, w in 1:2imgW
        if w != 2imgW C_ccAcc[h, w+1] += C_ccAcc[h, w] end
    end
    for h in 1:2imgH, w in 1:2imgW
        if h != 2imgH C_ccAcc[h+1, w] += C_ccAcc[h, w] end
    end

    C_coMat = real.(ifft(fft(calculated).*conj.(fft(observed))))

    for h in 1:imgH, w in 1:imgW
        C_cc = getRangeSum(C_ccAcc, h, w, h+imgH-1, w+imgW-1)
        C_co = C_coMat[h, w]

        plausibleSum += exp(- (C_oo - 2 * C_co + C_cc) / (2 * sigma^2) )
    end

    return plausibleSum
end

"""
function bayesEstimationWithFFTConvolution(
        observed::Array{Float64, 2},
        calculated::Array{Float64, 2},
        imgH::Int,
        imgW::Int,
        extendImgH::Int,
        extendImgW::Int,
        sigma::Float64)::Float64

    N_pix = imgH * imgW
    plausibleSum = Float64(0)

    C_o = sum(observed)
    C_oo = sum(observed .^ 2)

    C_cAcc = [calculated calculated; calculated calculated]
    C_ccAcc = C_cAcc .^ 2

    for h in 1:2extendImgH, w in 1:2extendImgW
        if w != 2extendImgW
            C_cAcc[h, w+1] += C_cAcc[h, w]
            C_ccAcc[h, w+1] += C_ccAcc[h, w]
        end
    end
    for h in 1:2extendImgH, w in 1:2extendImgW
        if h != 2extendImgH
            C_cAcc[h+1, w] += C_cAcc[h, w]
            C_ccAcc[h+1, w] += C_ccAcc[h, w]
        end
    end

    C_coLis = abs.(ifft(fft(calculated).*conj.(fft(observed))))

    basePlausible = sigma^(2-N_pix) * (2pi)^(1-N_pix/2)

    for h in 1:extendImgH, w in 1:extendImgW
        C_c = getRangeSum(C_cAcc, h, w, h+imgH-1, w+imgW-1)
        C_cc = getRangeSum(C_ccAcc, h, w, h+imgH-1, w+imgW-1)
        C_co = C_coLis[h, w]

        if N_pix * C_cc - C_c^2 < 1e-8 continue end
        if C_co < 1e-8 continue end

        thisPlausible = basePlausible
        index = N_pix * (C_cc * C_oo - C_co^2)
        index += 2C_o * C_co * C_c
        index -= C_cc * C_o^2
        index -= C_oo * C_c^2
        index /= 2 * sigma^2 * (N_pix * C_cc - C_c^2)

        if index < 0 continue end

        thisPlausible *= exp(-index)
        thisPlausible /= sqrt(N_pix * C_cc - C_c^2)

        plausibleSum += thisPlausible
    end

    return plausibleSum
end
"""

ext: jl

max_stars_repo_name: matsunagalab/BioAFM

max_issues_repo_name: matsunagalab/BioAFM

max_forks_repo_name: matsunagalab/BioAFM


text:
using jInv.Utils
using Test

# test checkDerivative
f = x-> x^2
dfRight = x -> 2*x
p,e,o = checkDerivative(f,dfRight,randn(),out=true)
@test p==true

F = x->x.^2
dFRight = x-> Diagonal(2*x)

x0 = randn(10).+2
p,e,o = checkDerivative(F,dFRight,x0,out=true)
@test p==true

dfWrong = x-> x
dFWrong = x->randn(length(x),length(x))
p,e,o = checkDerivative(f,dfWrong,randn(),out=false)
@test p==false

x0 = randn(10)
p,e,o = checkDerivative(F,dFWrong,x0,out=false)
@test p==false


text:
mutable struct Point
	x::Float64
	y::Float64
end

Base.show(io::IO, p::Point) = print(io, "Point($(p.x), $(p.y))")

getx(p::Point) = p.x
gety(p::Point) = p.y

setx(p::Point, x) = (p.x = x)
sety(p::Point, y) = (p.y = y)

mutable struct Circle
	x::Float64
	y::Float64
	r::Float64
end

getx(c::Circle) = c.x
gety(c::Circle) = c.y
getr(c::Circle) = c.r

setx(c::Circle, x) = (c.x = x)
sety(c::Circle, y) = (c.y = y)
setr(c::Circle, r) = (c.r = r)

Base.show(io::IO, c::Circle) = print(io, "Circle($(c.x), $(c.y), $(c.r))")

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# This file is a part of ValueShapes.jl, licensed under the MIT License (MIT).

using ValueShapes
using Test

using Random
using ArraysOfArrays
import TypedTables

@testset "scalar_shape" begin
    @test @inferred(ValueShapes._valshapeoftype(Int)) == ScalarShape{Int}()
    @test @inferred(ValueShapes._valshapeoftype(Complex{Float64})) == ScalarShape{Complex{Float64}}()

    @test @inferred(valshape(3)) == ScalarShape{Int}()

    @test @inferred(size(ScalarShape{Real}())) == ()
    @test @inferred(size(ScalarShape{Complex}())) == ()

    @test @inferred(size(ScalarShape{Real}())) == ()
    @test @inferred(size(ScalarShape{Complex}())) == ()

    @test @inferred(ValueShapes.default_unshaped_eltype(ScalarShape{Real}())) == Float64
    @test @inferred(ValueShapes.default_unshaped_eltype(ScalarShape{Complex{Real}}())) == Float64
    @test @inferred(ValueShapes.default_unshaped_eltype(ScalarShape{Complex{Int32}}())) == Int32

    @test @inferred(ValueShapes.shaped_type(ScalarShape{Real}())) == Float64
    @test @inferred(ValueShapes.shaped_type(ScalarShape{Real}(), Float32)) == Float32
    @test @inferred(ValueShapes.shaped_type(ScalarShape{Complex}())) == Complex{Float64}
    @test @inferred(ValueShapes.shaped_type(ScalarShape{Complex{Real}}(), Float32)) == Complex{Float32}
    @test @inferred(ValueShapes.shaped_type(ScalarShape{Complex{Int16}}())) == Complex{Int16}

    @test @inferred(totalndof(ScalarShape{Int}())) == 1
    @test @inferred(totalndof(ScalarShape{Complex{Float64}}())) == 2

    @test @inferred(ScalarShape{Real}()(undef)) === zero(Float64)
    @test @inferred(ScalarShape{Complex}()(undef)) === zero(Complex{Float64})

    @test typeof(@inferred(ScalarShape{Real}()([42]))) <: SubArray{Int,0}
    @test @inferred(ScalarShape{Real}()([42])[]) == 42

    @test @inferred (ScalarShape{Int}() <= ScalarShape{Real}()) == true
    @test @inferred (ScalarShape{Real}() <= ScalarShape{Int}()) == false
    @test @inferred (ScalarShape{Real}() >= ScalarShape{Int}()) == true

    let shape = ScalarShape{Real}(), data = [4.2]
        @test @inferred(unshaped(shape(data), shape)) == data
        @test @inferred(unshaped(shape(data)[], shape)) == data
        @test_throws ArgumentError unshaped(shape([3, 4]), shape)
        @test_throws ArgumentError unshaped(shape(data), ScalarShape{Integer}())
    end

    @test let
        shape = ScalarShape{Real}()
        valshape(shape.(push!(@inferred(VectorOfSimilarVectors{Float64}(shape)), @inferred(Vector{Float64}(undef, shape))))[1]) == valshape(shape(undef))
    end

    let
        A = collect(1:6)
        ac = ValueAccessor(ScalarShape{Real}(), 2)
        @test @inferred(getindex(A, ac)) == 3
        @test @inferred(view(A, ac))[] == 3
        @test view(A, ac) isa AbstractArray{<:Real,0}
        @test @inferred(setindex!(A, 7, ac)) === A
        @test A[ac] == 7
    end

    let
        A = collect(reshape(1:6*6, 6, 6))
        ac1 = ValueAccessor(ScalarShape{Real}(), 2)
        ac2 = ValueAccessor(ScalarShape{Real}(), 4)
        @test @inferred(getindex(A, ac1, ac2)) == 27
        @test @inferred(view(A, ac1, ac2))[] == 27
        @test view(A, ac1, ac2) isa AbstractArray{<:Real,0}
        @test @inferred(setindex!(A, 42, ac1, ac2)) === A
        @test A[ac1, ac2] == 42
    end
end


text:
@doc """
    improved_estimator(model::QuantumXXZ, T::Real, Js::AbstractArray, uf::UnionFind)

Returns the following observables as `Dict{String, Any}` using loop information `uf`

# Observables
- `"Sign"`
    - Sign of the weight function
- `"Sign * Energy"`
    - Energy per spin (site)
- `"Sign * Energy^2"`
- `"Sign * Magnetization"`
    - Total magnetization (Sz) per spin (site)
- `"Sign * |Magnetization|"`
- `"Sign * Magnetization^2"`
- `"Sign * Magnetization^4"`
"""
function improved_estimator(model::QuantumXXZ, T::Real, Jzs::AbstractArray, Jxys::AbstractArray, Gs::AbstractArray, uf::UnionFind)
    S2 = model.S2
    nsites = numsites(model)
    nspins = nsites*S2
    nbonds = numbonds(model)
    nc = numclusters(uf)

    ms = zeros(nc)
    for s in 1:nspins
        i = clusterid(uf, s)
        ms[i] += model.spins[s]
    end
    ms .*= 0.5/nsites

    E0 = 0.0
    for st in 1:numsitetypes(model)
        E0 += 0.5 * numsites(model, st) * Gs[st] * S2
    end
    for bt in 1:numbondtypes(model)
        nb = numbonds(model,bt)*S2*S2
        z = nb*Jzs[bt]
        x = nb*abs(Jxys[bt])
        if z > x
            ## AntiFerroIsing like
            E0 += 0.25z
        elseif z < -x
            ## FerroIsing like
            E0 -= 0.25z
        else
            ## XY like
            E0 += 0.25x
        end
    end
    nops = length(model.ops)
    E = E0-nops*T
    E2 = nops*(nops-1)*T^2 - 2*E0*T*nops + E0^2
    E /= nsites
    E2 /= nsites^2

    M = 0.0
    M2 = 0.0
    M4 = 0.0
    for m in ms
        M += m
        m2 = m*m
        M4 += m2*m2 + 6M2*m2
        M2 += m2
    end

    sgn = 1.0
    for op in model.ops
        if !op.isdiagonal
            if op.let_type == LET_Vertex || op.let_type == LET_Cross
                b = op.space-nsites
                bt = bondtype(model, b)
                sgn *= ifelse(Jxys[bt] > 0.0, -1.0, 1.0)
            end
        end
    end

    res = Measurement()
    res["Sign * Magnetization"] = sgn*M
    res["Sign * |Magnetization|"] = sgn*abs(M)
    res["Sign * Magnetization^2"] = sgn*M2
    res["Sign * Magnetization^4"] = sgn*M4
    res["Sign * Energy"] = sgn*E
    res["Sign * Energy^2"] = sgn*E2
    res["Sign"] = sgn

    return res
end

default_estimator(model::QuantumXXZ, update) = improved_estimator


text:
function direct_deps!(reach::Vector{Variable}, p, var, ir)
    for (v, st) in ir
        st.expr isa Expr && begin
            if var in st.expr.args
                push!(reach, v)
                if st.expr.head == :call &&
                    st.expr.args[1] isa GlobalRef && st.expr.args[1].name == :trace
                    # blocked
                else
                    direct_deps!(reach, p, v, ir)
                end
            end
        end
    end
end
function direct_deps(var::Variable, ir)
    reach = Variable[]
    for (v, st) in ir
        st.expr isa Expr && begin
            if var in st.expr.args
                push!(reach, v)
                direct_deps!(reach, var, v, ir)
            end
        end
    end
    return reach
end
function direct_dep_analysis(ir)
    sites = Variable[]
    addrs = QuoteNode[]
    var_sym_map = Dict{Variable, QuoteNode}()
    reach = Dict{Variable, Any}()
    for (v, st) in ir
        Jaynes.MacroTools.postwalk(st) do e
            Jaynes.@capture(e, call_(sym_, args__))
            call isa GlobalRef && call.name == :trace && begin
                push!(sites, v)
                push!(addrs, sym)
                var_sym_map[v] = sym
                r = direct_deps(v, ir)
                !isempty(r) && begin
                    reach[v] = r
                end
            end
            e
        end
    end
    return StaticAnalysis(reach, sites, addrs, var_sym_map)
end

function markov_blanket(ir::IRTools.IR)
    ds = dependency(direct_dep_analysis(ir))

    addrs = ds.addresses
    G = copy(ds.dependencies)

    # add empty diagonal
    for p in addrs
        if !haskey(G, p)
            G[p] = Set{Address}()
        end
    end

    # moralize
    parents = [[i for (i, js) in pairs(G) if k in js] for k in addrs]
    for pars in parents
        for p in pars
            for q in pars
                p == q && continue
                push!(G[p], q)
            end
        end
    end

    # symmetrize and add diagonal
    for (p, qs) in G
        for q in qs
            p == q && continue
            push!(G[q], p)
        end
        push!(G[p], p)
    end
    G
end

function markov_blanket(call::Function, type...)
    ir = lower_to_ir(call, type...)
    markov_blanket(ir)
end


text:
function construct_device!(
    psi_container::PSIContainer,
    sys::PSY.System,
    model::DeviceModel{St, D},
    ::Type{S},
) where {St <: PSY.Storage, D <: AbstractStorageFormulation, S <: PM.AbstractPowerModel}
    devices = get_available_components(St, sys)

    if !validate_available_devices(St, devices)
        return
    end

    #Variables
    active_power_variables!(psi_container, devices)
    reactive_power_variables!(psi_container, devices)
    energy_storage_variables!(psi_container, devices)

    #Initial Conditions
    initial_conditions!(psi_container, devices, D)

    #Constraints
    active_power_constraints!(psi_container, devices, model, S, get_feedforward(model))
    reactive_power_constraints!(psi_container, devices, model, S, get_feedforward(model))
    energy_capacity_constraints!(psi_container, devices, model, S, get_feedforward(model))
    feedforward!(psi_container, St, get_feedforward(model))

    # Energy Balanace limits
    energy_balance_constraint!(psi_container, devices, D, S, get_feedforward(model))

    return
end

function construct_device!(
    psi_container::PSIContainer,
    sys::PSY.System,
    model::DeviceModel{St, D},
    ::Type{S},
) where {
    St <: PSY.Storage,
    D <: AbstractStorageFormulation,
    S <: PM.AbstractActivePowerModel,
}

    devices = get_available_components(St, sys)

    if !validate_available_devices(St, devices)
        return
    end

    #Variables
    active_power_variables!(psi_container, devices)
    energy_storage_variables!(psi_container, devices)

    #Initial Conditions
    initial_conditions!(psi_container, devices, D)

    #Constraints
    active_power_constraints!(psi_container, devices, model, S, get_feedforward(model))
    energy_capacity_constraints!(psi_container, devices, model, S, get_feedforward(model))
    feedforward!(psi_container, St, get_feedforward(model))

    # Energy Balanace limits
    energy_balance_constraint!(psi_container, devices, D, S, get_feedforward(model))

    return
end

function construct_device!(
    psi_container::PSIContainer,
    sys::PSY.System,
    model::DeviceModel{St, BookKeepingwReservation},
    ::Type{S},
) where {St <: PSY.Storage, S <: PM.AbstractPowerModel}
    devices = get_available_components(St, sys)

    if !validate_available_devices(St, devices)
        return
    end

    #Variables
    active_power_variables!(psi_container, devices)
    reactive_power_variables!(psi_container, devices)
    energy_storage_variables!(psi_container, devices)
    storage_reservation_variables!(psi_container, devices)

    #Initial Conditions
    initial_conditions!(psi_container, devices, model.formulation)

    #Constraints
    active_power_constraints!(psi_container, devices, model, S, get_feedforward(model))
    reactive_power_constraints!(psi_container, devices, model, S, get_feedforward(model))
    energy_capacity_constraints!(psi_container, devices, model, S, get_feedforward(model))
    feedforward!(psi_container, St, get_feedforward(model))

    # Energy Balanace limits
    energy_balance_constraint!(
        psi_container,
        devices,
        model.formulation,
        S,
        get_feedforward(model),
    )

    return
end

function construct_device!(
    psi_container::PSIContainer,
    sys::PSY.System,
    model::DeviceModel{St, BookKeepingwReservation},
    ::Type{S},
) where {St <: PSY.Storage, S <: PM.AbstractActivePowerModel}
    devices = get_available_components(St, sys)

    if !validate_available_devices(St, devices)
        return
    end

    #Variables
    active_power_variables!(psi_container, devices)
    energy_storage_variables!(psi_container, devices)
    storage_reservation_variables!(psi_container, devices)

    #Initial Conditions
    initial_conditions!(psi_container, devices, model.formulation)

    #Constraints
    active_power_constraints!(psi_container, devices, model, S, get_feedforward(model))
    energy_capacity_constraints!(psi_container, devices, model, S, get_feedforward(model))
    feedforward!(psi_container, St, get_feedforward(model))

    # Energy Balanace limits
    energy_balance_constraint!(
        psi_container,
        devices,
        model.formulation,
        S,
        get_feedforward(model),
    )

    return
end


text:
struct Start <: EdgeModifier
    s
end

@testset "Unmeta Node" begin
    @test unmeta(Node(:a)) == Node(:a)
    @test unmeta(Node(:a) * Start(4)) == Node(:a)
end

@testset "Unmeta Edge" begin
    @test unmeta(Edge(Node(:a), Node(:b))) == Edge(Node(:a), Node(:b))
    @test unmeta(Edge(Node(:a), Node(:b)) * Start(4)) == Edge(Node(:a), Node(:b))
end


text:
using Test
# Test constraint stuff
n,m = 3,2
cE(x,u) = [2x[1:2]+u;
          x'x + 5]
pE = 3
cE(x) = [cos(x[1]) + x[2]*x[3]; x[1]*x[2]^2]
pE_N = 2
cI(x,u) = [x[3]-x[2]; u[1]*x[1]]
pI = 2
pI_N = 0

model, obj = Dynamics.dubinscar
obj.tf = 3
count_inplace_output(cI,n,m)
obj_con = ConstrainedObjective(obj,cE=cE,cI=cI)
@test obj_con.p == pE + pI
@test obj_con.pI == pI
@test obj_con.pI_N == pI_N
@test obj_con.p_N == n + pE_N + pI_N

N = 5
solver = Solver(model, obj, N=N)
@test solver.state.constrained == true
@test get_num_constraints(solver) == (5,2,3)
@test original_constraint_inds(solver) == trues(5)
@test get_constraint_labels(solver) == ["custom inequality", "custom inequality", "custom equality", "custom equality", "custom equality"]

# Add state and control bounds
obj_con = update_objective(obj_con, u_min=[-10,-Inf], u_max=10, x_min=[-Inf,-10,-Inf], x_max=[10,12,10])
pI_bnd = 1 + m + 1 + n
@test obj_con.pI == pI + pI_bnd
@test obj_con.p == pI + pI_bnd + pE
p = obj_con.p
pI = obj_con.pI
pE = p-pI

N = 5
solver = Solver(model, obj_con, N=N)
@test solver.state.constrained == true
@test get_num_constraints(solver) == (5+pI_bnd,2+pI_bnd,3)
@test original_constraint_inds(solver) == trues(5+pI_bnd)
@test get_constraint_labels(solver) == ["control (upper bound)", "control (upper bound)", "control (lower bound)", "state (upper bound)", "state (upper bound)", "state (upper bound)", "state (lower bound)",
    "custom inequality", "custom inequality", "custom equality", "custom equality", "custom equality"]

# Infeasible controls
solver_inf = Solver(model, obj_con, N=N)
solver_inf.opts.infeasible = true
@test solver_inf.opts.constrained == true
@test get_num_constraints(solver_inf) == (p+n,pI,pE+n)
@test original_constraint_inds(solver_inf) == [trues(p); falses(n)]
@test get_constraint_labels(solver_inf) == ["control (upper bound)", "control (upper bound)", "control (lower bound)", "state (upper bound)", "state (upper bound)", "state (upper bound)", "state (lower bound)",
    "custom inequality", "custom inequality", "custom equality", "custom equality", "custom equality",
    "* infeasible control","* infeasible control","* infeasible control"]

# Minimum time
obj_mintime = update_objective(obj_con, tf=:min)
solver_min = Solver(model, obj_mintime, N=N)
@test solver_min.opts.constrained == true
@test get_num_constraints(solver_min) == (p+3,pI+2,pE+1)
@test original_constraint_inds(solver_min) == [true; true; false; true; false; trues(4); trues(5); false]
@test get_constraint_labels(solver_min) == ["control (upper bound)", "control (upper bound)", "* √dt (upper bound)", "control (lower bound)", "* √dt (lower bound)", "state (upper bound)", "state (upper bound)", "state (upper bound)", "state (lower bound)",
    "custom inequality", "custom inequality", "custom equality", "custom equality", "custom equality",
    "* √dt (equality)"]

# Minimum time and infeasible
obj_mintime = update_objective(obj_con, tf=:min)
solver_min = Solver(model, obj_mintime, N=N)
solver_min.opts.infeasible = true
@test solver_min.opts.constrained == true
@test get_num_constraints(solver_min) == (p+3+n,pI+2,pE+1+n)
@test original_constraint_inds(solver_min) == [true; true; false; true; false; trues(4); trues(5); falses(4)]
@test get_constraint_labels(solver_min) == ["control (upper bound)", "control (upper bound)", "* √dt (upper bound)", "control (lower bound)", "* √dt (lower bound)", "state (upper bound)", "state (upper bound)", "state (upper bound)", "state (lower bound)",
    "custom inequality", "custom inequality", "custom equality", "custom equality", "custom equality",
    "* infeasible control","* infeasible control","* infeasible control","* √dt (equality)"]
get_constraint_labels(solver_min)


text:
@doc """ Driver for compare_walks_parallel.jl.
This file sets the parameters to be run.
A command line for running this program:
julia -p 8 -L "compare_walks_driver.jl" -e "run_all_jobs()"
"""
include("../src/NKLandscapes.jl")
using NKLandscapes
using Base.Test
include("compare_walks_parallel.jl")

const BOUNDARIES = ["Nearest-neighbor", "Random"]

#const N_VALUES = [60, 100, 140, 180]
const N_VALUES = [220, 240]
const K_VALUE =  8
const Q_VALUE = 0  # q for NKq landscape, set to 0 for NK landscape
const B_VALUE = 2
const RUNS = 100  # number of runs for each setting of the parameters

function run_all_jobs()
  print_header()
  trial = 1
  for n_value in N_VALUES
    jobs = [Job(trial, n_value, K_VALUE, Q_VALUE, B_VALUE) for trial = 1:RUNS ]
    run_jobs(jobs)
    trial += 1
  end
end


text:
module AutomotiveDrivingModels

using Printf
using LinearAlgebra
using Parameters
using StaticArrays
using Distributions
using Reexport
using Random
@reexport using Vec
@reexport using Records

# Roadways

export StraightRoadway,
       mod_position_to_roadway,
       get_headway

include("roadways/straight_1d_roadways.jl")

export CurvePt,
       Curve,
       CurveIndex,
       CurveProjection,
       CURVEINDEX_START,
       get_lerp_time,
       index_closest_to_point,
       get_curve_index,
       curveindex_end

include("roadways/curves.jl")

export
    LaneTag,
    LaneBoundary,
    Lane,
    LaneConnection,
    SpeedLimit,
    RoadSegment,
    Roadway,
    RoadProjection,
    RoadIndex,

    NULL_BOUNDARY,
    NULL_LANETAG,
    NULL_ROADINDEX,

    DEFAULT_SPEED_LIMIT,
    DEFAULT_LANE_WIDTH,

    is_in_exits,
    is_in_entrances,
    is_at_curve_end,
    is_between_segments_lo,
    is_between_segments_hi,
    is_between_segments,
    has_segment,
    has_lanetag,
    has_next,
    has_prev,
    next_lane,
    prev_lane,
    next_lane_point,
    prev_lane_point,
    connect!,
    move_along,
    n_lanes_left,
    n_lanes_right,
    lanes,
    lanetags,
    read

include("roadways/roadways.jl")

export
    Frenet,
    NULL_FRENET

include("roadways/frenet.jl")

export
        gen_straight_curve,
        gen_straight_segment,
        gen_straight_roadway,
        gen_stadium_roadway,
        gen_bezier_curve

include("roadways/roadway_generation.jl")

## agent definitions

export
    AbstractAgentDefinition,
    AgentClass,
    VehicleDef,
    BicycleModel,
    length,
    width,
    class

include("agent-definitions/agent_definitions.jl")

## states

export
    posf,
    posg,
    vel,
    velf,
    velg,
    VehicleState,
    Vehicle,
    get_vel_s,
    get_vel_t,
    get_center,
    get_footpoint,
    get_front,
    get_rear,
    get_lane

include("states/interface.jl")
include("states/vehicle_state.jl")

export
    TrajdataFrame,
    TrajdataState,
    Trajdata,
    TrajdataVehicleIterator

include("states/trajdatas.jl")

export
    Scene,
    SceneRecord

include("states/scenes.jl")

## Collision Checkers

export
    ConvexPolygon,
    CPAMemory,
    CollisionCheckResult,
    to_oriented_bounding_box!,
    get_oriented_bounding_box,
    is_colliding,
    is_potentially_colliding,
    get_collision_time,
    get_first_collision,
    is_collision_free,
    get_distance,
    get_edge,
    collision_checker,
    polygon

include("collision-checkers/minkowski.jl")
include("collision-checkers/parallel_axis.jl")

## Feature Extraction

export
    VehicleTargetPoint,
    VehicleTargetPointFront,
    VehicleTargetPointCenter,
    VehicleTargetPointRear,
    get_targetpoint_delta,
    NeighborLongitudinalResult,
    get_neighbor_fore,
    get_neighbor_rear,
    get_headway,
    get_neighbor_fore_along_lane,
    get_neighbor_fore_along_left_lane,
    get_neighbor_fore_along_right_lane,
    get_neighbor_rear_along_lane,
    get_neighbor_rear_along_left_lane,
    get_neighbor_rear_along_right_lane,
    FrenetRelativePosition,
    get_frenet_relative_position,
    get_lane_width,
    get_markerdist_left,
    get_markerdist_right

include("feature-extraction/neighbors_features.jl")
include("feature-extraction/lane_features.jl")

export
    AbstractFeature,
    FeatureValue,
    FeatureState,
    is_feature_valid,
    is_symbol_a_feature,
    allfeatures,
    symbol2feature,
    AbstractFeatureExtractor,
    FeatureExtractor,
    SubsetExtractor,
    StandardizingExtractor,
    pull_features!,
    rec_length

include("feature-extraction/interface.jl")
include("feature-extraction/features.jl")
include("feature-extraction/features_extractors.jl")

export
    LidarSensor,
    nbeams,
    observe!,
    RoadlineLidarSensor,
    nlanes,
    LanePortion,
    RoadwayLidarCulling,
    ensure_leaf_in_rlc!,
    get_lane_portions

include("feature-extraction/lidar_sensor.jl")

## Actions

export
    propagate,
    LaneFollowingAccel,
    AccelTurnrate,
    AccelDesang,
    LatLonAccel,
    AccelSteeringAngle,
    PedestrianLatLonAccel

include("actions/interface.jl")
include("actions/lane_following_accel.jl")
include("actions/accel_turnrate.jl")
include("actions/accel_desang.jl")
include("actions/lat_lon_accel.jl")
include("actions/accel_steering.jl")
include("actions/pedestrian_lat_lon_accel.jl")

export
    DriverModel,
    StaticDriver,
    get_name,
    action_type,
    set_desired_speed!,
    observe!,
    reset_hidden_state!,
    prime_with_history!

include("behaviors/interface.jl")

export
    LaneFollowingDriver,
    StaticLaneFollowingDriver,
    PrincetonDriver,
    IntelligentDriverModel,
    ProportionalSpeedTracker,
    track_longitudinal!,
    LateralDriverModel,
    ProportionalLaneTracker,
    LatLonSeparableDriver,
    Tim2DDriver,
    track_lane!,
    SidewalkPedestrianModel,
    LaneChangeChoice,
    LaneChangeModel,
    get_lane_offset,
    DIR_RIGHT,
    DIR_MIDDLE,
    DIR_LEFT,
    MOBIL,
    TimLaneChanger

include("behaviors/lane_following_drivers.jl")
include("behaviors/princeton_driver.jl")
include("behaviors/speed_trackers.jl")
include("behaviors/intelligent_driver_model.jl")
include("behaviors/lateral_driver_models.jl")
include("behaviors/lane_change_models.jl")
include("behaviors/MOBIL.jl")
include("behaviors/tim_lane_changer.jl")
include("behaviors/lat_lon_separable_driver.jl")
include("behaviors/tim_2d_driver.jl")
include("behaviors/sidewalk_pedestrian_model.jl")

export
    get_actions!,
    tick!,
    reset_hidden_states!,
    simulate,
    simulate!,
    EntityAction,
    run_callback,
    CollisionCallback

include("simulation/simulation.jl")
include("simulation/callbacks.jl")

export
    State1D,
    Vehicle1D,
    Scene1D

include("deprecated.jl")

end # AutomotiveDrivingModels


text:
using TullioStatistics
using Test
using Statistics
using StableRNGs

rng = StableRNG(111)

A = rand(rng,100);

@test TullioStatistics.mean(A) ≈ mean(A)


text:
## 3D cases
#==========================================================================================#

"""
    solve_case(components :: Vector{Horseshoe}, fs :: Freestream, refs :: References;
               name = :aircraft :: Symbol,
               print = false :: Boolean,
               print_components = false :: Boolean)

Perform a vortex lattice analysis given a vector of `Horseshoe`s, a `Freestream` condition, and `Reference` values.
"""
function solve_case(components, freestream :: Freestream, refs :: References; name = :aircraft, print = false, print_components = false)
    system = solve_system(components, freestream, refs)

    # Printing if needed
    if print_components
        print_coefficients(system, name, components = true)
    elseif print
        print_coefficients(system, name)
    end

    system
end

## Placeholder for functions I'm not sure where to put
#==========================================================================================#

spanwise_loading(CXs, weighted_areas) = vec(sum(CXs, dims = 1)) .* weighted_areas

## Span-loading
"""
    spanwise_loading(panels, CFs, S)
    spanwise_loading(mesh :: WingMesh, CFs, S)

Compute the spanwise loading of the forces given panels, associated force coefficients, and the reference area.

Alternatively, provide a `WingMesh` type which gets the panels for you.
"""
function spanwise_loading(panels, CFs, S)
    # Get y-coordinates of spanwise strips
    ys = @views vec(mean(x -> midpoint(x)[2], panels, dims = 1))

    # Compute weighted areas for spanwise strips
    area_scale = S ./ vec(sum(panel_area, panels, dims = 1))
    # chords     = vec(sum(average_chord, panels, dims = 1))

    # Compute spanwise coefficients
    span_CFs = permutedims(combinedimsview(spanwise_loading(CFs, area_scale)))

    [ ys span_CFs ]
end

spanwise_loading(wing :: WingMesh, CFs, S) = spanwise_loading(chord_panels(wing), CFs, S)

## Mesh connectivities
triangle_connectivities(inds) = @views [
                                         vec(inds[1:end-1,1:end-1]) vec(inds[1:end-1,2:end]) vec(inds[2:end,2:end])    ;
                                         vec(inds[2:end,2:end])     vec(inds[2:end,1:end-1]) vec(inds[1:end-1,1:end-1])
                                       ]

## Extrapolating surface values to neighbouring points
function extrapolate_point_mesh(mesh, weight = 0.75)
    m, n   = size(mesh)
    points = zeros(eltype(mesh), m + 1, n + 1)

    # The quantities are measured at the bound leg of a Horseshoe by default (0.25×)
    @views @. points[1:end-1,1:end-1] += weight       / 2 * mesh
    @views @. points[1:end-1,2:end]   += weight       / 2 * mesh
    @views @. points[2:end,1:end-1]   += (1 - weight) / 2 * mesh
    @views @. points[2:end,2:end]     += (1 - weight) / 2 * mesh

    points
end


text:
using Unitful
using Unitful: uconvert, NoUnits, cm, g, eV, MeV, u
using PeriodicTable: Element, elements

function datapath(args...)
    datadir = joinpath(@__DIR__, "..", "data")
    abspath(joinpath(datadir, args...))
end

_ucoalesce(unit, q::Quantity) = uconvert(unit, q)
_ucoalesce(unit, x::Number  ) = x*unit

_ustrip(x::Quantity) = uconvert(UF.NoUnits, x)
_ustrip(x) = x

function density(mat)
    el = aselement(mat)
    _ucoalesce(g*cm^-3, el.density)
end
aselement(el::Element) = el
aselement(key) = elements[key]

function linterpol(x, (x_lo, y_lo), (x_hi, y_hi))
    @assert x_lo <= x <= x_hi
    if x_lo == x_hi
        @assert y_lo == y_hi
        return y_lo
    end
    w = x_hi - x_lo
    w_hi = _ustrip((x - x_lo) / w)
    w_lo = _ustrip((x_hi - x) / w)
    @assert w_lo + w_hi ≈ 1
    y_lo * w_lo + y_hi * w_hi
end

lower_material(s::DataSource, Z::Int)       = Z
lower_material(s::DataSource, mat) = aselement(mat).number

lower_particle(s::DataSource, pt::Particle) = pt.energy
lower_particle(s::DataSource, E::UF.Energy) = E

function lookup(s::DataSource,
                           mat, pt, p)
    E = lower_particle(s, pt)
    _lookup(s, mat, E, p)
end

function get_table(s::DataSource, key)
    get_element_table(s, key)
end
function get_element_table(s::DataSource, key)
    Z = lower_material(s, key)
    s.element_tables[Z]
end

function _lookup(s::DataSource,mat,E::UF.Energy, p::Process)
    table = get_table(s, mat)
    Es = table.E
    Cs = getcol(table, p)
    E_min = first(Es)
    E_max = last(Es)
    @assert E_min <= E <= E_max
    E == E_min && return first(Cs)
    index_hi = searchsortedfirst(Es, E)
    index_lo = index_hi - 1
    E_lo = Es[index_lo]
    E_hi = Es[index_hi]
    # check for absorption edge
    absorption_edge = (E == E_hi) &&
        (get(Es, index_hi+1, nothing) == E_hi)
    if absorption_edge
        msg = "Material $mat has an absorption edge at E=$(E)."
        throw(ArgumentError(msg))
    end
    c_lo = Cs[index_lo]
    c_hi = Cs[index_hi]
    linterpol(E, E_lo => c_lo, E_hi => c_hi)
end

function tabletype(pairs...)
    symbols = Symbol[]
    types   = Type[]
    for (s, u) in pairs
        push!(symbols, s)
        Q = typeof(1.0*u)
        @assert (Q <: Quantity) || Q == Float64
        push!(types, Vector{Q})
    end
    NamedTuple{tuple(symbols...), Tuple{types...}}
end

function columnnames(::Type{S}) where {S<:DataSource}
    NT = eltype(first(S.types))
    fieldnames(NT)
end

function get_process_symbol(p)::Symbol
    Symbol(last(split(string(p), ".")))
end
@generated function getcol(table, process::P) where {P}
    s = get_process_symbol(process)
    :(table.$s)
end

empty_col(u) = typeof(1.0*u)[]

function load(::Type{S}; Zs, compounds=nothing, dir)::S where {S<:DataSource}
    header = map(string, columnnames(S))
    element_tables = asyncmap(Zs, ntasks=length(Zs)) do Z
        path = datapath(dir,"Z$Z.csv")
        out = emptytable(S)
        readcsv!(path, out, header=header)
    end
    if compounds === nothing
        S(element_tables)
    else
        compound_pairs = asyncmap(compounds) do compound
            path = datapath(dir, "$compound.csv")
            compound => readcsv!(path, emptytable(S), header=header)
        end
        compound_tables = Dict(compound_pairs...)
        S(element_tables, compound_tables)
    end
end


text:
using Kinetic, Plots, LinearAlgebra, JLD2
using ProgressMeter: @showprogress

cd(@__DIR__)
D = Dict{Symbol,Any}()
begin
    D[:matter] = "gas"
    D[:case] = "cavity"
    D[:space] = "2d1f2v"
    D[:flux] = "kfvs"
    D[:collision] = "bgk"
    D[:nSpecies] = 1
    D[:interpOrder] = 2
    D[:limiter] = "vanleer"
    D[:boundary] = "maxwell"
    D[:cfl] = 0.5
    D[:maxTime] = 1.0

    D[:x0] = -1.5
    D[:x1] = 1.5
    D[:nx] = 40
    D[:y0] = -1.5
    D[:y1] = 1.5
    D[:ny] = 40
    D[:pMeshType] = "uniform"
    D[:nxg] = 0
    D[:nyg] = 0

    D[:umin] = -5.0
    D[:umax] = 5.0
    D[:nu] = 24
    D[:vmin] = -5.0
    D[:vmax] = 5.0
    D[:nv] = 24
    D[:vMeshType] = "rectangle"
    D[:nug] = 0
    D[:nvg] = 0

    D[:knudsen] = 0.5
    D[:mach] = 0.0
    D[:prandtl] = 1.0
    D[:inK] = 0.0
    D[:omega] = 0.81
    D[:alphaRef] = 1.0
    D[:omegaRef] = 0.5

    D[:uLid] = 0.15
    D[:vLid] = 0.0
    D[:tLid] = 1.0
end

ks = SolverSet(D)
ctr, a1face, a2face = init_fvm(ks, ks.ps, :dynamic_array; structarray = true)

s2 = 0.03^2
s1 = 0.1#4π * s2
flr = 0.0001
init_density(x, y) = max(flr, 1.0 / s1 * exp(-(x^2 + y^2) / 4.0 / s2))

function init_field!(ks, ctr, a1face, a2face)
    for i = 1:ks.ps.nx, j = 1:ks.ps.ny
        rho0 = init_density(ks.ps.x[i, j], ks.ps.y[i, j])

        ctr[i, j].prim .= [rho0, 0.0, 0.0, 1.0]
        ctr[i, j].w .= prim_conserve(ctr[i, j].prim, ks.gas.γ)
        ctr[i, j].f .= maxwellian(ks.vs.u, ks.vs.v, ctr[i, j].prim)
    end
    for i in eachindex(a1face)
        a1face[i].fw .= 0.0
        a1face[i].ff .= 0.0
    end
    for i in eachindex(a2face)
        a2face[i].fw .= 0.0
        a2face[i].ff .= 0.0
    end
end

init_field!(ks, ctr, a1face, a2face)

t = 0.0
dt = timestep(ks, ctr, t)
nt = Int(ks.set.maxTime ÷ dt) + 1
res = zero(ks.ib.wL)

@showprogress for iter = 1:20#nt
    reconstruct!(ks, ctr)
    #evolve!(ks, ctr, a1face, a2face, dt; mode = Symbol(ks.set.flux), bc = Symbol(ks.set.boundary))

    # horizontal flux
    @inbounds Threads.@threads for j = 1:ks.pSpace.ny
        for i = 2:ks.pSpace.nx
            flux_kfvs!(
                a1face[i, j].fw,
                a1face[i, j].ff,
                ctr[i-1, j].f,
                ctr[i, j].f,
                ks.vSpace.u,
                ks.vSpace.v,
                ks.vSpace.weights,
                dt,
                a1face[i, j].len,
            )
        end
    end

    # vertical flux
    vn = ks.vSpace.v
    vt = -ks.vSpace.u
    @inbounds Threads.@threads for j = 2:ks.pSpace.ny
        for i = 1:ks.pSpace.nx
            wL = KitBase.local_frame(ctr[i, j-1].w, 0., 1.)
            wR = KitBase.local_frame(ctr[i, j].w, 0., 1.)
            swL = KitBase.local_frame(ctr[i, j-1].sw[:, 2], 0., 1.)
            swR = KitBase.local_frame(ctr[i, j].sw[:, 2], 0., 1.)

            KitBase.flux_kfvs!(
                a2face[i, j].fw,
                a2face[i, j].ff,
                ctr[i, j-1].f,
                ctr[i, j].f,
                vn,
                vt,
                ks.vSpace.weights,
                dt,
                a2face[i, j].len,
            )

            a2face[i, j].fw .= KitBase.global_frame(a2face[i, j].fw, 0., 1.)
        end
    end

    # boundary flux
    @inbounds Threads.@threads for j = 1:ks.pSpace.ny
        KitBase.flux_boundary_maxwell!(
            a1face[1, j].fw,
            a1face[1, j].ff,
            ks.ib.bcL,
            ctr[1, j].f,
            ks.vSpace.u,
            ks.vSpace.v,
            ks.vSpace.weights,
            dt,
            ctr[1, j].dy,
            1.,
        )

        KitBase.flux_boundary_maxwell!(
            a1face[ks.pSpace.nx+1, j].fw,
            a1face[ks.pSpace.nx+1, j].ff,
            ks.ib.bcL,
            ctr[ks.pSpace.nx, j].f,
            ks.vSpace.u,
            ks.vSpace.v,
            ks.vSpace.weights,
            dt,
            ctr[ks.pSpace.nx, j].dy,
            -1.,
        )
    end

    @inbounds Threads.@threads for i = 1:ks.pSpace.nx
        KitBase.flux_boundary_maxwell!(
            a2face[i, 1].fw,
            a2face[i, 1].ff,
            ks.ib.bcL,
            ctr[i, 1].f,
            vn,
            vt,
            ks.vSpace.weights,
            dt,
            ctr[i, 1].dx,
            1,
        )
        a2face[i, 1].fw .= KitBase.global_frame(a2face[i, 1].fw, 0., 1.)

        KitBase.flux_boundary_maxwell!(
            a2face[i, ks.pSpace.ny+1].fw,
            a2face[i, ks.pSpace.ny+1].ff,
            ks.ib.bcL,
            ctr[i, ks.pSpace.ny].f,
            vn,
            vt,
            ks.vSpace.weights,
            dt,
            ctr[i, ks.pSpace.ny].dy,
            -1,
        )
        a2face[i, ks.pSpace.ny+1].fw .= KitBase.global_frame(
            a2face[i, ks.pSpace.ny+1].fw,
            0.,
            1.,
        )
    end

    update!(ks, ctr, a1face, a2face, dt, res; coll = Symbol(ks.set.collision), bc = Symbol(ks.set.boundary))

    t += dt
end

plot_contour(ks, ctr)

max_stars_repo_name: CSMMLab/Flowmachine

max_issues_repo_name: CSMMLab/Flowmachine

max_forks_repo_name: CSMMLab/Flowmachine


text:
# This example shows how to use custom datatypes and reduction operators
# It computes the variance in parallel in a numerically stable way

using MPI, Statistics

MPI.Init()
const comm = MPI.COMM_WORLD
const root = 0

# Define a custom struct
# This contains the summary statistics (mean, variance, length) of a vector
struct SummaryStat
    mean::Float64
    var::Float64
    n::Float64
end
function SummaryStat(X::AbstractArray)
    m = mean(X)
    v = varm(X,m, corrected=false)
    n = length(X)
    SummaryStat(m,v,n)
end

# Define a custom reduction operator
# this computes the pooled mean, pooled variance and total length
function pool(S1::SummaryStat, S2::SummaryStat)
    n = S1.n + S2.n
    m = (S1.mean*S1.n + S2.mean*S2.n) / n
    v = (S1.n * (S1.var + S1.mean * (S1.mean-m)) +
         S2.n * (S2.var + S2.mean * (S2.mean-m)))/n
    SummaryStat(m,v,n)
end

X = randn(10,3) .* [1,3,7]'

# Perform a scalar reduction
summ = MPI.Reduce(SummaryStat(X), pool, root, comm)

if MPI.Comm_rank(comm) == root
    @show summ.var
end

# Perform a vector reduction:
# the reduction operator is applied elementwise
col_summ = MPI.Reduce(mapslices(SummaryStat,X,dims=1), pool, root, comm)

if MPI.Comm_rank(comm) == root
    col_var = map(summ -> summ.var, col_summ)
    @show col_var
end


text:
function and_combine(::Type{BooleanDetector}, arg0::Collection)
    return jcall(BooleanDetector, "andCombine", BooleanDetector, (Collection,), arg0)
end

function and_combine(::Type{BooleanDetector}, arg0::Vector{EventDetector})
    return jcall(BooleanDetector, "andCombine", BooleanDetector, (Vector{EventDetector},), arg0)
end

function g(obj::BooleanDetector, arg0::SpacecraftState)
    return jcall(obj, "g", jdouble, (SpacecraftState,), arg0)
end

function get_detectors(obj::BooleanDetector)
    return jcall(obj, "getDetectors", List, ())
end

function init(obj::BooleanDetector, arg0::SpacecraftState, arg1::AbsoluteDate)
    return jcall(obj, "init", void, (SpacecraftState, AbsoluteDate), arg0, arg1)
end

function not_combine(::Type{BooleanDetector}, arg0::EventDetector)
    return jcall(BooleanDetector, "notCombine", NegateDetector, (EventDetector,), arg0)
end

function or_combine(::Type{BooleanDetector}, arg0::Collection)
    return jcall(BooleanDetector, "orCombine", BooleanDetector, (Collection,), arg0)
end

function or_combine(::Type{BooleanDetector}, arg0::Vector{EventDetector})
    return jcall(BooleanDetector, "orCombine", BooleanDetector, (Vector{EventDetector},), arg0)
end


text:
abstract type FluxLayer end

layertype(l::T) where T = T

# Types for layers with parameters (e.g. weights and biases) and with similar handling
# w.r.t what shape of parameters means in terms of number of inputs and number of outputs
abstract type FluxParLayer <: FluxLayer end
struct FluxDense <: FluxParLayer end
layertype(l::Dense) = FluxDense()

abstract type FluxRecurrent <:FluxParLayer end
struct FluxRnn <: FluxRecurrent end
struct FluxLstm <: FluxRecurrent end
struct FluxGru <: FluxRecurrent end
layertype(l::Flux.Recur{<:Flux.RNNCell}) = FluxRnn()
layertype(l::Flux.Recur{<:Flux.LSTMCell}) = FluxLstm()
layertype(l::Flux.Recur{<:Flux.GRUCell}) = FluxGru()

abstract type FluxConvolutional{N} <: FluxParLayer end
struct FluxConv{N} <: FluxConvolutional{N} end
struct FluxConvTranspose{N}  <: FluxConvolutional{N} end
struct FluxDepthwiseConv{N} <: FluxConvolutional{N}  end
struct FluxCrossCor{N} <: FluxConvolutional{N}  end
layertype(l::Conv{N}) where N = FluxConv{N}()
layertype(l::ConvTranspose{N}) where N = FluxConvTranspose{N}()
layertype(l::DepthwiseConv{N}) where N = FluxDepthwiseConv{N}()
layertype(l::CrossCor{N}) where N = FluxCrossCor{N}()

abstract type FluxTransparentLayer <: FluxLayer end
# Invariant layers with parameters, i.e nin == nout always and parameter selection must
# be performed
abstract type FluxParInvLayer <: FluxTransparentLayer end
struct FluxDiagonal <: FluxParInvLayer end
struct FluxLayerNorm <: FluxParInvLayer end
abstract type FluxParNorm <: FluxParInvLayer end
struct FluxBatchNorm <: FluxParNorm end
struct FluxInstanceNorm <: FluxParNorm end
struct FluxGroupNorm <: FluxParNorm end

layertype(l::Flux.Diagonal) = FluxDiagonal()
layertype(l::LayerNorm) = FluxLayerNorm()
layertype(l::BatchNorm) = FluxBatchNorm()
layertype(l::InstanceNorm) = FluxInstanceNorm()
layertype(l::GroupNorm) = FluxGroupNorm()

# Transparent layers, i.e nin == nout always and there are no parameters
struct FluxNoParLayer <: FluxTransparentLayer end
layertype(l::MaxPool) = FluxNoParLayer()
layertype(l::MeanPool) = FluxNoParLayer()
layertype(l::Dropout) = FluxNoParLayer()
layertype(l::AlphaDropout) = FluxNoParLayer()

# Compositions? Might not have any common methods...
# MaxOut, Chain?


text:
# This is an implementation of point mass dynamics
# PointMass is a mutable structure which only contains data.
# We have functions that work on this data structure, updating the values of the member variables.

using LinearAlgebra

mutable struct PointMass
    m::Float64; # Mass of the rigid body.

    # Internal variables
    x::Vector{Float64};    # State vector

    function PointMass(mass::Float64)
        this = new();
        this.m = mass;
        this.In = Inertia;
        return this;
    end
end

function pmDyn(PM::PointMass,
    extF::extForces, Fc::Vector{Float64}, GravityInInertial::Vector{Float64})

    # States are: 6
    # x = [x,y,z]:inertial, [u,v,w]:inertial
    # extF and Fc are described in the inertial frame.

    xdot = Vector{Float64}(undef,6)

    unconstrainedF = genExtF_PM(PM,extF,GravityInInertial)
    TotalForce = unconstrainedF[1:3] + Fc[1:3]

    # Udwadia's formulation using 6 States
    xdot[1:3] = x[4:6]
    xdot[4:6] = TotalForce/PM.m
    return xdot
end

function initialisePointMass(p::PointMass,x0::Vector{Float64})
    # Initialise point mass with x0
    p.x = x0
    return RB
end

function updatePointMass(p::PointMass,x::Vector{Float64})
    p.x = x
    return p
end

max_stars_repo_name: isrlab/Multibody-Dynamics

max_issues_repo_name: isrlab/Multibody-Dynamics

max_forks_repo_name: isrlab/Multibody-Dynamics


text:
export Constant

"""
    Constant(x)

Constant sequence. Always returns `x`.
"""
struct Constant <: Sequence
  x
end

sample(seq::Constant, t::AbstractVector) = fill(seq.x, length(t))


text:
using CompScienceMeshes
using BEAST
using Test

sphere = readmesh(joinpath(dirname(@__FILE__),"assets","sphere5.in"))
numcells(sphere)

κ = 2π
direction = point(0,0,1)
f = BEAST.HH3DPlaneWave(direction, κ)

v1 = f(point(0,0,0))
v2 = f(point(0,0,0.5))

@test v1 ≈ +1
@test v2 ≈ -1

import BEAST.∂n
p = ∂n(f)

s = chart(m,first(cells(m)))
c = neighborhood(s, [1,1]/3)

r = cartesian(c)
n = normal(s)

w1 = p(c)
w2 = -im*κ*dot(direction, n)*f(r)

w1 ≈ w2

N = BEAST.HH3DHyperSingularFDBIO(im*κ)
X = BEAST.lagrangec0d1(sphere)

numfunctions(X)

Nxx = assemble(N, X, X)

@test size(Nxx) == (numfunctions(X), numfunctions(X))


text:
# This example shows a simple use for a channel and a feeder

using Distributed
using Logging
using Base.CoreLogging

# Set Debug for Logging
global_logger(ConsoleLogger(stdout, Logging.Debug))

function read_and_print(channel::Channel{Tuple{Vector{Int}, Int}})

    @debug("Running updater.")

    while true

        x, i = try
            take!(channel);
        catch e

            if isa(e, InvalidStateException)
                break
            end

            @warn("Something wrong when reading from channel. Will skip.", e)

            continue

        end

        @debug("Updater has read values from channel. x = $(x) and i = $(i).")

    end

    @debug("Update channel closed. Exiting thread.")

end

c = Channel{Tuple{Vector{Int}, Int}}(0)

# Task waits for the channel to be fed
task1 = @async read_and_print(c)

# Feeding the channel
task2 = @async for i = 1:10
    put!(c, (ones(Int, 3), i))
end

sleep(1)

# Wrong feeding
put!(c, 1)

# This should close `task1`
close(c)

# Just to be sure
sleep(1)

println("Task1 done = $(istaskdone(task1))")

println("Task2 done = $(istaskdone(task2))")


text:
#=
MIT License

Copyright (c) 2019 Jagadeesh Kotra

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
=#

using ArgParse

function parse_commandline_args()
    s = ArgParseSettings()
    @add_arg_table s begin

        "--bank"
           help = "bank / balance"
           required = true
           arg_type = Float64

        "--stake"
           help = "stake/Bet Amount to start with, also considered as base"
           required = true
           arg_type = Float64

        "--rolls"
           help = "N no of rolls"
           required = true
           arg_type = Int
           default = 1
        "--multiplier"
           help = "If won then, stake*multiplier"
           required = true
           arg_type = Float16

        "--report-freq"
           help = "Gap of rolls between each line printed to STDOUT"
           required = true
           arg_type = Int
           default = 1
    end
    return parse_args(s)
end

parsed_args = parse_commandline_args()

bank = parsed_args["bank"]

stake_base = parsed_args["stake"]
stake = 0

#for our 1st roll of Dice
stake = stake_base

multiplier = parsed_args["multiplier"]

w = 0
l = 0

for i = 1:parsed_args["rolls"]
    global bank -= stake

    if rand()*100 > 50.0

        global w += 1

        global bank += stake*multiplier
        global stake = stake_base

    else
        global l += 1
        global stake = stake*2
    end

    if i % parsed_args["report-freq"] == 0
        println(i," ",bank)
    end

    if stake > bank
        println(i," ",bank, " Bankrupt! Stake:",stake)
        break
    end

end

println("---END---")
println(w,"-",l," Win %:", 100*(w/(w+l)))

max_stars_repo_name: jkotra/LuckSimulations

max_issues_repo_name: jkotra/LuckSimulations

max_forks_repo_name: jkotra/LuckSimulations


text:
module HyperRectangles
using Compat

export AbstractHyperRectangle, HyperRectangle

abstract AbstractHyperRectangle{T, N}

type HyperRectangle{T, N} <: AbstractHyperRectangle{T, N}
    min::Vector{T}
    max::Vector{T}
end

function HyperRectangle{T}(min::Vector{T}, max::Vector{T})
    n = length(max)
    m = length(min)
    n == m || error("min and max vector lengths are different!")
    HyperRectangle{T, n}(min, max)
end

@inline Base.max(b::HyperRectangle) = b.max
@inline Base.min(b::HyperRectangle) = b.min

function (==){T1, T2, N}(b1::HyperRectangle{T1, N}, b2::HyperRectangle{T2, N})
    b1.min == b2.min && b1.max == b2.max
end

@inline isequal(b1::HyperRectangle, b2::HyperRectangle) = b1 == b2

function Base.contains{T1, T2, N}(b1::HyperRectangle{T1,N}, b2::HyperRectangle{T2, N})
    for i = 1:N
        b2.max[i] <= b1.max[i] && b2.min[i] >= b1.min[i] || return false
    end
    true
end

@inline Base.in(b1::HyperRectangle, b2::HyperRectangle) = contains(b2, b1)

# Splits an HyperRectangle into two new ones along an axis
# at a given axis value
function Base.split{T, N}(b::HyperRectangle{T,N}, axis::Int, value::T)
    b1max = copy(b.max)
    b1max[axis] = value

    b2min = copy(b.min)
    b2min[axis] = value

    return HyperRectangle{T, N}(b.min, b1max),
           HyperRectangle{T, N}(b2min, b.max)
end

#
# Based on:
# Near real-time CSG rendering using tree normalization and geometric pruning
# Computer Graphics and Applications, IEEE
#
function Base.union{T,N}(h1::HyperRectangle{T,N}, h2::HyperRectangle{T,N})
    mins = zeros(T, N)
    maxs = zeros(T, N)
    for i = 1:N
        mins[i] = min(h1.min[i], h2.min[i])
        maxs[i] = max(h1.max[i], h2.max[i])
    end
    HyperRectangle{T,N}(mins, maxs)
end

function Base.diff(h1::HyperRectangle, h2::HyperRectangle)
    h1
end

function Base.intersect{T,N}(h1::HyperRectangle{T,N}, h2::HyperRectangle{T,N})
    mins = zeros(T, N)
    maxs = zeros(T, N)
    for i = 1:N
        mins[i] = max(h1.min[i], h2.min[i])
        maxs[i] = min(h1.max[i], h2.max[i])
    end
    HyperRectangle{T,N}(mins, maxs)
end

if VERSION >= v"0.4.0-"
    function Base.call{T,N}(::Type{HyperRectangle{T,N}})
        max = fill(typemin(T), N)
        min = fill(typemax(T), N)
        HyperRectangle{T,N}(min, max)
    end
end

# submodules
include("Relations.jl")
include("Operations.jl")

end # module


text:
using DFControl, Test

testassetspath = joinpath(testdir, "testassets")
testjobpath = joinpath(testassetspath, "test_job")

@testset "Structure manipulation" begin
    job = Job(testjobpath, "localhost_test")

    @test length(Structures.symmetry_operators(job.structure)[1]) == 48

    job.structure = Structures.create_supercell(job.structure, 1, 0, 0, make_afm = true)
    @test length(Structures.symmetry_operators(job.structure)[1]) == 12
    @test length(job.structure) == 2
    @test isapprox(job.structure[:Ni1][1].magnetization, [0,0,-0.1])
    @test isapprox(job.structure[:Ni][1].magnetization, [0,0,0.1])
    @test isapprox(map(x->x.magnetization, job.structure[Structures.element(:Ni)]), [[0,0,0.1], [0,0,-0.1]])
    ngl = Structures.niggli_reduce(job.structure)
    @test Structures.volume(job.structure) == Structures.volume(ngl)
    @test job.structure.cell != ngl.cell
    prev = job.structure[Structures.element(:Ni)][1].position_cryst[1]
    prev_vol = Structures.volume(job.structure)
    Structures.scale_cell!(job.structure, diagm(0=>[2,1,1]))
    @test job.structure[Structures.element(:Ni)][1].position_cryst[1] == 0.5prev
    @test 2*prev_vol == Structures.volume(job.structure)

    out_str = outputdata(job)["scf"][:initial_structure]
    Structures.update_geometry!(job.structure, out_str)
    @test out_str == job.structure
end

@testset "supercell" begin
    job = Job(testjobpath, "localhost_test")
    struct2 = Structures.create_supercell(job.structure, 1, 2, 1)
    newpositions = [at.position_cart for at in struct2.atoms]
    oldposition = job.structure.atoms[1].position_cart
    cell_ = job.structure.cell

    @test job.structure.atoms[1].position_cart ==
          job.structure.cell' * job.structure.atoms[1].position_cryst
    @test oldposition == newpositions[1]
    @test oldposition + cell_[:, 1] ∈ newpositions
    @test oldposition + 2 * cell_[:, 2] ∈ newpositions
    @test oldposition + 1 * cell_[:, 3] ∈ newpositions
    @test oldposition + (1 * cell_[:, 3] + 1 * cell_[:, 2]) ∈ newpositions
end

@testset "registry" begin
    job = Job(testjobpath,  "localhost_test")
    rm(testjobpath, recursive=true)
    prevlen = length(Jobs.registered_jobs())
    save(job)
    @test length(Jobs.registered_jobs()) == prevlen + 1
    @test Jobs.registered_jobs()[end][1] == job.dir
end

# @testset "execs" begin
#     job = Job(testjobpath,  "localhost_test")
#     for c in job.calculations
#         e = c.exec
#         if e.exec == "pw.x"
#             Calculations.set_flags!(e, :nk => 4)
#             e.dir = "test/test"
#         elseif e.exec == "mpirun"
#             Calculations.set_flags!(e, :np => 4)
#             @test_throws ErrorException Calculations.set_flags!(e, :nasdf => 4)
#         end
#     end
#     ex = job.calculations[1].execs[2]
#     mpi = job.calculations[1].execs[1]
#     @test ex.dir == "test/test"
#     @test ex.flags[1].symbol == :nk
#     @test ex.flags[1].value == 4
#     @test mpi.flags[1].symbol == :np
#     @test mpi.flags[1].value == 4

#     save(job)
#     job = Job(testjobpath, "localhost_test")
#     ex = job.calculations[1].execs[2]
#     mpi = job.calculations[1].execs[1]
#     @test ex.flags[1].symbol == :nk
#     @test ex.flags[1].value == 4
#     @test mpi.flags[1].symbol == :np
#     @test mpi.flags[1].value == 4

#     for c in job.calculations
#         for e in c.execs
#             if e.exec == "pw.x"
#                 Calculations.rm_flags!(e, :nk)
#                 e.dir = "test/test"
#             elseif e.exec == "mpirun"
#                 Calculations.rm_flags!(e, :np)
#             end
#         end
#     end

#     @test isempty(ex.flags)
#     @test isempty(mpi.flags)
#     save(job)
# end

@testset "calculation management" begin
    job = Job(testjobpath, "localhost_test")
    ncalcs = length(job.calculations)
    t = pop!(job, "scf")
    @test length(job.calculations) == ncalcs - 1
    @test job[2].name == "bands"
    insert!(job, 2, t)
    @test job.calculations[2].name == "scf"

    Jobs.set_flow!(job, "" => false, "scf" => true)
    @test job["scf"].run == true
    @test job["nscf"].run == true
    @test job["bands"].run == false
    push!(job.header, "module load testmod")
    Jobs.set_headerword!(job, "testmod" => "testmod2")
    @test job.header[end] == "module load testmod2"

end

@testset "versioning" begin
    job = Job(testjobpath, "localhost_test")
    job[:nbnd] = 30
    curver = Client.last_version(job)
    save(job)
    @test job.version == curver + 1
    save(job)
    @test ispath(joinpath(job, Jobs.VERSION_DIR_NAME))
    @test ispath(joinpath(job, Jobs.VERSION_DIR_NAME, "$(curver+1)"))
    job[:nbnd] = 40
    save(job)
    @test job.version == curver + 3
    @test job["scf"][:nbnd] == 40
    switch_version!(job, curver + 1)
    @test job.version == curver + 1
    # @test DFControl.last_version(job) == 2
    @test job["scf"][:nbnd] == 30
    rm_version!(job, curver + 1)
    @test !ispath(joinpath(Jobs.main_job_dir(job), Jobs.VERSION_DIR_NAME,
                           "$(curver + 1)"))
    @test ispath(joinpath(Jobs.main_job_dir(job), Jobs.VERSION_DIR_NAME,
                          "$(curver + 2)"))
end

rm(testjobpath, recursive=true)
# function copy_outfiles()
#     for f in readdir(joinpath(testassetspath, "outputs"))
#         cp(joinpath(testassetspath, "outputs", f), joinpath(testjobpath, f); force = true)
#     end
# end

# copy_outfiles()

# job = Job(testjobpath);
# job4 = Job(testjobpath)

# #output stuff
# out = outputdata(job; print = false, onlynew = false);
# @test haskey(out, "nscf")
# @test haskey(out["nscf"], :fermi)

# @test isconverged(job["scf"])
# nscf = DFControl.calculation(job, "nscf")
# nscf2 = Calculation(nscf, "nscf2";
#                       data = [:testdata => (:testoption, "test"),
#                               :k_points => (:blabla, [1, 1, 1, 1, 1, 1])])

# @test data(nscf2, :testdata).option == :testoption
# @test data(nscf2, :testdata).data == "test"
# @test data(nscf2, :testdata).name == :testdata
# @test data(nscf2, :k_points).option == :blabla
# @test data(nscf2, :k_points).data == [1, 1, 1, 1, 1, 1]
# @test data(job["scf"], :k_points).option == :automatic
# @test_throws ErrorException job["bladbsflk"]
# @test_throws ErrorException nscf2[:bladkfj]

# set_kpoints!(nscf2, (3, 3, 3); print = false)
# @test data(nscf2, :k_points).data == DFControl.kgrid(3, 3, 3, :nscf)

# set_kpoints!(nscf2, [(3.0, 3.0, 3.0, 1.0), (3.0, 3.0, 3.0, 1.0)]; print = false)
# @test data(nscf2, :k_points).option == :crystal_b
# @test data(nscf2, :k_points).data == [(3.0, 3.0, 3.0, 1.0), (3.0, 3.0, 3.0, 1.0)]

# set_kpoints!(nscf2, (3, 3, 3, 0, 0, 1); print = false)
# @test data(nscf2, :k_points).option == :automatic
# @test data(nscf2, :k_points).data == [3, 3, 3, 0, 0, 1]

# fermi = outputdata(job, "nscf")[:fermi]
# @test fermi == 17.4572

# tcalc = gencalc_bands(job["bands"],
#                       [(0.5, 0.0, 0.5, 10.0), (0.0, 0.0, 0.0, 10.0), (0.5, 0.5, 0.5, 1.0)];
#                       name = "bands2")

# push!(job, tcalc)
# @test job["bands2"][:calculation] == "bands"
# @test data(job["bands2"], :k_points).data ==
#       [(0.5, 0.0, 0.5, 10.0), (0.0, 0.0, 0.0, 10.0), (0.5, 0.5, 0.5, 1.0)]

# tcalc = gencalc_nscf(job["nscf"], (10, 10, 10); name = "nscf2")
# push!(job, tcalc)
# @test job["nscf2"][:calculation] == "nscf"
# push!(job, gencalc_scf(job["scf"], (5, 5, 5, 1, 1, 1); name = "1scf2"))
# @test job["1scf2"][:calculation] == "scf"

# wanflags = [:write_hr => true, :wannier_plot => true]

# set_projections!(job, :Pt => ["s", "p", "d"])
# wancalc = gencalc_wan(job["nscf"], structure(job), fermi - 7.0, wanflags...; Epad = 5.0)
# append!(job, wancalc)
# @test job["wan"][:write_hr] == job["wan"][:wannier_plot] == true

# copy_outfiles()
# wanout = outputdata(job, "wan")
# @test length(wanout[:final_state]) == 20
# @test length(wanout[:wannierise]) == 203
# job.calculations = job.calculations[1:4]

# job[:nbnd] = 300
# @test job["scf"][:nbnd] == 300
# @test job["nscf"][:nbnd] == 300
# rm_flags!(job, :nbnd)
# @test !haskey(job["scf"].flags, :nbnd)

# job["scf"][:nbnd] = 300
# @test job["scf"][:nbnd] == 300
# rm_flags!(job, :nbnd)

# job["nscf"][:Hubbard_U] = [1.0]
# @test job["nscf"][:Hubbard_U] == [1.0]

# job["nscf"][:starting_magnetization] = [[1.0, 2.0, 1.0]]
# @test job["nscf"][:starting_magnetization] == [[1.0, 2.0, 1.0]]

# job["nscf"][:Hubbard_J] = [0 1 2]
# @test job["nscf"][:Hubbard_J] == [0 1 2]
# rm_flags!(job, :Hubbard_J, :starting_magnetization, :Hubbard_U)
# set_magnetization!(job, :Pt => [0.0, 0.0, 1.0])

# push!.((job,), gencalc_wan(nscf, structure(job), fermi - 7.0, wanflags...; Epad = 5.0))
# #TODO: add test with the new wancalc from projwfc

# set_flow!(job, "" => false)
# @test job.calculations[1].run == false
# set_flow!(job, "nscf" => true, "bands" => true)
# @test job.calculations[3].run

# save(job)
# job2 = Job(dir)
# begin
#     for (calc, calc2) in zip(job.calculations, job2.calculations)
#         for (f, v) in calc.flags
#             if f in (:Hubbard_J, :pseudo_dir, :wannier_plot)
#                 continue
#             else
#                 @test v == calc2.flags[f]
#             end
#         end
#         for (b1, b2) in zip(calc.data, calc2.data)
#             for n in fieldnames(typeof(b1))
#                 @test getfield(b1, n) == getfield(b2, n)
#             end
#         end
#     end
# end
# copy_outfiles()
# set_cutoffs!(job)
# @test job["scf"][:ecutwfc] == 32.0

# set_pseudos!(job, :Pt => Pseudo("Pt.UPF", joinpath(testdir, "testassets/pseudos")))
# @test job.structure.atoms[1].pseudo ==
#       Pseudo("Pt.UPF", joinpath(testdir, "testassets/pseudos"))
# set_pseudos!(job, :Pt, :test)
# @test job.structure.atoms[1].pseudo ==
#       Pseudo("Pt.UPF", joinpath(testdir, "testassets/pseudos"))

# set_pseudos!(job, :test)
# @test job.structure.atoms[1].pseudo ==
#       Pseudo("Pt.UPF", joinpath(testdir, "testassets/pseudos"))

# testorbs = ["s", "p"]
# set_projections!(job, :Pt => testorbs)
# @test convert.(String, [p.orb for p in Iterators.flatten(projections.(job.structure.atoms :Pt)))]) ==
#       testorbs
# set_wanenergies!(job, nscf, fermi - 7.0; Epad = 3.0)

# @test job["wanup"][:dis_froz_max] == 13.2921
# @test job["wanup"][:dis_win_max] == 16.292099999999998

# set_execflags!(job, "pw.x", :nk => 230, :ndiag => 3)
# @test DFControl.getfirst(x -> x.symbol == :nk, execs(job["nscf"])[2].flags).value == 230
# @test DFControl.getfirst(x -> x.symbol == :ndiag, execs(job["nscf"])[2].flags).value == 3
# rm_execflags!(job, "pw.x", :nk)
# @test isempty(filter(x -> x.symbol == :nk, execs(job["nscf"])[2].flags))

# set_execdir!(job, "pw.x", joinpath(homedir(), "bin"))
# @test execs(job["nscf"])[2].dir == joinpath(homedir(), "bin")

# set_name!(job["nscf"], "test")
# set_name!(job["test"], "nscf")

# set_serverdir!(job, "localhost")
# @test job.server_dir == "localhost"

# set_headerword!(job, "defpart" => "frontend"; print = false)
# @test any(occursin.("frontend", job.header))

# report = progressreport(job; onlynew = false, print = false)
# @test report[:fermi] == 17.4572
# @test length(report[:accuracy]) == 9
# oldat = job.structure.atoms[1]
# copy_outfiles()
# newatompos = outputdata(job, "vc_relax"; onlynew = false)[:final_structure]
# job.structure = newatompos
# push!(job.structure.atoms, oldat)

# cp(joinpath(testdir, "testassets", "pseudos"),
#    joinpath(testdir, "testassets", "pseudos_copy"); force = true)
# set_pseudos!(job, :Si => Pseudo("Si.UPF", joinpath(testdir, "testassets", "pseudos_copy")))
# save(job)
# @test ispath(joinpath(job.dir, "Si.UPF"))
# @test job.structure.atoms :Si)[1].pseudo == Pseudo("Si.UPF", job.dir)
# rm(joinpath(testdir, "testassets", "pseudos_copy"); recursive = true)

# job["nscf"][:occupations] = "smearing"
# job["nscf"][:degauss] = 2.0
# job["nscf"][:smearing] = "mp"
# projwfc = gencalc_projwfc(job["nscf"], -20, 10, 0.05)
# @test projwfc[:Emin] == -20
# @test projwfc[:degauss] == job["nscf"][:degauss]
# @test projwfc[:ngauss] == 1

# set_Hubbard_U!(job, :Si => 1.0)
# @test job.structure.atoms :Si)[1].dftu.U == 1.0

# prev_a = job.structure.cell[1, :]
# prev_b = job.structure.cell[2, :]
# prev_c = job.structure.cell[3, :]
# prev_pos = position_cart.(job.structure.atoms)
# scale_cell!(job, [2 0 0; 0 2 0; 0 0 2])
# @test prev_a .* 2 == job.structure.cell[1, :]
# @test prev_b .* 2 == job.structure.cell[2, :]
# @test prev_c .* 2 == job.structure.cell[3, :]
# for (p, at) in zip(prev_pos, job.structure.atoms)
#     @test round.(DFControl.ustrip.(p * 2), digits = 3) ==
#           round.(DFControl.ustrip.(at.position_cart), digits = 3)
# end

# set_magnetization!(job, :Pt => [1.0, 0.0, 0.0])
# @test magnetization(job.structure.atoms :Pt)[1]) == DFControl.Vec3(1.0, 0.0, 0.0)

# at             = job.structure.atoms[1]
# c              = job.structure.cell
# orig_pos_cart  = at.position_cart
# orig_pos_cryst = at.position_cryst

# set_position!(at, orig_pos_cart .+ c * [0.1, 0.1, 0.1], c)
# @test isapprox(at.position_cryst, orig_pos_cryst .+ [0.1, 0.1, 0.1])

# set_position!(at, orig_pos_cryst .+ [0.1, 0.1, 0.1], c)
# @test at.position_cart == c * at.position_cryst

# at2 = job.structure.atoms[2]
# p1, p2 = at.position_cart, at.position_cart)
# mid = (p1 + p2) / 2
# bondlength = distance(at, at2)
# scale_bondlength!(at, at2, 0.5, c)

# @test isapprox((at.position_cart + at.position_cart)) / 2, mid)
# @test isapprox(distance(at, at2), bondlength / 2)
# copy_outfiles()
# @test isapprox(bandgap(job), 2.6701999999999995)

# t = job["nscf"]
# curlen = length(job.calculations)
# id = findfirst(x -> x.name == "nscf", job.calculations)
# n = pop!(job, "nscf")
# @test n == t
# @test length(job.calculations) == curlen - 1

# nscf = job["scf"]
# rm_flags!(job, :nspin, :lda_plus_u, :noncolin)
# set_magnetization!(job, :Pt => [0.2, 1.0, 0.2])
# DFControl.sanitize_flags!(job)
# set_Hubbard_U!(job, :Pt => 2.3)
# DFControl.sanitize_magnetization!(job)
# DFControl.set_hubbard_flags!.(filter(x -> eltype(x) == QE,
#                                      job.calculations), (job.structure,))
# DFControl.set_starting_magnetization_flags!.(filter(x -> eltype(x) == QE,
#                                                     job.calculations),
#                                              (job.structure,))
# @test job["scf"][:lda_plus_u]
# @test job["scf"][:noncolin]
# @test job["scf"][:lda_plus_u_kind] == 1

# rm_flags!(job, :nspin, :lda_plus_u, :noncolin)
# set_magnetization!(job, :Pt => [0.0, 0.0, 0.5])
# DFControl.sanitize_flags!(job)
# DFControl.sanitize_magnetization!(job)
# DFControl.set_hubbard_flags!.(filter(x -> eltype(x) == QE,
#                                      job.calculations), (job.structure,))
# DFControl.set_starting_magnetization_flags!.(filter(x -> eltype(x) == QE,
#                                                     job.calculations),
#                                              (job.structure,))
# @test job["scf"][:nspin] == 2

# using LinearAlgebra
# new_str = create_supercell(structure(job), 1, 0, 0)
# prevcell = job.structure.cell
# @test norm(new_str.cell[:, 1]) == norm(prevcell[:, 1]) * 2
# prevlen = length(job.structure.atoms)
# @test length(new_str.atoms) == 2 * prevlen
# prevlen_Pt = length(job.structure.atoms :Pt))

# set_magnetization!(job, :Pt => [0, 0, 1])
# orig_projs = projections(job.structure.atoms :Pt)[1])

# job.structure = create_supercell(structure(job), 1, 0, 0; make_afm = true)
# DFControl.sanitize_magnetization!(job)
# DFControl.set_starting_magnetization_flags!.(filter(x -> eltype(x) == QE,
#                                                     job.calculations),
#                                              (job.structure,))
# @test length(job.structure.atoms :Pt)) == prevlen_Pt
# @test length(job.structure.atoms :Pt1)) == prevlen_Pt

# @test magnetization(job.structure.atoms :Pt)[1]) == [0, 0, 1]
# @test magnetization(job.structure.atoms :Pt1)[1]) == [0, 0, -1]

# DFControl.sanitize_projections!(job)
# @test projections(job.structure.atoms :Pt)[1]) != projections(job.structure.atoms :Pt1)[1])

# # job4.server_dir = "/tmp"
# # save(job4)
# # job3 = Job(job4.dir)

# # @test job3.server_dir == "/tmp"
# set_data_option!(job["scf"], :k_points, :blabla; print = false)
# @test data(job["scf"], :k_points).option == :blabla
# set_data_option!(job, :k_points, :test; print = false)
# @test data(job["scf"], :k_points).option == :test

# rm.(DFControl.inpath.(job.calculations))
# set_kpoints!(job["scf"], (6, 6, 6, 1, 1, 1))
# rm(joinpath(job, DFControl.VERSION_DIR_NAME); recursive = true)

# rm(joinpath(DFControl.main_job_dir(job), DFControl.VERSION_DIR_NAME); recursive = true)
# rm.(DFControl.inpath.(job.calculations))

# rm(joinpath(job, ".metadata.jld2"))
# rm(joinpath(job, "pw2wan_wandn.in"))
# rm(joinpath(job, "pw2wan_wanup.in"))
# rm.(joinpath.((job.dir,), filter(x -> occursin("UPF", x), readdir(job.dir))))


text:
module TestLister

using Test
using Transducers: TransducerLister

@testset "show" begin
    @test show(devnull, "text/markdown", TransducerLister()) isa Any
    @test show(devnull, "text/plain", TransducerLister()) isa Any
    @test show(devnull, "text/html", TransducerLister()) isa Any
end

end  # module


text:
function contour_mesh(f::Function,
                      bounds::HyperRectangle,
                      isosurface_value::Number=0.0,
                      resolution::Real=maximum(widths(bounds))/20.0)

    # Extract an isosurface from a vector-input function in 3 dimensions and
    # return it as a mesh.
    window_width = widths(bounds)
    sdf = SignedDistanceField(x -> f(SVector{3, Float64}(x[1], x[2], x[3])),
                              bounds,
                              resolution)
    mesh = HomogenousMesh(sdf, MarchingTetrahedra(isosurface_value))
end

contour_mesh(f::Function,
             lower_bound::Vec,
             upper_bound::Vec,
             isosurface_value::Number=0.0,
             resolution::Real=maximum(upper_bound - lower_bound)/20.0) = (
    contour_mesh(f,
                 HyperRectangle(lower_bound, upper_bound - lower_bound),
                 isosurface_value,
                 resolution))

contour_mesh(f::Function,
             lower_bound::AbstractVector,
             upper_bound::AbstractVector,
             isosurface_value::Number=0.0,
             resolution::Real=maximum(upper_bound - lower_bound)/20.0) = (
    contour_mesh(f,
                 Vec{3, Float64}(lower_bound[1], lower_bound[2], lower_bound[3]),
                 Vec{3, Float64}(upper_bound[1], upper_bound[2], upper_bound[3]),
                 isosurface_value,
                 resolution))


text:
Text provided under a Creative Commons Attribution license, CC-BY.  All code is made available under the FSF-approved BSD-3 license.  (c) Lorena A. Barba, Gilbert F. Forsyth 2017. Thanks to NSF for support via CAREER award #1149784.
[@LorenaABarba](https://twitter.com/LorenaABarba)

12 steps to Navier–Stokes
=====
***

The final two steps in this interactive module teaching beginning [CFD with Python](https://bitbucket.org/cfdpython/cfd-python-class) will both solve the Navier–Stokes equations in two dimensions, but with different boundary conditions.

The momentum equation in vector form for a velocity field $\vec{v}$ is:

$$\frac{\partial \vec{v}}{\partial t}+(\vec{v}\cdot\nabla)\vec{v}=-\frac{1}{\rho}\nabla p + \nu \nabla^2\vec{v}$$

This represents three scalar equations, one for each velocity component $(u,v,w)$. But we will solve it in two dimensions, so there will be two scalar equations.

Remember the continuity equation? This is where the [Poisson equation](./13_Step_10.ipynb) for pressure comes in!

Step 11: Cavity Flow with Navier–Stokes
----
***

Here is the system of differential equations: two equations for the velocity components $u,v$ and one equation for pressure:

$$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y} = -\frac{1}{\rho}\frac{\partial p}{\partial x}+\nu \left(\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2} \right) $$

$$\frac{\partial v}{\partial t}+u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y} = -\frac{1}{\rho}\frac{\partial p}{\partial y}+\nu\left(\frac{\partial^2 v}{\partial x^2}+\frac{\partial^2 v}{\partial y^2}\right) $$

$$\frac{\partial^2 p}{\partial x^2}+\frac{\partial^2 p}{\partial y^2} = -\rho\left(\frac{\partial u}{\partial x}\frac{\partial u}{\partial x}+2\frac{\partial u}{\partial y}\frac{\partial v}{\partial x}+\frac{\partial v}{\partial y}\frac{\partial v}{\partial y} \right)$$

From the previous steps, we already know how to discretize all these terms. Only the last equation is a little unfamiliar. But with a little patience, it will not be hard!

### Discretized equations

First, let's discretize the $u$-momentum equation, as follows:

$$
\begin{split}
& \frac{u_{i,j}^{n+1}-u_{i,j}^{n}}{\Delta t}+u_{i,j}^{n}\frac{u_{i,j}^{n}-u_{i-1,j}^{n}}{\Delta x}+v_{i,j}^{n}\frac{u_{i,j}^{n}-u_{i,j-1}^{n}}{\Delta y} = \\
& \qquad -\frac{1}{\rho}\frac{p_{i+1,j}^{n}-p_{i-1,j}^{n}}{2\Delta x}+\nu\left(\frac{u_{i+1,j}^{n}-2u_{i,j}^{n}+u_{i-1,j}^{n}}{\Delta x^2}+\frac{u_{i,j+1}^{n}-2u_{i,j}^{n}+u_{i,j-1}^{n}}{\Delta y^2}\right)
\end{split}
$$

Similarly for the $v$-momentum equation:

$$
\begin{split}
&\frac{v_{i,j}^{n+1}-v_{i,j}^{n}}{\Delta t}+u_{i,j}^{n}\frac{v_{i,j}^{n}-v_{i-1,j}^{n}}{\Delta x}+v_{i,j}^{n}\frac{v_{i,j}^{n}-v_{i,j-1}^{n}}{\Delta y} = \\
& \qquad -\frac{1}{\rho}\frac{p_{i,j+1}^{n}-p_{i,j-1}^{n}}{2\Delta y}
+\nu\left(\frac{v_{i+1,j}^{n}-2v_{i,j}^{n}+v_{i-1,j}^{n}}{\Delta x^2}+\frac{v_{i,j+1}^{n}-2v_{i,j}^{n}+v_{i,j-1}^{n}}{\Delta y^2}\right)
\end{split}
$$

Finally, the discretized pressure-Poisson equation can be written thus:

$$
\begin{split}
& \frac{p_{i+1,j}^{n}-2p_{i,j}^{n}+p_{i-1,j}^{n}}{\Delta x^2}+\frac{p_{i,j+1}^{n}-2p_{i,j}^{n}+p_{i,j-1}^{n}}{\Delta y^2} = \\
& \qquad \rho \left[ \frac{1}{\Delta t}\left(\frac{u_{i+1,j}-u_{i-1,j}}{2\Delta x}+\frac{v_{i,j+1}-v_{i,j-1}}{2\Delta y}\right) -\frac{u_{i+1,j}-u_{i-1,j}}{2\Delta x}\frac{u_{i+1,j}-u_{i-1,j}}{2\Delta x} - 2\frac{u_{i,j+1}-u_{i,j-1}}{2\Delta y}\frac{v_{i+1,j}-v_{i-1,j}}{2\Delta x} - \frac{v_{i,j+1}-v_{i,j-1}}{2\Delta y}\frac{v_{i,j+1}-v_{i,j-1}}{2\Delta y}\right]
\end{split}
$$

You should write these equations down on your own notes, by hand, following each term mentally as you write it.

As before, let's rearrange the equations in the way that the iterations need to proceed in the code. First, the momentum equations for the velocity at the next time step.

The momentum equation in the $u$ direction:

$$
\begin{split}
u_{i,j}^{n+1} = u_{i,j}^{n} & - u_{i,j}^{n} \frac{\Delta t}{\Delta x} \left(u_{i,j}^{n}-u_{i-1,j}^{n}\right) - v_{i,j}^{n} \frac{\Delta t}{\Delta y} \left(u_{i,j}^{n}-u_{i,j-1}^{n}\right) \\
& - \frac{\Delta t}{\rho 2\Delta x} \left(p_{i+1,j}^{n}-p_{i-1,j}^{n}\right) \\
& + \nu \left(\frac{\Delta t}{\Delta x^2} \left(u_{i+1,j}^{n}-2u_{i,j}^{n}+u_{i-1,j}^{n}\right) + \frac{\Delta t}{\Delta y^2} \left(u_{i,j+1}^{n}-2u_{i,j}^{n}+u_{i,j-1}^{n}\right)\right)
\end{split}
$$

The momentum equation in the $v$ direction:

$$
\begin{split}
v_{i,j}^{n+1} = v_{i,j}^{n} & - u_{i,j}^{n} \frac{\Delta t}{\Delta x} \left(v_{i,j}^{n}-v_{i-1,j}^{n}\right) - v_{i,j}^{n} \frac{\Delta t}{\Delta y} \left(v_{i,j}^{n}-v_{i,j-1}^{n})\right) \\
& - \frac{\Delta t}{\rho 2\Delta y} \left(p_{i,j+1}^{n}-p_{i,j-1}^{n}\right) \\
& + \nu \left(\frac{\Delta t}{\Delta x^2} \left(v_{i+1,j}^{n}-2v_{i,j}^{n}+v_{i-1,j}^{n}\right) + \frac{\Delta t}{\Delta y^2} \left(v_{i,j+1}^{n}-2v_{i,j}^{n}+v_{i,j-1}^{n}\right)\right)
\end{split}
$$

Almost there! Now, we rearrange the pressure-Poisson equation:

$$
\begin{split}
p_{i,j}^{n} = & \frac{\left(p_{i+1,j}^{n}+p_{i-1,j}^{n}\right) \Delta y^2 + \left(p_{i,j+1}^{n}+p_{i,j-1}^{n}\right) \Delta x^2}{2\left(\Delta x^2+\Delta y^2\right)} \\
& -\frac{\rho\Delta x^2\Delta y^2}{2\left(\Delta x^2+\Delta y^2\right)} \\
& \times \left[\frac{1}{\Delta t}\left(\frac{u_{i+1,j}-u_{i-1,j}}{2\Delta x}+\frac{v_{i,j+1}-v_{i,j-1}}{2\Delta y}\right)-\frac{u_{i+1,j}-u_{i-1,j}}{2\Delta x}\frac{u_{i+1,j}-u_{i-1,j}}{2\Delta x} -2\frac{u_{i,j+1}-u_{i,j-1}}{2\Delta y}\frac{v_{i+1,j}-v_{i-1,j}}{2\Delta x}-\frac{v_{i,j+1}-v_{i,j-1}}{2\Delta y}\frac{v_{i,j+1}-v_{i,j-1}}{2\Delta y}\right]
\end{split}
$$

The initial condition is $u, v, p = 0$ everywhere, and the boundary conditions are:

$u=1$ at $y=2$ (the "lid");

$u, v=0$ on the other boundaries;

$\frac{\partial p}{\partial y}=0$ at $y=0$;

$p=0$ at $y=2$

$\frac{\partial p}{\partial x}=0$ at $x=0,2$

Implementing Cavity Flow
----

```python
import numpy
from matplotlib import pyplot, cm
from mpl_toolkits.mplot3d import Axes3D
%matplotlib inline
```

```python
nx = 41
ny = 41
nt = 500
nit = 50
c = 1
dx = 2 / (nx - 1)
dy = 2 / (ny - 1)
x = numpy.linspace(0, 2, nx)
y = numpy.linspace(0, 2, ny)
X, Y = numpy.meshgrid(x, y)

rho = 1
nu = .1
dt = .001

u = numpy.zeros((ny, nx))
v = numpy.zeros((ny, nx))
p = numpy.zeros((ny, nx))
b = numpy.zeros((ny, nx))
```

The pressure Poisson equation that's written above can be hard to write out without typos.  The function `build_up_b` below represents the contents of the square brackets, so that the entirety of the PPE is slightly more manageable.

```python
def build_up_b(b, rho, dt, u, v, dx, dy):

    b[1:-1, 1:-1] = (rho * (1 / dt *
                    ((u[1:-1, 2:] - u[1:-1, 0:-2]) /
                     (2 * dx) + (v[2:, 1:-1] - v[0:-2, 1:-1]) / (2 * dy)) -
                    ((u[1:-1, 2:] - u[1:-1, 0:-2]) / (2 * dx))**2 -
                      2 * ((u[2:, 1:-1] - u[0:-2, 1:-1]) / (2 * dy) *
                           (v[1:-1, 2:] - v[1:-1, 0:-2]) / (2 * dx))-
                          ((v[2:, 1:-1] - v[0:-2, 1:-1]) / (2 * dy))**2))

    return b
```

The function `pressure_poisson` is also defined to help segregate the different rounds of calculations.  Note the presence of the pseudo-time variable `nit`.  This sub-iteration in the Poisson calculation helps ensure a divergence-free field.

```python
def pressure_poisson(p, dx, dy, b):
    pn = numpy.empty_like(p)
    pn = p.copy()

    for q in range(nit):
        pn = p.copy()
        p[1:-1, 1:-1] = (((pn[1:-1, 2:] + pn[1:-1, 0:-2]) * dy**2 +
                          (pn[2:, 1:-1] + pn[0:-2, 1:-1]) * dx**2) /
                          (2 * (dx**2 + dy**2)) -
                          dx**2 * dy**2 / (2 * (dx**2 + dy**2)) *
                          b[1:-1,1:-1])

        p[:, -1] = p[:, -2] # dp/dx = 0 at x = 2
        p[0, :] = p[1, :]   # dp/dy = 0 at y = 0
        p[:, 0] = p[:, 1]   # dp/dx = 0 at x = 0
        p[-1, :] = 0        # p = 0 at y = 2

    return p
```

Finally, the rest of the cavity flow equations are wrapped inside the function `cavity_flow`, allowing us to easily plot the results of the cavity flow solver for different lengths of time.

```python
def cavity_flow(nt, u, v, dt, dx, dy, p, rho, nu):
    un = numpy.empty_like(u)
    vn = numpy.empty_like(v)
    b = numpy.zeros((ny, nx))

    for n in range(nt):
        un = u.copy()
        vn = v.copy()

        b = build_up_b(b, rho, dt, u, v, dx, dy)
        p = pressure_poisson(p, dx, dy, b)

        u[1:-1, 1:-1] = (un[1:-1, 1:-1]-
                         un[1:-1, 1:-1] * dt / dx *
                        (un[1:-1, 1:-1] - un[1:-1, 0:-2]) -
                         vn[1:-1, 1:-1] * dt / dy *
                        (un[1:-1, 1:-1] - un[0:-2, 1:-1]) -
                         dt / (2 * rho * dx) * (p[1:-1, 2:] - p[1:-1, 0:-2]) +
                         nu * (dt / dx**2 *
                        (un[1:-1, 2:] - 2 * un[1:-1, 1:-1] + un[1:-1, 0:-2]) +
                         dt / dy**2 *
                        (un[2:, 1:-1] - 2 * un[1:-1, 1:-1] + un[0:-2, 1:-1])))

        v[1:-1,1:-1] = (vn[1:-1, 1:-1] -
                        un[1:-1, 1:-1] * dt / dx *
                       (vn[1:-1, 1:-1] - vn[1:-1, 0:-2]) -
                        vn[1:-1, 1:-1] * dt / dy *
                       (vn[1:-1, 1:-1] - vn[0:-2, 1:-1]) -
                        dt / (2 * rho * dy) * (p[2:, 1:-1] - p[0:-2, 1:-1]) +
                        nu * (dt / dx**2 *
                       (vn[1:-1, 2:] - 2 * vn[1:-1, 1:-1] + vn[1:-1, 0:-2]) +
                        dt / dy**2 *
                       (vn[2:, 1:-1] - 2 * vn[1:-1, 1:-1] + vn[0:-2, 1:-1])))

        u[0, :]  = 0
        u[:, 0]  = 0
        u[:, -1] = 0
        u[-1, :] = 1    # set velocity on cavity lid equal to 1
        v[0, :]  = 0
        v[-1, :] = 0
        v[:, 0]  = 0
        v[:, -1] = 0

    return u, v, p
```

Let's start with `nt = 100` and see what the solver gives us:

```python
u = numpy.zeros((ny, nx))
v = numpy.zeros((ny, nx))
p = numpy.zeros((ny, nx))
b = numpy.zeros((ny, nx))
nt = 100
u, v, p = cavity_flow(nt, u, v, dt, dx, dy, p, rho, nu)
```

```python
fig = pyplot.figure(figsize=(11,7), dpi=100)
# plotting the pressure field as a contour
pyplot.contourf(X, Y, p, alpha=0.5, cmap=cm.viridis)
pyplot.colorbar()
# plotting the pressure field outlines
pyplot.contour(X, Y, p, cmap=cm.viridis)
# plotting velocity field
pyplot.quiver(X[::2, ::2], Y[::2, ::2], u[::2, ::2], v[::2, ::2])
pyplot.xlabel('X')
pyplot.ylabel('Y');
```

![png](output_30_0.png)

You can see that two distinct pressure zones are forming and that the spiral pattern expected from lid-driven cavity flow is beginning to form.  Experiment with different values of `nt` to see how long the system takes to stabilize.

```python
u = numpy.zeros((ny, nx))
v = numpy.zeros((ny, nx))
p = numpy.zeros((ny, nx))
b = numpy.zeros((ny, nx))
nt = 700
u, v, p = cavity_flow(nt, u, v, dt, dx, dy, p, rho, nu)
```

```python
fig = pyplot.figure(figsize=(11, 7), dpi=100)
pyplot.contourf(X, Y, p, alpha=0.5, cmap=cm.viridis)
pyplot.colorbar()
pyplot.contour(X, Y, p, cmap=cm.viridis)
pyplot.quiver(X[::2, ::2], Y[::2, ::2], u[::2, ::2], v[::2, ::2])
pyplot.xlabel('X')
pyplot.ylabel('Y');
```

![png](output_33_0.png)

The quiver plot shows the magnitude of the velocity at the discrete points in the mesh grid we created.
(We're actually only showing half of the points because otherwise it's a bit of a mess.  The `X[::2, ::2]` syntax above is a convenient way to ask for every other point.)

Another way to visualize the flow in the cavity is to use a `streamplot`:

```python
fig = pyplot.figure(figsize=(11, 7), dpi=100)
pyplot.contourf(X, Y, p, alpha=0.5, cmap=cm.viridis)
pyplot.colorbar()
pyplot.contour(X, Y, p, cmap=cm.viridis)
pyplot.streamplot(X, Y, u, v)
pyplot.xlabel('X')
pyplot.ylabel('Y');
```

![png](output_35_0.png)

## Learn More

The interactive module **12 steps to Navier–Stokes** is one of several components of the Computational Fluid Dynamics class taught by Prof. Lorena A. Barba in Boston University between 2009 and 2013.

For a sample of what the other components of this class are, you can explore the **Resources** section of the Spring 2013 version of [the course's Piazza site](https://piazza.com/bu/spring2013/me702/resources).

***

```python
from IPython.core.display import HTML
def css_styling():
    styles = open("../styles/custom.css", "r").read()
    return HTML(styles)
css_styling()
```

<link href='http://fonts.googleapis.com/css?family=Fenix' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=Alegreya+Sans:100,300,400,500,700,800,900,100italic,300italic,400italic,500italic,700italic,800italic,900italic' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=Source+Code+Pro:300,400' rel='stylesheet' type='text/css'>
<style>
    @font-face {
        font-family: "Computer Modern";
        src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');
    }
    div.cell{
        width:800px;
        margin-left:16% !important;
        margin-right:auto;
    }
    h1 {
        font-family: 'Alegreya Sans', sans-serif;
    }
    h2 {
        font-family: 'Fenix', serif;
    }
    h3{
		font-family: 'Fenix', serif;
        margin-top:12px;
        margin-bottom: 3px;
       }
	h4{
		font-family: 'Fenix', serif;
       }
    h5 {
        font-family: 'Alegreya Sans', sans-serif;
    }
    div.text_cell_render{
        font-family: 'Alegreya Sans',Computer Modern, "Helvetica Neue", Arial, Helvetica, Geneva, sans-serif;
        line-height: 135%;
        font-size: 120%;
        width:600px;
        margin-left:auto;
        margin-right:auto;
    }
    .CodeMirror{
            font-family: "Source Code Pro";
			font-size: 90%;
    }
/*    .prompt{
        display: None;
    }*/
    .text_cell_render h1 {
        font-weight: 200;
        font-size: 50pt;
		line-height: 100%;
        color:#CD2305;
        margin-bottom: 0.5em;
        margin-top: 0.5em;
        display: block;
    }
    .text_cell_render h5 {
        font-weight: 300;
        font-size: 16pt;
        color: #CD2305;
        font-style: italic;
        margin-bottom: .5em;
        margin-top: 0.5em;
        display: block;
    }

    .warning{
        color: rgb( 240, 20, 20 )
        }
</style>
<script>
    MathJax.Hub.Config({
                        TeX: {
                           extensions: ["AMSmath.js"]
                           },
                tex2jax: {
                    inlineMath: [ ['$','$'], ["\\(","\\)"] ],
                    displayMath: [ ['$$','$$'], ["\\[","\\]"] ]
                },
                displayAlign: 'center', // Change this to 'center' to center equations.
                "HTML-CSS": {
                    styles: {'.MathJax_Display': {"margin": 4}}
                }
        });
</script>

> (The cell above executes the style for this notebook.)

max_stars_repo_name: mkitti/CFDJulia

max_issues_repo_name: mkitti/CFDJulia

max_forks_repo_name: mkitti/CFDJulia


text:
function get_all_paths_in_tree!(
    tree    :: Union{LeafOrNodeCausalH{S}, LeafOrNodeCausalNH{S}},
    result  = [],
    path    = []
    ) where {S}

    push!(path, tree.featid)
    test = true

    if !is_leaf(tree.left)
        get_all_paths_in_tree!(tree.left, result, path)
        test = false
    end
    if !is_leaf(tree.right)
        get_all_paths_in_tree!(tree.right, result, path)
        test = false
    end

    if test
        push!(result, copy(path))
    end

    pop!(path)
    return result
end

function get_occurence_frequencies(forest :: EnsembleCausal{S}) where {S}
    n_trees = length(forest.trees)
    dico = Dict{Set{Int}, Float64}()
    full_set = Set(1:size(forest.X, 2))
    empty_set = Set()

    for i in 1:n_trees
        tree = forest.trees[i].tree
        paths_list = get_all_paths_in_tree!(tree)
        for j in 1:length(paths_list)
            full_path = paths_list[j]
            for k in 1:length(full_path)
                set_u = Set(full_path[1:k])
                if set_u == full_set
                    break
                end
                if haskey(dico, set_u)
                    dico[set_u] += 1
                else
                    dico[set_u] = 1
                end
            end
        end
    end

    delete!(dico, full_set)
    delete!(dico, empty_set)

    nb_sets = sum(values(dico))

    for key in keys(dico)
        dico[key] = dico[key] / nb_sets
    end

    return dico

end

function sample_U(dico :: Dict{Set{Int}, Float64}, n :: Int)
    return StatsBase.sample(collect(keys(dico)), StatsBase.pweights(collect(values(dico))), n)
end