text:
using MathProgBase, Gurobi
using Distributions, DataFrames

"""
Computes the optimal q value via the cross validation procedure.
Args:
    data (array):       vector of binary outcomes.
    tr_ind (array):     Indices for training sample.
    test_ind (array):   Indices for validation sample.
    focus_ind (array):  Indices for focus covariates.
    aux_ind (array):    Indices for auxiliary covariates.
    beta0 (int):        Coefficient taking value either 1 or -1 to normalize the scale for the first
                        first covariate in x_foc.
    q (int):            Cardinality constraint for the covariate selection.
    T (int):            Time limit specified for the MIO solver.
    bnd (array):        (((k-1)+d) by 2) matrix where the first and second columns respectively
                        store the lower and upper bounds of the unknown coefficients. First (k-1)
                        rows correspond to the bounds of the focused covariates excluding the first one.
                        Remaining d rows correspond to the bounds of the auxiliary covariates.
    mio (int) :         MIO formulation of the best subset maximum score problem
                        set mio = 1 for Method 1 of the paper
                        set mio = 2 for Method 2 of the paper
    tau (int):          Tuning parameter for enlarging the estimated bounds.
Yields:
    best_q:             Best q value.
    bhat (array):       Maximum score estimates for the unknown coefficients.
    score (array):      Value of maximum score objective function.
    gap   (array):      MIO optimization gap value in case of early termination
                        gap = 0 ==> optimal solution is found within the time limit.
    rtime (array):      Time used by the MIO solver in the estimation procedure.
    ncount (array):     Node count.
"""
function CVMaxScore(tr_ind::AbstractArray,test_ind::AbstractArray,data::AbstractArray,focus_ind::AbstractArray,aux_ind::AbstractArray,beta0::Integer,q_range::AbstractMatrix,T::Real,tol::Real,bnd::AbstractArray,mio::Integer)
    q_num = len(q_range)
    fold=size(tr_ind)[1]
    score=zeros(fold,q_num)
    gap=zeros(fold,q_num)
    rtime=zeros(fold,q_num)
    ncount=zeros(fold,q_num)
    bhat=zeros(len(focus_ind)+len(aux_ind)-1,fold,q_num)
    val_score=zeros(q_num,1)

    for q=1:q_num
        for i=1:fold
            println("(q,fold) : ",q_range[q]," ",i+1)
            y=data[tr_ind[:,i],1]
            datax=data[tr_ind[:,i],2:size(data)[2]]

            bhat[:,i,q],score[i,q],gap[i,q],rtime[i,q],ncount[i,q]  = MaxScore(y,datax[:,focus_ind],datax[:,aux_ind],beta0,q_range[q],T,tol[q],bnd,mio)

            y_v=data[test_ind[:,i],1]
            datax_v=data[test_ind[:,i],2:size(data)[2]]
            val_score[q] = val_score[q] + mean(Int.(y_v .== (datax_v*[beta0;bhat[:,i,q]].>=0)))
        end
        val_score[q]=val_score[q]/fold
    end
    best_q = maximum(val_score)
    return best_q, bhat,score,gap,rtime,ncount
end

export CVMaxScore


text:
"""
    print the err text
"""
printerr(io::IO, err) = print(io, sprint(showerror, err, catch_backtrace()))
printerr(err) = printerr(stdout, err)

ext: jl


text:
abstract type Bound{T} end

Base.show(io::IO, bound::Bound) = summary(io, bound)

@inline getT(bound::Bound{T}) where T = T

@inline g(bound::Bound{T}) where T = zero(T)

@inline ∂g∂pos(bound::Bound{T}) where T = @SVector zeros(T, 6)
@inline ∂g∂vel(bound::Bound{T}) where T = @SVector zeros(T, 6)
@inline schurf(bound::Bound{T}) where T = @SVector zeros(T, 6)
@inline schurD(bound::Bound{T}) where T = @SMatrix zeros(T, 6, 6)


text:
#=
read_nrn4_and_extract_voxels.jl:
- Julia version: 1.1
- Author: Jirik
- Date: 2019-03-05
=#

# include("../src/LarSurf.jl")
using Revise
using Plasm
# using LinearAlgebraicRepresentation
# Lar = LinearAlgebraicRepresentation
using LarSurf

pth = Pio3d.datasets_join_path("medical/orig/sample-data/nrn4.pklz")
datap = Pio3d.read3d(pth)

data3d = datap["data3d"]
voxelsize_mm = datap["voxelsize_mm"]
threshold = 4000

vertices, nvoxels = LarSurf.vertices_and_count_voxels_as_the_square_face_preprocessing(data3d, voxelsize_mm, threshold)
voxels = LarSurf.create_cuboid_voxels(data3d, nvoxels, threshold)

voxels_linear01 = LarSurf.ind_to_sparse(hcat(voxels), prod(size(data3d)), 1)


text:
module TestExamples2DEulerMulticomponent

using Test
using Trixi

include("test_trixi.jl")

# pathof(Trixi) returns /path/to/Trixi/src/Trixi.jl, dirname gives the parent directory
EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "tree_2d_dgsem")

@testset "Compressible Euler Multicomponent" begin
  # NOTE: Some of the L2/Linf errors are comparably large. This is due to the fact that some of the
  #       simulations are set up with dimensional states. For example, the reference pressure in SI
  #       units is 101325 Pa, i.e., pressure has values of O(10^5)

  @trixi_testset "elixir_eulermulti_shock_bubble.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_shock_bubble.jl"),
      l2   = [73.78467629094177, 0.9174752929795251, 57942.83587826468, 0.1828847253029943, 0.011127037850925347],
      linf = [196.81051991521073, 7.8456811648529605, 158891.88930113698, 0.811379581519794, 0.08011973559187913],
      tspan = (0.0, 0.001))
  end

  @trixi_testset "elixir_eulermulti_ec.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_ec.jl"),
      l2   = [0.050182236154087095, 0.050189894464434635, 0.2258715597305131, 0.06175171559771687],
      linf = [0.3108124923284472, 0.3107380389947733, 1.054035804988521, 0.29347582879608936])
  end

  @trixi_testset "elixir_eulermulti_es.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_es.jl"),
      l2   = [0.0496546258404055, 0.04965550099933263, 0.22425206549856372, 0.004087155041747821, 0.008174310083495642, 0.016348620166991283, 0.032697240333982566],
      linf = [0.2488251110766228, 0.24832493304479406, 0.9310354690058298, 0.017452870465607374, 0.03490574093121475, 0.0698114818624295, 0.139622963724859])
  end

  @trixi_testset "elixir_eulermulti_convergence_ec.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_convergence_ec.jl"),
      l2   = [0.00012290225488326508, 0.00012290225488321876, 0.00018867397906337653, 4.8542321753649044e-5, 9.708464350729809e-5],
      linf = [0.0006722819239133315, 0.0006722819239128874, 0.0012662292789555885, 0.0002843844182700561, 0.0005687688365401122])
  end

  @trixi_testset "elixir_eulermulti_convergence_es.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_convergence_es.jl"),
      l2   = [2.2661773867001696e-6, 2.266177386666318e-6, 6.593514692980009e-6, 8.836308667348217e-7, 1.7672617334696433e-6],
      linf = [1.4713170997993075e-5, 1.4713170997104896e-5, 5.115618808515521e-5, 5.3639516094383666e-6, 1.0727903218876733e-5])
  end

  @trixi_testset "elixir_eulermulti_convergence_es.jl with flux_chandrashekar" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_convergence_es.jl"),
      l2   = [1.8621737639352465e-6, 1.862173764098385e-6, 5.942585713809631e-6, 6.216263279534722e-7, 1.2432526559069443e-6],
      linf = [1.6235495582606063e-5, 1.6235495576388814e-5, 5.854523678827661e-5, 5.790274858807898e-6, 1.1580549717615796e-5],
      volume_flux = flux_chandrashekar)
  end
end

end # module


text:
export Signature, MemorySettings, MemoryWrapper, Sig

Sig(name::Symbol) = Signature(name)

struct Signature
    name::Symbol
    Base.string(x::Signature) = "Function:$(x.name)"
end

"Contains the settings for a function to utilize in compilation"
struct MemorySettings
    storageType::Type
    inputType::Type
    outputType::Type
    device::Device
    deviceArray::Type

    MemorySettings() = MemorySettings(Double, Double, Double, CPUDev, Array)

    MemorySettings(storageType::Type, inputType::Type, outputType::Type, device::Device, deviceArray::Type) =
        MemorySettings(storageType, inputType, outputType, device, deviceArray)

    MemorySettings(settings::MemorySettings; storageType=settings.storageType,
                   inputType=settings.inputType, outputType=settings.outputType,
                   device=settings.device) =
            MemorySettings(storageType, inputType, outputType, device, arraytype(settings.device))
end

"Wrapper designed to take in memory settings and return some form of object with those settings builtin"
struct MemoryWrapper{ReturnType}
    sig::Signature
    fun::Function

    MemoryWrapper{R}(sig::Signature, fun::Function) where R = new{R}(sig, fun)

    function (wrapper::MemoryWrapper)(settings::MemorySettings)
        return wrapper.fun(settings)
    end
end


text:
include("libavutil_h.jl")

#include("audio_fifo.jl")
#include("avutil.jl")
#include("channel_layout.jl")
#include("dict.jl")
#include("fifo.jl")
#include("file.jl")
#include("imgutils.jl")
include("log.jl")
include("mem.jl")
include("opt.jl")
include("pixdesc.jl")
#include("rational.jl")
#include("samplefmt.jl")
#include("xtea.jl")


text:
Pkg.add("MathematicalSystems");
repo = LibGit2.GitRepo(joinpath(Pkg.dir(), "MathematicalSystems"));
LibGit2.branch!(repo, "Pin_3451c", "3451ce51bfa15078279610d5d8c3eb62f23ca3ea");

Pkg.add("LazySets");
repo = LibGit2.GitRepo(joinpath(Pkg.dir(), "LazySets"));
LibGit2.branch!(repo, "Pin_08150", "08150a1307f9c25c5e02fe7fe9a74e102049bd6e");

Pkg.add("Reexport");

Pkg.clone("https://github.com/JuliaReach/SX.jl.git");
repo = LibGit2.GitRepo(joinpath(Pkg.dir(), "SX"));
LibGit2.branch!(repo, "Pin_0bee1", "0bee143c44c1e6d9210e6754910f03bca7ad1a11");

try
    Pkg.clone("https://github.com/JuliaReach/Reachability.jl.git");
catch
    repo = LibGit2.GitRepo(joinpath(Pkg.dir(), "Reachability"));
    LibGit2.branch!(repo, "Pin_73f9e", "73f9ee92342e4d46418589d34200c5682bcb8609");
end

Pkg.clone("https://github.com/JuliaReach/ARCH2018_RE.git")


text:
module SimTools

    import Statistics
    import BSON
    import Serialization: serialize, deserialize
    import ProgressMeter: ProgressThresh, Progress, next!, finish!, update!

    include("bisection_search.jl")
    include("gradient_descent.jl")
    include("is_stationary.jl")
    include("OnDiskMonitor.jl")
    include("utils.jl")

end


text:
# RELOG: Reverse Logistics Optimization
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.

module RELOG

include("instance/structs.jl")

include("graph/structs.jl")

include("graph/build.jl")
include("graph/csv.jl")
include("instance/compress.jl")
include("instance/geodb.jl")
include("instance/parse.jl")
include("instance/validate.jl")
include("model/build.jl")
include("model/getsol.jl")
include("model/solve.jl")
include("model/resolve.jl")
include("reports/plant_emissions.jl")
include("reports/plant_outputs.jl")
include("reports/plants.jl")
include("reports/products.jl")
include("reports/tr_emissions.jl")
include("reports/tr.jl")
include("reports/write.jl")
end

max_stars_repo_name: ANL-CEEESA/RELOG

max_issues_repo_name: ANL-CEEESA/RELOG

max_forks_repo_name: ANL-CEEESA/RELOG


text:
# Generates a series of circles that will form the surface
function generate_circles{T<:AbstractFloat}(
        center::Tuple{T,T},
        radius::T,
        n_points::Int64,   # Points on a contour
        n_contours::Int64  # Number of contours
    )
    n_contours >= 2 || error("NO")
    X = Array(T, n_points*n_contours, 2)
    angle_step = 2π/n_points
    radius_step = radius/(n_contours-1)
    for j = 1:n_contours
        r = (n_contours-j) * radius_step
        for i = 1:n_points
            idx = (j-1)*n_points + i
            θ = (i-1)*angle_step
            X[idx, 1] = center[1] + r * cos(θ)
            X[idx, 2] = center[2] + r * sin(θ)
        end
    end
    X
end

# Returns an n_rows+1 by n_cols Array with the last row being a copy of the first row
# This completes the surface (or else you're missing a slice of pie)
function circular_reshape{T<:AbstractFloat}(A::Array{T}, n_rows::Int64, n_cols::Int64)
    length(A) == n_rows*n_cols || throw(DimensionMismatch("Dimensions do not conform"))
    T[A[(j-1)*n_rows + mod(i-1, n_rows) + 1] for i = 1:n_rows+1, j = 1:n_cols]
end

function coordinate_matrices{T<:AbstractFloat}(W::Matrix{T}, n_points::Int64, n_contours::Int64)
    size(W,2) == 3 || error("Matrix must be n by 3")
    X = circular_reshape(W[:,1], n_points, n_contours)
    Y = circular_reshape(W[:,2], n_points, n_contours)
    Z = circular_reshape(W[:,3], n_points, n_contours)
    (X, Y, Z)
end

function generate_point_cluster{T<:AbstractFloat}(n, center::Tuple{T,T}, radius::T)
    u = radius * sqrt(rand(n))
    v = 2π*rand(n)
    (u .* hcat(cos(v), sin(v))) .+ [center[1] center[2]]
end

function generate_point_ring{T<:AbstractFloat}(n, center::Tuple{T,T}, band_center::T, width::T)
    u = (rand(n) * width) .+ (band_center - width/2)
    v = rand(n) * 2band_center
    r = Array(T, n)
    for i = 1:n
        r[i] = v[i] < u[i] ? u[i] : 2band_center - u[i]
    end
    θ = 2π*rand(n)
    r .* hcat(cos(θ), sin(θ)) .+ [center[1] center[2]]
end

# Performs a form of KPCA for 3 dimensions
function kpca_3d{T<:Base.LinAlg.BlasReal}(K::Matrix{T})
    Λ, V = LAPACK.syev!('V', 'U', centerkernelmatrix!(copy(K)))
    W = V[:,end:-1:end-2] * 10  # eigenvalues are in ascending order
end

using MLKernels, PyPlot

#===================================================================================================
  Kernel Trick Visualization - Linearly Separable Classes
===================================================================================================#

n_c1 = 60
n_c2 = 180

cluster = 0.55
ring    = (0.725, 0.975)  # radius 1 & 2

xy_mid = (cluster^2 + ring[1]^2)/2
z_lim  = ring[2] * cluster * sqrt(2)

# Class Sampling
srand(1234)
c1 = generate_point_cluster(n_c1, (0.0, 0.0), cluster)
c2 = generate_point_ring(n_c2, (0.0, 0.0), sum(ring)/2, ring[2]-ring[1])

# Plot points in original feature space
PyPlot.figure("FeatureSpace")
PyPlot.scatter(vec(c1[:,1]), vec(c1[:,2]), c="r", marker="D")
PyPlot.scatter(vec(c2[:,1]), vec(c2[:,2]), c="m", marker="o")
PyPlot.xlim([-1.05,1.05])
PyPlot.ylim([-1.05,1.05])
PyPlot.title("Original Feature Space")
PyPlot.xlabel("Feature 1")
PyPlot.ylabel("Feature 2")

# Explicit Feature Map (Polynomial with d = 2, c = 0)
#   (x,y) -> (x^2, y^2, √2*x*y)
ϕ(x,y) = [x^2 y^2 √2*x*y]
c1_ϕ = vcat([ϕ(c1[i,:]...) for i = 1:size(c1,1)]...)
c2_ϕ = vcat([ϕ(c2[i,:]...) for i = 1:size(c2,1)]...)

# Plot points in kernel Hilbert space
PyPlot.figure("HilbertSpace")
PyPlot.scatter3D(vec(c1_ϕ[:,1]), vec(c1_ϕ[:,2]), vec(c1_ϕ[:,3]), c="r", marker="D")
PyPlot.scatter3D(vec(c2_ϕ[:,1]), vec(c2_ϕ[:,2]), vec(c2_ϕ[:,3]), c="m", marker="o")
PyPlot.xlim([-0.05,1.05])
PyPlot.ylim([-0.05,1.05])
PyPlot.zlim([-z_lim - 0.05, z_lim + 0.05])
PyPlot.title("Kernel Hilbert Space")
PyPlot.xlabel("Dimension 1")
PyPlot.ylabel("Dimension 2")
PyPlot.zlabel("Dimension 3")

# Add Hyperplane
X = [ 0.0    0.0   ;  xy_mid xy_mid]
Y = [ xy_mid xy_mid;  0.0    0.0   ]
Z = [-z_lim  z_lim ; -z_lim  z_lim ]
PyPlot.plot_surface(X, Y, Z, rstride=1, cstride=1, alpha=0.5, color="b")

# Add Circular Intersection
v = linspace(0, 2π, 360)
X = xy_mid*(cos(v).+1)/2
Y = xy_mid .- X
Z = √2 * sqrt(X) .* sqrt(Y) .* sign(sin(v))
PyPlot.plot(X, Y, Z, color="k", ls="-")

# Kernel Geometry
U = generate_circles((0.0,0.0),ring[2],60,2)
W = vcat([ϕ(U[i,:]...) for i = 1:size(U,1)]...)
X, Y, Z = coordinate_matrices(W, 60, 2)
PyPlot.figure("KernelGeometry")
PyPlot.plot_surface(X, Y, Z, rstride=1, cstride=1, alpha=0.5, linewidth=0,
                    cmap=ColorMap("coolwarm"))
PyPlot.xlim([-0.05,1.05])
PyPlot.ylim([-0.05,1.05])
PyPlot.zlim([-z_lim - 0.05, z_lim + 0.05])
PyPlot.title("Kernel Geometry")
PyPlot.xlabel("Dimension 1")
PyPlot.ylabel("Dimension 2")
PyPlot.zlabel("Dimension 3")

# Add Circular Intersection
v = linspace(0, 2π, 360)
X = xy_mid*(cos(v).+1)/2
Y = xy_mid .- X
Z = √2 * sqrt(X) .* sqrt(Y) .* sign(sin(v))
PyPlot.plot(X, Y, Z, color="k", ls="-")

# Add Intersection to original space
PyPlot.figure("FeatureSpace")
PyPlot.plot(sqrt(xy_mid)*sin(v), sqrt(xy_mid)*cos(v), c="b", ls="-")

#===================================================================================================
  Kernel Trick Visualization
===================================================================================================#

for (κ, title, n_points, n_contours, center, radius) in (
        (SquaredDistanceKernel(0.5), "Squared-Distance Kernel", 40, 15, (0.0, 0.0), 1.0),
        (SineSquaredKernel(), "Sine-Squared Kernel", 40, 15, (0.0, 0.0), 0.5),
        (ChiSquaredKernel(), "Chi-Squared Kernel", 40, 15, (0.5, 0.5), 0.5),
        (GaussianKernel(), "Gaussian Kernel", 40, 15, (0.0, 0.0), 1.0),
        (LaplacianKernel(), "Laplacian Kernel", 40, 15, (0.0, 0.0), 1.0),
        (PeriodicKernel(), "Periodic Kernel", 40, 15, (0.0, 0.0), 0.5),
        (RationalQuadraticKernel(), "Rational-Quadratic Kernel", 40, 15, (0.0, 0.0), 1.0),
        (MaternKernel(), "Matern Kernel", 40, 15, (0.0, 0.0), 1.0),
        (PolynomialKernel(), "Polynomial Kernel", 40, 15, (0.0, 0.0), 1.0),
        (SigmoidKernel(), "Sigmoid Kernel", 40, 15, (0.0, 0.0), 1.0)
    )
    U = generate_circles(center, radius, n_points, n_contours)
    W = kpca_3d((isnegdef(κ) ? -1 : 1) * kernelmatrix(κ, U))
    X, Y, Z = coordinate_matrices(W, n_points, n_contours)
    PyPlot.figure(title)
    PyPlot.plot_surface(X, Y, Z, rstride=1, cstride=1, alpha=1.0, linewidth=0,
                        cmap=ColorMap("coolwarm"))
    PyPlot.title(title)
    PyPlot.xlabel("Component 1")
    PyPlot.ylabel("Component 2")
    PyPlot.zlabel("Component 3")
end


text:
module AWSExceptions

using HTTP
using JSON
using XMLDict
using XMLDict: XMLDictElement

export AWSException, ProtocolNotDefined, InvalidFileName, NoCredentials

struct ProtocolNotDefined <: Exception
    message::String
end
Base.show(io::IO, e::ProtocolNotDefined) = println(io, e.message)

struct InvalidFileName <: Exception
    message::String
end
Base.show(io::IO, e::InvalidFileName) = println(io, e.message)

struct NoCredentials <: Exception
    message::String
end
Base.show(io::IO, e::NoCredentials) = println(io, e.message)

struct AWSException <: Exception
    code::String
    message::String
    info::Union{XMLDictElement, Dict, String, Nothing}
    cause::HTTP.StatusError
end
function Base.show(io::IO, e::AWSException)
    message = isempty(e.message) ? "" : (" -- " * e.message)
    println(io, "AWSException: ", string(e.code, message, "\n", e.cause))
end

http_message(e::HTTP.StatusError) = String(copy(e.response.body))
http_status(e::HTTP.StatusError) = e.status
content_type(e::HTTP.StatusError) = HTTP.header(e.response, "Content-Type")
is_valid_xml_string(str) = startswith(str, '<')

function AWSException(e::HTTP.StatusError)
    code = string(http_status(e))
    message = "AWSException"
    info = Dict{String, Dict}()

    # Extract API error code from Lambda-style JSON error message...
    if endswith(content_type(e), "json")
        info = JSON.parse(http_message(e))
    end

    # Extract API error code from JSON error message...
    if occursin(r"^application/x-amz-json-1\.[01]$", content_type(e))
        info = JSON.parse(http_message(e))
        if haskey(info, "__type")
            code = rsplit(info["__type"], '#', limit=2)[end]
        end
    end

    # Extract API error code from XML error message...
    error_content_types = ["", "application/xml", "text/xml"]
    if content_type(e) in error_content_types && is_valid_xml_string(http_message(e))
        info = parse_xml(http_message(e))
    end

    # There are times when Errors or Error are returned back
    info = get(info, "Errors", info)
    info = get(info, "Error", info)

    code = get(info, "Code", code)

    # There are also times when the response back is (M|m)essage
    message = get(info, "Message", message)
    message = get(info, "message", message)

    return AWSException(code, message, info, e)
end

end


text:
using Test
include("../src/UTriMatrix.jl")

@testset "UTriMatrix.jl" begin
    n = 4
    nentries = (n*(n+1)) ÷ 2
    m = UTriMatrix(n, BigInt(1))
    @test size(m) == (n,n)
    @test length(m) == nentries

    for i in 1:n
        for j in i:n
            m[i,j] = i+j
        end
    end

    @test m[1,1] == 2
    @test m[1,2] == 3
    @test m[2,3] == 5
    @test m[4,4] == 8
    @test_throws BoundsError m[2,1]
end

@testset "BitsUTriMatrix.jl" begin
    n = 4
    nentries = (n*(n+1)) ÷ 2
    m = BitsUTriMatrix(n, 1)
    @test size(m) == (n,n)
    @test length(m) == nentries

    for i in 1:n
        for j in i:n
            m[i,j] = i+j
        end
    end

    @test m[1,1] == 2
    @test m[1,2] == 3
    @test m[2,3] == 5
    @test m[4,4] == 8
    @test_throws BoundsError m[2,1]
end


text:
using DSGE, ModelConstructors, Dates, CSV, DataFrames, Plots, Test, HDF5
include("tvcred_parameterize.jl")

regenerate_reference_forecasts = false

# Forecast settings
fcast_date = DSGE.quartertodate("2020-Q4")
date_fcast_end  = iterate_quarters(fcast_date, 60)

# Flexible AIT Rule
imperfect_cred_flexait = .33
imperfect_cred_taylor  = 1. - imperfect_cred_flexait
φ_π                    = 10.
φ_y                    = 6.
Thalf                  = 10.
ρ_smooth               = 0.7
pgap_init              = 0.125
ygap_init              = 12.
pgap_ygap_init_date    = Date(2020, 6, 30)
start_zlb_date         = Date(2021, 3, 31)
end_zlb_date           = Date(2022, 12, 31)

# additional settings to implement Flexible AIT rule
flexait_custom = Dict{Symbol, Setting}(:add_initialize_pgap_ygap_pseudoobs => Setting(:add_initialize_pgap_ygap_pseudoobs, true),
                                       :add_pgap => Setting(:add_pgap, true), :add_ygap => Setting(:add_ygap, true),
                                       :add_altpolicy_pgap => Setting(:add_altpolicy_pgap, true),
                                       :add_altpolicy_ygap => Setting(:add_altpolicy_ygap, true),
                                       :forecast_horizons => Setting(:forecast_horizons,
                                                                     subtract_quarters(date_fcast_end, fcast_date)),
                                       :forecast_smoother => Setting(:forecast_smoother, :durbin_koopman),
                                       :contemporaneous_and_proportional_antshocks =>
                                       Setting(:contemporaneous_and_proportional_antshocks, Symbol[:biidc]),
                                       :antshocks =>
                                       Setting(:antshocks, Dict{Symbol, Int}(:biidc => 1,
                                                                             :φ => 1, :ziid => 1)),
                                       :ant_eq_mapping =>
                                       Setting(:ant_eq_mapping, Dict{Symbol, Symbol}(:biidc => :biidc,
                                                                                :φ => :φ,
                                                                                :ziid => :ziid)),
                                       :ant_eq_E_mapping =>
                                       Setting(:ant_eq_E_mapping, Dict{Symbol, Symbol}(:φ => :Eφ)),
                                       :proportional_antshocks =>
                                       Setting(:proportional_antshocks, [:biidc, :φ, :ziid]),
                                       :n_anticipated_obs_gdp => Setting(:n_anticipated_obs_gdp, 1),
                                       :add_anticipated_obs_gdp => Setting(:add_anticipated_obs_gdp, true),
                                        :meas_err_anticipated_obs_gdp =>
                                        Setting(:meas_err_anticipated_obs_gdp, 1.),
                                       :flexible_ait_policy_change =>
                                       Setting(:flexible_ait_policy_change, false))

# Initialize model object
m = Model1002("ss59"; custom_settings = flexait_custom)
usual_model_settings!(m, "200001", cdvt = "200001", fcast_date = fcast_date)
m <= Setting(:time_varying_trends, true)
get_setting(m, :regime_dates)[5] = Date(2020, 12, 31)
setup_regime_switching_inds!(m, cond_type = :full)

m10 = Model1002("ss10") # for help initializing parameters of m
θ10 = h5read(joinpath(dirname(@__FILE__), "..", "reference", "tvcred_reference_forecast.h5"), "para")
DSGE.update!(m10, θ10)
update_vals = zeros(Float64, length(m.parameters))
for k in map(x-> x.key, m10.parameters)
    update_vals[m.keys[k]] = m10[k].value
end
DSGE.update!(m, update_vals) # make sure parameters match Model 1002 ss10, w/zeros for other parameters

# Parameterize
tvcred_parameterize!(m)

# Load data
df = CSV.read(joinpath(dirname(@__FILE__), "..", "reference", "timevaryingcred_df.csv"), DataFrame)
df[!, :obs_pgap] .= NaN # NaN out just to be safe
df[!, :obs_ygap] .= NaN
ind_init = findfirst(df[!, :date] .== pgap_ygap_init_date)
df[end, :obs_gdp] = df[end - 1, :obs_gdp1] # To match reference forecast

# Set up Flexible AIT with a temporary ZLB. Note that df will be altered
θ = (φ_π = φ_π, φ_y = φ_y, Thalf = Thalf, pgap = pgap_init,
     ygap = ygap_init, ρ_smooth = ρ_smooth,
     cred = imperfect_cred_flexait,
     historical_policy = DSGE.taylor_rule()) # NamedTuple with parameters of the Flexible AIT policy
df[ind_init, :obs_pgap] = -θ[:pgap] # initialize pgap and ygap in the date
df[ind_init, :obs_ygap] = -θ[:ygap] # pgap_ygap_init_date, default to 2020:Q2

m <= Setting(:flexible_ait_φ_π, θ[:φ_π])
m <= Setting(:flexible_ait_φ_y, θ[:φ_y])
m <= Setting(:ait_Thalf, θ[:Thalf])
m <= Setting(:gdp_Thalf, θ[:Thalf])
m <= Setting(:pgap_value, θ[:pgap])
m <= Setting(:pgap_type, :flexible_ait)
m <= Setting(:ygap_value, θ[:ygap])
m <= Setting(:ygap_type, :flexible_ait)
m <= Setting(:flexible_ait_ρ_smooth, θ[:ρ_smooth])
m <= Setting(:alternative_policies, AltPolicy[θ[:historical_policy]])
m <= Setting(:skip_altpolicy_state_init, true)

## Set up temporary ZLB
m <= Setting(:gensys2, true)

# Set up credibility for ZLB and alternative policy, if applicable
m <= Setting(:uncertain_temporary_altpolicy, true)
m <= Setting(:temporary_altpolicy_names, [:zero_rate])
m <= Setting(:zero_rate_zlb_value, 0.)
m <= Setting(:uncertain_altpolicy, true)
gensys2_first_regime = get_setting(m, :n_regimes) + 1
get_setting(m, :regime_dates)[gensys2_first_regime] = start_zlb_date
n_zlb_reg = DSGE.subtract_quarters(end_zlb_date, start_zlb_date) + 1

# Set up regime_eqcond_info
m <= Setting(:replace_eqcond, true)
reg_dates = deepcopy(get_setting(m, :regime_dates))
regime_eqcond_info = Dict{Int, DSGE.EqcondEntry}()
for (regind, date) in zip(gensys2_first_regime:(n_zlb_reg - 1 + gensys2_first_regime), # See comments starting at line 57
                          DSGE.quarter_range(reg_dates[gensys2_first_regime],
                                             DSGE.iterate_quarters(reg_dates[gensys2_first_regime], n_zlb_reg - 1)))
    reg_dates[regind] = date
    regime_eqcond_info[regind] = DSGE.EqcondEntry(DSGE.zero_rate(), [0., 1.])
end
reg_dates[n_zlb_reg + gensys2_first_regime] = DSGE.iterate_quarters(reg_dates[gensys2_first_regime], n_zlb_reg)
regime_eqcond_info[n_zlb_reg + gensys2_first_regime] = DSGE.EqcondEntry(DSGE.flexible_ait(), [0., 1.])
nreg0 = length(reg_dates)
m <= Setting(:regime_dates,             reg_dates)
m <= Setting(:regime_eqcond_info, regime_eqcond_info)

# Set up regime indices and extra regimes for parameters
setup_regime_switching_inds!(m; cond_type = :full)
set_regime_vals_fnct(m, get_setting(m, :n_regimes))

# Set up TVIS information set
m <= Setting(:tvis_information_set, vcat([i:i for i in 1:(gensys2_first_regime - 1)],
                                         [i:get_setting(m, :n_regimes) for i in
                                          gensys2_first_regime:get_setting(m, :n_regimes)]))

# output_vars = [:histobs, :histpseudo, :forecastobs, :forecastpseudo, :forecaststates]
output_vars = [:forecastobs, :forecastpseudo]
modal_params = map(x -> x.value, m.parameters)
@assert !haskey(m.settings, :temporary_altpolicy_length)
outp0 = DSGE.forecast_one_draw(m, :mode, :full, output_vars, modal_params, df,
                                     regime_switching = true, n_regimes = get_setting(m, :n_regimes))
for i in keys(get_setting(m, :regime_eqcond_info))
    get_setting(m, :regime_eqcond_info)[i].weights = [θ[:cred], 1. - θ[:cred]]
end

# Compare fixed vs. time-varying credibility
outp33 = DSGE.forecast_one_draw(m, :mode, :full, output_vars, modal_params, df,
                                regime_switching = true, n_regimes = get_setting(m, :n_regimes))
sysp33 = compute_system(m; tvis = true)
for i in nreg0:(nreg0 + 15)
    reg_dates[i] = DSGE.iterate_quarters(reg_dates[nreg0], i - nreg0)
    regime_eqcond_info[i] = DSGE.EqcondEntry(DSGE.flexible_ait(), [.33, 1-.33])
end
m <= Setting(:regime_dates,             reg_dates)
m <= Setting(:regime_eqcond_info, regime_eqcond_info)
setup_regime_switching_inds!(m; cond_type = :full)
m <= Setting(:tvis_information_set, vcat([i:i for i in 1:(gensys2_first_regime - 1)],
                                         [i:get_setting(m, :n_regimes) for i in
                                          gensys2_first_regime:get_setting(m, :n_regimes)]))
set_regime_vals_fnct(m, get_setting(m, :n_regimes))
m <= Setting(:temporary_altpolicy_length, n_zlb_reg)
for i in keys(get_setting(m, :regime_eqcond_info))
    get_setting(m, :regime_eqcond_info)[i].weights = [.33, 1-.33]
end
sysp33_tv = compute_system(m; tvis = true)
outp33_tv = DSGE.forecast_one_draw(m, :mode, :full, output_vars, modal_params, df,
                                   regime_switching = true, n_regimes = get_setting(m, :n_regimes))

@testset "Compare Fixed to Time-Varying Credibility" begin
    for reg in 1:n_regimes(sysp33)
        @test sysp33[reg, :TTT] ≈ sysp33_tv[reg, :TTT]
        @test sysp33[reg, :RRR] ≈ sysp33_tv[reg, :RRR]
        @test sysp33[reg, :CCC] ≈ sysp33_tv[reg, :CCC]
        @test sysp33[reg, :QQ] ≈ sysp33_tv[reg, :QQ]
        @test sysp33[reg, :ZZ] ≈ sysp33_tv[reg, :ZZ]
        @test sysp33[reg, :DD] ≈ sysp33_tv[reg, :DD]
    end
    for k in keys(outp33)
        @test outp33[k] ≈ outp33_tv[k]
    end

    if !regenerate_reference_forecasts
        testfcast = h5read(joinpath(dirname(@__FILE__), "..", "reference", "tvcred_reference_meas_fix_forecast.h5"), "forecastobs")
        inds = vcat(1:9, 12:13, 20:21)  # ignore k-periods ahead observables, reference data generated when a bug existed
        @test maximum(abs.(testfcast[inds, :] -
                           outp33[:forecastobs][inds, :])) < 1.1e-3#5e-4 # some numerical differences b/c fixes to calculations
    end

end

# Compare perfect credibility to permanent case
for i in keys(get_setting(m, :regime_eqcond_info))
    get_setting(m, :regime_eqcond_info)[i].weights = [1., 0.]
end
credvec = ones(17)
for (i, k) in enumerate(sort!(collect(keys(regime_eqcond_info))))
    if (get_setting(m, :temporary_altpolicy_length) - 1) < i
        get_setting(m, :regime_eqcond_info)[k].weights = [credvec[i - (get_setting(m, :temporary_altpolicy_length) - 1)],
                                                          1. - credvec[i - (get_setting(m, :temporary_altpolicy_length) - 1)]]
    end
end
out_temp = DSGE.forecast_one_draw(m, :mode, :full, output_vars, modal_params, df,
                                  regime_switching = true, n_regimes = get_setting(m, :n_regimes))

m <= Setting(:uncertain_temporary_altpolicy, false)
m <= Setting(:uncertain_altpolicy, false)

credvec = collect(range(0., stop = 1., length = 17))
for (i, k) in enumerate(sort!(collect(keys(regime_eqcond_info))))
    if i <= length(credvec)
        get_setting(m, :regime_eqcond_info)[k].weights = [credvec[i], 1.0 - credvec[i]]
    end
end

out_perm = DSGE.forecast_one_draw(m, :mode, :full, output_vars, modal_params, df,
                              regime_switching = true, n_regimes = get_setting(m, :n_regimes))

# Check perfectly credible ZLB match permanent equivalent
m <= Setting(:uncertain_temporary_altpolicy, false)
m <= Setting(:uncertain_altpolicy, true)
credvec = collect(range(0., stop = 1., length = 17))
for (i, k) in enumerate(sort!(collect(keys(regime_eqcond_info))))
    if (get_setting(m, :temporary_altpolicy_length) ) < i
        get_setting(m, :regime_eqcond_info)[k].weights = [credvec[i - (get_setting(m, :temporary_altpolicy_length))],
                                                          1. - credvec[i - (get_setting(m, :temporary_altpolicy_length))]]
    end
end

out_credzlb = DSGE.forecast_one_draw(m, :mode, :full, output_vars, modal_params, df,
                                     regime_switching = true, n_regimes = get_setting(m, :n_regimes))

m <= Setting(:uncertain_altpolicy, true)
m <= Setting(:uncertain_temporary_altpolicy, true)
for i in keys(get_setting(m, :regime_eqcond_info))
    get_setting(m, :regime_eqcond_info)[i].weights = [1., 0.]
end
credvec = collect(range(0., stop = 1., length = 17))

for (i, k) in enumerate(sort!(collect(keys(regime_eqcond_info))))
    if (get_setting(m, :temporary_altpolicy_length)) < i
        get_setting(m, :regime_eqcond_info)[k].weights = [credvec[i - (get_setting(m, :temporary_altpolicy_length))],
                                                          1. - credvec[i - (get_setting(m, :temporary_altpolicy_length))]]
    else
        firsted = 1.0
        get_setting(m, :regime_eqcond_info)[k].weights = [firsted, 1.0 - firsted]
    end
end

out_temp_credzlb = DSGE.forecast_one_draw(m, :mode, :full, output_vars, modal_params, df,
                                          regime_switching = true, n_regimes = get_setting(m, :n_regimes))

if !regenerate_reference_forecasts
    @testset "Compare Perfect Credibility to Permanent Alternative Policy" begin
        @test out_temp[:forecastobs] ≈ out_perm[:forecastobs]
        @test out_temp_credzlb[:forecastobs] ≈ out_credzlb[:forecastobs]
    end
end

for i in keys(get_setting(m, :regime_eqcond_info))
    get_setting(m, :regime_eqcond_info)[i].weights = [0., 1.]
end
credvec = collect(range(0., stop = 1., length = 17))
for (i, k) in enumerate(sort!(collect(keys(regime_eqcond_info))))
    if (get_setting(m, :temporary_altpolicy_length) - 1) < i
        get_setting(m, :regime_eqcond_info)[k].weights = [credvec[i - (get_setting(m, :temporary_altpolicy_length) - 1)],
                                                          1. - credvec[i - (get_setting(m, :temporary_altpolicy_length) - 1)]]
    end
end
out1 = DSGE.forecast_one_draw(m, :mode, :full, output_vars, modal_params, df,
                              regime_switching = true, n_regimes = get_setting(m, :n_regimes))

if !regenerate_reference_forecasts
    inds = vcat(1:9, 12:13)
    @testset "Compare TV Credibility to Reference Forecast" begin
        tvtestfcast = h5read(joinpath(dirname(@__FILE__), "..", "reference", "tvcred_reference_forecast.h5"),
                             VERSION >= v"1.5" ? "tvforecastobs_1p5" : "tvforecastobs")
        @test maximum(abs.(out1[:forecastobs][inds, :] - tvtestfcast[inds, :])) < 1e-3
    end
end

if regenerate_reference_forecasts
    h5open(joinpath(dirname(@__FILE__), "..", "reference", "tvcred_reference_meas_fix_forecast.h5"), "w") do file
        write(file, "para", θ10)
        write(file, "forecastobs", outp33[:forecastobs])
        if VERSION >= v"1.5"
            write(file, "tvforecastobs", tvtestfcast_othervers)
            write(file, "tvforecastobs_1p5", out1[:forecastobs])
        else
            write(file, "forecastobs_tv0to1", out1[:forecastobs])
            write(file, "forecastobs_fixed33", outp33[:forecastobs])
            write(file, "forecastobs_fixed0", outp0[:forecastobs])
            write(file, "forecastobs_fixed100", out_temp[:forecastobs])
            write(file, "forecastobs_tv0to1_ZLBcred_1", out_credzlb[:forecastobs])
        end
    end
end

nothing


text:
# sample from weights

struct AliasingScalarSampler
  domain::Vector{Float64}
  weights::StatsBase.ProbabilityWeights
  AliasingScalarSampler() = new()
  AliasingScalarSampler(x::Vector{<:Real}, p_x::Vector{<:Real}; SNRfloor::Float64=0.0) = begin
    pxf = Float64.(p_x)
    pxf .-= quantile(pxf, SNRfloor)
    pxf[pxf .< 0.0] = 0.0
    pxf ./= norm(pxf)
    wim = StatsBase.ProbabilityWeights(pxf)
    new(x, wim)
  end
end

function rand!(ass::AliasingScalarSampler, smpls::Array{Float64})
    StatsBase.alias_sample!(ass.domain, ass.weights, smpls)
    nothing
end

function rand(ass::AliasingScalarSampler, N::Int=1)
  smpls = zeros(N)
  rand!(ass,smpls)
  smpls
end


text:
# Laplacian Eigenmaps
# -------------------
# Laplacian Eigenmaps for Dimensionality Reduction and Data Representation,
# M. Belkin, P. Niyogi, Neural Computation, June 2003; 15 (6):1373-1396

#### LEM type
immutable LEM{T <: Real} <: SpectralResult
    k::Int
    λ::AbstractVector{T}
    t::T
    proj::Projection{T}
    component::Vector{Int}

    LEM{T}(k::Int, λ::AbstractVector{T}, t::T, proj::Projection{T}) = new(k, λ, t, proj)
    LEM{T}(k::Int, λ::AbstractVector{T}, t::Float64, proj::Projection{T}, cc::Vector{Int}) = new(k, λ, t, proj, cc)
end

## properties
outdim(M::LEM) = size(M.proj, 1)
projection(M::LEM) = M.proj

eigvals(M::LEM) = M.λ
neighbors(M::LEM) = M.k
ccomponent(M::LEM) = M.component

## show & dump
function show(io::IO, M::LEM)
    print(io, "Laplacian Eigenmaps(outdim = $(outdim(M)), neighbors = $(neighbors(M)), t = $(M.t))")
end

function dump(io::IO, M::Isomap)
    show(io, M)
    println(io, "eigenvalues: ")
    Base.showarray(io, eigvals(M)', header=false, repr=false)
    println(io)
    println(io, "projection:")
    Base.showarray(io, projection(M), header=false, repr=false)
end

## interface functions
function transform{T<:Real}(::Type{LEM}, X::DenseMatrix{T}; d::Int=2, k::Int=12, t::T=1.0)
    n = size(X, 2)

    # Construct NN graph
    D, E = find_nn(X, k)

    W = zeros(T,n,n)
    for i = 1 : n
        jj = E[:, i]
        W[i,jj] = D[:, i]
    end
    W .^= 2
    W ./= maximum(W)

    # Select largest connected component
    CC = components(E)
    C = length(CC) == 1 ? CC[1] : CC[indmax(map(size, CC))]

    # Compute weights
    W[W .> eps(T)] = exp(-W[W .> eps(T)] ./ t)
    D = diagm(sum(W,2)[:])
    L = D - W

    λ, V = decompose(L, D, d)
    return LEM{T}(k, λ, t, V', C)
end


text:
using jInv.Mesh
using DivSigGrad
using Multigrid

println("************************************************* Example 2D ******************************************************");

domain = [0.0, 1.0, 0.0, 1.0];
n = [50,50];
Mr = getRegularMesh(domain,n)
m  = exp.(randn(prod(n)));
Ar = getDivSigGradMatrix(vec(m),Mr);
Ar = Ar + 1e-8*norm(Ar,1)*speye(size(Ar,2));

levels      = 3;
numCores 	= 8;
maxIter     = 3;
relativeTol = 1e-2;
relaxType   = "SPAI";
relaxParam  = 1.0;
relaxPre 	= 2;
relaxPost   = 2;
cycleType   ='V';
coarseSolveType = "NoMUMPS";

MG = getMGparam(levels,numCores,maxIter,relativeTol,relaxType,relaxParam,relaxPre,relaxPost,cycleType,coarseSolveType);

N = size(Ar,2);

b = Ar*rand(N,3);
x = zeros(N,3);
SA_AMGsetup(Ar,MG,Float64,true,size(b,2),true);

println("****************************** Stand-alone AMG: ******************************")
solveMG(MG,b,x,true);

x[:] = 0.0;
println("****************************** CG preconditioned with AMG: ******************************")
solveCG_MG(Ar,MG,b,x,true)

println("****************************** CG preconditioned with BiCGSTAB: ******************************")
solveBiCGSTAB_MG(Ar,MG,b,x,true);

println("****************************** transposeHierarchy *********************************");
transposeHierarchy(MG,true);

println("****************************** GMRES preconditioned with AMG: (only one rhs...) ******************************")
x[:] = 0.0;
b1 = vec(b[:,1]);
x1 = zeros(N);
MG = getMGparam(levels,numCores,maxIter,relativeTol,relaxType,relaxParam,relaxPre,relaxPost,cycleType,coarseSolveType);
SA_AMGsetup(Ar,MG,Float64,true,1,true);
solveGMRES_MG(Ar,MG,b1,x1,true,2)

x1 = zeros(N);
MG.maxOuterIter = 10;
solveBiCGSTAB_MG(Ar,MG,b1,x1,true)

MG = getMGparam(levels,numCores,maxIter,relativeTol,"Jac-GMRES",relaxParam,relaxPre,relaxPost,'K',coarseSolveType);
SA_AMGsetup(Ar,MG,Float64,true,size(b,2),true);
x[:] = 0.0;
solveGMRES_MG(Ar,MG,b,x,true,2)

Ar = 0;
b = 0;
x = 0;
Mr = 0;


text:
module jobs

end

max_stars_repo_name: tmaone/xsh

max_issues_repo_name: turing-complete/xsh

max_forks_repo_name: turing-complete/xsh


text:
__precompile__()

module GeneralizedChartParsing

# from standard library
using Distributed
using SpecialFunctions: digamma
using SparseArrays:     SparseMatrixCSC, sparse, spzeros, findnz
using LinearAlgebra:    normalize

# other packages
using LogProbs
using FastClosures
using DataFrames
using Distributions: Dirichlet

# imports
import Base: zero, one, +, *
import Base: IteratorSize, eltype, iterate
import Base: parse, rand, merge, convert, show

# exports
export NamedFunction
export Grammar, grammar_from_string, train_grammar
export add_score, set_scores
export add_random_prob_score, add_expected_count_score, add_forest_score
export count_score, enum_forest_score, forest_score
export tree_struct, tree_structs, expected_counts_dataframe

include("utils.jl")
include("LazyDicts.jl")
include("newtrees.jl")
include("Grammars.jl")
include("scores.jl")
include("grammar_from_string.jl")
include("variational_bayes.jl")

using .Trees
using .LazyDicts

end # module


text:
@testset "KeyedDataset" begin
    @testset "Construction" begin
        @testset "Invalid" begin
            # Test that you can construct an invalid dataset if you really want to
            @test_throws KeyAlignmentError KeyedDataset(
                OrderedSet(Pattern[(:__, :time)]),
                LittleDict(
                    (:val1,) => KeyedArray(rand(4); time=1:4),
                    (:val2,) => KeyedArray(rand(4); time=2:5),
                ),
            )

            ds = KeyedDataset(
                OrderedSet(Pattern[(:__, :time)]),
                LittleDict(
                    (:val1,) => KeyedArray(rand(4); time=1:4),
                    (:val2,) => KeyedArray(rand(4); time=2:5),
                ),
                false
            )
            @test_throws KeyAlignmentError validate(ds)
        end

        @testset "Empty" begin
            expected = KeyedDataset(OrderedSet{Pattern}(), LittleDict{Tuple, KeyedArray}())
            @test KeyedDataset() == expected

            patterns = Pattern[(:train, :_, :target)]
            expected = KeyedDataset(OrderedSet(patterns), LittleDict{Tuple, KeyedArray}())
            @test KeyedDataset(; constraints=patterns) == expected
        end

        @testset "KeyedArrays" begin
            ds = KeyedDataset(
                :val1 => KeyedArray(
                    rand(4, 3, 2);
                    time=DateTime(2021, 1, 1, 11):Hour(1):DateTime(2021, 1, 1, 14),
                    loc=1:3,
                    obj=[:a, :b],
                ),
                :val2 => KeyedArray(
                    rand(4, 3, 2) .+ 1.0;
                    time=DateTime(2021, 1, 1, 11):Hour(1):DateTime(2021, 1, 1, 14),
                    loc=1:3,
                    obj=[:a, :b],
                ),
            )

            # Test that we successfully extracted the dims
            @test issetequal([:time, :loc, :obj], dimnames(ds))

            # Test generated constraints
            @test constraintmap(ds) == LittleDict(
                Pattern((:__, :time)) => Set([(:val1, :time), (:val2, :time)]),
                Pattern((:__, :loc)) => Set([(:val1, :loc), (:val2, :loc)]),
                Pattern((:__, :obj)) => Set([(:val1, :obj), (:val2, :obj)]),
            )

            # Test that we have a data dict entry for each value column
            @test issetequal([(:val1,), (:val2,)], keys(ds.data))

            @testset "Variable Keys" begin
                # Test construction with different key types and lengths
                ds = KeyedDataset(
                    "val1" => KeyedArray(
                        rand(4, 3, 2);
                        time=DateTime(2021, 1, 1, 11):Hour(1):DateTime(2021, 1, 1, 14),
                        loc=1:3,
                        obj=[:a, :b],
                    ),
                    (:group1, 2) => KeyedArray(
                        rand(4, 3, 2) .+ 1.0;
                        time=DateTime(2021, 1, 1, 11):Hour(1):DateTime(2021, 1, 1, 14),
                        loc=1:3,
                        obj=[:a, :b],
                    ),
                )

                # Test that we successfully extracted the dims
                @test issetequal([:time, :loc, :obj], dimnames(ds))
                # Test that we have a data dict entry for each value column
                @test issetequal([("val1",), (:group1, 2)], keys(ds.data))
            end
        end

        @testset "Flatten" begin
            # Technically this test is just showing how you can use the `flatten` function
            # to construct a KeyedDatasets from an existing nested structure of KeyedArrays.
            @testset "NamedTuples" begin
                # In this case, the resulting keys from flatten need to be symbols with
                # an `:_` delimiter
                data = (
                    group1 = (
                        a = KeyedArray(
                            rand(4, 3, 2);
                            time=DateTime(2021, 1, 1, 11):Hour(1):DateTime(2021, 1, 1, 14),
                            loc=1:3,
                            obj=[:a, :b],
                        ),
                        b = KeyedArray(
                            rand(4, 2);
                            time=DateTime(2021, 1, 1, 11):Hour(1):DateTime(2021, 1, 1, 14),
                            label=["x", "y"],
                        )
                    ),
                    group2 = (
                        a = KeyedArray(
                            rand(4, 3, 2) .+ 1.0;
                            time=DateTime(2021, 1, 1, 11):Hour(1):DateTime(2021, 1, 1, 14),
                            loc=1:3,
                            obj=[:a, :b],
                        ),
                        b = KeyedArray(
                            rand(4, 2) .+ 1.0;
                            time=DateTime(2021, 1, 1, 11):Hour(1):DateTime(2021, 1, 1, 14),
                            label=["x", "y"],
                        )
                    )
                )
                ds = KeyedDataset(flatten(data)...)

                # Test that we successfully extracted the dims
                @test issetequal([:time, :loc, :obj, :label], dimnames(ds))

                # Test generated constraints
                @test constraintmap(ds) == LittleDict(
                    Pattern((:__, :time)) => Set([
                        (:group1, :a, :time),
                        (:group1, :b, :time),
                        (:group2, :a, :time),
                        (:group2, :b, :time),
                    ]),
                    Pattern((:__, :loc)) => Set([(:group1, :a, :loc), (:group2, :a, :loc)]),
                    Pattern((:__, :obj)) => Set([(:group1, :a, :obj), (:group2, :a, :obj)]),
                    Pattern((:__, :label)) => Set([
                        (:group1, :b, :label),
                        (:group2, :b, :label),
                    ]),
                )

                # Test that we successfully extracted the flattened kwargs.
                @test issetequal(
                    [(:group1, :a), (:group1, :b), (:group2, :a), (:group2, :b)], keys(ds.data)
                )

                # Test an example where we pass a delimiter to flatten to a named tuple,
                # that we can use as kwargs...
                ds = KeyedDataset(; flatten(data, :_)...)

                @test constraintmap(ds) == LittleDict(
                    Pattern((:__, :time)) => Set([
                        (:group1_a, :time),
                        (:group1_b, :time),
                        (:group2_a, :time),
                        (:group2_b, :time),
                    ]),
                    Pattern((:__, :loc)) => Set([(:group1_a, :loc), (:group2_a, :loc)]),
                    Pattern((:__, :obj)) => Set([(:group1_a, :obj), (:group2_a, :obj)]),
                    Pattern((:__, :label)) => Set([(:group1_b, :label), (:group2_b, :label)]),
                )

                @test issetequal(
                    [(:group1_a,), (:group1_b,), (:group2_a,), (:group2_b,)], keys(ds.data)
                )
            end
            @testset "Pairs" begin
                # In this case, we have the option to having tuple keys which allows for
                # multi-dimensional style indexing.
                data = [
                    :group1 => [
                        :a => KeyedArray(
                            rand(4, 3, 2);
                            time=DateTime(2021, 1, 1, 11):Hour(1):DateTime(2021, 1, 1, 14),
                            loc=1:3,
                            obj=[:a, :b],
                        ),
                        :b => KeyedArray(
                            rand(4, 2);
                            time=DateTime(2021, 1, 1, 11):Hour(1):DateTime(2021, 1, 1, 14),
                            label=["x", "y"],
                        )
                    ],
                    :group2 => [
                        :a => KeyedArray(
                            rand(4, 3, 2) .+ 1.0;
                            time=DateTime(2021, 1, 1, 11):Hour(1):DateTime(2021, 1, 1, 14),
                            loc=1:3,
                            obj=[:a, :b],
                        ),
                        :b => KeyedArray(
                            rand(4, 2) .+ 1.0;
                            time=DateTime(2021, 1, 1, 11):Hour(1):DateTime(2021, 1, 1, 14),
                            label=["x", "y"],
                        )
                    ]
                ]
                ds = KeyedDataset(flatten(data)...)

                # Test that we successfully extracted the dims
                @test issetequal([:time, :loc, :obj, :label], dimnames(ds))

                # Test that we successfully extracted the flattened pairs as tuples.
                @test issetequal(
                    [(:group1, :a), (:group1, :b), (:group2, :a), (:group2, :b)], keys(ds.data)
                )

                # We can still choose to flatten to a single symbol if we want.
                ds = KeyedDataset(flatten(data, :_)...)

                @test issetequal(
                    [(:group1_a,), (:group1_b,), (:group2_a,), (:group2_b,)], keys(ds.data)
                )

                @test constraintmap(ds) == LittleDict(
                    Pattern((:__, :time)) => Set([
                        (:group1_a, :time),
                        (:group1_b, :time),
                        (:group2_a, :time),
                        (:group2_b, :time),
                    ]),
                    Pattern((:__, :loc)) => Set([(:group1_a, :loc), (:group2_a, :loc)]),
                    Pattern((:__, :obj)) => Set([(:group1_a, :obj), (:group2_a, :obj)]),
                    Pattern((:__, :label)) => Set([(:group1_b, :label), (:group2_b, :label)]),
                )
            end

            @testset "Mixed key types" begin
                # In this case, we have the option to having tuple keys which allows for
                # multi-dimensional style indexing.
                data = [
                    :group1 => [
                        "a" => KeyedArray(
                            rand(4, 3, 2);
                            time=DateTime(2021, 1, 1, 11):Hour(1):DateTime(2021, 1, 1, 14),
                            loc=1:3,
                            obj=[:a, :b],
                        ),
                        "b" => KeyedArray(
                            rand(4, 2);
                            time=DateTime(2021, 1, 1, 11):Hour(1):DateTime(2021, 1, 1, 14),
                            label=["x", "y"],
                        )
                    ],
                    :group2 => [
                        "a" => KeyedArray(
                            rand(4, 3, 2) .+ 1.0;
                            time=DateTime(2021, 1, 1, 11):Hour(1):DateTime(2021, 1, 1, 14),
                            loc=1:3,
                            obj=[:a, :b],
                        ),
                        "b" => KeyedArray(
                            rand(4, 2) .+ 1.0;
                            time=DateTime(2021, 1, 1, 11):Hour(1):DateTime(2021, 1, 1, 14),
                            label=["x", "y"],
                        )
                    ]
                ]
                ds = KeyedDataset(flatten(data)...)

                # Test that we successfully extracted the dims
                @test issetequal([:time, :loc, :obj, :label], dimnames(ds))

                # Test that we successfully extracted the flattened pairs as tuples.
                @test issetequal(
                    [(:group1, "a"), (:group1, "b"), (:group2, "a"), (:group2, "b")], keys(ds.data)
                )
            end
        end
    end
    @testset "show" begin
        ds = KeyedDataset(
            :val1 => KeyedArray(
                rand(4, 3, 2);
                time=DateTime(2021, 1, 1, 11):Hour(1):DateTime(2021, 1, 1, 14),
                loc=1:3,
                obj=[:a, :b],
            ),
            :val2 => KeyedArray(
                rand(4, 3, 2) .+ 1.0;
                time=DateTime(2021, 1, 1, 11):Hour(1):DateTime(2021, 1, 1, 14),
                loc=1:3,
                obj=[:a, :b],
            ),
        )

        # This is likely to change in the future, so we keep this test simple
        s = sprint(show, ds)
        @test startswith(s, "KeyedDataset")
        @test occursin("2 components", s)
        @test occursin("3 constraints", s)

        # Having unused constraints is fine since we might add them later.
        s = sprint(show, KeyedDataset(pairs(ds.data)...; constraints=Pattern[(:__, :foo)]))
        @test occursin("2 components", s)
        @test occursin("1 constraints", s)

        # While possible to construct it's rather pointless to have a dimension path
        # that can map to multiple constraints.
        bad_ds = KeyedDataset(
            pairs(ds.data)...;
            constraints=Pattern[
                (:__, :time),
                (:val1, :time),
            ],
        )
        @test_throws ArgumentError sprint(show, bad_ds)
    end

    @testset "dimpaths" begin
        ds = KeyedDataset(
            :val1 => KeyedArray(
                rand(4, 3, 2);
                time=DateTime(2021, 1, 1, 11):Hour(1):DateTime(2021, 1, 1, 14),
                loc=1:3,
                obj=[:a, :b],
            ),
            :val2 => KeyedArray(
                rand(4, 3, 2) .+ 1.0;
                time=DateTime(2021, 1, 1, 11):Hour(1):DateTime(2021, 1, 1, 14),
                loc=1:3,
                obj=[:a, :b],
            ),
        )

        @test issetequal(
            dimpaths(ds),
            [
                (:val1, :time),
                (:val1, :loc),
                (:val1, :obj),
                (:val2, :time),
                (:val2, :loc),
                (:val2, :obj),
            ],
        )

        @test issetequal(dimpaths(ds, Pattern(:__, :time)), [(:val1, :time), (:val2, :time)])
    end

    @testset "dimnames" begin
        ds = KeyedDataset(
            :val1 => KeyedArray(
                rand(4, 3, 2);
                time=DateTime(2021, 1, 1, 11):Hour(1):DateTime(2021, 1, 1, 14),
                loc=1:3,
                obj=[:a, :b],
            ),
            :val2 => KeyedArray(
                rand(4, 3, 2) .+ 1.0;
                time=DateTime(2021, 1, 1, 11):Hour(1):DateTime(2021, 1, 1, 14),
                id=1:3,
                obj=[:a, :b],
            ),
        )
        @test issetequal(dimnames(ds), [:time, :loc, :obj, :id])
    end

    @testset "axiskeys" begin
        ds = KeyedDataset(
            :val1 => KeyedArray(
                rand(4, 3, 2);
                time=DateTime(2021, 1, 1, 11):Hour(1):DateTime(2021, 1, 1, 14),
                loc=1:3,
                obj=[:a, :b],
            ),
            :val2 => KeyedArray(
                rand(4, 3, 2) .+ 1.0;
                time=DateTime(2021, 1, 1, 11):Hour(1):DateTime(2021, 1, 1, 14),
                id=1:3,
                obj=[:a, :b],
            ),
        )

        @test issetequal(
            axiskeys(ds),
            [
                DateTime(2021, 1, 1, 11):Hour(1):DateTime(2021, 1, 1, 14),
                1:3,
                [:a, :b],
            ],
        )

        @test axiskeys(ds, (:val1, :time)) == DateTime(2021, 1, 1, 11):Hour(1):DateTime(2021, 1, 1, 14)
        @test unique(axiskeys(ds, :time)) == [DateTime(2021, 1, 1, 11):Hour(1):DateTime(2021, 1, 1, 14)]
    end

    @testset "validate" begin
        ds = KeyedDataset(
            :val1 => KeyedArray(
                rand(4, 3, 2);
                time=DateTime(2021, 1, 1, 11):Hour(1):DateTime(2021, 1, 1, 14),
                loc=1:3,
                obj=[:a, :b],
            ),
            :val2 => KeyedArray(
                rand(4, 3, 2) .+ 1.0;
                time=DateTime(2021, 1, 1, 11):Hour(1):DateTime(2021, 1, 1, 14),
                id=1:3,
                obj=[:a, :b],
            ),
        )

        @test validate(ds)
        @test validate(ds, Pattern(:__, :time))

        # Intentionally break the internal keys
        ds.data[(:val1,)] = KeyedArray(
            rand(4, 3, 2);
            time=DateTime(2021, 1, 2, 11):Hour(1):DateTime(2021, 1, 2, 14),
            loc=1:3,
            obj=[:a, :b],
        )

        @test_throws KeyAlignmentError validate(ds)
        @test_throws KeyAlignmentError validate(ds, Pattern(:__, :time))
        @test validate(ds, Pattern(:__, :obj))
    end
end


text:
using HalideCall
const libaccum = "./libaccum.so"
A = rand(Float32,2000,3000);
out = similar(A);
Abuf = HalideBuffer(A)
outbuf = HalideBuffer(out, (1,1))
fill!(out, 0);
ccall((:accum_vectorized,libaccum), Void, (Ptr{Void}, Ptr{Void}), [Abuf], [outbuf])
@time (for i = 1:10; ccall((:accum_vectorized,libaccum), Void, (Ptr{Void}, Ptr{Void}), [Abuf], [outbuf]); end)
@time (for i = 1:10; ccall((:accum_tiled,libaccum), Void, (Ptr{Void}, Ptr{Void}), [Abuf], [outbuf]); end)

using KernelTools

function accum_vectorized!(out, A)
    for y = 1:size(A, 2), x = 1:size(A, 1)
        @inbounds out[x, y] += A[x, y]
    end
    out
end

function accum_tiled_byhand!(out, A)
    for y = 1:size(A, 2), x = 1:32:size(A, 1)
        for xi = x:min(x+31, size(A, 1))
            @inbounds out[xi, y] += A[xi, y]
        end
    end
    out
end

function accum_tiled!(out, A)
    KernelTools.@tile (y,32,x,256) begin
        for y = 1:size(A, 2), x = 1:size(A, 1)
            @inbounds out[x, y] += A[x, y]
        end
    end
    out
end

function accum_tiledx!(out, A)
    for y = 1:size(A, 2)
        KernelTools.@tile (x,32) begin
            for x = 1:size(A, 1)
                @inbounds out[x, y] += A[x, y]
            end
        end
    end
    out
end

accum_vectorized!(out, A)
gc()
@time (for i = 1:10; accum_vectorized!(out, A); end)
accum_tiled_byhand!(out, A)
@time (for i = 1:10; accum_tiled_byhand!(out, A); end)
accum_tiled!(out, A)
@time (for i = 1:10; accum_tiled!(out, A); end)
accum_tiledx!(out, A)
@time (for i = 1:10; accum_tiledx!(out, A); end)


text:
module MobileWithLogTiming

# For profiling

using ModiaLang

# ModiaLang models
include("$(ModiaLang.path)/models/Blocks.jl")
include("$(ModiaLang.path)/models/Electric.jl")
include("$(ModiaLang.path)/models/Rotational.jl")

import Modia3D
using  Modia3D.ModiaInterface

const depthMax = 6
const enableVisualization = false

const rodLength       = 0.1
const rodDiameter     = 0.05*rodLength
const rodMaterial     = "DryWood"
const rodColor        = VisualMaterial(color="Brown", transparency=0.5)
const rodTranslation1 = [0.0,  rodLength/2, 0.0]
const rodTranslation2 = [0.0, -rodLength/2, 0.0]

const sphereDiameter  = 0.3*rodLength
const sphereMaterial  = "DryWood"
const sphereColor     = VisualMaterial(color="Blue")

const barWidth       = 0.3*rodLength
const barThickness   = 0.3*rodLength
const barMaterial    = "DryWood"
const barColor       = VisualMaterial(color="DeepSkyBlue3", transparency=0.5)
const leaveBarLength = 0.8*rodLength

const revDiameter = rodDiameter
const revLength   = 1.2*barThickness
const revColor    = VisualMaterial(color="Red")

const damping = 0.006

Rod = Model(
    frame0 = Object3D(feature = Solid(shape = Cylinder(axis=2, length=rodLength, diameter=rodDiameter),
                                    solidMaterial  = rodMaterial,
                                    visualMaterial = rodColor)),
    frame1 = Object3D(parent = :frame0, translation = rodTranslation1),
    frame2 = Object3D(parent = :frame0, translation = rodTranslation2)
)
Bar = Model(
    L = missing,
    frame0 = Object3D(feature = Solid(shape = Beam(axis = 1, length = :L, width=barWidth, thickness=barThickness),
                                    solidMaterial  = barMaterial,
                                    visualMaterial = barColor)),
    frame1 = Object3D(parent = :frame0, translation = :[-L/2, 0.0, 0.0]),
    frame2 = Object3D(parent = :frame0, translation = :[ L/2, 0.0, 0.0])
)
RevoluteWithDamping(;obj1, obj2, axis=3, phi_start=0.0) = Model(
    rev     = RevoluteWithFlange(obj1=obj1, obj2=obj2, axis=axis, phi=Var(init=phi_start)),
    cyl     = Object3D(parent  = obj1,
                    feature = Visual(shape = Cylinder(axis=axis, diameter=revDiameter, length=revLength),
                                        visualMaterial = revColor)),
    damper  = Damper | Map(d=damping),
    fixed   = Fixed,
    connect = :[(damper.flange_b, rev.flange),
                (damper.flange_a, fixed.flange)]
)
function barLength(depth)
    if depth<=2
        return leaveBarLength
    end
    return 2*barLength(depth-1)
end
function createMobile(depth)
    if depth == 1
        Model(
            rod    = Rod,
            sphere = Object3D(parent = :(rod.frame0), translation = rodTranslation2,
                              feature = Solid(shape = Sphere(diameter=sphereDiameter),
                                              solidMaterial  = sphereMaterial,
                                              visualMaterial = sphereColor))
        )
    else
        Model(
            rod  = Rod,
            bar  = Bar | Map(L=barLength(depth)),
            sub1 = createMobile(depth-1),
            sub2 = createMobile(depth-1),
            rev0 = RevoluteWithDamping(obj1=:(rod.frame2), obj2=:(bar.frame0)),
            rev1 = RevoluteWithDamping(obj1=:(bar.frame1), obj2=:(sub1.rod.frame1)),
            rev2 = RevoluteWithDamping(obj1=:(bar.frame2), obj2=:(sub2.rod.frame1))
        )
    end
end
Mobile = Model(
    world      = Object3D(feature=Scene(gravityField=UniformGravityField(g=9.81, n=[0, -1, 0]),
                                        enableContactDetection=false,
                                        enableVisualization=enableVisualization)),
    worldFrame = Object3D(parent = :world, feature = Visual(shape = CoordinateSystem(length=rodLength))),
    top        = createMobile(depthMax),
    rev0       = RevoluteWithDamping(obj1 = :world, obj2 = :(top.rod.frame1), phi_start=0.2)
)

println("... @instantiateModel:")
@time mobile = @instantiateModel(buildModia3D(Mobile), unitless=true, log=false, logDetails=false, logModel=false, logStateSelection=false,
                                 logCode=false, logExecution=true, logTiming=false, evaluateParameters=true)

const stopTime = 5.0
const tolerance = 1e-4
const requiredFinalStates = missing

println("... first simulation:")
@time simulate!(mobile, stopTime=stopTime, tolerance=tolerance, logTiming=false, requiredFinalStates=requiredFinalStates)

println("... second simulation:")
@time simulate!(mobile, stopTime=stopTime, tolerance=tolerance, logTiming=true, log=true, requiredFinalStates=requiredFinalStates)

#@usingModiaPlot
#plot(mobile, "rev0.rev.phi")

end


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

using Markdown
using InteractiveUtils

# ╔═╡ 15b8afb9-ace9-4b5b-894e-dce9dda5bf5f
begin
	using Pkg
	using Random
end

# ╔═╡ 45d72214-919c-11eb-2e69-23c0786baacb
md"""
### Functions as variable

Src: [Julia 1.0 Cookbook -  Bogumił Kamiński , Przemysław Szufel - 2018](https://www.packtpub.com/product/julia-1-0-programming-cookbook/9781788998369)
"""

# ╔═╡ 9608fdc8-807e-4be2-91f8-950c0ceec2dd
Random.seed!(0)

# ╔═╡ c15cd8d2-3178-4490-bcdb-7ec54ec41318
md"""
In this example, we will consider a simplified agent-based simulation model, where a group of agents moves randomly over a two-dimensional surface. Additionally, we assume that, at each step of the simulation, the closer the agent is to the starting point, the lower the probability of movement.

the probability of a change in an agent's location is equal to $\frac{1}{k}$, where $k$ is the rank of an agent (the agent that managed to travel the longest distance from its origin has the highest rank).
"""

# ╔═╡ 1911333e-be43-4294-a374-ce15a7538a42
mutable struct Agent
	id::Int
	x::Float64
	y::Float64
	times_moved::Int
end

# ╔═╡ e33a47c8-ce78-4aad-8deb-a97c9cb47b25
function move!(agent::Agent)
	angle = rand() * 2π
	agent.x += cos(angle)
	agent.y += sin(angle)
	agent.times_moved += 1
end

# ╔═╡ 03409d5c-a6d2-43c4-9367-c52b87a9da01
function step!(pop::Vector{Agent})
	sort!(pop, by=a -> (a.x * a.x + a.y * a.y), rev=true)
	foreach(ix -> (rand() < (1/ix) && move!(pop[ix])), 1:length(pop))
end

# ╔═╡ afd70dd1-a99a-4003-a57f-e037855c7ad1
md"""
Now let us define a population of 30 agents:
"""

# ╔═╡ c2ce27b3-64f3-4a9f-a72d-8c56a8492fcc
pop = Agent.(1:30, 0., 0., 0)

# ╔═╡ ca503785-20e7-4ed0-bb4a-6330d710c27b
md"""
Now we define what happens at each simulation time step:
"""

# ╔═╡ d9f6a75e-c588-42e7-a815-318c5e33f706
md"""
Let us run the simulation:
"""

# ╔═╡ 476e9c3b-805b-4e89-a537-238c95ed546c
foreach(_a -> step!(pop), 1:1000)

# ╔═╡ d86c15dc-c0a1-4834-9328-055ef838f756
md"""
Finally, let us filter out the agents that are at least 25 units away from the origin:
"""

# ╔═╡ 30907ce3-a6c4-47f7-9aea-c40fb4368599
filter(a -> √(a.x * a.x + a.y * a.y) ≥ 25., pop)

# ╔═╡ 20ffa0af-d9ac-46a1-a6bd-597f2f8f1bed
md"""
In Julia, a function is a first-class object.
We start by defining `struct` that we will operate on. We also define a `move!` function that mutates an agent's state.

The `map` function takes two parameters: a function and a collection. The function given as the parameter is executed over each element of the collection, and the result is returned as Array.

In the `step!` function, when sorting with `sort!`, we provide the `by=` parameter. The value of this parameter is an anonymous function that calculates a weight for each Array's element. The `foreach` function allows us to execute the `move!` function for each agent. Once again, we provide an anonymous function as a parameter to the `move!` function.

Finally, we again use the `foreach` function to actually run the simulation. The `filter` function can be used to filter out those Array elements having a desired set of values.
"""

# ╔═╡ Cell order:
# ╠═45d72214-919c-11eb-2e69-23c0786baacb
# ╠═15b8afb9-ace9-4b5b-894e-dce9dda5bf5f
# ╠═9608fdc8-807e-4be2-91f8-950c0ceec2dd
# ╟─c15cd8d2-3178-4490-bcdb-7ec54ec41318
# ╠═1911333e-be43-4294-a374-ce15a7538a42
# ╠═e33a47c8-ce78-4aad-8deb-a97c9cb47b25
# ╠═03409d5c-a6d2-43c4-9367-c52b87a9da01
# ╟─afd70dd1-a99a-4003-a57f-e037855c7ad1
# ╠═c2ce27b3-64f3-4a9f-a72d-8c56a8492fcc
# ╟─ca503785-20e7-4ed0-bb4a-6330d710c27b
# ╟─d9f6a75e-c588-42e7-a815-318c5e33f706
# ╠═476e9c3b-805b-4e89-a537-238c95ed546c
# ╟─d86c15dc-c0a1-4834-9328-055ef838f756
# ╠═30907ce3-a6c4-47f7-9aea-c40fb4368599
# ╟─20ffa0af-d9ac-46a1-a6bd-597f2f8f1bed

max_stars_repo_name: pascal-p/julia-notebooks

max_issues_repo_name: pascal-p/julia-notebooks

max_forks_repo_name: pascal-p/julia-notebooks


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

name = "gmmlib"
version = v"20.2.5"

# Collection of sources required to build this package
sources = [
    GitSource("https://github.com/intel/gmmlib.git",
              "c1b7e077b5b0b6b97057562aaec574f0da97b37b"),
]

# Bash recipe for building across all platforms
script = raw"""
cd gmmlib
install_license LICENSE.md

CMAKE_FLAGS=()

if [[ ${target} == *musl* ]]; then
    # https://www.openwall.com/lists/musl/2018/09/11/2
    sed -i '/-fstack-protector/d' Source/GmmLib/Linux.cmake
fi

# Release build for best performance
# https://github.com/intel/gmmlib/issues/70
CMAKE_FLAGS+=(-DBUILD_TYPE=Release)

# Install things into $prefix
CMAKE_FLAGS+=(-DCMAKE_INSTALL_PREFIX=${prefix})

# Explicitly use our cmake toolchain file and tell CMake we're cross-compiling
CMAKE_FLAGS+=(-DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN})
CMAKE_FLAGS+=(-DCMAKE_CROSSCOMPILING:BOOL=ON)

# Don't run tests
CMAKE_FLAGS+=(-DRUN_TEST_SUITE:Bool=OFF)

cmake -B build -S . -GNinja ${CMAKE_FLAGS[@]}
ninja -C build -j ${nproc} install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Platform("i686", "linux"; libc="glibc"),
    Platform("x86_64", "linux"; libc="glibc"),
    Platform("i686", "linux"; libc="musl"),
    Platform("x86_64", "linux"; libc="musl"),
]

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

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

build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
               preferred_gcc_version=v"5", lock_microarchitecture=false)

max_stars_repo_name: mortenpi/Yggdrasil

max_issues_repo_name: mortenpi/Yggdrasil

max_forks_repo_name: mortenpi/Yggdrasil


text:
"""
    FixedRateDescent

Descent method which minimizes the objective function in the direction
of antigradient at each step.
"""
mutable struct FixedRateDescent{T<:AbstractFloat,V<:AbstractVector{T}} <: OptBuffer
    x::V
    xpre::V
    g::V
    α::T
    y::T
    ypre::T
end

function FixedRateDescent(x::AbstractVector, α::Real)
    F = float(eltype(x))
    return FixedRateDescent(
        similar(x, F),
        similar(x, F),
        similar(x, F),
        convert(F, α),
        F(NaN),
        F(NaN),
    )
end

FixedRateDescent(x::AbstractVector) = FixedRateDescent(x, 0.01)

function init!(optfn!, M::FixedRateDescent{T}, x0; kw...) where {T}
    optfn!(x0, zero(T), x0)
    return
end

@inline reset!(::FixedRateDescent) = nothing

function reset!(M::FixedRateDescent, x0, α = M.α)
    if length(M.x) != length(x0)
        foreach((M.x, M.g)) do v
            resize!(v, length(x0))
        end
    end
    M.α = α
    return
end

@inline function callfn!(fdf, M::FixedRateDescent, x, α, d)
    __update_arg!(M, x, α, d)
    y, g = fdf(M.x, M.g)
    __update_grad!(M, g)
    return y, g
end

function step!(optfn!, M::FixedRateDescent; constrain_step = infstep)
    d = rmul!(M.g, -1)
    maxstep = constrain_step(M.x, d)
    s = M.α <= maxstep ? M.α : maxstep / 2
    optfn!(M.x, s, d)
    return s
end

@inline function __update_arg!(M::FixedRateDescent, x, α, d)
    map!(M.x, d, x) do a, b
        muladd(α, a, b)
    end
    return
end

@inline function __update_arg!(M::FixedRateDescent, x)
    if x !== M.x
        copy!(M.x, x)
    end
    return
end

@inline function __update_grad!(M::FixedRateDescent, g)
    if M.g !== g
        copy!(M.g, g)
    end
    return
end

"""
    MomentumDescent

Descent method which minimizes the objective function in the direction
of antigradient at each step.
"""
mutable struct MomentumDescent{T<:AbstractFloat,V<:AbstractVector{T}} <: OptBuffer
    x::V
    xpre::V
    g::V
    v::V
    learn_rate::T
    decay_rate::T
    y::T
    ypre::T
end

function MomentumDescent(x::AbstractVector; learn_rate::Real=0.01, decay_rate::Real=0.9)
    F = float(eltype(x))
    return MomentumDescent(
        similar(x, F),
        similar(x, F),
        similar(x, F),
        similar(x, F),
        convert(F, learn_rate),
        convert(F, decay_rate),
        F(NaN),
        F(NaN),
    )
end

function init!(optfn!, M::MomentumDescent{T}, x0; kw...) where {T}
    optfn!(x0, zero(T), x0)
    fill!(M.v, zero(T))
    return M
end

@inline function reset!(M::MomentumDescent)
    fill!(M.v, false)
    return M
end

function reset!(
    M::MomentumDescent,
    x0,
    learn_rate = M.learn_rate,
    decay_rate = M.decay_rate,
)
    n = length(x0)
    if length(M.x) != n
        for v in (M.x, M.g, M.v)
            resize!(v, n)
        end
    end
    M.v .= false
    M.x .= x0
    M.learn_rate = learn_rate
    M.decay_rate = decay_rate
    return M
end

@inline function callfn!(fdf, M::MomentumDescent, x, α, d)
    __update_arg!(M, x, α, d)
    y, g = fdf(M.x, M.g)
    __update_grad!(M, g)
    return y, g
end

function step!(optfn!, M::MomentumDescent; constrain_step = infstep)
    map!(M.v, M.v, M.g) do v, g
        M.decay_rate * v - M.learn_rate * g
    end
    maxstep = constrain_step(M.x, M.v)
    if maxstep <= 1
        M.v *= maxstep / 2
    end
    optfn!(M.x, one(maxstep), M.v)
    return M.learn_rate
end

@inline function __update_arg!(M::MomentumDescent, x, α, d)
    map!(M.x, d, x) do a, b
        muladd(α, a, b)
    end
    return
end

@inline function __update_arg!(M::MomentumDescent, x)
    if x !== M.x
        copy!(M.x, x)
    end
    return
end

@inline function __update_grad!(M::MomentumDescent, g)
    if M.g !== g
        copy!(M.g, g)
    end
    return
end

"""
    NesterovMomentum

Descent method which minimizes the objective function in the direction
of antigradient at each step.
"""
mutable struct NesterovMomentum{T<:AbstractFloat,V<:AbstractVector{T}} <: OptBuffer
    x::V
    xpre::V
    g::V
    v::V
    α::T # learning rate
    β::T # decay rate
    y::T
    ypre::T
end

function NesterovMomentum(x::AbstractVector; learn_rate::Real=0.01, decay_rate::Real=0.9)
    F = float(eltype(x))

    # α and β are parameters of "half step" v <- βv - α∇f
    # applied twice, it must yield v <- decay_rate × v - learn_rate × ∇f
    α = learn_rate
    β = decay_rate
    return NesterovMomentum(
        similar(x, F),
        similar(x, F),
        similar(x, F),
        similar(x, F),
        convert(F, α),
        convert(F, β),
        F(NaN),
        F(NaN),
    )
end

function init!(optfn!, M::NesterovMomentum{T}, x0; kw...) where {T}
    optfn!(x0, zero(T), x0)
    fill!(M.v, zero(T))
    return M
end

@inline function reset!(M::NesterovMomentum)
    fill!(M.v, false)
    return M
end

function reset!(M::NesterovMomentum, x0, learn_rate=M.learn_rate, decay_rate=M.decay_rate)
    n = length(x0)
    if length(M.x) != n
        for v in (M.x, M.g, M.v)
            resize!(v, n)
        end
    end
    M.α = learn_rate / (1 + decay_rate)
    M.β = sqrt(decay_rate)
    return M
end

@inline function callfn!(fdf, M::NesterovMomentum, x, α, d)
    __update_arg!(M, x, α, d)
    y, g = fdf(M.x, M.g)
    __update_grad!(M, g)
    return y, g
end

function step!(optfn!, M::NesterovMomentum; constrain_step = infstep)
    α, β = M.α, M.β
    M.v *= β
    maxstep = constrain_step(M.x, M.v)
    if maxstep <= 1
        M.v *= maxstep / 2
    end
    optfn!(M.x, one(maxstep), M.v)
    M.v .-= α * M.g
    M.x .-= α * M.g
    return M.α
end

@inline function __update_arg!(M::NesterovMomentum, x, α, d)
    map!(M.x, d, x) do a, b
        muladd(α, a, b)
    end
    return
end

@inline function __update_arg!(M::NesterovMomentum, x)
    if x !== M.x
        copy!(M.x, x)
    end
    return
end

@inline function __update_grad!(M::NesterovMomentum, g)
    if M.g !== g
        copy!(M.g, g)
    end
    return
end


text:
module PoisonAgFEMTests

using IterativeSolvers: cg
using AlgebraicMultigrid

using STLCutters
using Gridap
using GridapEmbedded
using Test

using STLCutters: compute_stl_model
using STLCutters: read_stl, merge_nodes, get_bounding_box

# Manufactured solution
u(x) = x[1] + x[2] - x[3]
f(x) = - Δ(u)(x)
ud(x) = u(x)

stlpath = joinpath(@__DIR__,"data/Bunny-LowPoly.stl")

@time geo = STLGeometry( stlpath )

n = 28
δ = 0.2
pmin,pmax = get_bounding_box(geo)
diagonal = pmax-pmin
pmin = pmin - diagonal*δ
pmax = pmax + diagonal*δ
partition = (n,n,n)

bgmodel = CartesianDiscreteModel(pmin,pmax,partition)

# Cut the background model
cutter = STLCutter()
@time cutgeo,facet_to_inoutcut = cut(cutter,bgmodel,geo)

strategy = AggregateAllCutCells()
aggregates = aggregate(strategy,cutgeo,facet_to_inoutcut)

# Setup integration meshes
Ω_bg = Triangulation(bgmodel)
Ω = Triangulation(cutgeo,PHYSICAL)
Γd = EmbeddedBoundary(cutgeo)

# Setup normal vectors
n_Γd = get_normal_vector(Γd)

#writevtk(Ω,"trian_O")
#writevtk(Γd,"trian_Gd",cellfields=["normal"=>n_Γd])
#writevtk(Triangulation(bgmodel),"bgtrian")

# Setup Lebesgue measures
order = 1
degree = 2*order
dΩ = Measure(Ω,degree)
dΓd = Measure(Γd,degree)
#dΓg = Measure(Γg,degree)

vol = sum( ∫(1)*dΩ  )
surf = sum( ∫(1)*dΓd )
@show vol,surf

# Setup FESpace
Ω_act = Triangulation(cutgeo,ACTIVE)
Vstd = FESpace(Ω_act,ReferenceFE(lagrangian,Float64,order),conformity=:H1)

@time V = AgFEMSpace(Vstd,aggregates) #,Vser)
U = TrialFESpace(V)
# Weak form
γd = 10.0
h = (pmax - pmin)[1] / partition[1]

a(u,v) =
  ∫( ∇(v)⋅∇(u) ) * dΩ  +
  ∫( (γd/h)*v*u  - v*(n_Γd⋅∇(u)) - (n_Γd⋅∇(v))*u ) * dΓd

l(v) =
  ∫( v*f ) * dΩ +
  ∫( (γd/h)*v*ud - (n_Γd⋅∇(v))*ud ) * dΓd

# FE problem
@time op = AffineFEOperator(a,l,U,V)

A = get_matrix(op)
b = get_vector(op)

@time p = aspreconditioner(smoothed_aggregation(A))

@time  x = cg(A,b,verbose=true,Pl=p,reltol=1e-10)

@time uh = FEFunction(U,x)

@time e = u - uh

# Postprocess
l2(u) = sqrt(sum( ∫( u*u )*dΩ ))
h1(u) = sqrt(sum( ∫( u*u + ∇(u)⋅∇(u) )*dΩ ))

el2 = l2(e)
eh1 = h1(e)
ul2 = l2(uh)
uh1 = h1(uh)

#colors = color_aggregates(aggregates,bgmodel)
#writevtk(Ω_bg,"trian",celldata=["aggregate"=>aggregates,"color"=>colors],cellfields=["uh"=>uh])
#writevtk(Ω,"results",cellfields=["uh"=>uh])

@test el2/ul2 < 1.e-9
@test eh1/uh1 < 1.e-9

end # module


text:
"""
    time_loop(sys, state)

Time stepping procedure for the lattice Boltzmann state `state` given parameters `sys`
"""
function time_loop(sys::SysConst, state::LBM_state_2D; verbose=false)
    for t in 1:sys.Tmax
        if t % sys.tdump == 0
            mass = 0.0
            mass = sum(state.height)
            if verbose
                println("Time step $t mass is $(round(mass, digits=3))")
            end
        end
        Swalbe.filmpressure!(state, sys)
        Swalbe.h∇p!(state)
        Swalbe.slippage!(state, sys)
        state.Fx .= -state.h∇px .- state.slipx
        state.Fy .= -state.h∇py .- state.slipy
        Swalbe.equilibrium!(state)
        Swalbe.BGKandStream!(state)
        Swalbe.moments!(state)
    end
    return state
end
function time_loop(sys::SysConst, state::State, θ; verbose=false)
    for t in 1:sys.Tmax
        if t % sys.tdump == 0
            mass = 0.0
            mass = sum(state.height)
            if verbose
                println("Time step $t mass is $(round(mass, digits=3))")
            end
        end
        Swalbe.filmpressure!(state, sys, θ)
        Swalbe.h∇p!(state)
        Swalbe.slippage!(state, sys)
        state.Fx .= -state.h∇px .- state.slipx
        state.Fy .= -state.h∇py .- state.slipy
        Swalbe.equilibrium!(state)
        Swalbe.BGKandStream!(state)
        Swalbe.moments!(state)
    end
    return state
end

function time_loop(sys::SysConst, state::State, Δh::Vector; verbose=false)
    for t in 1:sys.Tmax
        if t % sys.tdump == 0
            mass = 0.0
            mass = sum(state.height)
            if verbose
                println("Time step $t mass is $(round(mass, digits=3))")
            end
        end
        push!(Δh, maximum(state.height) - minimum(state.height))
        Swalbe.filmpressure!(state, sys)
        Swalbe.h∇p!(state)
        Swalbe.slippage!(state, sys)
        state.Fx .= -state.h∇px .- state.slipx
        state.Fy .= -state.h∇py .- state.slipy
        Swalbe.equilibrium!(state, sys)
        Swalbe.BGKandStream!(state)
        Swalbe.moments!(state)
    end
    return state
end

function time_loop(sys::SysConst, state::State, f::Function, measure::Vector; verbose=false)
    for t in 1:sys.Tmax
        if t % sys.tdump == 0
            mass = 0.0
            mass = sum(state.height)
            if verbose
                println("Time step $t mass is $(round(mass, digits=3))")
            end
        end
        Swalbe.filmpressure!(state, sys)
        Swalbe.h∇p!(state)
        Swalbe.slippage!(state, sys)
        state.Fx .= -state.h∇px .- state.slipx
        state.Fy .= -state.h∇py .- state.slipy
        f(measure, state)
        Swalbe.equilibrium!(state)
        Swalbe.BGKandStream!(state)
        Swalbe.moments!(state)

    end
    return state, measure
end
# Time loop with snapshots saved
function time_loop(sys::SysConst, state::State, θ, data::Matrix; verbose=false)
    for t in 1:sys.Tmax
        if t % sys.tdump == 0
            mass = 0.0
            mass = sum(state.height)
            if verbose
                println("Time step $t mass is $(round(mass, digits=3))")
            end
        end
        Swalbe.filmpressure!(state, sys, θ)
        Swalbe.h∇p!(state)
        Swalbe.slippage!(state, sys)
        state.Fx .= -state.h∇px .- state.slipx
        state.Fy .= state.h∇py .+ state.slipy
        Swalbe.equilibrium!(state)
        Swalbe.BGKandStream!(state)
        Swalbe.moments!(state)

        Swalbe.snapshot!(data, state.height, t, dumping = sys.tdump)

    end
    return state
end

function time_loop(sys::SysConst_1D, state::State_1D; verbose=false)
    for t in 1:sys.Tmax
        if t % sys.tdump == 0
            mass = 0.0
            mass = sum(state.height)
            if verbose
                println("Time step $t mass is $(round(mass, digits=3))")
            end
        end
        Swalbe.filmpressure!(state, sys)
        Swalbe.h∇p!(state)
        Swalbe.slippage!(state, sys)
        state.F .= -state.h∇p .- state.slip
        Swalbe.equilibrium!(state)
        Swalbe.BGKandStream!(state)
        Swalbe.moments!(state)
    end
    return state
end

function time_loop(sys::SysConst_1D, state::State_1D, θ; verbose=false)
    for t in 1:sys.Tmax
        if t % sys.tdump == 0
            mass = 0.0
            mass = sum(state.height)
            if verbose
                println("Time step $t mass is $(round(mass, digits=3))")
            end
        end
        Swalbe.filmpressure!(state, sys, θ)
        Swalbe.h∇p!(state)
        Swalbe.slippage!(state, sys)
        state.F .= -state.h∇p .- state.slip
        Swalbe.equilibrium!(state)
        Swalbe.BGKandStream!(state)
        Swalbe.moments!(state)
    end
    return state
end

function time_loop(sys::SysConst_1D, state::State_1D, Δh::Vector; verbose=false)
    for t in 1:sys.Tmax
        if t % sys.tdump == 0
            mass = 0.0
            mass = sum(state.height)
            if verbose
                println("Time step $t mass is $(round(mass, digits=3))")
            end
        end
        push!(Δh, maximum(state.height) - minimum(state.height))
        Swalbe.filmpressure!(state, sys)
        Swalbe.h∇p!(state)
        Swalbe.slippage!(state, sys)
        state.F .= -state.h∇p .- state.slip
        Swalbe.equilibrium!(state, sys)
        Swalbe.BGKandStream!(state)
        Swalbe.moments!(state)
    end
    return state
end

function time_loop(sys::SysConst_1D, state::State_1D, f::Function, measure::Vector; verbose=false)
    for t in 1:sys.Tmax
        if t % sys.tdump == 0
            mass = 0.0
            mass = sum(state.height)
            if verbose
                println("Time step $t mass is $(round(mass, digits=3))")
            end
        end
        Swalbe.filmpressure!(state, sys)
        Swalbe.h∇p!(state)
        Swalbe.slippage!(state, sys)
        state.F .= -state.h∇p .- state.slip
        f(measure, state)
        Swalbe.equilibrium!(state)
        Swalbe.BGKandStream!(state)
        Swalbe.moments!(state)

    end
    return state, measure
end

function time_loop(sys::SysConst_1D, state::State_1D, data::Matrix; verbose=false)
    for t in 1:sys.Tmax
        if t % sys.tdump == 0
            mass = 0.0
            mass = sum(state.height)
            if verbose
                println("Time step $t mass is $(round(mass, digits=3))")
            end
        end
        Swalbe.filmpressure!(state, sys)
        Swalbe.h∇p!(state)
        Swalbe.slippage!(state, sys)
        state.F .= -state.h∇p .- state.slip
        Swalbe.equilibrium!(state)
        Swalbe.BGKandStream!(state)
        Swalbe.moments!(state)

        Swalbe.snapshot!(data, state.height, t, dumping = sys.tdump)

    end
    return data
end
"""
    run_flat(Sys::SysConst, device::String)

Performs a simulation of an flat interface without forces.

# Arguments

- `verbos :: Bool`: Enables consol output
- `T :: AbstractFloat`: Precision of output, default `Float64`

# Theory

Nothing at all should happen.
As the initial state is falt and no force is applied the fluid has no way to flow.
The equality
`` h(\\mathbf{x},0) = h(\\mathbf{x},\\infty), ``
should be satisfied for arbitrary many time steps.

# Examples
```julia
julia> using Swalbe, Test

julia> sys = Swalbe.SysConst(Lx=100, Ly=100, Tmax=5000);

julia> h = Swalbe.run_flat(sys, "CPU", verbos=false);

julia> @test all(h.height .== 1.0) # Check if all height values are identical to 1.0 (initial condition)
Test passed
```
"""
function run_flat(sys::SysConst, device::String; verbos=true)
    println("Simulating a flat interface without driving forces (nothing should happen) in two dimensions")
    state = Swalbe.Sys(sys, device)
    state.height .= 1.0
    time_loop(sys, state, verbose=verbos)

    return state.height
end

function run_flat(sys::SysConst_1D; verbos=true, T=Float64)
    println("Simulating a flat interface without driving forces (nothing should happen) in one dimension")
    state = Swalbe.Sys(sys)
    state.height .= 1.0
    time_loop(sys, state, verbose=verbos)

    return state.height
end

"""
    run_random(sys::SysConst, device::String)

Simulation of an random undulated interface

# Arguments
- `h₀ :: Float` : Average initial height
- `ϵ :: Float` : Amplitude of the flucutation
- `verbos :: Bool`: Enables consol output
- `T :: AbstractFloat`: Precision of output, default `Float64`

# Theory

Initial randomly perturbed fluid surface.
Unstable wavemodes should grow while wavemodes larger q₀ should be damped out.
Measuring this is on the TODO list.

# Examples
```julia
julia> using Swalbe, Test

julia> sys = Swalbe.SysConst(Lx=100, Ly=100, Tmax=5000);

julia> Swalbe.randinterface!(height, h₀, ϵ)

julia> h = Swalbe.run_random(sys, "CPU", h₀=10, ϵ=0.1, verbos=false);
```
"""
function run_random(sys::SysConst, device::String; h₀=1.0, ϵ=0.01, verbos=true)
    println("Simulating a random undulated interface in two dimensions")
    state = Swalbe.Sys(sys, device)
    Swalbe.randinterface!(state.height, h₀, ϵ)
    Swalbe.equilibrium!(state)
    time_loop(sys, state, 1/9, verbose=verbos)

    return state.height
end
# 1D case
function run_random(sys::SysConst_1D; h₀=1.0, ϵ=0.01, verbos=true)
    println("Simulating a random undulated interface in one dimension")
    state = Swalbe.Sys(sys)
    Swalbe.randinterface!(state.height, h₀, ϵ)
    Swalbe.equilibrium!(state)
    time_loop(sys, state, 1/9, verbose=verbos)
    return state.height
end

"""
    run_rayleightaylor(sys::SysConst, device::String)

Simulation of an random undulated interface and a gravitanional pull.

# Arguments
- `kx :: Int` : wavemode in x-direction, kx=18 -> 18 sine waves fitting into the domain
- `ky :: Int` : wavemode in y-direction
- `h₀ :: Float` : Average initial height
- `ϵ :: Float` : Amplitude of the flucutation
- `verbos :: Bool`: Enables consol output
- `T :: AbstractFloat`: Precision of output, default `Float64`

# Theory

Initial randomly perturbed fluid surface hanging from a substrate.
Here we have an interplay between gravity and surface tension.
The critical wavemode can be computed according to
`` q_0 =  ``
Measuring this is on the TODO list.

# Examples
```julia
julia> using Swalbe, Test

julia> sys = Swalbe.SysConst(Lx=100, Ly=100, Tmax=5000);

julia> Swalbe.randinterface!(height, h₀, ϵ)

julia> h = Swalbe.run_random(sys, "CPU", h₀=10, ϵ=0.1, verbos=false);
```
"""
function run_rayleightaylor(sys::SysConst, device::String; kx=15, ky=18, h₀=1.0, ϵ=0.001, verbos=true, T=Float64)
    println("Simulating the Rayleigh Taylor instability in two dimensions")
    state = Swalbe.Sys(sys, device)
    for i in 1:sys.Lx, j in 1:sys.Ly
        state.height[i,j] = h₀ * (1 + ϵ * sin(2π*kx*i/(sys.Lx-1)) * sin(2π*ky*j/(sys.Ly-1)))
    end
    diff = []
    Swalbe.equilibrium!(state)
    Swalbe.time_loop(sys, state, diff, verbose=verbos)
    return state.height, diff
end
# 1D case
function run_rayleightaylor(sys::SysConst_1D; k=15, h₀=1.0, ϵ=0.001, verbos=true, T=Float64)
    println("Simulating the Rayleigh Taylor instability in one dimension")
    state = Swalbe.Sys(sys)
    for i in 1:sys.L
        state.height[i] = h₀ * (1 + ϵ * sin(2π*k*i/(sys.L-1)))
    end
    diff = []
    Swalbe.equilibrium!(state)
    Swalbe.time_loop(sys, state, diff, verbose=verbos)
    return state.height, diff
end

"""
    run_dropletrelax()

Simulates an out of equilibrium droplet
"""
function run_dropletrelax(sys::SysConst, device::String; radius=20, θ₀=1/6, center=(sys.Lx÷2, sys.Ly÷2), verbos=true)
    println("Simulating an out of equilibrium droplet in two dimensions")
    state = Swalbe.Sys(sys, device)
    Swalbe.singledroplet(state.height, radius, θ₀, center)
    Swalbe.equilibrium!(state)
    area = []
    time_loop(sys, state, Swalbe.wetted!, area, verbose=verbos)
    return state.height, area
end
# 1D case
function run_dropletrelax(sys::SysConst_1D; radius=20, θ₀=1/6, center=(sys.L÷2), verbos=true )
    println("Simulating an out of equilibrium droplet in one dimensions")
    area = []
    state = Swalbe.Sys(sys)
    Swalbe.singledroplet(state.height, radius, θ₀, center)
    Swalbe.equilibrium!(state)
    time_loop(sys, state, Swalbe.wetted!, area, verbose=verbos)
    return state.height, area
end

"""
    run_dropletpatterned()

Simulates an droplet on a patterned substrate
"""
function run_dropletpatterned(sys::SysConst, device::String; radius=20, θ₀=1/6, center=(sys.Lx÷2, sys.Ly÷2), θₛ=fill(1/9, sys.Lx, sys.Ly), verbos=true)
    println("Simulating a droplet on a patterned substrate in two dimensions")
    state = Swalbe.Sys(sys, device)
    Swalbe.singledroplet(state.height, radius, θ₀, center)
    Swalbe.equilibrium!(state)
    time_loop(sys, state, θₛ, verbose=verbos)

    return state.height

end
# 1D case
function run_dropletpatterned(sys::SysConst_1D; radius=20, θ₀=1/6, center=sys.L÷2, θₛ=fill(1/2, sys.L), verbos=true)
    println("Simulating a droplet on a patterned substrate in one dimension")
    state = Swalbe.Sys(sys)
    Swalbe.singledroplet(state.height, radius, θ₀, center)
    Swalbe.equilibrium!(state)
    time_loop(sys, state, θₛ, verbose=verbos)
    return state.height
end

"""
    run_dropletforced()

Simulates an droplet on a patterned substrate
"""
function run_dropletforced(sys::SysConst, device::String; radius=20, θ₀=1/6, center=(sys.Lx÷2, sys.Ly÷2), fx=0.0, fy=0.0, verbos=true)
    bodyforce = zeros(2)
    bodyforce[1] = fx
    bodyforce[2] = fy
    println("Simulating a sliding droplet in two dimensions")
    state = Swalbe.Sys(sys, device)
    Swalbe.singledroplet(state.height, radius, θ₀, center)
    Swalbe.equilibrium!(state)
    println("Starting the lattice Boltzmann time loop")
    time_loop(sys, state, Swalbe.inclination!, bodyforce, verbose=verbos)

    return state.height, state.velx, state.vely

end
# 1D case
function run_dropletforced(sys::SysConst_1D; radius=20, θ₀=1/6, center=(sys.L÷2), θₛ=fill(1/9, sys.L), f=0.0, verbos=true)
    bodyforce = zeros(1)
    bodyforce[1] = f
    println("Simulating a sliding droplet in one dimension")
    state = Swalbe.Sys(sys)
    Swalbe.singledroplet(state.height, radius, θ₀, center)
    Swalbe.equilibrium!(state)
    println("Starting the lattice Boltzmann time loop")
    time_loop(sys, state, Swalbe.inclination!, bodyforce, verbose=verbos)
    return state.height, state.vel

end

# function run_active_thin_film(
#     sys::SysConst_1D;
#     h₀=1,
#     ϵ=0.001,
#     θ₀=1/9,
#     verbos=true,
#     rho = zeros(sys.L),
#     lap_rho = zeros(sys.L),
#     grad_rho = zeros(sys.L),
#     grad_h = zeros(sys.L),
#     lap_h = zeros(sys.L),
#     Gam = 0.5,
#     mobl = 1.0,
#     difu = 0.1,
#     T=Float64
# )
#     println("Simulating a active thin film")
#     fout, ftemp, feq, height, vel, pressure, dgrad, F, slip, h∇p, fluc = Swalbe.Sys(sys, true, T)
#     Swalbe.randinterface!(height, h₀, ϵ)
#     # Swalbe.sinewave1D!(height, h₀, 1, ϵ, 1)
#     rho .= 0.1
#     Swalbe.equilibrium!(feq, height, vel)
#     ftemp .= feq
#     println("Starting the lattice Boltzmann time loop")
#     for t in 1:sys.Tmax
#         if t % sys.tdump == 0
#             mass = 0.0
#             mass = sum(height)
#             if verbos
#                 println("Time step $t mass is $(round(mass, digits=3))")
#             end
#         end
#         Swalbe.update_rho(fluc, rho, dgrad, height, lap_rho, grad_rho, grad_h, lap_h, M=mobl, D=difu)
#         # Swalbe.filmpressure!(pressure, height, dgrad, sys.γ, θ₀, sys.n, sys.m, sys.hmin, sys.hcrit)
#         Swalbe.filmpressure!(pressure, height, dgrad, rho, sys.γ, θ₀, sys.n, sys.m, sys.hmin, sys.hcrit, Gamma=Gam)
#         Swalbe.∇f!(h∇p, pressure, dgrad, height)
#         Swalbe.slippage!(slip, height, vel, sys.δ, sys.μ)
#         # HereSwalbe.thermal!(fluc, height, sys.kbt, sys.μ, sys.δ)
#         # Here we a force that is like pull of an inclined plane
#         F .= h∇p .+ slip # .+ fluc
#         Swalbe.equilibrium!(feq, height, vel)
#         Swalbe.BGKandStream!(fout, feq, ftemp, -F)
#         Swalbe.moments!(height, vel, fout)
#     end

#     return height, vel, rho
#
# end