text:
"""
    methodlearner(method, traindata, validdata[; callbacks=[], kwargs...]) -> Learner
    methodlearner(method, data; pctgval = 0.2, kwargs...)

Create a [`Learner`](#) to train a model for learning method `method` using
`data`.

## Keyword arguments

- `callbacks = []`: [`Callback`](#)s to use during training.
- `batchsize = 16`: Batch size used for the training data loader.
- `backbone = nothing`: Backbone model to construct task-specific model from using
   [`methodmodel`](#)`(method, backbone)`.
- `model = nothing`: Complete model to use. If given, the `backbone` argument is ignored.
- `optimizer = ADAM()`: Optimizer passed to `Learner`.
- `lossfn = `[`methodlossfn`](#)`(method)`: Loss function passed to `Learner`.

Any other keyword arguments will be passed to [`methoddataloaders`](#).

## Examples

Full example:

```julia
data, blocks = loaddataset("imagenette2-160", (Image, Label))
method = ImageClassificationSingle(blocks)
learner = methodlearner(method, data)
fitonecycle!(learner, 10)
```

Custom training and validation split:

```julia
learner = methodlearner(method, traindata, validdata)
```

Using callbacks:

```julia
learner = methodlearner(method, data; callbacks=[
    ToGPU(), Checkpointer(), LogMetrics(TensorboardBackend())
])
```
"""
function methodlearner(
        method::LearningMethod,
        traindata,
        validdata;
        backbone=nothing,
        model=nothing,
        callbacks=[],
        pctgval=0.2,
        batchsize=16,
        optimizer=ADAM(),
        lossfn=methodlossfn(method),
        kwargs...,
    )
    if isnothing(model)
        model = isnothing(backbone) ? methodmodel(method) : methodmodel(method, backbone)
    end
    dls = methoddataloaders(traindata, validdata, method, batchsize; kwargs...)
    return Learner(model, dls, optimizer, lossfn, callbacks...)
end

function methodlearner(method, data; pctgval=0.2, kwargs...)
    traindata, validdata = splitobs(shuffleobs(data), at=1 - pctgval)
    return methodlearner(method, traindata, validdata; kwargs...)
end

function plotpredictions(method::LearningMethod, learner::Learner; n=4)
    cb = FluxTraining.getcallback(learner, ToDevice)
    devicefn = isnothing(cb) ? identity : cb.movedatafn
    backfn = isnothing(cb) ? identity : cpu
    xs, ys = first(learner.data.validation)
    b = last(size(xs))
    xs, ys = DataLoaders.collate([s for (s, _) in zip(DataLoaders.obsslices((xs, ys)), 1:min(n, b))])
    ŷs = learner.model(devicefn(xs)) |> backfn
    return plotpredictions(method, xs, ŷs, ys)
end


text:
# A dummy struct used to easily collect alpha values for paletted images. It allows us to simply
# call `pointer_from_objref` on a `Vector{_AlphaBuffer}`
struct _AlphaBuffer
    val::N0f8
end

function _inspect_png_read(fpath, gamma::Union{Nothing,Float64}=nothing)
    fp = open_png(fpath)
    png_ptr = create_read_struct()
    info_ptr = create_info_struct(png_ptr)
    png_init_io(png_ptr, fp)
    png_set_sig_bytes(png_ptr, PNG_BYTES_TO_CHECK)
    png_read_info(png_ptr, info_ptr)

    width = png_get_image_width(png_ptr, info_ptr)
    height = png_get_image_height(png_ptr, info_ptr)
    color_type_orig = png_get_color_type(png_ptr, info_ptr)
    color_type = color_type_orig
    bit_depth_orig = png_get_bit_depth(png_ptr, info_ptr)
    bit_depth = bit_depth_orig
    num_channels = png_get_channels(png_ptr, info_ptr)
    interlace_type = png_get_interlace_type(png_ptr, info_ptr)

    backgroundp = png_color_16p()
    if png_get_bKGD(png_ptr, info_ptr, Ref(backgroundp)[]) != 0
        png_set_background(png_ptr, backgroundp, PNG_BACKGROUND_GAMMA_FILE, 1, 1.0)
    end

    if color_type == PNG_COLOR_TYPE_PALETTE
        png_set_palette_to_rgb(png_ptr)
        color_type = PNG_COLOR_TYPE_RGB
    end

    if color_type == PNG_COLOR_TYPE_GRAY && bit_depth < 8
        png_set_expand_gray_1_2_4_to_8(png_ptr)
        png_set_packing(png_ptr)
        bit_depth = UInt8(8)
    end

    if png_get_valid(png_ptr, info_ptr, PNG_INFO_tRNS) != 0
        png_set_tRNS_to_alpha(png_ptr)
        if color_type == PNG_COLOR_TYPE_GRAY || color_type == PNG_COLOR_TYPE_RGB
            color_type |= PNG_COLOR_MASK_ALPHA
        end
    end

    screen_gamma = PNG_DEFAULT_sRGB
    image_gamma = Ref{Cdouble}(0.0)
    intent = Ref{Cint}(-1)
    white_x = Ref{Cdouble}(0.0)
    white_y = Ref{Cdouble}(0.0)
    red_x = Ref{Cdouble}(0.0)
    red_y = Ref{Cdouble}(0.0)
    green_x = Ref{Cdouble}(0.0)
    green_y = Ref{Cdouble}(0.0)
    blue_x = Ref{Cdouble}(0.0)
    blue_y = Ref{Cdouble}(0.0)
    if isnothing(gamma)
        if png_get_valid(png_ptr, info_ptr, PNG_INFO_sRGB) != 0
            if png_get_sRGB(png_ptr, info_ptr, intent) != 0
                png_set_gamma(png_ptr, screen_gamma, PNG_DEFAULT_sRGB);
            else
                if png_get_gAMA(png_ptr, info_ptr, image_gamma) != 0
                    png_set_gamma(png_ptr, screen_gamma, image_gamma[])
                else
                    image_gamma[] = 0.45455
                    png_set_gamma(png_ptr, screen_gamma, image_gamma[])
                end
            end
        elseif png_get_valid(png_ptr, info_ptr, PNG_INFO_gAMA) != 0
            if png_get_gAMA(png_ptr, info_ptr, image_gamma) != 0
                png_set_gamma(png_ptr, screen_gamma, image_gamma[])
            else
                image_gamma[] = 0.45455
                png_set_gamma(png_ptr, screen_gamma, image_gamma[])
            end
        end
    elseif gamma != 1
        image_gamma[] = gamma
        png_set_gamma(png_ptr, screen_gamma, image_gamma[])
    end

    buffer_eltype = _buffer_color_type(color_type, bit_depth)
    bit_depth == 16 && png_set_swap(png_ptr)
    n_passes = png_set_interlace_handling(png_ptr)
    png_read_update_info(png_ptr, info_ptr)
    chunk_info(chunk) = png_get_valid(png_ptr, info_ptr, chunk)

    @show(fpath,
        chunk_info(PNG_INFO_sRGB),
        chunk_info(PNG_INFO_iCCP),
        chunk_info(PNG_INFO_gAMA),
        chunk_info(PNG_INFO_tRNS),
        chunk_info(PNG_INFO_cHRM),
        chunk_info(PNG_INFO_PLTE),
        chunk_info(PNG_INFO_sPLT),
        chunk_info(PNG_INFO_hIST),
        chunk_info(PNG_INFO_sBIT),
        chunk_info(PNG_INFO_bKGD),
        chunk_info(PNG_INFO_pCAL),
        (height, width, num_channels),
        (color_type_orig, color_type),
        (bit_depth_orig, bit_depth),
        (image_gamma[], screen_gamma),
        intent[],
        (white_x[], white_y[], red_x[], red_y[], green_x[], green_y[], blue_x[], blue_y[]),
        (interlace_type, n_passes),
        buffer_eltype,
    )

    png_destroy_read_struct(Ref{Ptr{Cvoid}}(png_ptr), Ref{Ptr{Cvoid}}(info_ptr), C_NULL)
    close_png(fp)
end

ext: jl


text:
using StochasticBlockModel
using Test

@testset "Random assignemnt" begin

    A = [0 1 1 0 0 0 0 0;
         1 0 0 1 0 0 0 0;
         1 0 0 1 0 0 0 0;
         0 1 1 0 0 0 0 0;
         0 0 0 0 0 1 1 0;
         0 0 0 0 1 0 0 1;
         0 0 0 0 1 0 0 1;
         0 0 0 0 0 1 1 0]
    dataset = StochasticBlockModel.Dataset(A, 2)
    x = StochasticBlockModel.randomAssignments(dataset, seed=1)
    @test x == [1 0; 1 0; 0 1; 0 1; 0 1; 1 0; 1 0; 0 1]
    x = StochasticBlockModel.randomAssignments(dataset, seed=2)
    @test x == [1 0; 0 1; 0 1; 0 1; 0 1; 1 0; 0 1; 1 0]

end

@testset "Local Search Assignments" begin

    A = [0 1 1 0 0 0 0 0;
         1 0 0 1 0 0 0 0;
         1 0 0 1 0 0 0 0;
         0 1 1 0 0 0 0 0;
         0 0 0 0 0 1 1 0;
         0 0 0 0 1 0 0 1;
         0 0 0 0 1 0 0 1;
         0 0 0 0 0 1 1 0]
    dataset = StochasticBlockModel.Dataset(A, 2)
    x = StochasticBlockModel.randomAssignments(dataset, seed=1)
    estimator = StochasticBlockModel.SBMEstimator("ls1", 10.0, false, 0, false)
    w = StochasticBlockModel.optimalProbMatrix(dataset, x)
    obj_before = StochasticBlockModel.calculateObjective(dataset, w, x)
    # Run local search
    x_ = StochasticBlockModel.localSearchAssignments(estimator, dataset, w, x)
    w_ = StochasticBlockModel.optimalProbMatrix(dataset, x_)
    obj_after = StochasticBlockModel.calculateObjective(dataset, w_, x_)
    @test x_ == [0 1; 1 0; 1 0; 0 1; 0 1; 1 0; 1 0; 0 1]
    @test obj_after ≈ 2.454822555520437 rtol = 1e-5

end

@testset "Local Search Methods" begin

    A = [0 1 1 0 0 0 0 0;
         1 0 0 1 0 0 0 0;
         1 0 0 1 0 0 0 0;
         0 1 1 0 0 0 0 0;
         0 0 0 0 0 1 1 0;
         0 0 0 0 1 0 0 1;
         0 0 0 0 1 0 0 1;
         0 0 0 0 0 1 1 0]
    dataset = StochasticBlockModel.Dataset(A, 2)
    estimator = StochasticBlockModel.SBMEstimator("ls1", 10.0, false, 0, false)
    sbm, x, opt_results = StochasticBlockModel.localSearch1(estimator, dataset)
    @test opt_results.UB ≈ 2.454822555520437 rtol = 1e-5
    @test opt_results.LB == -Inf

    estimator = StochasticBlockModel.SBMEstimator("ls2", 10.0, false, 0, false)
    sbm, x, opt_results = StochasticBlockModel.localSearch2(estimator, dataset)
    @test opt_results.UB ≈ 2.454822555520437 rtol = 1e-5
    @test opt_results.LB == -Inf

    dataset = StochasticBlockModel.Dataset("../instances/RW/zachary.in")
    estimator = StochasticBlockModel.SBMEstimator("ls1", 10.0, false, 0, false)
    sbm, x, opt_results = StochasticBlockModel.localSearch1(estimator, dataset)
    time_ls1 = opt_results.solvetime
    @test opt_results.UB ≈ 67.88545219325196 rtol = 1e-5
    @test opt_results.LB == -Inf
    @test opt_results.status == :LocalOptimum_LS1
    @test sbm.w ≈ [0.57483 1.49603; 1.49603 0.421296] rtol = 1e-5

    estimator = StochasticBlockModel.SBMEstimator("ls2", 10.0, false, 0, false)
    sbm, x, opt_results = StochasticBlockModel.localSearch2(estimator, dataset)
    @test opt_results.UB ≈ 63.973982362428586 rtol = 1e-5
    @test opt_results.LB == -Inf
    @test opt_results.status == :LocalOptimum_LS2
    # @test time_ls1 < opt_results.solvetime
    @test sbm.w ≈ [0.518201 1.5768; 1.5768 0.309462] rtol = 1e-5

end

@testset "Local Search 2" begin

    A = [0 1 1 0 0 0 0 0;
         1 0 0 1 0 0 0 0;
         1 0 0 1 0 0 0 0;
         0 1 1 0 0 0 0 0;
         0 0 0 0 0 1 1 0;
         0 0 0 0 1 0 0 1;
         0 0 0 0 1 0 0 1;
         0 0 0 0 0 1 1 0]
    dataset = StochasticBlockModel.Dataset(A, 2)
    estimator = StochasticBlockModel.SBMEstimator("ls2", 10.0, false, 0, false)
    sbm, x, opt_results = StochasticBlockModel.localSearch2(estimator, dataset)
    @test opt_results.UB ≈ 2.454822555520437 rtol = 1e-5
    @test opt_results.LB == -Inf
end


text:
using ADatasets
using EvalCurves

"""
    gridSearch(f, parameters...)
Maps `f` to product of `parameters`.
"""
gridsearch(f, parameters...) = map(p -> f(p), Base.product(parameters...))

# function printandrun(f, p)
#     println(p)
#     (p, f(p))
# end

global mainfolder = "D:/dev/julia/"
if !isdir(mainfolder)
    mainfolder = "/home/bimjan/dev/julia/"
    if !isdir(mainfolder)
        error("The main folder is unknown")
    end
end
const datafolder = mainfolder * "data/loda/public/datasets/numerical"

loaddata(dataset, difficulty) =  ADatasets.makeset(ADatasets.loaddataset(dataset, difficulty, datafolder)..., 0.8, "low")

computeauc(score, labels) = EvalCurves.auc(EvalCurves.roccurve(score, labels)...)


text:
################################################################################
#
#   fista : Accelerated Forward Backward splitting
#
################################################################################

function fista!{T<:FloatingPoint}(
    rmodel::SupervisedRiskModel,     # model definition (prediction model and loss function)
    θ::StridedArray{T},              # initial parameter guess
    x::StridedArray{T},              # inputs
    y::StridedArray{T};              # responses
    reg::Regularizer = ZeroReg(),    # regularizer
    λ::T             = one(T),       # step size
    cbinterval::Int  = 0,            # callback interval
    callback         = fista_trace,    # callback function
    maxsteps::Int    = 100,          # max number of steps
    maxtime::Int     = 0         )   # max calculation time (0 = no limit)

    g      = similar(θ)  # holds gradient
    dualθ  = copy(θ)     # dual θ
    prevθ  = copy(θ)     # θ of previous step
    deltaθ = similar(θ)  # variation of theta

    function fbstep!(θ1, g, θ0, λ)
        copy!(θ1, θ0)
        axpy!(-λ, g, θ1)
        prox!(reg, θ1, θ1, λ)
    end

    α    = 1.      # acceleration parameter
    starttime = time()
    for t in 1:maxsteps
        if (maxtime != 0) && (time() - starttime > maxtime)
            break
        end

        v_risk, _ = value_and_addgrad!(rmodel, 0., g, 1., dualθ, x, y)
        fbstep!(θ, g, dualθ, λ)
        deltaθ = θ - prevθ
        if dot(g, deltaθ) >= 0.
            α = 1.
            μ = 0.
        else
            prevα = α
            α = 1. + sqrt(1+4α*α)/2.
            μ = (prevα - 1.) / α
        end
        copy!(dualθ, θ)
        axpy!(μ, deltaθ, dualθ)
        copy!(prevθ, θ)

        # callback
        if cbinterval > 0 && t % cbinterval == 0
            callback(t, θ, v_risk, value(reg, θ), α)
        end
    end

    θ
end

# Non mutating version
function fista{T<:FloatingPoint}(
    rmodel::SupervisedRiskModel,     # model definition
    θ::StridedArray{T},              # initial parameter guess
    x::StridedArray{T},              # inputs
    y::StridedArray{T};              # responses
    vargs... )

    fista!(rmodel, copy(θ), x, y; vargs...)
end


text:
using MCMCBenchmarks,CSV

ProjDir = @__DIR__
cd(ProjDir)
folder = "2019_05_30T19_40_00"
results = CSV.read(string(ProjDir,"/results/",folder*"/results.csv"))

# Make plots
include("plotting.jl")


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

name = "Attr"
version = v"2.4.48"

# Collection of sources required to build attr
sources = [
    "https://download.savannah.gnu.org/releases/attr/attr-$(version).tar.gz" =>
    "5ead72b358ec709ed00bbf7a9eaef1654baad937c001c044fe8b74c57f5324e7",
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/attr-*/
./configure --prefix=${prefix} --host=${target}
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line.  We are manually disabling
# many platforms that do not seem to work.
platforms = [p for p in supported_platforms() if Sys.islinux(p)]

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

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

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: blenessy/Yggdrasil

max_issues_repo_name: blenessy/Yggdrasil

max_forks_repo_name: blenessy/Yggdrasil


text:
# Doesn't work for colinear/parallel lines. (𝘂 × 𝘃 = 𝟬).
# For 𝗹₁(r) = 𝘅₁ + r𝘂 and 𝗹₂(s) = 𝘅₂ + s𝘃
# 1) 𝘅₁ + r𝘂 = 𝘅₂ + s𝘃                  subtracting 𝘅₁ from both sides
# 2) r𝘂 = (𝘅₂-𝘅₁) + s𝘃                  𝘄 = 𝘅₂-𝘅₁
# 3) r𝘂 = 𝘄 + s𝘃                        cross product with 𝘃 (distributive)
# 4) r(𝘂 × 𝘃) = 𝘄 × 𝘃 + s(𝘃 × 𝘃)        𝘃 × 𝘃 = 𝟬
# 5) r(𝘂 × 𝘃) = 𝘄 × 𝘃                   let 𝘄 × 𝘃 = 𝘅 and 𝘂 × 𝘃 = 𝘇
# 6) r𝘇 = 𝘅                             dot product 𝘇 to each side
# 7) r𝘇 ⋅ 𝘇 = 𝘅 ⋅ 𝘇                     divide by 𝘇 ⋅ 𝘇
# 8) r = (𝘅 ⋅ 𝘇)/(𝘇 ⋅ 𝘇)
# We need to ensure r, s ∈ [0, 1], hence we need to solve for s too.
# 1) 𝘅₂ + s𝘃 = 𝘅₁ + r𝘂                     subtracting 𝘅₂ from both sides
# 2) s𝘃 = -𝘄 + r𝘂                          cross product with 𝘄
# 3) s(𝘃 × 𝘄) = -𝘄 × 𝘄 + r(𝘂 × 𝘄)          𝘄 × 𝘄 = 𝟬
# 4) s(𝘃 × 𝘄) = r(𝘂 × 𝘄)                   using 𝘂 × 𝘄 = -(𝘄 × 𝘂), likewise for 𝘃 × 𝘄
# 5) s(𝘄 × 𝘃) = r(𝘄 × 𝘂)                   let 𝘄 × 𝘂 = 𝘆. use 𝘄 × 𝘃 = 𝘅
# 6) s𝘅 = r𝘆                               dot product 𝘅 to each side
# 7) s(𝘅 ⋅ 𝘅) = r(𝘆 ⋅ 𝘅)                   divide by (𝘅 ⋅ 𝘅)
# 9) s = r(𝘅 ⋅ 𝘆)/(𝘅 ⋅ 𝘅)
# The cross product of two vectors in the plane is a vector of the form (0, 0, k),
# hence, in 2D:
# r = (𝘅 ⋅ 𝘇)/(𝘇 ⋅ 𝘇) = x₃/z₃
# s = r(𝘅 ⋅ 𝘆)/(𝘅 ⋅ 𝘅) = y₃/z₃
function Base.intersect(l₁::LineSegment{Point{2,T}},
                        l₂::LineSegment{Point{2,T}}) where {T}
    𝘄 = l₂[1] - l₁[1]
    𝘂₁= l₁[2] - l₁[1]
    𝘂₂= l₂[2] - l₂[1]
    z = 𝘂₁ × 𝘂₂
    r = (𝘄 × 𝘂₂)/z
    s = (𝘄 × 𝘂₁)/z
    valid = 0 ≤ r && r ≤ 1 && 0 ≤ s && s ≤ 1
    return valid ? l₂(s) : Point(T(1e6),T(1e6))
end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
module TestPaperSelfgravitatingGasDynamics

using Test
using Trixi

include("test_trixi.jl")

# Start with a clean environment: remove Trixi output directory if it exists
outdir = "out"
isdir(outdir) && rm(outdir, recursive=true)

# pathof(Trixi) returns /path/to/Trixi/src/Trixi.jl, dirname gives the parent directory
const EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "paper-self-gravitating-gas-dynamics")

# Numerical examples from the Euler-gravity paper
@testset "paper-self-gravitating-gas-dynamics" begin
  @testset "elixir_euler_eoc_test.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_eoc_test.jl"),
      l2   = [0.0001740977055972079, 0.0003369355182519592, 0.0003369355182518708, 0.0006099171220334989],
      linf = [0.001079347149189669, 0.0018836938381321389, 0.001883693838132583, 0.003971575376718217])
  end

  @testset "elixir_euler_eoc_test.jl with polydeg=4" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_eoc_test.jl"),
      l2   = [1.7187201161597772e-5, 2.678065111772951e-5, 2.678065111783027e-5, 4.952504160091526e-5],
      linf = [0.0001501749544159381, 0.00016549482504535362, 0.00016549482504601976, 0.0004372960291432193],
      polydeg = 4)
  end

  @testset "elixir_hypdiff_eoc_test.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_hypdiff_eoc_test.jl"),
      l2   = [0.003154024896093942, 0.012394432074951856, 0.02185973823794725],
      linf = [0.01731850928579215, 0.07843510773347553, 0.11242300176349201])
  end

  @testset "elixir_hypdiff_eoc_test.jl with polydeg=4" begin
  @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_hypdiff_eoc_test.jl"),
    l2   = [0.0002511283012128458, 0.0008808243846610255, 0.0016313343228567005],
    linf = [0.0017290715087938668, 0.003129184465704738, 0.01000728849316701],
    polydeg = 4)
  end

  @testset "elixir_eulergravity_eoc_test.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulergravity_eoc_test.jl"),
      l2   = [0.00024871265138964204, 0.0003370077102132591, 0.0003370077102131964, 0.0007231525513793697],
      linf = [0.0015813032944647087, 0.0020494288423820173, 0.0020494288423824614, 0.004793821195083758],
      tspan = (0.0, 0.1))
  end

  @testset "elixir_eulergravity_eoc_test.jl with polydeg=4" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulergravity_eoc_test.jl"),
      l2   = [1.9537712148648045e-5, 2.7564396197947587e-5, 2.7564396197967635e-5, 5.688838772067586e-5],
      linf = [0.00012335710672761735, 0.00020086268350816283, 0.00020086268350727465, 0.0004962155455632278],
      tspan = (0.0, 0.1), polydeg = 4)
  end

  @testset "elixir_eulergravity_eoc_test.jl with 1st order RK3S*" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulergravity_eoc_test.jl"),
      l2   = [0.00024871265138959434, 0.000337007710281087, 0.0003370077102811394, 0.0007231525515231289],
      linf = [0.0015813032941613958, 0.002049428843978518, 0.0020494288439798503, 0.004793821198143977],
      tspan = (0.0, 0.1), timestep_gravity=Trixi.timestep_gravity_erk51_3Sstar!)
  end

  @testset "elixir_eulergravity_eoc_test.jl with 3rd order RK3S*" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulergravity_eoc_test.jl"),
      l2   = [0.0002487126513894034, 0.00033700771023049785, 0.00033700771023048245, 0.0007231525514158737],
      linf = [0.0015813032943847727, 0.002049428842844314, 0.0020494288428452023, 0.004793821195971937],
      tspan = (0.0, 0.1), timestep_gravity=Trixi.timestep_gravity_erk53_3Sstar!)
  end

  @testset "elixir_eulergravity_jeans_instability.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulergravity_jeans_instability.jl"),
      l2   = [10733.63378538114, 13356.780607423452, 1.6721408750434403e-6, 26834.076821148774],
      linf = [15194.296424901113, 18881.481685044182, 6.809717915492584e-6, 37972.99700513482],
      tspan = (0.0, 0.1))
  end

  @testset "elixir_eulergravity_jeans_instability.jl with RK3S*" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulergravity_jeans_instability.jl"),
      l2   = [10734.598193238024, 13358.217234481384, 1.9105446886104614e-6, 26836.487841241516],
      linf = [15195.661004798487, 18883.512035906537, 7.864715319982429e-6, 37976.408478975296],
      tspan = (0.0, 0.1),
      parameters=ParametersEulerGravity(background_density=1.5e7,
                                        gravitational_constant=6.674e-8,
                                        cfl=2.4,
                                        resid_tol=1.0e-4,
                                        n_iterations_max=1000,
                                        timestep_gravity=timestep_gravity_erk52_3Sstar!))
  end

  @testset "Printing" begin
    # These tests rely on the previous tests above and re-use some of their variables.
    show(stdout, parameters)
    show(stdout, semi)
    show(stdout, semi_euler.boundary_conditions)
    show(stdout, TrivialCallback())
    show(stdout, equations_euler)
  end

  @testset "elixir_eulergravity_sedov_blast_wave.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulergravity_sedov_blast_wave.jl"),
      l2   = [0.046315994852653024, 0.0650818006233669, 0.06508180062336677, 0.4896707211656037],
      linf = [2.3874843337593776, 4.07876384374792, 4.07876384374792, 16.23914384809855],
      tspan = (0.0, 0.05))
  end

  @testset "elixir_eulergravity_sedov_blast_wave.jl with ref-level=8 and no AMR" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulergravity_sedov_blast_wave.jl"),
      l2   = [0.00289222135995042, 0.013724813590853825, 0.013724813590853832, 0.05822904710548214],
      linf = [0.26361780693997594, 1.3908873830688688, 1.3908873830688688, 4.066701303607613],
      tspan = (0.0, 0.005), initial_refinement_level=8, amr_callback=TrivialCallback())
  end
end

# Clean up afterwards: delete Trixi output directory
@test_nowarn rm(outdir, recursive=true)

end #module


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""Local represents directly-attached storage with node affinity (Beta feature)

    IoK8sApiCoreV1LocalVolumeSource(;
        fsType=nothing,
        path=nothing,
    )

    - fsType::String : Filesystem type to mount. It applies only when the Path is a block device. Must be a filesystem type supported by the host operating system. Ex. \"ext4\", \"xfs\", \"ntfs\". The default value is to auto-select a fileystem if unspecified.
    - path::String : The full path to the volume on the node. It can be either a directory or block device (disk, partition, ...).
"""
mutable struct IoK8sApiCoreV1LocalVolumeSource <: SwaggerModel
    fsType::Any # spec type: Union{ Nothing, String } # spec name: fsType
    path::Any # spec type: Union{ Nothing, String } # spec name: path

    function IoK8sApiCoreV1LocalVolumeSource(;fsType=nothing, path=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1LocalVolumeSource, Symbol("fsType"), fsType)
        setfield!(o, Symbol("fsType"), fsType)
        validate_property(IoK8sApiCoreV1LocalVolumeSource, Symbol("path"), path)
        setfield!(o, Symbol("path"), path)
        o
    end
end # type IoK8sApiCoreV1LocalVolumeSource

const _property_map_IoK8sApiCoreV1LocalVolumeSource = Dict{Symbol,Symbol}(Symbol("fsType")=>Symbol("fsType"), Symbol("path")=>Symbol("path"))
const _property_types_IoK8sApiCoreV1LocalVolumeSource = Dict{Symbol,String}(Symbol("fsType")=>"String", Symbol("path")=>"String")
Base.propertynames(::Type{ IoK8sApiCoreV1LocalVolumeSource }) = collect(keys(_property_map_IoK8sApiCoreV1LocalVolumeSource))
Swagger.property_type(::Type{ IoK8sApiCoreV1LocalVolumeSource }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1LocalVolumeSource[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1LocalVolumeSource }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1LocalVolumeSource[property_name]

function check_required(o::IoK8sApiCoreV1LocalVolumeSource)
    (getproperty(o, Symbol("path")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1LocalVolumeSource }, name::Symbol, val)
end


text:
module SemanticReleaseExample

f() = "Hello World"

h() = "Hello, again!"

end


text:
using Distributions: ProbabilityWeights, sample

function execpolicy(model, current_state)
    probs = model(current_state)
    action = sample(1:2, ProbabilityWeights(probs)) # Living that stochastic policy life...
    log_prob = log(probs[action])
    return action, log_prob
end

#=
function getaction(current_state)
    x = relu.(W1 * current_state .+ b1)
    probs = softmax(W2 * x .+ b2)
    action = argmax(probs)
    log_probab = log(probs[action])
    return action, log_probab
end
=#


text:
module Fields
import Luna: Grid, Maths, PhysData, Modes
import Luna.PhysData: wlfreq, ε_0, μ_0
import StaticArrays: SVector
import Cubature: hcubature
import NumericalIntegration: integrate, SimpsonEven
import Random: AbstractRNG, GLOBAL_RNG
import Statistics: mean
import Hankel
import LinearAlgebra: dot, norm
import FFTW
import BlackBoxOptim
import Optim
import DelimitedFiles: readdlm
import HCubature: hquadrature
import DSP: unwrap
import Logging: @warn

abstract type AbstractField end
abstract type TimeField <: AbstractField end

"""
    PulseField(λ0, energy, ϕ, τ0, Itshape)

Represents a temporal pulse with shape defined by `Itshape`.

# Fields
- `λ0::Float64`: the central field wavelength
- `energy::Float64`: the pulse energy
- `power::Float64`: the pulse peak power (**after** applying any spectral phases)
- `ϕ::Vector{Float64}`: spectral phases (CEP, group delay, GDD, TOD, ...)
- `Itshape`: a callable `f(t)` to get the shape of the intensity/power in the time domain
"""
struct PulseField{eT, pT, iT} <: TimeField
    λ0::Float64
    energy::eT
    power::pT
    ϕ::Vector{Float64}
    Itshape::iT
end

function PulseField(;λ0, Itshape, energy=nothing, power=nothing, ϕ=Float64[])
    if !isnothing(power)
        if !isnothing(energy)
            error("only one of `energy` or `power` can be specified")
        end
    elseif isnothing(energy)
        error("one of `energy` or `power` must be specified")
    end
    PulseField(λ0, energy, power, ϕ, Itshape)
end

"""
    GaussField(;λ0, τfwhm, energy, ϕ, m=1)

Construct a (super)Gaussian shaped pulse with intensity/power FWHM `τfwhm`, either
`energy` or peak `power` specified, superGaussian parameter `m=1` and other parameters
as defined for [`PulseField`](@ref).
"""
function GaussField(;λ0, τfwhm, energy=nothing, power=nothing, ϕ=Float64[], m=1)
    if !isnothing(power)
        if !isnothing(energy)
            error("only one of `energy` or `power` can be specified")
        end
    elseif isnothing(energy)
        error("one of `energy` or `power` must be specified")
    end
    PulseField(λ0, energy, power, ϕ, t -> Maths.gauss(t, fwhm=τfwhm, power=2*m))
end

"""
    SechField(;λ0, energy, τw=nothing, τfwhm=nothing, ϕ=0.0, τ0=0.0)

Construct a Sech^2(t/τw) shaped pulse, specifying either the
natural width `τw`, or the intensity/power FWHM `τfwhm`, and either
`energy` or peak `power` specified.
Other parameters are as defined for [`PulseField`](@ref).
"""
function SechField(;λ0, energy=nothing, power=nothing, τw=nothing, τfwhm=nothing,
                    ϕ=Float64[])
    if !isnothing(τfwhm)
        if !isnothing(τw)
            error("only one of `τw` or `τfwhm` can be specified")
        else
            τw = τfwhm/(2*log(1 + sqrt(2)))
        end
    elseif isnothing(τw)
        error("one of `τw` or `τfwhm` must be specified")
    end
    if !isnothing(power)
        if !isnothing(energy)
            error("only one of `energy` or `power` can be specified")
        end
    elseif isnothing(energy)
        error("one of `energy` or `power` must be specified")
    end
    PulseField(λ0, energy, power, ϕ, t -> sech(t/τw)^2)
end

"""
    make_Et(p::PulseField, grid)

Create electric field for `PulseField`, either the field (for `RealGrid`) or
the envelope (for `EnvGrid`)
"""
function make_Et(p::PulseField, grid::Grid.RealGrid)
    t = grid.t
    ω0 = PhysData.wlfreq(p.λ0)
    @. sqrt(p.Itshape(t))*cos(ω0*t)
end

function make_Et(p::PulseField, grid::Grid.EnvGrid)
    t = grid.t
    Δω = PhysData.wlfreq(p.λ0) - grid.ω0
    @. sqrt(p.Itshape(t))*exp(im*(Δω*t))
end

"""
    (p::PulseField)(Eω, grid, energy_t, FT)

Add the field to `Eω` for the provided `grid`, `energy_t` function and Fourier transform `FT`
"""
function (p::PulseField)(grid, FT)
    Et = make_Et(p, grid)
    if length(p.ϕ) >= 1
        Et = FT \ prop_taylor!(FT * Et, grid, p.ϕ, p.λ0)
    end
    if !isnothing(p.energy)
        energy_t = Fields.energyfuncs(grid)[1]
        Et .*= sqrt(p.energy)/sqrt(energy_t(Et))
    else
        Pt = It(Et, grid)
        Et .*= sqrt(p.power)/sqrt(maximum(Pt))
    end

    FT * Et
end

"""
    CWField(Pavg, Aωfunc)

Represents a continuous-wave field with spectral phase/amplitude defined by `Aωfunc`.

# Fields
- `Pavg::Float64`: the average power
- `Aωfunc`: a callable `f(ω)` to get the amplitude/phase of the field in the frequency domain
"""
struct CWField{aT} <: TimeField
    Pavg::Float64
    Aωfunc::aT
end

"""
    CWSech(;λ0, Pavg, Δλ)

Construct a CW field with Sech^2 spectral power density and random phase, with spectral
full-width half-maximim of `Δλ` and other parameters as defined for [`CWField`](@ref).
"""
function CWSech(;λ0, Pavg, Δλ, rng=GLOBAL_RNG)
    ωw = PhysData.ΔλΔω(Δλ, λ0)/(2*log(1 + sqrt(2)))
    ω0 = PhysData.wlfreq(λ0)
    Aωfunc(ω) = let rng=rng, ωw=ωw, ω0=ω0
        sech((ω - ω0)/ωw)*exp(1im*2π*rand(rng))
    end
    CWField(Pavg, Aωfunc)
end

"""
    (c::CWField)(grid, FT)

Get the field for the provided `grid` and Fourier transform `FT`
"""
function (c::CWField)(grid::Grid.EnvGrid, FT)
    Eω = c.Aωfunc.(grid.ω)
    istart = findfirst(isequal(1.0), grid.twin)
    iend = findlast(isequal(1.0), grid.twin)
    shape = abs.(Eω)
    Eω′ = Eω
    while true
        Eω = shape .* exp.(1im .* angle.(Eω))
        Et = FT \ Eω
        Et .*= sqrt(c.Pavg) / sqrt(mean(It(Et, grid)[istart:iend]))
        Et .*= grid.twin
        Eω = FT * Et
        err = norm(Eω .- Eω′)/maximum(abs.(Eω))
        err < 1e-2 && break
        Eω′ = Eω
    end
    Eω
end

struct DataField <: TimeField
    ω::Vector{Float64}
    Iω::Vector{Float64}
    ϕω::Vector{Float64}
    energy::Float64
    ϕ::Vector{Float64}
    λ0::Float64
end

"""
    DataField(ω, Iω, ϕω; energy, ϕ=Float64[], λ0=NaN)

Represents a field with spectral power density `Iω` and spectral phase `ϕω`, sampled on
radial frequency axis `ω`.
"""
DataField(ω, Iω, ϕω; energy, ϕ=Float64[], λ0=NaN) = DataField(ω, Iω, ϕω, energy, ϕ, λ0)

"""
    DataField(ω, Eω; energy, ϕ=Float64[], λ0=NaN)

Create a `DataField` from the complex frequency-domain field `Eω` sampled on radial
frequency grid `ω`.
"""
DataField(ω, Eω; energy, ϕ=Float64[], λ0=NaN) = DataField(ω, abs2.(Eω), unwrap(angle.(Eω)),
                                                          energy, ϕ, λ0)

"""
    DataField(fpath; energy, ϕ=Float64[], λ0=NaN)

Create a `DataField` by loading `ω`, `Iω`, and `ϕω` from the file at `fpath`. The file must
contain 3 columns:

- frequency in Hz
- spectral power density (arbitrary units)
- unwrapped spectral phase
"""
function DataField(fpath; energy, ϕ=Float64[], λ0=NaN)
    dat = readdlm(fpath, ' ')
    DataField(dat[:, 1]*2π, dat[:, 2], dat[:, 3]; energy, ϕ)
end

"""
    (d::DataField)(grid, FT)

Interpolate the `DataField` onto the provided `grid` (note the argument `FT` is unused).
"""
function (d::DataField)(grid::Grid.AbstractGrid, FT)
    if maximum(grid.ω) < maximum(d.ω)
        @warn("Interpolating onto a coarser grid may clip the input spectrum.")
    end
    energy_ω = Fields.energyfuncs(grid)[2]
    ϕg = Maths.BSpline(d.ω, d.ϕω).(grid.ω)
    Ig = Maths.BSpline(d.ω, d.Iω).(grid.ω)
    Ig[Ig .< 0] .= 0
    Ig[.!(minimum(d.ω) .< grid.ω .< maximum(d.ω))] .= 0
    Ig .*= grid.ωwin
    Eω = sqrt.(Ig) .* exp.(1im.*ϕg)
    Eω .*= sqrt(d.energy/energy_ω(Eω))
    τ = length(grid.t) * (grid.t[2] - grid.t[1])/2
    Eω .*= exp.(-1im .* grid.ω .* τ)
    if length(d.ϕ) >= 1
        λ0 = isnan(d.λ0) ? wlfreq(Maths.moment(d.ω, d.Iω)) : d.λ0
        prop_taylor!(Eω, grid, d.ϕ, λ0)
    end
    Eω
end

"""
    PropagatedField(propagator!, field)

A wrapper around a previously defined `TimeField` which applies the mutating propagation
function `propagator!(Eω, grid)` to the field `Eω`.
"""
struct PropagatedField{pT, fT<:TimeField} <: TimeField
    propagator!::pT
    field::fT
end

PropagatedField(::Nothing, field::fT) where fT <:TimeField = field

function (pf::PropagatedField)(grid::Grid.AbstractGrid, FT)
    Eω = pf.field(grid, FT)
    pf.propagator!(Eω, grid)
    Eω
end

"""
    ShotNoise(rng=GLOBAL_RNG)

Creates one photon per mode quantum noise (shot noise) to add to an input field.
If no random number generator `rng` is provided, it defaults to `GLOBAL_RNG`
"""
struct ShotNoise{rT<:AbstractRNG} <: TimeField
    rng::rT
end

function ShotNoise(rng=GLOBAL_RNG)
    ShotNoise(rng)
end

"""
    (s::ShotNoise)(Eω, grid)

Get shotnoise for the provided `grid`. The optional parameter `FT`
is unused and is present for interface compatibility with [`TimeField`](@ref).
"""
function (s::ShotNoise)(grid::Grid.RealGrid, FT=nothing)
    δω = grid.ω[2] - grid.ω[1]
    δt = grid.t[2] - grid.t[1]
    amp = @. sqrt(PhysData.ħ*grid.ω/δω)
    rFFTamp = sqrt(2π)/2δt*amp
    φ = 2π*rand(s.rng, size(grid.ω)...)
    @. rFFTamp * exp(1im*φ)
end

function (s::ShotNoise)(grid::Grid.EnvGrid, FT=nothing)
    δω = grid.ω[2] - grid.ω[1]
    δt = grid.t[2] - grid.t[1]
    amp = zero(grid.ω)
    amp[grid.sidx] = @. sqrt(PhysData.ħ*grid.ω[grid.sidx]/δω)
    FFTamp = sqrt(2π)/δt*amp
    φ = 2π*rand(s.rng, size(grid.ω)...)
    @. FFTamp * exp(1im*φ)
end

"""
    SpatioTemporalField(λ0, energy, ϕ, τ0, Ishape)

Represents a spatiotemporal pulse with shape defined by `Ishape`.

# Fields
- `λ0::Float64`: the central field wavelength
- `energy::Float64`: the pulse energy
- `ϕ::Float64`: the CEO phase
- `τ0::Float64`: the temproal shift from grid time 0
- `Ishape`: a callable `f(t, xs)` to get the shape of the intensity/power in the time-space domain
"""
struct SpatioTemporalField{iT} <: AbstractField
    λ0::Float64
    energy::Float64
    ϕ::Float64
    τ0::Float64
    Ishape::iT
    propz::Float64
end

"Gaussian temporal-spatial field defined radially"
function GaussGauss(t, r::AbstractVector, fwhm, m, w0)
    Maths.gauss.(t, fwhm=fwhm, power=2*m) .* Maths.gauss.(r, w0/2)'
end

"Gaussian temporal-spatial field defined on x-y grid"
function GaussGauss(t, r::AbstractArray{T,3} where T, fwhm, m, w0)
    Maths.gauss.(t, fwhm=fwhm, power=2*m) .* Maths.gauss.(r, w0/2)
end

"""
    GaussGaussField(;λ0, τfwhm, energy, w0, ϕ=0.0, τ0=0.0, m=1)

Construct a (super)Gaussian shaped pulse with intensity/power FWHM `τfwhm`,
superGaussian parameter `m=1` and Gaussian shaped spatial profile with waist `w0`,
propagation distance from the waist of `propz`,
and other parameters as defined for [`TimeField`](@ref).
"""
function GaussGaussField(;λ0, τfwhm, energy, w0, ϕ=0.0, τ0=0.0, m=1, propz=0.0)
    SpatioTemporalField(λ0, energy, ϕ, τ0,
                        (t, xs) -> GaussGauss(t, xs, τfwhm, m, w0),
                        propz)
end

function make_Etr(s::SpatioTemporalField, grid::Grid.RealGrid, spacegrid)
    t = grid.t .- s.τ0
    ω0 = PhysData.wlfreq(s.λ0)
    sqrt.(s.Ishape(t, spacegrid.r)) .* cos.(ω0.*t .+ s.ϕ)
end

function make_Etr(s::SpatioTemporalField, grid::Grid.EnvGrid, spacegrid)
    t = grid.t .- s.τ0
    Δω = PhysData.wlfreq(s.λ0) - grid.ω0
    sqrt.(s.Ishape(t, spacegrid.r)) .* exp.(im .* (s.ϕ .+ Δω.*t))
end

transform(spacegrid::Hankel.QDHT, FT, Etr) = spacegrid * (FT * Etr)

transform(spacegrid::Grid.FreeGrid, FT, Etr) = FT * Etr

"""
    (s::SpatioTemporalField)(grid, spacegrid, FT)

Get the field for the provided `grid`, `spacegrid` function
and Fourier transform `FT`
"""
function (s::SpatioTemporalField)(grid, spacegrid, FT)
    Etr = make_Etr(s, grid, spacegrid)
    energy_t = Fields.energyfuncs(grid, spacegrid)[1]
    Etr .*= sqrt(s.energy)/sqrt(energy_t(Etr))
    Eωk = transform(spacegrid, FT, Etr)
    if s.propz != 0.0
        prop!(Eωk, s.propz, grid, spacegrid)
    end
    Eωk
end

# TODO: this for FreeGrid
function prop!(Eωk, z, grid, q::Hankel.QDHT)
    kzsq = @. (grid.ω/PhysData.c)^2 - (q.k^2)'
    kzsq[kzsq .< 0] .= 0
    kz = sqrt.(kzsq)
    @. Eωk *= exp(-1im * z * (kz - grid.ω/PhysData.c))
end

"""
    gauss_beam(k, ω0; z=0.0, pol=:y)

Gaussian beam field distribution with waist radius `ω0` and wavenumber `k`,
at position `z` from focus. `pol` describes the polarisation direction,
one of `:x` or `:y`.
"""
function gauss_beam(k, ω0; z=0.0, pol=:y)
    let k=k, ω0=ω0, z=z, pol=pol
        function fieldfunc(xs)
            zr = k*ω0^2/2
            ω = ω0*sqrt(1 + (z/zr)^2)
            R1 = z/(z^2 + zr^2)
            ψ = atan(z/zr)
            phase = exp(-1im * (k*z + k*xs[1]^2*R1/2 - ψ))
            E = ω0/ω * exp(-xs[1]^2/ω^2) * phase
            if pol==:x
                return SVector(E, 0.0)
            else
                return SVector(0.0, E)
            end
        end
    end
end

function int2D(field1, field2, lowerlim, upperlim)
    Ifunc(xs) = 0.5*sqrt(ε_0/μ_0)*dot(conj(field1(xs)), field2(xs))*xs[1]
    abs(hcubature(Ifunc, lowerlim, upperlim)[1])
end

function normalised_field(fieldfunc, rmax)
    scale = 1.0/sqrt(int2D(fieldfunc, fieldfunc, (0.0,0.0), (rmax, 2π)))
    return let scale=scale, fieldfunc=fieldfunc
        (xs) -> fieldfunc(xs) .* scale
    end
end

normalised_gauss_beam(k, ω0; pol=:y) = normalised_field(gauss_beam(k, ω0, pol=pol), 6*ω0)

"""
    coupled_field(i, mode, E, fieldfunc; energy, kwargs...)

Create an element of an input field tuple (for use in `Luna.setup`) based on coupling
field `E` into a `mode`. The index `i` species the mode index. The temporal fields are
initialised using `fieldfunc` (e.g. one of `GaussField`, `SechField` etc.) with the
same keyword arguments.
"""
function coupled_field(i, mode, E, fieldfunc; energy, kwargs...)
    ei = energy * Modes.overlap(mode, E)^2
    (mode=i, fields=(fieldfunc(;energy=ei, kwargs...),))
end

"""
    gauss_beam_init(modes, k, ω0, fieldfunc; energy, kwargs...)

Create an input field tuple (for use in `Luna.setup`) based on coupling a focused
Gaussian beam with focused spot size `ω0` and wavenumber `k` into `modes`.
The temporal fields are initialised using `fieldfunc` (e.g. one of `GaussField`,
`SechField` etc.) with the same keyword arguments.

```jldoctest
julia> a = 125e-6;
julia> energy = 1e-3;
julia> λ0 = 800e-9;
julia> modes = (Capillary.MarcatiliMode(a, :He, 1.0, m=1), Capillary.MarcatiliMode(a, :He, 1.0, m=2));
julia> fields = Fields.gauss_beam_init(modes, 2*pi/λ0, a*0.64, Fields.GaussField; λ0=λ0, τfwhm=30e-15, energy=energy);
julia> fields[1].fields[1].energy/energy ≈ 0.98071312
true
julia> fields[2].fields[1].energy/energy ≈ 0.0061826217
true
"""
function gauss_beam_init(modes, k, ω0, fieldfunc; energy, kwargs...)
    gauss = normalised_gauss_beam(k, ω0)
    tuple(collect(coupled_field(i, mode, gauss, fieldfunc; energy=energy, kwargs...) for (i,mode) in enumerate(modes))...)
end

It(Et, grid::Grid.RealGrid) = abs2.(Maths.hilbert(Et))
It(Et, grid::Grid.EnvGrid) = abs2.(Et)

iFT(Eω, grid::Grid.RealGrid) = FFTW.irfft(Eω, length(grid.t), 1)
iFT(Eω, grid::Grid.EnvGrid) = FFTW.ifft(Eω, 1)

"Calculate energy from modal field E(t)"
function energyfuncs(grid::Grid.RealGrid)
    function energy_t(Et)
        return integrate(grid.t, It(Et, grid), SimpsonEven())
    end

    prefac = 2π/(grid.ω[end]^2)
    function energy_ω(Eω)
        prefac*integrate(grid.ω, abs2.(Eω), SimpsonEven())
    end
    return energy_t, energy_ω
end

function energyfuncs(grid::Grid.EnvGrid)
    function energy_t(Et)
        return integrate(grid.t, It(Et, grid), SimpsonEven())
    end

    δω = grid.ω[2] - grid.ω[1]
    Δω = length(grid.ω)*δω
    prefac = 2π*δω/(Δω^2)
    function energy_ω(Eω)
        prefac*sum(abs2.(Eω))
    end
    return energy_t, energy_ω
end

function energyfuncs(grid::Grid.RealGrid, q::Hankel.QDHT)
    function energy_t(Et)
        Eta = Maths.hilbert(Et)
        tintg = integrate(grid.t, abs2.(Eta), SimpsonEven())
        return 2π*PhysData.c*PhysData.ε_0/2 * Hankel.integrateR(tintg, q)
    end

    prefac = 2π*PhysData.c*PhysData.ε_0/2 * 2π/(grid.ω[end]^2)
    function energy_ω(Eω)
        ωintg = integrate(grid.ω, abs2.(Eω), SimpsonEven())
        return prefac*Hankel.integrateK(ωintg, q)
    end
    return energy_t, energy_ω
end

function energyfuncs(grid::Grid.EnvGrid, q::Hankel.QDHT)
    function energy_t(Et)
        tintg = integrate(grid.t, abs2.(Et), SimpsonEven())
        return 2π*PhysData.c*PhysData.ε_0/2 * Hankel.integrateR(tintg, q)
    end

    δω = grid.ω[2] - grid.ω[1]
    Δω = length(grid.ω)*δω
    prefac = 2π*PhysData.c*PhysData.ε_0/2 * 2π*δω/(Δω^2)
    function energy_ω(Eω)
        ωintg = dropdims(sum(abs2.(Eω); dims=1), dims=1)
        return prefac*Hankel.integrateK(ωintg, q)
    end
    return energy_t, energy_ω
end

function energyfuncs(grid::Grid.RealGrid, xygrid::Grid.FreeGrid)
    δx = xygrid.x[2] - xygrid.x[1]
    δy = xygrid.y[2] - xygrid.y[1]
    δt = grid.t[2] - grid.t[1]
    prefac_t = PhysData.c*PhysData.ε_0/2 * δx * δy * δt
    function energy_t(Et)
        Eta = Maths.hilbert(Et)
        return  prefac_t * sum(abs2.(Eta))
    end

    δω = grid.ω[2] - grid.ω[1]
    Δω = grid.ω[end]
    δkx = xygrid.kx[2] - xygrid.kx[1]
    Δkx = length(xygrid.kx)*δkx
    δky = xygrid.ky[2] - xygrid.ky[1]
    Δky = length(xygrid.ky)*δky
    prefac = PhysData.c*PhysData.ε_0/2 * 2π*δω/(Δω^2) * 2π*δkx/(Δkx^2) * 2π*δky/(Δky^2)
    energy_ω(Eω) = prefac * sum(abs2.(Eω))

    return energy_t, energy_ω
end

function energyfuncs(grid::Grid.EnvGrid, xygrid::Grid.FreeGrid)
    δx = xygrid.x[2] - xygrid.x[1]
    δy = xygrid.y[2] - xygrid.y[1]
    δt = grid.t[2] - grid.t[1]
    prefac_t = PhysData.c*PhysData.ε_0/2 * δx * δy * δt
    function energy_t(Et)
        return  prefac_t * sum(abs2.(Et))
    end

    δω = grid.ω[2] - grid.ω[1]
    Δω = length(grid.ω)*δω
    δkx = xygrid.kx[2] - xygrid.kx[1]
    Δkx = length(xygrid.kx)*δkx
    δky = xygrid.ky[2] - xygrid.ky[1]
    Δky = length(xygrid.ky)*δky
    prefac = PhysData.c*PhysData.ε_0/2 * 2π*δω/(Δω^2) * 2π*δkx/(Δkx^2) * 2π*δky/(Δky^2)
    energy_ω(Eω) = prefac * sum(abs2.(Eω))

    return energy_t, energy_ω
end

"""
    prop_taylor!(Eω, grid, ϕs, λ0)
    prop_taylor!(Eω, grid::Grid.AbstractGrid, ϕs, λ0)

Add spectral phase, given as Taylor-expansion coefficients `ϕs` around central wavelength
`λ0`, to the frequency-domain field `Eω`. Sampling axis of `Eω` can be given either as an
`AbstractGrid` or the frequency axis `ω`.
"""
function prop_taylor!(Eω, ω, ϕs, λ0)
    Δω = ω .- wlfreq(λ0)
    ϕ = zeros(length(ω))
    for (n, ϕi) in enumerate(ϕs)
        ϕ .+= Δω.^(n-1)./factorial(n-1) * ϕi
    end
    Eω .*= exp.(-1im.*ϕ)
end

prop_taylor!(Eω, grid::Grid.AbstractGrid, ϕs, λ0) = prop_taylor!(Eω, grid.ω, ϕs, λ0)

"""
    prop_taylor(Eω, grid, ϕs, λ0)
    prop_taylor(Eω, grid::Grid.AbstractGrid, ϕs, λ0)

Return a copy of the frequency-domain field `Eω` with added spectral phase, given as Taylor-expansion coefficients `ϕs` around central wavelength `λ0`.
Sampling axis of `Eω` can be given either as an `AbstractGrid` or the frequency axis `ω`.
"""
prop_taylor(Eω, args...) = prop_taylor!(copy(Eω), args...)

"""
    prop_material!(Eω, ω, material, thickness, λ0=nothing;
                   P=1, T=PhysData.roomtemp, lookup=nothing)

Linearly propagate the frequency-domain field `Eω` through a certain `thickness` of a `material`.
If the central wavelength `λ0` is given, remove the group delay at this wavelength.
Keyword arguments `P` (pressure), `T` (temperature)
and `lookup` (whether to use lookup table instead of Sellmeier expansion).
"""
function prop_material!(Eω, ω, material, thickness, λ0=nothing; kwargs...)
    propagator_material(material; kwargs...)(Eω, ω, thickness, λ0)
end

prop_material!(Eω, grid::Grid.AbstractGrid, args...; kwargs...) = prop_material!(
    Eω, grid.ω, args...; kwargs...)

"""
    prop_material(Eω, ω, material, thickness, λ0=nothing;
                   P=1, T=PhysData.roomtemp, lookup=nothing)

Return a copy of the frequency-domain field `Eω` after linear propagation through a certain
`thickness` of a `material`. If the central wavelength `λ0` is given, remove the group
delay at this wavelength. Keyword arguments `P` (pressure), `T` (temperature)
and `lookup` (whether to use lookup table instead of Sellmeier expansion).
"""
prop_material(Eω, args...; kwargs...) = prop_material!(copy(Eω), args...; kwargs...)

"""
    propagator_material(material; P=1, T=PhysData.roomtemp, lookup=nothing)

Create a function `prop!(Eω, ω, thickness, λ0)` which propagates the field `Eω` through
a certain `thickness` of a `material`. If the central wavelength `λ0` is given, remove the group
delay at this wavelength.  Keyword arguments `P` (pressure), `T` (temperature)
and `lookup` (whether to use lookup table instead of Sellmeier expansion).
"""
function propagator_material(material; P=1, T=PhysData.roomtemp, lookup=nothing)
    n = PhysData.ref_index_fun(material, P, T; lookup=lookup)
    β1 = PhysData.dispersion_func(1, n)
    function prop!(Eω, ω, thickness, λ0=nothing)
        β = ω./PhysData.c .* n.(wlfreq.(ω))
        if !isnothing(λ0)
            β .-= β1(λ0) .* (ω .- wlfreq(λ0))
        end
        β[.!isfinite.(β)] .= 0
        Eω .*= exp.(-1im.*real(β).*thickness)
    end
    prop!
end

"""
    prop_mirror!(Eω, ω, mirror, reflections)

Propagate the field `Eω` linearly by adding a number of `reflections` from the `mirror` type.
"""
function prop_mirror!(Eω, ω, mirror, reflections)
    λ = wlfreq.(ω)
    t = PhysData.lookup_mirror(mirror).(λ) # transfer function
    tn = t.^reflections
    tn[.!isfinite.(tn)] .= 0
    if reflections < 0
        tn = exp.(1im .* angle.(tn))
    end
    Eω .*= tn
end

prop_mirror!(Eω, grid::Grid.AbstractGrid, args...) = prop_mirror!(Eω, grid.ω, args...)

prop_mirror(Eω, args...) = prop_mirror!(copy(Eω), args...)

"""
    prop_mode!(Eω, ω, mode, distance, λ0=nothing)

Propagate the field `Eω` linearly by a certain `distance` in the given `mode`. If the
central wavelength `λ0` is given, remove the group delay at this wavelength. Propagation
includes both dispersion and loss.
"""
function prop_mode!(Eω, ω, mode, distance, λ0=nothing)
    β(z) = ω./PhysData.c .* Modes.neff.(mode, ω; z=z)
    β1(z) = Modes.dispersion(mode, 1, wlfreq(λ0); z=z)
    βint, err = hquadrature(β, 0, abs(distance))
    if !isnothing(λ0)
        β1int, err = hquadrature(β1, 0, abs(distance))
        βint .-= β1int .* (ω .- wlfreq(λ0))
    end
    expφ = exp.(-1im.*sign(distance).*conj(βint))
    expφ[.!isfinite.(expφ)] .= 0
    Eω .*= expφ
end

prop_mode!(Eω, grid::Grid.AbstractGrid, args...) = prop_mode!(Eω, grid.ω, args...)

prop_mode(Eω, args...) = prop_mode!(copy(Eω), args...)

"""
    optcomp_taylor(Eω, grid, λ0; order=2)

Maximise the peak power of the field `Eω` by adding Taylor-expanded spectral phases up to
order `order`.
"""
function optcomp_taylor(Eω::AbstractVecOrMat, grid, λ0; order=2, boundfac=8)
    τ = length(grid.t) * (grid.t[2] - grid.t[1])/2
    EωFTL = abs.(Eω) .* exp.(-1im .* grid.ω .* τ)
    ItFTL = _It(iFT(EωFTL, grid), grid)
    target = 1/maximum(ItFTL)

    Eωnorm = Eω ./ sqrt(maximum(ItFTL))

    function f(disp)
        # disp here is just the dispersion terms (2nd order and higher)
        ϕs = [0, 0, disp...]
        Eωp = prop_taylor(Eωnorm, grid, ϕs, λ0)
        Itp = _It(iFT(Eωp, grid), grid)
        1/maximum(Itp)
    end

    τ0FTL = Maths.fwhm(grid.t, ItFTL)/(2*sqrt(log(2)))
    τ0 = Maths.fwhm(grid.t, _It(iFT(Eω, grid), grid))/(2*sqrt(log(2)))

    ϕ2_0 = τ0FTL*sqrt(τ0^2 - τ0FTL^2) # GDD to stretch Gaussian from FTL to actual duration

    # for Gaussian with pure GDD, sqrt(ϕ2_0) is the FTL duration, so use that as guide
    bounds = boundfac*(sqrt(ϕ2_0) .^(2:order))
    srange = [(-bi, bi) for bi in bounds]
    res = BlackBoxOptim.bboptimize(f; SearchRange=srange,
                                   TraceMode=:silent, TargetFitness=target)
    ϕs = [0, 0, BlackBoxOptim.best_candidate(res)...]
    ϕs, prop_taylor(Eω, grid, ϕs, λ0)
end

function optcomp_taylor(Eω, grid, λ0; order=2)
    out = similar(Eω)
    cidcs = CartesianIndices(size(Eω)[3:end])
    ϕsout = zeros(order+1, size(cidcs)...)
    for ci in cidcs
        ϕsi, Eωi = optcomp_taylor(Eω[:, :, ci], grid, λ0; order=order)
        out[:, :, ci] .= Eωi
        ϕsout[:, ci] .= ϕsi
    end
    ϕsout, out
end

_It(Et::AbstractVector, grid) = It(Et, grid)
_It(Et::AbstractMatrix, grid) = dropdims(sum(It(Et, grid); dims=2); dims=2)

"""
    optcomp_material(Eω, grid, material, λ0; kwargs...)

Maximise the peak power of the field `Eω` by linear propagation through the `material`.
Keyword arguments `kwargs` are the same as for [`prop_material`](@ref).
"""
function optcomp_material(Eω::AbstractVecOrMat, grid, material, λ0,
                          min_thickness, max_thickness; kwargs...)
    τ = length(grid.t) * (grid.t[2] - grid.t[1])/2
    EωFTL = abs.(Eω) .* exp.(-1im .* grid.ω .* τ)
    ItFTL = _It(iFT(EωFTL, grid), grid)
    target = 1/maximum(ItFTL)

    Eωnorm = Eω ./ sqrt(maximum(ItFTL))

    prop! = propagator_material(material; kwargs...)

    function f(d)
        # d is the material insertion
        Eωp = copy(Eωnorm)
        prop!(Eωp, grid.ω, d, λ0)
        Itp = _It(iFT(Eωp, grid), grid)
        1/maximum(Itp)
    end

    # res = BlackBoxOptim.bboptimize(f; SearchRange=[srange],
    #                                TraceMode=:silent, TargetFitness=target)
    # prop_material(Eω, grid, material, BlackBoxOptim.best_candidate(res), λ0; kwargs...)
    res = Optim.optimize(f, min_thickness, max_thickness)
    res.minimizer, prop_material(Eω, grid, material, res.minimizer, λ0; kwargs...)
end

function optcomp_material(Eω, args...; kwargs...)
    out = similar(Eω)
    cidcs = CartesianIndices(size(Eω)[3:end])
    dout = zeros(size(cidcs))
    for ci in cidcs
        di, Eωi = optcomp_material(Eω[:, :, ci], args...; kwargs...)
        out[:, :, ci] .= Eωi
        dout[ci] = di
    end
    dout, out
end

"""
    optfield_cep(Eω, grid)

Find the value of the absolute phase which produces the maximal field strength in the time
domain.
"""
function optfield_cep(Eω::AbstractVector, grid)
    res = Optim.optimize(-π, π) do ϕ
        Et = real(iFT(Eω*exp(1im*ϕ), grid))
        1/maximum(Et)
    end

    res.minimizer, Eω*exp(1im*ϕ)
end

function optfield_cep(Eω::AbstractMatrix, grid; mode=1)
    res = Optim.optimize(-π, π) do ϕ
        Et = real(iFT(Eω[:, mode]*exp(1im*ϕ), grid))
        1/maximum(Et)
    end

    res.minimizer, Eω*exp(1im*res.minimizer)
end

end


text:
using Pkg
pkg"activate ."
using DrWatsonSim
const ds = DrWatsonSim
using Test
using BSON
using DrWatson
using Dates

include("helper_functions.jl")

@testset "DrWatsonSim.jl" begin
    @testset "Metadata" begin
        @eval ds max_lock_retries = 10

        dummy_project() do folder
            @testset "Locking functions" begin
                @test_logs (:info, r"creating") ds.assert_metadata_directory()
                # Check if index file was created
                ds.lock("metadata")
                @test isdir(ds.metadatadir("metadata.lck"))
                ds.lock("foo")
                @test isdir(ds.metadatadir("foo.lck"))
                @test_throws ErrorException ds.lock("foo")
                ds.unlock("foo")
                ds.lock("foo")
                ds.unlock("foo")
                @test_throws ErrorException ds.lock("metadata")
                @test_throws ErrorException ds.lock("metadata")
                ds.unlock("metadata")
                @test_throws ErrorException ds.unlock("metadata")
                ds.lock("metadata")
                ds.unlock("metadata")
                @test ds.semaphore_status("bar") == 0
                ds.semaphore_enter("bar")
                @test isfile(ds.metadatadir("bar.sem"))
                ds.semaphore_enter("bar")
                ds.semaphore_enter("bar")
                @test ds.semaphore_status("bar") == 3
                ds.semaphore_exit("bar")
                @test ds.semaphore_status("bar") == 2
                ds.semaphore_exit("bar")
                @test ds.semaphore_status("bar") == 1
                ds.semaphore_exit("bar")
                @test ds.semaphore_status("bar") == 0
                @test !isfile(ds.metadatadir("bar.sem"))
                @test_throws ErrorException ds.semaphore_exit("bar")
                function sem_test()
                    function blocked_worker(v)
                        ds.lock("foo", wait_for_semaphore="bar")
                        v[1]=1
                        ds.unlock("foo")
                    end
                    v = [0]
                    ds.semaphore_enter("bar")
                    @async blocked_worker(v)
                    @test v[1] == 0
                    ds.semaphore_enter("bar")
                    @test v[1] == 0
                    ds.semaphore_exit("bar")
                    ds.semaphore_exit("bar")
                    yield()
                    @test v[1] == 1
                end
                @sync sem_test()
            end
        end

        dummy_project() do folder
            @testset "Identifer Creation" begin
                ds.assert_metadata_directory()
                @test ds.hash_path(datadir("sims","a.bson")) == hash("data/sims/a.bson")
            end
        end

        dummy_project() do folder
            @testset "Metadata creation" begin
                m = Metadata(datadir("fileA"))
                @test m.path == datadir("fileA")
                @test isfile(ds.metadatadir(ds.hash_path(m.path)|>ds.to_file_name))
                mb = Metadata(datadir("fileB"))
                @test isfile(ds.metadatadir(ds.hash_path(mb.path)|>ds.to_file_name))
                @test m.mtime == 0
                touch(datadir("fileA"))
                @test_logs (:warn, r"changed") Metadata(datadir("fileA"))
                @test_nowarn m = Metadata!(datadir("fileA"))
                @test m.mtime > 0
                A = rand(3,3)
                m["some_data"] = A
                raw_loaded = BSON.load(ds.metadatadir(ds.to_file_name(ds.hash_path(m.path))))
                @test raw_loaded["data"]["some_data"] == A
                rename!(m, datadir("fileC"))
                @test "some_data" in keys(Metadata(datadir("fileC")))
            end
        end

        dummy_project() do folder
            @testset "Metadata Pentest" begin
                ds.assert_metadata_directory()
                @sync for i in 1:500
                    @async begin
                        m = Metadata(datadir("file$i"))
                        s = rand(1:100)
                        m["data"] = rand(s,s)
                        for j in 1:10
                            rename!(m, datadir("file$(i)_$j"))
                        end
                    end
                end
                index = filter(x->endswith(x,".bson"),readdir(ds.metadatadir()))
                @test length(index) == 500
                files = Set{String}()
                for f in index
                    d = BSON.load(ds.metadatadir(f))
                    @test endswith(d["path"],"_10")
                    push!(files,d["path"])
                end
                @test length(files) == 500
            end
        end

    end
    @testset "Simulation" begin
        @eval ds max_lock_retries = 10000

        dummy_project() do folder
            @testset "id generation" begin
                mkdir(datadir("sims2"))
                @test ds.get_next_simulation_id(datadir("sims")) == 1
                @test ds.get_next_simulation_id(datadir("sims2")) == 1
                @test ds.get_next_simulation_id(datadir("sims2")) == 2
                rm(datadir("sims2","1"))
                @test ds.get_next_simulation_id(datadir("sims2")) == 1
                for i in 1:10
                    @test ds.get_next_simulation_id(datadir("sims")) == 1+i
                end
                for i in 1:10
                    rm(datadir("sims","$i"))
                end
                @test ds.get_next_simulation_id(datadir("sims")) == 1
            end
        end

        dummy_project() do folder
            @testset "long running computation" begin
                Pkg.develop(PackageSpec(url=joinpath(@__DIR__,"..")))
                pkg"add BSON"
                pkg"add Dates"
                file = scriptsdir("long_running_script.jl")
                cp(joinpath(@__DIR__, "long_running_script.jl"), file)
                run(`julia $file`)
                for i in 1:4
                    folder = datadir("sims","$i")
                    file = datadir("sims","$i","output.bson")
                    @test isfile(file)
                    result = BSON.load(file)[:result]
                    m = Metadata(folder)
                    p = m["parameters"]
                    @test p[:a]^p[:b] == result
                    @test m["type"] == "Simple Computation"
                    @test m["started at"] < now()
                    m_new = Metadata(joinpath(folder,"newfile"))
                    @test m_new["extra"] == "This should be blocked"
                    @test m["simulation_submit_group"] == ["data/sims/$j" for j in 1:4]
                end
            end
        end

        dummy_project() do folder
            @testset "Rerun simulation" begin
                Pkg.develop(PackageSpec(url=joinpath(@__DIR__,"..")))
                pkg"add BSON"
                pkg"add Dates"
                file = scriptsdir("long_re_running_script.jl")
                cp(joinpath(@__DIR__, "long_re_running_script.jl"), file)
                run(`julia $file`)
                for i in 1:4
                    folder = datadir("sims","$i")
                    fileA = datadir("sims","$i","output.bson")
                    fileB = datadir("sims","$i","output_first_run.bson")
                    resultA = BSON.load(fileA)[:result]
                    resultB = BSON.load(fileB)[:result]
                    @test resultA == resultB
                end
            end
        end
    end

    @testset "Searching" begin
        dummy_project() do folder
            @testset "get by path" begin
                ds.assert_metadata_directory()
                m = Metadata(datadir("sims","1"))
                m["Foo"] = "Bar"
                m2 = get_metadata(datadir("sims","1","111"))
                @test m2["Foo"] == m["Foo"]
                @test get_metadata(datadir("sims","1","111"),include_parents=false) == nothing
                m = Metadata(datadir("sims","2"))
                m["Foo"] = "Baz"
                @test length(get_metadata()) == 2
                @test length(get_metadata("Foo","Baz")) == 1
                @test get_metadata("Foo","Baz")[1].path == m.path
            end
        end
    end

end

max_stars_repo_name: sebastianpech/DrWatsonSim

max_issues_repo_name: sebastianpech/DrWatsonSim

max_forks_repo_name: sebastianpech/DrWatsonSim


text:
struct JobGraph
	scheduler::Scheduler
	sampleStatus::Ref{String}
	fileIOLock::ReentrantLock
	paramIDs::Dict{String,JobID}
	loadSampleID::JobID
	normalizeID::JobID
	setupSimplicesID::JobID
	dimreductionID::JobID
	makeplotID::JobID
	exportSingleID::JobID
	exportMultipleID::JobID
end

struct SampleData
	data::Matrix
	sa::DataFrame
	va::DataFrame
end
SampleData() = SampleData(zeros(0,0),DataFrame(),DataFrame(),"")

struct ReducedSampleData
	F::Factorization
	sa::DataFrame
	va::DataFrame
end

guesslastsampleannot(df::DataFrame) =
	something(findlast(j->!(eltype(df[!,j])<:Union{Real,Missing}), 1:size(df,2)), 1) # a reasonable guess for which column to use as the last sample annotation

getdelim(filepath::String) = lowercase(splitext(filepath)[2])==".csv" ? ',' : '\t'

loadcsv(filepath::String; delim, transpose::Bool=false) =
	CSV.read(filepath, DataFrame; delim=delim, transpose=transpose, ntasks=1) # ntasks=1 might prevent crashes on some systems???

function loadsample(st, input::Dict{String,Any})::DataFrame
	@assert length(input)==2
	filepath = input["filepath"]::String
	rowsAsSamples   = parse(Bool,input["rowsassamples"])
	filepath == :__INVALID__ && return Nothing
	filepath::String
	isempty(filepath) && return Nothing
	@assert isfile(filepath) "Sample file not found: \"$filepath\""
	df = loadcsv(filepath; delim=getdelim(filepath), transpose=!rowsAsSamples)
	@assert size(df,2)>1 "Invalid data set. Must contain at least one sample annotation and one variable."
	@assert guesslastsampleannot(df)<size(df,2) "Invalid data set. Numerical data (variables) must come after sample annotations."
	df
end

# callback function
function showsampleannotnames(df::DataFrame, toGUI)
	indLastSampleAnnot = guesslastsampleannot(df)
	put!(toGUI, :displaysampleannotnames=>(names(df)[1:min(max(indLastSampleAnnot+10,40),end)], indLastSampleAnnot))
end

function normalizesample(st, input::Dict{String,Any})
	@assert length(input)==5
	df = input["dataframe"]
	method             = input["method"]
	lastSampleAnnot    = input["lastsampleannot"]
	logTransform       = parse(Bool,input["logtransform"])
	logTransformOffset = parse(Float64, input["logtransformoffset"])
	@assert method in ("None", "Mean=0", "Mean=0,Std=1")

	nbrSampleAnnots = findfirst(x->string(x)==lastSampleAnnot, names(df))
	@assert nbrSampleAnnots != nothing "Couldn't find sample annotation: \"$lastSampleAnnot\""

	sa = df[:, 1:nbrSampleAnnots]
	va = DataFrame(VariableId=propertynames(df)[nbrSampleAnnots+1:end])
	originalData = Matrix(df[!,nbrSampleAnnots+1:end])'
	@assert eltype(originalData) <: Union{Number,Missing}

	if logTransform
		@assert all(x->x>-logTransformOffset, skipmissing(originalData)) "Data contains negative values, cannot log-transform."
		originalData .= log2.(originalData.+logTransformOffset)
	end

	data = zeros(size(originalData))
	if any(ismissing,originalData)
		# Replace missing values with mean over samples with nonmissing data
		@info "Reconstructing missing values (taking the mean over all nonmissing samples)"
		for i=1:size(data,1)
			m = ismissing.(originalData[i,:])
			data[i,.!m] .= originalData[i,.!m]
			reconstructed = mean(originalData[i,.!m])
			data[i,m] .= isnan(reconstructed) ? 0.0 : reconstructed
		end
	else
		data .= originalData # just copy
	end
	@assert all(!isnan, data) "Input data cannot contain NaNs. Use empty cells for missing values."

	X = data
	if method == "Mean=0,Std=1"
		X = normalizemeanstd!(X)
	elseif method == "Mean=0"
		X = normalizemean!(X)
	end
	SampleData(X,sa,va)
end

function setupsimplices(st, input::Dict{String,Any})
	@assert length(input)==6
	sampleData  = input["sampledata"]
	method      = Symbol(input["method"])
	sampleAnnot = Symbol(input["sampleannot"])
	timeAnnot   = Symbol(input["timeannot"])
	kNN         = parse(Int,input["knearestneighbors"])
	distNN      = parse(Float64, input["distnearestneighbors"])
	@assert method in (:SA,:Time,:NN,:NNSA)

	local G, description
	if method == :SA
		G = groupsimplices(sampleData.sa[!,sampleAnnot])
		description = "GroupBy=$sampleAnnot"
	elseif method == :Time
		eltype(sampleData.sa[!,timeAnnot])<:Number || @warn "Expected time annotation to contain numbers, got $(eltype(sampleData.sa[!,timeAnnot])). Fallback to default sorting."
		G = timeseriessimplices(sampleData.sa[!,timeAnnot], groupby=sampleData.sa[!,sampleAnnot])
		description = "Time=$timeAnnot, GroupBy=$sampleAnnot"
	elseif method == :NN
		G = neighborsimplices(sampleData.data; k=kNN, r=distNN, dim=50)
		description = "Nearest Neighbors k=$kNN, r=$distNN"
	elseif method == :NNSA
		G = neighborsimplices(sampleData.data; k=kNN, r=distNN, dim=50, groupby=sampleData.sa[!,sampleAnnot])
		description = "Nearest Neighbors k=$kNN, r=$distNN, GroupBy=$sampleAnnot"
	end
	G, description
end

function dimreduction(st, input::Dict{String,Any})
	@assert length(input)==3
	sampleData      = input["sampledata"]
	sampleSimplices = input["samplesimplices"][1]
	method          = Symbol(input["method"])

	X = sampleData.data::Matrix{Float64}

	# dim = 3
	dim = min(10, size(X)...)

	if method==:PMA
		F = pma(X, sampleSimplices, nsv=dim)
	elseif method==:PCA
		F = svdbyeigen(X,nsv=dim)
	end
	ReducedSampleData(F, sampleData.sa, sampleData.va)
end

function makeplot(st, input::Dict{String,Any})
	@assert length(input)==17
	reduced            = input["reduced"]::ReducedSampleData
	dimReductionMethod = Symbol(input["dimreductionmethod"])
	sampleAnnot        = Symbol(input["sampleannot"])
	sampleSimplices, sampleSimplexDescription = input["samplesimplices"]
	xaxis              = parse(Int,input["xaxis"])
	yaxis              = parse(Int,input["yaxis"])
	zaxis              = parse(Int,input["zaxis"])
	plotWidth          = parse(Int,input["plotwidth"])
	plotHeight         = parse(Int,input["plotheight"])
	showPoints         = parse(Bool,input["showpoints"])
	showLines          = parse(Bool,input["showlines"])
	showTriangles      = parse(Bool,input["showtriangles"])
	markerSize         = parse(Float64,input["markersize"])
	lineWidth          = parse(Float64,input["linewidth"])
	triangleOpacity    = parse(Float64,input["triangleopacity"])
	colorByMethod      = Symbol(input["colorby"])
	colorAnnot         = Symbol(input["colorannot"])
	# shapeByMethod      = Symbol(input["shapeby"])
	# shapeAnnot         = Symbol(input["shapeannot"])

	title = "$dimReductionMethod - $sampleSimplexDescription"

	colorByMethod == :Auto && (colorAnnot=sampleAnnot)
	colorBy = colorByMethod == :None ? repeat([""],size(reduced.sa,1)) : reduced.sa[!,colorAnnot]
	colorDict = (eltype(colorBy) <: Real) ? nothing : colordict(colorBy)

	shapeBy, shapeDict = nothing, nothing
	# if shapeByMethod == :Custom
	# 	shapeBy = reduced.sa[!,shapeAnnot]
	# 	shapeDict = shapedict(shapeBy)
	# end

	# TODO: handle missing values in sample annotations?

	dims = [xaxis, yaxis, zaxis]
	plotArgs = plotsimplices(reduced.F.V[:,dims],sampleSimplices,colorBy,colorDict, title=title,
	                         drawPoints=showPoints, drawLines=showLines, drawTriangles=showTriangles,
	                         opacity=triangleOpacity, markerSize=markerSize, lineWidth=lineWidth,
	                         shapeBy=shapeBy, shapeDict=shapeDict,
	                         xLabel=string("PMA",xaxis), yLabel=string("PMA",yaxis), zLabel=string("PMA",zaxis))
	plotArgs, plotWidth, plotHeight
end

showplot(data, toGUI::Channel) = put!(toGUI, :displayplot=>data)

function exportsingle(st, input::Dict{String,Any})
	@assert length(input)==5
	reduced  = input["reduced"]::ReducedSampleData
	mode     = Symbol(input["mode"])
	filepath = input["filepath"]
	sortMode = Symbol(input["sort"])
	dim      = parse(Int,input["dim"])
	@assert mode in (:Variables, :Samples)
	@assert sortMode in (:Abs, :Descending, :Ascending, :Original)

	colName = Symbol("PMA",dim)

	if mode == :Variables
		df = copy(reduced.va)
		df[!,colName] = reduced.F.U[:,dim]
	else
		df = copy(reduced.sa[!,1:1])
		df[!,colName] = reduced.F.V[:,dim]
	end

	if sortMode==:Abs
		sort!(df, colName, by=abs, rev=true)
	elseif sortMode==:Descending
		sort!(df, colName, rev=true)
	elseif sortMode==:Ascending
		sort!(df, colName)
	end

	CSV.write(filepath, df, delim=getdelim(filepath))
end

function exportmultiple(st, input::Dict{String,Any})
	@assert length(input)==4
	reduced  = input["reduced"]::ReducedSampleData
	mode     = Symbol(input["mode"])
	filepath = input["filepath"]
	dim      = parse(Int,input["dim"])
	@assert mode in (:Variables, :Samples)

	if mode == :Variables
		df = copy(reduced.va)
		PMAs = reduced.F.U
	else
		df = copy(reduced.sa[!,1:1])
		PMAs = reduced.F.V
	end

	for d in 1:dim
		df[!,Symbol("PMA",d)] = PMAs[:,d]
	end

	CSV.write(filepath, df, delim=getdelim(filepath))
end

function samplestatus(jg::JobGraph)
	status = jobstatus(jg.scheduler, jg.loadSampleID)
	status==:done && return "Sample loaded."
	status in (:waiting,:running) && return "Loading sample."
	"Please load sample."
end

function setsamplestatus(jg::JobGraph, toGUI::Channel)
	sampleStatus = samplestatus(jg)
	sampleStatus!=jg.sampleStatus[] && put!(toGUI, :samplestatus=>sampleStatus)
	jg.sampleStatus[] = sampleStatus
end

function JobGraph(;verbose=false)
	scheduler = Scheduler(verbose=verbose)
	# scheduler = Scheduler(threaded=false, verbose=verbose) # For DEBUG

	# Data Nodes (i.e. parameters chosen in the GUI)
	paramIDs = Dict{String,JobID}()

	loadSampleID = createjob!(loadsample, scheduler,
	                          getparamjobid(scheduler,paramIDs,"samplefilepath")=>"filepath",
	                          getparamjobid(scheduler,paramIDs,"loadrowsassamples")=>"rowsassamples")

	normalizeID = createjob!(normalizesample, scheduler,
	                         loadSampleID=>"dataframe",
	                         getparamjobid(scheduler,paramIDs,"lastsampleannot")=>"lastsampleannot",
	                         getparamjobid(scheduler,paramIDs,"normalizemethod")=>"method",
	                         getparamjobid(scheduler,paramIDs,"logtransform")=>"logtransform",
	                         getparamjobid(scheduler,paramIDs,"logtransformoffset")=>"logtransformoffset")

	setupSimplicesID = createjob!(setupsimplices, scheduler,
	                              normalizeID=>"sampledata",
	                              getparamjobid(scheduler,paramIDs,"samplesimplexmethod")=>"method",
	                              getparamjobid(scheduler,paramIDs,"sampleannot")=>"sampleannot",
	                              getparamjobid(scheduler,paramIDs,"timeannot")=>"timeannot",
	                              getparamjobid(scheduler,paramIDs,"knearestneighbors")=>"knearestneighbors",
	                              getparamjobid(scheduler,paramIDs,"distnearestneighbors")=>"distnearestneighbors")

	dimreductionID = createjob!(dimreduction, scheduler,
	                            normalizeID=>"sampledata",
	                            setupSimplicesID=>"samplesimplices",
	                            getparamjobid(scheduler,paramIDs,"dimreductionmethod")=>"method")

	makeplotID = createjob!(makeplot, scheduler,
	                        dimreductionID=>"reduced",
	                        getparamjobid(scheduler,paramIDs,"dimreductionmethod")=>"dimreductionmethod",
	                        getparamjobid(scheduler,paramIDs,"sampleannot")=>"sampleannot",
	                        setupSimplicesID=>"samplesimplices",
	                        getparamjobid(scheduler,paramIDs,"xaxis")=>"xaxis",
	                        getparamjobid(scheduler,paramIDs,"yaxis")=>"yaxis",
	                        getparamjobid(scheduler,paramIDs,"zaxis")=>"zaxis",
	                        getparamjobid(scheduler,paramIDs,"plotwidth")=>"plotwidth",
	                        getparamjobid(scheduler,paramIDs,"plotheight")=>"plotheight",
	                        getparamjobid(scheduler,paramIDs,"showpoints")=>"showpoints",
	                        getparamjobid(scheduler,paramIDs,"showlines")=>"showlines",
	                        getparamjobid(scheduler,paramIDs,"showtriangles")=>"showtriangles",
	                        getparamjobid(scheduler,paramIDs,"markersize")=>"markersize",
	                        getparamjobid(scheduler,paramIDs,"linewidth")=>"linewidth",
	                        getparamjobid(scheduler,paramIDs,"triangleopacity")=>"triangleopacity",
	                        getparamjobid(scheduler,paramIDs,"colorby")=>"colorby",
	                        getparamjobid(scheduler,paramIDs,"colorannot")=>"colorannot")
	                        # getparamjobid(scheduler,paramIDs,"shapeby")=>"shapeby",
	                        # getparamjobid(scheduler,paramIDs,"shapeannot")=>"shapeannot")

	exportSingleID = createjob!(exportsingle, scheduler,
	                            dimreductionID=>"reduced",
	                            getparamjobid(scheduler,paramIDs,"exportmode")=>"mode",
	                            getparamjobid(scheduler,paramIDs,"exportsinglepath")=>"filepath",
	                            getparamjobid(scheduler,paramIDs,"exportsingledim")=>"dim",
	                            getparamjobid(scheduler,paramIDs,"exportsinglesort")=>"sort")

	exportMultipleID = createjob!(exportmultiple, scheduler,
	                              dimreductionID=>"reduced",
	                              getparamjobid(scheduler,paramIDs,"exportmode")=>"mode",
	                              getparamjobid(scheduler,paramIDs,"exportmultiplepath")=>"filepath",
	                              getparamjobid(scheduler,paramIDs,"exportmultipledim")=>"dim")

	JobGraph(scheduler,
	         Ref(""),
	         ReentrantLock(),
	         paramIDs,
	         loadSampleID,
	         normalizeID,
	         setupSimplicesID,
	         dimreductionID,
	         makeplotID,
	         exportSingleID,
	         exportMultipleID)
end

getparamjobid(s::Scheduler, paramIDs::Dict{String,JobID}, name::String, create::Bool=true) = create ? get!(paramIDs,name,createjob!(s, :__INVALID__, name=name)) : paramIDs[name]
getparamjobid(jg::JobGraph, name::String, args...) = getparamjobid(jg.scheduler, jg.paramIDs, name, args...)
function setparam(jg::JobGraph, name::String, value)
	if haskey(jg.paramIDs, name)
		setresult!(jg.scheduler, getparamjobid(jg,name), value)
	else
		@warn "Unknown variable name: $name"
	end
end

function process_step(jg::JobGraph, fromGUI::Channel, toGUI::Channel, lastSchedulerTime::Ref{UInt64}; verbosityLevel)
	scheduler = jg.scheduler

	verbosityLevel>=3 && @info "[Processing] tick"
	timeNow = time_ns()
	if isready(fromGUI)
		try
			msg = take!(fromGUI)
			msgName = msg.first
			msgArgs = msg.second
			verbosityLevel>=2 && @info "[Processing] Got message: $msgName $msgArgs"

			if msgName == :exit
				return false
			elseif msgName == :cancel
				verbosityLevel>=2 && @info "[Processing] Cancelling all future events."
				cancelall!(scheduler)
			elseif msgName == :setvalue
				setparam(jg, msgArgs[1], msgArgs[2])
			elseif msgName == :loadsample
				schedule!(x->showsampleannotnames(x,toGUI), scheduler, jg.loadSampleID)
			elseif msgName == :normalize
				schedule!(scheduler, jg.normalizeID)
			elseif msgName == :dimreduction
				schedule!(scheduler, jg.dimreductionID)
			elseif msgName == :showplot
				schedule!(x->showplot(x,toGUI), scheduler, jg.makeplotID)
			elseif msgName == :exportsingle
				schedule!(scheduler, jg.exportSingleID)
			elseif msgName == :exportmultiple
				schedule!(scheduler, jg.exportMultipleID)
			else
				@warn "Unknown message type: $(msgName)"
			end
		catch e
			@warn "[Processing] Error processing GUI message."
			showerror(stdout, e, catch_backtrace())
		end
		setsamplestatus(jg, toGUI)
	elseif wantstorun(scheduler) || (isactive(scheduler) && (timeNow-lastSchedulerTime[])/1e9 > 5.0)
		lastSchedulerTime[] = timeNow
		try
			step!(scheduler)
			status = statusstring(scheduler)
			!isempty(status) && verbosityLevel>=2 && @info "Job status: $status"
		catch e
			@warn "[Processing] Error processing event: $e"
			#showerror(stdout, e, catch_backtrace())
		end
		setsamplestatus(jg, toGUI)
	else
		sleep(0.05)#yield()
	end
	true
end

function process_thread(jg::JobGraph, fromGUI::Channel, toGUI::Channel; verbosityLevel)
	lastSchedulerTime = Ref{UInt64}(0)
	try
		while process_step(jg, fromGUI, toGUI, lastSchedulerTime; verbosityLevel=verbosityLevel)
		end
	catch e
		@warn "[Processing] Fatal error."
		showerror(stdout, e, catch_backtrace())
	end
	verbosityLevel>=2 && @info "[Processing] Exiting thread."
	put!(toGUI, :exited=>nothing)
end

"""
	pmaapp()

Start the Principal Moment Analysis App.
"""
function pmaapp(; verbosityLevel=1, return_job_graph=false)
	# This is the GUI thread

	jg = JobGraph(verbose=verbosityLevel>=2)

	verbosityLevel>=1 && @info "[PMAGUI] Using $(Threads.nthreads()) of $(Sys.CPU_THREADS) available threads."
	Threads.nthreads() == 1 && verbosityLevel>=1 && @warn "[PMAGUI] Threading not enabled, please set the environment variable JULIA_NUM_THREADS to the desired number of threads."

	# init
	fromGUI = Channel{Pair{Symbol,Vector}}(Inf)
	toGUI   = Channel{Pair{Symbol,Any}}(Inf)

	# start processing thread
	processingThreadRunning = true
	Threads.@spawn process_thread(jg, fromGUI, toGUI; verbosityLevel=verbosityLevel)

	# setup gui
	w = Window(Dict(:width=>512,:height=>768))

	# event listeners
	handle(w, "msg") do args
		msgName = Symbol(args[1])
		msgArgs = args[2:end]
		verbosityLevel>=2 && @info "[GUI] sending message: $msgName $(join(msgArgs,", "))"
		processingThreadRunning && put!(fromGUI, msgName=>msgArgs)
	end

	doc = read(joinpath(@__DIR__,"content.html"),String)
	body!(w,doc,async=false)
	js(w, js"init()")

	while isopen(w.content.sock) # is there a better way to check if the window is still open?
		verbosityLevel>=3 && @info "[GUI] tick"

		if isready(toGUI)
			msg = take!(toGUI)
			msgName = msg.first
			msgArgs = msg.second
			verbosityLevel>=2 && @info "[GUI] got message: $msgName"

			if msgName == :samplestatus
				js(w, js"""setSampleStatus($msgArgs)""")
			elseif msgName == :displayplot
				plotArgs,plotWidth,plotHeight = msgArgs
				pl = plot(plotArgs...)
				relayout!(pl, Dict(:margin=>attr(l=1, r=1, b=1, t=50))) # by doing this here - we seem to avoid a bug that causes axis ticks to disappear from subsequent Volcano(!) plots.

				display(pl)
				# Slight hack to set window size properly
				plSize = size(pl)
				size(pl.window, floor(Int,plotWidth*1.05), floor(Int,plotHeight*1.05))
			elseif msgName == :displaysampleannotnames
				js(w, js"""setSampleAnnotNames($(msgArgs[1]), $(msgArgs[2]-1))""")
			elseif msgName == :exited
				processingThreadRunning = false
			end
		end
		sleep(0.05)#yield() # Allow GUI to run
	end

	verbosityLevel>=2 && @info "[GUI] Waiting for scheduler thread to finish."
	processingThreadRunning && put!(fromGUI, :exit=>[])
	# wait until all threads have exited
	while processingThreadRunning
		msg = take!(toGUI)
		msgName = msg.first
		msgArgs = msg.second
		verbosityLevel>=2 && @info "[GUI] got message: $msgName"
		msgName == :exited && (processingThreadRunning = false)
		sleep(0.05)
	end
	verbosityLevel>=2 && @info "[GUI] Scheduler thread finished."

	return_job_graph ? jg : nothing
end


text:
# ---
# title: 64. Minimum Path Sum
# id: problem64
# author: zhwang
# date: 2022-01-17
# difficulty: Medium
# categories: Array, Dynamic Programming
# link: <https://leetcode.com/problems/minimum-path-sum/description/>
# hidden: true
# ---
#
# Given a `m x n` `grid` filled with non-negative numbers, find a path from top
# left to bottom right, which minimizes the sum of all numbers along its path.
#
# **Note:** You can only move either down or right at any point in time.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/11/05/minpath.jpg)
#
#
#
#     Input: grid = [[1,3,1],[1,5,1],[4,2,1]]
#     Output: 7
#     Explanation: Because the path 1 -> 3 -> 1 -> 1 -> 1 minimizes the sum.
#
#
# **Example 2:**
#
#
#
#     Input: grid = [[1,2,3],[4,5,6]]
#     Output: 12
#
#
#
#
# **Constraints:**
#
#   * `m == grid.length`
#   * `n == grid[i].length`
#   * `1 <= m, n <= 200`
#   * `0 <= grid[i][j] <= 100`
#
#
## @lc code=start
using LeetCode

function min_path_sum!(grid::Vector{Vector{Int}})::Int
    m, n = length(grid), length(grid[1])
    for i in 2:n
        grid[1][i] += grid[1][i - 1]
    end
    for i in 2:m
        grid[i][1] += grid[i - 1][1]
    end
    for i in 2:m
        for j in 2:n
            grid[i][j] += min(grid[i - 1][j], grid[i][j - 1])
        end
    end
    return grid[end][end]
end
min_path_sum(grid::Vector{Vector{Int}}) = min_path_sum!(copy(grid))
## @lc code=end


text:
module blFeats

    import LightGraphs
    import JLD2
    import FileIO
    import YAML
    import DataFrames
    const df = DataFrames
    const lg = LightGraphs
    include("../offline-laptop/mykavosh/kavosh.jl")

    macro p(s)
        quote
            $(esc(parse("fetch.([@spawn " * "$s"[2:end] * ")")))
        end
    end

    # Todo
    # - sensible names / categories
    # - sensible directory structure
    # - sensible way of loading / generating from source files
    # - convert graphholder etc. to new system
    # - datamat should accept a list of features
    # - should parallelise only if feature doesn't exist?
    # - make all keys case insensitive
    # - should probably use "get(dict,key,default)" instead of accessing dict directly

    const CACHE_DIR = "/home/olie/.dropbox-work/julia_graph_cache/"

    type Feature
        name::String
        padfunction::Function # one arg, graph und meta[]
        getfunction::Function # with one argument, graph und meta
    end

    FEATURE_DICT = Dict()

    type GraphAndMeta
        graph::LightGraphs.SimpleGraphs.AbstractSimpleGraph
        meta::Dict{String,Any}
        feature::Dict{String,Any}
    end

    function GraphAndMeta(graph; name=randstring(32))
        meta = Dict(
            "name" => name,
            "location" => CACHE_DIR,
            "category" => "Uncategorised",
        )
        return GraphAndMeta(graph,meta,Dict())
    end

    function konect_loader(directory)
        meta_file = split(readstring(`fish -c "ls $(directory)/**meta.*"`))[1]
        graph_file = split(readstring(`fish -c "ls $(directory)/**out.*"`))[1]
        meta = Dict()
        # Handy script for tidying up metadata:
        # sed 's/:/-/2' $(filename)
        # should probably do it in place here rather than messing with the source files

        try
            meta = YAML.load_file(meta_file)
            if haskey(meta,"Name") meta["name"] = meta["Name"] end
        catch(error)
            if isa(error,YAML.ParserError) || isa(error,YAML.ScannerError)
                meta["name"] = split(graph_file,"/")[end]
            else
                println(typeof(error))
                throw(error)
            end
        end
        meta["location"]=CACHE_DIR

        # Should change type of graph we make based on this
        meta["type"] = join(split(strip(readstring(`head -n1 $(graph_file)`)), " ")[2:end], " ")

        cachedfile = meta["location"]*meta["name"]*".jld2"
        try
            candidate = loadGraphAndMeta(cachedfile)
            merge!(candidate.meta,meta)
            return candidate
        catch(e)
            if isa(e,SystemError)
                print(e)
                graph = graph_from_konect(graph_file)
                return GraphAndMeta(graph,meta,Dict())
            else
                throw(e)
            end
        end
    end

    function graph_from_konect(filename)
        edgeframe= df.readtable(filename,separator=' ',allowcomments=true,commentmark='%')
        edgelist = Array(edgeframe[:,1:2])
        numnodes = maximum(Array(edgelist[:,1:2])[:])
        a = lg.Graph(numnodes)
        for i in 1:size(edgelist)[1]
            edge = tuple(edgelist[i,1:2]...)
            lg.add_edge!(a,edge)
        end
        return a
    end

    function thomps_loader(filename,directed=true)
        meta = Dict()
        meta["location"]=CACHE_DIR
        meta["name"] = split(split(filename,'/')[end],"txt.")[1]
        meta["name"] = directed ? meta["name"] : meta["name"]*"_undirected"
        meta["category"] = "Animal"
        meta["cite"] = "Thompson & Townsend"
        meta["type"] = "sym unweighted"
        cachedfile = meta["location"]*meta["name"]*".jld2"
        try
            candidate = loadGraphAndMeta(cachedfile)
            merge!(candidate.meta,meta)
            return candidate
        catch(e)
            if isa(e,SystemError)
                print(e)
                graph = graph_from_thomps(filename)
                graph = directed ? graph : lg.Graph(graph)
                return GraphAndMeta(graph,meta,Dict())
            else
                throw(e)
            end
        end
    end

    function graph_from_thomps(filename)
        raw = df.readtable(filename,separator='\t')
        adj = Matrix(raw[:,2:end])
        #adj = adj'*adj # not a real undirected graph, want binary or with transpose
        lg.DiGraph(adj)
    end

    #konect_graphs = []
    #konect_strings = split(readstring(`fish -c "ls /home/olie/Dropbox_work/big_data/data_oct16/konect/**/out.* "`))

    function generategraph(graphfunc,args)
        graph = GraphAndMeta(graphfunc(args...),name=string(graphfunc)*join(string.(args),"_")*"_"*randstring(32))
        graph.meta["function"] = graphfunc
        graph.meta["args"] = args
        graph.meta["category"] = "Generated"
        return graph
    end

    function getfeature!(gm,feature::String;)
        if haskey(gm.feature,feature) return gm.feature[feature] end
        gm.feature[feature] = FEATURE_DICT[feature].getfunction(gm)
        saveGraphAndMeta(gm)
        return gm.feature[feature]
    end

    function saveGraphAndMeta(gm;location="")
        location = location == "" ? gm.meta["location"]*gm.meta["name"]*".jld2" : location
        FileIO.save(location,gm.meta["name"],gm)
    end

    function loadGraphAndMeta(location)
        return [v for v in values(FileIO.load(location))][1]
    end

    function loadbyname(name)
        location = CACHE_DIR*name
        return loadGraphAndMeta(location)
    end

    type MotifFeature
        deets::Dict{Array{UInt64,1},Float64}
    end

    type AdjacencyFeature
        deets::SparseMatrixCSC{Int64,Int64}
    end

    function array_from_motif_dicts(motif_dicts)
        function mapping_from_motif_count_array(motifs)
            karr = []
            for ks in keys.(motifs)
                karr = vcat(karr, collect(ks))
            end
            #nauty.label_to_adj.(collect(Set(karr)),4)
            dmap = Dict()
            [dmap[j]=i for (i,j) in enumerate(collect(Set(karr)))]
            return dmap
        end

        dmap = mapping_from_motif_count_array(motif_dicts);
        result = Array{Float64,2}(length(dmap),0)
        for example in motif_dicts
            temparr = zeros(Float64,length(dmap),1)
            for (key, value) in example
                temparr[dmap[key]] = value
            end
            result = hcat(result,temparr)
        end
        rmap = Dict()
        for (k,v) in dmap
            rmap[v] = k
        end
        (result',dmap,rmap)
    end

    function getdatamat(graphs,feature::String)
        return FEATURE_DICT[feature].padfunction(graphs)
    end

    function loadall(filter="*")
        loadGraphAndMeta.(split(readstring(`fish -c "ls \"$(CACHE_DIR)\"$(filter)"`),"\n")[1:end-1])
    end

    function dictcat(dicts)
        ansdict = Dict()
        for dict in dicts
            for (k,v) in dict
                ansdict[k] = push!(get(ansdict,k,[]),v)
            end
        end
        ansdict
    end

    function dictstat(dicts)
        ansdict = Dict()
        catdict = dictcat(dicts)
        dictdot(catdict,v->(mean(v),std(Float64.(v))))
    end

    function zscore(x,μ,σ)
        (x-μ)/σ
    end

    function dictdot!(dict,func)
        for (k,v) in dict
            dict[k] = func(v)
        end
        dict
    end

    function dictdot(dict,func)
        ansdict = Dict()
        for (k,v) in dict
            ansdict[k] = func(v)
        end
        ansdict
    end

    function zmotif(motifdict,normdict)
        ansdict = Dict()
        for (k,v) in motifdict
            # rule of succession, sort of:
            # if motif does not show up at all, assume 1% prevelance
            # with std of 1pp.
            # really we need number of motifs to do this properly
            meanandstd = get(normdict,k,(0.01,0.01))
            ansdict[k] = zscore(v,meanandstd[1],meanandstd[2])
        end
        ansdict
    end

    function getclosedense(g,graphs)
        tol = 0.0001
        firsttry = filter(x-> isapprox(lg.density(x.graph),lg.density(g.graph),atol=tol), graphs)
        while length(firsttry) < 10
            tol *= 2
            firsttry = filter(x-> isapprox(lg.density(x.graph),lg.density(g.graph),atol=tol), graphs)
        end
        firsttry
    end

    function normmotifs(g,graphs)
        motifs = dictstat([blFeats.getfeature!(x,"motif_4") for x in getclosedense(g,graphs)])
        zmotif(blFeats.getfeature!(g,"motif_4"),motifs)
    end

    benchgraphs = loadall("*erdos*")

    function normmotifs(g)
        normmotifs(g, benchgraphs)
    end

    append_to_feature_dict!(feat::Feature) = FEATURE_DICT[feat.name] = feat

    append_to_feature_dict!(Feature(
            "motif_3",
            (graphs)->array_from_motif_dicts(
                [getfeature!(graph,"motif_3") for graph in graphs]),
            (g)->kavosh.getsubgraphs(g.graph,3))
    )

    append_to_feature_dict!(Feature(
            "motif_4",
            (graphs)->array_from_motif_dicts(
                [getfeature!(graph,"motif_4") for graph in graphs]),
            (g)->kavosh.getsubgraphs(g.graph,4))
    )

    append_to_feature_dict!(Feature(
            "motif_4_norm",
            (graphs)->array_from_motif_dicts(
                [getfeature!(graph,"motif_4_norm") for graph in graphs]),
            normmotifs)
    )

    append_to_feature_dict!(Feature(
            "motif_5",
            (graphs)->array_from_motif_dicts(
                [getfeature!(graph,"motif_5") for graph in graphs]),
            (g)->kavosh.getsubgraphs(g.graph,5))
    )

    append_to_feature_dict!(Feature(
            "motif_6",
            (graphs)->array_from_motif_dicts(
                [getfeature!(graph,"motif_6") for graph in graphs]),
            (g)->kavosh.getsubgraphs(g.graph,6))
    )
    # Todo: graph mixing matrices / sorted degree matrices
end


text:
function solve(
    i::MinimumBudget{ElementaryPathInstance{T, Maximise}, T},
    ::DynamicProgramming,
) where {T}
    return solve(i, BellmanFordAlgorithm())
end

function solve(
    i::MinimumBudget{ElementaryPathInstance{T, Maximise}, T},
    ::BellmanFordAlgorithm,
) where {T}
    V = Dict{Tuple{T, Int}, Float64}()
    S = Dict{Tuple{T, Int}, Vector{Edge{T}}}()

    # Initialise. For β = 0, this is exactly Bellman-Ford algorithm.
    # Otherwise, use the same initialisation as Bellman-Ford.
    β0 = solve(i.instance, BellmanFordAlgorithm())
    for v in vertices(i.instance.graph)
        S[v, 0] = β0.solutions[v]
        V[v, 0] =
            length(S[v, 0]) == 0 ? 0 :
            sum(i.instance.rewards[e] for e in S[v, 0])
    end

    for β in 1:(i.min_budget)
        for v in vertices(i.instance.graph)
            V[v, β] = -Inf
            S[v, β] = Edge{T}[]
        end

        V[i.instance.src, β] = 0.0
    end

    # Dynamic part.
    for β in 1:(i.min_budget)
        # Loop needed when at least a weight is equal to zero. TODO: remove it when all weights are nonzero?
        while true
            changed = false

            for v in vertices(i.instance.graph)
                for w in inneighbors(i.instance.graph, v)
                    # Compute the remaining part of the budget still to use.
                    remaining_budget = max(0, β - i.weights[Edge(w, v)])

                    # If the explored subproblem has no solution, skip it.
                    if V[w, remaining_budget] == -Inf
                        continue
                    end

                    # If using the solution to the currently explored subproblem would
                    # lead to a cycle, skip it.
                    if any(src(e) == v for e in S[w, remaining_budget])
                        continue
                    end

                    # Compute the amount of budget already used by the solution to the currently explored subproblem.
                    used_budget = 0
                    if length(S[w, remaining_budget]) > 0
                        used_budget =
                            sum(i.weights[e] for e in S[w, remaining_budget])
                    end
                    if used_budget < remaining_budget
                        continue
                    end

                    # Compute the maximum: is passing through w advantageous?
                    if V[w, remaining_budget] + i.instance.rewards[Edge(w, v)] >
                       V[v, β] && used_budget >= remaining_budget
                        changed = true

                        V[v, β] =
                            V[w, remaining_budget] +
                            i.instance.rewards[Edge(w, v)]
                        S[v, β] = vcat(S[w, remaining_budget], Edge(w, v))
                    end
                end
            end

            if !changed
                break
            end
        end
    end

    S_refined = Dict(β => S[i.instance.dst, β] for β in 0:(i.min_budget))
    return BudgetedElementaryPathSolution(
        i,
        S[i.instance.dst, i.min_budget],
        S_refined,
    )
end