text:
using JuLIPMaterials
using Base.Test
using JuLIP, JuLIP.Potentials
using PyPlot
MST = JuLIPMaterials
BCC = MST.BCC

# data = Pkg.dir("JuLIP") * "/data/"
# eam_W = JuLIP.Potentials.FinnisSinclair(data*"W-pair-Wang-2014.plt", data*"W-e-dens-Wang-2014.plt")

morse = Morse(A = 3.5, r0 = rnn("W")) * C2Shift(2.3 * rnn("W"))
at = bulk("W")
variablecell!(at)
set_calculator!(at, morse)

s = "W"
at1 = BCC.screw_111.(s, 70.1, bsign = 1, soln =:anisotropic, calc=morse)
x, y, z1 = xyz(at1)
scatter(x,y,c=z1)

set_calculator!(at1,morse)
f = norm.(forces(at1))
R = norm.(positions(at1))
loglog(R,f,".")
loglog(R,R.^-2,".")
loglog(R,R.^-1,".")

cell(at1)
BCC.lattice_constants_111("W", calc=morse)


text:
module OSQP

export OSQPMathProgBaseInterface
using SparseArrays
using LinearAlgebra

if isfile(joinpath(dirname(@__FILE__), "..", "deps", "deps.jl"))
    include("../deps/deps.jl")
else
    error("OSQP not properly installed. Please run Pkg.build(\"OSQP\")")
end

function __init__()
    # Get version
    ver_array = split(version(), ".")
    ver_string = string(ver_array[1], ".", ver_array[2], ".", ver_array[3])  # Get string without dev vers
    vnum = VersionNumber(ver_string)

    depsdir = realpath(joinpath(dirname(@__FILE__), "..", "deps"))
    if (vnum.major != 0 && vnum.minor != 6)
        error("Current OSQP version installed is $(osqp_version()), but we require version 0.6.*. Delete the contents of the `$depsdir` directory except for the files `build.jl` and `.gitignore`, then rerun Pkg.build(\"OSQP\").")
    end
end

include("constants.jl")
include("types.jl")
include("interface.jl")
include("MPB_wrapper.jl")
include("MOI_wrapper.jl")
const Optimizer = MathOptInterfaceOSQP.Optimizer

end # module

ext: jl


text:
using Logging

@testset "Travis CI deploy configuration" begin; with_logger(NullLogger()) do
    # Regular tag build
    withenv("TRAVIS_CI" => "true",
            "TRAVIS_PULL_REQUEST" => "false",
            "TRAVIS_REPO_SLUG" => "JuliaDocs/Documenter.jl",
            "TRAVIS_BRANCH" => "master",
            "TRAVIS_TAG" => "v1.2.3",
            "DOCUMENTER_KEY" => "SGVsbG8sIHdvcmxkLg==",
        ) do
        cfg = Documenter.Travis()
        @test Documenter.deploy_folder(cfg; repo="github.com/JuliaDocs/Documenter.jl.git",
                  devbranch="master", devurl="dev", push_preview=true) == "v1.2.3"
        @test Documenter.documenter_key(cfg) === "SGVsbG8sIHdvcmxkLg=="
        @test Documenter.authentication_method(cfg) === Documenter.SSH
    end
    # Broken tag build
    withenv("TRAVIS_CI" => "true",
            "TRAVIS_PULL_REQUEST" => "false",
            "TRAVIS_REPO_SLUG" => "JuliaDocs/Documenter.jl",
            "TRAVIS_BRANCH" => "master",
            "TRAVIS_TAG" => "not-a-version",
            "DOCUMENTER_KEY" => "SGVsbG8sIHdvcmxkLg==",
        ) do
        cfg = Documenter.Travis()
        @test Documenter.deploy_folder(cfg; repo="github.com/JuliaDocs/Documenter.jl.git",
                  devbranch="master", devurl="dev", push_preview=true) === nothing
    end
    # Regular/broken devbranch build
    withenv("TRAVIS_CI" => "true",
            "TRAVIS_PULL_REQUEST" => "false",
            "TRAVIS_REPO_SLUG" => "JuliaDocs/Documenter.jl",
            "TRAVIS_BRANCH" => "master",
            "TRAVIS_TAG" => nothing,
            "DOCUMENTER_KEY" => "SGVsbG8sIHdvcmxkLg==",
        ) do
        cfg = Documenter.Travis()
        @test Documenter.deploy_folder(cfg; repo="github.com/JuliaDocs/Documenter.jl.git",
                  devbranch="master", devurl="hello-world", push_preview=true) == "hello-world"
        @test Documenter.deploy_folder(cfg; repo="github.com/JuliaDocs/Documenter.jl.git",
                  devbranch="not-master", devurl="hello-world", push_preview=true) === nothing
        @test Documenter.documenter_key(cfg) === "SGVsbG8sIHdvcmxkLg=="
    end
    # Regular pull request build
    withenv("TRAVIS_CI" => "true",
            "TRAVIS_PULL_REQUEST" => "42",
            "TRAVIS_REPO_SLUG" => "JuliaDocs/Documenter.jl",
            "TRAVIS_BRANCH" => "something",
            "TRAVIS_TAG" => nothing,
            "DOCUMENTER_KEY" => "SGVsbG8sIHdvcmxkLg==",
        ) do
        cfg = Documenter.Travis()
        @test Documenter.deploy_folder(cfg; repo="github.com/JuliaDocs/Documenter.jl.git",
                  devbranch="master", devurl="hello-world", push_preview=true) == "previews/PR42"
        @test Documenter.deploy_folder(cfg; repo="github.com/JuliaDocs/Documenter.jl.git",
                  devbranch="master", devurl="hello-world", push_preview=false) === nothing
        @test Documenter.documenter_key(cfg) === "SGVsbG8sIHdvcmxkLg=="
    end
    # Missing/broken environment variables
    withenv("TRAVIS_CI" => "true",
            "TRAVIS_PULL_REQUEST" => "false",
            "TRAVIS_REPO_SLUG" => "JuliaDocs/Documenter.jl",
            "TRAVIS_BRANCH" => "master",
            "TRAVIS_TAG" => "v1.2.3",
            "DOCUMENTER_KEY" => nothing,
        ) do
        cfg = Documenter.Travis()
        @test Documenter.deploy_folder(cfg; repo="github.com/JuliaDocs/Documenter.jl.git",
                  devbranch="master", devurl="hello-world", push_preview=false) === nothing
    end
end end

@testset "GitHub Actions deploy configuration" begin; with_logger(NullLogger()) do
    # Regular tag build with GITHUB_TOKEN
    withenv("GITHUB_EVENT_NAME" => "push",
            "GITHUB_REPOSITORY" => "JuliaDocs/Documenter.jl",
            "GITHUB_REF" => "refs/tags/v1.2.3",
            "GITHUB_ACTOR" => "github-actions",
            "GITHUB_TOKEN" => "SGVsbG8sIHdvcmxkLg==",
            "DOCUMENTER_KEY" => nothing,
        ) do
        cfg = Documenter.GitHubActions()
        @test Documenter.deploy_folder(cfg; repo="github.com/JuliaDocs/Documenter.jl.git",
                  devbranch="master", devurl="dev", push_preview=true) == "v1.2.3"
        @test Documenter.authentication_method(cfg) === Documenter.HTTPS
        @test Documenter.authenticated_repo_url(cfg) === "https://github-actions:SGVsbG8sIHdvcmxkLg==@github.com/JuliaDocs/Documenter.jl.git"
    end
    # Regular tag build with SSH deploy key (SSH key prioritized)
    withenv("GITHUB_EVENT_NAME" => "push",
            "GITHUB_REPOSITORY" => "JuliaDocs/Documenter.jl",
            "GITHUB_REF" => "refs/tags/v1.2.3",
            "GITHUB_ACTOR" => "github-actions",
            "GITHUB_TOKEN" => "SGVsbG8sIHdvcmxkLg==",
            "DOCUMENTER_KEY" => "SGVsbG8sIHdvcmxkLg==",
        ) do
        cfg = Documenter.GitHubActions()
        @test Documenter.deploy_folder(cfg; repo="github.com/JuliaDocs/Documenter.jl.git",
                  devbranch="master", devurl="dev", push_preview=true) == "v1.2.3"
        @test Documenter.authentication_method(cfg) === Documenter.SSH
        @test Documenter.documenter_key(cfg) === "SGVsbG8sIHdvcmxkLg=="
    end
    # Broken tag build
    withenv("GITHUB_EVENT_NAME" => "push",
            "GITHUB_REPOSITORY" => "JuliaDocs/Documenter.jl",
            "GITHUB_REF" => "refs/tags/not-a-version",
            "GITHUB_ACTOR" => "github-actions",
            "GITHUB_TOKEN" => "SGVsbG8sIHdvcmxkLg==",
            "DOCUMENTER_KEY" => nothing,
        ) do
        cfg = Documenter.GitHubActions()
        @test Documenter.deploy_folder(cfg; repo="github.com/JuliaDocs/Documenter.jl.git",
                  devbranch="master", devurl="dev", push_preview=true) === nothing
    end
    # Regular devbranch build with GITHUB_TOKEN
    withenv("GITHUB_EVENT_NAME" => "push",
            "GITHUB_REPOSITORY" => "JuliaDocs/Documenter.jl",
            "GITHUB_REF" => "refs/heads/master",
            "GITHUB_ACTOR" => "github-actions",
            "GITHUB_TOKEN" => "SGVsbG8sIHdvcmxkLg==",
            "DOCUMENTER_KEY" => nothing,
        ) do
        cfg = Documenter.GitHubActions()
        @test Documenter.deploy_folder(cfg; repo="github.com/JuliaDocs/Documenter.jl.git",
                  devbranch="master", devurl="hello-world", push_preview=true) == "hello-world"
        @test Documenter.deploy_folder(cfg; repo="github.com/JuliaDocs/Documenter.jl.git",
                  devbranch="not-master", devurl="hello-world", push_preview=true) === nothing
        @test Documenter.authentication_method(cfg) === Documenter.HTTPS
        @test Documenter.authenticated_repo_url(cfg) === "https://github-actions:SGVsbG8sIHdvcmxkLg==@github.com/JuliaDocs/Documenter.jl.git"
    end
    # Regular devbranch build with SSH deploy key (SSH key prioritized)
    withenv("GITHUB_EVENT_NAME" => "push",
            "GITHUB_REPOSITORY" => "JuliaDocs/Documenter.jl",
            "GITHUB_REF" => "refs/heads/master",
            "GITHUB_ACTOR" => "github-actions",
            "GITHUB_TOKEN" => "SGVsbG8sIHdvcmxkLg==",
            "DOCUMENTER_KEY" => "SGVsbG8sIHdvcmxkLg==",
        ) do
        cfg = Documenter.GitHubActions()
        @test Documenter.deploy_folder(cfg; repo="github.com/JuliaDocs/Documenter.jl.git",
                  devbranch="master", devurl="hello-world", push_preview=true) == "hello-world"
        @test Documenter.deploy_folder(cfg; repo="github.com/JuliaDocs/Documenter.jl.git",
                  devbranch="not-master", devurl="hello-world", push_preview=true) === nothing
        @test Documenter.authentication_method(cfg) === Documenter.SSH
        @test Documenter.documenter_key(cfg) === "SGVsbG8sIHdvcmxkLg=="
    end
    # Regular pull request build with GITHUB_TOKEN
    withenv("GITHUB_EVENT_NAME" => "pull_request",
            "GITHUB_REPOSITORY" => "JuliaDocs/Documenter.jl",
            "GITHUB_REF" => "refs/pull/42/merge",
            "GITHUB_ACTOR" => "github-actions",
            "GITHUB_TOKEN" => "SGVsbG8sIHdvcmxkLg==",
            "DOCUMENTER_KEY" => nothing,
        ) do
        cfg = Documenter.GitHubActions()
        @test Documenter.deploy_folder(cfg; repo="github.com/JuliaDocs/Documenter.jl.git",
                  devbranch="master", devurl="hello-world", push_preview=true) == "previews/PR42"
        @test Documenter.deploy_folder(cfg; repo="github.com/JuliaDocs/Documenter.jl.git",
                  devbranch="not-master", devurl="hello-world", push_preview=false) === nothing
        @test Documenter.authentication_method(cfg) === Documenter.HTTPS
        @test Documenter.authenticated_repo_url(cfg) === "https://github-actions:SGVsbG8sIHdvcmxkLg==@github.com/JuliaDocs/Documenter.jl.git"
    end
    # Regular pull request build with SSH deploy key (SSH key prioritized)
    withenv("GITHUB_EVENT_NAME" => "pull_request",
            "GITHUB_REPOSITORY" => "JuliaDocs/Documenter.jl",
            "GITHUB_REF" => "refs/pull/42/merge",
            "GITHUB_ACTOR" => "github-actions",
            "GITHUB_TOKEN" => "SGVsbG8sIHdvcmxkLg==",
            "DOCUMENTER_KEY" => "SGVsbG8sIHdvcmxkLg==",
        ) do
        cfg = Documenter.GitHubActions()
        @test Documenter.deploy_folder(cfg; repo="github.com/JuliaDocs/Documenter.jl.git",
                  devbranch="master", devurl="hello-world", push_preview=true) == "previews/PR42"
        @test Documenter.deploy_folder(cfg; repo="github.com/JuliaDocs/Documenter.jl.git",
                  devbranch="not-master", devurl="hello-world", push_preview=false) === nothing
        @test Documenter.authentication_method(cfg) === Documenter.SSH
        @test Documenter.documenter_key(cfg) === "SGVsbG8sIHdvcmxkLg=="
    end
    # Missing environment variables
    withenv("GITHUB_EVENT_NAME" => "push",
            "GITHUB_REPOSITORY" => "JuliaDocs/Documenter.jl",
            "GITHUB_REF" => "refs/tags/v1.2.3",
            "GITHUB_ACTOR" => "github-actions",
            "GITHUB_TOKEN" => nothing,
            "DOCUMENTER_KEY" => nothing,
        ) do
        cfg = Documenter.GitHubActions()
        @test Documenter.deploy_folder(cfg; repo="github.com/JuliaDocs/Documenter.jl.git",
                  devbranch="master", devurl="hello-world", push_preview=true) === nothing
    end
end end

@testset "Autodetection of deploy system" begin
    withenv("TRAVIS_REPO_SLUG" => "JuliaDocs/Documenter.jl",
            "GITHUB_REPOSITORY" => nothing,
        ) do
        cfg = Documenter.auto_detect_deploy_system()
        @test cfg isa Documenter.Travis
    end
    withenv("TRAVIS_REPO_SLUG" => nothing,
            "GITHUB_REPOSITORY" => "JuliaDocs/Documenter.jl",
        ) do
        cfg = Documenter.auto_detect_deploy_system()
        @test cfg isa Documenter.GitHubActions
    end
    withenv("TRAVIS_REPO_SLUG" => nothing,
            "GITHUB_REPOSITORY" => nothing,
        ) do
        cfg = Documenter.auto_detect_deploy_system()
        @test cfg === nothing
    end
end


text:
using PyPlot, Interpolations, DelimitedFiles

function boxcarn(A::AbstractArray, ntosmooth::Int)
           out = zeros(Float64, size(A))
           R = CartesianIndices(A)
           Ifirst, Ilast = first(R), last(R)
           I1 = oneunit(Ifirst)
           for I in R
               n, s = 0, zero(eltype(out))
               for J in max(Ifirst, I-ntosmooth*I1):min(Ilast, I+ntosmooth*I1)
                   s += A[J]
                   n += 1
               end
               out[I] = s/n
           end
           out
end

xy = [-0.0018781493357764578 -1.695087036188732
       0.03190563444800737 -1.2727897388914347
       0.08197434722858449 -0.7486257443884563
       0.13212322491983508 -0.3421037562986724
       0.1756527714154833 -0.22171323866239145
       0.22207970682546951 -0.3532409528172247
       0.2646243701328447 -0.7875343563902888
       0.29499541914796146 -1.3570487860742109
       0.3199267063673842 -1.9437127805771874
       0.3505611543747137 -2.8997652313330287
       0.3673156207054512 -3.4869445716903362
       0.3827187356848373 -4.091015803939533
       0.3940219880897847 -4.678538708199725
       0.3963582226294091 -4.10696289509849
       0.39715987173614276 -1.2833829592304165
       0.3995190105359596 1.2545808520384787
       0.4016720109940449 4.095052679798442
       0.40242785158039396 4.985856619331195
       0.4385364177737059 3.9965357306459
       0.4565964269354099 3.4934722858451672
       0.4774278515803939 2.9233566193311957
       0.5190563444800732 1.8335432890517636
       0.5522789738891433 1.0793346312414114
       0.5854557947778286 0.3923499770957397
       0.624083829592304 -0.2942911131470449
       0.6817796610169491 -0.9629237288135588
       0.7350893266147505 -1.1948579935868064
       0.7800160329821344 -1.1247995877233166
       0.823442510306917 -0.8531550618415018
       0.8654145671094824 -0.44714841960604623
       0.9046037562986715 0.04271644525881957
       0.9424072377462207 0.5661074209803028
       0.97882501145213 1.123024507558406
       1.0030920751259733 1.5111085661933128
       ];

x2 = LinRange(0, 1, 201) # new x-grid
figure()
plot(xy[:,1],xy[:,2])

itp = interpolate((xy[:,1],), xy[:,2], Gridded(Linear()))
yf = boxcarn(itp(x2), 1)
plot(x2, yf)
writedlm("func.txt", [x2 yf], ',')


text:
module SynchronousExamples

println("\nSynchronousExamples: Demonstrating the ability to simulate models with synchronous semantics")

using Modia
using Modia.Synchronous: sample, Clock, previous, hold

# Desired:
#   using ModiaMath: plot
#   using Test
#
# In order that these packages need not to be defined in the user environment, they are included via Modia:
using Modia.ModiaMath: plot

@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Modia.Test
end

@testset "Synchronous" begin

@model SynchronousOperators begin
  a = Var(size=())
  b = Var(size=())
  c = Var(start=0.0)
  obs=Float(start=0.0)
@equations begin
  a = Clock(0.1)
  b = sample(time, a)
  c = previous(c, a) + 0.5
  der(obs) = 1000000*(c-obs)
  end
end

#=
result = simulate(SynchronousOperators, 5.0, storeEliminated=false, logSimulation=false)
plot(result, ("obs"), heading="SynchronousOperators",figure=15)
@test result["obs"][end] == 25
=#

@model MassWithSpringDamper begin
  m=1
  k=1 # spring constant
  d=0.1 # damping coefficient
  x=Float(start=0.0) # Position
  v=Float(start=0.0) # Velocity
  f=Float() # Force
@equations begin
  der(x) = v
  m*der(v) = f - k*x - d*v
  end
end

@model SpeedControl begin
  @extends MassWithSpringDamper(k=0)
  @inherits v, f
  K = 5 # Gain of speed P controller
  vref = 100 # Speed ref.
  vd=Float(start=0.0)
  u=Float(start=0.0)
  fobs=Float(start=0.0)
@equations begin
  # speed sensor
  vd = sample(v, Clock(0.1))

  # P controller for speed
  u = K*(vref-vd)

  # force actuator
  f = hold(u)
  der(fobs) = 1000000*(f-fobs)
  end
end

result = simulate(SpeedControl, 5.0, storeEliminated=false, logSimulation=false)
plot(result, ("v", "fobs"), heading="SpeedControl", figure=15)
@test result["v"][end] == 98.03921568627446

@model SpeedControlPI begin
  @extends MassWithSpringDamper(k=0)
  @inherits v, f
  K = 2 # Gain of speed P controller
  Ti = 2 # Integral time
  vref = 100 # Speed reference
  dt=0.1 # sampling interval
  vd=Float()
  u=Float(start=0)
  e=Float()
  intE=Float(start=0)

  fobs=Float(start=0.0)
@equations begin
  # speed sensor
  vd = sample(v, Clock(dt))

  # PI controller for speed
  e = vref-vd
  intE = previous(intE, Clock(dt)) + e
  u = K*(e + intE/Ti)

  # force actuator
  f = hold(u)
  der(fobs) = 1000000*(f-fobs)
  end
end

result = simulate(SpeedControlPI, 5.0, storeEliminated=false, logSimulation=false)
plot(result, ("v", "fobs"), heading="SpeedControlPI", figure=16)
@show result["v"][end]
@test result["v"][end] == 100.2849917097788

@model ControlledMassBasic begin
  @extends MassWithSpringDamper(k=0) # k=100
  @inherits x, v, f
  KOuter = 10 # Gain of position PI controller
  KInner = 20 # Gain of speed P controller
  Ti = 10 # Integral time for pos. PI controller
  xref = 10 # Position reference
  xd=Float(size=())
  eOuter=Float(size=())
  intE=Float(start=0.0)
  uOuter=Float(size=())
  vd=Float(size=())
  vref=Float(size=())
  uInner=Float(start=0.0)
  fobs=Float(start=1990)
@equations begin
  # position sensor
  xd = sample(x, Clock(0.02))
  # outer PI controller for position
  eOuter = xref-xd
  intE = previous(intE, Clock(0.02)) + eOuter
  uOuter = KOuter*(eOuter + intE/Ti)
  # speed sensor
  vd = sample(v, Clock(0.02))
  # inner P controller for speed
  vref = uOuter
  uInner = KInner*(vref-vd)
  # force actuator
  f = hold(uInner)
  der(fobs) = 100000*(f-fobs)
  end
end

#result = simulate(ControlledMassBasic, 5.0, storeEliminated=false, logSimulation=true)
#plot(result, ("v", "fobs"), figure=4)

@model DiscretePIController begin
  K=0.1 # Gain
  Ti=1E10 # Integral time
  dt=1.0 # sampling interval
  ref=1 # set point
  u=Float(); ud=Float(size=())
  y=Float(); yd=Float(size=())
  e=Float(size=())
  intE=Float(start=0)

  fobs=Float(start=0)
@equations begin
  # sensor
  ud = sample(u, Clock(dt))
  # PI controller
  e = ref-ud
  intE = previous(intE, Clock(dt)) + e
  yd = K*(e + intE/Ti)
  # actuator
  y = hold(yd)

  der(fobs) = 100000*(y-fobs)
  end
end

@model SpeedControl2 begin
  m=MassWithSpringDamper(d=0.5, k=0)
  c=DiscretePIController(K=2, dt=0.1)
@equations begin
  connect(m.v, c.u)
  connect(c.y, m.f)
  end
end

#result = simulate(SpeedControl2, 5.0, storeEliminated=false, logSimulation=true)
#plot(result, ("v", "fobs"), figure=5)

@model SpeedControl3 begin
  m=MassWithSpringDamper(d=0.5, k=0)
  @extends DiscretePIController(K=2, dt=0.1)
  @inherits u, y, dt, e, intE
@equations begin
  u = m.v
  m.f = y
  end
end

#result = simulate(SpeedControl3, 5.0, storeEliminated=false, logSimulation=true)
#plot(result, ("v", "fobs"), figure=6)

end # testset
end


text:
"""

```
finalresults!(R::Result, options::Options, trajectory::Trajectory, samples::Samples)
```

Function that computes the final results of all the data computed by averaging
according to the sampling of each type of data, and converts to common units.

Computes also the final distribution functions and KB integrals

This function modified the values contained in the R data structure

"""
function finalresults!(
    R::Result,
    options::Options,
    trajectory::Trajectory,
    samples::Samples,
)

    # Setup the distance vector
    for i = 1:R.nbins
        R.d[i] = shellradius(i, options.binstep)
    end

    # Counters
    @. R.md_count = R.md_count / (samples.md * R.nframes_read)
    @. R.solute_atom = R.solute_atom / (samples.md * R.nframes_read)
    @. R.solvent_atom = R.solvent_atom / (samples.md * R.nframes_read)
    @. R.md_count_random = R.md_count_random / (samples.random * R.nframes_read)
    @. R.rdf_count = R.rdf_count / (samples.md * R.nframes_read)
    @. R.rdf_count_random = R.rdf_count_random / (samples.random * R.nframes_read)

    # Volumes and Densities
    R.volume.total = R.volume.total / R.nframes_read
    R.density.solvent = R.density.solvent / R.nframes_read
    R.density.solute = R.density.solute / R.nframes_read

    R.volume.shell = R.volume.shell / R.nframes_read
    R.volume.domain = R.volume.domain / R.nframes_read
    R.volume.bulk = R.volume.bulk / R.nframes_read

    R.density.solvent_bulk = R.density.solvent_bulk / R.nframes_read

    #
    # Computing the distribution functions and KB integrals, from the MDDF
    # and from the RDF
    #

    for ibin = 1:R.nbins

        # For the MDDF

        if R.md_count_random[ibin] > 0.0
            R.mddf[ibin] = R.md_count[ibin] / R.md_count_random[ibin]
            for i = 1:trajectory.solute.natomspermol
                R.solute_atom[ibin, i] = R.solute_atom[ibin, i] / R.md_count_random[ibin]
            end
            for j = 1:trajectory.solvent.natomspermol
                R.solvent_atom[ibin, j] = R.solvent_atom[ibin, j] / R.md_count_random[ibin]
            end
        end
        if ibin == 1
            R.sum_md_count[ibin] = R.md_count[ibin]
            R.sum_md_count_random[ibin] = R.md_count_random[ibin]
        else
            R.sum_md_count[ibin] = R.sum_md_count[ibin-1] + R.md_count[ibin]
            R.sum_md_count_random[ibin] =
                R.sum_md_count_random[ibin-1] + R.md_count_random[ibin]
        end
        R.kb[ibin] =
            units.Angs3tocm3permol *
            (1 / R.density.solvent_bulk) *
            (R.sum_md_count[ibin] - R.sum_md_count_random[ibin])

        # For the RDF

        if R.rdf_count_random[ibin] > 0.0
            R.rdf[ibin] = R.rdf_count[ibin] / R.rdf_count_random[ibin]
        end
        if ibin == 1
            R.sum_rdf_count[ibin] = R.rdf_count[ibin]
            R.sum_rdf_count_random[ibin] = R.rdf_count_random[ibin]
        else
            R.sum_rdf_count[ibin] = R.sum_rdf_count[ibin-1] + R.rdf_count[ibin]
            R.sum_rdf_count_random[ibin] =
                R.sum_rdf_count_random[ibin-1] + R.rdf_count_random[ibin]
        end
        R.kb_rdf[ibin] =
            units.Angs3tocm3permol *
            (1 / R.density.solvent_bulk) *
            (R.sum_rdf_count[ibin] - R.sum_rdf_count_random[ibin])

    end

    nothing
end


text:
using Pkg
pkg"activate ."

using Genie, HTTP
import Genie.Router: route, POST, @params

route("/jsonpayload", method = POST) do
  @show @params(:JSON_PAYLOAD)
end
Genie.AppServer.startup()

HTTP.request("POST", "http://localhost:8000/jsonpayload", [("Content-Type", "application/json; charset=utf-8")], """{"greeting":"hello"}""")

exit(0)


text:
type Pixel
x::Int64
y::Int64
color::Int64
end

max_stars_repo_path: Module 2/Chapter03/pixel.jl

max_issues_repo_path: Module 2/Chapter03/pixel.jl

max_forks_repo_path: Module 2/Chapter03/pixel.jl


text:
using BenchmarkTools
using Aligator

include("singlepath.jl")
include("multipath.jl")

# singlepath = [:cohencu, :freire1, :freire2, :(petter(1)), :(petter(2)), :(petter(3)), :(petter(4))]
# multipath = [:divbin, :euclidex, :fermat, :knuth, :lcm, :mannadiv, :wensley]

singlepath = [:cubes, :eucliddiv, :freire1, :freire2, :(petter(2)), :(petter(3)), :(petter(4)), :(petter(22)), :(petter(23)), :psolv1s, :psolv2s]
multipath = [:euclidex, :fermat, :knuth, :lcm, :wensley, :divbin, :psolv2m, :psolv3m, :psolv4m, :psolv10m]

group = BenchmarkGroup()
group["singlepath"] = BenchmarkGroup()
group["multipath"] = BenchmarkGroup()

function create(suite, instances)
    for loop in instances
        suite[string(loop)] = @benchmarkable aligator($loop) samples=1 seconds=60 evals=5
    end
end

create(group["singlepath"], singlepath)
create(group["multipath"], multipath)

# tune!(group)
# results = run(group)

# res = show(results)

# results["singlepath"]
# results["multipath"]


text:
module ConvHMM

using BSON
using ExpFamilyDistributions
using LinearAlgebra
using PaddedViews
using MarkovModels

export Regressors1D

include("regressor.jl")

export ARNormal1D
export ARNormal1DSet
export DARNormal1D

export accstats_λ
export accstats_ξ
export accstats_h
export elbo
export loglikelihood
export predict
export save
export update_λ!
export update_h!

include("arnormal.jl")
include("darnormal.jl")

end

max_stars_repo_name: BUTSpeechFIT/ConvHMM

max_issues_repo_name: BUTSpeechFIT/ConvHMM

max_forks_repo_name: BUTSpeechFIT/ConvHMM


text:
using DelimitedFiles

Random.seed!(1234)

δ = 1.0
φ = 0.95
ση = 1.0
σε = 1.0

θ0 = [δ;atanh(φ);ση;σε]

β = [0.0;1.0]

N = 1000

@testset "Kalman filter regression tests" begin
    arh = ARH(β)

    X,Y = Filtering.simulate(arh,θ0,N,fill([1.0],N))

    kf = kalman_filter(arh,θ0,Y,fill([1.0],length(Y)))

    @testset "Filter" begin

        xf0 = readdlm("regressions/ar/kf_xf.csv")[:,1]
        Pf0 = readdlm("regressions/ar/kf_Pf.csv")[:,1]
        xp0 = readdlm("regressions/ar/kf_xp.csv")[:,1]
        Pp0 = readdlm("regressions/ar/kf_Pp.csv")[:,1]
        Σ0  = readdlm("regressions/ar/kf_Σ.csv")[:,1]
        ε0  = readdlm("regressions/ar/kf_ε.csv")[:,1]

        @test all(vcat(kf.xf...) .== xf0)
        @test all(vcat(kf.Pf...) .== Pf0)
        @test all(vcat(kf.xp...) .== xp0)
        @test all(vcat(kf.Pp...) .== Pp0)
        @test all(vcat(kf.Σ...)  .== Σ0)
        @test all(vcat(kf.ε...)  .== ε0)
    end

    ks = kalman_smoother(arh,θ0,kf)

    @testset "Smoother" begin

        x0  = readdlm("regressions/ar/ks_x.csv")[:,1]
        P0  = readdlm("regressions/ar/ks_P.csv")[:,1]

        @test all(vcat(ks.x...)  .== x0)
        @test all(vcat(ks.P...)  .== P0)
    end

end


text:
@reexport module ModelKit

export @var,
    @unique_var,
    AbstractHomotopy,
    AbstractSystem,
    CompiledHomotopy,
    CompiledSystem,
    Expression,
    Homotopy,
    InterpretedHomotopy,
    InterpretedSystem,
    System,
    TaylorVector,
    Variable,
    coefficients,
    coeffs_as_dense_poly,
    degree,
    degrees,
    differentiate,
    dense_poly,
    evaluate,
    evaluate!,
    evaluate_and_jacobian!,
    expand,
    exponents_coefficients,
    expressions,
    horner,
    is_homogeneous,
    jacobian,
    jacobian!,
    multi_degrees,
    parameters,
    nparameters,
    nvariables,
    monomials,
    optimize,
    parameters,
    subs,
    support_coefficients,
    rand_poly,
    taylor!,
    to_dict,
    to_number,
    variables,
    variable_groups,
    vectors

using ..DoubleDouble: ComplexDF64

using Arblib:
    Arblib,
    Arb,
    ArbRef,
    ArbVector,
    ArbRefVector,
    Acb,
    AcbRef,
    AcbVector,
    AcbRefVector,
    AcbRefMatrix
import LinearAlgebra
import MultivariatePolynomials:
    MultivariatePolynomials,
    coefficients,
    degree,
    differentiate,
    nvariables,
    monomials,
    subs,
    variables
using Parameters: @unpack
const MP = MultivariatePolynomials

include("model_kit/symengine.jl")
include("model_kit/symbolic.jl")

include("model_kit/taylor_vector.jl")

include("model_kit/instructions.jl")
include("model_kit/slp_interpreter.jl")
include("model_kit/slp_compiler.jl")
#
include("model_kit/abstract_system_homotopy.jl")
include("model_kit/compiled_system_homotopy.jl")
include("model_kit/interpreted_system_homotopy.jl")

end # module


text:
InputFolder = './Images/FISH/';
OutputFolder = './Results/Images/FISH_spt/';

@iA = '*_C01*.tif';

@fxg_mLapThrBinWatTiles [iA] > [L];
params.GRad = 2;
params.Thr = 1;
params.GaussianD = -1;
params.MinArea = 9;
/endf

@fxm_sMarkObjCentroids [L] > [S];
params = [];
/endf

/show iA > S;
/keep S > tif;

max_stars_repo_name: SebastienTs/LOBSTER

max_issues_repo_name: SebastienTs/LOBSTER

max_forks_repo_name: SebastienTs/LOBSTER


text:
# basisfunctions.jl

"""
    fourier(k::Int)

Implements the Fourier basis of functions ϕ_k, defined, if k is odd, by
ϕ_k(x)=sqrt(2)sin((k+1)π*x) and if k is even, by ϕ_k(x)=sqrt(2)cos(k*x).
This function is defined for k≥0 and ϕ_0≡1.

# Examples

```julia
x = -1.0:0.01:2.0
y = fourier(3).(x)
```
"""
function fourier(k::Int)
  k ≥ 0 || throw(AssertionError("k should be non-negative"))
  sqrttwo = sqrt(2.0)
  if k == 0
    return x -> 1.0
  elseif isodd(k)
    return x -> sqrttwo * sinpi(float(k+1)*x)
  else
    return x -> sqrttwo * cospi(float(k)*x)
  end
end

"""
    faberschauderone(x::Float64)

Implements the first Faber-Schauder function defined by 1-2x for 0≤x≤1/2 and
-1+2x for 0.5≤x≤1, and is 1-periodically extended to all x∈R.

#See also: faberschauder

#Warning

Note the difference between faberschauderone and faberschauder. The first is a
function that takes a Float64 and returns a Float64, the second takes (j,k) and
returns a anonymous function that takes a Float64 and returns a Float64.

#Examples

```julia
x=-2.0:0.001:2.0
y=faberschauderone.(x)
```
"""
function faberschauderone(x::Float64)
    y = mod(x, 1.0)
    if 0≤y≤0.5
        return 1.0 - 2*y
    else
        return -1.0 + 2*y
    end
end

"""
    faberschauder(j::Int, k::Int)

Implements the k-th Faber-Schauder function of level j. Here, j≥0 and 1≤k≤2^j.
It is a one periodic function and defined on [0,1] by 2^(j+1)(x-(k-1)2^(-j)) on
(k-1)2^(-j)≤x≤(k-1/2)2^(-j) and 1 - 2^(j+1)(x-(k-1/2)2^(-j)) on
[(k-1/2)2^(-j), k2^(-j)] and zero outside these intervals.

# See also: faberschauderone.

#Warning

Note the difference between faberschauderone and faberschauder. The first is a
function that takes a Float64 and returns a Float64, the second takes (j,k) and
returns a anonymous function that takes a Float64 and returns a Float64.

#Example using PyPlot
```julia
using PyPlot
clf()
j=3
x=-2.0:0.001:2.0
for k in 1:2^j
    y = faberschauder(j,k).(x)
    plot(x,y)
end
```
"""
function faberschauder(j::Int,k::Int)
  j ≥ 0 || throw(AssertionError("j should be a nonnegative integer."))
  1≤k≤2^j || throw(AssertionError("k should be an integer between 1 and 2^j."))
  twotothepowerminusj = 2.0^-j
  k_float = float(k)
  j_float = float(j)
  return function(x::Float64)
    y = mod(x,1.0)
    if y ≤ (k_float-1.0)*twotothepowerminusj || y ≥ k_float*twotothepowerminusj
      return 0.0
    elseif (k_float-1.0)*twotothepowerminusj < y ≤ (k_float-0.5)*twotothepowerminusj
      return 2.0^(j_float+1.0)*y-2*k_float+2.0
    else
      return 2.0*k_float - 2.0^(j_float+1.0) * y
    end
  end
end


text:
##############################################################################################################################
#															     #
# J4PR - a small library and package wrapper written in Julia by Cornel Cofaru at 0x0α Research                              #
# 															     #
##############################################################################################################################
VERSION >= v"0.6.0" && __precompile__(true)

"""
**J4PR**

A small library and package wrapper written in Julia for A.I. hacking.
It relies on various packages from the Julia ecosystem and a quite simple
constructs to provide a consistent way of solving problems.
Documentation for the main types and functions is provided.
Digging into the code will yield more low-level documentation
in the form of comments :) .

To begin the documentation process regarding the library, type
`?DataCell`, `?FunctionCell`, `?PipeCell` and `?datacell`.
Also, start exploring **src/j4pr.jl** as well as the included files.
The code is commented as much as possible and should at least partially
make up for the lack in documentation.

For bugs, remarks and other related issues, e-mail at **j4pr@oxoaresearch.com**.
Please note that support is not oficially available however any feedback is welcomed.
"""
module j4pr

    	##############################################################################################################################
    	# Dependencies														     #
    	##############################################################################################################################

	# Imports section
	import Base: Test, convert, promote_rule, show, +, -, *, |>, ~, vcat, hcat, getindex, setindex!, size, endof, ndims,
		     isnan, deleteat!, values, start, next, done, length, eltype, strip, interrupt
	import StatsBase: sample
	import MLDataPattern: nobs, getobs, datasubset, targets, gettargets
	import UnicodePlots, LossFunctions, Distances, MultivariateStats, Clustering, DecisionTree, NetworkLearning, MLKernels

	# Using section
	using DataStructures: SortedDict
	using StatsBase: countmap, fit, Histogram
	using Reexport, StaticArrays, DataArrays, Compat, LearnBase, MLLabelUtils, MLLabelUtils.LabelEncoding

	# If Julia > 0.7, skip some packages
	if v"0.6" <= VERSION < v"0.7-"
		import LIBSVM
		using RDatasets, Images, ImageInTerminal
	end

    	##############################################################################################################################
    	# Global variables and logging configuration										     #
	##############################################################################################################################
	global const j4pr_version = "0.1.3-alpha"                      				# The current version of j4pr

	oinfoglobal = j4pr_version*Dates.format(Dates.now(), " dd-mm-YYYY")			# Define information string

    	##############################################################################################################################
    	# Export														     #
    	##############################################################################################################################
	@reexport using MLDataPattern

	# Some packages do not work for Julia > 0.6
	if v"0.6" <= VERSION < v"0.7-"
    		export
		rdataset,									# R datasets
		im2targets, targets2im,								# From Images.jl Arrays to DataCells and back
		libsvm,										# LIBSVM classifier/regressor
		kNNClassifier, knn, knnr,							# k-nearest neighbours classifier and regressor
		ROC, findop, changeop!, simpleop,						# Operating point optimization
		rocplot
	end

	export
		# [core]
		@async_cell, @remote_cell,							# Macros
	   	AbstractCell, DataCell, FunctionCell, PipeCell,					# The basic types and constructors
		Model, ModelProperties,								# Model related types
	   	DataGenerator,									# Data generator submodule
		functioncell, datacell,                                                         # Functions to create cells
	   	getx, getx!, gety, gety!, getf, getf!,                               		# Functions the get various cell fields
	   	nvars, size, start, next, done, eltype, length, ndims, endof,			# Various useful functions for data cells
	   	classsizes, nclasssizes, classnames, countmapn,
		nclass, deleteat!, deleteat, idx,
		varsubset, labelencn,								# Lazy subset of variables
	   	strip, 										# Return tuple from data cell
	   	pipestack, pipeparallel, pipeserial,                                          	# Create pipes
	   	flatten, mat, 		                                                       	# Flatten a pipe made out of data
	   	addlabels, unlabel, labelize,							# Label-related functionality
	   	interrupt,

		# [lib]
		countapp, countapp!, countappw, countappw!, 					# Counting utils
	  	pintgen,									# Generate numbers based on priors
		gini, misclassification,							# Purity utility functions
		linearsplit, densitysplit,							# Split vectors according to different criteria
		confusionmatrix,								# Confusion matrix

		# [lib/data]
		DataGenerator, 									# Small sub-module that generates some datasets
		cslice,										# Class slicing
	   	sample,										# Data sampling (sub/super sampling)
	   	scaler!,									# Data scaling
	   	filterg,									# Data processing (for missing values)
	   	filterdomain!,									# Data domain filtering
	   	ohenc, ohenc_integer, ohenc_binary,						# One hot encoder
	   	lineplot, scatterplot, densityplot1d, 						# Plots
		densityplot2d,

		# [lib/ext]
		### fa, far,									# Factor Analysis and reconstruction
		### kpca, kpcar,								# Kernel-PCA and reconstruction
		dist,										# Distances
		whiten,										# Data Whitening
		pca, pcar,									# PCA and reconstruction
		ppca, ppcar,									# Probabilistic PCA and reconstruction
		ica,										# ICA
		mds,										# MDS
		lda, ldasub,									# Multiclass LDA, Subspace LDA
		lr,										# Linear regression (Linear+Ridge)
		kmeans, kmeans!,								# K-means clustering
		kmedoids, kmedoids!,								# K-medoids clustering
		affinityprop,									# Affinity propagation clustering
		dbscan,										# DBSCAN clustering
		tree, randomforest,								# Decision tree, random forest classifiers
		treer, randomforestr,								# Decisiontree, random forest regression
		aboostump,									# Adaptively boosted stump classifier
		mlkernel, 									# Machine Learning kernels from MLKernels.jl

		# [lib/sup]
		loss,										# Calculate losses based on MLLabelUtils.jl and LossFunctions.jl
		LinDiscClassifier, lindisc,							# Linear discriminant classifier
		QuadDiscClassifier, quaddisc,							# Quadratic discriminant classifier
		ParzenClassifier, parzen,							# Parzen window density estimator/regressor/classifier
		ClassifierCombiner, votecombiner, wvotecombiner,				# Label combiners
		naivebayescombiner, bkscombiner,
		meancombiner, wmeancombiner, productcombiner,					# Continuous output combiners
		mediancombiner, dtcombiner,
		RandomSubspace, randomsubspace,							# Random sub-space ensemble
		AdaBoost, adaboost,								# AdaBoost ensemble
		DecisionStump, stump, stumpr,							# Decision stump classifier and regressor
		networklearner,									# Network learning

		# [lib/unsup]
		kernelize!, kernelize, kernel							# Construct kernels

    	##############################################################################################################################
    	# Include														     #
    	##############################################################################################################################

    	# [core]
    	include("core/abstractcell.jl")                                                         # AbstractCell related
    	include("core/datacell.jl")                                                         	# DataCell related
    	include("core/functioncell.jl")                                                         # FunctionCell related
    	include("core/pipecell.jl")                                                         	# PipeCell related
	include("core/calls.jl")                                                                # < call > methods
    	include("core/convert_promote.jl")                                                      # Conversion and Promotion rules in J4PR
    	include("core/printers.jl")                                                             # Text output (e.g. in REPL) for J4PR objects
    	include("core/coreutils.jl")                                                            # Utility functions (manipulate data and pipes)
    	include("core/macros.jl") 								# Macros
    	include("core/parallel.jl") 								# Code associated with Cell parallelism
	include("core/checks.jl")								# Various checks needed for the library
	include("core/version.jl")								# Fancy J4PR info

	# [lib]
	include("lib/libutils.jl")								# Utility functions for learning
	include("lib/mldatapattern.jl")								# MLDataPattern interface

	# [lib/data] e.g. data manipulation
		include("lib/data/cslice.jl")							# Class slicing (e.g. select classes)
		include("lib/data/labelutils.jl")                                               # Manipulation of datacell labels
		include("lib/data/datagen.jl")							# Datasets in datacell format
		include("lib/data/sample.jl")							# Data sampling
		include("lib/data/scale.jl")							# Data scaling
		include("lib/data/filter.jl")							# Data filter: generic, can be used for missing values
		include("lib/data/filterdomain.jl")						# Data domain filtering
		include("lib/data/ohenc.jl")							# One-hot encoding
		include("lib/data/textanalysis_interface.jl")					# TextAnalysis.jl interface (so far empty)
		include("lib/data/videoio_interface.jl")					# VideoIO.jl interface (so far empty)

	# [lib/ext] e.g external algorithms
		### include("lib/ext/kpca.jl")							# Kernel-PCA (MultivariateStats.jl)
		### include("lib/ext/fa.jl")							# Factor Analysis (MultivariateStats.jl)
		include("lib/ext/dist.jl")							# Distances (Distances.jl)
		include("lib/ext/whiten.jl")							# Data whitening (MultivariateStats.jl)
		include("lib/ext/pca.jl")							# PCA (MultivariateStats.jl)
		include("lib/ext/ppca.jl")							# Probabilistic PCA (MultivariateStats.jl)
		include("lib/ext/ica.jl")							# ICA (MultivariateStats.jl)
		include("lib/ext/mds.jl")							# MDS (MultivariateStats.jl)
		include("lib/ext/lda.jl")							# Multiclass LDA, Subspace LDA (MultivariateStats.jl)
		include("lib/ext/lr.jl")							# Linear regression (MultivariateStats.jl)
		include("lib/ext/kmeans.jl")							# K-means clustering (Clustering.jl)
		include("lib/ext/kmedoids.jl")							# K-medoids clustering (Clustering.jl)
		include("lib/ext/affinityprop.jl")						# Affinity propagation clustering (Clustering.jl)
		include("lib/ext/dbscan.jl")							# DBSCAN clustering (Clustering.jl)
		include("lib/ext/decisiontree.jl")						# DecisionTree, Random Forest, Boosted stumps (DecisionTree.jl)
		include("lib/ext/mlkernel.jl") 							# Machine Learning kernels (MLKernels.jl)

	# [lib/unsup] e.g. unsupervised learning
		include("lib/unsup/kernel.jl")							# Lightweight way of constructing kernels

	# [lib/sup] e.g. supervised learning
		include("lib/sup/lindisc.jl")							# Linear discriminant classifier
		include("lib/sup/quaddisc.jl")							# Quadratic discriminant classifier
		include("lib/sup/parzen.jl")							# Parzen window density estimator/regressor/classifier
		include("lib/sup/loss.jl")							# Calculate losses based on MLLabelUtils.jl and LossFunctions.jl
		include("lib/sup/combiner.jl")							# Combine classifier outputs (labels, class posteriors, etc.)
		include("lib/sup/randomsubspace.jl")						# Random subspace ensemble framework
		include("lib/sup/adaboost.jl")							# AdaBoost ensemble framework
		include("lib/sup/stump.jl")							# Decision stump classifier and regressor
		include("lib/sup/networklearner.jl")						# Network learning

	# Some packages do not work for Julia > 0.6
	if v"0.6" <= VERSION < v"0.7-"
		include("lib/data/images_interface.jl")						# Functions to transform to/from Images.jl representations from/to DataCells
		include("lib/data/rdataset.jl")							# Loads different R datasets
		include("lib/ext/libsvm.jl")							# LIBSVM classifier/regressor (LIBSVM.jl)
		include("lib/sup/knn.jl")							# kNN classifier/regressor/density estimator
		include("lib/sup/roc.jl")							# ROC Analysis and operating point related functionality
	end

	# [exp] Experimental stuff
		include("exp/plotting.jl")							# Plots for labeled/unlabeled datasets (UnicodePlots.jl)
												# limited and biased view on the library's performance)

    	##############################################################################################################################
    	# Post-loading steps 													     #
    	##############################################################################################################################

	# Print welcome message
    	version(j4pr_version)

end


text:
import Base.fetch
abstract Router

macro namespace(path, block)
    block.head != :block && error("Second arg is not a block")
    bargs = map(block.args) do line
        if line.args[1] == :handle!
            line.args[5] = path * line.args[5]
        end
        esc(line)
    end

    block.args = bargs
    block
end

"""
# Mount
takes one routing tree `r2` and attaches it to another routing tree `r1`
the mount will applied according to the specified path `to`

return returns the leaf node
"""
function mount!(r1::Router, r2::Router; to="/")
  error("mount! not implemented.")
end

"""
# handle!

attaches the handler::Function to the handlemap contained on the leaf endpoint
of the path.
"""
function handle!(
  fn::Handler,
  ep::Endpoint,
  method::Verb
)
  error("handle!(::Handler, ::Endpoint, ::Verb) implemented")
end

function handle!(
  fn::Handler,
  router::Router,
  method::Verb,
  path::Vector
)
  error("handle!(::Handler, ::Router, ::Verb, ::Vector) implemented")
end

handle!(
  fn::Handler,
  router::Router,
  method::Verb,
  path::AbstractString
) = handle!(fn, router, method, parseurl(path))

handle!(
  router::Router,
  fn::Handler,
  method::Verb,
  path::AbstractString
) = handle!(fn, router, method, parseurl(path))

function use!(router::Router, fn::Function)
  error("use! is not implemented.")
end

use!(fn::Function, r::Router) = use!(r, fn)
fetch(r::Router, tokens::Vector) =
finalize!(r::TRouter) = error("finalize! not implemented.")

# Router specific helper functions


text:
# This file contains all utility functions that can be used for testing
# Modelica models with the built-in unit testing support of Julia

"""
    regressionTest(
        omc:: OMCSession, name:: String, refdir:: String;
        relTol:: Real = 1e-6, variableFilter:: String = "", outputFormat="csv",
        resample=true
    )

Performs a regression test that ensures that variable values in the simulation
output are approximately equal to variables in the reference file in the
directory given by `refdir`.
Both the simulation output and the reference file must have the standard name
`"\$(name)_res.\$outputFormat"`.
This function also assumes that the simulation a simulation of the model named
`name` has already been run with
[`simulate(omc:: OMCSession, name::String, settings:: Dict{String, Any})`](@ref).

The test consists of the following checks performed with `@test`:

* Are there no variables in the simulation output that have no corresponding
    variable in the reference file? For this check, variables are ignored that
    occur in the simulation file but should not have been selected by the
    regex `variableFilter`. This can happen in OpenModelica, because sometimes
    alias variables are added to the output even though they do not match the
    `variableFilter`.
* Are there no variables in the reference file that have no corresponding
    variable in the simulation output?
* Is the intersection of common variables in both files nonempty?
* Is the length of the simulation output equal to the length of the reference
    file?
* If the length is not equal and `resample` is `true`, resample the data to
    equal length.
* Are there any variables in the simulation file that do not satisfy
    `isapprox(act, ref, rtol=relTol)` for all values?

NOTE: There is a OM scripting function `compareSimulationResults()` that could
be used for this task, but it is not documented and does not react to changes
of its `relTol` parameter in a predictable way.
"""
function regressionTest(
        omc:: OMCSession, name:: String, refdir:: String; relTol:: Real = 1e-6,
        variableFilter:: String = "", outputFormat="csv", resample=true
    )
    actname = "$(name)_res.$outputFormat"
    # make refdir relative to CWD of OMC
    omcrefdir = relpath(refdir, sendExpression(omc, "cd()"))
    refname = joinpath(omcrefdir, actname)
    actvars = sendExpression(omc, "readSimulationResultVars(\"$actname\")")
    refvars = sendExpression(omc, "readSimulationResultVars(\"$refname\")")
    missingRef = setdiff(Set(actvars), Set(refvars))
    # ignore variables without reference if they should not have been selected in the first place
    if length(variableFilter) > 0
        varEx = Regex("^($(variableFilter))\$")
        missingRef = filter(x -> occursin(varEx,x), missingRef)
    end
    @test isempty(missingRef)
    missingAct = setdiff(Set(refvars), Set(actvars))
    @test isempty(missingAct)
    # if variable sets differ, we should only check the variables that are present in both files
    vars = collect(intersect(Set(actvars), Set(refvars)))
    @test !isempty(vars)

    wd = sendExpression(omc, "cd()")
    actpath = joinpath(wd, actname)
    refpath = joinpath(wd, refname)
    function readSimulationResultMat(fn:: String)
        # TODO can be replaced by DataFrame(MAT.matread(actpath)) once
        # https://github.com/JuliaIO/MAT.jl/pull/132 is merged
        fnrel = relpath(fn, sendExpression(omc, "cd()"))
        data = sendExpression(omc, "readSimulationResult(\"$fnrel\", {$(join(vars, ", "))})")
        df = DataFrame(Dict(zip(vars, data)))
        return df
    end
    if outputFormat == "csv"
        actdata = DataFrame(CSV.File(actpath))
        refdata = DataFrame(CSV.File(refpath))
    elseif outputFormat == "mat"
        actdata = readSimulationResultMat(actpath)
        refdata = readSimulationResultMat(refpath)
    else
        throw(MoSTError("unknown output format $outputFormat", ""))
    end
    # check if length is equal
    @test size(actdata, 1) == size(refdata, 1)
    n = min(size(actdata, 1), size(refdata, 1))
    equalsize = size(actdata, 1) == size(refdata, 1)
    # check if start and end time is equal
    @test actdata[1, "time"] == refdata[1, "time"]
    @test actdata[end, "time"] == refdata[end, "time"]
    equalstart = actdata[1, "time"] == refdata[1, "time"]
    equalend = actdata[end, "time"] == refdata[end, "time"]
    # resample if required
    if equalstart && equalend && !equalsize && resample
        actdata = resamplequi(actdata, n)
        refdata = resamplequi(refdata, n)
    end
    # find unequal variables
    unequalVars = filter(x -> !isapprox(actdata[1:n,Symbol(x)], refdata[1:n,Symbol(x)]; rtol=relTol), vars)
    @test isempty(unequalVars)
end

"""
    resamplequi!(data, n)

Performs equidistant resampling to `n` data points according to the `time`
column in `data`.

For each point on the resampled time axis, the first value from the original
data is selected whose time stamp is greater or equal to that point in time.
"""
function resamplequi(data, n)
    # equidistant resampling according to time axis info
    time = data[!, "time"]
    newtime = range(time[1], time[end], length=n)
    # compute permutation array
    perm = fill(1, n)
    idx = 1
    for i in 1:n
        # advance index in old time array until we reach the next index
        # whose time stamp lies in the present or future
        while time[idx] < newtime[i]
            idx += 1
        end
        perm[i] = idx
    end
    # apply permutation to rows
    res = DataFrame([v => (v == "time" ? newtime : data[!, Symbol(v)][perm]) for v in names(data)])
    return res
end

"""
    testmodel(omc, name; override=Dict(), refdir="regRefData", regRelTol:: Real= 1e-6)

Performs a full test of the model named `name` with the following steps:

* Load the model using [`loadModel(omc:: OMCSession, name:: String)`](@ref) (called inside `@test`)
* Simulate the model using [`simulate(omc:: OMCSession, name::String, settings:: Dict{String, Any})`](@ref) (called inside `@test`)
* If a reference file exists in `refdir`, perform a regression test with
    [`regressionTest(omc:: OMCSession, name:: String, refdir:: String; relTol:: Real = 1e-6, variableFilter:: String = "", outputFormat="csv")`](@ref).
"""
function testmodel(omc, name; override=Dict(), refdir="regRefData", regRelTol:: Real= 1e-6)
    @test isnothing(loadModel(omc, name))
    settings = getSimulationSettings(omc, name; override=override)
    outputFormat = settings["outputFormat"]
    varfilter = settings["variableFilter"]
    @test isnothing(simulate(omc, name, settings))

    # compare simulation results to regression data
    if isfile("$(joinpath(refdir, name))_res.$outputFormat")
        regressionTest(omc, name, refdir; relTol=regRelTol, variableFilter=varfilter, outputFormat=outputFormat)
    else
        write(Base.stderr, "WARNING: no reference data for regression test of $name\n")
    end
end

max_stars_repo_name: THM-MoTE/MoST

max_issues_repo_name: THM-MoTE/MoST

max_forks_repo_name: THM-MoTE/MoST


text:
export Mixture, sample_mixture, log_pdf_mixture

struct Mixture
    "Number of Gaussian kernels"
    Nψ::Int64

    "Dimension of the state"
    Nx::Int64

    "Array of Nψ dimensions"
    dist::Array{MvNormal,1}

    "Vector of weights for each mode"
    w::Array{Float64,1}
end

function Mixture(Nx::Int64)

    @assert Nx == 3 "Method is not implemented for Nx != 3 "
    # Define number of modes
    Nψ = 2^Nx

    # Define weights, mean and covariance
    w   = rand(Nψ)
    w ./= sum(w)

    dist = []

    loc = [-4 -4.0 -4; -4 4 -4; 4 -4 -4; 4 4 -4; -4 -4 4; -4 4 4; 4 -4 4; 4 4 4]

    @assert size(loc) == (Nψ, Nx)
    collect(Iterators.product(-4:4, -4:4))

    @inbounds for i=1:Nψ
        # Use identity for the covariance matric
        push!(dist, MvNormal(loc[i,:], 1.0))
    end

    return Mixture(Nψ, Nx, dist, w)
end

function sample_mixture(M::Mixture, Ne::Int64)
    out = zeros(M.Nx, Ne)
    @inbounds for (i,di) in enumerate(M.dist)
        out += M.w[i]*rand(di,Ne)
    end
    return out
end

function log_pdf_mixture(M::Mixture, X::Array{Float64,2})
    NxX, Ne = size(X)
    @assert NxX == M.Nx "Dimension of the sample is wrong"
    logpdfX = zeros(Ne)
    @inbounds for (i,di) in enumerate(M.dist)
        logpdfX += M.w[i]*pdf(di, X)
    end
    @. logpdfX = log(logpdfX)

    return logpdfX
end


text:
"""
    gmax = gaugemax(gauge::VisClaw.Gauge)

Maximal values and their occurrence times in a gauge
"""
function gaugemax(gauge::VisClaw.Gauge)

    max_AMRlevel = findmax(gauge.AMRlevel)[1]
    max_eta = NaN
    t_eta = NaN
    if !isempty(gauge.eta)
        max_eta, tind = findmax(gauge.eta)
        t_eta = gauge.time[tind]
    end

    max_vel = NaN
    t_vel = NaN
    if !isempty(gauge.u)
        vel = sqrt.( (gauge.u).^2 + (gauge.v).^2 )
        vel[isnan.(vel)] .= 0.0
        max_vel, tind = findmax(vel)
        t_vel = gauge.time[tind]
    end

    gmax = VisClaw.Gaugemax(gauge.label, gauge.id, gauge.loc,
                            max_AMRlevel, max_eta, max_vel,
                            t_eta, t_vel,
                            gauge.unittime)
    return  gmax
end


text:
module variable
include("./helper.jl")

import .helper: ensure_float

mutable struct Variable

    value::Float64
    requires_grad::Bool
    parent_nodes::Any
    grad::Float64

    function Variable(value, requires_grad=False, parent_nodes=Any[], grad=0)
        value = ensure_float(value)

        return new(value, requires_grad, parent_nodes, grad)
    end

    function Base.:+(v1::Variable, v2::Variable)
        return add(v1, v2)
    end

    function Base.:*(v1::Variable, v2::Variable)
        return mul(v1, v2)
    end

end

function add(v1::Variable, v2::Variable)
    requires_grad = v1.requires_grad || v2.requires_grad
    result = v1.value + v2.value

    parent_nodes = Any[]
    if v1.requires_grad
        push!(parent_nodes, Node(v1, 1))
    end

    if v2.requires_grad
        push!(parent_nodes, Node(v2, 1))
    end

    return Variable(result, requires_grad, parent_nodes)
end

function mul(v1::Variable, v2::Variable)
    requires_grad = v1.requires_grad || v2.requires_grad
    result = v1.value * v2.value

    parent_nodes = Any[]
    if v1.requires_grad
        push!(parent_nodes, Node(v1, v2.value))
    end

    if v2.requires_grad
        push!(parent_nodes, Node(v2, v1.value))
    end

    return Variable(result, requires_grad, parent_nodes)
end

mutable struct Node
    # Node stores the dependencies of a
    # variable. This is essential during
    # the backpropagation.

    variable::Variable
    grad::Float64

    function Node(variable::Variable, grad::Any)

        grad = ensure_float(grad)
        return new(variable, grad)
    end
end

export Variable, Node

end

max_stars_repo_name: brymer-meneses/Autograd

max_issues_repo_name: brymer-meneses/Autograd

max_forks_repo_name: brymer-meneses/Autograd


text:
# Container data structure (structure-of-arrays style) for DG elements
# TODO: Taal refactor, remove NVARS, POLYDEG?
mutable struct ElementContainer2D{RealT<:Real, uEltype<:Real, NVARS, POLYDEG} <: AbstractContainer
  inverse_jacobian::Vector{RealT}        # [elements]
  node_coordinates::Array{RealT, 4}      # [orientation, i, j, elements]
  surface_flux_values::Array{uEltype, 4} # [variables, i, direction, elements]
  cell_ids::Vector{Int}                  # [elements]
  # internal `resize!`able storage
  _node_coordinates::Vector{RealT}
  _surface_flux_values::Vector{uEltype}
end

nvariables(::ElementContainer2D{RealT, uEltype, NVARS, POLYDEG}) where {RealT, uEltype, NVARS, POLYDEG} = NVARS
polydeg(::ElementContainer2D{RealT, uEltype, NVARS, POLYDEG}) where {RealT, uEltype, NVARS, POLYDEG} = POLYDEG
Base.eltype(::ElementContainer2D{RealT, uEltype}) where {RealT, uEltype} = uEltype

# Only one-dimensional `Array`s are `resize!`able in Julia.
# Hence, we use `Vector`s as internal storage and `resize!`
# them whenever needed. Then, we reuse the same memory by
# `unsafe_wrap`ping multi-dimensional `Array`s around the
# internal storage.
function Base.resize!(elements::ElementContainer2D, capacity)
  n_nodes = polydeg(elements) + 1
  nvars = nvariables(elements)
  @unpack _node_coordinates, _surface_flux_values,
          inverse_jacobian, cell_ids = elements

  resize!(inverse_jacobian, capacity)

  resize!(_node_coordinates, 2 * n_nodes * n_nodes * capacity)
  elements.node_coordinates = unsafe_wrap(Array, pointer(_node_coordinates),
                                          (2, n_nodes, n_nodes, capacity))

  resize!(_surface_flux_values, nvars * n_nodes * 2 * 2 * capacity)
  elements.surface_flux_values = unsafe_wrap(Array, pointer(_surface_flux_values),
                                             (nvars, n_nodes, 2 * 2, capacity))

  resize!(cell_ids, capacity)

  return nothing
end

function ElementContainer2D{RealT, uEltype, NVARS, POLYDEG}(capacity::Integer) where {RealT<:Real, uEltype<:Real, NVARS, POLYDEG}
  n_nodes = POLYDEG + 1
  nan_RealT = convert(RealT, NaN)
  nan_uEltype = convert(uEltype, NaN)

  # Initialize fields with defaults
  inverse_jacobian = fill(nan_RealT, capacity)

  _node_coordinates = fill(nan_RealT, 2 * n_nodes * n_nodes * capacity)
  node_coordinates = unsafe_wrap(Array, pointer(_node_coordinates),
                                 (2, n_nodes, n_nodes, capacity))

  _surface_flux_values = fill(nan_uEltype, NVARS * n_nodes * 2 * 2 * capacity)
  surface_flux_values = unsafe_wrap(Array, pointer(_surface_flux_values),
                                    (NVARS, n_nodes, 2 * 2, capacity))

  cell_ids = fill(typemin(Int), capacity)

  return ElementContainer2D{RealT, uEltype, NVARS, POLYDEG}(
    inverse_jacobian, node_coordinates, surface_flux_values, cell_ids,
    _node_coordinates, _surface_flux_values)
end

# Return number of elements
@inline nelements(elements::ElementContainer2D) = length(elements.cell_ids)
# TODO: Taal performance, 1:nelements(elements) vs. Base.OneTo(nelements(elements))
@inline eachelement(elements::ElementContainer2D) = Base.OneTo(nelements(elements))

# Create element container and initialize element data
function init_elements(cell_ids, mesh::TreeMesh2D,
                       equations::AbstractEquations{2, NVARS},
                       basis::LobattoLegendreBasis{T, NNODES}, ::Type{RealT}, ::Type{uEltype}) where {RealT<:Real, uEltype<:Real, NVARS, T, NNODES}
  # Initialize container
  n_elements = length(cell_ids)
  elements = ElementContainer2D{RealT, uEltype, NVARS, NNODES-1}(n_elements)

  init_elements!(elements, cell_ids, mesh, basis.nodes)
  return elements
end

function init_elements!(elements, cell_ids, mesh::TreeMesh2D, nodes)
  n_nodes = length(nodes)

  # Store cell ids
  elements.cell_ids .= cell_ids

  # Calculate inverse Jacobian and node coordinates
  for element in eachelement(elements)
    # Get cell id
    cell_id = cell_ids[element]

    # Get cell length
    dx = length_at_cell(mesh.tree, cell_id)

    # Calculate inverse Jacobian as 1/(h/2)
    elements.inverse_jacobian[element] = 2/dx

    # Calculate node coordinates
    for j in 1:n_nodes, i in 1:n_nodes
      elements.node_coordinates[1, i, j, element] = (
          mesh.tree.coordinates[1, cell_id] + dx/2 * nodes[i])
      elements.node_coordinates[2, i, j, element] = (
          mesh.tree.coordinates[2, cell_id] + dx/2 * nodes[j])
    end
  end

  return elements
end

# Container data structure (structure-of-arrays style) for DG interfaces
# TODO: Taal refactor, remove NVARS, POLYDEG?
mutable struct InterfaceContainer2D{uEltype<:Real, NVARS, POLYDEG} <: AbstractContainer
  u::Array{uEltype, 4}      # [leftright, variables, i, interfaces]
  neighbor_ids::Matrix{Int} # [leftright, interfaces]
  orientations::Vector{Int} # [interfaces]
  # internal `resize!`able storage
  _u::Vector{uEltype}
  _neighbor_ids::Vector{Int}
end

nvariables(::InterfaceContainer2D{uEltype, NVARS, POLYDEG}) where {uEltype, NVARS, POLYDEG} = NVARS
polydeg(::InterfaceContainer2D{uEltype, NVARS, POLYDEG}) where {uEltype, NVARS, POLYDEG} = POLYDEG
Base.eltype(::InterfaceContainer2D{uEltype}) where {uEltype} = uEltype

# See explanation of Base.resize! for the element container
function Base.resize!(interfaces::InterfaceContainer2D, capacity)
  n_nodes = polydeg(interfaces) + 1
  nvars = nvariables(interfaces)
  @unpack _u, _neighbor_ids, orientations = interfaces

  resize!(_u, 2 * nvars * n_nodes * capacity)
  interfaces.u = unsafe_wrap(Array, pointer(_u),
                             (2, nvars, n_nodes, capacity))

  resize!(_neighbor_ids, 2 * capacity)
  interfaces.neighbor_ids = unsafe_wrap(Array, pointer(_neighbor_ids ),
                                        (2, capacity))

  resize!(orientations, capacity)

  return nothing
end

function InterfaceContainer2D{uEltype, NVARS, POLYDEG}(capacity::Integer) where {uEltype<:Real, NVARS, POLYDEG}
  n_nodes = POLYDEG + 1
  nan = convert(uEltype, NaN)

  # Initialize fields with defaults
  _u = fill(nan, 2 * NVARS * n_nodes * capacity)
  u = unsafe_wrap(Array, pointer(_u),
                  (2, NVARS, n_nodes, capacity))

  _neighbor_ids = fill(typemin(Int), 2 * capacity)
  neighbor_ids = unsafe_wrap(Array, pointer(_neighbor_ids ),
                             (2, capacity))

  orientations = fill(typemin(Int), capacity)

  return InterfaceContainer2D{uEltype, NVARS, POLYDEG}(
    u, neighbor_ids, orientations,
    _u, _neighbor_ids)
end

# Return number of interfaces
@inline ninterfaces(interfaces::InterfaceContainer2D) = length(interfaces.orientations)

# Create interface container and initialize interface data in `elements`.
function init_interfaces(cell_ids, mesh::TreeMesh2D,
                         elements::ElementContainer2D{RealT, uEltype, NVARS, POLYDEG}) where {RealT<:Real, uEltype<:Real, NVARS, POLYDEG}
  # Initialize container
  n_interfaces = count_required_interfaces(mesh, cell_ids)
  interfaces = InterfaceContainer2D{uEltype, NVARS, POLYDEG}(n_interfaces)

  # Connect elements with interfaces
  init_interfaces!(interfaces, elements, mesh)
  return interfaces
end

# Count the number of interfaces that need to be created
function count_required_interfaces(mesh::TreeMesh2D, cell_ids)
  count = 0

  # Iterate over all cells
  for cell_id in cell_ids
    for direction in eachdirection(mesh.tree)
      # Only count interfaces in positive direction to avoid double counting
      if direction % 2 == 1
        continue
      end

      # If no neighbor exists, current cell is small or at boundary and thus we need a mortar
      if !has_neighbor(mesh.tree, cell_id, direction)
        continue
      end

      # Skip if neighbor has children
      neighbor_cell_id = mesh.tree.neighbor_ids[direction, cell_id]
      if has_children(mesh.tree, neighbor_cell_id)
        continue
      end

      # Skip if neighbor is on different rank -> create MPI interface instead
      if mpi_isparallel() && !is_own_cell(mesh.tree, neighbor_cell_id)
        continue
      end

      count += 1
    end
  end

  return count
end

# Initialize connectivity between elements and interfaces
function init_interfaces!(interfaces, elements, mesh::TreeMesh2D)
  # Construct cell -> element mapping for easier algorithm implementation
  tree = mesh.tree
  c2e = zeros(Int, length(tree))
  for element in eachelement(elements)
    c2e[elements.cell_ids[element]] = element
  end

  # Reset interface count
  count = 0

  # Iterate over all elements to find neighbors and to connect via interfaces
  for element in eachelement(elements)
    # Get cell id
    cell_id = elements.cell_ids[element]

    # Loop over directions
    for direction in eachdirection(mesh.tree)
      # Only create interfaces in positive direction
      if direction % 2 == 1
        continue
      end

      # If no neighbor exists, current cell is small and thus we need a mortar
      if !has_neighbor(mesh.tree, cell_id, direction)
        continue
      end

      # Skip if neighbor has children
      neighbor_cell_id = mesh.tree.neighbor_ids[direction, cell_id]
      if has_children(mesh.tree, neighbor_cell_id)
        continue
      end

      # Skip if neighbor is on different rank -> create MPI interface instead
      if mpi_isparallel() && !is_own_cell(mesh.tree, neighbor_cell_id)
        continue
      end

      # Create interface between elements (1 -> "left" of interface, 2 -> "right" of interface)
      count += 1
      interfaces.neighbor_ids[2, count] = c2e[neighbor_cell_id]
      interfaces.neighbor_ids[1, count] = element

      # Set orientation (x -> 1, y -> 2)
      interfaces.orientations[count] = div(direction, 2)
    end
  end

  @assert count == ninterfaces(interfaces) ("Actual interface count ($count) does not match " *
                                            "expectations $(ninterfaces(interfaces))")
end

# Container data structure (structure-of-arrays style) for DG boundaries
# TODO: Taal refactor, remove NVARS, POLYDEG?
mutable struct BoundaryContainer2D{RealT<:Real, uEltype<:Real, NVARS, POLYDEG} <: AbstractContainer
  u::Array{uEltype, 4}              # [leftright, variables, i, boundaries]
  neighbor_ids::Vector{Int}         # [boundaries]
  orientations::Vector{Int}         # [boundaries]
  neighbor_sides::Vector{Int}       # [boundaries]
  node_coordinates::Array{RealT, 3} # [orientation, i, elements]
  n_boundaries_per_direction::SVector{4, Int} # [direction]
  # internal `resize!`able storage
  _u::Vector{uEltype}
  _node_coordinates::Vector{RealT}
end

nvariables(::BoundaryContainer2D{RealT, uEltype, NVARS, POLYDEG}) where {RealT, uEltype, NVARS, POLYDEG} = NVARS
polydeg(::BoundaryContainer2D{RealT, uEltype, NVARS, POLYDEG}) where {RealT, uEltype, NVARS, POLYDEG} = POLYDEG
Base.eltype(::BoundaryContainer2D{RealT, uEltype}) where {RealT, uEltype} = uEltype

# See explanation of Base.resize! for the element container
function Base.resize!(boundaries::BoundaryContainer2D, capacity)
  n_nodes = polydeg(boundaries) + 1
  nvars = nvariables(boundaries)
  @unpack _u, _node_coordinates,
          neighbor_ids, orientations, neighbor_sides = boundaries

  resize!(_u, 2 * nvars * n_nodes * capacity)
  boundaries.u = unsafe_wrap(Array, pointer(_u),
                             (2, nvars, n_nodes, capacity))

  resize!(_node_coordinates, 2 * n_nodes * capacity)
  boundaries.node_coordinates = unsafe_wrap(Array, pointer(_node_coordinates),
                                            (2, n_nodes, capacity))

  resize!(neighbor_ids, capacity)

  resize!(orientations, capacity)

  resize!(neighbor_sides, capacity)

  return nothing
end

function BoundaryContainer2D{RealT, uEltype, NVARS, POLYDEG}(capacity::Integer) where {RealT<:Real, uEltype<:Real, NVARS, POLYDEG}
  n_nodes = POLYDEG + 1
  nan_RealT = convert(RealT, NaN)
  nan_uEltype = convert(uEltype, NaN)

  # Initialize fields with defaults
  _u = fill(nan_uEltype, 2 * NVARS * n_nodes * capacity)
  u = unsafe_wrap(Array, pointer(_u),
                  (2, NVARS, n_nodes, capacity))

  neighbor_ids = fill(typemin(Int), capacity)

  orientations = fill(typemin(Int), capacity)

  neighbor_sides = fill(typemin(Int), capacity)

  _node_coordinates = fill(nan_RealT, 2 * n_nodes * capacity)
  node_coordinates = unsafe_wrap(Array, pointer(_node_coordinates),
                                 (2, n_nodes, capacity))

  n_boundaries_per_direction = SVector(0, 0, 0, 0)

  return BoundaryContainer2D{RealT, uEltype, NVARS, POLYDEG}(
    u, neighbor_ids, orientations, neighbor_sides,
    node_coordinates, n_boundaries_per_direction,
    _u, _node_coordinates)
end

# Return number of boundaries
@inline nboundaries(boundaries::BoundaryContainer2D) = length(boundaries.orientations)

# Create boundaries container and initialize boundary data in `elements`.
function init_boundaries(cell_ids, mesh::TreeMesh2D,
                         elements::ElementContainer2D{RealT, uEltype, NVARS, POLYDEG}) where {RealT<:Real, uEltype<:Real, NVARS, POLYDEG}
  # Initialize container
  n_boundaries = count_required_boundaries(mesh, cell_ids)
  boundaries = BoundaryContainer2D{RealT, uEltype, NVARS, POLYDEG}(n_boundaries)

  # Connect elements with boundaries
  init_boundaries!(boundaries, elements, mesh)
  return boundaries
end

# Count the number of boundaries that need to be created
function count_required_boundaries(mesh::TreeMesh2D, cell_ids)
  count = 0

  # Iterate over all cells
  for cell_id in cell_ids
    for direction in eachdirection(mesh.tree)
      # If neighbor exists, current cell is not at a boundary
      if has_neighbor(mesh.tree, cell_id, direction)
        continue
      end

      # If coarse neighbor exists, current cell is not at a boundary
      if has_coarse_neighbor(mesh.tree, cell_id, direction)
        continue
      end

      # No neighbor exists in this direction -> must be a boundary
      count += 1
    end
  end

  return count
end

# Initialize connectivity between elements and boundaries
function init_boundaries!(boundaries, elements, mesh::TreeMesh2D)
  # Reset boundaries count
  count = 0

  # Initialize boundary counts
  counts_per_direction = MVector(0, 0, 0, 0)

  # OBS! Iterate over directions first, then over elements, and count boundaries in each direction
  # Rationale: This way the boundaries are internally sorted by the directions -x, +x, -y etc.,
  #            obviating the need to store the boundary condition to be applied explicitly.
  # Loop over directions
  for direction in eachdirection(mesh.tree)
    # Iterate over all elements to find missing neighbors and to connect to boundaries
    for element in eachelement(elements)
      # Get cell id
      cell_id = elements.cell_ids[element]

      # If neighbor exists, current cell is not at a boundary
      if has_neighbor(mesh.tree, cell_id, direction)
        continue
      end

      # If coarse neighbor exists, current cell is not at a boundary
      if has_coarse_neighbor(mesh.tree, cell_id, direction)
        continue
      end

      # Create boundary
      count += 1
      counts_per_direction[direction] += 1

      # Set neighbor element id
      boundaries.neighbor_ids[count] = element

      # Set neighbor side, which denotes the direction (1 -> negative, 2 -> positive) of the element
      if direction in (2, 4)
        boundaries.neighbor_sides[count] = 1
      else
        boundaries.neighbor_sides[count] = 2
      end

      # Set orientation (x -> 1, y -> 2)
      if direction in (1, 2)
        boundaries.orientations[count] = 1
      else
        boundaries.orientations[count] = 2
      end

      # Store node coordinates
      enc = elements.node_coordinates
      if direction == 1 # -x direction
        boundaries.node_coordinates[:, :, count] .= enc[:, 1,   :,   element]
      elseif direction == 2 # +x direction
        boundaries.node_coordinates[:, :, count] .= enc[:, end, :,   element]
      elseif direction == 3 # -y direction
        boundaries.node_coordinates[:, :, count] .= enc[:, :,   1,   element]
      elseif direction == 4 # +y direction
        boundaries.node_coordinates[:, :, count] .= enc[:, :,   end, element]
      else
        error("should not happen")
      end
    end
  end

  @assert count == nboundaries(boundaries) ("Actual boundaries count ($count) does not match " *
                                            "expectations $(nboundaries(boundaries))")
  @assert sum(counts_per_direction) == count

  boundaries.n_boundaries_per_direction = SVector(counts_per_direction)

  return boundaries.n_boundaries_per_direction
end

# Container data structure (structure-of-arrays style) for DG L2 mortars
# Positions/directions for orientations = 1, large_sides = 2:
# mortar is orthogonal to x-axis, large side is in positive coordinate direction wrt mortar
#           |    |
# upper = 2 |    |
#           |    |
#                | 3 = large side
#           |    |
# lower = 1 |    |
#           |    |
# TODO: Taal refactor, remove NVARS, POLYDEG?
# TODO: Taal refactor, mutable struct or resize! for AMR?
mutable struct L2MortarContainer2D{uEltype<:Real, NVARS, POLYDEG} <: AbstractContainer
  u_upper::Array{uEltype, 4}  # [leftright, variables, i, mortars]
  u_lower::Array{uEltype, 4}  # [leftright, variables, i, mortars]
  neighbor_ids::Matrix{Int}   # [position, mortars]
  # Large sides: left -> 1, right -> 2
  large_sides::Vector{Int}  # [mortars]
  orientations::Vector{Int} # [mortars]
  # internal `resize!`able storage
  _u_upper::Vector{uEltype}
  _u_lower::Vector{uEltype}
  _neighbor_ids::Vector{Int}
end

nvariables(::L2MortarContainer2D{uEltype, NVARS, POLYDEG}) where {uEltype, NVARS, POLYDEG} = NVARS
polydeg(::L2MortarContainer2D{uEltype, NVARS, POLYDEG}) where {uEltype, NVARS, POLYDEG} = POLYDEG
Base.eltype(::L2MortarContainer2D{uEltype}) where {uEltype} = uEltype

# See explanation of Base.resize! for the element container
function Base.resize!(mortars::L2MortarContainer2D, capacity)
  n_nodes = polydeg(mortars) + 1
  nvars = nvariables(mortars)
  @unpack _u_upper, _u_lower, _neighbor_ids,
          large_sides, orientations = mortars

  resize!(_u_upper, 2 * nvars * n_nodes * capacity)
  mortars.u_upper = unsafe_wrap(Array, pointer(_u_upper),
                                (2, nvars, n_nodes, capacity))

  resize!(_u_lower, 2 * nvars * n_nodes * capacity)
  mortars.u_lower = unsafe_wrap(Array, pointer(_u_lower),
                                (2, nvars, n_nodes, capacity))

  resize!(_neighbor_ids, 3 * capacity)
  mortars.neighbor_ids = unsafe_wrap(Array, pointer(_neighbor_ids ),
                                        (3, capacity))

  resize!(large_sides, capacity)

  resize!(orientations, capacity)

  return nothing
end

function L2MortarContainer2D{uEltype, NVARS, POLYDEG}(capacity::Integer) where {uEltype<:Real, NVARS, POLYDEG}
  n_nodes = POLYDEG + 1
  nan = convert(uEltype, NaN)

  # Initialize fields with defaults
  _u_upper = fill(nan, 2 * NVARS * n_nodes * capacity)
  u_upper = unsafe_wrap(Array, pointer(_u_upper),
                        (2, NVARS, n_nodes, capacity))

  _u_lower = fill(nan, 2 * NVARS * n_nodes * capacity)
  u_lower = unsafe_wrap(Array, pointer(_u_lower),
                        (2, NVARS, n_nodes, capacity))

  _neighbor_ids = fill(typemin(Int), 3 * capacity)
  neighbor_ids = unsafe_wrap(Array, pointer(_neighbor_ids),
                             (3, capacity))

  large_sides  = fill(typemin(Int), capacity)

  orientations = fill(typemin(Int), capacity)

  return L2MortarContainer2D{uEltype, NVARS, POLYDEG}(
    u_upper, u_lower, neighbor_ids, large_sides, orientations,
    _u_upper, _u_lower, _neighbor_ids)
end

# Return number of L2 mortars
@inline nmortars(l2mortars::L2MortarContainer2D) = length(l2mortars.orientations)

# Allow printing container contents
function Base.show(io::IO, ::MIME"text/plain", c::L2MortarContainer2D)
  @nospecialize c # reduce precompilation time

  println(io, '*'^20)
  for idx in CartesianIndices(c.u_upper)
    println(io, "c.u_upper[$idx] = $(c.u_upper[idx])")
  end
  for idx in CartesianIndices(c.u_lower)
    println(io, "c.u_lower[$idx] = $(c.u_lower[idx])")
  end
  println(io, "transpose(c.neighbor_ids) = $(transpose(c.neighbor_ids))")
  println(io, "c.large_sides = $(c.large_sides)")
  println(io, "c.orientations = $(c.orientations)")
  print(io,   '*'^20)
end

# Create mortar container and initialize mortar data in `elements`.
function init_mortars(cell_ids, mesh::TreeMesh2D,
                      elements::ElementContainer2D{RealT, uEltype, NVARS, POLYDEG},
                      ::LobattoLegendreMortarL2) where {RealT<:Real, uEltype<:Real, NVARS, POLYDEG}
  # Initialize containers
  n_mortars = count_required_mortars(mesh, cell_ids)
  mortars = L2MortarContainer2D{uEltype, NVARS, POLYDEG}(n_mortars)

  # Connect elements with mortars
  init_mortars!(mortars, elements, mesh)
  return mortars
end

# Count the number of mortars that need to be created
function count_required_mortars(mesh::TreeMesh2D, cell_ids)
  count = 0

  # Iterate over all cells and count mortars from perspective of coarse cells
  for cell_id in cell_ids
    for direction in eachdirection(mesh.tree)
      # If no neighbor exists, cell is small with large neighbor or at boundary -> do nothing
      if !has_neighbor(mesh.tree, cell_id, direction)
        continue
      end

      # If neighbor has no children, this is a conforming interface -> do nothing
      neighbor_id = mesh.tree.neighbor_ids[direction, cell_id]
      if !has_children(mesh.tree, neighbor_id)
        continue
      end

      count +=1
    end
  end

  return count
end

# Initialize connectivity between elements and mortars
function init_mortars!(mortars, elements, mesh::TreeMesh2D)
  # Construct cell -> element mapping for easier algorithm implementation
  tree = mesh.tree
  c2e = zeros(Int, length(tree))
  for element in eachelement(elements)
    c2e[elements.cell_ids[element]] = element
  end

  # Reset interface count
  count = 0

  # Iterate over all elements to find neighbors and to connect via interfaces
  for element in eachelement(elements)
    # Get cell id
    cell_id = elements.cell_ids[element]

    for direction in eachdirection(mesh.tree)
      # If no neighbor exists, cell is small with large neighbor -> do nothing
      if !has_neighbor(mesh.tree, cell_id, direction)
        continue
      end

      # If neighbor has no children, this is a conforming interface -> do nothing
      neighbor_cell_id = mesh.tree.neighbor_ids[direction, cell_id]
      if !has_children(mesh.tree, neighbor_cell_id)
        continue
      end

      # Create mortar between elements:
      # 1 -> small element in negative coordinate direction
      # 2 -> small element in positive coordinate direction
      # 3 -> large element
      count += 1
      mortars.neighbor_ids[3, count] = element
      if direction == 1
        mortars.neighbor_ids[1, count] = c2e[mesh.tree.child_ids[2, neighbor_cell_id]]
        mortars.neighbor_ids[2, count] = c2e[mesh.tree.child_ids[4, neighbor_cell_id]]
      elseif direction == 2
        mortars.neighbor_ids[1, count] = c2e[mesh.tree.child_ids[1, neighbor_cell_id]]
        mortars.neighbor_ids[2, count] = c2e[mesh.tree.child_ids[3, neighbor_cell_id]]
      elseif direction == 3
        mortars.neighbor_ids[1, count] = c2e[mesh.tree.child_ids[3, neighbor_cell_id]]
        mortars.neighbor_ids[2, count] = c2e[mesh.tree.child_ids[4, neighbor_cell_id]]
      elseif direction == 4
        mortars.neighbor_ids[1, count] = c2e[mesh.tree.child_ids[1, neighbor_cell_id]]
        mortars.neighbor_ids[2, count] = c2e[mesh.tree.child_ids[2, neighbor_cell_id]]
      else
        error("should not happen")
      end

      # Set large side, which denotes the direction (1 -> negative, 2 -> positive) of the large side
      if direction in [2, 4]
        mortars.large_sides[count] = 1
      else
        mortars.large_sides[count] = 2
      end

      # Set orientation (x -> 1, y -> 2)
      if direction in [1, 2]
        mortars.orientations[count] = 1
      else
        mortars.orientations[count] = 2
      end
    end
  end

  @assert count == nmortars(mortars) ("Actual mortar count ($count) does not match " *
                                      "expectations $(nmortars(mortars))")
end

# Container data structure (structure-of-arrays style) for DG MPI interfaces
mutable struct MPIInterfaceContainer2D{uEltype<:Real, NVARS, POLYDEG} <: AbstractContainer
  u::Array{uEltype, 4}           # [leftright, variables, i, interfaces]
  local_element_ids::Vector{Int} # [interfaces]
  orientations::Vector{Int}      # [interfaces]
  remote_sides::Vector{Int}      # [interfaces]
  # internal `resize!`able storage
  _u::Vector{uEltype}
end

nvariables(::MPIInterfaceContainer2D{uEltype, NVARS, POLYDEG}) where {uEltype, NVARS, POLYDEG} = NVARS
polydeg(::MPIInterfaceContainer2D{uEltype, NVARS, POLYDEG}) where {uEltype, NVARS, POLYDEG} = POLYDEG
Base.eltype(::MPIInterfaceContainer2D{uEltype}) where {uEltype} = uEltype

# See explanation of Base.resize! for the element container
function Base.resize!(mpi_interfaces::MPIInterfaceContainer2D, capacity)
  n_nodes = polydeg(mpi_interfaces) + 1
  nvars = nvariables(mpi_interfaces)
  @unpack _u, local_element_ids, orientations, remote_sides = mpi_interfaces

  resize!(_u, 2 * nvars * n_nodes * capacity)
  mpi_interfaces.u = unsafe_wrap(Array, pointer(_u),
                                 (2, nvars, n_nodes, capacity))

  resize!(local_element_ids, capacity)

  resize!(orientations, capacity)

  resize!(remote_sides, capacity)

  return nothing
end

function MPIInterfaceContainer2D{uEltype, NVARS, POLYDEG}(capacity::Integer) where {uEltype<:Real, NVARS, POLYDEG}
  n_nodes = POLYDEG + 1
  nan = convert(uEltype, NaN)

  # Initialize fields with defaults
  _u = fill(nan, 2 * NVARS * n_nodes * capacity)
  u = unsafe_wrap(Array, pointer(_u),
                  (2, NVARS, n_nodes, capacity))

  local_element_ids = fill(typemin(Int), capacity)

  orientations = fill(typemin(Int), capacity)

  remote_sides = fill(typemin(Int), capacity)

  return MPIInterfaceContainer2D{uEltype, NVARS, POLYDEG}(
    u, local_element_ids, orientations, remote_sides,
    _u)
end

# TODO: Taal, rename to ninterfaces?
# Return number of interfaces
nmpiinterfaces(mpi_interfaces::MPIInterfaceContainer2D) = length(mpi_interfaces.orientations)

# Create MPI interface container and initialize MPI interface data in `elements`.
function init_mpi_interfaces(cell_ids, mesh::TreeMesh2D,
                             elements::ElementContainer2D{RealT, uEltype, NVARS, POLYDEG}) where {RealT<:Real, uEltype<:Real, NVARS, POLYDEG}
  # Initialize container
  n_mpi_interfaces = count_required_mpi_interfaces(mesh, cell_ids)
  mpi_interfaces = MPIInterfaceContainer2D{uEltype, NVARS, POLYDEG}(n_mpi_interfaces)

  # Connect elements with interfaces
  init_mpi_interfaces!(mpi_interfaces, elements, mesh)
  return mpi_interfaces
end

# Count the number of MPI interfaces that need to be created
function count_required_mpi_interfaces(mesh::TreeMesh2D, cell_ids)
  count = 0

  # Iterate over all cells
  for cell_id in cell_ids
    for direction in eachdirection(mesh.tree)
      # If no neighbor exists, current cell is small or at boundary and thus we need a mortar
      if !has_neighbor(mesh.tree, cell_id, direction)
        continue
      end

      # Skip if neighbor has children
      neighbor_cell_id = mesh.tree.neighbor_ids[direction, cell_id]
      if has_children(mesh.tree, neighbor_cell_id)
        continue
      end

      # Skip if neighbor is on this rank -> create regular interface instead
      if mpi_isparallel() && is_own_cell(mesh.tree, neighbor_cell_id)
        continue
      end

      count += 1
    end
  end

  return count
end

# Initialize connectivity between elements and interfaces
function init_mpi_interfaces!(mpi_interfaces, elements, mesh::TreeMesh2D)
  # Reset interface count
  count = 0

  # Iterate over all elements to find neighbors and to connect via mpi_interfaces
  for element in eachelement(elements)
    # Get cell id
    cell_id = elements.cell_ids[element]

    # Loop over directions
    for direction in eachdirection(mesh.tree)
      # If no neighbor exists, current cell is small and thus we need a mortar
      if !has_neighbor(mesh.tree, cell_id, direction)
        continue
      end

      # Skip if neighbor has children
      neighbor_cell_id = mesh.tree.neighbor_ids[direction, cell_id]
      if has_children(mesh.tree, neighbor_cell_id)
        continue
      end

      # Skip if neighbor is on this MPI rank -> create regular interface instead
      if mpi_isparallel() && is_own_cell(mesh.tree, neighbor_cell_id)
        continue
      end

      # Create interface between elements
      count += 1
      mpi_interfaces.local_element_ids[count] = element

      if direction in (2, 4) # element is "left" of interface, remote cell is "right" of interface
        mpi_interfaces.remote_sides[count] = 2
      else
        mpi_interfaces.remote_sides[count] = 1
      end

      # Set orientation (x -> 1, y -> 2)
      if direction in (1, 2) # x-direction
        mpi_interfaces.orientations[count] = 1
      else # y-direction
        mpi_interfaces.orientations[count] = 2
      end
    end
  end

  @assert count == nmpiinterfaces(mpi_interfaces) ("Actual interface count ($count) does not match "
                                                   * "expectations $(nmpiinterfaces(mpi_interfaces))")
end


text:
module Nexus

export parse_taxon_block

const NTAXA_DIM_START = "Dimensions ntax="

function parse_taxon_block(path::String)

    taxa = String[]

    is_ntaxa_line = false
    current_taxon = 1

    for line in eachline(path)
        stripped_line = strip(line)
        if isempty(line)
            continue
        elseif stripped_line == "Begin taxa;"
            is_ntaxa_line = true
        elseif is_ntaxa_line
            @assert startswith(stripped_line, NTAXA_DIM_START) && endswith(line, ';')
            str_ntaxa = stripped_line[(length(NTAXA_DIM_START) + 1):(end - 1)]
            ntaxa = parse(Int, str_ntaxa)
            taxa = fill("", ntaxa)
            is_ntaxa_line = false
        elseif !isempty(taxa)
            if stripped_line == "Taxlabels"
                continue
            elseif length(taxa) >= current_taxon
                taxa[current_taxon] = stripped_line
                current_taxon += 1
            else # all taxa done
                break;
            end
        end
    end

    return taxa
end

end


text:
module MyGenie

using Genie, Genie.Router, Genie.Renderer, Genie.AppServer

function main()
  Base.eval(Main, :(const UserApp = MyGenie))

  include("genie.jl")

  Base.eval(Main, :(const Genie = MyGenie.Genie))
  Base.eval(Main, :(using Genie))
end

main()

end


text:
"""
    generateidealdata(numberofgridpoints, kind)

numberofgridpoint
kind:
* 'nested' for array with n gridpoint arrays
* 'flat' for array with n gridpoints
""" generateidealdata
function generateidealdata(numberofgridpoints::Integer, kind)
    #   Arg-Check
    n  = numberofgridpoints
    n < 2 && throw(ArgumentError("number of gridpoints has to be greater than 1"))
    #   Gridpoints erstellen
    if kind == "flat"
        gridpoints = [collect(1:n) ones(Int, n)]'
    elseif kind == "nested"
        gridpoints = [[i,1] for i in 1:n]
    elseif kind == "static-nested"
        gridpoints = [SVector{2, Int}(i,1) for i in 1:n]
    else
        throw(ArgumentError("kind has to be either 'flat', 'nested' or 'static-nested'"))
    end

    return gridpoints
end

max_stars_repo_name: Affos/Bamberg

max_issues_repo_name: Affos/Bamberg

max_forks_repo_name: Affos/Bamberg


text:
# This file contains code to convert serialized data from old formats to work with newer
# code.

#
# In version 1.0.1 the DeterministicMetadata struct added field time_series type.
# Deserialization needs to add this field and value.
#

function pre_deserialize_conversion!(raw, sys::System)
    old = raw["data_format_version"]
    if old == "1.0.0"
        for component in raw["data"]["components"]
            for ts_metadata in get(component, "time_series_container", [])
                if ts_metadata["__metadata__"]["type"] == "DeterministicMetadata" &&
                   !haskey(ts_metadata, "time_series_type")
                    # This will allow deserialization to work.
                    # post_deserialize_conversion will fix the type.
                    ts_metadata["time_series_type"] = Dict(
                        "__metadata__" => Dict(
                            "module" => "InfrastructureSystems",
                            "type" => "AbstractDeterministic",
                        ),
                    )
                end
            end
        end
    else
        error("conversion of data from $old to $DATA_FORMAT_VERSION is not supported")
    end
end

function post_deserialize_conversion!(sys::System, raw)
    old = raw["data_format_version"]
    if old == "1.0.0"
        for component in IS.iterate_components_with_time_series(sys.data.components)
            ts_container = IS.get_time_series_container(component)
            for key in keys(ts_container.data)
                if key.time_series_type == IS.DeterministicMetadata
                    ts_metadata = ts_container.data[key]
                    ts = get_time_series(
                        AbstractDeterministic,
                        component,
                        get_name(ts_metadata),
                        len = get_horizon(ts_metadata),
                        count = 1,
                    )
                    ts_metadata.time_series_type = typeof(ts)
                end
            end
        end
    else
        error("conversion of data from $old to $DATA_FORMAT_VERSION is not supported")
    end
end


text:
@testset "Base" begin
    @testset "RefValue" begin
        x = Ref(1)
        p, re = Functors.functor(x)
        @test p == (x = 1,)
        @test re(p) isa Base.RefValue{Int}
    end
end


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

commit_hash = "475890c3a760300f5b088c0c308d2b3b95b2acbb"
# sha256sum of the zip file
sha256_sum = "42e25b9ddd245fc2c09318269dda67c55bba77b0c5304a7f11c13d5715ae7b4f"

# Collection of sources required to build cddlibBuilder
sources = [
    "https://github.com/cddlib/cddlib/archive/$commit_hash.zip" =>
    sha256_sum,
]
name = "cddlib"
version = v"0.94j+"

# Bash recipe for building across all platforms
script = raw"""
cd cddlib-475890c3a760300f5b088c0c308d2b3b95b2acbb
sed -i '/^libcddgmp_la_LDFLAGS/ s/$/  $(CDD_LDFLAGS)/' lib-src/Makefile.gmp.am
./bootstrap # needed because we are building a commit and not a release
CPPFLAGS=-I$prefix/include ./configure --prefix=$prefix --host=$target
cd lib-src # Call make in the lib-src subfolder to avoid building the doc folder since pdflatex is not installed
make -j
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

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

# Dependencies that must be installed before this package can be built
dependencies = [
    "https://github.com/JuliaPackaging/Yggdrasil/releases/download/GMP-v6.1.2-1/build_GMP.v6.1.2.jl"
]

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


text:
# Fit particles to evenly spaced grid. There might be unused spots.
function mdplace!( r::MolecularDynamicsTrial )
    # Make sure this isn't a trial in progress
    if r.currentStep != 0
        error("Cannot place particles: currentStep = ", r.currentStep);
    end

    # We need at least enough grid notches for the number of bodies; round up
    gridNotches = ceil( r.numBodies ^ (1/3) );

    # We need to know how long each grid unit is
    gridUnit = r.L / gridNotches;

    # Iterate through the bodies and lay them down in neat little rows
    println("On step 1\nPlacing particles")
    for i in 1:r.numBodies
        num = i - 1; # Because Julia and MATLAB start at 1, not 0

        # Provide jth coordinate value and shift radix representation right
        for j in 1:3
            grid = ( num % gridNotches );       # find radix^j place
            num = ( num - grid ) / gridNotches; # shift right
            grid += 0.1rand() + 0.45;           # middle of cell, some variation
            r.y[j,i,1] = grid * gridUnit;       # scale by box size
        end
    end

    # Calculate energies for 1st step
    potentialenergyandforce!( r, 1 )
    kineticenergy!( r, 1 )
    r.e[1] = r.ket[1] + r.pet[1]

    # Particles are now placed, so say we are on step 1
    r.currentStep = 1;
end


text:
# # Iter Data Interface

# Types (and instances of those types) that implement the `AbstractIterData`
# interface are used to have a unified access to relevant iteration data in
# all functions of the main algorithm.

# First of all, we make it so that any `AbstractResult` is broadcasted wholly:
Base.broadcastable( id :: AbstractIterData ) = Ref( id )

# ## Mandatory Methods

# Any subtype of `AbstractIterData` should implement these methods.

# First, we need a constructor.
# This constructor should take a base type implementing `AbstractIterData` and
# its first states as arguments, where `x` and `fx` are vectors of floats and `Δ`
# is either a float or a vector thereof:
function _init_iter_data( T :: Type{<:AbstractIterData}, x :: VecF, fx :: VecF, Δ :: NumOrVecF )
    return nothing
end

# ### Getters

# There are Getters for the mathematical objects relevant during optimzation:
"Return current iteration site vector ``xᵗ``."
function get_x( :: AbstractIterData{XT,YT,DT} ) where {XT,YT,DT}
    return XT()
end

"Return current value vector ``f(xᵗ)``."
function get_fx(  :: AbstractIterData{XT,YT,DT} ) where {XT,YT,DT}
    return YT()
end

"Return current trust region radius (vector) ``Δᵗ``."
function get_delta( :: AbstractIterData{XT,YT,DT} ) where {XT,YT,DT}
    return zero(eltype(DT))
end

# We also need the iteration result index for our database.
# This should be implemented but works as is if only `MockDB` is used.
"Index (or `id`) of current iterate in database."
get_x_index( :: AbstractIterData ) :: Int = -1

# For printing information and for stopping we need the iteration count:
"Return number of iterations so far."
get_num_iterations( :: AbstractIterData ) :: Int = 0
"Return the number of model improvement iterations so far."
get_num_model_improvements( :: AbstractIterData ) :: Int = 0

# Finally, we also want to display the last iteration classification:
"Return the iteration classification of `ITER_TYPE`."
it_stat( :: AbstractIterData ) :: ITER_TYPE = SUCCESSFULL

# ### Setters

# Implement the setters and note the leading underscore!
# The algorithm will actually call `set_x!` instead of `_set_x!` etc.
# These derived methods are implemented below and ensure that we actually
# store copies of the input!

"Set current iteration site to `x̂`."
_set_x!( :: AbstractIterData, x̂ :: Vec ) :: Nothing = nothing

"Set current iteration value vector to `x̂`."
_set_fx!( :: AbstractIterData, ŷ :: Vec ) :: Nothing = nothing

"Set current trust region radius (vector?) to `Δ`."
_set_delta!( :: AbstractIterData, Δ :: NumOrVec ) :: Nothing = nothing

"Set the current iteration database id to `val`."
set_x_index!( :: AbstractIterData, val :: Int ) :: Nothing = nothing

# The itaration counters are modified by the following methods:
"Set the iteration counter to `N`."
set_num_iterations!( :: AbstractIterData, N :: Int = 0 ) :: Nothing = nothing

"Set the improvement counter to `N`."
set_num_model_improvements!( :: AbstractIterData, N :: Int = 0 ) :: Nothing = nothing ;

# Of course, we also need a way to set the iteration classification:
"Set the iteration classification."
it_stat!( :: AbstractIterData, :: ITER_TYPE ) :: Nothing = nothing;

# ## Derived Methods

# The actual setters simply ensure proper copying:
set_x!( id :: AbstractIterData, x :: Vec ) = _set_x!(id, copy(x))
set_fx!( id :: AbstractIterData, fx :: Vec ) = _set_fx!(id, copy(fx))
set_delta!( id :: AbstractIterData, Δ :: NumOrVec ) = _set_delta!(id, copy(Δ))

# These two auxillary methods are derived from the above definitions:
"Increase the iteration counter by `N`."
function inc_num_iterations!( id :: AbstractIterData, N :: Int = 1 )
    current_num_iter = get_num_iterations( id )
    return set_num_iterations!( id, current_num_iter + N)
end

"Increase the model improvement counter by `N`."
function inc_num_model_improvements!( id :: AbstractIterData, N :: Int = 1 )
    current_num_model_improvements = get_num_model_improvements( id )
    return set_num_iterations!( id, current_num_model_improvements + N )
end

# The actual constructor ensures that all arguments have eltype<:AbstractFloat:
"""
    init_iter_data( T , x, fx, Δ )
Return an instance of "base" type `T` implementing `AbstractIterData` with
correct type parameters for `x`, `fx` and `Δ`.
`x` and `fx` should be vectors of floats and `Δ` can either be a float or
a vector of floats.
"""
function init_iter_data( T :: Type{<:AbstractIterData}, x :: Vec, fx :: Vec, Δ :: NumOrVec )
    global MIN_PRECISION
    base_type = @eval $(T.name.name)    # strip any type parameters from T
    XTe = Base.promote_eltype( x, MIN_PRECISION )
    YTe = Base.promote_eltype( fx, MIN_PRECISION )
    DTe = Base.promote_eltype( Δ, MIN_PRECISION )
	return _init_iter_data( base_type, XTe.(x), YTe.(y), DTe.(Δ) )
end

# # `AbstractIterSaveable`

# The only thing we require for an `AbstractIterSaveable` is a constructor (think "extractor")
# that gets us a saveable object to store in the database.
# It can have arbritrary keyword arguments but should accept `kwargs...` in any
# case:

get_saveable( :: Type{<:AbstractIterSaveable}, id :: AbstractIterData; kwargs... ) = nothing
get_saveable( :: Type{<:Nothing}, id :: AbstractIterData; kwargs... ) = nothing


text:
module MusicXML

using EzXML

end


text:
using Test

@testset "CombinatorialMaps" begin
  include("CombinatorialMaps.jl")
end

@testset "SimplicialSets" begin
  include("SimplicialSets.jl")
  include("DualSimplicialSets.jl")
end

@testset "Meshes" begin
  include("MeshInterop.jl")
  include("MeshGraphics.jl")
end