text:
using QuasiGeostrophy, CUDA, BenchmarkTools, FFTW
# 1D Tests file: include("sandbox/gpu_tests.jl")
Ω = S¹(0,2π)
Nx = 2^8;
arraytype = CuArray
grid = FourierGrid(Nx, Ω, arraytype = arraytype)
ϕ, ϕ2 = create_fields(names = ("ϕ", "ϕ2"), grid = grid, arraytype = arraytype)
∂x = create_operators(grid, arraytype = arraytype)
# initialize fields with nontrivial data
x = grid.grid[1]
k = grid.wavenumbers[1]

P = ϕ.metadata.transform.forward
c_P = plan_fft(ϕ.data)
scale = round(Int, 1/(x[2]-x[1]))
ν = 10^4
κ = 1e-4
iv = @. tanh( ν * (x - π/2) ) * tanh( ν * (x - 3π/2) )
ϕ(iv)
ϕ2(iv)
Δ =  FourierOperator(∂x.op .* ∂x.op,FourierOperatorMetaData("Δ"))
Δt = 0.01 / scale * 2π
t = randn(1) .* 0

for i in 1:100*scale
    f = ∂x(ϕ) + κ * Δ(ϕ)
    nϕ = ϕ + Δt * f
    g  = ∂x(nϕ)  + κ * Δ(nϕ)
    nϕ = ϕ + Δt/2 * ( f +  g)
    ϕ.data .= nϕ.data
    t[1] += Δt
end
println(norm(ϕ-ϕ2)/norm(ϕ2))

## 2D Tests
# 2D Test
Ω = S¹(0,2π) × S¹(0,2π)
Nx = Ny =  2^7;
grid = FourierGrid((Nx, Ny), Ω, arraytype = arraytype)
fieldnames = ("ϕ", "ϕ2", "u", "v", "ψ")
ϕ, ϕ2, u, v, ψ = create_fields(names = fieldnames, grid = grid, arraytype = arraytype)
∂x, ∂y =  create_operators(grid, arraytype = arraytype)
Δ =  FourierOperator(∂x.op .* ∂x.op .+ ∂y.op .* ∂y.op,FourierOperatorMetaData("Δ"))
# initialize fields with nontrivial data
#  Heuns Method
x, y = grid.grid
scale = round(Int, 1/(x[2]-x[1]))
ν = 10000
κ = 0e-2
iv = @. tanh( ν*(x - π/2)) * tanh(ν*(x - 3π/2) ) * tanh( ν*(y - π/2)) * tanh(ν*(y - 3π/2) )
ϕ(iv)
ϕ2(iv)
Δt = 0.01 / scale * 2π
t = randn(1) .* 0
u( (x .+ y) .* 0 .+1 )
v( (x .+ y) .* 0 .+1 )
ψ(@. sin(x/2) * sin(y/2))
norm(∂x(u) + ∂y(v))
for i in 1:100*scale
    f = u * ∂x(ϕ) + v * ∂y(ϕ) + κ * Δ(ϕ)
    nϕ = ϕ + Δt * f
    nϕ = ϕ + Δt/2 * ( f +  u * ∂x(nϕ) + v * ∂y(nϕ)  + κ * Δ(nϕ))
    ϕ.data .= nϕ.data
    t[1] += Δt
end
for i in 1:100*scale
    f =  -u * ∂x(ϕ) + -v * ∂y(ϕ) + κ * Δ(ϕ)
    nϕ = ϕ + Δt * f
    nϕ = ϕ + Δt/2 * ( f +  -u * ∂x(nϕ) + -v * ∂y(nϕ)  + κ * Δ(nϕ))
    ϕ.data .= nϕ.data
    t[1] += Δt
end
println(norm(ϕ-ϕ2) / norm(ϕ2))

max_stars_repo_name: upiterbarg/QuasiGeostrophy

max_issues_repo_name: upiterbarg/QuasiGeostrophy

max_forks_repo_name: upiterbarg/QuasiGeostrophy


text:
using BenchmarkTools, PoincareInvariants

const SUITE = BenchmarkGroup()

SUITE["FirstPoincareInvariants"] = BenchmarkGroup()

SUITE["SecondPoincareInvariants"] = BenchmarkGroup()

SUITE["SecondPoincareInvariants"]["canonical"] = BenchmarkGroup()
for N in [100 * 100, 1000 * 1000], D in [2, 12, 100]
    name = "SecondPoincareInvariant{Float64}(Ω, $D, $N)"

    SUITE["SecondPoincareInvariants"]["canonical"][name] =
        @benchmarkable compute!(pinv, phasepoints, 0, nothing) setup=begin
            Ω(z, t, p) = CanonicalSymplecticMatrix($D)
            pinv = SecondPoincareInvariant{Float64}(Ω, $D, $N, Val(false))
            phasepoints = rand(getpointnum(pinv), $D)
        end
end

results = run(SUITE)

ext: jl


text:
# Learnings:
# - the ! operator
    # - push with ! operator means it will overwrite list with the new value
    # - push without ! operator means it will return the new value, but not modify list

function computeGroup(rows::Array{String})
    # println("Rows:\n", rows)
    seen = Set{Char}()
    for r in rows
        for c in r
            push!(seen, c)
        end
    end
    return length(seen)
end

# example
exampleGroup = [
"abcx"
"abcy"
"abcz"
]
@assert computeGroup(exampleGroup) == 6

function computeGroupV2(rows::Array{String})
    # println("Rows:\n", rows)
    sets = []
    for r in rows
        s = Set{Char}()
        for c in r
            push!(s, c)
        end
        push!(sets, s)
    end
    common = intersect(sets...)
    return length(common)
end

function parseLines(lines)::Array{Array{String}}
    groups = []
    g = []
    for l in lines
        if l == ""
            # reset
            push!(groups, g)
            g = []
        else
            push!(g, l)
        end
    end

    # add last group
    push!(groups, g)

    return groups
end

function run(fname, version)
    # read file
    f = open(fname, "r")
    lines = readlines(f)
    close(f)

    groups = parseLines(lines)
    f = version == 1 ? computeGroup : computeGroupV2
    counts = map(f, groups)
    result = sum(counts)
    println("Ran ", fname, " got ", result)
    return result
end

@assert run("6ex.txt", 1) == 11
run("6.txt", 1)

@assert run("6ex.txt", 2) == 6
run("6.txt", 2)


text:
# This build includes components from LLDBuilder and WABTBuilder

using BinaryProvider # requires BinaryProvider 0.3.0 or later

# Parse some basic command-line arguments
const verbose = "--verbose" in ARGS
const prefix = Prefix(get([a for a in ARGS if a != "--verbose"], 1, joinpath(@__DIR__, "usr")))
products = [
    LibraryProduct(prefix, String["libLLVM"], :libLLVM),
    ExecutableProduct(prefix, "llc", :llc),
    ExecutableProduct(prefix, "lld", :lld),
    ExecutableProduct(prefix, "wast2wasm", :wast2wasm),
    ExecutableProduct(prefix, "wasm2wast", :wasm2wast),
]

# Download binaries from hosted location
bin_prefix = "https://github.com/tshort/LLDBuilder/releases/download/v6.0.0-3"

# Listing of files generated by BinaryBuilder:
download_info = Dict(
    Linux(:aarch64, :glibc) => ("$bin_prefix/lld.aarch64-linux-gnu.tar.gz", "539650c04e9acb9e78568f7753da63cdee6692ae796196eb9e409d733de721c1"),
    Linux(:armv7l, :glibc, :eabihf) => ("$bin_prefix/lld.arm-linux-gnueabihf.tar.gz", "d4ef33968c21e3987a08642c93cca020ce60745b7d99cd280cf5045006b4a6e0"),
    Linux(:i686, :glibc) => ("$bin_prefix/lld.i686-linux-gnu.tar.gz", "11d835bcdca411d457a35e3018801ed338492947f44d3af964b22f13646a5ada"),
    Windows(:i686) => ("$bin_prefix/lld.i686-w64-mingw32.tar.gz", "36f4d6cf50118238746f440d22da3fcc9e63635a10453138b5ed5688c55142bb"),
    MacOS(:x86_64) => ("$bin_prefix/lld.x86_64-apple-darwin14.tar.gz", "d34f798a2a70295a1be115ffc1d19db9f1d947492137a274461b7dee49710e29"),
    Linux(:x86_64, :glibc) => ("$bin_prefix/lld.x86_64-linux-gnu.tar.gz", "9036170580a90654740bc0058e88f5d7a8e4ed353c477beeb3301dcfee6a5d70"),
    Windows(:x86_64) => ("$bin_prefix/lld.x86_64-w64-mingw32.tar.gz", "1e3c17bde002a072013e5e14fad9597f1155456aa8b79459c5dfc3ddc5ba21b0"),
)
bin_prefix = "https://github.com/tshort/WABTBuilder/releases/download/wabt1.0.0-0"
download_info2 = Dict(
    BinaryProvider.Linux(:aarch64, :glibc, :blank_abi) => ("$bin_prefix/WABT.aarch64-linux-gnu.tar.gz", "d979e0b71799da6fe64dff2e31743ae1ce37aca4a474f6cbd6d1e258dfc69966"),
    BinaryProvider.Linux(:aarch64, :musl, :blank_abi) => ("$bin_prefix/WABT.aarch64-linux-musl.tar.gz", "447d13319b4e869f07e68ac225176e859cc748a0f8f2d977c189809b0571ae99"),
    BinaryProvider.Linux(:armv7l, :glibc, :eabihf) => ("$bin_prefix/WABT.arm-linux-gnueabihf.tar.gz", "d4fb9d6c6f63b69b352fab1a750feec51a38f6d76787af29e5887a92b2278a9f"),
    BinaryProvider.Linux(:armv7l, :musl, :eabihf) => ("$bin_prefix/WABT.arm-linux-musleabihf.tar.gz", "c6fac837327a2ac64f578b9d43163019058f45caf3181cc7307325b19e7147e7"),
    BinaryProvider.Linux(:i686, :glibc, :blank_abi) => ("$bin_prefix/WABT.i686-linux-gnu.tar.gz", "06f246047e4ac8a534d515d6deeb0f84141a5e09ad16ad238fb0312cd9559e61"),
    BinaryProvider.Linux(:i686, :musl, :blank_abi) => ("$bin_prefix/WABT.i686-linux-musl.tar.gz", "319e87b6ed1386f13f5d14e7c6cd5e1ea1cc1f5f926a43b3502578736e80a0bc"),
    BinaryProvider.Windows(:i686, :blank_libc, :blank_abi) => ("$bin_prefix/WABT.i686-w64-mingw32.tar.gz", "b10d7b7abad86da5222da5e46d0b0d868423d883820bd7b9e1bd0b01c574eb7e"),
    BinaryProvider.Linux(:powerpc64le, :glibc, :blank_abi) => ("$bin_prefix/WABT.powerpc64le-linux-gnu.tar.gz", "70e5b51b41d85e041379ea28259e12e45d129c36bb490ce5b007eda866ff13af"),
    BinaryProvider.MacOS(:x86_64, :blank_libc, :blank_abi) => ("$bin_prefix/WABT.x86_64-apple-darwin14.tar.gz", "9f437ff428c777f24687f2408f5ca5faea18da8ef0955637164b691f0e44be5f"),
    BinaryProvider.Linux(:x86_64, :glibc, :blank_abi) => ("$bin_prefix/WABT.x86_64-linux-gnu.tar.gz", "ee24ce9f5f9b732fcba5e8ab3311a318579846ebd3f5961a0c0812c0432a5f7a"),
    BinaryProvider.Linux(:x86_64, :musl, :blank_abi) => ("$bin_prefix/WABT.x86_64-linux-musl.tar.gz", "0d9ef975f70abc997d98785b067940bcf7d483b2705298513adcba844c6f63aa"),
    BinaryProvider.Windows(:x86_64, :blank_libc, :blank_abi) => ("$bin_prefix/WABT.x86_64-w64-mingw32.tar.gz", "0bc40a171839d4f73f85bdfdc526f8e5f41262b58c65094ad499466a0f313cd8"),
)

# Install unsatisfied or updated dependencies:
unsatisfied = any(!satisfied(p; verbose=verbose) for p in products)
if haskey(download_info, platform_key())
    url, tarball_hash = download_info[platform_key()]
    if unsatisfied || !isinstalled(url, tarball_hash; prefix=prefix)
        # Download and install binaries
        install(url, tarball_hash; prefix=prefix, force=true, verbose=verbose)
    end
    url, tarball_hash = download_info2[platform_key()]
    if unsatisfied || !isinstalled(url, tarball_hash; prefix=prefix)
        # Download and install binaries
        install(url, tarball_hash; prefix=prefix, force=true, verbose=verbose)
    end
elseif unsatisfied
    # If we don't have a BinaryProvider-compatible .tar.gz to download, complain.
    # Alternatively, you could attempt to install from a separate provider,
    # build from source or something even more ambitious here.
    error("Your platform $(triplet(platform_key())) is not supported by this package!")
end

# Write out a deps.jl file that will contain mappings for our products
write_deps_file(joinpath(@__DIR__, "deps.jl"), products)


text:
# is for scheduling robot reference motion (not a must at the moment (no channels))
function scheduleReferenceMotion(ref::ReferencePath,
                                 referenceMotion!::Function)
    error("von scheduleReferenceMotion Zeile 4")
    referenceMotion!(ref)
end

scheduleReferenceMotion(; referencePath, referenceMotion) = scheduleReferenceMotion(referencePath, referenceMotion)

#-------------- Commands that can be used in referenceMotion programs -----
### ------------------------- ptpJointSpace -------------------------------

function ptpJointSpace(ref::ReferencePath, positions::Union{Matrix{Float64}, Vector{Float64}, Float64}, iargs...; kwargs...)
    checkErrorPtpJointSpace(iargs, kwargs)
    return push!(ref.refMotion, ArbitraryMotion(ptpJointSpace!, positions, nothing, nothing))
end
ptpJointSpace(;referencePath, positions) = ptpJointSpace(referencePath, positions)
ptpJointSpace(positions::Union{Matrix{Float64}, Vector{Float64}, Float64}, ref::ReferencePath, iargs...; kwargs...) =
ptpJointSpace(ref, positions, iargs...; kwargs...)
ptpJointSpace(ref::ReferencePath, iargs...; kwargs...) =
    checkErrorPtpJointSpace(iargs, kwargs)

#
# --------- check errors with might occur by the user -------------------
function checkErrorPtpJointSpace(iargs, kwargs)
    if length(iargs) > 0
        error("... from ptpJointSpace: first input argument must be a referencePath, second input argument must be position. No further or less input arguments are allowed.")
    end
    if length(kwargs) > 0
        error("... from ptpJointSpace: no keyword arguments are allowed. You used $kwargs.")
    end
    return nothing
end

### ------------------- calculateRobotMovement -----------------------
# is called from computeSignal or computeForceElement

calculateRobotMovement(; referencePath, instantiatedModel) = calculateRobotMovement(referencePath, instantiatedModel)

function calculateRobotMovement(referencePath, instantiatedModel)
    referencePath.instantiatedModel = instantiatedModel
    if ModiaLang.isEvent(instantiatedModel)
        if ModiaLang.getTime(instantiatedModel) < referencePath.nextEventTime
            # Reschedule next time event
            ModiaLang.setNextEvent!(instantiatedModel, referencePath.nextEventTime)
        else
            # At the current event instant a new command has to be processed
            # Execute next command in referenceMotion
            if !isempty(referencePath.refMotion)
                item = popfirst!(referencePath.refMotion)
                item.func(referencePath, item.inputArg1, item.inputArg2, item.inputArg3)
    end; end; end

    # only if a ptp path is defined getPosition! can be executed
    t = ModiaLang.getTime(instantiatedModel)
    #if t >= 0.1 && t <= 0.2
    #    println("... time = ", t, ", ptpPath = ", referencePath.ptpPath)
    #end
    if !isnothing(referencePath.ptpPath)
        getPosition!(referencePath.ptpPath, ModiaLang.getTime(instantiatedModel),
            referencePath.position,
            referencePath.velocity,
            referencePath.acceleration)
        #println("... time = ", t, ", position = ", referencePath.position)
    end
    return referencePath
end


text:
module VisualRegressionTests

using Requires
using FileIO
using ColorTypes
using ColorVectorSpace
using ImageFiltering

# ---------------------------------------------

struct VisualTest
    testFilename::AbstractString
    testFunction::Function # function taking (output_filename, args...; kw...) and writing to output_filename
    args
    kw

    refFilename::AbstractString
end

function VisualTest(testFunction::Function, refFilename::AbstractString, args...; kw...)
    VisualTest(tempname() * ".png", testFunction, args, kw, refFilename)
end

# ---------------------------------------------

@enum VisualTestStatus EXACT_MATCH CLOSE_MATCH DOES_NOT_MATCH PROCESSING_ERROR

mutable struct VisualTestResult
    testFilename::AbstractString
    testImage
    refFilename::AbstractString
    refImage
    status::VisualTestStatus
    diff::Float64
    err
end

Base.success(result::VisualTestResult) = result.status in (EXACT_MATCH, CLOSE_MATCH)

# ---------------------------------------------

export
    compare_images,
    test_images,
    VisualTest,
    VisualTestStatus,
    VisualTestResult,
    EXACT_MATCH,
    CLOSE_MATCH,
    DOES_NOT_MATCH,
    PROCESSING_ERROR,
    @visualtest,
    @plottest

include("utils.jl")
include("imgcomp.jl")
include("macros.jl")
function __init__()
    @require Gtk = "4c0ca9eb-093a-5379-98c5-f87ac0bbbf44" include("gui.jl")
end

end # module


text:
######################################################################
# Distributions Package UnivariateDistribution
######################################################################

#################### Categorical ####################

Categorical{T<:Real}(p::AbstractVector{T}) =
  Categorical(convert(Vector{Float64}, p))

######################################################################
# Distributions Package MultivariateDistribution
######################################################################

#################### Dirichlet ####################

Dirichlet{T<:Real}(alpha::AbstractVector{T}) =
  Dirichlet(convert(Vector{Float64}, alpha))

#################### Multinomial ####################

Multinomial{T<:Real}(n::Real, p::AbstractVector{T}) =
  Multinomial(convert(Int, n), convert(Vector{Float64}, p))

#################### MvNormal ####################

MvNormal{T<:Real, U<:Real}(μ::AbstractVector{T}, Σ::AbstractMatrix{U}) =
  MvNormal(convert(Vector{Float64}, μ), convert(Matrix{Float64}, Σ))

MvNormal{T<:Real, U<:Real}(μ::AbstractVector{T}, σ::AbstractVector{U}) =
  MvNormal(convert(Vector{Float64}, μ), convert(Vector{Float64}, σ))

MvNormal{T<:Real}(μ::AbstractVector{T}, σ::Real) =
  MvNormal(convert(Vector{Float64}, μ), convert(Float64, σ))

MvNormal{T<:Real}(Σ::AbstractMatrix{T}) =
  MvNormal(convert(Matrix{Float64}, Σ))

MvNormal{T<:Real}(σ::AbstractVector{T}) =
  MvNormal(convert(Vector{Float64}, σ))

#################### MvNormalCanon ####################

MvNormalCanon{T<:Real, U<:AbstractMatrix}(h::AbstractVector{T}, J::U) =
  MvNormalCanon(convert(Vector{Float64}, h), convert(Matrix{Float64}, J))

MvNormalCanon{T<:Real, U<:AbstractVector}(h::AbstractVector{T}, prec::U) =
  MvNormalCanon(convert(Vector{Float64}, h), convert(Vector{Float64}, prec))

MvNormalCanon{T<:Real, U<:Real}(h::AbstractVector{T}, prec::U) =
  MvNormalCanon(convert(Vector{Float64}, h), convert(Float64, prec))

MvNormalCanon{T<:Real}(J::AbstractMatrix{T}) =
  MvNormalCanon(convert(Matrix{Float64}, J))

MvNormalCanon{T<:Real}(prec::AbstractVector{T}) =
  MvNormalCanon(convert(Vector{Float64}, prec))

MvNormalCanon{T<:Real}(d::Real, prec::T) =
  MvNormalCanon(convert(Int, d), convert(Float64, prec))

#################### MvTDist ####################

MvTDist{T<:Real}(df::Real, μ::AbstractVector{T}, C::PDMat) =
  MvTDist(convert(Float64, df), convert(Vector{Float64}, μ), C)

MvTDist{T<:Real, U<:Real}(df::Real, μ::AbstractVector{T}, Σ::AbstractMatrix{U}) =
  MvTDist(convert(Float64, df), convert(Vector{Float64}, μ),
          convert(Matrix{Float64}, Σ))

MvTDist{T<:Real}(df::Real, Σ::AbstractMatrix{T}) =
  MvTDist(convert(Float64, df), convert(Matrix{Float64}, Σ))

#################### VonMisesFisher ####################

VonMisesFisher{T<:Real}(μ::AbstractVector{T}, κ::Real) =
  VonMisesFisher(convert(Vector{Float64}, μ), convert(Float64, κ))

######################################################################
# Distributions Package MatrixDistribution
######################################################################

#################### InverseWishart ####################

InverseWishart{T<:Real}(df::Real, Ψ::AbstractMatrix{T}) =
  InverseWishart(convert(Float64, df), convert(Matrix{Float64}, Ψ))

#################### Wishart ####################

Wishart{T<:Real}(df::Real, S::AbstractMatrix{T}) =
  Wishart(convert(Float64, df), convert(Matrix{Float64}, S))


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

using Markdown
using InteractiveUtils

# ╔═╡ 50ec2550-7256-11ec-1908-fbece8adaf13
begin
	import Pkg
	Pkg.activate(@__DIR__)
	Pkg.instantiate()
	using ThreeBodyDecay
	using Plots
	import Plots.PlotMeasures: mm
	#
	using PlutoUI
end

# ╔═╡ 602a1012-cf64-425c-a347-c582253255dc
md"""
# Comparison of the Dalitz-plot models
"""

# ╔═╡ c5c79442-02d4-4008-a76d-b7e255cce202
theme(:wong2, grid=false, frame=:box,
	xlims=(:auto,:auto), ylim=(:auto,:auto), lab="")

# ╔═╡ 72076d51-bb4e-4eb8-91fe-66f2515fa78c
begin
	const mπ = 0.14
	const m0 = 1.04 #0.782
	const mρ = 0.77
	const Γρ = 0.15
end;

# ╔═╡ 08b66d3d-7ce6-4bf9-9971-ec799d8e5a2e
md"""
## The case  $\omega \to 3\pi$: small dalitz, tails of $\rho$
"""

# ╔═╡ f6bc93cc-6660-4021-b6d6-d75e02ad439c
const ms = ThreeBodyMasses(mπ,mπ,mπ; m0)

# ╔═╡ d23809a4-1174-4150-adb6-875ca20c567a
"""
	ϕ(σs)

Kibble function of Mandelstam invatiant variables
"""
function ϕ(σs)
	ms²=ms^2
	return	λ(λ(ms²[4],σs[1],ms²[1]),
			  λ(ms²[4],σs[2],ms²[2]),
			  λ(ms²[4],σs[3],ms²[3]))
end

# ╔═╡ 37d00c66-5784-42fc-8d05-2fbb6ee00ceb
"""
	ρ(k,σs)

Three-body phase-space factor projected to one dimension
"""
function ρ(k,σs)
	ms²=ms^2
	i,j,_ = ijk(k)
	return	sqrt(λ(ms²[4],σs[k],ms²[k])*λ(σs[k],ms²[i],ms²[j]))
end

# ╔═╡ e4147784-1a90-40c2-868a-aa41412a760b
begin
	abstract type Model end
	struct BW <: Model
		m::Float64
		Γ::Float64
	end
	amplitude(ξ::BW, σ) = ξ.m*ξ.Γ/(ξ.m^2-σ-1im*ξ.m*ξ.Γ)
end

# ╔═╡ a1c461c6-108e-4ab3-8e2b-ab5286dc26f5
function A(ξ::Model, σs)
	ϕ0=ϕ(σs)
	return ϕ0>0 ? 0.0 : sqrt(-ϕ0) *
		  (amplitude(ξ,σs[1])+amplitude(ξ,σs[2])+amplitude(ξ,σs[3]))
end

# ╔═╡ 336e370e-66f8-4ed2-b3b1-e67c922f5c25
md"""
## Two model: regular $\rho$ and modified $\rho$
"""

# ╔═╡ d006f277-1673-4817-a934-a1aa9a0f7b91
begin
	ξ₁ = BW(mρ,Γρ) # nominal ρ
	ξ₂ = BW(mρ, Γρ+50e-3) # modified ρ
end;

# ╔═╡ 05477d32-20cb-450b-b891-85aefeb6ca42
begin
	plot(size=(700,300), layout=grid(1,2),
		leg=:left, ylab=["|A|²" "|A|²ρ"], xlab="m(ππ) [GeV]",
		bottom_margin=5mm, left_margin=5mm)
	#
	plot!(sp=1, m->abs2(amplitude(ξ₁,m^2)), sqrt.(lims1(ms))..., lab="model1")
	plot!(sp=1, m->abs2(amplitude(ξ₂,m^2)), sqrt.(lims1(ms))..., lab="model2")
	#
	plot!(sp=2, m->abs2(amplitude(ξ₁,m^2))*ρ(1,m^2), sqrt.(lims1(ms))...)
	plot!(sp=2, m->abs2(amplitude(ξ₂,m^2))*ρ(1,m^2), sqrt.(lims1(ms))...)
end

# ╔═╡ 002b72b4-8276-4054-83ac-a27cf3d66eba
I(ξ::Model, σs) = abs2(A(ξ,σs))

# ╔═╡ 8b0c8a14-1eb8-4843-b0ea-a5e7eb833b37
plot(size=(700,300), layout=grid(1,2),
	plot(Base.Fix1(I, ξ₁), ms),
	plot(Base.Fix1(I, ξ₂), ms))

# ╔═╡ 245b618d-0393-479b-abc6-070939ee6421
md"""
## MC toys: likelihood ratio on pseudodata
"""

# ╔═╡ 44e36584-dc7c-49ce-bd37-258e3e92dfa3
md"""
$\mathcal{L}_i(D) = \sum_\text{D}\log |A_i(s,t)|^2$

$\Delta\mathcal{L}(D) = \mathcal{L}_1(D) - \mathcal{L}_2(D)$
"""

# ╔═╡ 74f40b1c-b447-44cb-9d7b-753b1084a03f
binneddensity(ξ::Model) = getbinned2dDensity(
	(σ1,σ3)->I(ξ, Invariants(ms;σ1,σ3)), lims1(ms), lims3(ms),  50, 50);

# ╔═╡ 81f8ef1a-f4a4-4310-be39-b1fab4de2c55
function experiment(ξ₀,ξₜ,N)
	gds = binneddensity(ξ₀)
	pseudodata = [rand(gds) for _ in 1:N]
	σsv = [Invariants(ms; σ1, σ3) for (σ1,σ3) in pseudodata]
	Δllh = sum(σs->log(I(ξ₀, σs))-log(I(ξₜ, σs)), σsv)
	return Δllh/N
end

# ╔═╡ d0bcedb6-e1bc-4a2d-ada9-d0df93a09212
begin
	Nev = 500
	Nans = 200
end

# ╔═╡ 97af4995-8b3f-4c20-a358-d435b87fc3a0
begin
	ensample1 = [experiment(ξ₁,ξ₂,Nev) for _ in 1:Nans]
	ensample2 = -[experiment(ξ₂,ξ₁,Nev) for _ in 1:Nans]
end;

# ╔═╡ 3374ea60-45ea-4a02-aa3f-60c4f1487325
begin
	bins = range(extrema(vcat(ensample1,ensample2))..., length=30)
	plot(xlab="LLH₁-LLH₂")
	stephist!(ensample1; bins, lab="data(model1)")
	stephist!(ensample2; bins, lab="data(model2)")
end

# ╔═╡ 270bdb41-75f9-4e26-9008-5dc897136757
md"""
## Analytic expectation for toys
"""

# ╔═╡ ec2e5405-682f-44e5-924f-905e8477be06
function moment(ξ₀,ξₜ,n::Int)
	t(σs) = log(I(ξ₀,σs)) - log(I(ξₜ,σs))
	integrand(σs) = t(σs)^n * I(ξ₀,σs)
	cvalue = three_body_phase_space_integral(integrand, ms)
	return real(cvalue)
end

# ╔═╡ 308c6a64-fc93-4e69-ba25-9918c42046eb
function meansigma(ξ₀,ξₜ)
	N = moment(ξ₀,ξₜ,0)
	μN = moment(ξ₀,ξₜ,1)
	μ²N = moment(ξ₀,ξₜ,2)
	#
	μ = μN/N
	σ = sqrt(μ²N/N - (μN/N)^2)
	(; μ, σ)
end

# ╔═╡ 6620abf3-dd33-46ec-974d-ce2d608567a5
begin
	μ₁,σ₁ = meansigma(ξ₁,ξ₂)
	μ₂,σ₂ = meansigma(ξ₂,ξ₁)
end;

# ╔═╡ 45a8e3d5-20fe-48d1-b433-3d298a3a689e
let
	σₑ(N) = sqrt((σ₁^2+σ₂^2)/N)
	#
	separation(N) = (μ₁+μ₂)/σₑ(N) # = sqrt(N)*cosnt
	#
	plot(separation, 100:1_000_000, xscale=:log10,
		xlab="#events", ylab="separation in #σ")
end

# ╔═╡ Cell order:
# ╟─602a1012-cf64-425c-a347-c582253255dc
# ╠═50ec2550-7256-11ec-1908-fbece8adaf13
# ╠═c5c79442-02d4-4008-a76d-b7e255cce202
# ╠═72076d51-bb4e-4eb8-91fe-66f2515fa78c
# ╟─08b66d3d-7ce6-4bf9-9971-ec799d8e5a2e
# ╠═f6bc93cc-6660-4021-b6d6-d75e02ad439c
# ╠═d23809a4-1174-4150-adb6-875ca20c567a
# ╠═37d00c66-5784-42fc-8d05-2fbb6ee00ceb
# ╠═e4147784-1a90-40c2-868a-aa41412a760b
# ╠═a1c461c6-108e-4ab3-8e2b-ab5286dc26f5
# ╟─336e370e-66f8-4ed2-b3b1-e67c922f5c25
# ╠═d006f277-1673-4817-a934-a1aa9a0f7b91
# ╠═05477d32-20cb-450b-b891-85aefeb6ca42
# ╠═002b72b4-8276-4054-83ac-a27cf3d66eba
# ╠═8b0c8a14-1eb8-4843-b0ea-a5e7eb833b37
# ╟─245b618d-0393-479b-abc6-070939ee6421
# ╟─44e36584-dc7c-49ce-bd37-258e3e92dfa3
# ╠═74f40b1c-b447-44cb-9d7b-753b1084a03f
# ╠═81f8ef1a-f4a4-4310-be39-b1fab4de2c55
# ╠═d0bcedb6-e1bc-4a2d-ada9-d0df93a09212
# ╠═97af4995-8b3f-4c20-a358-d435b87fc3a0
# ╠═3374ea60-45ea-4a02-aa3f-60c4f1487325
# ╟─270bdb41-75f9-4e26-9008-5dc897136757
# ╠═ec2e5405-682f-44e5-924f-905e8477be06
# ╠═308c6a64-fc93-4e69-ba25-9918c42046eb
# ╠═6620abf3-dd33-46ec-974d-ce2d608567a5
# ╠═45a8e3d5-20fe-48d1-b433-3d298a3a689e

max_stars_repo_name: mmikhasenko/github-pages-test

max_issues_repo_name: mmikhasenko/github-pages-test

max_forks_repo_name: mmikhasenko/github-pages-test


text:
using Base.Test
using CGP
CGP.Config.init("cfg/test.yaml")

function test_mutate(c::Chromosome, child::Chromosome)
      @test child != c
      @test child.genes != c.genes
      @test child.nin == c.nin
      @test child.nout == c.nout
      @test distance(c, child) > 0
end

@testset "Mutation tests" begin
    for ct in CGP.chromosomes
        println(ct)
        nin = rand(1:100); nout = rand(1:100)
        c = ct(nin, nout)
        cgenes = deepcopy(c.genes)
        @testset "Clone $ct" begin
            copy = clone(c)
            @test copy.genes == c.genes
        end
        @testset "Constructor $ct" begin
            child = ct(c)
            test_mutate(c, child)
        end
        @testset "Mutate genes $ct" begin
            child = gene_mutate(c)
            test_mutate(c, child)
            @test length(child.genes) == length(c.genes)
        end
        @testset "Active genes $ct" begin
            child = active_gene_mutate(c)
            test_mutate(c, child)
            @test length(child.genes) == length(c.genes)
        end
        @testset "Add nodes $ct" begin
            child = add_nodes(c)
            test_mutate(c, child)
            @test length(child.genes) > length(c.genes)
            @test forward_connections(c) != forward_connections(child)
            @test issubset([n.p for n in c.nodes], [n.p for n in child.nodes])
            @test issubset([n.f for n in c.nodes], [n.f for n in child.nodes])
        end
        @testset "Delete nodes $ct" begin
            child = delete_nodes(c)
            test_mutate(c, child)
            @test length(child.genes) < length(c.genes)
            @test forward_connections(c) != forward_connections(child)
            @test issubset([n.p for n in child.nodes], [n.p for n in c.nodes])
            @test issubset([n.f for n in child.nodes], [n.f for n in c.nodes])
        end
        @testset "Mixed node mutate $ct" begin
            child = mixed_node_mutate(c)
            test_mutate(c, child)
        end
        @testset "Add subtree $ct" begin
            child = add_nodes(c)
            test_mutate(c, child)
            @test length(child.genes) > length(c.genes)
            @test forward_connections(c) != forward_connections(child)
            @test issubset([n.p for n in c.nodes], [n.p for n in child.nodes])
            @test issubset([n.f for n in c.nodes], [n.f for n in child.nodes])
        end
        @testset "Delete subtree $ct" begin
            child = delete_nodes(c)
            test_mutate(c, child)
            @test length(child.genes) < length(c.genes)
            @test forward_connections(c) != forward_connections(child)
            @test issubset([n.p for n in child.nodes], [n.p for n in c.nodes])
            @test issubset([n.f for n in child.nodes], [n.f for n in c.nodes])
        end
        @testset "Mixed subtree mutate $ct" begin
            child = mixed_subtree_mutate(c)
            test_mutate(c, child)
        end
    end
end

nothing


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule flux_core_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("flux_core")
JLLWrappers.@generate_main_file("flux_core", UUID("381ff186-7949-5def-80e3-c2e748080dc1"))
end  # module flux_core_jll


text:
@testset "Contrast Stretching" begin

    for T in (Gray{N0f8}, Gray{N0f16}, Gray{Float32}, Gray{Float64})
        #=
        Create an image that spans a narrow graylevel range.  Then
        quantize the 256 bins down to 16 and determine how many bins have
        non-zero counts.
        =#
        img = T.([i/255.0 for i = 64:128, j = 1:10])
        _, counts_before = build_histogram(img,16, minval = 0, maxval = 1)
        nonzero_before = sum(counts_before .!= 0)

        #=
        Stretch the histogram. Then quantize the 256 bins down to 16 and
        verify that all 16 bins have non-zero counts. This will confirm that
        the dynamic range of the original image has been increased.
        =#
        ret = adjust_histogram(img, ContrastStretching(t = 0.4, slope = 17))
        edges, counts_after = build_histogram(ret,16, minval = 0, maxval = 1)
        nonzero_after = sum(counts_after .!= 0)
        @test nonzero_before < nonzero_after
        @test nonzero_after == 16
        @test eltype(img) == eltype(ret)

        #=
        Verify that the function can cope with a NaN value.
        =#
        if T <: Gray{Float32} || T <: Gray{Float64}
            img[1] = NaN
            ret = adjust_histogram(img, ContrastStretching(t = 0.4, slope = 17))
            edges, counts_after = build_histogram(ret,16, minval = 0, maxval = 1)
            nonzero_after = sum(counts_after .!= 0)
            @test nonzero_before < nonzero_after
            @test nonzero_after == 16
            @test eltype(img) == eltype(ret)
        end

        #=
        Verify that when the slope is set to a very large value the contrast
        streching behaves like a thresholding function.
        =#
        ret = adjust_histogram(img, ContrastStretching(t = 0.37, slope = 1000))
        edges, counts_after = build_histogram(ret,16, minval = 0, maxval = 1)
        @test sum(counts_after .!= 0) == 2
    end

    for T in (RGB{N0f8}, RGB{N0f16}, RGB{Float32}, RGB{Float64})
        #=
        Create a color image that spans a narrow graylevel range.  Then
        quantize the 256 bins down to 16 and determine how many bins have
        non-zero counts.
        =#
        imgg = Gray{Float32}.([i/255.0 for i = 64:128, j = 1:10])
        img = colorview(RGB,imgg,imgg,imgg)
        img = T.(img)
        _, counts_before = build_histogram(img, 16, minval = 0, maxval = 1)
        nonzero_before = sum(counts_before .!= 0)

        #=
        Stretch the histogram. Then quantize the 256 bins down to 16 and
        verify that all 16 bins have non-zero counts. This will confirm that
        the dynamic range of the original image has been increased.
        =#
        ret = adjust_histogram(img, ContrastStretching(t = 0.4, slope = 17))
        edges, counts_after = build_histogram(ret, 16, minval = 0, maxval = 1)
        nonzero_after = sum(counts_after .!= 0)
        @test nonzero_before < nonzero_after
        @test nonzero_after == 16
        @test eltype(img) == eltype(ret)

        #=
        Verify that when the slope is set to a very large value the contrast
        streching behaves like a thresholding function.
        =#
        ret = adjust_histogram(img, ContrastStretching(t = 0.37, slope = 1000))
        edges, counts_after = build_histogram(ret, 16, minval = 0, maxval = 1)
        @test sum(counts_after .!= 0) == 2
    end
end


text:
using DimensionalData, Statistics, Test, BenchmarkTools, Unitful

using DimensionalData: val, basetype, slicedims, dims2indices, formatdims,
      @dim, reducedims, dimnum, basetype, X, Y, Z, Time, Forward

include("dimension.jl")
include("primitives.jl")
include("array.jl")
include("selector.jl")
include("methods.jl")
include("benchmarks.jl")


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

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

### data directories ###
data_dirs = ["ALOI", "Annthyroid", "Arrhythmia", "Cardiotocography", "Glass",
             "HeartDisease", "Hepatitis", "Ionosphere", "KDDCup99", "Lymphography",
             "PageBlocks", "Parkinson", "PenDigits", "Pima", "Shuttle", "SpamBase",
             "Stamps", "WBC", "WDBC", "WPBC", "Waveform"]

### learning scenario ###
initial_pool_strategy = [("Pu", Dict())]
split_strategy = [("Sf", Dict())]

init_strategies = [SimpleCombinedStrategy(RuleOfThumbScott(), BoundedTaxErrorEstimate(0.05, 0.02, 0.98)),
                   WangCombinedInitializationStrategy(solver, 2.0.^range(-4, stop=4, step=1.0), BoundedTaxErrorEstimate(0.05, 0.02, 0.98))]

#### models ####
models = [Dict(:type => :SVDDneg, :param => Dict{Symbol, Any}()),
          Dict(:type => :SSAD, :param => Dict{Symbol, Any}(:κ => 0.1))]
classify_precision = SVDD.OPT_PRECISION

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

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

num_al_iterations = 100
num_resamples_initial_pool = 5
exp_name = "evaluation_part2"

max_stars_repo_name: kit-dbis/ocal-evaluation

max_issues_repo_name: kit-dbis/ocal-evaluation

max_forks_repo_name: kit-dbis/ocal-evaluation


text:
#= Run as:
cd(".../Modia/test")
include("runsymbolictests.jl")
=#

using Base.Test

@testset "SymbolicTests" begin

include("symbolic/BLTandPantelides/setup.jl")
include("symbolic/BLTandPantelides/testBLTandPantelides.jl")

include("symbolic/DAEquations/setup.jl")
include("symbolic/DAEquations/testSymbolicTransform.jl")

include("symbolic/TestExactlyRemoveSingularities.jl")

include("symbolic/TestStateSelectionAlgorithm.jl")

include("symbolic/TestTearingAlgorithm.jl")

end
nothing


text:
# PREAMBLE

# PKG_SETUP

# ## Setup

using ContactImplicitMPC
using LinearAlgebra

# ## Simulation
s = get_simulation("hopper_2D", "flat_2D_lc", "flat");
model = s.model
env = s.env
nq = model.nq
nu = model.nu
nc = model.nc

# ## Run Policy
function run_policy(s::Simulation; H_sim::Int = 4000, verbose = verbose,
		offset::AbstractVector = zeros(s.model.nq))
	model = s.model
	env = s.env
	ref_traj = deepcopy(get_trajectory(model, env,
	    joinpath(module_dir(), "src/dynamics/hopper_2D/gaits/gait_in_place.jld2"),
	    load_type = :joint_traj))

	H = ref_traj.H
	h = ref_traj.h
	N_sample = 5
	H_mpc = 10
	h_sim = h / N_sample
	κ_mpc = 1.0e-4

	obj = TrackingVelocityObjective(model, env, H_mpc,
		q = [[Diagonal(1.0e-2 * [0.1,3,1,3])   for t = 1:H_mpc-2]; [Diagonal(1.0e-1 * [0.1,3,1,3])   for t = 1:2]],
		v = [[Diagonal(1.0e+1 * [0.1,3,1,3])   for t = 1:2]; [Diagonal(1.0e-3 * [0.1,3,1,3])   for t = 1:H_mpc-2]],
	    u = [Diagonal(1.0e-0 * [3e-3, 1e0]) for t = 1:H_mpc],
	    γ = [Diagonal(1.0e-100 * ones(model.nc)) for t = 1:H_mpc],
	    b = [Diagonal(1.0e-100 * ones(model.nc * friction_dim(env))) for t = 1:H_mpc])

	p = ci_mpc_policy(ref_traj, s, obj,
	    H_mpc = H_mpc,
	    N_sample = N_sample,
	    κ_mpc = κ_mpc,
		mode = :configuration,
	    n_opts = NewtonOptions(
			r_tol = 3e-4,
			max_iter = 5,
			max_time = ref_traj.h, # HARD REAL TIME
			),
	    mpc_opts = CIMPCOptions(
	        ),
		ip_opts = InteriorPointOptions(
			undercut = 5.0,
			κ_tol = κ_mpc,
			r_tol = 1.0e-8,
			diff_sol = true,
			solver = :empty_solver,
			max_time = 1e5),
	    )

	q1_ref = copy(ref_traj.q[2]) + offset
	q0_ref = copy(ref_traj.q[1]) + offset
	q1_sim = ContactImplicitMPC.SVector{model.nq}(q1_ref)
	q0_sim = ContactImplicitMPC.SVector{model.nq}(copy(q1_sim - (q1_ref - q0_ref) / N_sample))
	@assert norm((q1_sim - q0_sim) / h_sim - (q1_ref - q0_ref) / h) < 1.0e-8

	v1_sim = (ref_traj.q[2] - ref_traj.q[1]) ./ h

	sim = simulator(s, H_sim, h=h_sim, policy=p);

	status = simulate!(sim, q1_sim, v1_sim);

	return deepcopy(sim.traj)
end

# ## Collect runs
function collect_runs(s::Simulation; n::Int = 1, H_sim::Int = 4000, verbose::Bool = false)
	trajs = []
	ContactImplicitMPC.Random.seed!(100)
	offset_max = [+0.50, 0.30, +0.30, +0.00]
	offset_min = [-0.50, 0.00, -0.30, -0.30]
	offset_Δ = offset_max .- offset_min
	for i = 1:n
		verbose && println("sample = $i/$n")
		offset = offset_min .+ offset_Δ .* ContactImplicitMPC.rand(s.model.nq)
		traj = run_policy(s, H_sim = H_sim, offset = offset, verbose = verbose)
		push!(trajs, traj)
	end
	return trajs
end

# ## Visualize runs
function visualize_runs!(vis::ContactImplicitMPC.Visualizer, model::Model, trajs::AbstractVector;
		sample=max(1, Int(floor(trajs[1].H / 100))), h=trajs[1].h*sample,  α=1.0,
		anim::ContactImplicitMPC.MeshCat.Animation=ContactImplicitMPC.MeshCat.Animation(Int(floor(1/h))),
		name::Symbol=ContactImplicitMPC.model_name(model))

	for (i, traj) in enumerate(trajs)
		name_i = Symbol(string(name) * string(i))
		anim = visualize_robot!(vis, model, traj, sample = sample, h = h, α = α, anim = anim, name = name_i)
	end
	return anim
end

# ## Experiment
H_sim = 1000
trajs = collect_runs(s, n = 100, H_sim = H_sim, verbose = true);

# ## Visualizer
vis = ContactImplicitMPC.Visualizer()
ContactImplicitMPC.render(vis)

# ## Visualize
anim = visualize_runs!(vis, s.model, trajs, α=0.1, sample=5)


text:
using DrWatson
@quickactivate :FractalDimension # uses DynamicalSystems, PyPlot
include(srcdir("style.jl"))

ds = Systems.lorenz()
integ = integrator(ds)

function trajectory(ds::CDS{IIP, S, D}, integ, T, u, dt, Ttr) where {IIP, S, D}
    t0 = ds.t0
    tvec = (t0+Ttr):dt:(T+t0+Ttr)
    D = isnothing(a) ? D : length(a)
    sol = Vector{SVector{D, eltype(S)}}(undef, length(tvec))
    step!(integ, Ttr)
    for (i, t) in enumerate(tvec)
        while t > integ.t
            step!(integ)
        end
        sol[i] = SVector{D, eltype(S)}(obtain_access(integ(t), a))
    end
    return Dataset(sol)
end

max_stars_repo_name: Datseris/FractalDimension

max_issues_repo_name: Datseris/FractalDimension

max_forks_repo_name: Datseris/FractalDimension


text:
using PackageCompiler

import Pkg
Pkg.add("JLD")
Pkg.add("DataStructures")

create_sysimage([:JLD], sysimage_path="ext/EQ-JL/emews/worker.so")

max_stars_repo_name: emews/EQ-JL

max_issues_repo_name: emews/EQ-JL

max_forks_repo_name: emews/EQ-JL


text:
"""
```
get_scenario_input_file(m, scen::Scenario)
```

Get file name of raw scenario targets from `inpath(m, \"scenarios\")`.
"""
function get_scenario_input_file(m::AbstractDSGEModel, scen::Scenario)
    basename = string(scen.key) * "_" * scen.vintage * ".jld2"
    return inpath(m, "scenarios", basename)
end

"""
```
count_scenario_draws!(m, scen::Scenario)
```

Return the number of draws for `scen`, determined using
`get_scenario_input_file(m, scen)`, and update the
`n_draws` field of `scen` with this count.
"""
function count_scenario_draws!(m::AbstractDSGEModel, scen::Scenario)
    input_file = get_scenario_input_file(m, scen)
    draws = JLD2.jldopen(input_file, "r") do file
        dataset = read(file, "arr")
        size(dataset, 1)
    end
    scen.n_draws = draws
    return draws
end

"""
```
load_scenario_targets!(m, scen::Scenario, draw_index)
```

Add the targets from the `draw_index`th draw of the raw scenario targets to
`scen.targets`.
"""
function load_scenario_targets!(m::AbstractDSGEModel, scen::Scenario, draw_index::Int)
    if isempty(scen.target_names)
        n_periods = forecast_horizons(m)
        for target in keys(m.observables)
            scen.targets[target] = fill(NaN, n_periods)
        end
    else
        path = get_scenario_input_file(m, scen)
        raw_targets, target_inds = JLD2.jldopen(path, "r") do file
            arr = read(file, "arr")
            inds = read(file, "target_indices")
            arr[draw_index, :, :], inds
        end

        @assert collect(keys(target_inds)) == scen.target_names "Target indices in $path do not match target names in $(scen.key)"

        for (target_name, target_index) in target_inds
            scen.targets[!, target_name] = raw_targets[target_index, :]
        end
    end

    @assert collect(keys(target_inds)) == scen.target_names "Target indices in $path do not match target names in $(scen.key)"

    for (target_name, target_index) in target_inds
        scen.targets[!, target_name] = raw_targets[target_index, :]
    end

    return scen.targets
end

"""
```
get_scenario_filename(m, scen::AbstractScenario, output_var;
    pathfcn = rawpath, fileformat = :jld2, directory = "")
```

Get scenario file name of the form
`pathfcn(m, \"scenarios\", output_var * filestring * string(fileformat))`. If
`directory` is provided (nonempty), then the same file name in that directory
will be returned instead.
"""
function get_scenario_filename(m::AbstractDSGEModel, scen::AbstractScenario, output_var::Symbol;
                               pathfcn::Function = rawpath,
                               fileformat::Symbol = :jld2,
                               directory::String = "")
    filestring_addl = Vector{String}()
    if isa(scen, SingleScenario)
        push!(filestring_addl, "scen=" * string(scen.key))
    elseif isa(scen, ScenarioAggregate)
        push!(filestring_addl, "sagg=" * string(scen.key))
    end
    push!(filestring_addl, "svin=" * scen.vintage)

    base = string(output_var) * "." * string(fileformat)
    path = pathfcn(m, "scenarios", base, filestring_addl)
    if !isempty(directory)
        path = joinpath(directory, basename(path))
    end
    return path
end

"""
```
get_scenario_output_files(m, scen::SingleScenario, output_vars)
```

Return a `Dict{Symbol, String}` mapping `output_vars` to the raw simulated
scenario outputs for `scen`.
"""
function get_scenario_output_files(m::AbstractDSGEModel, scen::SingleScenario,
                                   output_vars::Vector{Symbol})
    output_files = Dict{Symbol, String}()
    for var in output_vars
        output_files[var] = get_scenario_filename(m, scen, var)
    end
    return output_files
end

"""
```
get_scenario_mb_input_file(m, scen::AbstractScenario, output_var)
```

Call `get_scenario_filename` while replacing `forecastut` and `forecast4q` in
`output_var` with `forecast`.
"""
function get_scenario_mb_input_file(m::AbstractDSGEModel, scen::AbstractScenario, output_var::Symbol)
    input_file = get_scenario_filename(m, scen, output_var)
    input_file = replace(input_file, "forecastut" => "forecast")
    input_file = replace(input_file, "forecast4q" => "forecast")
    return input_file
end

"""
```
get_scenario_mb_output_file(m, scen::AbstractScenario, output_var;
    directory = "")
```

Call `get_scenario_filename` while tacking on `\"mb\"` to the front of the base
file name.
"""
function get_scenario_mb_output_file(m::AbstractDSGEModel, scen::AbstractScenario, output_var::Symbol;
                                     directory::String = "")
    fullfile = get_scenario_filename(m, scen, output_var, pathfcn = workpath, directory = directory)
    joinpath(dirname(fullfile), "mb" * basename(fullfile))
end

"""
```
get_scenario_mb_metadata(m, scen::SingleScenario, output_var)

get_scenario_mb_metadata(m, agg::ScenarioAggregate, output_var)
```

Return the `MeansBands` metadata dictionary for `scen`.
"""
function get_scenario_mb_metadata(m::AbstractDSGEModel, scen::SingleScenario, output_var::Symbol)
    forecast_output_file = get_scenario_mb_input_file(m, scen, output_var)
    metadata = get_mb_metadata(m, :mode, :none, output_var, forecast_output_file)
    metadata[:scenario_key] = scen.key
    metadata[:scenario_vint] = scen.vintage

    return metadata
end

function get_scenario_mb_metadata(m::AbstractDSGEModel, agg::ScenarioAggregate, output_var::Symbol)
    forecast_output_file = get_scenario_mb_input_file(m, agg.scenarios[1], output_var)
    metadata = get_mb_metadata(m, :mode, :none, output_var, forecast_output_file)

    # Initialize start and end date
    start_date = date_forecast_start(m)
    end_date   = maximum(keys(metadata[:date_inds]))

    for scen in agg.scenarios
        forecast_output_file = get_scenario_mb_input_file(m, scen, output_var)
        scen_dates = JLD2.jldopen(forecast_output_file, "r") do file
            read(file, "date_indices")
        end

        # Throw error if start date for this scenario doesn't match
        if Dict([reverse(scen_date) for scen_date = pairs(scen_dates)])[1] != start_date
            error("All scenarios in agg must start from the same date")
        end

        # Update end date if necessary
        end_date = max(end_date, maximum(keys(metadata[:date_inds])))
    end

    dates = quarter_range(start_date, end_date)
    metadata[:date_inds] = OrderedDict{Date, Int}(d => i for (i, d) in enumerate(dates))
    metadata[:scenario_key] = agg.key
    metadata[:scenario_vint] = agg.vintage

    return metadata
end

"""
```
read_scenario_output(m, scen::SingleScenario, class, product, var_name)

read_scenario_output(m, agg::ScenarioAggregate, class, product, var_name)
```

Given either `scen` or `agg`, read in and return all draws of and the
appropriate reverse transform for `var_name`.

The third function that takes in two models is used for when we have scenarios from two different models.
"""
function read_scenario_output(m::AbstractDSGEModel, scen::SingleScenario, class::Symbol, product::Symbol,
                              var_name::Symbol)
    # Get filename
    filename = get_scenario_mb_input_file(m, scen, Symbol(product, class))

    fcast_dates_dict = load(filename, "date_indices")
    fcast_dates = map(x -> x[1], sort(collect(fcast_dates_dict), by = x->x[2]))

    jldopen(filename, "r") do file
        # Get index corresponding to var_name
        class_long = get_class_longname(class)
        indices = FileIO.load(filename, "$(class_long)_indices")
        var_ind = indices[var_name]

        # Read forecast outputs
        fcast_series = read_forecast_series(filename, product, var_ind)

        # Parse transform
        class_long = get_class_longname(class)
        transforms = load(filename, string(class_long) * "_revtransforms")
        transform = parse_transform(transforms[var_name])

        fcast_series, transform, fcast_dates
    end
end

function read_scenario_output(m::AbstractDSGEModel, m904::AbstractDSGEModel, agg::ScenarioAggregate, class::Symbol,
                              product::Symbol, var_name::Symbol)
    # Aggregate scenarios
    nscens = length(agg.scenarios)
    agg_draws = Vector{Matrix{Float64}}(undef, nscens)
    scen_dates = Vector{Date}(undef, nscens)

    # If not sampling, initialize vector to record number of draws in each
    # scenario in order to update `agg.proportions` and `agg.total_draws` at the
    # end
    if !agg.sample
        n_scen_draws = zeros(Int, nscens)
    end

    # Initialize transform so it can be assigned from within the following for
    # loop. Each transform read in from read_scenario_output will be the
    # same. We just want to delegate the transform parsing to the recursive
    # read_scenario_output call.
    transform = identity

    for (i, scen) in enumerate(agg.scenarios)
        if in(:scenarios, fieldnames(typeof(scen))) #length(scen.scenarios)>1
            scen_draws, transform, scen_dates = read_scenario_output(m, m904, scen, class, product, var_name)
        else
            if scen.key==:bor8 || scen.key==:bor9 || scen.key==:bor8_02 || scen.key==:bor9_02
                if var_name==:obs_corepce
                    var_name = :obs_gdpdeflator
                end
                # Recursively read in scenario draws
                scen_draws, transform, scen_dates = read_scenario_output(m904, scen, class, product, var_name)
            else
                # Recursively read in scenario draws
                scen_draws, transform, scen_dates = read_scenario_output(m, scen, class, product, var_name)
            end
        end
        # Sample if desired
        agg_draws[i]  = if agg.sample
            pct = agg.proportions[i]
            actual_ndraws = size(scen_draws, 1)
            desired_ndraws = convert(Int, round(pct * agg.total_draws))

            sampled_inds = if agg.replace
                sample(1:actual_ndraws, desired_ndraws, replace = true)
            else
                if desired_ndraws == 0
                    Int[]
                else
                    quotient  = convert(Int, floor(actual_ndraws / desired_ndraws))
                    remainder = actual_ndraws % desired_ndraws
                    vcat(repeat(1:actual_ndraws, quotient),
                         sample(1:actual_ndraws, remainder, replace = false))
                end
            end
            sort!(sampled_inds)
            scen_draws[sampled_inds, :]
        else
            # Record number of draws in this scenario
            n_scen_draws[i] = size(scen_draws, 1)
            scen_draws
        end
    end

    # Stack draws from all component scenarios
    fcast_series = cat(agg_draws..., dims = 1)

    # If not sampling, update `agg.proportions` and `agg.total_draws`
    if !agg.sample
        agg.total_draws = sum(n_scen_draws)
        agg.proportions = n_scen_draws ./ agg.total_draws
    end

    return fcast_series, transform, scen_dates
end

#=function read_scenario_output(m::AbstractDSGEModel, m904::AbstractDSGEModel, agg::ScenarioAggregate, class::Symbol,
                              product::Symbol, var_name::Symbol)
    # Aggregate scenarios
    nscens = length(agg.scenarios)
    agg_draws = Vector{Matrix{Float64}}(nscens)

    # If not sampling, initialize vector to record number of draws in each
    # scenario in order to update `agg.proportions` and `agg.total_draws` at the
    # end
    if !agg.sample
        n_scen_draws = zeros(Int, nscens)
    end

    # Initialize transform so it can be assigned from within the following for
    # loop. Each transform read in from read_scenario_output will be the
    # same. We just want to delegate the transform parsing to the recursive
    # read_scenario_output call.
    transform = identity

    for (i, scen) in enumerate(agg.scenarios)
        @show scen
        @show scen.key
        if in(:scenarios, fieldnames(scen)) #length(scen.scenarios)>1
            scen_draws, transform = read_scenario_output(m, m904, scen, class, product, var_name)
        else
            # If BOR8 or BOR9 (switching or non-switching), read GDP Deflator instead of core PCE
            if scen.key==:bor8 || scen.key==:bor9 || scen.key==:bor8_02 || scen.key==:bor9_02
                if var_name==:obs_corepce
                    var_name = :obs_gdpdeflator
                end
                # Recursively read in scenario draws
                scen_draws, transform = read_scenario_output(m904, scen, class, product, var_name)
            else
                # Recursively read in scenario draws
                scen_draws, transform = read_scenario_output(m, scen, class, product, var_name)
            end
        end
        # Sample if desired
        agg_draws[i] = if agg.sample
            pct = agg.proportions[i]
            actual_ndraws = size(scen_draws, 1)
            desired_ndraws = convert(Int, round(pct * agg.total_draws))

            sampled_inds = if agg.replace
                sample(1:actual_ndraws, desired_ndraws, replace = true)
            else
                if desired_ndraws == 0
                    Int[]
                else
                    quotient  = convert(Int, floor(actual_ndraws / desired_ndraws))
                    remainder = actual_ndraws % desired_ndraws
                    vcat(repmat(1:actual_ndraws, quotient),
                         sample(1:actual_ndraws, remainder, replace = false))
                end
            end
            sort!(sampled_inds)
            scen_draws[sampled_inds, :]
        else
            # Record number of draws in this scenario
            n_scen_draws[i] = size(scen_draws, 1)
            scen_draws
        end
    end

    # Stack draws from all component scenarios
    fcast_series = cat(1, agg_draws...)

    # If not sampling, update `agg.proportions` and `agg.total_draws`
    if !agg.sample
        agg.total_draws = sum(n_scen_draws)
        agg.proportions = n_scen_draws ./ agg.total_draws
    end

    return fcast_series, transform
end=#

function read_scenario_output(m::AbstractDSGEModel, agg::ScenarioAggregate, class::Symbol,
                              product::Symbol, var_name::Symbol)
    # Aggregate scenarios
    nscens = length(agg.scenarios)
    agg_draws = Vector{Matrix{Float64}}(undef, nscens)
    agg_dates = Vector{Vector{Date}}(undef, nscens)

    # If not sampling, initialize vector to record number of draws in each
    # scenario in order to update `agg.proportions` and `agg.total_draws` at the
    # end
    if !agg.sample
        n_scen_draws = zeros(Int, nscens)
    end

    # Initialize transform so it can be assigned from within the following for
    # loop. Each transform read in from read_scenario_output will be the
    # same. We just want to delegate the transform parsing to the recursive
    # read_scenario_output call.
    transform = identity

    for (i, scen) in enumerate(agg.scenarios)
        # Recursively read in scenario draws
        scen_draws, transform, scen_dates = read_scenario_output(m, scen, class, product, var_name)

        # Sample if desired
        agg_draws[i], agg_dates[i] = if agg.sample
            pct = agg.proportions[i]
            actual_ndraws = size(scen_draws, 1)
            desired_ndraws = convert(Int, round(pct * agg.total_draws))

            sampled_inds = if agg.replace
                sample(1:actual_ndraws, desired_ndraws, replace = true)
            else
                if desired_ndraws == 0
                    Int[]
                else
                    quotient  = convert(Int, floor(actual_ndraws / desired_ndraws))
                    remainder = actual_ndraws % desired_ndraws
                    vcat(repeat(1:actual_ndraws, outer = quotient),
                         sample(1:actual_ndraws, remainder, replace = false))
                end
            end
            sort!(sampled_inds)
            scen_draws[sampled_inds, :], scen_dates
        else
            # Record number of draws in this scenario
            n_scen_draws[i] = size(scen_draws, 1)
            scen_draws, scen_dates
        end
    end

    # Stack draws from all component scenarios
    fcast_series = cat(agg_draws..., dims = 1)

   #= fcast_dates_dict = load(get_scenario_mb_input_file(m, scen, Symbol(:forecast, class)), "date_indices")
    fcast_dates = map(x -> x[1], sort(collect(fcast_draws_dates_dict), by = x->x[2]))=#

    # If not sampling, update `agg.proportions` and `agg.total_draws`
    if !agg.sample
        agg.total_draws = sum(n_scen_draws)
        agg.proportions = n_scen_draws ./ agg.total_draws
    end

    return fcast_series, transform, agg_dates
end

"""
```
read_scenario_mb(m, scen::AbstractScenario, output_var; directory = "")
```

Read in an alternative scenario `MeansBands` object.
"""
function read_scenario_mb(m::AbstractDSGEModel, scen::AbstractScenario, output_var::Symbol;
                          directory::String = "")
    filepath = get_scenario_mb_output_file(m, scen, output_var, directory = directory)
    read_mb(filepath)
end


text:
using Meshes
# jit compile
mesh("cube.stl")
Profile.clear() # in case we have any previous profiling data
@profile mesh("cube.stl")
using ProfileView
ProfileView.view()


text:
export HamiltonianSplittingSpin

"""
    HamiltonianSplittingSpin( maxwell_solver,
                              kernel_smoother_0, kernel_smoother_1,
                              particle_group, e_dofs, a_dofs )

Hamiltonian splitting type for Vlasov-Maxwell

- Integral over the spline function on each interval (order p+1)
- Integral over the spline function on each interval (order p)
- `e_dofs` describing the two components of the electric field
- `a_dofs` describing the potential vector
"""
struct HamiltonianSplittingSpin
    maxwell_solver::AbstractMaxwellSolver
    kernel_smoother_0::ParticleMeshCoupling1D
    kernel_smoother_1::ParticleMeshCoupling1D
    kernel_smoother_2::ParticleMeshCoupling1D
    particle_group::ParticleGroup

    spline_degree::Int64
    Lx::Float64
    x_min::Float64
    delta_x::Float64

    cell_integrals_0::SVector
    cell_integrals_1::SVector

    e_dofs::Array{Array{Float64,1}}
    a_dofs::Array{Array{Float64,1}}
    j_dofs::Array{Array{Float64,1}}
    part::Array{Array{Float64,1}}

    HH::Float64

    function HamiltonianSplittingSpin(
        maxwell_solver,
        kernel_smoother_0,
        kernel_smoother_1,
        kernel_smoother_2,
        particle_group,
        e_dofs,
        a_dofs,
        HH=0.0,
    ) #0.00022980575)

        # Check that n_dofs is the same for both kernel smoothers.
        @assert kernel_smoother_0.n_dofs == kernel_smoother_1.n_dofs

        j_dofs = [zeros(Float64, kernel_smoother_0.n_dofs) for i in 1:2]

        nx = maxwell_solver.n_dofs
        part = [zeros(Float64, nx) for i in 1:4]
        x_min = maxwell_solver.xmin
        Lx = maxwell_solver.Lx
        spline_degree = 3
        delta_x = Lx / kernel_smoother_1.n_dofs

        cell_integrals_1 = SVector{3}([0.5, 2.0, 0.5] ./ 3.0)
        cell_integrals_0 = SVector{4}([1.0, 11.0, 11.0, 1.0] ./ 24.0)

        return new(
            maxwell_solver,
            kernel_smoother_0,
            kernel_smoother_1,
            kernel_smoother_2,
            particle_group,
            spline_degree,
            Lx,
            x_min,
            delta_x,
            cell_integrals_0,
            cell_integrals_1,
            e_dofs,
            a_dofs,
            j_dofs,
            part,
            HH,
        )
    end
end

export strang_splitting!

"""
    strang_splitting( h, dt, number_steps)

Strang splitting
- time splitting object
- time step
- number of time steps
"""
function strang_splitting!(h::HamiltonianSplittingSpin, dt::Float64, number_steps::Int64)
    for i_step in 1:number_steps
        operatorHE(h, 0.5dt)
        operatorHp(h, 0.5dt)
        operatorHA(h, 0.5dt)
        operatorHs(h, dt)
        operatorHA(h, 0.5dt)
        operatorHp(h, 0.5dt)
        operatorHE(h, 0.5dt)
    end
end

"""

    operatorHp(h, dt)

```math
\\begin{aligned}
\\dot{x} & =p \\\\
\\dot{E}_x & = - \\int (p f ) dp ds
\\end{aligned}
```
"""
function operatorHp(h::HamiltonianSplittingSpin, dt::Float64)
    n_cells = h.kernel_smoother_0.n_dofs

    fill!(h.j_dofs[1], 0.0)
    fill!(h.j_dofs[2], 0.0)

    for i_part in 1:(h.particle_group.n_particles)

        # Read out particle position and velocity
        x_old = h.particle_group.array[1, i_part]
        v_old = h.particle_group.array[2, i_part]

        # Then update particle position:  X_new = X_old + dt * V
        x_new = x_old + dt * v_old

        # Get charge for accumulation of j
        w = get_charge(h.particle_group, i_part)

        add_current_update_v!(h.j_dofs[1], h.kernel_smoother_1, x_old, x_new, w)

        x_new = mod(x_new, h.Lx)
        h.particle_group.array[1, i_part] = x_new
    end

    # Update the electric field.
    return compute_e_from_j!(h.e_dofs[1], h.maxwell_solver, h.j_dofs[1], 1)
end

"""
    operatorHA(h, dt)

```math
\\begin{aligned}
\\dot{p} = (A_y, A_z) \\cdot \\partial_x (A_y, A_z)   \\\\
\\dot{Ey} = -\\partial_x^2 A_y + A_y \\rho \\\\
\\dot{Ez} = -\\partial_x^2 A_z + A_z \\rho \\\\
\\end{aligned}
```
"""
function operatorHA(h::HamiltonianSplittingSpin, dt::Float64)
    n_cells = h.kernel_smoother_0.n_dofs

    fill!.(h.part, 0.0)

    qm = h.particle_group.q_over_m
    # Update velocities
    aa = zeros(Float64, n_cells)

    tmp_jdofs1 = zeros(h.kernel_smoother_1.n_span + 1)
    tmp_jdofs2 = zeros(h.kernel_smoother_0.n_span)

    for i_part in 1:(h.particle_group.n_particles)

        # fill!(h.j_dofs[1], 0.0)
        # fill!(h.j_dofs[2], 0.0)

        # Evaluate b at particle position (splines of order p)
        xi = h.particle_group.array[1, i_part]
        vi = h.particle_group.array[2, i_part]
        wi = get_charge(h.particle_group, i_part)

        # add_charge!( h.j_dofs[2], h.kernel_smoother_0, xi, 1.0)

        xn = (xi - h.kernel_smoother_0.xmin) / h.kernel_smoother_0.delta_x
        index = trunc(Int, xn)
        xn = xn - index
        index = index - h.kernel_smoother_0.spline_degree

        uniform_bsplines_eval_basis!(
            h.kernel_smoother_0.spline_val, h.kernel_smoother_0.spline_degree, xn
        )

        for i in 1:(h.kernel_smoother_0.n_span)
            tmp_jdofs2[i] = h.kernel_smoother_0.spline_val[i] * h.kernel_smoother_0.scaling
        end

        # add_charge!( h.j_dofs[1], h.kernel_smoother_1, xi, 1.0)

        xn = (xi - h.kernel_smoother_1.xmin) / h.kernel_smoother_1.delta_x
        index = trunc(Int, xn)
        xn = xn - index
        index = index - h.kernel_smoother_1.spline_degree

        uniform_bsplines_eval_basis!(
            h.kernel_smoother_1.spline_val, h.kernel_smoother_1.spline_degree, xn
        )

        for i in 1:(h.kernel_smoother_1.n_span)
            tmp_jdofs1[i] = h.kernel_smoother_1.spline_val[i] * h.kernel_smoother_1.scaling
        end

        # values of the derivatives of basis function
        # compute_rderivatives_from_basis!(aa, h.maxwell_solver, h.j_dofs[1])

        coef = 1 / h.maxwell_solver.delta_x
        # relation betwen spline coefficients for strong Ampere
        for i in 1:(h.kernel_smoother_1.n_span)
            index1d = mod1(index + i, n_cells)
            aa[index1d] = coef * (tmp_jdofs1[i] - tmp_jdofs1[i + 1])
        end
        ## treat Periodic point
        #aa[end] =  coef * ( field_in[end] - field_in[1] )
        # h.j_dofs[1] .= aa

        for i in 1:(h.kernel_smoother_1.n_span)
            index1d = mod1(index + i, n_cells)
            tmp_jdofs1[i] = aa[index1d]
        end

        #p11 = h.a_dofs[1]'h.j_dofs[1]
        #p21 = h.a_dofs[2]'h.j_dofs[1]

        #p12 = h.a_dofs[1]'h.j_dofs[2]
        #p22 = h.a_dofs[2]'h.j_dofs[2]

        p11 = 0.0
        p21 = 0.0
        p12 = 0.0
        p22 = 0.0
        for i in 1:(h.kernel_smoother_0.n_span)
            index1d = mod1(index + i, n_cells)
            p11 += h.a_dofs[1][index1d] * tmp_jdofs1[i]
            p21 += h.a_dofs[2][index1d] * tmp_jdofs1[i]
            p12 += h.a_dofs[1][index1d] * tmp_jdofs2[i]
            p22 += h.a_dofs[2][index1d] * tmp_jdofs2[i]
        end

        vi = vi - dt * (p11 * p12 + p22 * p21)

        h.particle_group.array[2, i_part] = vi

        # below we solve electric field
        # first define part1 and part2 to be 0 vector
        for i in 1:(h.kernel_smoother_0.n_span)
            index1d = mod1(index + i, n_cells)
            h.part[1][index1d] = h.part[1][index1d] - tmp_jdofs2[i] * dt * wi * p12
            h.part[2][index1d] = h.part[2][index1d] - tmp_jdofs2[i] * dt * wi * p22
        end
    end

    # Update the electric field. Also, we still need to scale with 1/Lx

    compute_e_from_j!(h.e_dofs[2], h.maxwell_solver, h.part[1], 2)
    compute_e_from_j!(h.e_dofs[3], h.maxwell_solver, h.part[2], 2)

    # define part3 and part4
    compute_lderivatives_from_basis!(h.part[3], h.maxwell_solver, h.a_dofs[1])
    compute_lderivatives_from_basis!(h.part[4], h.maxwell_solver, h.a_dofs[2])

    compute_e_from_b!(h.e_dofs[2], h.maxwell_solver, dt, h.part[3])
    return compute_e_from_b!(h.e_dofs[3], h.maxwell_solver, dt, h.part[4])
end

"""
    operatorHE(h, dt)

```math
\\begin{aligned}
\\dot{v}   & =  E_x \\\\
\\dot{A}_y & = -E_y \\\\
\\dot{A}_z & = -E_z
\\end{aligned}
```
"""
function operatorHE(h::HamiltonianSplittingSpin, dt::Float64)
    np = h.particle_group.n_particles

    if np > 2
        chunks = Iterators.partition(1:np, np ÷ nthreads())
    else
        chunks = Iterators.partition(1:np, nthreads())
    end

    @sync for i_chunk in chunks
        @spawn begin
            for i_part in i_chunk
                xi = h.particle_group.array[1, i_part]
                vi = h.particle_group.array[2, i_part]

                e1 = evaluate(h.kernel_smoother_1, xi, h.e_dofs[1])

                h.particle_group.array[2, i_part] = vi + dt * e1
            end
        end
    end

    h.a_dofs[1] .-= dt * h.e_dofs[2]
    return h.a_dofs[2] .-= dt * h.e_dofs[3]
end

"""
    operatorHs(h, dt)

Push H_s: Equations to be solved
```math
\\begin{aligned}
\\dot{s} &= s x B = (s_y \\partial_x A_y +s_z \\partial_x Az, -s_x \\partial_x A_y, -s_x \\partial_x A_z)  \\\\
\\dot{p} &= s \\cdot \\partial_x B = -s_y \\partial^2_{x} A_z + s_z \\partial^2_{x} A_y \\\\
\\dot{E}_y &=   \\int (s_z \\partial_x f) dp ds \\\\
\\dot{E}_z &= - \\int (s_y \\partial_x f) dp ds
\\end{aligned}
```
"""
function operatorHs(h::HamiltonianSplittingSpin, dt::Float64)
    n_cells = h.kernel_smoother_0.n_dofs

    fill!.(h.part, 0.0)

    hat_v = zeros(Float64, 3, 3)

    S = zeros(Float64, 3)
    St = zeros(Float64, 3)
    V = zeros(Float64, 3)
    aa = zeros(Float64, n_cells)

    @inbounds for i_part in 1:(h.particle_group.n_particles)
        xi = h.particle_group.array[1, i_part]
        vi = h.particle_group.array[2, i_part]

        # Evaluate efields at particle position
        fill!(h.j_dofs[1], 0.0)
        fill!(h.j_dofs[2], 0.0)

        add_charge!(h.j_dofs[2], h.kernel_smoother_1, xi, 1.0)

        compute_rderivatives_from_basis!(h.j_dofs[1], h.maxwell_solver, h.j_dofs[2])

        Y = h.a_dofs[1]'h.j_dofs[1]
        Z = h.a_dofs[2]'h.j_dofs[1]
        V .= [0, Z, -Y]

        hat_v[1, 2] = Y
        hat_v[1, 3] = Z
        hat_v[2, 1] = -Y
        hat_v[3, 1] = -Z

        s1 = h.particle_group.array[4, i_part]
        s2 = h.particle_group.array[5, i_part]
        s3 = h.particle_group.array[6, i_part]

        vnorm = norm(V)

        if vnorm > 1e-14
            S .= (
                [s1, s2, s3] .+
                (
                    sin(dt * vnorm) / vnorm * hat_v +
                    0.5 * (sin(dt / 2 * vnorm) / (vnorm / 2))^2 * hat_v^2
                ) * [s1, s2, s3]
            )
        else
            S .= [s1, s2, s3]
        end

        h.particle_group.array[4, i_part] = S[1]
        h.particle_group.array[5, i_part] = S[2]
        h.particle_group.array[6, i_part] = S[3]

        if vnorm > 1e-14
            St .= (
                dt .* [s1, s2, s3] +
                (
                    2 * (sin(dt * vnorm / 2) / vnorm)^2 * hat_v +
                    2.0 / (vnorm^2) * (dt / 2 - sin(dt * vnorm) / 2 / vnorm) * hat_v^2
                ) * [s1, s2, s3]
            )
        else
            St .= dt .* [s1, s2, s3]
        end

        wi = get_charge(h.particle_group, i_part)

        # define part1 and part2
        h.part[1] .+= wi[1] * St[3] * h.j_dofs[2]
        h.part[2] .+= wi[1] * St[2] * h.j_dofs[2]

        # update velocity
        fill!(h.j_dofs[2], 0.0)
        add_charge!(h.j_dofs[2], h.kernel_smoother_2, xi, 1.0)
        compute_rderivatives_from_basis!(h.j_dofs[1], h.maxwell_solver, h.j_dofs[2])
        compute_rderivatives_from_basis!(aa, h.maxwell_solver, h.j_dofs[1])
        h.j_dofs[1] .= aa
        vi =
            vi -
            h.HH * (h.a_dofs[2]' * h.j_dofs[1] * St[2] + h.a_dofs[1]' * h.j_dofs[1] * St[3])

        set_v!(h.particle_group, i_part, vi)
    end

    # Update bfield
    compute_rderivatives_from_basis!(h.j_dofs[1], h.maxwell_solver, h.part[1])
    compute_rderivatives_from_basis!(h.j_dofs[2], h.maxwell_solver, -h.part[2])

    h.j_dofs[1] .*= h.HH
    h.j_dofs[2] .*= h.HH

    compute_e_from_j!(h.e_dofs[2], h.maxwell_solver, h.j_dofs[1], 2)
    return compute_e_from_j!(h.e_dofs[3], h.maxwell_solver, h.j_dofs[2], 2)
end


text:
"""
    precalculateSOPRRB(model::SequencingModel)

A method that preloads and returns the calculations required to find a rough relaxed bound for SOPs

# Arguments
- `model::SequencingModel`: The SOP
"""
function precalculateSOPRRB(model::SequencingModel)
    paramType = valtype(model.objective.f.components)
    values = Dict{Int, Vector{Tuple{Int, paramType}}}()
    lastArc = Vector{Tuple{Int,paramType}}()
    for (states,objValue) in model.objective.f.components
        startState = states[1]
        endState = states[2]
        if haskey(values,startState)
            insert!(
                values[startState],
                searchsorted(values[startState], (endState, objValue), by = x -> x[2]).start,
                (endState, objValue)
            )
        else
            values[startState] = [(endState, objValue)]
        end
        if endState == model.seqEnd
            insert!(
                lastArc,
                searchsorted(lastArc, (startState, objValue), by = x -> x[2]).start,
                (startState, objValue)
            )
        end
    end
    return values, lastArc
end

"""
    sopRRB(visited::BitVector, state::U, startingValue::T, values::Dict{Int, Vector{Tuple{U, T}}},lastArc::Vector{Tuple{U, T}})where{T<:Real,U<:Int}

A method that takes in the current state of the sop and the precalculated values and returns a rough relaxe bound

# Arguments
- `visited::BitVector`: The visited nodes
- `state::U`: The current node
- `values::Dict{U, Vector{Tuple{U, T}}}`: The presorted out arc values
- `lastArc::Vector{Tuple{Int, T}}`: The presorted in arc values for the last edge
"""
function sopRRB(visited::BitVector, state::U, startingValue::T, values::Dict{Int, Vector{Tuple{U, T}}},lastArc::Vector{Tuple{U, T}})where{T<:Real,U<:Int}
    rrb = startingValue
    #add the first arc
    for pair in values[state]
        if !visited[pair[1]] && pair[1]!=length(visited)
            rrb += pair[2]
            break
        end
    end
    #add the last arc
    for pair in lastArc
        if !visited[pair[1]]
            rrb += pair[2]
            break
        end
    end

    worstValue = 0
    for i in 1:length(visited)-1
        if !visited[i]
            for pair in values[i]
                if !visited[pair[1]] && pair[1]!=length(visited)
                    rrb += pair[2]
                    if pair[2] > worstValue
                        worstValue = pair[2]
                    end
                    break
                end
            end
        end
    end
    rrb -= worstValue
    return rrb
end

"""
    sopRRAB(arc::Union{ExactArc,SimpleArc},values::Dict{Int, Vector{Tuple{U, T}}},lastArc::Vector{Tuple{U, T}})where{T<:Real,U<:Int}

A method that takes in an arc from a relaxed DD of the sop, and the precalculated values, then returns a rough relaxed bound

# Arguments
- `arc::Union{ExactArc,SimpleArc}`: The arc to value
- `values::Dict{U, Vector{Tuple{U, T}}}`: The presorted out arc values
- `lastArc::Vector{Tuple{Int, T}}`: The presorted in arc values for the last edge
"""
function sopRRAB(arc::Union{ExactArc,SimpleArc},values::Dict{Int, Vector{Tuple{U, T}}},lastArc::Vector{Tuple{U, T}})where{T<:Real,U<:Int}
    rrb = arc.value
    visited = copy(arc.origin.allDown)
    visited[arc.assignment] = true
    #add the first arc
    for pair in values[arc.assignment]
        if !visited[pair[1]] && pair[1]!=length(visited)
            rrb += pair[2]
            break
        end
    end
    #add the last arc
    for pair in lastArc
        if !visited[pair[1]]
            rrb += pair[2]
            break
        end
    end

    worstValues = Vector{Int}()
    for i in 1:length(visited)-1
        if !visited[i]
            valueFound = false
            for pair in values[i]
                if !visited[pair[1]] && pair[1]!=length(visited)
                    if isempty(worstValues)
                        push!(worstValues, pair[2])
                    else
                        insert!(
                            worstValues,
                            searchsorted(worstValues, pair[2]).start,
                            pair[2]
                        )
                    end
                    valueFound = true
                    break
                end
            end
            if !valueFound
                if isempty(worstValues)
                    push!(worstValues, 0)
                else
                    insert!(
                        worstValues,
                        searchsorted(worstValues, 0).start,
                        0
                    )
                end
            end
        end
    end
    worstValues = worstValues[1:(length(visited)-length(arc.destination)-2)]
    rrb += sum(worstValues)
    return rrb
end


text:
module QGL
using Compat

include("yaml.jl")
include("config.jl")
include("channels.jl")
include("pulses.jl")
include("pulseshapes.jl")
include("controlflow.jl")
include("APS2.jl")
include("compiler.jl")
include("sequences.jl")

end


text:
@testset "linsolve" begin
    function rosenbrock()
        function f!(fvec, x)
            fvec[1] = 1 - x[1]
            fvec[2] = 10(x[2]-x[1]^2)
        end
        function j!(fjac, x)
            fjac[1,1] = -1
            fjac[1,2] = 0
            fjac[2,1] = -20x[1]
            fjac[2,2] = 10
        end
        (OnceDifferentiable(f!, j!, [-1.2, 1.0], [-1.2, 1.0]), [-1.2; 1.0], "Rosenbrock")
    end
    prob = rosenbrock()
    df = prob[1]
    res_default_solve = nlsolve(df, prob[2]; method=:newton)
    prob = rosenbrock()
    df = prob[1]
    res_gmres_solve = nlsolve(df, prob[2]; method=:newton, linsolve = (x, A, b)->gmres!(fill!(x, eltype(x)(0)), A, b))
    @test res_gmres_solve.zero ≈ res_default_solve.zero
    prob = rosenbrock()
    df = prob[1]
    res_idrs_solve = nlsolve(df, prob[2]; method=:newton, linsolve = (x, A, b)->idrs!(fill!(x, eltype(x)(0)), A, b))
    @test res_idrs_solve.zero ≈ res_default_solve.zero
    testy = [false]
    function loud_solve(x, A, b)
        testy[1] = true
        gmres!(x, A, b)
    end
    prob = rosenbrock()
    df = prob[1]
    nlsolve(df, prob[2]; method=:newton, linsolve = loud_solve)
    @test testy[1] == true
end


text:
module Dungalyze

# package code goes here
# using StaticArrays, Dierckx, CSV, DataFrames, MAT, MATLAB, CoordinateTransformations, Plots, Unitful, Rematch, Distributions, Missings
# include(joinpath(Pkg.dir("Dungalyze"), "src", "calibrations.jl"))
# include(joinpath(Pkg.dir("Dungalyze"), "src", "turning_point.jl"))

using LinearAlgebra, SparseArrays
include("/home/yakir/.juliaold/v0.6/Dungalyze/src/types.jl")
include("/home/yakir/.juliaold/v0.6/Dungalyze/src/plots.jl")
include("/home/yakir/.juliaold/v0.6/Dungalyze/src/checks.jl")

# include(joinpath(Pkg.dir("Dungalyze"), "src", "types.jl"))
# include(joinpath(Pkg.dir("Dungalyze"), "src", "plots.jl"))
# include(joinpath(Pkg.dir("Dungalyze"), "src", "checks.jl"))

home = "/home/yakir/google-drive/.shared/Coffee beetle experiments"
figfolder = "dungalyze figures"
if !isdir(figfolder)
    mkdir(figfolder)
end

const CHECK = true
for experiment in readdir(home)
    # if !r"elin"i(experiment)
        path = joinpath(home, experiment)
        a = Experiment(path)
        plot(a)
    # end
end

#=df = DataFrame(Experiment = "", Run = "", Feeder_y = 0.0, Turning_point_x = 0.0, Turning_point_y = 0.0, Pellet_x = 0.0, Pellet_y = 0.0, Track_end_x = 0.0, Track_end_y = 0.0, ORG_nest_x = 0.0, ORG_nest_y = 0.0, Gravity_turn2pellet_x = 0.0, Gravity_turn2pellet_y = 0.0, Gravity_pellet2end_x = 0.0, Gravity_pellet2end_y = 0.0, Gravity_x = 0.0, Gravity_y = 0.0, Dist2turn = 0.0, Dis2pellet = 0.0, Dis2end = 0.0, Dis2org_nest = 0.0, Dis2grav_turn2pellet = 0.0, Dis2grav_pellet2end = 0.0, Dis2gravity = 0.0)
allowmissing!(df)
function doall(experiment::String; bad::Vector{String} = String[])
    _, _experiment = splitdir(experiment)
    mkdir(_experiment)
    cd(_experiment)
    ds = Final[]
    runs = get_runs(experiment, bad=bad)
    for r in runs
        d = Final(r)
        push!(df, row(d))
        push!(ds, d)
    end
    plot(ds)
    cd("..")
end
for experiment in readdir(home)
    bad = problems[problems[:Experiment] .== experiment, :Run]
    doall(joinpath(home, experiment), bad = bad)
end
deleterows!(df, 1)
CSV.write("all.csv", df)

experiment = "/home/yakir/google-drive/.shared/Coffee beetle experiments/Transfer pellet moved Therese"

for experiment in readdir(home)
    a = Experiment(joinpath(home, experiment))
    plot(a)
end

for experiment in readdir(home)
    if r"elin"i(experiment)
        f = joinpath(home, experiment)
        # println(f)
        a = Experiment(f)
        plot(a)
    end
end

runs = get_runs(experiment)
experiment = get_experiment(runs[1])
ds = Final{typeof(experiment), <:FinalPellet}[]

for r in runs
    d = Final(r)
    push!(ds, d)
end

plot(ds)

=##=function getpois(path::String)
    b = CSV.read(joinpath(path, poi_filename), header=["File", "Number", "POI"], types=[String, Int, String])
    b[:POI] .= strip.(b[:POI])
    pois = Dict{String, Dict{String, Int}}()
    by(b, :File) do df
        pois[df[1, :File]] = Dict(df[i,:POI] => df[i,:Number] for i in 1:nrow(df))
        end
        return pois
    end=##=
    =##=function extract_raw!(data::Dict{Symbol,Tuple{Matrix{Float64},Vector{Float64}}}, pois::Dict{String, Dict{String, Int}}, path::String, file::String)
        v = pois[file]
        matopen(joinpath(path, file), "r") do i
            X, Y, T, tmp = read(i, "xdata", "ydata", "tdata", "status")
            fr = tmp["FrameRate"]
            for (poi, column) in v
                data[Symbol(poi)] = _extract_raw(X[:, column], Y[:, column], T[:, column], fr)
            end
        end
    end=##=
    function doall(experiment::String; bad::Vector{String} = String[])
        runs = get_runs(experiment, bad=bad)
        ds = Final[]
        for r in runs
            a = sort_experiment(r)
            # d = Final(a)
            b = Calibrated(a)
            # c = Standardized(b)
            d = get_final(b)
            # p1 = plot(d)
            # p2 = plotsmooth(c, d)
            m = match(r"/.shared/(.*)$", r)
            path = joinpath("/home/yakir/tmp", m.captures[1])
            # mkpath(path)
            # savefig(p1, joinpath(path, "summary.pdf"))
            # savefig(p2, joinpath(path, "smooth.pdf"))
            push!(ds, d)
        end
        =##=if length(ds) > 1
            p = plot(ds)
            m = match(r"/.shared/(.*)$", experiment)
            path = joinpath("/home/yakir/tmp", m.captures[1])
            savefig(p, joinpath(path, "summary.pdf"))
        end=##=
    end

    # end # module

    home = "/home/yakir/google-drive/.shared/Coffee beetle experiments"
    a = Dict("Awayd_closed_nest_Belén" => ["64"],
             "Closed nest_Belén" => ["24"],
             "Daway_Belén" => [""],
             "Dright_Belén" => [""],
             "Dleft_closednest_Belén" => [""],
             "Transfer_Belén" => [""])
    for (experiment, bad) in a
        println(experiment)
        doall(joinpath(home, experiment), bad = bad)
    end

    doall.(joinpath.(home, readdir(home)))

    doall(joinpath(home, "Transfer pellet moved Therese"))

    experiments = filter(x -> !(r"^Transfer pellet moved Therese$"(x)), readdir(home))
    doall.(joinpath.(home, experiments))

    home = "/home/yakir/google-drive/.shared/Coffee beetle experiments"
    doall(joinpath(home, "Transfer_Belén"))

    home = "/home/yakir/google-drive/.shared/Coffee beetle experiments"
    doall(joinpath(home, "Transfer nest Therese"))

    experiments = filter(r"Bel[e,é]n", readdir(home))
    doall.(joinpath.(home, experiments))

    experiments = filter(r"therese"i, readdir(home))
    doall.(joinpath.(home, experiments))

    experiment = joinpath(home, "Closed nest_Belén")

    home = "/home/yakir/google-drive/.shared/Coffee beetle experiments"
    experiment = joinpath(home, filter(r"therese"i, readdir(home))[6])
    runs = get_runs(experiment)
    r = runs[1]
    a = sort_experiment(r)

    b = Calibrated(a)
    d = Final(b)

    home = "/home/yakir/google-drive/.shared/Coffee beetle experiments"
    include(joinpath(Pkg.dir("Dungalyze"), "src", "checks.jl"))

    for experiment in readdir(home)
        bad = problems[problems[:Experiment] .== experiment, :Run]
        doall(joinpath(home, experiment), bad = bad)
    end

    home = "/home/yakir/google-drive/.shared/Coffee beetle experiments"
    experiment = joinpath(home, "Closed nest Therese")
    bad = [""]

    experiment_path = joinpath(home, experiment)
    all_runs = get_runs(experiment_path)
    bad_runs = joinpath.(experiment_path, problems[problems[:Experiment] .== experiment, :Run])
    runs = setdiff(all_runs, bad_runs)
    for r in runs
        a = sort_experiment(r)
        b = Calibrated(a)
        d = Final(b)=#

end # module


text:
mutable struct Layer
    name::String
    surface::Cairo.CairoSurface
    redraw::Bool
    visible::Bool
    selectable::Bool
    dummy::Bool
    clearOnReset::Bool

    function Layer(name::String, surface::Union{Cairo.CairoSurface, Nothing}=nothing; redraw::Bool=true, visible::Bool=true, selectable::Bool=true, dummy::Bool=false, clearOnReset::Bool=true)
        return new(name, surface === nothing ? CairoARGBSurface(0, 0) : surface, redraw, visible, selectable, dummy, clearOnReset)
    end
end

include("drawable_room.jl")

Base.:(==)(lhs::Layer, rhs::Layer) = lhs.name == rhs.name

const redrawingFuncs = Dict{String, Function}()

function redrawCanvas!()
    draw(canvas)
end

function getLayerByName(layers::Array{Layer, 1}, name::String, default::String="")
    for layer in layers
        if layer.name == name
            return layer
        end
    end

    if default !== nothing || default != ""
        return getLayerByName(layers, default)
    end

    return false
end

function redrawLayer!(layer::Layer)
    layer.redraw = true

    redrawCanvas!()
end

# Void redraw because layers might not be set yet in tools
redrawLayer!(layers::Array{Layer, 1}, name::String) = redrawLayer!(getLayerByName(layers, name))
redrawLayer!(::Nothing) = false

layerName(l::Layer) = l.name
layerName(l::Nothing) = ""

function layerSurfaceInvalid(layer::Layer)
    return layer.surface.ptr == C_NULL
end

function resetLayer!(layer::Layer, room::Room)
    canvasWidth, canvasHeight = Int(width(layer.surface)), Int(height(layer.surface))

    if layerSurfaceInvalid(layer) || canvasWidth != room.size[1] || canvasHeight != room.size[2]
        deleteSurface(layer.surface)
        layer.surface = CairoARGBSurface(room.size[1], room.size[2])
    end

    if layer.clearOnReset
        clearSurface(getSurfaceContext(layer.surface))
    end
end

resetLayer!(layer::Layer, room::DrawableRoom) = resetLayer!(layer, room.room)

const globalLayerVisibility = Dict{String, Bool}()

function setGlobalLayerVisiblity(canvas::Gtk.GtkCanvas, name::String, visible::Bool=true)
    if get!(globalLayerVisibility, name, true) != visible
        globalLayerVisibility[name] = visible

        # Make sure all relevant layers are marked for redraw
        for (map, data) in drawableRooms
            for (room, dr) in data
                layer = getLayerByName(dr.layers, name)

                layer.redraw = true
            end
        end

        draw(canvas)
    end
end

function setGlobalLayerVisiblity(name::String, widget::Gtk.GtkCheckMenuItem)
    setGlobalLayerVisiblity(canvas, name, get_gtk_property(widget, :active, Bool))
end

function combineLayers!(ctx::Cairo.CairoContext, layers::Array{Layer, 1}, camera::Camera, room::DrawableRoom; alpha=nothing)
    for layer in layers
        if !layer.dummy && (layerSurfaceInvalid(layer) || layer.redraw)
            resetLayer!(layer, room)

            redrawFunc = get(redrawingFuncs, layer.name, (layer, room, camera) -> true)

            @catchall begin
                redrawFunc(layer, room, camera)

                layer.redraw = false
            end
        end

        if layer.visible && !layer.dummy && get!(globalLayerVisibility, layer.name, true)
            applyLayer!(ctx, layer, alpha=alpha)
        end
    end
end

function applyLayer!(ctx::Cairo.CairoContext, layer::Layer, x=0, y=0; alpha=nothing)
    drawImage(ctx, layer.surface, x, y, alpha=alpha)
end

max_stars_repo_name: asdfsdfadfasfasd/Ahorn

max_issues_repo_name: asdfsdfadfasfasd/Ahorn

max_forks_repo_name: asdfsdfadfasfasd/Ahorn


text:
using ArgParse,Match
using DataFrames,CSV,Pipe,CategoricalArrays

struct Yvariable
    label::String
    valtype::Type
end

function load_data(xfile, yfile, yvars::Array{Yvariable,1};idpair::Pair=1 => :id) where T <: String
    xrawdf = CSV.File(xfile) |> DataFrame |> df -> rename(df, idpair)
    yrawdf = @pipe CSV.File(yfile) |> DataFrame |> df -> rename(df, idpair) |> df -> select(df, ["id",getfield.(yvars, :label)...])

    merged_df = innerjoin(xrawdf, yrawdf, on=:id)

    noncatlabels::Array{String,1} = @pipe yvars |> filter(yc -> yc.valtype != CategoricalArray, _) |> getfield.(_, :label)

    xdf = select(merged_df,names(xrawdf)) |> selectNumerical
    yrawdf = select(merged_df,names(yrawdf))

    categorical!(yrawdf, Not(["id",noncatlabels...]))

    continousNames = selectNumerical(yrawdf) |> names
    ycontinous = select(yrawdf, ["id",continousNames...])

    categoricalNames = select(yrawdf, eltype.(eachcol(yrawdf)) .<: Union{Missing, CategoricalValue}) |> names

    ycategoricals::Array{DataFrame,1} = @pipe select(yrawdf, categoricalNames) |> # create a new dataframe only containing categorical columns
        onehot.(eachcol(_)) |>  # transform each column to a dataframe with onehot values
        insertcols!.(_, 1, :id => yrawdf[:,:id]) # insert column with row identifierto each dataframe

    ydf = undef
    if hasCategorical(yvars) && hasContinous(yvars)
        ydf = innerjoin(ycontinous, ycategoricals..., on=:id)
    elseif hasCategorical(yvars)
        ydf = ycategoricals |> first
    else
        ydf = ycontinous
    end

    (xdf = xdf, ydf = ydf)
end

function hasContinous(ycols)
    filter(ycol -> ycol.valtype <: Number, ycols) |> !isempty
end

function hasCategorical(ycols)
    filter(ycol -> ycol.valtype == CategoricalArray, ycols) |> !isempty
end

function todaframe(df, colname)
    @pipe onehot(df[:,colname]) |> insertcols!(_, 1, :id => df[:,:id])
end

"""
$(FUNCTIONNAME)(modelfile::String,xfile::String, outfile::String)

    Loads model from jld2 file, predicts using xfile and exports residual matrix into .csv file
"""
function predict_xres(modelfile::String, xfile::String, outfile::String)::DataFrame
    pls, stdevs, means, modelvariables, transformations = loadmodel(modelfile)

    #----------------------------------------------------
    #TODO: implement proper handling of transformations
    #= transformations = getTransformations(pls)
    push!(transformations,docenter(means))
    push!(transformations,douvscale(stdevs))

    [transform(preddataset) for transform in transformations] =#
    #----------------------------------------------------

    douv = "uv" in transformations

    outfileext = splitext(outfile)|> last
    delim = outfileext == ".csv" ? "," : "\t"

    xdf = CSV.File(xfile) |> DataFrame
    preddataset = @pipe xdf |>
        select(_,Symbol.(modelvariables))|>
        parseDataFrame(_) |>
        normalize!(_, doscale=douv, stdevs=stdevs, means=means)

    YpredVar, Xres = predictY(preddataset, pls)

    outdf::DataFrame = @pipe DataFrame(Xres) |> rename!(_, Symbol.(preddataset.value_columns)) |> hcat(xdf[:,[1]],_)

    CSV.write(outfile, outdf, delim = delim)

    outdf
end

function docenter(means)
    dataset->dataset.X .-= means'
end

function douvscale(stdevs)
    dataset->dataset.X ./= stdevs'
end

function splitArrayArgs(parsed_args, field)
    if !isnothing(parsed_args[field])
        return split(parsed_args[field], ";") |> a -> convert(Array{String,1}, a) |> strarr -> filter(str -> length(str) > 1, strarr)
    end

    return []
end

"""
$(FUNCTIONNAME)(xdf::DataFrame, ydf::DataFrame, A::Int64, modelfile::String; filters = [dataset -> dataset.mvs .< 0.25])

    Calibrates PLS model based on datatypes in DataFrame for y

    Columns of type CategoricalArray is handled by one-hot precedure

    The calibrated model is saved to specified locations

    Variables with more than 25% missing values are excluded by default
"""
function calibrate_model(xdf::DataFrame, ydf::DataFrame, A::Int64, modelfile::String; filters = [dataset -> dataset.mvs .< 0.25], transformations=["center","uv"])

    uvscale = "uv" in transformations

    # allow some missingness in x data
    xobsfilters = [ds -> calcObsMissingValues(ds) .<= 0.5]

    # do not allow any missing values observations in y dataset
    yobsfilters = [ds -> calcObsMissingValues(ds) .== 0]

    # parse datasets to dataframes
    xdataset = xdf |> parseDataFrame
    ydataset = ydf |> parseDataFrame

    # filter datasets for sample missingness
    xobsmask = filterObservations(xdataset,filters=xobsfilters)
    yobsmask = filterObservations(ydataset,filters=yobsfilters)

    # check overlap between datasets for inclusion of observations
    obsmask = xobsmask .& yobsmask

    # filter variables for both datasets using the same filter list
    xvarmask = filterVariables(xdataset,filters=filters)
    yvarmask = filterVariables(ydataset,filters=filters)

    # copy datasets from filter masks
    xdataset = copydataset(xdataset,obsmask,xvarmask) |> ds -> normalize(ds,doscale=uvscale)
    ydataset = copydataset(ydataset,obsmask,yvarmask) |> ds -> normalize(ds,doscale=true)

    # calculate PLS model
    pls = calcPLS(xdataset, ydataset, A)

    # save PLS model to disk
    savemodel(pls, xdataset, modelfile,transformations)

    pls
end

"""
$(FUNCTIONNAME)(parsed_args::Dict{String,Any})

    Calibrates PLS model based on datatypes in DataFrame for y

    Columns of type CategoricalArray is handled by one-hot precedure

    The calibrated model is saved to specified locations
"""
function calibrate_model(parsed_args::Dict{String,Any})::MultivariateModel

    mvcutoff = parsed_args["mvcutoff"]
    minvalues = parsed_args["minvalues"]
    uvscale = parsed_args["uvscale"]

    transformations = ["center"]
    if uvscale
        push!(transformations,"uv")
    end

    varfilters = [dataset -> dataset.mvs .<= mvcutoff, dataset -> dataset.varvalues .>= minvalues]

    yvars::Array{Yvariable,1} = Array{Yvariable,1}()

    @pipe splitArrayArgs(parsed_args, "ycategorical") |> Yvariable.(_, CategoricalArray) |> append!(yvars, _)

    @pipe splitArrayArgs(parsed_args, "ycontinous") |> Yvariable.(_, Float64) |> append!(yvars, _)

    yvars |> println

    if !isempty(yvars)
        xfile = parsed_args["xfile"]
        yfile = parsed_args["yfile"]
        modelfile = parsed_args["modelfile"]

        A = parsed_args["components"]

        xdf, ydf = load_data(xfile, yfile, yvars)

        calibrate_model(xdf, ydf, A, modelfile,filters=varfilters, transformations = transformations)
    end
end

"""
$(FUNCTIONNAME)(model::T, dataset::Dataset, name::String) where T <: MultivariateModel

    Save PCA or PLS model as JLD2 file
"""
function correct(parsed_args)
    xfile = parsed_args["xfile"]
    modelfile = parsed_args["modelfile"]
    outfile = parsed_args["outfile"]

    predict_xres(modelfile, xfile, outfile)
end


text:
"Diffusion MRI modelling"
module SpinDoctor

using Arpack: eigs
using DiffEqCallbacks: PresetTimeCallback
using Expokit: expmv!
# using GLMakie
using CairoMakie
using GLPK: Optimizer
using LinearAlgebra
using MiniQhull: delaunay
using OrdinaryDiffEq
using Polynomials: fit
using Printf
using QuadGK
using Roots: find_zeros
using SparseArrays
using SpecialFunctions
using StaticArrays
using Statistics: mean
using TetGen: RawTetGenIO, facetlist!, tetrahedralize
using Triangulate: TriangulateIO, triangulate
using UnPack
using WriteVTK

export TimeProfile,
    PGSE,
    CosOGSE,
    SinOGSE,
    DoublePGSE,
    integral,
    bvalue_no_q,
    intervals,
    get_interval,
    constant_intervals,
    is_constant,
    echotime
export Setup, CylinderSetup, SphereSetup, NeuronSetup, Experiment
export get_coefficients
export assemble_matrices
export create_cells
export create_surfaces_cylinder
export create_surfaces_sphere
export create_surfaces_neuron
export create_geometry
export create_model
export split_mesh
export create_directions
export get_cmpt_volumes, get_mesh_volumes
export compute_laplace_eig
export length2eig
export eig2length
export Model
export savefield
export savefield_time
export save_btpde_results
export solve_analytical
export solve_btpde
export solve_btpde_midpoint
export solve_hadc
export solve_mf
export solve_karger
export fit_adc
export fit_tensors
export plot_mesh
export get_values

# Rexport default ODE solver
export QNDF, MagnusGL6

# Datatypes
include("datatypes/sequences.jl")
include("datatypes/experiment.jl")
include("datatypes/femesh.jl")
include("datatypes/gradient.jl")
include("datatypes/model.jl")
include("datatypes/create_model.jl")
include("datatypes/setup.jl")
include("datatypes/get_coefficients.jl")

# Geometry
include("geometry/call_tetgen.jl")
include("geometry/convexhull.jl")
include("geometry/create_cells.jl")
include("geometry/create_directions.jl")
include("geometry/create_fibonacci_sphere.jl")
include("geometry/create_geometry.jl")
include("geometry/create_surfaces_cylinder.jl")
include("geometry/create_surfaces_neuron.jl")
include("geometry/create_surfaces_sphere.jl")
include("geometry/deform_domain.jl")
include("geometry/get_volumes.jl")
include("geometry/get_mesh_surface.jl")
include("geometry/get_mesh_surfacenormals.jl")
include("geometry/gmesh2fem.jl")
include("geometry/read_cells.jl")
include("geometry/read_surfaces.jl")
include("geometry/read_tetgen.jl")
include("geometry/save_cells.jl")
include("geometry/save_surfaces.jl")
include("geometry/save_tetgen.jl")
include("geometry/split_mesh.jl")

# Matrix assembly
include("matrix_assembly/assemble_mass_matrix.jl")
include("matrix_assembly/assemble_stiffness_matrix.jl")
include("matrix_assembly/assemble_flux_matrices.jl")
include("matrix_assembly/assemble_flux_matrix.jl")
include("matrix_assembly/couple_flux_matrix.jl")
include("matrix_assembly/assemble_matrices.jl")

# Matrix formalism
include("matrix_formalism/eig2length.jl")
include("matrix_formalism/compute_laplace_eig.jl")
include("matrix_formalism/solve_mf.jl")

# Analytical
include("analytical/alpha_func.jl")
include("analytical/compute_bc.jl")
include("analytical/compute_beta.jl")
include("analytical/compute_int_I.jl")
include("analytical/compute_int_J.jl")
include("analytical/compute_int_K.jl")
include("analytical/compute_JY.jl")
include("analytical/compute_nu.jl")
include("analytical/compute_v.jl")
include("analytical/find_alpha.jl")
include("analytical/solve_analytical.jl")

# BTPDE
include("btpde/solve_btpde.jl")
include("btpde/solve_btpde_midpoint.jl")
include("btpde/solve_karger.jl")

# ADC
include("adc/solve_hadc.jl")

# Postprocessing
include("postprocess/fit_adc.jl")
include("postprocess/fit_tensors.jl")
include("postprocess/savefield.jl")

# Plot
include("plot/plot_mesh.jl")

# Utils
include("utils/get_values.jl")

end


text:
function Base.:^(a::UncertainDataset, b::UncertainDataset; n = 10000)
    N = length(a)
    n_vals_b = length(b)

    if N != n_vals_b
        throw(ArgumentError("Dataset lengths do not match ($N, $n_vals_b)"))
    end

    UncertainDataset([a[i] ^ b[i] for i = 1:N])
end

function Base.:^(a::T, b::UncertainDataset; n = 10000) where {T <: Number}
    N = length(a)
    UncertainDataset([a ^ b[i] for i = 1:N])
end

function Base.:^(a::Vector{T}, b::UncertainDataset; n = 10000) where {T <: Number}
    N = length(a)
    n_vals_b = length(b)

    if N != n_vals_b
        throw(ArgumentError("Dataset lengths do not match ($N, $n_vals_b)"))
    end

    UncertainDataset([a[i] ^ b[i] for i = 1:N])
end

function Base.:^(a::UncertainDataset, b::T; n = 10000) where {T <: Number}
    N = length(a)

    UncertainDataset([a[i] ^ b for i = 1:N])
end

function Base.:^(a::UncertainDataset, b::Vector{T}; n = 10000) where {T <: Number}
    N = length(a)
    n_vals_b = length(b)

    if N != n_vals_b
        throw(ArgumentError("Dataset lengths do not match ($N, $n_vals_b)"))
    end

    UncertainDataset([a[i] ^ b[i] for i = 1:N])
end


text:
using TimerOutputs

@params struct LinearElasticityResult
    comp
    u
end
GPUUtils.whichdevice(r::LinearElasticityResult) = whichdevice(r.u)

function simulate(problem::StiffnessTopOptProblem, topology = ones(getncells(TopOptProblems.getdh(problem).grid)); round = true, hard = true, xmin = 0.001)
    if round
        if hard
            solver = FEASolver(Displacement, Direct, problem, xmin = 0.0)
        else
            solver = FEASolver(Displacement, Direct, problem, xmin = xmin)
        end
        if solver.vars isa GPUArrays.GPUArray
            solver.vars = map(round, typeof(solver.vars)(topology))
        else
            solver.vars .= Base.round.(topology)
        end
    else
        solver = FEASolver(Displacement, Direct, problem, xmin = xmin)
        copyto!(solver.vars, topology)
    end

    solver(Val{true})
    comp = dot(solver.u, solver.globalinfo.f)
    return LinearElasticityResult(comp, copy(solver.u))
end


text:
using OwnTime
using Profile

function profile_me()
    A = rand(200, 200, 400)
    maximum(A)
end

Profile.init(1_000_000, 0.001)
Profile.clear()
@profile profile_me()
owntime()
t = totaltime()
framecounts(t)
frametotal(t)
frames(t)

max_stars_repo_name: Arkoniak/OwnTime

max_issues_repo_name: Arkoniak/OwnTime

max_forks_repo_name: Arkoniak/OwnTime


text:
function _test_plain_print_methods(list::Array)
    for object in list
        normal = repr(object)
        io = IOBuffer()
        show(io, "text/plain", object)
        grabbed = String(take!(io))
        @test grabbed !== nothing
    end
end

function _test_html_print_methods(list::Array)
    for object in list
        normal = repr(object)
        io = IOBuffer()
        show(io, "text/html", object)
        grabbed = String(take!(io))
        @test grabbed !== nothing
    end
end

@testset "Test Model Print Methods" begin
    template = get_thermal_dispatch_template_network()
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5")

    dm_model = DecisionModel(template, c_sys5; optimizer=GLPK_optimizer)
    @test build!(dm_model; output_dir=mktempdir(cleanup=true)) == PSI.BuildStatus.BUILT
    @test solve!(dm_model; optimizer=GLPK_optimizer) == RunStatus.SUCCESSFUL
    results = ProblemResults(dm_model)
    variables = read_variables(results)

    list = [
        template,
        dm_model,
        PSI.get_model(template, ThermalStandard),
        PSI.get_network_model(template),
        results,
    ]

    _test_plain_print_methods(list)
    _test_html_print_methods(list)
end

@testset "Test Simulation Print Methods" begin
    template_uc = get_template_basic_uc_simulation()
    template_ed = get_template_nomin_ed_simulation()
    set_device_model!(template_ed, InterruptibleLoad, StaticPowerLoad)
    set_device_model!(template_ed, HydroEnergyReservoir, HydroDispatchReservoirBudget)
    set_network_model!(template_uc, NetworkModel(
        CopperPlatePowerModel,
        # MILP "duals" not supported with free solvers
        # duals = [CopperPlateBalanceConstraint],
    ))
    set_network_model!(
        template_ed,
        NetworkModel(
            CopperPlatePowerModel,
            duals=[CopperPlateBalanceConstraint],
            use_slacks=true,
        ),
    )
    c_sys5_hy_uc = PSB.build_system(PSITestSystems, "c_sys5_hy_uc")
    c_sys5_hy_ed = PSB.build_system(PSITestSystems, "c_sys5_hy_ed")

    models = SimulationModels(
        decision_models=[
            DecisionModel(template_uc, c_sys5_hy_uc; name="UC", optimizer=HiGHS_optimizer),
            DecisionModel(template_ed, c_sys5_hy_ed; name="ED", optimizer=HiGHS_optimizer),
        ],
    )

    sequence = SimulationSequence(
        models=models,
        feedforwards=Dict(
            "ED" => [
                SemiContinuousFeedforward(
                    component_type=ThermalStandard,
                    source=OnVariable,
                    affected_values=[ActivePowerVariable],
                ),
                EnergyLimitFeedforward(
                    component_type=HydroEnergyReservoir,
                    source=ActivePowerVariable,
                    affected_values=[ActivePowerVariable],
                    number_of_periods=12,
                ),
            ],
        ),
        ini_cond_chronology=InterProblemChronology(),
    )

    sim_not_built = Simulation(
        name="printing_sim",
        steps=2,
        models=models,
        sequence=sequence,
        simulation_folder=mktempdir(cleanup=true),
    )

    sim = Simulation(
        name="printing_sim",
        steps=2,
        models=models,
        sequence=sequence,
        simulation_folder=mktempdir(cleanup=true),
    )

    build!(sim)
    execute!(sim)
    results = SimulationResults(sim)
    results_uc = get_decision_problem_results(results, "UC")
    list = [models, sequence, template_uc, template_ed, sim, sim_not_built, results_uc]
    _test_plain_print_methods(list)
    _test_html_print_methods(list)
end