text:
module GtkTextUtils

    using Gtk, GtkExtensions, JuliaWordsUtils
    import Gtk: GtkTextIter, GLib.mutable, hasselection
    import JuliaWordsUtils: select_word_backward

    export select_word_double_click, select_word, get_line_text, cursor_position,
    get_text_iter_at_cursor, get_current_line_text,
    get_text_left_of_iter, get_text_right_of_iter,
    move_cursor_to_sentence_start, move_cursor_to_sentence_end,
    hasselection,replace_text

    function hasselection(b::GtkTextBuffer)
        (found,it_start,it_end) = selection_bounds(b)
        found
    end
    function replace_text(buffer::GtkTextBuffer,itstart::T,itend::T,str::AbstractString) where {T<:GtkTextIters}
        pos = offset(itstart)+1
        splice!(buffer,itstart:itend)
        insert!(buffer,GtkTextIter(buffer,pos),str)
    end

    cursor_position(b::GtkTextBuffer) = get_gtk_property(b,:cursor_position,Int)

    get_text_iter_at_cursor(b::GtkTextBuffer) =
        GtkTextIter(b,cursor_position(b)+1) #+1 because there's a -1 in gtk.jl

    function get_current_line_text(buffer::GtkTextBuffer)
        it = get_text_iter_at_cursor(buffer)
        return get_line_text(buffer,it)
    end

    function select_word(it::GtkTextIter,buffer::GtkTextBuffer,stop_at_dot::Bool)

        (txt, line_start, line_end) = get_line_text(buffer,it)

        pos = offset(it) - offset(line_start) +1#not sure about the +1 but it feels better
        if pos <= 0
            return ("",GtkTextIter(buffer,offset(it)),
            GtkTextIter(buffer,offset(it)))
        end

        word,i,j = extend_word(pos, txt, stop_at_dot)

        its = GtkTextIter(buffer, i + offset(line_start) )
        ite = GtkTextIter(buffer, j + offset(line_start) + 1)

        return (word,its,ite)
    end
    select_word(it::GtkTextIter,buffer::GtkTextBuffer) = select_word(it,buffer,true)

    function select_word_backward(it::GtkTextIter,buffer::GtkTextBuffer,stop_at_dot::Bool)

        (txt, line_start, line_end) = get_line_text(buffer,it)
        pos = offset(it) - offset(line_start) #position of cursor in txt

        if pos <= 0 || length(txt) == 0
            return ("",GtkTextIter(buffer,offset(it)),
            GtkTextIter(buffer,offset(it)))
        end

        txt = CharArray(txt,pos)
        (i,j) = select_word_backward(pos,txt,stop_at_dot)

        its = GtkTextIter(buffer, i + offset(line_start) )
        ite = GtkTextIter(buffer, offset(it))

        return (txt[i:j],its,it)
    end

    function get_line_text(buffer::GtkTextBuffer, it::GtkTextIter)

        itstart, itend = mutable(it), mutable(it)
        li = get_gtk_property(itstart, :line, Integer)

        skip(itstart, :backward_line) #seems there's no skip to line start
        li != get_gtk_property(itstart, :line, Integer) && skip(itstart, 1, :line)#for fist line
        !get_gtk_property(itend, :ends_line, Bool) && skip(itend, :forward_to_line_end)

        return ( (itstart:itend).text[String], itstart, itend)
    end

    function select_word_double_click(textview::GtkTextView,buffer::GtkTextBuffer,x::Integer,y::Integer)

        (x,y) = text_view_window_to_buffer_coords(textview,x,y)
        iter_end = get_iter_at_position(textview,x,y)
        #iter_end = mutable( get_text_iter_at_cursor(buffer) ) #not using this because the cursor position is modified somewhere

        (w, iter_start, iter_end) = select_word(iter_end,buffer)
        select_range(buffer,iter_start,iter_end)
    end

    function get_text_right_of_cursor(buffer::GtkTextBuffer)
        it = mutable(get_text_iter_at_cursor(buffer))
        return (it:(it+1)).text[String]
    end
    function get_text_left_of_cursor(buffer::GtkTextBuffer)
        it = mutable(get_text_iter_at_cursor(buffer))
        return ((it-1):it).text[String]
    end

    get_text_left_of_iter(it::MutableGtkTextIter) = ((it-1):it).text[String]
    get_text_right_of_iter(it::MutableGtkTextIter) = (it:(it+1)).text[String]

    get_text_left_of_iter(it::GtkTextIter) = ((mutable(it)-1):mutable(it)).text[String]
    get_text_right_of_iter(it::GtkTextIter) = (mutable(it):(mutable(it)+1)).text[String]

    nlines(it_start,it_end) = abs(line(it_end)-line(it_start))+1

    function move_cursor_to_sentence_start(buffer::GtkTextBuffer)
        it = mutable( get_text_iter_at_cursor(buffer) )
        skip(it, :backward_sentence_start)
        place_cursor(buffer,it)
    end
    function move_cursor_to_sentence_end(buffer::GtkTextBuffer)
        it = mutable( get_text_iter_at_cursor(buffer) )
        skip(it, :forward_sentence_end)
        place_cursor(buffer,it)
    end

    function toggle_wrap_mode(v::GtkTextView)
        wm = get_gtk_property(v, :wrap_mode, Int)
        wm = convert(Bool, wm)
        setproperty!(v, :wrap_mode, !wm)
        nothing
    end

end # module


text:
@testset "Interval arithmetic" begin
    a = -1.5 ± 0.5
    b = -1..1
    @test a + b == -3.0..0.0
    @test a * b == -2.0..2.0
    @test a * b + a == -4.0..1.0
    @test a * (b + 1) == -4.0..0.0
end

@testset "Interval matrix arithmetic" begin
    a = 1.0..1.3; b = 2.0..3.5; c = -0.5 ± 0.5; d = 0.0 ± 0.1
    a₊ = 2.0..2.6; b₊ = 4.0..7.0; c₊ = -2.0..0.0; d₊ = -0.2..0.2
    a₋ = -0.3..0.3; b₋ = -1.5..1.5; c₋ = -1.0..1.0; d₋ = -0.2..0.2
    A = IntervalMatrix([a b; c d])

    B = A + A
    @test B isa IntervalMatrix && B == IntervalMatrix([a₊ b₊; c₊ d₊])

    B = A - A
    @test B isa IntervalMatrix && B == IntervalMatrix([a₋ b₋; c₋ d₋])

    # multiply interval and interval matrix
    x = 0.0..2.0
    @test x * A == A * x ==
          IntervalMatrix([0.0..2.6 0.0..7.0; -2.0..0.0 -0.2..0.2])
    # multiply scalar and interval matrix
    x = 1.0
    for A2 in [x * A, A * x, A / x]
        @test A2 == A && typeof(A2) == typeof(A)
    end

    # arithmetic closure using interval matrices and non-interval matrices
    Ainf = inf(A)
    @test A + Ainf isa IntervalMatrix
    @test Ainf + A isa IntervalMatrix
    @test A - Ainf isa IntervalMatrix
    @test Ainf - A isa IntervalMatrix
    @test A * Ainf isa IntervalMatrix
    @test Ainf * A isa IntervalMatrix
end

@testset "Matrix multiplication" begin

    # test default settings
    @test get_multiplication_mode() == Dict(:multiplication => :fast)

    A = IntervalMatrix([2..4 -2..1; -1..2 2..4])
    set_multiplication_mode(:slow)
    @test A * A == IntervalMatrix([0..18 -16..8; -8..16 0..18])
    @test A * mid.(A) == IntervalMatrix([5..12.5 -8..2; -2..8 5..12.5])
    @test mid.(A) * A == IntervalMatrix([5..12.5 -8..2; -2..8 5..12.5])

    # set_multiplication_mode(:rank1)
    # @test A * A == [0..18 -16..8; -8..16 0..18]

    set_multiplication_mode(:fast)
    @test A * A == IntervalMatrix([-2..19.5 -16..10; -10..16 -2..19.5])
    @test A * mid.(A) == IntervalMatrix([5..12.5 -8..2; -2..8 5..12.5])
    @test mid.(A) * A == IntervalMatrix([5..12.5 -8..2; -2..8 5..12.5])
end

ext: jl


text:
module test_julia_dates

using Test
using Dates # DateTime DateFormat DatePeriod TimePeriod year month day hour minute second

const RFC3339 = DateFormat("yyyy-mm-dd\\THH:MM:SSZ")

created_at = "2019-06-15T09:49:47Z"
dt = parse(DateTime, created_at, RFC3339)
@test Dates.format(dt, RFC3339) == created_at

@test year(dt) == 2019
@test month(dt) == 6
@test day(dt) == 15
@test hour(dt) == 9
@test minute(dt) == 49
@test second(dt) == 47

@test Year(2019) isa DatePeriod
@test Hour(10) isa TimePeriod

@test dayname(dt) == dayname(dayofweek(dt)) == "Saturday"
@test join(dayname.(Dates.Mon:Dates.Sun), ", ") == "Monday, Tuesday, Wednesday, Thursday, Friday, Saturday, Sunday"

@test length(January:December) == 12

@test Date(0) == Date(0, 1, 1)
@test Date(1) == Date(1, 1, 1)
@test Time(0) == Time(0, 0, 0)
@test Time(1) == Time(1, 0, 0)

end # module test_julia_dates


text:
using Variography
using Meshes
using GeoStatsBase
using GeoStatsImages
using Unitful
using Distances
using LinearAlgebra
using CSV, DelimitedFiles
using Plots; gr(size=(600,400))
using ReferenceTests, ImageIO
using Test, Random
import CoDa: Composition

# workaround for GR warnings
ENV["GKSwstype"] = "100"

# environment settings
isCI = "CI" ∈ keys(ENV)
islinux = Sys.islinux()
visualtests = !isCI || (isCI && islinux)
datadir = joinpath(@__DIR__,"data")

# list of tests
testfiles = [
  "empirical.jl",
  "partition.jl",
  "varioplane.jl",
  "theoretical.jl",
  "nesting.jl",
  "sampling.jl",
  "pairwise.jl",
  "fitting.jl"
]

@testset "Variography.jl" begin
  for testfile in testfiles
    include(testfile)
  end
end


text:
"""
    Git(;
        ignore=String[],
        name=nothing,
        email=nothing,
        branch=nothing,
        ssh=false,
        jl=true,
        manifest=false,
        gpgsign=false,
    )

Creates a Git repository and a `.gitignore` file.

## Keyword Arguments
- `ignore::Vector{<:AbstractString}`: Patterns to add to the `.gitignore`.
  See also: [`gitignore`](@ref).
- `name::AbstractString`: Your real name, if you have not set `user.name` with Git.
- `email::AbstractString`: Your email address, if you have not set `user.email` with Git.
- `branch::AbstractString`: The desired name of the repository's default branch.
- `ssh::Bool`: Whether or not to use SSH for the remote.
  If left unset, HTTPS is used.
- `jl::Bool`: Whether or not to add a `.jl` suffix to the remote URL.
- `manifest::Bool`: Whether or not to commit `Manifest.toml`.
- `gpgsign::Bool`: Whether or not to sign commits with your GPG key.
  This option requires that the Git CLI is installed,
  and for you to have a GPG key associated with your committer identity.
"""
@plugin struct Git <: Plugin
    ignore::Vector{String} = String[]
    name::Union{String, Nothing} = nothing
    email::Union{String, Nothing} = nothing
    branch::Union{String, Nothing} = nothing
    ssh::Bool = false
    jl::Bool = true
    manifest::Bool = false
    gpgsign::Bool = false
end

# Try to make sure that no files are created after we commit.
priority(::Git, ::typeof(posthook)) = 5
gitignore(p::Git) = p.ignore

function validate(p::Git, t::Template)
    if p.gpgsign && !git_is_installed()
        throw(ArgumentError("Git: gpgsign is set but the Git CLI is not installed"))
    end

    foreach((:name, :email)) do k
        user_k = "user.$k"
        if getproperty(p, k) === nothing && isempty(LibGit2.getconfig(user_k, ""))
            throw(ArgumentError("Git: Global Git config is missing required value '$user_k'"))
        end
    end
end

# Set up the Git repository.
function prehook(p::Git, t::Template, pkg_dir::AbstractString)
    LibGit2.with(LibGit2.init(pkg_dir)) do repo
        LibGit2.with(GitConfig(repo)) do config
            foreach((:name, :email)) do k
                v = getproperty(p, k)
                v === nothing || LibGit2.set!(config, "user.$k", v)
            end
        end
        commit(p, repo, pkg_dir, "Initial commit")
        pkg = basename(pkg_dir)
        suffix = p.jl ? ".jl" : ""
        url = if p.ssh
            "git@$(t.host):$(t.user)/$pkg$suffix.git"
        else
            "https://$(t.host)/$(t.user)/$pkg$suffix"
        end
        default = LibGit2.branch(repo)
        branch = something(p.branch, default)
        if branch != default
            LibGit2.branch!(repo, branch)
            delete_branch(GitReference(repo, "refs/heads/$default"))
        end
        close(GitRemote(repo, "origin", url))
    end
end

# Create the .gitignore.
function hook(p::Git, t::Template, pkg_dir::AbstractString)
    ignore = mapreduce(gitignore, vcat, t.plugins)
    # Only ignore manifests at the repo root.
    p.manifest || "Manifest.toml" in ignore || push!(ignore, "/Manifest.toml")
    unique!(sort!(ignore))
    gen_file(joinpath(pkg_dir, ".gitignore"), join(ignore, "\n"))
end

# Commit the files.
function posthook(p::Git, ::Template, pkg_dir::AbstractString)
    # Special case for issue 211.
    if Sys.iswindows()
        files = filter(f -> startswith(f, "_git2_"), readdir(pkg_dir))
        foreach(f -> rm(joinpath(pkg_dir, f)), files)
    end

    # Ensure that the manifest exists if it's going to be committed.
    manifest = joinpath(pkg_dir, "Manifest.toml")
    if p.manifest && !isfile(manifest)
        touch(manifest)
        with_project(Pkg.update, pkg_dir)
    end

    LibGit2.with(GitRepo(pkg_dir)) do repo
        LibGit2.add!(repo, ".")
        msg = "Files generated by PkgTemplates"
        v = version_of("PkgTemplates")
        v === nothing || (msg *= "\n\nPkgTemplates version: $v")
        # TODO: Put the template config in the message too?
        commit(p, repo, pkg_dir, msg)
    end
end

function commit(p::Git, repo::GitRepo, pkg_dir::AbstractString, msg::AbstractString)
    if p.gpgsign
        run(pipeline(`git -C $pkg_dir commit -S --allow-empty -m $msg`; stdout=devnull))
    else
        LibGit2.commit(repo, msg)
    end
end

needs_username(::Git) = true

function git_is_installed()
    return try
        run(pipeline(`git --version`; stdout=devnull))
        true
    catch
        false
    end
end

if isdefined(Pkg, :dependencies)
    function version_of(pkg::AbstractString)
        for p in values(Pkg.dependencies())
            p.name == pkg && return p.version
        end
        return nothing
    end
else
    version_of(pkg::AbstractString) = get(Pkg.installed(), pkg, nothing)
end


text:
# The following iterative sequence is defined for the set of positive integers:
#
# n → n/2 (n is even)
# n → 3n + 1 (n is odd)
#
# Using the rule above and starting with 13, we generate the following sequence:
#
# 13 → 40 → 20 → 10 → 5 → 16 → 8 → 4 → 2 → 1
#
# It can be seen that this sequence (starting at 13 and finishing at 1) contains
# 10 terms. Although it has not been proved yet (Collatz Problem), it is thought
# that all starting numbers finish at 1.
#
# Which starting number, under one million, produces the longest chain?
#
# NOTE: Once the chain starts the terms are allowed to go above one million.

using ProjectEulerSolutions

# Memoization type solution to store intermediate results without recursion.
# It does end up recalculating some solutions multiple times.
function p014solution_memoize(num::Integer=3)::Integer
    cache = Dict{Integer, Integer}(1 => 1)
    for k = 2:num
        n = k
        i = 1
        while !haskey(cache, n)
            if n % 2 == 0
                n = fld(n, 2)
            else
                n = 3 * n + 1
            end
            i += 1
        end
        @inbounds cache[k] = get(cache, n, 0) + i - 1
    end

    # Find max value
    return maxdictkey(cache)
end

# Return maximum value in dictionary
function maxdictkey(cache::Dict{Integer, Integer})::Integer
    # Find max value
    maxkey, maxvalue = 0, 0
    for (key, value) in cache
        if value > maxvalue
            maxkey, maxvalue = key, value
        end
    end
    return maxkey
end

# Dynamic programming type solution to store intermediate results
function collatz!(cache::Dict{Integer, Integer}, n::Integer)::Integer
    if haskey(cache, n)
        return cache[n]
    else
        if n % 2 == 0
            @inbounds cache[n] = collatz!(cache, fld(n, 2)) + 1
        else
            @inbounds cache[n] = collatz!(cache, 3 * n + 1) + 1
        end
        return cache[n]
    end
end

# Same type of answer, but with no caching.  Sometimes the "smartest" solution
# isn't always the fastest.
function p014solution_nocache(num::Integer=3)::Integer
    longest = 0
    max_terms = 0
    for i in 1:num
        n = i
        nterms = 1
        while n != 1
            nterms += 1
            if nterms > max_terms
                max_terms = nterms
                longest = i
            end
            if n % 2 == 0
                n = fld(n, 2)
            else
                n = 3 * n + 1
            end
        end
    end

    return longest
end

# Recursive type solution to store intermediate results.  Theoretically
# recalculates no numbers, but recursion adds an overhead.
function p014solution_recurse(num::Integer=3)::Integer
    cache = Dict{Integer, Integer}(1 => 1)
    for k = num:-1:2
        if !haskey(cache, k)
            collatz!(cache, k)
        end
    end

    # Find max value
    return maxdictkey(cache)
end

p014 = Problems.Problem(Dict("Recurse" => p014solution_recurse,
                             "Memoize" => p014solution_memoize,
                             "No Cache" => p014solution_nocache))

Problems.benchmark(p014, 1_000_000)

max_stars_repo_name: gnujosh/julia-euler-project

max_issues_repo_name: gnujosh/julia-euler-project

max_forks_repo_name: gnujosh/julia-euler-project


text:
const _smallbuf_size = 256
#@show _smallbuf_size

mutable struct SpaceTokenizer{T}
    io::T
    buf::Array{UInt8,1}
    smallbuf::Array{UInt8,1}
    smallbufpos::Int
    smallbuflen::Int

    SpaceTokenizer(stream::T) where {T <: IO} = EachToken(stream, 2^10)
    SpaceTokenizer(stream::T, maxbuf::Int) where {T <: IO} = begin
        smallbuf = Array{UInt8,1}(_smallbuf_size)
        smallbuflen = readbytes!(stream, smallbuf, _smallbuf_size)
        new{T}(stream, Array{UInt8,1}(maxbuf), smallbuf, 1, smallbuflen)
    end
end

"""
This could be executed after the underyling IOBuffer has been
reset.
"""
function reset(itr::SpaceTokenizer)
  itr.smallbuflen = readbytes!(itr.io, itr.smallbuf, _smallbuf_size)
  itr.smallbufpos = 1
end

@inline function _read_spaces(itr)
    #@show String(smallbuf), smallbufpos, smallbuflen
    @inbounds while itr.smallbuflen > 0
        if isspacecode(itr.smallbuf[itr.smallbufpos])
            itr.smallbufpos += 1 # move the position
            if itr.smallbufpos > itr.smallbuflen
                # need to refill the buffer
                itr.smallbuflen = readbytes!(itr.io, itr.smallbuf, _smallbuf_size)
                itr.smallbufpos = 1
            end
        else
            break
        end
    end
    return itr.smallbuflen > 0
end

@inline function _buffer_nonspaces(itr)
    curbuf = 0
    @inbounds while itr.smallbuflen > 0
        if isspacecode(itr.smallbuf[itr.smallbufpos])
            break
        else
            curbuf += 1
            if curbuf > length(itr.buf)
                throw(ArgumentError("tokensize exceeded buffer"))
            end
            itr.buf[curbuf] = itr.smallbuf[itr.smallbufpos] # save the current value
            itr.smallbufpos += 1
            if itr.smallbufpos > itr.smallbuflen
                # need to refill the buffer
                itr.smallbuflen = readbytes!(itr.io, itr.smallbuf, _smallbuf_size)
                itr.smallbufpos = 1
            end
        end
    end
    return curbuf
end

@inline function step!(itr::SpaceTokenizer)
    # The step function has to read through spaces
    # then store non-spaces into buf
    # Then read through the rest of the spaces
    if _read_spaces(itr)
        return _buffer_nonspaces(itr)
    else
        return -1
    end
end

@inline function end_of_stream(itr::SpaceTokenizer)
    return itr.smallbuflen == 0
end

@inline function next(::Type{T}, itr::SpaceTokenizer) where {T <: Union{Float16,Float32,Float64,Integer}}
    if (curlen = step!(itr)) == -1
        return zero(T)
    else
        return myparse(T, itr.buf, 1, curlen)
    end
end


text:
# * DSL for calculation of quantities

# The somewhat complicated setup with walking the expression tree in
# =make_field= (twice implemented :) allows for very clean expression
# of the formulaic dependencies between different quantities, almost
# as if it were pure math. All ingoing quantities are either unitful,
# or made unitful using the set base units. Furthermore, outgoing
# quantities are transformed to the base units, even though the
# expression may result in a different (but equivalent) unit
# expression. This way, even if a period time is provided in
# femtoseconds, the wavenumber will always be returned in Kaysers, for
# instance.

# ** Testing presence of quantity in the namespace
# These methods are used for the =if x= constructs in the quantity
# conversion DSL. We return =true= if the symbol =x= is present in
# the =params= dict. Since we implement it as an expression walking
# algorithm, we can have constructs as =if x || y=, which will
# expand to =if :x in keys(params) || :y in keys(params)=. At first,
# however, we check if the symbols =isdefined= in the =Main= scope,
# which is generally the case for operators and such, to avoid
# generating tentative key lookups for every symbol that is
# encountered.

test_symbol_walk(node, params) = node
test_symbol_walk(node::Symbol, params) =
    isdefined(Main, node) ? node : Expr(:call, :in, Expr(:quote, node),
                                        Expr(:call, :keys, params))
test_symbol_walk(node::Expr, params) =
    Expr(node.head, test_symbol_walk.(node.args, Ref(params))...)#

# ** Expression walker
# This walks the quantity conversions block, replacing symbols with
# references to dictionary items and converting quantities to
# applicable base units, before assignment. As above, before
# generating dictionary lookups for all symbols encountered, we
# first check if the symbol =isdefined= in the =Main= scope.

get_reference(r::Symbol,params) = Expr(:ref, params, Expr(:quote, r))
get_reference(r,params) = r

function get_symbol(s::Symbol, params, escape)
    if isdefined(Main, s)
        escape ? esc(s) : s
    else
        test_expression = test_symbol_walk(s, params)
        ref_expression = get_reference(s, params)
        Expr(:if, test_expression, ref_expression, escape ? esc(s) : s)
    end
end
get_symbol(s, params, escape) = s

walk(node, params, escape) = node
walk(node::Symbol, params, escape) = get_symbol(node, params, escape)

function walk(node::Expr, params, escape)
    if node.head ∈ [:line, :quote, :macrocall]
        node
    elseif node.head ∈ [:if, :elseif]
        # Dispatch "if x" to test_symbol_walk, which checks if :x is
        # present as key in params.
        args = Any[test_symbol_walk(node.args[1], params)]
        append!(args, [walk(a, params, escape) for a in node.args[2:end]])
        Expr(node.head, args...)
    else
        if node.head == :(=)
            target,expr = node.args
            r = get_reference(target, params)
            result = walk(expr, params, escape)
            Expr(node.head, r, result)
        else
            args = walk.(node.args, Ref(params), Ref(escape))
            Expr(node.head, args...)
        end
    end
end

# ** Namespace macro

"""
    @namespace!(exprs, params)

This macro uses the dictionary `params` as a "namespace", i.e. all
symbols are assumed to be keys in this dictionary. We need to escape
the generated expression tree to modify in the scope of the caller,
not the global scope.
"""
macro namespace!(exprs, params)
    local tree = walk(exprs, params, false)
    quote
        $(esc(tree))()
    end
end

# ** Test of "competing quantities"
"""
    test_field_parameters(field_params, set)

This function ensures that one and only one of \"competing\"
quantities is specified."""
function test_field_parameters(field_params, set)
    info = set ∩ keys(field_params)
    set_string = join(set, ", ", " and ")

    length(info) == 0 &&
        (length(set) > 1 && error("Need to provide one of $(set_string)") ||
         error("Required parameter $(set_string) missing"))
    length(info) > 1 &&
        error("Can only specify one of $(set_string)")

    info
end


text:
"""
    StackVector([itr])

Construct a `StackVector` containing up to 64 `Bool` values.
A `StackVector` is stored as an integer in memory and is immutable.
"""
struct StackVector <: AbstractVector{Bool}
    x::UInt
    len::Int

    function StackVector(x::UInt, len::Int, ::Unsafe)
        new(x, len)
    end
end

function StackVector(x::UInt, len::Int)
    (len % UInt) > Sys.WORD_SIZE && throw(DomainError(len, "len must be 0:$(Sys.WORD_SIZE)"))
    return StackVector(x & mask(len), len, unsafe)
end

mask(L) = ifelse(L == 64, typemax(UInt), UInt(1) << (L & 63) - 1)
Base.size(s::StackVector) = (length(s),)
Base.length(s::StackVector) = s.len

function Base.hash(x::StackVector, h::UInt)
    base = (0x93774c8a392b33bb % UInt) * length(x)
    return hash(x.x, h ⊻ base)
end

StackVector() = StackVector(UInt(0), 0, unsafe)

@noinline function throw_stackvec_err()
    throw(DomainError("StackVector can only contain $(Sys.WORD_SIZE) elements"))
end

function StackVector(itr)
    bits = zero(UInt)
    len = 0
    for i in itr
        len += 1
        len > Sys.WORD_SIZE && throw_stackvec_err()
        val = convert(UInt, convert(Bool, i))
        bits |= (val << ((len-1) & 63))
    end
    return StackVector(bits, len, unsafe)
end

function Base.getindex(s::StackVector, i::Int)
    @boundscheck checkbounds(s, i)
    return isodd(s.x >>> unsigned(i-1))
end

"""
    setindex(collection, v, i::Int)

Return a copy of `collection` with the value at index `i` set to `v`.

# Examples

```jldoctest
julia> x = StackVector([true, false])
2-element StackVector:
 1
 0

julia> setindex(x, false, 1)
2-element StackVector:
 0
 0
```
"""
function setindex(s::StackVector, v::Bool, i::Int)
    @boundscheck checkbounds(s, i)
    u = UInt(1) << ((i-1) & 63)
    typeof(s)(ifelse(v, s.x | u, s.x & ~u), s.len, unsafe)
end

function push(s::StackVector, v::Bool, ::Unsafe)
    return StackVector(s.x | UInt(v) << (length(s) & 63), length(s)+1, unsafe)
end

function push(s::StackVector, v)
    length(s) == 64 && throw_stackvec_err()
    v_ = convert(Bool, v)
    return push(s, v_, unsafe)
end

function pop(s::StackVector, ::Unsafe)
    return StackVector(s.x & ~(UInt(1) << (length(s) & 63)), length(s) - 1, unsafe)
end

pop(s::StackVector) = isempty(s) ? throw_empty_err() : pop(s, unsafe)

function Base.iterate(s::StackVector, i::Int=0)
    i+1 > length(s) && return nothing
    isodd(s.x >>> (i&63)), i+1
end

Base.in(v::Bool, s::StackVector) = !iszero(ifelse(v, s.x, s.x ⊻ mask(length(s))))

Base.maximum(s::StackVector) = !minimum(~s)
function Base.minimum(s::StackVector)
    isempty(s) && throw_empty_err()
    return s.x == mask(length(s))
end

Base.sum(s::StackVector) = count_ones(s.x)

function Base.convert(::Type{BitVector}, s::StackVector)
    b = BitVector(undef, length(s))
    !isempty(s) && @inbounds b.chunks[1] = s.x
    b
end

Base.:~(s::StackVector) = StackVector(s.x ⊻ mask(length(s)), s.len, unsafe)

function Base.findfirst(s::StackVector)
    isempty(s) && return nothing
    i = trailing_zeros(s.x) + 1
    i > length(s) && return nothing
    return i
end

function Base.findfirst(f::Function, s::StackVector)
    isempty(s) && return nothing
    ft::Bool = f(true)
    ff::Bool = f(false)
    ft & ff && return 1
    !(ft | ff) && return nothing
    return findfirst(ifelse(ft, s, ~s))
end

Base.argmin(s::StackVector) = argmax(~s)
function Base.argmax(s::StackVector)
    isempty(s) && throw_empty_err()
    f = findfirst(s)
    return f === nothing ? 1 : f
end

function Base.reverse(s::StackVector)
    x = s.x
    x = ((x & 0xaaaaaaaaaaaaaaaa) >>> 1)  | ((x & 0x5555555555555555) << 1)
    x = ((x & 0xcccccccccccccccc) >>> 2)  | ((x & 0x3333333333333333) << 2)
    x = ((x & 0xf0f0f0f0f0f0f0f0) >>> 4)  | ((x & 0x0f0f0f0f0f0f0f0f) << 4)
    x = bswap(x)
    x >>>= sizeof(UInt) << 3 - length(s)
    typeof(s)(x, s.len, unsafe)
end

function Base.circshift(s::StackVector, k::Int)
    isempty(s) && return s
    shift = k % length(s)
    left = ifelse(k < 0, length(s)+shift, shift) & 63
    right = (length(s) - left) & 63
    bits = ((s.x << left) | (s.x >>> right)) & mask(length(s))
    return StackVector(bits, length(s), unsafe)
end


text:
#} functions for parsing Standard Algebraic Notation (SAN) and converting
# movetext to a array of Move objects

import Base: start, done, next

# Examples:
# e4 exd4 d8=Q d8=Q+ d8=Q#
# Ne4 Nfe4 Nxe4 Nfxe4 Nxe4+ Nxe4# Nfxe4+ Nfxe4#
# O-O O-O-O O-O-O+ O-O+ O-O#  O-O-O#

PIECE_CODES = Dict("P"=>1, "N"=>2, "B" =>3, "R" =>4, "Q"=>5, "K"=>6, "O-O"=>7, "O-O-O"=>8)

UTF8_PIECE_SYMBOLS = Dict(
    'K' => '♔', 'Q' => '♕', 'R' => '♖', 'B' => '♗', 'N' => '♘', 'P' => '♙',
    'k' => '♚', 'q' => '♛', 'r' => '♜', 'b' => '♝', 'n' => '♞', 'p' => '♟')

COLOR_WHITE = true
COLOR_BLACK = false

type Move
  # bitfield: [ color | piece piece piece | iscap | ispromotion | iscastle | ischeck | ismate]
  color::Bool
  capture::Bool
  check::Bool
  mate::Bool
  piece::Int8
  promotion::Int8  # piece code of promotion, or 0 if not
  fromsquare::Int8
  tosquare::Int8
  #comment::String
end

MoveList = Array{Move,1}
MoveText = String

function square2subs(square::String)
  file = Int(square[1]) - Int('a') + 1
  rank = Int(square[2]) - Int('1') + 1
  (rank, file)
end
function index(square::String)
  # linear index of a square given in algebraic notation
  rank, file = square2subs(square)
  8*file + rank
end

capture(m::SubString) = 'x' in m
capture(m::Move) = m.capture
promotion(m::String) = '=' in m
promotion(m::Move) = m.ispromotion
check(m::String) = '+' in m
check(m::Move) = m.ischeck
mate(m::String) = '#' in m
mate(m::Move) = m.ismate

function chompspecial!(m::SubString, c::Char)
  if m[end] == c
    m = m[1:end-1]
    return true
  else
    return false
  end
end

function chomppiece!(m::SubString)
  if m[1] == 'O'
    if length(m) >= 4 && m[4] == '-'
      m = m[6:end]
      return PIECE_CODES["O-O-O"]
    else
      m = m[4:end]
      return PIECE_CODES["O-O"]
    end
  elseif isupper(m[1])
    p = m[1]
    m = m[2:end]
    return PIECE_CODES[string(p)]
  else
    return PIECE_CODES["P"]
  end
end

function chomppromotion!(m::SubString)
  # assumes no check or mate at end anymore
  if m[end-1] == '='
    c = PIECE_CODES[m[end]]
    m = m[1:end-2]
    return c
  else
    return 0
  end
end

"""
tosquare(s)

Parse string `s` and return the destination square.
"""
function tosquare(m::String)

end

"""
fromsquare(s)

Parse string `s` and return the origin square.
"""
function fromsquare(m::String)

end

function parse(movetext::String)
  moves = Array{Move}(1)
  col = COLOR_WHITE
  for m in split(movetext)
    code= chomppiece!(m)
    ma = chompspecial!(m, '#')
    ch = chompspecial!(m, '#')
    prom = chomppromotion!(m)
    cap = capture(m)
    to = 1
    fr = 1
    push!(moves, Move(col,cap,ch,ma,code,prom,to,fr))
    col = !col
  end
  return moves
end

println(parse("d4 Nf6 c4 g6"))


text:
@testset "Axiom Network Tests" begin

    axiom_network = ALifeBenchmark.get_axiom_node()

    @test axiom_network.string == "001"
    @test axiom_network.out_excitatory[1].string == "000"
    @test axiom_network.out_excitatory[1].out_excitatory[1].string == "01"

    @test axiom_network.out_excitatory[1].in_excitatory[1] == axiom_network
    @test axiom_network.out_excitatory[1].out_excitatory[1].in_excitatory[1] == axiom_network.out_excitatory[1]

end

@testset "Development Helper Test" begin

    node = ALifeBenchmark.Node("1", [], [], [], [ALifeBenchmark.Node("0")])
    rule = ALifeBenchmark.Rule("1", "11", "", false, false, false, false, false, false, 0, 0)

    ALifeBenchmark.apply_to_all(ALifeBenchmark.fill_temp_links!, node)
    ALifeBenchmark._apply_rule!(rule, node)
    ALifeBenchmark.apply_to_all(ALifeBenchmark.apply_temp_links!, node)

    @test node.string == "11"
    @test node.out_excitatory[1].string == "0"

    node = ALifeBenchmark.Node("1", [], [], [], [ALifeBenchmark.Node("0")])
    rule = ALifeBenchmark.Rule("1", "11", "010", false, false, false, true, false, false, 0, 0)

    ALifeBenchmark.apply_to_all(ALifeBenchmark.fill_temp_links!, node)
    ALifeBenchmark._apply_rule!(rule, node)
    ALifeBenchmark.apply_to_all(ALifeBenchmark.apply_temp_links!, node)

    @test node.string == "11"
    @test length(node.out_excitatory) == 2
    @test node.out_excitatory[2].in_excitatory[1] == node
    @test node.out_excitatory[2].out_excitatory[1] == node.out_excitatory[1]

    node = ALifeBenchmark.Node("1", [], [ALifeBenchmark.Node("101")], [ALifeBenchmark.Node("0")], [])
    rule = ALifeBenchmark.Rule("1", "11", "010", true, false, false, false, true, false, 0, 0)

    ALifeBenchmark.apply_to_all(ALifeBenchmark.fill_temp_links!, node)
    ALifeBenchmark._apply_rule!(rule, node)
    ALifeBenchmark.apply_to_all(ALifeBenchmark.apply_temp_links!, node)

    @test node.string == "11"
    @test length(node.out_excitatory) == 0
    @test length(node.out_inhibitory) == 1
    @test length(node.in_excitatory) == 1
    @test length(node.in_inhibitory) == 1

    new = node.in_inhibitory[1]
    @test new.string == "010"
    @test length(new.out_excitatory) == 0
    @test length(new.out_inhibitory) == 1
    @test length(new.in_excitatory) == 0
    @test length(new.in_inhibitory) == 1
    @test new.out_inhibitory[1] == node
    @test new.in_inhibitory[1] == node.in_excitatory[1]

end

@testset "Complete Development Test 1" begin

    node1 = ALifeBenchmark.get_axiom_node()
    node2 = node1.out_excitatory[1]
    node3 = node2.out_excitatory[1]

    rules = [
        ALifeBenchmark.Rule("0", "011", "1101", false, true, false, false, false, true, 1, 1),
        ALifeBenchmark.Rule("00", "101", "", false, true, true, true, false, false, 1, 1),
        ALifeBenchmark.Rule("011", "", "0110", false, false, true, true, false, false, 1, 1),
        ALifeBenchmark.Rule("1101", "110", "1010", false, true, true, false, false, true, 1, 1)
    ]

    @test node1.string == "001"
    @test node2.string == "000"
    @test node3.string == "01"

    network = ALifeBenchmark.Network(node1)

    ALifeBenchmark.develop_nodes!(network, rules)

    node4 = node2.out_excitatory[2]

    @test node1.string == "101"
    @test node2.string == "101"
    @test node3.string == "011"
    @test node4.string == "1101"

    @test length(node1.in_inhibitory) == 0
    @test length(node1.in_excitatory) == 0
    @test length(node1.out_inhibitory) == 0
    @test length(node1.out_excitatory) == 1

    @test length(node2.in_inhibitory) == 0
    @test length(node2.in_excitatory) == 1
    @test length(node2.out_inhibitory) == 0
    @test length(node2.out_excitatory) == 2

    @test length(node3.in_inhibitory) == 0
    @test length(node3.in_excitatory) == 2
    @test length(node3.out_inhibitory) == 0
    @test length(node3.out_excitatory) == 0

    @test length(node4.in_inhibitory) == 0
    @test length(node4.in_excitatory) == 1
    @test length(node4.out_inhibitory) == 0
    @test length(node4.out_excitatory) == 1

    @test node1.out_excitatory[1] == node2
    @test node2.in_excitatory[1] == node1

    @test node2.out_excitatory[1] == node3
    @test node3.in_excitatory[1] == node2

    @test node2.out_excitatory[2] == node4
    @test node4.in_excitatory[1] == node2

    @test node4.out_excitatory[1] == node3
    @test node3.in_excitatory[2] == node4

    ALifeBenchmark.develop_nodes!(network, rules)

    node3 = node4.out_excitatory[1]
    node5 = node4.out_inhibitory[1]

    @test length(node1.in_inhibitory) == 0
    @test length(node1.in_excitatory) == 0
    @test length(node1.out_inhibitory) == 0
    @test length(node1.out_excitatory) == 1

    @test length(node2.in_inhibitory) == 0
    @test length(node2.in_excitatory) == 1
    @test length(node2.out_inhibitory) == 0
    @test length(node2.out_excitatory) == 3

    @test length(node3.in_inhibitory) == 0
    @test length(node3.in_excitatory) == 2
    @test length(node3.out_inhibitory) == 0
    @test length(node3.out_excitatory) == 0

    @test length(node4.in_inhibitory) == 0
    @test length(node4.in_excitatory) == 2
    @test length(node4.out_inhibitory) == 1
    @test length(node4.out_excitatory) == 1

    @test length(node5.in_inhibitory) == 1
    @test length(node5.in_excitatory) == 1
    @test length(node5.out_inhibitory) == 0
    @test length(node5.out_excitatory) == 1

    @test node1.out_excitatory[1] == node2
    @test node2.in_excitatory[1] == node1

    @test node2.out_excitatory[2] == node3
    @test node3.in_excitatory[1] == node2

    @test node2.out_excitatory[1] == node4
    @test node4.in_excitatory[1] == node2

    @test node4.out_excitatory[1] == node3
    @test node3.in_excitatory[2] == node4

    @test node2.out_excitatory[3] == node5
    @test node5.in_excitatory[1] == node2

    @test node5.out_excitatory[1] == node4
    @test node4.in_excitatory[2] == node5

    @test node4.out_inhibitory[1] == node5
    @test node5.in_inhibitory[1] == node4
end

@testset "Complete Development Test 2" begin

    node1 = ALifeBenchmark.get_axiom_node()
    node2 = node1.out_excitatory[1]
    node3 = node2.out_excitatory[1]

    rules = [
        ALifeBenchmark.Rule("001", "0", "01", false, false, false, false, false, false, 1, 1),
        ALifeBenchmark.Rule("0", "01", "000", false, true, false, false, false, true, 1, 1),
        ALifeBenchmark.Rule("000", "101", "00", false, false, false, false, false, false, 1, 1),
        ALifeBenchmark.Rule("", "10", "000", true, false, true, true, true, true, 1, 1),
        ALifeBenchmark.Rule("101", "100", "", false, true, false, false, true, true, 1, 1)
    ]

    @test node1.string == "001"
    @test node2.string == "000"
    @test node3.string == "01"

    network = ALifeBenchmark.Network(node1)

    ALifeBenchmark.develop_nodes!(network, rules)

    node4 = node2.in_excitatory[2]
    node5 = node1.out_excitatory[2]
    node6 = node2.out_excitatory[2]

    @test node1.string == "0"
    @test node2.string == "101"
    @test node3.string == "01"
    @test node4.string == "01"
    @test node5.string == "00"
    @test node6.string == "000"

    @test length(node1.in_inhibitory) == 0
    @test length(node1.in_excitatory) == 0
    @test length(node1.out_inhibitory) == 0
    @test length(node1.out_excitatory) == 2

    @test length(node2.in_inhibitory) == 0
    @test length(node2.in_excitatory) == 2
    @test length(node2.out_inhibitory) == 0
    @test length(node2.out_excitatory) == 2

    @test length(node3.in_inhibitory) == 0
    @test length(node3.in_excitatory) == 3
    @test length(node3.out_inhibitory) == 0
    @test length(node3.out_excitatory) == 0

    @test length(node4.in_inhibitory) == 0
    @test length(node4.in_excitatory) == 0
    @test length(node4.out_inhibitory) == 0
    @test length(node4.out_excitatory) == 1

    @test length(node5.in_inhibitory) == 0
    @test length(node5.in_excitatory) == 1
    @test length(node5.out_inhibitory) == 0
    @test length(node5.out_excitatory) == 1

    @test length(node6.in_inhibitory) == 0
    @test length(node6.in_excitatory) == 1
    @test length(node6.out_inhibitory) == 0
    @test length(node6.out_excitatory) == 1

    @test node1.out_excitatory[1] == node2
    @test node2.in_excitatory[1] == node1

    @test node1.out_excitatory[2] == node5
    @test node5.in_excitatory[1] == node1

    @test node2.out_excitatory[1] == node3
    @test node3.in_excitatory[1] == node2

    @test node2.out_excitatory[2] == node6
    @test node6.in_excitatory[1] == node2

    @test node4.out_excitatory[1] == node2
    @test node2.in_excitatory[2] == node4

    @test node5.out_excitatory[1] == node3
    @test node3.in_excitatory[2] == node5

    @test node6.out_excitatory[1] == node3
    @test node3.in_excitatory[3] == node6
end

max_stars_repo_name: tochsner/ALifeBenchmark

max_issues_repo_name: tochsner/ALifeBenchmark

max_forks_repo_name: tochsner/ALifeBenchmark


text:
using TrajectoryOptimization
const TO = TrajectoryOptimization
using BenchmarkTools
using TrajOptPlots
using LinearAlgebra
using StaticArrays
using FileIO
using MeshCat
using GeometryTypes
using CoordinateTransformations

# Start visualizer
Rot = MRP{Float64}
model = Dynamics.YakPlane(Rot)
if !isdefined(Main,:vis)
    vis = Visualizer(); open(vis);
    set_mesh!(vis, model)
end

# Generate problem
prob = Problems.YakProblems(UnitQuaternion{Float64,IdentityMap})
U0 = deepcopy(controls(prob))

# Solve
solver = iLQRSolver(prob)

initial_controls!(solver, U0)
# benchmark_solve!(solver)
rollout!(solver)
TO.state_diff_size(solver.model)

solver.opts.verbose = true
cost(solver)
solve!(solver)
iterations(solver)
cost(solver)
visualize!(vis, solver)
X = states(solver)
findfirst(isnan.(X))
X[15]

solver.obj.J


text:
include("saha.jl")
include("constants.jl")

function bremsstrahlung(N, T, ν)

    A = 0

    gaunt = 1.0

    #gaunt= (sqrt(3) / π)* log((4.0 / ζ)*((keV * T) / (heV* ν)))

    Ne , HI, HII, HeI, HeII, HeIII = saha(T, N)

    if T <= 2e5
        A = 18.2 + 3/2 * log(T) - log(ν)
    else
        A = 24.5 + log(T) - log(ν)
    end

    κ = 9.78e-3 * (Ne^2 / (ν^2 * T^(3/2))) * A

    #κ = 1.3705e-35 * T^(0.5) * Ne * Ne * ν^(-3) * (1 - exp(-heV * ν / (keV * T) )) * gaunt
    #κ = 3.7e8 * T^(0.5) * Ne * HeIII * ν^(-3) * (1 - exp(-heV * ν / (keV * T) )) * gaunt

    return κ
end

export bremsstrahlung


text:
struct GeneratorWithElType{T, S}
    source::S
end

function IteratorInterfaceExtensions.getiterator(source::Base.Generator)
    TS = eltype(source.iter)
    T = Base.return_types(source.f, (TS,))[1]
    S = typeof(source)

    if isconcretetype(T)
        return GeneratorWithElType{T,S}(source)
    else
        return source
    end
end

function TableTraits.isiterabletable(source::Base.Generator)
    TS = eltype(source.iter)
    T = Base.return_types(source.f, (TS,))[1]

    if isconcretetype(T)
        if T <: NamedTuple
            return true
        else
            return false
        end
    else
        return missing
    end
end

Base.IteratorSize(::Type{GeneratorWithElType{T,S}}) where {T,S} = Base.IteratorSize(S)

Base.size(it::GeneratorWithElType, dims...) = size(it.source, dims...)

Base.length(it::GeneratorWithElType) = length(it.source)

Base.eltype(::Type{GeneratorWithElType{T,TS}}) where {T,TS} = T

Base.iterate(it::GeneratorWithElType) = iterate(it.source)

Base.iterate(it::GeneratorWithElType, state) = iterate(it.source, state)


text:
module TensorGames

using PATHSolver
using SparseArrays
using ChainRulesCore: ChainRulesCore, NoTangent, ZeroTangent, @non_differentiable
using LinearAlgebra: qr, diag, I
using ForwardDiff: ForwardDiff

include("nash.jl")

export compute_equilibrium, expected_cost

end # module


text:
# Index manipulations
"""
    permute(t::AbstractTensorMap{S}, p1::NTuple{N₁}, p2::NTuple{N₂} = ();
                copy::Bool = false)

Permute the indices of tensor map `t` to a new tensor with indices in `p1` playing the role
of the codomain of the map, and indices in `p2` indicating the domain.
If `copy == false`, keep the original tensor if possible.
"""
function permute(t::TensorMap{S}, p1::IndexTuple,  p2::IndexTuple=();
                    copy::Bool = false) where {S}
    cod = ProductSpace{S}(map(n->space(t, n), p1))
    dom = ProductSpace{S}(map(n->dual(space(t, n)), p2))
    # share data if possible
    if !copy
        if p1 === codomainind(t) && p2 === domainind(t)
            return t
        elseif has_shared_permute(t, p1, p2)
            return TensorMap(reshape(t.data, dim(cod), dim(dom)), cod, dom)
        end
    end
    # general case
    @inbounds begin
        return add!(true, t, false, similar(t, cod←dom), p1, p2)
    end
end

"""
    permute(t::AdjointTensorMap{S}, p1::IndexTuple, p2::IndexTuple=();
                        copy::Bool = false) where {S}

Permute the indices of tensor map `t` to a new tensor with indices in `p1` playing the role
of the codomain of the map, and indices in `p2` indicating the domain.
If `copy == false`, keep the original tensor if possible.
"""
function permute(t::AdjointTensorMap{S}, p1::IndexTuple, p2::IndexTuple=();
                    copy::Bool = false) where {S}
    p1′ = map(n->adjointtensorindex(t, n), p2)
    p2′ = map(n->adjointtensorindex(t, n), p1)
    adjoint(permute(adjoint(t), p1′, p2′; copy = copy))
end

"""
    has_shared_permute(t::TensorMap, p1, p2)

Return true if `p1 == (1,...,length(codomain))` and
`p2 == length(codomain)+(1,...,length(domain))`;
Or in the trivial case, the matrix with `p1` and `p2` can be fused into a vector with
continuous memory with stride = 1.
"""
function has_shared_permute(t::TensorMap, p1, p2)
    if p1 === codomainind(t) && p2 === domainind(t)
        return true
    elseif sectortype(t) === Trivial
        stridet = i->stride(t[], i)
        sizet = i->size(t[], i)
        canfuse1, d1, s1 = TO._canfuse(sizet.(p1), stridet.(p1))
        canfuse2, d2, s2 = TO._canfuse(sizet.(p2), stridet.(p2))
        return canfuse1 && canfuse2 && s1 == 1 && (d2 == 1 || s2 == d1)
    else
        return false
    end
end

function has_shared_permute(t::AdjointTensorMap, p1, p2)
    p1′ = adjointtensorindices(t, p2)
    p2′ = adjointtensorindices(t, p1)
    return has_shared_permute(t', p1′, p2′)
end

"""
    Base.permute!(tdst::AbstractTensorMap{S, N₁, N₂}, tsrc::AbstractTensorMap{S},
                    p1::IndexTuple{N₁}, p2::IndexTuple{N₂}=()) where {S, N₁, N₂}

Permute the indices of `tsrc` and write the result into `tdst`, with indices in `p1` playing
the role of the codomain or range of the map `tdst`, and indices in `p2` indicating the
domain of `tdst`.

For more details see [`add!`](@ref), of which this is a special case.
"""
@propagate_inbounds Base.permute!(tdst::AbstractTensorMap{S, N₁, N₂},
                                    tsrc::AbstractTensorMap{S},
                                    p1::IndexTuple{N₁},
                                    p2::IndexTuple{N₂}=()) where {S, N₁, N₂} =
    add_permute!(true, tsrc, false, tdst, p1, p2)

"""
    braid(t::TensorMap{S}, levels::IndexTuple, p1::IndexTuple, p2::IndexTuple=();
            copy::Bool = false) where {S}

Braid the indices of tensor map `t` to a new tensor with indices in `p1` playing the role
of the codomain of the map, and indices in `p2` indicating the domain, according to the
`levels`. If `copy == false`, keep the original tensor if possible.
"""
function braid(t::TensorMap{S}, levels::IndexTuple, p1::IndexTuple, p2::IndexTuple=();
                copy::Bool = false) where {S}
    @assert length(levels) == numind(t)
    if BraidingStyle(sectortype(S)) isa SymmetricBraiding
        return permute(t, p1, p2; copy = copy)
    end
    if !copy && p1 == codomainind(t) && p2 == domainind(t)
        return t
    end
    # general case
    cod = ProductSpace{S}(map(n->space(t, n), p1))
    dom = ProductSpace{S}(map(n->dual(space(t, n)), p2))
    @inbounds begin
        return add_braid!(true, t, false, similar(t, cod←dom), p1, p2, levels)
    end
end

"""
    braid!(tdst::AbstractTensorMap{S, N₁, N₂}, tsrc::AbstractTensorMap{S},
            levels::IndexTuple, p1::IndexTuple{N₁}, p2::IndexTuple{N₂}=()) where {S, N₁, N₂}

Braid the indices of tensor map `tsrc` and write into `tdst` with indices in `p1` playing
the role of the codomain of the map, and indices in `p2` indicating the domain, according to
the `levels`.
"""
@propagate_inbounds braid!(tdst::AbstractTensorMap{S, N₁, N₂}, tsrc::AbstractTensorMap{S},
                            levels::IndexTuple, p1::IndexTuple{N₁},
                            p2::IndexTuple{N₂}=()) where {S, N₁, N₂} =
    add!(true, tsrc, false, tdst, p1, p2, levels)

"""
    LinearAlgebra.transpose(t::TensorMap{S}, p1::IndexTuple = reverse(domainind(t)),
                            p2::IndexTuple = reverse(codomainind(t));
                            copy::Bool = false) where {S}

Transpose the indices of tensor map `t` to a new tensor with indices in `p1` playing the
role of the codomain of the map, and indices in `p2` indicating the domain.
If `copy == false`, keep the original tensor if possible.
"""
function LinearAlgebra.transpose(t::TensorMap{S}, p1::IndexTuple = reverse(domainind(t)),
                                    p2::IndexTuple = reverse(codomainind(t));
                                    copy::Bool = false) where {S}
    if sectortype(S) === Trivial
        return permute(t, p1, p2; copy = copy)
    end
    if !copy && p1 == codomainind(t) && p2 == domainind(t)
        return t
    end
    # general case
    cod = ProductSpace{S}(map(n->space(t, n), p1))
    dom = ProductSpace{S}(map(n->dual(space(t, n)), p2))
    @inbounds begin
        return add_transpose!(true, t, false, similar(t, cod←dom), p1, p2)
    end
end

"""
    LinearAlgebra.transpose(t::AdjointTensorMap{S},
                            p1::IndexTuple = reverse(domainind(t)),
                            p2::IndexTuple = reverse(codomainind(t));
                            copy::Bool = false) where {S}

Transpose the indices of tensor map `t` to a new tensor with indices in `p1` playing the
role of the codomain of the map, and indices in `p2` indicating the domain.
If `copy == false`, keep the original tensor if possible.
"""
function LinearAlgebra.transpose(t::AdjointTensorMap{S},
                                    p1::IndexTuple = reverse(domainind(t)),
                                    p2::IndexTuple = reverse(codomainind(t));
                                    copy::Bool = false) where {S}
    p1′ = map(n->adjointtensorindex(t, n), p2)
    p2′ = map(n->adjointtensorindex(t, n), p1)
    adjoint(transpose(adjoint(t), p1′, p2′; copy = copy))
end

"""
    LinearAlgebra.transpose!(tdst::AbstractTensorMap, tsrc::AbstractTensorMap,
                                p1::IndexTuple, p2::IndexTuple)

Transpose the indices of tensor map `tsrc` and write into `tdst` with indices in `p1`
playing the role of the codomain of the map, and indices in `p2` indicating the domain.
"""
LinearAlgebra.transpose!(tdst::AbstractTensorMap, tsrc::AbstractTensorMap,
                            p1::IndexTuple, p2::IndexTuple) =
    add_transpose!(true, tsrc, false, tdst, p1, p2)

"""
    twist(t::AbstractTensorMap, i::Int; inv::Bool = false)

Twist the `i`th index of the tensor map `t` and returned as a new tensor.
Inverse the twist if `inv == true`.
"""
twist(t::AbstractTensorMap, i::Int; inv::Bool = false) = twist!(copy(t), i; inv = inv)

"""
    twist!(t::AbstractTensorMap, i::Int; inv::Bool = false)

Twist the `i`th index of the tensor map `t` and replace it.
Inverse the twist if `inv == true`.
"""
function twist!(t::AbstractTensorMap, i::Int; inv::Bool = false)
    if i > numind(t)
        msg = "Can't twist index $i of a tensor with only $(numind(t)) indices."
        throw(ArgumentError(msg))
    end
    BraidingStyle(sectortype(t)) == Bosonic() && return t
    N₁ = numout(t)
    for (f1, f2) in fusiontrees(t)
        θ = i <= N₁ ? twist(f1.uncoupled[i]) : twist(f2.uncoupled[i-N₁])
        inv && (θ = θ')
        rmul!(t[f1, f2], θ)
    end
    return t
end


text:
"""
    StabilizerReg{Tzv, Tm} <: AbstractRegister{1}

A data struct for stabilizer states in Yao. To initialize a `StabilizerReg`,
one can use [`zero_stabilizer_state`](@ref) and [`rand_stabilizer_state`](@ref).
"""
mutable struct StabilizerReg{Tzv, Tm} <: AbstractRegister{1}
    st::Stabilizer{Tzv, Tm}
end

Base.copy(reg::StabilizerReg{Tzv, Tm}) where {Tzv, Tm} = StabilizerReg{Tzv, Tm}(copy(reg.st))
function Base.show(io::IO, reg::StabilizerReg)
    println(io, typeof(reg), " with $(Yao.nqubits(reg)) qubits:")
    show(io, reg.st)
end

"""
    zero_stabilizer_state(n::Integer) -> StabilizerReg

Initialize a stabilizer state |0...0> of n-qubits.
"""
zero_stabilizer_state(n::Integer) = StabilizerReg(one(Stabilizer, n))

"""
    rand_stabilizer_state(n::Integer) -> StabilizerReg

Initialize a random stabilizer state of n-qubits.
"""
rand_stabilizer_state(n::Integer) = StabilizerReg(random_stabilizer(n))

function Yao.addbits!(reg::StabilizerReg, n::Integer)
    reg.st = reg.st ⊗ one(Stabilizer, n)
    return reg
end
Yao.fidelity(r1::StabilizerReg, r2::StabilizerReg) = LinearAlgebra.dot(r1.st, r2.st)
Yao.instruct!(reg::StabilizerReg, op...) = throw(MethodError(Yao.instruct!, (reg, op...)))
Yao.nactive(reg::StabilizerReg) = Yao.nqubits(reg)
Yao.nbatch(reg::StabilizerReg{Tzv, Tm}) where {Tzv, Tm} = 1
Yao.nqubits(reg::StabilizerReg) = QuantumClifford.nqubits(reg.st)
Yao.nremain(reg::StabilizerReg) = 0
function Yao.partial_tr(reg::StabilizerReg, locs)
    n = Yao.nqubits(reg)
    qubit_map = sort!(setdiff(1:n, locs))
    new_st = zero(Stabilizer, length(qubit_map))
    rank = 0
    temp_st = traceout!(copy(reg.st), locs)
    for r = 1:length(temp_st)
        rank += 1
        for i = 1:length(qubit_map)
            if !iszero(temp_st[r])
                new_st[rank, i] = temp_st[r, qubit_map[i]]
            end
        end
    end
    return StabilizerReg(new_st)
end
function Yao.reorder!(reg::StabilizerReg, orders)
    reg2 = copy(reg)
    for i = 1:length(orders), r = 1:length(reg.st)
        reg.st[r, orders[i]] = reg2.st[r, i]
    end
    return reg
end
Yao.invorder!(reg::StabilizerReg) = Yao.reorder!(reg, Yao.nqubits(reg):-1:1)
Yao.invorder(reg::StabilizerReg) = Yao.invorder!(copy(reg))

PauliString{N} = KronBlock{N, N, <:NTuple{N,Yao.YaoBlocks.PauliGate}}
"""
    PauliString(paulis...)

Generate a `KronBlock` to represent a Pauli operator.
"""
PauliString(args::YaoBlocks.PauliGate...) = kron(args...)

# focus!
# insert_qubits!
# collapseto!
# density_matrix
# probs
# purify
# relax!
# state
# statevec
# collapseto!
# density_matrix
# von_neumann_entropy
# ρ
# basis


text:
using MinimaxCalibratedEBayes
using Test
using Plots
pgfplotsx()

xs = range(-3,3;step=0.5)
ws = Float64.(collect(1:14))
tmp_hist = MCEBHistogram(xs)

@test StatsBase.midpoints(tmp_hist.grid) == StatsBase.midpoints(tmp_hist)

calib_normal = DiscretizedAffineEstimator(tmp_hist, Normal(0, 0.5))
calib_sin = DiscretizedAffineEstimator(tmp_hist, sin)

calib_sinc = DiscretizedAffineEstimator(tmp_hist, SincKernel(m=10_000))

calib_array = [calib_sin, calib_normal, calib_sinc]
@test isa(calib_array, AbstractVector{<:DiscretizedAffineEstimator})

plot(calib_normal)
plot(calib_sin)

plot(calib_array, label=["bla" "blabla" "sinc"])

z = DiscretizedStandardNormalSample(0.5, tmp_hist)
zs = DiscretizedStandardNormalSample.(randn(100), tmp_hist)

@test var(marginalize(Normal(0,0), StandardNormalSample(1.0))) == 1.0
@test var(marginalize(Normal(0,1), StandardNormalSample(1.0))) ≈ 2.0

tmp_marginalize = marginalize(Normal(), z)
marginalize_probs = pdf(tmp_marginalize)
@test sum(pdf(tmp_marginalize)) == 1

midpoints(marginalize_probs)

plot(marginalize_probs)

first(tmp_hist)

zs[1].mhist
tmp_calib = DiscretizedAffineEstimator(tmp_hist, ws, 1.0)

@test tmp_calib(-4) == 2.00
@test tmp_calib.([-4, 4]) == [2.0, 15.0]
@test tmp_calib.([-4, -3, -2.9, 4]) == [2.0, 2.0, 3.0, 15.0]

@test step(tmp_hist) == step(xs)
@test first(tmp_hist) == first(xs)
@test last(tmp_hist) == last(xs)


text:
export PDEfunc!, pdefunc_ops
pdefunc_ops = [:objective,:field,:forw_model,:jacob_forw,:jacob_adj,:hess,:hess_gn]

function PDEfunc!(op::Symbol,
                  v::Union{AbstractArray{F,1},AbstractArray{Complex{F},1}},
                  Q::Union{AbstractArray{F,1},AbstractArray{Complex{F},1}},
                  Dobs::Union{AbstractArray{F,1},AbstractArray{Complex{F},1}},
                  input::Union{AbstractArray{F,1},AbstractArray{Complex{F},1}},
                  model::Model{I,F},
                  opts::PDEopts{I,F};
                  grad::Union{AbstractArray{F,1},AbstractArray{Complex{F},1}}=Array{F}(0),
                  srcfreqmask::AbstractArray{Bool,2}=Array{Bool}(0,0)) where {F<:AbstractFloat,I<:Integer}

    op in pdefunc_ops || error("Unrecognized op $op, must be one of $pdefunc_ops")
    nsrc = length(model.xsrc)*length(model.ysrc)*length(model.zsrc)
    nfreq = length(model.freq)
    nrec = length(model.xrec)*length(model.yrec)*length(model.zrec)
    nmodel = prod(model.n)
    if isempty(srcfreqmask)
        srcfreqmask = trues(nsrc,nfreq)
    end
    size(srcfreqmask,1)==nsrc && size(srcfreqmask,2)==nfreq || error("srcfreqmask must be a nsrc x nfreq matrix")
    Iactive = find(srcfreqmask)
    npde_out = length(Iactive)
    (iS,iF) = ind2sub((nsrc,nfreq),Iactive)
    freqsxsy = Array{Array{I,1},1}(nfreq)

    numcompsrc = (length(model.n)==2||model.n[3]==1) ? nsrc : 1
    length(Dobs)==0 || length(Dobs)==nrec*npde_out || error("Dobs must be a nrec*npde_out length vector")

    if !isempty(Dobs)
        Dobs = reshape(Dobs,nrec,npde_out)
    end

    if length(Q)==nsrc*nsrc
        Q = reshape(Q,(nsrc,nsrc))
    elseif length(Q) ==nsrc*nsrc*nfreq
        Q = reshape(Q,(nsrc,nsrc,nfreq))
    else
        throw(Exception("Q must be nsrc x nsrc or nsrc x nsrc x nfreq"))
    end
    for i in 1:nfreq
        J = find(iF.==i)
        if !isempty(J)
            freqsxsy[i] = iS[J]
        else
            freqsxsy[i] = Array{I,1}()
        end
    end

    # No work to do, just return the right outputs
    if npde_out==0
        z = zeros(eltype(v),(nmodel))
        if op==:objective
            return 0.0
        elseif op in [:forw_model,:jacob_forw]
            return (zeros{Complex{F}}(nrec,0))
        else
            return (z)
        end
    end

    # Set up outputs
    if op==:objective
        f = 0.0
        compute_grad = length(grad) > 0
        if compute_grad
            length(grad)==length(v) || error("length(g) $(length(g)) != length(v) $(length(v))")
            fill!(grad,0.0)
        end
        misfit = opts.misfit
    elseif op in [:forw_model,:jacob_forw]
        output = zeros(Complex{F},(nrec,npde_out))
    else
        output = zeros(F,(nmodel))
        if op==:jacob_adj
            input = reshape(input,nrec,npde_out)
        end
    end
    npdes = 1
    freqs = model.freq
    w = fwi_wavelet(freqs,model.t0,model.f0)

    freq_idx = 0
    for k in 1:nfreq

        if isempty(freqsxsy[k])
            continue
        end
        freq_idx = freq_idx+1
        freq = freqs[k]
        isrc = freqsxsy[k]

        src_batches = index_block(isrc,numcompsrc)

        (H,comp_grid,T,DT_adj) = helmholtz_system(v,model,freq,opts)
        phys_to_comp = comp_grid.phys_to_comp_grid
        comp_to_phys = comp_grid.comp_to_phys_grid
        (Ps,Pr) = src_rec_interp_operators(model,comp_grid)
        if op in [:jacob_forw :hess_gn :hess]
            δm = phys_to_comp*input;
        end

        for j in 1:length(src_batches)

            current_src_idx = src_batches[j]
            # Scaling so that the wavefield amplitudes are the same for different grid spacings
            if ndims(Q)==3
                q = Q[:,current_src_idx,k]
            else
                q = Q[:,current_src_idx]
            end

            q = jo_convert(Complex{F},Ps*(w[k]*q))
            q = q*prod(model.d)/prod(comp_grid.comp_d)
            data_idx = npdes:npdes+length(current_src_idx)-1

            # Wavefield solve
            U = H\q
            sum_srcs = x->real(comp_to_phys*sum(x,2))
            if op==:objective
                #figure()
                #imshow(real(Pr*U),aspect="auto")
                #title("Pred data at $(model.freq[k]) Hz")
                (ϕ,δϕ) = misfit(Pr*U,Dobs[:,data_idx])
                f += ϕ
                if compute_grad
                    V = H'\(-(Pr'*δϕ))
                    grad .+= squeeze(sum_srcs(T(U)'*V),2)
                end
            elseif op==:field
                output[:,data_idx] = U
            elseif op==:forw_model
                output[:,data_idx] = Pr*U
            elseif op==:jacob_forw
                δU = H\(-T(U)*δm)
                output[:,data_idx] = Pr*δU
            elseif op==:jacob_adj
                V = H'\(-(Pr'*input[:,data_idx]))
                output .+= sum_srcs(T(U)'*V)
            elseif op==:hess_gn
                δU = H\(-T(U)*δm)
                δU .= H'\(-(Pr'*(Pr*δU)))
                output .+= sum_srcs(T(U)*δU)
            elseif op==:hess
                (ϕ,δϕ,δ2ϕ) = misfit(Pr*U,Dobs[:,data_idx])
                δU = H\(-T(U)*δm)
                V = H'\(-Pr'*δϕ)
                δV = H'\(-T(V)*δm - Pr'*(δ2ϕ.* (Pr*δU) ) )
                output .+= sum_srcs(DT_adj(U,δm,δU)*V+T(U)*δV)
            end
            npdes += length(current_src_idx)
        end
    end
if(op==:objective)
    return f
else
    return output
end
end

function src_rec_interp_operators(model::Model{I,F},comp_grid::ComputationalGrid{I,F}) where {F<:AbstractFloat,I<:Integer}
    ndims = (length(model.n)==2||model.n[3]==1) ? 2 : 3
    kaiser_window_param = 4
    DT = Complex{F}
    RT = Complex{F}
    opS(x,y) = joSincInterp(x,y,DomainT=DT,RangeT=RT)
    if ndims==2
        (zt,xt) = odn_to_grid(comp_grid)
        Ps = joKron(opS(model.xsrc,xt),opS(model.zsrc,zt))
        Pr = joKron(opS(xt,model.xrec),opS(zt,model.zrec))
    else
        (xt,yt,zt) = odn_to_grid(comp_grid)
        Ps = joKron(opS(model.zsrc,zt),opS(model.ysrc,yt),opS(model.xsrc,xt))
        Pr = joKron(opS(zt,model.zrec),opS(yt,model.yrec),opS(xt,model.xrec))
    end
    return (Ps,Pr)
end


text:
@testset "geometry.coordsinternal" begin

@testset "internalcoords" begin
    emptycoords = InternalCoords(10)
    @test length(emptycoords) == 10
    coords = internalcoords([
        nothing nothing nothing;
        1 nothing nothing;
        2 1 nothing;
        3 2 1
    ], [
        nothing nothing nothing;
        1.0 nothing nothing;
        1.0 2/3 nothing;
        1.0 2/3 0.0
    ])
    @test length(coords) == 4
    @test label1(coords, 2) == 1
    @test label2(coords, 3) == 1
    @test label3(coords, 4) == 1
    @test distance(coords, 2) == 1.0
    @test angle(coords, 3) == 2 / 3
    @test dihedral(coords, 4) == 0.0
end

end # geometry.coordsinternal


text:
include("../aimajulia.jl");
using aimajulia;
using aimajulia.utils;

function beautify_node_args(n)
    if (typeof(n) <: Node)
        return beautify_node(n);
    elseif (typeof(n) <: Tuple)
        return map(beautify_node_args, n);
    else
        return n;
    end
end

function bfgs{T <: AbstractProblem}(problem::T, f::Function)
    mf = MemoizedFunction(f);
    local node = Node{typeof(problem.initial)}(problem.initial);
    if (goal_test(problem, node.state))
        return node;
    end
    local frontier = PQueue();
    push!(frontier, node, mf);
    local explored = Set{String}();
    while (length(frontier) != 0)
        node = pop!(frontier);
        if (goal_test(problem, node.state))
            return node;
        end
        push!(explored, node.state);
        for child_node in expand(node, problem)
            if (!(child_node.state in explored) &&
                !(child_node in collect(getindex(x, 2) for x in frontier.array)))
                push!(frontier, child_node, mf);
            elseif (child_node in [getindex(x, 2) for x in frontier.array])
            #Recall that Nodes can share the same state and different values for other fields.
                local existing_node = pop!(collect(getindex(x, 2)
                                                    for x in frontier.array
                                                    if (getindex(x, 2) == child_node)));

                eval_memoized_function(mf, child_node);
                eval_memoized_function(mf, existing_node);

                if (eval_memoized_function(mf, child_node) < eval_memoized_function(mf, existing_node))
                    delete!(frontier, existing_node);
                    push!(frontier, child_node, mf);
                end
            end
        end
        print("length of memoization dictionary: ", length(mf.values), " ");
        println(map(beautify_node_args, collect(keys(mf.values)))...);
    end
    return nothing;
end

function astar_memoized_search(problem::GraphProblem; h::Union{Void, Function}=nothing)
    local mh::MemoizedFunction; #memoized h(n) function
    if (!(typeof(h) <: Void))
        mh = h;
    else
        mh = problem.h;
    end
    return bfgs(problem,
                                    (function(node::Node; h::MemoizedFunction=mh, prob::GraphProblem=problem)
                                        return node.path_cost + eval_memoized_function(h, prob, node);end));
end

astar_str="";for i in 1:5
	if (Node{String}("P")==get(astar_memoized_search(GraphProblem("A", "B", aimajulia.romania)).parent))
		astar_str = astar_str * "1";
        println("Test: Passed!");
	else
		astar_str = astar_str * "0";
        println("Test: Failed!");
	end
end;println(count(i->(i=='1'), astar_str), " of 5 tries passed!");

max_stars_repo_name: mikhail-j/aima-julia

max_issues_repo_name: mikhail-j/aima-julia

max_forks_repo_name: mikhail-j/aima-julia


text:
"""
    struct ComputeShaderPipeline{PushConstsT, NBuffers}

An overarching structure that holds everything needed to run a single
compute-only pipeline: the shader, descriptor set and pipeline layouts, the
pipeline, a single small descriptor pool, and the descriptor sets.
"""
struct ComputeShaderPipeline{PushConstsT,NBuffers}
    shader::ShaderModule
    descriptor_set_layout::DescriptorSetLayout
    pipeline_layout::PipelineLayout
    pipeline::Pipeline
    descriptor_pool::DescriptorPool
    descriptor_sets::Vector{DescriptorSet}
end


text:
module SoftRobots

import DataStructures: OrderedDict
import Iterators: product
import Base: convert
import MultiPoly: MPoly

using SpatialFields
using Meshing
using GeometryTypes
using PyLCM
using PyCall
global spatial
global lcmdrake

function __init__()
    const global spatial = pyimport("scipy.spatial")
    const global lcmdrake = pyimport("drake")
end

abstract Object
abstract StaticObject <: Object
abstract DynamicObject <: Object

type DampedSpring
    parent::Int
    child::Int
    stiffness::Float64
    damping::Float64
    rest_length::Float64
end

type PointMass
    mass::Float64
end

type SoftRobot{T, IndexOffset} <: DynamicObject
    nodes::Vector{PointMass}
    edges::Vector{DampedSpring}
    faces::Vector{Face{3, T, IndexOffset}}
end

function SoftRobot{T, IndexOffset}(nodes::Vector{PointMass}, edges::Vector{DampedSpring}, faces::Vector{Face{3, T, IndexOffset}})
  SoftRobot{T, IndexOffset}(nodes, edges, faces)
end

immutable FixedObject <: StaticObject
end

type World
    objects::Array{Object}
    gravity::Vector
    friction::Number

    World(gravity::Vector) = new([], gravity, 1.0)
end

function World3D()
    World([0, 0, -9.81])
end

abstract ObjectState
abstract DynamicObjectState <: ObjectState
abstract StaticObjectState <: ObjectState

typealias WorldState Dict{Object, ObjectState}

type SoftRobotState{T} <: DynamicObjectState
    positions::Vector{Point{3, T}}
    velocities::Vector{Point{3, T}}
    barrier::ScalarField
    velocity_field::VectorField

    SoftRobotState(positions, velocities) = new(positions, velocities)
end

function SoftRobotState{T}(positions::Vector{Vector{T}}, velocities::Vector{Vector{T}})
    SoftRobotState{T}(map(x -> convert(Point{3, T}, x), positions), map(x -> convert(Point{3, T}, x), velocities))
end

function SoftRobotState{T}(positions::Vector{Point{3, T}}, velocities::Vector{Point{3, T}})
    SoftRobotState{T}(map(x -> convert(Point{3, T}, x), positions), map(x -> convert(Point{3, T}, x), velocities))
end

type FixedObjectState <: StaticObjectState
    barrier::ScalarField
    velocity_field::VectorField

    FixedObjectState(barrier) = new(barrier, [MPoly{Float64}(OrderedDict(zeros(Int64, length(barrier.vars)) => 0.0), barrier.vars) for i in 1:length(barrier.vars)])
end

function update!(world::World, states::Dict{Object, ObjectState}, dt::Number)
    for (object, state) in states
        update_barrier!(object, state)
        update_velocity_field!(object, state)
    end
    for object in keys(states)
        collisions!(object, states, world.friction, dt)
        if isa(object, SoftRobot)
            if any(x -> any(isnan, x), states[object].positions)
                @show states[object]
                error("nans in positions after collision")
            end
        end
    end
    for (object, state) in states
        dynamics!(object, state, world.gravity, dt)
    end
end

function update_barrier!(object::FixedObject, state::FixedObjectState)
    # Nothing to do here
end

function update_velocity_field!(object::FixedObject, args...)
    # Nothing to do here
end

function collisions!(object::FixedObject, args...)
    # Nothing to do here
end

function dynamics!(object::FixedObject, args...)
    # Nothing to do here
end

function normalize!(x)
    n = norm(x)
    @inbounds @simd for i = 1:length(x)
        x[i] /= n
    end
end

function update_barrier!{T}(robot::SoftRobot, state::SoftRobotState{T})
    node_is_on_face = falses(length(robot.nodes))
    for face in robot.faces
        for i in 1:length(face)
            node_is_on_face[onebased(face, i)] = true
        end
    end
    mesh = HomogenousMesh(state.positions, robot.faces)
    normals = decompose(GeometryTypes.Normal{3, T}, mesh)
    if any(x -> any(isnan, x), normals[node_is_on_face])
        @show robot.faces
        @show state.positions
        @show state.positions[node_is_on_face]
        @show normals[node_is_on_face]
    end
    # state.barrier = HermiteRadialField(state.positions[node_is_on_face], normals[node_is_on_face])
    values = zeros(length(state.positions))
    values[!node_is_on_face] = -1.0
    state.barrier = InterpolatingSurface(state.positions, values, SpatialFields.XSquaredLogX())
    # state.barrier = InterpolatingSurface(state.positions, values)
end

function update_velocity_field!{T}(robot::SoftRobot, state::SoftRobotState{T})
    # mean_velocity = mean(state.velocities)
    # dimension = length(state.positions[1])
    # state.velocity_field = PolynomialVectorField([MPoly{T}(OrderedDict(zeros(dimension) => v), [:x, :y, :z][1:dimension]) for v in mean_velocity])
    state.velocity_field = PolynomialVectorField(SpatialFields.linear_fit(state.positions, state.velocities))
end

function collisions!(robot::SoftRobot, states::Dict{Object, ObjectState}, friction::Number, dt::Number)
    my_state = states[robot]
    for (obj, other_state) in states
        if obj == robot
            # Don't collide with yourself
            continue
        end
        barrier_collisions!(robot, my_state, other_state, friction, dt)
    end
end

function dynamics!(robot::SoftRobot, state::SoftRobotState, gravity::Vector, dt::Number)
    dim = length(state.positions[1])
    total_forces = [robot.nodes[i].mass * gravity for i = 1:length(robot.nodes)]
    update_internal_forces!(robot, state, total_forces)
    euler_integrate!(robot, state, total_forces, dt)
end

function rotmat(theta)
    [cos(theta) -sin(theta);
     sin(theta) cos(theta)]
end

function barrier_collisions!{T}(robot::SoftRobot, my_state::SoftRobotState{T}, other_state::ObjectState, friction::Number, dt::Number)
    barrier_graident = grad(other_state.barrier)
    for i = 1:length(robot.nodes)
        barrier_velocity::Point{3, T} = evaluate(other_state.velocity_field, my_state.positions[i])
        relative_velocity::Point{3, T} = my_state.velocities[i] - barrier_velocity
        predicted_relative_position::Point{3, T} = my_state.positions[i] + relative_velocity * dt
        if evaluate(other_state.barrier, predicted_relative_position) <= 0
            gradient::Point{3, T} = evaluate(barrier_graident, predicted_relative_position)
            normal = gradient ./ norm(gradient)
            if any(isnan, normal)
                @show gradient
                @show normal
                error("nan in normal")
            end
            penetration_velocity = -dot(relative_velocity, normal)
            if penetration_velocity > 0
                # my_state.velocities[i] += normal * penetration_velocity
                my_state.velocities[i] = barrier_velocity
            end
            # gradient = rotmat(friction * pi*(rand(1) - 0.5)) * gradient
            # normal = gradient / norm(gradient)
            # penetration_velocity = -dot(relative_velocity, normal)
            # if penetration_velocity > 0
            #     my_state.velocities[i] += normal * penetration_velocity
            # end
        end
    end
end

function update_internal_forces!{T}(robot::SoftRobot, state::SoftRobotState{T}, total_forces)
    dim = length(total_forces[1])
    displacement = Array{T}(dim)
    direction = Array{T}(dim)
    relative_velocity = Array{T}(dim)
    force = Array{T}(dim)
    distance = 0.
    edge_velocity = 0.
    force = 0.
    for edge in robot.edges
        for j = 1:dim
            displacement[j] = state.positions[edge.child][j] - state.positions[edge.parent][j]
        end
        distance = norm(displacement)
        for j = 1:dim
            direction[j] = displacement[j] ./ distance
            relative_velocity[j] = state.velocities[edge.child][j] - state.velocities[edge.parent][j]
        end
        edge_velocity = dot(relative_velocity, direction)
        for j = 1:dim
            force = direction[j] * (edge.stiffness * (edge.rest_length - distance) - edge.damping * edge_velocity)
            total_forces[edge.child][j] += force
            total_forces[edge.parent][j] -= force
        end
    end
end

function euler_integrate!{T}(robot::SoftRobot, state::SoftRobotState{T}, total_forces, dt::Real)
    dim = length(state.positions[1])
    for i = 1:length(robot.nodes)
        delta_v = Point{3, T}(total_forces[i]) / robot.nodes[i].mass * dt
        # state.velocities[i] += half_delta_v
        # state.velocities[i] += Point(total_forces[i]) / robot.nodes[i].mass * dt
        state.positions[i] += state.velocities[i] * dt
        state.velocities[i] += delta_v
    end
end

function convex_hull{T}(nodes::Vector{Point{3, T}})
    hull = spatial[:ConvexHull](hcat(map(x -> convert(Vector, x), nodes)...)')
    simplices = hull[:simplices]
    simplices += 1

    faces = Face{3, Int, 0}[]
    # Reorient simplices so that normals always point out from the hull
    for i = 1:size(simplices, 1)
        v = Vector{T}[nodes[simplices[i, j]] for j in 1:size(simplices, 2)]

        # Given a face of a convex hull, all of the points in the body must be on the same side of that face. So in theory, we just need to pick one point not on the face and check which side it's on. But to attempt to be robust to numerical issues, we'll actually sum the dot product with the normal for every point in the body, and check whether that sum is positive or negative. If this becomes a performance bottleneck, we can revisit it later.
        n = Point{3, T}(cross(v[2] - v[1], v[3] - v[1]))
        b = dot(n, nodes[simplices[i, 1]])
        total = zero(T)
        for j = 1:length(nodes)
            total += dot(n, nodes[j]) - b
        end
        if total < 0
            # Then the normal is pointing the right way
            push!(faces, Face{3, Int, 0}((simplices[i,:])...))
        else
            # Otherwise the normal is pointing the wrong way, so we flip the face
            push!(faces, Face{3, Int, 0}((simplices[i,end:-1:1])...))
        end
    end
    return faces
end

function barrier_mesh(bounds::HyperRectangle, barrier::ScalarField)
    ub = maximum(bounds)
    lb = minimum(bounds)
    widths = maximum(bounds) - minimum(bounds)
    barrier_field = SignedDistanceField(x -> SoftRobots.evaluate(barrier, x), bounds, minimum(widths)/20);
    barrier_mesh = HomogenousMesh(barrier_field, 0.0)
    rescaled_points = Point{3,Float64}[Vec(v-1) ./ (Vec(size(barrier_field))-1) .* (ub - lb) + lb for v in vertices(barrier_mesh)]
     HomogenousMesh(rescaled_points, barrier_mesh.faces)
 end

 function barrier_mesh(field::HermiteRadialField)
    center = Vec(mean(field.points))
    widths = Vec(2*(maximum(field.points) - minimum(field.points)))
    lb = Vec(center - 0.5*widths)
    ub = Vec(center + 0.5*widths)
    bounds = HyperRectangle(lb, widths)
    barrier_mesh(bounds, field)
end

function barrier_mesh(field::InterpolatingSurface)
    center = Vec(mean(field.points))
    widths = Vec(2*(maximum(field.points) - minimum(field.points)))
    lb = Vec(center - 0.5*widths)
    ub = Vec(center + 0.5*widths)
    bounds = HyperRectangle(lb, widths)
    barrier_mesh(bounds, field)
end

barrier_mesh(state::SoftRobotState) = barrier_mesh(state.barrier)

include("drawing.jl")

function copy_dict_values(dict::Dict)
    Dict(zip(keys(dict), deepcopy(values(dict))))
end

function simulate{StateType <: WorldState}(world::World, initial_world_state::StateType, times; draw_callback::Function=(args...)->nothing, history_downsampling=1)
    state = copy_dict_values(initial_world_state)
    history = Tuple{Float64, StateType}[]
    dts = diff(times)

    for i in 2:length(times)
        SoftRobots.update!(world, state, times[i] - times[i-1])
        if mod(i, history_downsampling) == 0
            push!(history, (times[i], copy_dict_values(state)))
        end
        draw_callback(world, state)
    end
    history
end

include("robots.jl")

end


text:
##### Mass tendencies

#####
##### First order fluxes
#####

function flux(::Advect{Mass}, m, state, aux, t, ts, direction)
    return state.ρu
end

#####
##### Second order fluxes
#####

struct MoistureDiffusion{PV <: Mass} <: TendencyDef{Flux{SecondOrder}, PV} end
function flux(
    ::MoistureDiffusion{Mass},
    m,
    state,
    aux,
    t,
    ts,
    diffusive,
    hyperdiff,
)
    ν, D_t, τ = turbulence_tensors(m, state, diffusive, aux, t)
    d_q_tot = (-D_t) .* diffusive.moisture.∇q_tot
    return d_q_tot * state.ρ
end

#####
##### Sources
#####

function source(s::Subsidence{Mass}, m, args)
    @unpack state, aux, diffusive = args
    z = altitude(m, aux)
    w_sub = subsidence_velocity(s, z)
    k̂ = vertical_unit_vector(m, aux)
    return -state.ρ * w_sub * dot(k̂, diffusive.moisture.∇q_tot)
end

function source(s::RemovePrecipitation{Mass}, m, args)
    @unpack state = args
    @unpack ts = args.precomputed
    if has_condensate(ts)
        nt = remove_precipitation_sources(s, m, args)
        return nt.S_ρ_qt
    else
        FT = eltype(state)
        return FT(0)
    end
end

function source(s::WarmRain_1M{Mass}, m, args)
    nt = warm_rain_sources(m, args)
    return nt.S_ρ_qt
end

function source(s::RainSnow_1M{Mass}, m, args)
    nt = rain_snow_sources(m, args)
    return nt.S_ρ_qt
end


text:
"""
    f = initmetagauss(sx, st, nso, nto)

Generate the spectral amplitudes `f` for a metagaussian field of size `nso * nso * nto`
with slopes `sx` and `st`.
"""
function initmetagauss(sx, st, ns, nt)

    sx = abs(sx)
    st = abs(st)

#   This spares a lot of memory!
#   [kx,ky,kt]=ndgrid([0:ns/2 -ns/2+1:-1],[0:ns/2 -ns/2+1:-1],[0:nt/2 -nt/2+1:-1])
    kx = [collect(0:ns/2)' collect(-ns/2+1:-1)']'
    kx = repeat(kx*ones(1,ns), outer=[1,1,nt])
# kky=permdute(kkx,[2 1 3])

if (nt>1)
    kt = [collect(0:nt/2)' collect(-nt/2+1:-1)']'
    kt = repeat(kt*ones(1,ns),outer=[1 1 ns])
    kt = permute(kt,[3 2 1]).^2

    kt[:,:,1] .= 0.000001
end
    kx = kx.^2+permutedims(kx,[2 1 3]).^2

    kx[1,1,:] .= 0.000001
# f = (kx.^(-(sx+1)/4)).*kt.^(-st/4)
    if (nt>1)
        kx = (kx.^(-(sx+1)/4)).*kt.^(-st/4)
    else
        kx = (kx.^(-(sx+1)/4))
    end
    kx[1,1,1] = 0
    kx[1,1,:] .= 0
    if (nt>1)
       kx[:,:,1] .= 0
    end
   #  kx=kx./sqrt(sum(lin(abs(kx).^2)))*ns*ns*nt
   kx = kx./sqrt(sum(abs2.(kx[:])))*ns*ns*nt
end


text:
using LinearAlgebra
using ForwardDiff
using Test

export newton_raphson
export test_newton_raphson

∇(f, x) = ForwardDiff.jacobian(f, x)

function newton_raphson(f::Function, x0::Vector{Float64}, ϵ::Float64)
    xs = [x0]
    d = inv(-∇(f, x0)) * f(x0)
    while norm(d) > ϵ
        x = xs[end] + d
        d = inv(-∇(f, x)) * f(x)
        push!(xs, x)
    end
    return xs
end

function test_newton_raphson()
    f(x) = [x[1]^2 + x[2]^2 - 4*x[3] - 0;
            x[1]^2 + x[3]^2 + 0*x[3] - (1/4);
            x[1]^2 + x[2]^2 + x[3]^2 - 1]
    x0 = [1.0, 1.0, 1.0]
    ϵ = 0.001
    x = newton_raphson(f, x0, ϵ)
    @info "Newton-Raphson iterations:"
    @info x
    @test isapprox(norm(f(x[end])), 0; atol=ϵ)
end


text:
using KernelAbstractions, Test
include(joinpath(@__DIR__, "utils.jl")) # Load backend

@kernel function copy_kernel!(A, @Const(B))
    I = @index(Global)
    @inbounds A[I] = B[I]
end

function mycopy!(A::Array, B::Array)
    @assert size(A) == size(B)
    kernel = copy_kernel!(CPU(), 8)
    kernel(A, B, ndrange=length(A))
end

A = zeros(128, 128)
B = ones(128, 128)
event = mycopy!(A, B)
wait(event)
@test A == B

if has_cuda && has_cuda_gpu()
    function mycopy!(A::CuArray, B::CuArray)
        @assert size(A) == size(B)
        copy_kernel!(CUDADevice(), 256)(A, B, ndrange=length(A))
    end

    A = CuArray{Float32}(undef, 1024)
    B = CUDA.ones(Float32, 1024)
    event = mycopy!(A, B)
    wait(event)
    @test A == B
end

if has_rocm && has_rocm_gpu()
    function mycopy!(A::ROCArray, B::ROCArray)
        @assert size(A) == size(B)
        copy_kernel!(ROCDevice(), 256)(A, B, ndrange=length(A))
    end

    A = zeros(Float32, 1024) |> ROCArray
    B = ones(Float32, 1024) |> ROCArray
    event = mycopy!(A, B)
    wait(event)
    @test A == B
end


text:
# WARNING this file is only loaded if GBIF.jl is also active
# This all happens thanks to the Requires.jl package

import Base: getindex
import Base: setindex!
import SimpleSDMLayers: clip, latitudes, longitudes

"""
    Base.getindex(p::T, occurrence::GBIF.GBIFRecord) where {T <: SimpleSDMLayer}

Extracts the value of a layer at a given position for a `GBIFRecord`. If the
`GBIFRecord` has no latitude or longitude, this will return `nothing`.
"""
function Base.getindex(p::T, occurrence::GBIF.GBIFRecord) where {T <: SimpleSDMLayer}
   ismissing(occurrence.latitude) && return nothing
   ismissing(occurrence.longitude) && return nothing
   return p[occurrence.longitude, occurrence.latitude]
end

"""
    Base.setindex!(p::T, v, occurrence::GBIFRecord) where {T <: SimpleSDMResponse}

Changes the values of the cell including the point at the requested latitude and
longitude.
"""
function Base.setindex!(p::SimpleSDMResponse{T}, v::T, occurrence::GBIF.GBIFRecord) where {T}
   ismissing(occurrence.latitude) && return nothing
   ismissing(occurrence.longitude) && return nothing
   setindex!(p, v, occurrence.longitude, occurrence.latitude)
end

"""
    clip(p::T, r::GBIF.GBIFRecords)

Returns a clipped version (with a 10% margin) around all occurences in a
GBIFRecords collection.
"""
function SimpleSDMLayers.clip(p::T, r::GBIF.GBIFRecords) where {T <: SimpleSDMLayer}
   occ_latitudes = filter(!ismissing, [r[i].latitude for i in 1:length(r)])
   occ_longitudes = filter(!ismissing, [r[i].longitude for i in 1:length(r)])

   lat_min = minimum(occ_latitudes)
   lat_max = maximum(occ_latitudes)
   lon_min = minimum(occ_longitudes)
   lon_max = maximum(occ_longitudes)

   lat_Δ = abs(lat_max - lat_min)
   lon_Δ = abs(lon_max - lon_min)

   scaling = 0.1
   lon_s = scaling*lon_Δ
   lat_s = scaling*lat_Δ

   lat_max = min(p.top, lat_max+lat_s)
   lat_min = max(p.bottom, lat_min-lat_s)
   lon_max = min(p.right, lon_max+lon_s)
   lon_min = max(p.left, lon_min-lon_s)

   return p[left=lon_min, right=lon_max, bottom=lat_min, top=lat_max]
end

"""
    Base.getindex(p::T, r::GBIF.GBIFRecords) where {T <: SimpleSDMLayer}

Returns the values of a layer at all occurrences in a `GBIFRecords` collection.
"""
function Base.getindex(layer::T, records::GBIF.GBIFRecords) where {T <: SimpleSDMLayer}
   return convert(Vector{SimpleSDMLayers._inner_type(layer)}, filter(!isnothing, [layer[records[i]] for i in 1:length(records)]))
end

SimpleSDMLayers.longitudes(record::GBIF.GBIFRecord) = record.longitude
SimpleSDMLayers.latitudes(record::GBIF.GBIFRecord) = record.latitude

"""
    latitudes(records::GBIFRecords)

Returns the non-missing latitudes.
"""
SimpleSDMLayers.latitudes(records::GBIF.GBIFRecords) = filter(!ismissing, [latitudes(record) for i in 1:length(records)])

"""
    longitudes(records::GBIFRecords)

Returns the non-missing longitudes.
"""
SimpleSDMLayers.longitudes(records::GBIF.GBIFRecords) = filter(!ismissing, [longitudes(record) for i in 1:length(records)])


text:
abstract type QOperator end
abstract type QOperatorType end
struct BosonicDestroy <: QOperatorType end
struct BosonicCreate <: QOperatorType end
struct PauliSP <: QOperatorType end
struct PauliSM <: QOperatorType end

mutable struct QOperatorMeta
    type::QOperatorType
    h_idx::Int
end

Base.hash(op::T, h::UInt) where T<:Sym{QOperator} = hash(T, hash(getmetadata(op, QOperatorMeta).type, hash(getmetadata(op, QOperatorMeta).h_idx, h)))

# TermInterface.exprhead(::QOperator) = :call
# TermInterface.istree(::Type{<:QOperator}) = false
# TermInterface.arguments(::Type{<:QOperator}) = []

SymbolicUtils.promote_symtype(f, Ts::Type{<:QOperator}...) = promote_type(Ts...)
SymbolicUtils.promote_symtype(f, T::Type{<:QOperator}, Ts...) = T
SymbolicUtils.promote_symtype(f, T::Type{<:QOperator}, S::Type{<:Number}) = T
SymbolicUtils.promote_symtype(f, T::Type{<:Number}, S::Type{<:QOperator}) = S
SymbolicUtils.promote_symtype(f, T::Type{<:QOperator}, S::Type{<:QOperator}) = promote_type(T,S)

SymbolicUtils.promote_symtype(f, Ts::Type{<:QOperatorType}...) = promote_type(Ts...)
SymbolicUtils.promote_symtype(f, T::Type{<:QOperatorType}, Ts...) = T
SymbolicUtils.promote_symtype(f, T::Type{<:QOperatorType}, S::Type{<:Number}) = T
SymbolicUtils.promote_symtype(f, T::Type{<:Number}, S::Type{<:QOperatorType}) = S
SymbolicUtils.promote_symtype(f, T::Type{<:QOperatorType}, S::Type{<:QOperatorType}) = promote_type(T,S)

SymbolicUtils.islike(::Symbolic{QOperator}, ::Type{Number}) = true

for f in [+, -, *, \, /, ^]
    @eval SymbolicUtils.promote_symtype(::$(typeof(f)), Ts::Type{<:QOperator}...) = promote_type(Ts...)
    @eval SymbolicUtils.promote_symtype(::$(typeof(f)), T::Type{<:QOperator}, Ts...) = T
    @eval SymbolicUtils.promote_symtype(::$(typeof(f)), T::Type{<:QOperator}, S::Type{<:Number}) = T
    @eval SymbolicUtils.promote_symtype(::$(typeof(f)), T::Type{<:Number}, S::Type{<:QOperator}) = S
    @eval SymbolicUtils.promote_symtype(::$(typeof(f)), T::Type{<:QOperator}, S::Type{<:QOperator}) = promote_type(T,S)
end

function Base.isequal(a::Term{QOperator},b::Term{QOperator})
    length(arguments(a))==length(arguments(b)) || return false
    for (arg_a,arg_b) ∈ zip(arguments(a), arguments(b))
        isequal(arg_a, arg_b) || return false
    end
    return true
end

function commutator(a, b)
    return normal_order(normal_order(a * b) - normal_order(b * a))
end

macro boson(xs...)
    global hilbert_space_iterator
    defs = map(enumerate(xs)) do (i, x)
        n, t = _name_type(x)
        T = esc(t)
        nt = _name_type(x)
        n, t = nt.name, nt.type
        :($(esc(n)) = Sym{QOperator}($(Expr(:quote, n))); $(esc(n)) = setmetadata($(esc(n)), QOperatorMeta, QOperatorMeta(BosonicDestroy(), hilbert_space_iterator + $(esc(i)))))
    end
    hilbert_space_iterator += length(xs)
    Expr(:block, defs...,
         :(tuple($(map(x->esc(_name_type(x).name), xs)...))))
end

macro parameter(xs...)
    defs = map(enumerate(xs)) do (i, x)
        n, t = _name_type(x)
        T = esc(t)
        nt = _name_type(x)
        n, t = nt.name, nt.type
        :($(esc(n)) = Sym{$T}($(Expr(:quote, n))))
    end
    Expr(:block, defs...,
         :(tuple($(map(x->esc(_name_type(x).name), xs)...))))
end

# macro paulism(xs...)
#     global hilbert_space_iterator
#     defs = map(enumerate(xs)) do (i, x)
#         n, t = _name_type(x)
#         T = esc(t)
#         nt = _name_type(x)
#         n, t = nt.name, nt.type
#         :($(esc(n)) = Sym{QOperator}($(Expr(:quote, n)), metadata = QOperatorMeta(PauliSM(), hilbert_space_iterator + $(esc(i)))))
#     end
#     hilbert_space_iterator += length(xs)
#     Expr(:block, defs...,
#          :(tuple($(map(x->esc(_name_type(x).name), xs)...))))
# end


text:
ALL_TASK_OFFSETS[("linux-x86_64", v"1.5.4")] = Dict(
    :END => 608,
    :backtrace => 56,
    :bufsz => 296,
    :copy_stack => 304,
    :ctx => 88,
    :donenotify => 32,
    :eh => 312,
    :exception => 48,
    :excstack => 328,
    :gcstack => 320,
    :locks => 344,
    :logstate => 64,
    :next => 0,
    :prio => 338,
    :queue => 8,
    :result => 40,
    :sizeof_ctx => 200,
    :start => 72,
    :state => 24,
    :sticky => 80,
    :stkbuf => 288,
    :tid => 336,
    :timing_stack => 600,
    :tls => 16,
    :tls_base_context => 6648,
)


text:
function Simulation(
	ω₀::Real,
	domains::Vararg{Domain{1},N};
	k₂₀::Real=0,
	k₃₀::Real=0,
	k₁₀::Real=sqrt(ω₀^2-k₂₀^2-k₃₀^2)
	) where N

	lattices = map(z->z.lattice,domains)
	all(isequal.(Ref(lattices[1]),lattices)) || throw("in 1D, lattice must be shared by all domains")

	starts = map(d->d.boundary.shape.start,domains)
	stops =  map(d->d.boundary.shape.stop ,domains)
	start = min(starts...,stops...)
	stop  = max(starts...,stops...)
	a::Float64 = stop-start
	nx = ceil(Int,a/domains[1].lattice.dx)
	dx = a/nx

	lat = Lattice(dx;origin=start+dx/2)

	domains = map(d->Domain(d.boundary,lat,d.dielectric,d.pump,d.χ;type=d.type,name=d.name),domains)

	smoothed = Ref(false)

	domain_index = Int[]
	indices = CartesianIndex{1}[]
	nnm = (Int[],)
	nnp = (Int[],)
	name = Symbol[]
	x = Point{1}[]
	ε = ComplexF64[]
	F = Float64[]
	interior = BitArray([])
	exterior = falses(length(x))
	surface = BitArray([])
	bulk = BitArray([])
	corner = BitArray([])

	append_domains!(domain_index,nnm,nnp,name,x,indices,ε,F,interior,bulk,surface,corner,domains)

	removed = get_sites_to_remove(domains,domain_index,x)

	trim_nn!(nnm,nnp,removed,bulk)

	remove_sites!((domain_index,name,x,indices,ε,F,interior,bulk,surface,corner,nnm,nnp),removed)

	x_half = generate_half_x(domains,domain_index,x,indices)

	domain_wall = generate_and_domain_wall!(interior,bulk,surface,corner,domains,domain_index,x,indices)

	σ = boundary_layer(domains,x,domain_index)
	σ_half = boundary_layer(domains,x_half[1],vcat(domain_index,domain_index[end]))
	α = 1 .+ 1im*σ[1]/k₁₀
	α_half = 1 .+ 1im*σ_half[1]/k₁₀

	curlcurl = Curlcurl(domains[1].lattice,α,α_half,nnm,nnp,indices,interior,surface,domain_wall)
	Σ = SelfEnergy(domains,surface,domain_index,α_half,a,nnm,nnp,indices,interior,domain_wall)

	return Simulation{1,Symmetric,ComplexF64,typeof(domains),2,typeof(Σ.f)}(
		smoothed,
		ε,
		F,
		curlcurl,
		Σ,
		domains,
		domain_index,
		x,
		x_half,
		(α,),
		σ,
		σ_half,
		name,
		ω₀,
		k₁₀,
		k₂₀,
		k₃₀)
end

function generate_and_domain_wall!(interior,bulk,surface,corner,domains,domain_index,x::Array{Point{1}},indices)
	lattices = map(z->z.lattice,domains)
	shapes = map(z->z.boundary.shape,domains)
	DX = map(z->z.dx,lattices)
	domain_wall = falses(length(x))
	lattice_wall = falses(length(x))
	for i ∈ eachindex(x)
		if interior[i]
			ds = domain_index[i]
			xt = lattices[ds][indices[i]-CartesianIndex(1,)]
			dom1 = which_domain(domains,shapes,xt)
			xt = lattices[ds][indices[i]+CartesianIndex(1,)]
			dom2 = which_domain(domains,shapes,xt)

			domain_wall[i] = !(domain_index[i]==dom1==dom2)
			domain_wall[i] = domain_wall[i] && (dom1!==0 && dom2!==0)
			if domain_wall[i]
				corner[i] = false
				surface[i] = false
				bulk[i] = true
			end
		end
	end
	return domain_wall
end

function generate_half_x(domains,domain_index,x::Array{Point{1}},indices)
	half_x = (Array{Point{1},1}(undef,length(x)+1),)
	lattices = map(z->z.lattice,domains)
	for i ∈ 1:length(x)
		ds = domain_index[i]
		half_x[1][i] = lattices[ds][indices[i][1]-.5]
		if i==length(x)
			half_x[1][i+1] = lattices[ds][indices[i][1]+.5]
		end
	end
	return half_x
end

function smooth!(sim::Simulation{1},num_sub_pixel::Int=NUM_SUB_PIXEL)
    domains = sim.domains
    domain_wall = sim.domain_wall
	lattices = map(z->z.lattice,sim.domains)
    X = Array{Point{1}}(undef,num_sub_pixel)
    E = Array{ComplexF64}(undef,num_sub_pixel)
    F = Array{Float64}(undef,num_sub_pixel)
	dw_inds = findall(domain_wall)
    for dw ∈ eachindex(dw_inds)
		d = dw_inds[dw]
		ds = sim.domain_index[d]
	    I = Tuple(sim.indices[d])
        is = LinRange(I-.5,I+.5,num_sub_pixel)
        for i ∈ eachindex(is), j ∈ eachindex(js)
            X[i] = lattices[ds][is[i]]
        end
        simulation_dielectric!(E,domains,X)
        simulation_pump!(F,domains,X)
        sim.ε[d] = mean(E)
        sim.F[d] = mean(F)
    end
	sim.smoothed[] = true
    return nothing
end

function boundary_layer(domain,x::Point{1},dim::Int)
	typeof(domain.boundary.shape)<:Interval || throw("1D only defined for shape=Interval, but shape=",domain.boundary.shape)
	N = 2
	bl1 = domain.boundary.bls[1]
	bl2 = domain.boundary.bls[2]
	σx = bl1(x)+bl2(x)
	return σx
end

function Base.conj(sim::Simulation)
	return Simulation(sim.ω₀,map(conj,sim.domains)...)
end


text:
using DrWatson
@quickactivate
using ArgParse
using GenerativeAD
using StatsBase: fit!, predict
using StatsBase
using BSON

s = ArgParseSettings()
@add_arg_table! s begin
   "max_seed"
		default = 1
		arg_type = Int
		help = "seed"
	"dataset"
		default = "MNIST"
		arg_type = String
		help = "dataset"
	"anomaly_classes"
		arg_type = Int
		default = 10
		help = "number of anomaly classes"
	"method"
		arg_type = String
		default = "leave-one-out"
		help = "method for data creation -> \"leave-one-out\" or \"leave-one-in\" "
    "contamination"
    	arg_type = Float64
    	help = "contamination rate of training data"
    	default = 0.0
end
parsed_args = parse_args(ARGS, s)
@unpack dataset, max_seed, anomaly_classes, method, contamination = parsed_args

modelname = "fAnoGAN-GP"

"""
	sample_params()

Should return a named tuple that contains a sample of model parameters.
"""
function sample_params()
	# first sample the number of layers
	nlayers = rand(2:4)
	kernelsizes = reverse((3,5,7,9)[1:nlayers])
	channels = reverse((16,32,64,128)[1:nlayers])
	scalings = reverse((1,2,2,2)[1:nlayers])

	par_vec = (
		2 .^(3:8),
		10f0 .^(-4:-3),
		10f0 .^(-4:-3),
		2 .^ (5:7),
		["relu", "swish", "tanh"],
		1:Int(1e8),
	)
	argnames = (
		:zdim,
		:lr_gan,
		:lr_enc,
		:batch_size,
		:activation,
		:init_seed,
	)
	parameters = (;zip(argnames, map(x->sample(x, 1)[1], par_vec))...)
	return merge(parameters, (nlayers=nlayers, kernelsizes=kernelsizes,
		channels=channels, scalings=scalings))
end

"""
	fit(data, parameters)

This is the most important function - returns `training_info` and a tuple or a vector of tuples `(score_fun, final_parameters)`.
`training_info` contains additional information on the training process that should be saved, the same for all anomaly score functions.
Each element of the return vector contains a specific anomaly score function - there can be multiple for each trained model.
Final parameters is a named tuple of names and parameter values that are used for creation of the savefile name.
"""
function fit(data, parameters)
	# construct model - constructor should only accept kwargs
	default_params = (lr_gan = 0.00005, lr_enc = 0.001, batch_size=128, max_iters=10000, n_critic=5, usegpu=true)

	@info "idim -> $(size(data[1][1])[1:3]), flipped -> $(reverse(size(data[1][1])[1:3]))"
	# construct model - constructor should only accept kwargs
	model = GenerativeAD.Models.fAnoGAN_GP(;idim = reverse(size(data[1][1])[1:3]), usegpu=true, parameters...)

	#max_iter = 50 # this should be enough for testing purpouses
	params = merge(default_params, parameters)
	# fit train data

	try
		global info, fit_t, _, _, _ = @timed StatsBase.fit!(model, data[1][1],
			max_iters=params.max_iters, lr_gan=params.lr_gan, lr_enc=params.lr_enc,
			batch_size=params.batch_size, n_critic=params.n_critic)
	catch e
		# return an empty array if fit fails so nothing is computed
		@info "Failed training due to \n$e"
		return (fit_t = NaN, history=nothing, npars=nothing, model=nothing), []
	end

	#println("printing return form fit! -> ", info)
	model = info[1]
	#println(model)

	# construct return information - put e.g. the model structure here for generative models
	training_info = (
		fit_t = fit_t,
		history = info[2],
		model = model  # pytorchmodel
		)

	# now return the different scoring functions
	training_info, [(x -> predict(model, x), parameters)]
end

#___________________________________________________________________________________________________________________
try_counter = 0
max_tries = 10*max_seed
cont_string = (contamination == 0.0) ? "" : "_contamination-$contamination"
while try_counter < max_tries
	parameters = sample_params()

	for seed in 1:max_seed
		for i in 1:anomaly_classes
			savepath = datadir("experiments/images_$(method)$cont_string/$(modelname)/$(dataset)/ac=$(i)/seed=$(seed)")
			mkpath(savepath)

			data = GenerativeAD.load_data(dataset, seed=seed, method=method, anomaly_class_ind=i)

			@info "Trying to fit $modelname on $dataset with parameters $(parameters)..."
			@info "Train/validation/test splits: $(size(data[1][1], 4)) | $(size(data[2][1], 4)) | $(size(data[3][1], 4))"
			@info "Number of features: $(size(data[1][1])[1:3])"

			# here, check if a model with the same parameters was already tested
			if GenerativeAD.check_params(savepath, parameters) & GenerativeAD.Models.check_scaling(size(data[1][1])[1:3], parameters.scalings)
				# fit model
				training_info, results = fit(data, parameters)
				# saving model separately
				if training_info.model !== nothing
					# save model with build-in save funciton (torch.save) as .pt file (pt ~ pytorch)
					training_info.model.model.save_model(joinpath(savepath, savename("model", parameters, digits=5)))
					training_info = merge(training_info, (model = nothing,))
				end
				save_entries = merge(training_info, (modelname = modelname, seed = seed, dataset = dataset, anomaly_class = i,
					contamination=contamination))
				# now loop over all anomaly score funs
				for result in results
					GenerativeAD.experiment(result..., data, savepath; save_entries...)
				end
				global try_counter = max_tries + 1
			else
				@info "Model already present, sampling new hyperparameters..."
				global try_counter += 1
			end
		end
	end
end
(try_counter == max_tries) ? (@info "Reached $(max_tries) tries, giving up.") : nothing