text:
module ClientAI
  import Base: position

  include("config.jl")
  include("hexgrid.jl")
  include("bot.jl")
  include("actions.jl")
  include("events.jl")
  include("ai.jl")

  macro checked(ex)
    return quote
      try
        $(esc(ex))
      catch e
        showerror(STDERR, e)
        Base.show_backtrace(STDOUT, catch_backtrace())
      end
    end
  end
  export @checked

	export AbstractAI

  # general functions and types
  export AbstractAI, Position, Config
  export botid, position, to_dict

  # events
  export AbstractEvent, HitEvent, DeathEvent, SightEvent, RadarEvent, DetectionEvent, DamageEvent, MoveEvent, NoActionEvent
  export on_event, event_dispatch

  # actions
  export AbstractAction, MoveAction, RadarAction, CannonAction
  export make_action, plan_actions, best_action, randomize

  # bot functions
  export AbstractBot
  export is_alive, team, name, hitpoints, filter_valid

  # grid functions
  export distance, radius, center, circle
  export view_area, radar_area, move_area, damage_area, map_area
end

max_stars_repo_name: Condimentum/spacetyckiting

max_issues_repo_name: Condimentum/spacetyckiting

max_forks_repo_name: Condimentum/spacetyckiting


text:
n = 10 # rand(1:10)
@test matrixdepot("lehmer", Float64, n) == matrixdepot("lehmer", n)

A = ones(n, 1) * [1:n;]'
A = A./A'
A = tril(A) + tril(A)' - Matrix(1.0I, n, n)

@test matrixdepot("lehmer", n) ≈ A
println("'lehmer' passed test...")

ext: jl


text:
using BLASBenchmarksGPU
using Test

import CUDA
import InteractiveUtils
import LinearAlgebra
import Statistics

InteractiveUtils.versioninfo(stdout)
LinearAlgebra.versioninfo(stdout)
CUDA.versioninfo(stdout)
@info "" [CUDA.capability(dev) for dev in CUDA.devices()]

@testset "CUDA: Float32=Float16×Float16" begin
    backend = :CUDA
    TA = Float16
    TB = Float16
    TC = Float32
    sizes = [2^n for n in 7:9]
    @info "" backend TA TB TC sizes
    bench_result = BLASBenchmarksGPU.runbench(
        backend,
        TA,
        TB,
        TC;
        sizes = sizes,
    )
    @test bench_result isa BLASBenchmarksGPU.BenchmarkResult
    plot_directory = mktempdir()
    plot_filename = joinpath(plot_directory, "plot.png")
    BLASBenchmarksGPU.plotbench(
        bench_result,
        plot_filename,
    )
end


text:
using Test

@testset "Hypothesis tests" begin include("test_hypothesis_tests.jl") end
@testset "FSSD" begin include("fssd_tests.jl") end
@testset "KSD" begin include("ksd_tests.jl") end
@testset "Univariate distributions" begin include("test_univariate.jl") end


text:
# Manifold interface: every manifold (subtype of Manifold) defines the functions
# project_tangent!(m, g, x): project g on the tangent space to m at x
# retract!(m, x): map x back to a point on the manifold m
# out-of-place versions are defined by default

# For mathematical references, see e.g.

# The Geometry of Algorithms with Orthogonality Constraints
# Alan Edelman, Tomás A. Arias, and Steven T. Smith
# SIAM. J. Matrix Anal. & Appl., 20(2), 303–353. (51 pages)

# Optimization Algorithms on Matrix Manifolds
# P.-A. Absil, R. Mahony, R. Sepulchre
# Princeton University Press, 2008

module Manifolds

export Manifold,
       project_tangent, project_tangent!, retract, retract!,
       Flat, Sphere, Stiefel, PowerManifold, ProductManifold

abstract type Manifold
end

# fallback for out-of-place ops
project_tangent(M::Manifold,g,x) = project_tangent!(M, copy(x), x)
retract(M::Manifold,x) = retract!(M, copy(x))

"""Flat Euclidean space {R,C}^N, with projections equal to the identity."""
struct Flat <: Manifold
end
# all the functions below are no-ops, and therefore the generated code
# for the flat manifold should be exactly the same as the one with all
# the manifold stuff removed
retract(M::Flat, x) = x
retract!(M::Flat,x) = x
project_tangent(M::Flat, g, x) = g
project_tangent!(M::Flat, g, x) = g

"""Spherical manifold {|x| = 1}."""
struct Sphere <: Manifold
end
retract!(S::Sphere, x) = normalize!(x)
project_tangent!(S::Sphere,g,x) = (g .= g .- real(vecdot(x,g)).*x)

"""
N x n matrices with orthonormal columns, i.e. such that X'X = I.
Special cases: N x 1 = sphere, N x N = orthogonal/unitary group.
Stiefel() uses a SVD algorithm to compute the retraction. To use a Cholesky-based orthogonalization, use Stiefel_CholQR()
"""
abstract type Stiefel <: Manifold end
struct Stiefel_CholQR <: Stiefel end
struct Stiefel_SVD <: Stiefel end
function Stiefel(retraction=:SVD)
    if retraction == :CholQR
        Stiefel_CholQR()
    elseif retraction == :SVD
        Stiefel_SVD()
    else
        error("Unsupported argument.")
    end
end
function retract!(S::Stiefel_SVD, X)
    U,S,V = svd(X)
    X .= U*V'
end
function retract!(S::Stiefel_CholQR, X)
    overlap = X'X
    X .= X/chol(overlap)
end
project_tangent!(S::Stiefel, G, X) = (G .-= X*((X'G .+ G'X)./2))

"""
Multiple copies of the same manifold. Points are stored as inner_dims x outer_dims,
e.g. the product of 2x2 Stiefel manifolds of dimension N x n would be a N x n x 2 x 2 matrix.
"""
struct PowerManifold<:Manifold
    "Type of embedded manifold"
    inner_manifold::Manifold
    "Dimension of the embedded manifolds"
    inner_dims::Tuple
    "Number of embedded manifolds"
    outer_dims::Tuple
end
function retract!(m::PowerManifold, x)
    for i=1:prod(m.outer_dims)
        retract!(m.inner_manifold,get_inner(m, x, i))
    end
    x
end
function project_tangent!(m::PowerManifold, g, x)
    for i=1:prod(m.outer_dims)
        project_tangent!(m.inner_manifold,get_inner(m, g, i),get_inner(m, x, i))
    end
    g
end
@inline function get_inner(m::PowerManifold, x, i::Int)
    size_inner = prod(m.inner_dims)
    size_outer = prod(m.outer_dims)
    @assert 1 <= i <= size_outer
    return reshape(view(x, (i-1)*size_inner+1:i*size_inner), m.inner_dims)
end
@inline get_inner(m::PowerManifold, x, i::Tuple) = get_inner(m, x, ind2sub(m.outer_dims, i...))

"""
Product of two manifolds {P = (x1,x2), x1 ∈ m1, x2 ∈ m2}.
P is stored as a flat 1D array, and x1 is before x2 in memory.
Use get_inner(m, x, {1,2}) to access x1 and x2 in their original format.
"""
struct ProductManifold<:Manifold
    m1::Manifold
    m2::Manifold
    dims1::Tuple
    dims2::Tuple
end
function retract!(m::ProductManifold, x)
    retract!(m.m1, get_inner(m,x,1))
    retract!(m.m2, get_inner(m,x,2))
    x
end
function project_tangent!(m::ProductManifold, g, x)
    project_tangent!(m.m1, get_inner(m, g, 1), get_inner(m, x, 1))
    project_tangent!(m.m2, get_inner(m, g, 2), get_inner(m, x, 2))
    g
end
function get_inner(m::ProductManifold, x, i)
    N1 = prod(m.dims1)
    N2 = prod(m.dims2)
    @assert length(x) == N1+N2
    if i == 1
        return reshape(view(x, 1:N1),m.dims1)
    elseif i == 2
        return reshape(view(x, N1+1:N1+N2), m.dims2)
    else
        error("Only two components in a product manifold")
    end
end
end


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: mediapackage_vod
using AWS.Compat
using AWS.UUIDs

"""
    configure_logs(id)
    configure_logs(id, params::Dict{String,<:Any})

Changes the packaging group's properities to configure log subscription

# Arguments
- `id`: The ID of a MediaPackage VOD PackagingGroup resource.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"egressAccessLogs"`:
"""
configure_logs(id; aws_config::AbstractAWSConfig=global_aws_config()) = mediapackage_vod("PUT", "/packaging_groups/$(id)/configure_logs"; aws_config=aws_config)
configure_logs(id, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = mediapackage_vod("PUT", "/packaging_groups/$(id)/configure_logs", params; aws_config=aws_config)

"""
    create_asset(id, packaging_group_id, source_arn, source_role_arn)
    create_asset(id, packaging_group_id, source_arn, source_role_arn, params::Dict{String,<:Any})

Creates a new MediaPackage VOD Asset resource.

# Arguments
- `id`: The unique identifier for the Asset.
- `packaging_group_id`: The ID of the PackagingGroup for the Asset.
- `source_arn`: ARN of the source object in S3.
- `source_role_arn`: The IAM role ARN used to access the source S3 bucket.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"resourceId"`: The resource ID to include in SPEKE key requests.
- `"tags"`:
"""
create_asset(id, packagingGroupId, sourceArn, sourceRoleArn; aws_config::AbstractAWSConfig=global_aws_config()) = mediapackage_vod("POST", "/assets", Dict{String, Any}("id"=>id, "packagingGroupId"=>packagingGroupId, "sourceArn"=>sourceArn, "sourceRoleArn"=>sourceRoleArn); aws_config=aws_config)
create_asset(id, packagingGroupId, sourceArn, sourceRoleArn, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = mediapackage_vod("POST", "/assets", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("id"=>id, "packagingGroupId"=>packagingGroupId, "sourceArn"=>sourceArn, "sourceRoleArn"=>sourceRoleArn), params)); aws_config=aws_config)

"""
    create_packaging_configuration(id, packaging_group_id)
    create_packaging_configuration(id, packaging_group_id, params::Dict{String,<:Any})

Creates a new MediaPackage VOD PackagingConfiguration resource.

# Arguments
- `id`: The ID of the PackagingConfiguration.
- `packaging_group_id`: The ID of a PackagingGroup.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"cmafPackage"`:
- `"dashPackage"`:
- `"hlsPackage"`:
- `"mssPackage"`:
- `"tags"`:
"""
create_packaging_configuration(id, packagingGroupId; aws_config::AbstractAWSConfig=global_aws_config()) = mediapackage_vod("POST", "/packaging_configurations", Dict{String, Any}("id"=>id, "packagingGroupId"=>packagingGroupId); aws_config=aws_config)
create_packaging_configuration(id, packagingGroupId, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = mediapackage_vod("POST", "/packaging_configurations", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("id"=>id, "packagingGroupId"=>packagingGroupId), params)); aws_config=aws_config)

"""
    create_packaging_group(id)
    create_packaging_group(id, params::Dict{String,<:Any})

Creates a new MediaPackage VOD PackagingGroup resource.

# Arguments
- `id`: The ID of the PackagingGroup.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"authorization"`:
- `"egressAccessLogs"`:
- `"tags"`:
"""
create_packaging_group(id; aws_config::AbstractAWSConfig=global_aws_config()) = mediapackage_vod("POST", "/packaging_groups", Dict{String, Any}("id"=>id); aws_config=aws_config)
create_packaging_group(id, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = mediapackage_vod("POST", "/packaging_groups", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("id"=>id), params)); aws_config=aws_config)

"""
    delete_asset(id)
    delete_asset(id, params::Dict{String,<:Any})

Deletes an existing MediaPackage VOD Asset resource.

# Arguments
- `id`: The ID of the MediaPackage VOD Asset resource to delete.

"""
delete_asset(id; aws_config::AbstractAWSConfig=global_aws_config()) = mediapackage_vod("DELETE", "/assets/$(id)"; aws_config=aws_config)
delete_asset(id, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = mediapackage_vod("DELETE", "/assets/$(id)", params; aws_config=aws_config)

"""
    delete_packaging_configuration(id)
    delete_packaging_configuration(id, params::Dict{String,<:Any})

Deletes a MediaPackage VOD PackagingConfiguration resource.

# Arguments
- `id`: The ID of the MediaPackage VOD PackagingConfiguration resource to delete.

"""
delete_packaging_configuration(id; aws_config::AbstractAWSConfig=global_aws_config()) = mediapackage_vod("DELETE", "/packaging_configurations/$(id)"; aws_config=aws_config)
delete_packaging_configuration(id, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = mediapackage_vod("DELETE", "/packaging_configurations/$(id)", params; aws_config=aws_config)

"""
    delete_packaging_group(id)
    delete_packaging_group(id, params::Dict{String,<:Any})

Deletes a MediaPackage VOD PackagingGroup resource.

# Arguments
- `id`: The ID of the MediaPackage VOD PackagingGroup resource to delete.

"""
delete_packaging_group(id; aws_config::AbstractAWSConfig=global_aws_config()) = mediapackage_vod("DELETE", "/packaging_groups/$(id)"; aws_config=aws_config)
delete_packaging_group(id, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = mediapackage_vod("DELETE", "/packaging_groups/$(id)", params; aws_config=aws_config)

"""
    describe_asset(id)
    describe_asset(id, params::Dict{String,<:Any})

Returns a description of a MediaPackage VOD Asset resource.

# Arguments
- `id`: The ID of an MediaPackage VOD Asset resource.

"""
describe_asset(id; aws_config::AbstractAWSConfig=global_aws_config()) = mediapackage_vod("GET", "/assets/$(id)"; aws_config=aws_config)
describe_asset(id, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = mediapackage_vod("GET", "/assets/$(id)", params; aws_config=aws_config)

"""
    describe_packaging_configuration(id)
    describe_packaging_configuration(id, params::Dict{String,<:Any})

Returns a description of a MediaPackage VOD PackagingConfiguration resource.

# Arguments
- `id`: The ID of a MediaPackage VOD PackagingConfiguration resource.

"""
describe_packaging_configuration(id; aws_config::AbstractAWSConfig=global_aws_config()) = mediapackage_vod("GET", "/packaging_configurations/$(id)"; aws_config=aws_config)
describe_packaging_configuration(id, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = mediapackage_vod("GET", "/packaging_configurations/$(id)", params; aws_config=aws_config)

"""
    describe_packaging_group(id)
    describe_packaging_group(id, params::Dict{String,<:Any})

Returns a description of a MediaPackage VOD PackagingGroup resource.

# Arguments
- `id`: The ID of a MediaPackage VOD PackagingGroup resource.

"""
describe_packaging_group(id; aws_config::AbstractAWSConfig=global_aws_config()) = mediapackage_vod("GET", "/packaging_groups/$(id)"; aws_config=aws_config)
describe_packaging_group(id, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = mediapackage_vod("GET", "/packaging_groups/$(id)", params; aws_config=aws_config)

"""
    list_assets()
    list_assets(params::Dict{String,<:Any})

Returns a collection of MediaPackage VOD Asset resources.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`: Upper bound on number of records to return.
- `"nextToken"`: A token used to resume pagination from the end of a previous request.
- `"packagingGroupId"`: Returns Assets associated with the specified PackagingGroup.
"""
list_assets(; aws_config::AbstractAWSConfig=global_aws_config()) = mediapackage_vod("GET", "/assets"; aws_config=aws_config)
list_assets(params::AbstractDict{String, Any}; aws_config::AbstractAWSConfig=global_aws_config()) = mediapackage_vod("GET", "/assets", params; aws_config=aws_config)

"""
    list_packaging_configurations()
    list_packaging_configurations(params::Dict{String,<:Any})

Returns a collection of MediaPackage VOD PackagingConfiguration resources.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`: Upper bound on number of records to return.
- `"nextToken"`: A token used to resume pagination from the end of a previous request.
- `"packagingGroupId"`: Returns MediaPackage VOD PackagingConfigurations associated with
  the specified PackagingGroup.
"""
list_packaging_configurations(; aws_config::AbstractAWSConfig=global_aws_config()) = mediapackage_vod("GET", "/packaging_configurations"; aws_config=aws_config)
list_packaging_configurations(params::AbstractDict{String, Any}; aws_config::AbstractAWSConfig=global_aws_config()) = mediapackage_vod("GET", "/packaging_configurations", params; aws_config=aws_config)

"""
    list_packaging_groups()
    list_packaging_groups(params::Dict{String,<:Any})

Returns a collection of MediaPackage VOD PackagingGroup resources.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`: Upper bound on number of records to return.
- `"nextToken"`: A token used to resume pagination from the end of a previous request.
"""
list_packaging_groups(; aws_config::AbstractAWSConfig=global_aws_config()) = mediapackage_vod("GET", "/packaging_groups"; aws_config=aws_config)
list_packaging_groups(params::AbstractDict{String, Any}; aws_config::AbstractAWSConfig=global_aws_config()) = mediapackage_vod("GET", "/packaging_groups", params; aws_config=aws_config)

"""
    list_tags_for_resource(resource-arn)
    list_tags_for_resource(resource-arn, params::Dict{String,<:Any})

Returns a list of the tags assigned to the specified resource.

# Arguments
- `resource-arn`: The Amazon Resource Name (ARN) for the resource. You can get this from
  the response to any request to the resource.

"""
list_tags_for_resource(resource_arn; aws_config::AbstractAWSConfig=global_aws_config()) = mediapackage_vod("GET", "/tags/$(resource-arn)"; aws_config=aws_config)
list_tags_for_resource(resource_arn, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = mediapackage_vod("GET", "/tags/$(resource-arn)", params; aws_config=aws_config)

"""
    tag_resource(resource-arn, tags)
    tag_resource(resource-arn, tags, params::Dict{String,<:Any})

Adds tags to the specified resource. You can specify one or more tags to add.

# Arguments
- `resource-arn`: The Amazon Resource Name (ARN) for the resource. You can get this from
  the response to any request to the resource.
- `tags`: A collection of tags associated with a resource

"""
tag_resource(resource_arn, tags; aws_config::AbstractAWSConfig=global_aws_config()) = mediapackage_vod("POST", "/tags/$(resource-arn)", Dict{String, Any}("tags"=>tags); aws_config=aws_config)
tag_resource(resource_arn, tags, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = mediapackage_vod("POST", "/tags/$(resource-arn)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("tags"=>tags), params)); aws_config=aws_config)

"""
    untag_resource(resource-arn, tag_keys)
    untag_resource(resource-arn, tag_keys, params::Dict{String,<:Any})

Removes tags from the specified resource. You can specify one or more tags to remove.

# Arguments
- `resource-arn`: The Amazon Resource Name (ARN) for the resource. You can get this from
  the response to any request to the resource.
- `tag_keys`: A comma-separated list of the tag keys to remove from the resource.

"""
untag_resource(resource_arn, tagKeys; aws_config::AbstractAWSConfig=global_aws_config()) = mediapackage_vod("DELETE", "/tags/$(resource-arn)", Dict{String, Any}("tagKeys"=>tagKeys); aws_config=aws_config)
untag_resource(resource_arn, tagKeys, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = mediapackage_vod("DELETE", "/tags/$(resource-arn)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("tagKeys"=>tagKeys), params)); aws_config=aws_config)

"""
    update_packaging_group(id)
    update_packaging_group(id, params::Dict{String,<:Any})

Updates a specific packaging group. You can't change the id attribute or any other
system-generated attributes.

# Arguments
- `id`: The ID of a MediaPackage VOD PackagingGroup resource.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"authorization"`:
"""
update_packaging_group(id; aws_config::AbstractAWSConfig=global_aws_config()) = mediapackage_vod("PUT", "/packaging_groups/$(id)"; aws_config=aws_config)
update_packaging_group(id, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = mediapackage_vod("PUT", "/packaging_groups/$(id)", params; aws_config=aws_config)


text:
@testset "model2fit_sum.jl" begin

    @testset "Sum of 2 Gaussian with const centers" begin

        model = Gaussian_Peak() + Gaussian_Peak()
        model_θ = Float64[1,5,2,  # h1=1,μ1=5,σ1=2
                          2,10,3] # h2=2,μ2=10,σ2=3

        X = collect(1:1.0:10)
        Y = eval_y(model,X,model_θ)

        indices_μ = [2,5]
        values_μ = [5.0,10.0]
        model_const_μ = Model2Fit_Const_Parameters(model,[2,5],[5.0,10.0])
        model_const_μ_θ = deleteat!(copy(model_θ),indices_μ)

        @test parameter_size(model_const_μ) == length(model_const_μ_θ)

        @test get_model(model_const_μ) == model
        @test get_model_θ(model_const_μ, model_const_μ_θ) == model_θ

        Y_const = eval_y(model_const_μ,X,model_const_μ_θ)
        @test Y == Y_const

    end

end


text:
# using Revise
using Test

@testset "BifurcationKit" begin

	@testset "Linear Solvers" begin
		include("precond.jl")
		include("test_linear.jl")
	end

	@testset "Newton" begin
		include("test_newton.jl")
		include("test-bordered-problem.jl")
	end

	@testset "Continuation" begin
		include("test_bif_detection.jl")
		include("test-cont-non-vector.jl")
		include("simple_continuation.jl")
		include("testNF.jl")
	end

	@testset "Events / User function" begin
		include("event.jl")
	end

	@testset "Fold Codim 2" begin
		include("testJacobianFoldDeflation.jl")
		include("codim2.jl")
	end

	@testset "Hopf Codim 2" begin
		include("testHopfMA.jl")
		include("lorenz84.jl")
		include("COModel.jl")
	end

	@testset "Periodic orbits" begin
		include("test_potrap.jl")
		include("test_SS.jl")
		include("poincareMap.jl")
		include("stuartLandauSH.jl")
		include("stuartLandauTrap.jl")
		include("stuartLandauCollocation.jl")
		include("testLure.jl")
	end

	@testset "Wave" begin
		include("test_wave.jl")
	end
end


text:
abstract type AbstractLinearITRM <: InterpretableTimeseriesRegressionModel end

hasintercept(m::AbstractLinearITRM) = m.hasintercept
getintercept(m::AbstractLinearITRM) = hasintercept(m) ? m.params[end] : 0.0
ndigits(m::AbstractLinearITRM) = m.ndigits
nfactors(m::AbstractLinearITRM) = m.nfactors
nfeatures(m::AbstractLinearITRM) = m.nfeatures
maxlag(m::AbstractLinearITRM) = m.maxlag
actualmaxlag(m::AbstractLinearITRM, ps::Nothing) = m.maxlag
minlag(m::AbstractLinearITRM) = m.minlag
actualminlag(m::AbstractLinearITRM, ps::Nothing) = m.minlag
mincoef(m::AbstractLinearITRM) = m.mincoef
maxcoef(m::AbstractLinearITRM) = m.maxcoef
r(m::AbstractLinearITRM, v::Float64) = round(v, digits = ndigits(m))
hasparams(m::AbstractLinearITRM) = !isnothing(params(m))
isvalidparams(m::AbstractLinearITRM, ps::Vector{Float64}) = true
setparams!(m::AbstractLinearITRM, ps::Vector{Float64}) = isvalidparams(m, ps) && (m.params = ps)
params(m::AbstractLinearITRM) = m.params
setfeaturenames!(m::AbstractLinearITRM, fns::Vector{String}) = (m.featurenames = fns)
featurename(m::AbstractLinearITRM, fi::Int) =
    (length(m.featurenames) > 0) ? m.featurenames[fi] : "feature_$(fi)"

function StatsBase.predict(m::AbstractLinearITRM, X::AbstractMatrix{N},
    ps::Vector{Float64} = params(m)) where {N<:Number}
    preds = Array{Float64}(undef, size(X, 1)-actualmaxlag(m, ps))
    StatsBase.predict!(m, preds, X, ps)
end

function StatsBase.predict!(m::AbstractLinearITRM, preds::Vector{Float64},
    X::AbstractMatrix{N}, ps::Vector{Float64} = params(m)) where {N<:Number}
    previ = 0
    for row in (actualmaxlag(m, ps)+1):size(X, 1)
        preds[(previ += 1)] = predictforrow(m, X, row, ps)
    end
    preds # Note that due to lag this will be shorter than num rows of X
end

function Base.show(io::IO, m::AbstractLinearITRM)
    writedescription(io, m)
    write(io, "  model: ")
    writemodel(io, m)
end

function writemodel(io::IO, m::AbstractLinearITRM)
    if !hasparams(m)
        write(io, "<no params set yet>")
        return
    end
    hasintercept(m) && write(io, string(r(m, getintercept(m))))
    for fi in 1:nfactors(m)
        writefactor(io, m, fi)
    end
end

modelname(m::AbstractLinearITRM) = "Linear ITRM (Interpretable Timeseries Regression Model)"

function writedescription(io::IO, m::AbstractLinearITRM)
    write(io, modelname(m) * "\n")
    if mincoef(m) < 0.0 && maxcoef(m) > 0.0
        write(io, "  with positive and negative coefficients")
    elseif maxcoef(m) > 0.0
        write(io, "  with positive coefficients only")
    elseif mincoef(m) < 0.0
        write(io, "  with negative coefficients only")
    end
    write(io, " (allowed range: [$(mincoef(m)), $(maxcoef(m))])\n")
    write(io, "  num factors:  $(nfactors(m))\n")
    write(io, "  num features: $(nfeatures(m))\n")
    aminl, amaxl = actualminlag(m), actualmaxlag(m)
    ar = (aminl == amaxl) ? string(aminl) : "$(aminl)-$(amaxl)"
    ars = hasparams(m) ? ", $(ar) (actual)" : ""
    write(io, "  lags:         $(minlag(m))-$(maxlag(m)) (valid range)$(ars)\n")
end

modelasstring(m::AbstractLinearITRM, df::DataFrame, ps::Vector{Float64} = m.params) =
    modelasstring(m, ps, names(df))

function modelasstring(m::AbstractLinearITRM, params::Vector{Float64} = m.params,
    featurenames::AbstractVector = String[])
    setparams!(m, params)
    length(featurenames) > 0 && setfeaturenames!(m, map(string, featurenames))
    iob = IOBuffer()
    writemodel(iob, m)
    String(take!(iob))
end

function calcfitness(m::AbstractLinearITRM, X::AbstractMatrix{N}, y::AbstractVector{N},
    ps::Vector{Float64} = params(m), lossfn::Function = l2loss) where {N<:Number}
    yhat = predict(m, X, ps)
    lossfn(yhat, align(yhat, y))
end

align(yhat, y) = y[(end-length(yhat)+1):end]

function predictandcalcloss(m::AbstractLinearITRM, X::AbstractMatrix{Float64},
            y::AbstractVector{Float64}, ps::Vector{Float64}, lossfn::Function)
    yhat = predict(m, X, ps)
    return yhat, lossfn(yhat, align(yhat, y))
end

fit!(m::AbstractLinearITRM, X::AbstractMatrix{Float64}, y::AbstractVector{Float64}; kws...) =
    fit!(m, InterpretableModels.DefaultOptimizer(m), X, y; kws...)


text:
@symbol_func function K_DR(cur_reactor::AbstractReactor)
  cur_reactor.is_pulsed || return 0.0

  cur_denominator = cur_reactor.sigma_CS

  cur_denominator *= (12/5) * cur_reactor.pi

  cur_denominator -= cur_reactor.B_CS ^ 2

  cur_numerator = 3.0

  cur_numerator *= cur_reactor.B_CS ^ 2

  cur_K = cur_numerator

  cur_K /= cur_denominator

  cur_K
end


text:
module FTPServer

import Base: Process
import Base: close
using Conda
using Memento
using PyCall
using Random: randstring
const LOGGER = getlogger(@__MODULE__)

const pylogging = PyNULL()
const pyopenssl_crypto = PyNULL()
const pyopenssl_SSL = PyNULL()
const pyftpdlib_authorizers = PyNULL()
const pyftpdlib_handlers = PyNULL()
const pyftpdlib_servers = PyNULL()

# Defaults from pyftpdlib example
const USER = "user"
const PASSWD = "12345"
const HOST = "localhost"
const PORT = 2021
const PERM = "elradfmwM"
const DEBUG = false

const SCRIPT = abspath(dirname(@__FILE__), "server.py")
const ROOT = abspath(joinpath(dirname(dirname(@__FILE__)), "deps", "usr", "ftp"))
const HOMEDIR = joinpath(ROOT, "data")
const CERT = joinpath(ROOT, "test.crt")
const KEY = joinpath(ROOT, "test.key")
const PYTHON_CMD = joinpath(
    Conda.PYTHONDIR, Sys.iswindows() ? "python.exe" : "python"
)

function __init__()
    Memento.register(LOGGER)
    copy!(pyopenssl_crypto, pyimport_conda("OpenSSL.crypto", "OpenSSL"))
    copy!(pyopenssl_SSL, pyimport_conda("OpenSSL.SSL", "OpenSSL"))
    copy!(pyftpdlib_servers, pyimport_conda("pyftpdlib.servers", "pyftpdlib", "invenia"))

    DEBUG && pylogging[:basicConfig](level=pylogging[:DEBUG])
    mkpath(ROOT)
end

mutable struct Server
    homedir::AbstractString
    port::Int
    username::AbstractString
    password::AbstractString
    permissions::AbstractString
    security::Symbol
    process::Process
    io::IO
end

"""
    Server(homedir=$HOMEDIR; username=$USER, password=$PASSWD, permissions=$PERM, security=:none)

A Server stores settings for create an pyftpdlib server.

# Arguments
- `homedir::AbstractString`: Directory where you want store to store your data for the
  test server.

# Keywords
- `username`: Default login username
- `password`: Default login password
- `permission`: Default user read/write permissions
- `security`: Security method to use for connecting (options: `:none`, `:implicity`, `:explicity`).
  Passing in `:none` will use FTP and passing in `:implicity` or `:explicity` will use the appropriate
  FTPS connection.
"""
function Server(
    homedir::AbstractString=HOMEDIR; username="", password="", permissions="elradfmwM",
    security::Symbol=:none,
)
    if isempty(username)
        username = string("user", rand(1:9999))
    end
    if isempty(password)
        password = randstring(40)
    end

    cmd = `$PYTHON_CMD $SCRIPT $username $password $homedir --permissions $permissions`
    if security != :none
        cmd = `$cmd --tls $security --cert-file $CERT --key-file $KEY --gen-certs TRUE`
    end
    io = Pipe()

    # Note: open(::AbstractCmd, ...) won't work here as it doesn't allow us to capture STDERR.
    process = run(pipeline(cmd, stdout=io, stderr=io), wait=false)

    line = readline(io)
    m = match(r"starting FTP.* server on .*:(?<port>\d+)", line)
    if m !== nothing
        port = parse(Int, m[:port])
        Server(homedir, port, username, password, permissions, security, process, io)
    else
        kill(process)
        error(line, String(readavailable(io)))  # Display traceback
    end
end

"""
    serve(f, args...; kwargs...)

Passes `args` and `kwargs` to the `Server` constructor and runs the function `f` by passing
in the `server` instance. Upon completion the `server` will automatically be shutdown.
"""
function serve(f::Function, args...; kwargs...)
    server = Server(args...; kwargs...)

    try
        f(server)
    finally
        close(server)
    end
end

hostname(server::Server) = "localhost"
port(server::Server) = server.port
username(server::Server) = server.username
password(server::Server) = server.password
close(server::Server) = kill(server.process)

localpath(server::Server, path::AbstractString) = joinpath(server.homedir, split(path, '/')...)

function tempfile(path::AbstractString)
    content = randstring(rand(1:100))
    open(path, "w") do fp
        write(fp, content)
    end
    return content
end

function setup_home(dir::AbstractString)
    mkdir(dir)
    tempfile(joinpath(dir, "test_download.txt"))
    tempfile(joinpath(dir, "test_download2.txt"))
    mkdir(joinpath(dir, "test_directory"))
end

"""
    init()

Creates a test $HOMEDIR with a few sample files if one hasn't already been setup.

```
$HOMEDIR/test_download.txt
$HOMEDIR/test_download2.txt
$HOMEDIR/test_directory/
````
"""
init() = isdir(FTPServer.HOMEDIR) || setup_home(FTPServer.HOMEDIR)

"""
    cleanup()

Cleans up the default FTPServer.ROOT directory:

- $HOMEDIR
- $CERT
- $KEY
"""
function cleanup()
    rm(FTPServer.HOMEDIR, recursive=true)
    isfile(FTPServer.CERT) && rm(FTPServer.CERT)
    isfile(FTPServer.KEY) && rm(FTPServer.KEY)
end

end # module


text:
using Plexifications
using Base.Test: @test, @test_throws

@test_throws ErrorException begin
    inv(DualPerplex(0, 0.0))
end

@test begin
    !iszerodivisor(DualPerplex(1,2,3,4))
end

@test begin
    iszerodivisor(DualPerplex(1,1,3,4))
end

@test begin
    a = Plexifications.randomBigFloat()
    isreal(DualPerplex(a))
end

@test begin
    a = 1
    !isreal(DualPerplex(a, a))
end

@test begin
    z = random(DualPerplex{Int})
    z == +(z)
end

@test begin
    a = Plexifications.randomBigFloat()
    real(DualPerplex(a)) == a
end

@test begin
    io = IOBuffer()
    show(io, DualPerplex(1, 2, 3, 4))
    l = String(take!(io))
    r = "[1: 1, R: 2, A: 3, RA: 4]"
    l == r
end

@test begin
    x = random(DualPerplex{BigFloat})
    y = random(DualPerplex{BigFloat})
    x != y
end

@test begin
    a = Plexifications.randomBigInt()
    z = random(DualPerplex{BigInt})
    l = a + z
    r = z + a
    l == r
end

@test begin
    a = Plexifications.randomBigInt()
    z = random(DualPerplex{BigInt})
    l = a - z
    r = -(z - a)
    l == r
end

@test begin
    x = random(DualPerplex{BigInt})
    y = random(DualPerplex{BigInt})
    l = x + y
    r = y + x
    l == r
end

@test begin
    x = random(Prplex{BigInt})
    y = random(DualPerplex{BigInt})
    l = x + y
    r = y + x
    l == r
end

@test begin
    x = random(DualPerplex{BigInt})
    y = random(DualPerplex{BigInt})
    l = x * y
    r = y * x
    l == r
end

@test begin
    x = random(Prplex{BigInt})
    y = random(DualPerplex{BigInt})
    l = x * y
    r = y * x
    l == r
end

@test begin
    x = random(DualPerplex{BigInt})
    y = random(DualPerplex{BigInt})
    l = x - y
    r = -(y - x)
    l == r
end

@test begin
    x = random(Prplex{BigInt})
    y = random(DualPerplex{BigInt})
    l = x - y
    r = -(y - x)
    l == r
end

@test begin
    x = random(DualPerplex{BigInt})
    y = random(DualPerplex{BigInt})
    z = random(DualPerplex{BigInt})
    l = (x + y) + z
    r = x + (y + z)
    l == r
end

@test begin
    x = random(DualPerplex{BigInt})
    y = random(DualPerplex{BigInt})
    z = random(DualPerplex{BigInt})
    l = (x * y) * z
    r = x * (y * z)
    l == r
end

@test begin
    x = random(DualPerplex{BigInt})
    barA(barA(x)) == x
end

@test begin
    x = random(DualPerplex{BigInt})
    barR(barR(x)) == x
end

@test begin
    x = random(DualPerplex{BigInt})
    barR(barA(x)) == barA(barR(x))
end

@test begin
    x = random(DualPerplex{Rational{BigInt}})
    inv(inv(x)) == x
end

@test begin
    x = random(DualPerplex{BigInt})
    y = random(DualPerplex{BigInt})
    l = barA(x * y)
    r = barA(x) * barA(y)
    l == r
end

@test begin
    x = random(DualPerplex{BigInt})
    y = random(DualPerplex{BigInt})
    l = barA(x + y)
    r = barA(x) + barA(y)
    l == r
end

@test begin
    x = random(DualPerplex{BigInt})
    y = random(DualPerplex{BigInt})
    l = barA(x - y)
    r = barA(x) - barA(y)
    l == r
end

@test begin
    x = random(DualPerplex{BigInt})
    y = random(DualPerplex{BigInt})
    l = barR(x * y)
    r = barR(x) * barR(y)
    l == r
end

@test begin
    x = random(DualPerplex{BigInt})
    y = random(DualPerplex{BigInt})
    l = barR(x + y)
    r = barR(x) + barR(y)
    l == r
end

@test begin
    x = random(DualPerplex{BigInt})
    y = random(DualPerplex{BigInt})
    l = barR(x - y)
    r = barR(x) - barR(y)
    l == r
end

@test begin
    x = random(DualPerplex{Rational{BigInt}})
    y = random(DualPerplex{Rational{BigInt}})
    l = inv(x * y)
    r = inv(y) * inv(x)
    l == r
end

@test begin
    x = random(DualPerplex{Rational{BigInt}})
    y = random(DualPerplex{Rational{BigInt}})
    l = inv(x / y)
    r = y / x
    l == r
end

@test begin
    x = random(DualPerplex{Rational{BigInt}})
    y = random(DualPerplex{Rational{BigInt}})
    l = inv(x \ y)
    r = y \ x
    l == r
end

@test begin
    x = random(DualPerplex{BigInt})
    y = random(DualPerplex{BigInt})
    z = random(DualPerplex{BigInt})
    l = x * (y + z)
    r = (x * y) + (x * z)
    l == r
end

@test begin
    a = Plexifications.randomBigInt()
    x = random(DualPerplex{BigInt})
    y = random(DualPerplex{BigInt})
    l = a * (x + y)
    r = (a * x) + (a * y)
    l == r
end

@test begin
    a = Plexifications.randomBigInt()
    x = random(DualPerplex{BigInt})
    l = a * x
    r = x * a
    l == r
end

@test begin
    x = random(DualPerplex{BigInt})
    y = random(DualPerplex{BigInt})
    z = random(DualPerplex{BigInt})
    l = (x + y) * z
    r = (x * z) + (y * z)
    l == r
end

@test begin
    x = random(DualPerplex{BigInt})
    y = random(DualPerplex{BigInt})
    z = random(DualPerplex{BigInt})
    l = x * (y - z)
    r = (x * y) - (x * z)
    l == r
end

@test begin
    x = random(DualPerplex{BigInt})
    y = random(DualPerplex{BigInt})
    z = random(DualPerplex{BigInt})
    l = (x - y) * z
    r = (x * z) - (y * z)
    l == r
end

@test begin
    x = random(DualPerplex{BigInt})
    n = norm(x)
    n > big(0) || n == big(0)
end

@test begin
    x = random(DualPerplex{BigInt})
    y = random(DualPerplex{BigInt})
    l = norm(x * y)
    r = norm(x) * norm(y)
    l == r
end

@test_throws ErrorException begin
    inv(DualPerplex(0, 0.0))
end

@test_throws ErrorException begin
    random(DualPerplex{Int}) / 0
end

@test_throws ErrorException begin
    0 \ random(DualPerplex{Int})
end

@test begin
    a = Plexifications.randomBigRational()
    z = random(DualPerplex{Rational{BigInt}})
    l = z / a
    r = a \ z
    l == r
end

@test begin
    a = Plexifications.randomBigRational()
    z = random(DualPerplex{Rational{BigInt}})
    l = a / z
    r = z \ a
    l == r
end

@test begin
    z = random(DualPerplex{BigInt})
    z == z + zero(DualPerplex{BigInt})
end

@test begin
    z = random(DualPerplex{BigInt})
    z == z * one(DualPerplex{BigInt})
end

@test begin
    z = random(DualPerplex{BigInt})
    z == z + zero(z)
end

@test begin
    z = random(DualPerplex{BigInt})
    z == z * one(z)
end

@test begin
    x = random(DualPerplex{BigInt})
    y = random(DualPerplex{BigInt})
    l = commutator(x, y)
    r = zero(DualPerplex{BigInt})
    l == r
end

@test begin
    x = random(DualPerplex{BigInt})
    y = random(DualPerplex{BigInt})
    l = alternatorL(x, y)
    r = zero(DualPerplex{BigInt})
    l == r
end

@test begin
    x = random(DualPerplex{BigInt})
    y = random(DualPerplex{BigInt})
    l = alternatorR(x, y)
    r = zero(DualPerplex{BigInt})
    l == r
end

@test begin
    x = random(DualPerplex{BigInt})
    y = random(DualPerplex{BigInt})
    l = flexator(x, y)
    r = zero(DualPerplex{BigInt})
    l == r
end

@test begin
    x = random(DualPerplex{BigInt})
    y = random(DualPerplex{BigInt})
    z = random(DualPerplex{BigInt})
    l = associator(x, y, z)
    r = zero(DualPerplex{BigInt})
    l == r
end

@test begin
    x = random(DualPerplex{BigInt})
    y = random(DualPerplex{BigInt})
    z = random(DualPerplex{BigInt})
    l = jacobiator(x, y, z)
    r = 3 * x * y * z
    l == r
end

@test begin
    l = DualPerplex(Prplex(1,2))
    r = DualPerplex(1,2,0)
    l == r
end

@test begin
    l = DualPerplex(1, 2, 3.0)
    r = DualPerplex(1.0, 2.0, 3.0)
    l == r
end

@test begin
    l = DualPerplex(1, 2, 3, 4.0)
    r = DualPerplex(1.0, 2.0, 3.0, 4.0)
    l == r
end


text:
module Databula

using DungBase
using Dates, VideoIO, Combinatorics, JLSO, UUIDs, Observables, Diana, MAT, DelimitedFiles, SparseArrays, DataStructures
import REPL
using REPL.TerminalMenus

export register_video, register_calibration, add_resfile, test_integrity, test_duration, edit_video#, register_experiment, register_run, register_poi

# const coffeesource = joinpath("/home/yakir/mnt", "coffeesource")
const coffeesource = joinpath(homedir(), "coffeesource")
const sourcefolder = joinpath(coffeesource, "database")
const pixelfolder = joinpath(sourcefolder, "pixel")

if !isdir(pixelfolder)
    @info "creating the source folder" coffeesource
    mkpath(pixelfolder)
else
    @info "found existing source folder" coffeesource
end
f = "videos.jlso"
file = joinpath(sourcefolder, f)
if !isfile(file)
    @info "creating file" f
    JLSO.save(file, AbstractTimeLine[])
else
    @info "found existing file" file
end
f = "calibrations.jlso"
file = joinpath(sourcefolder, f)
if !isfile(file)
    @info "creating file" f
    JLSO.save(file, Calibration[])
else
    @info "found existing file" file
end
f = "pois.jlso"
file = joinpath(sourcefolder, f)
if !isfile(file)
    @info "creating file" f
    JLSO.save(file, Vector{Pair{String, Pair{Symbol, POI{Calibration, Temporal}}}})
else
    @info "found existing file" file
end
f = "runs.jlso"
file = joinpath(sourcefolder, f)
if !isfile(file)
    @info "creating file" f
    JLSO.save(file, Dict{String, Dict{Symbol, POI{Calibration, Temporal}}})
else
    @info "found existing file" file
end

include("terminalmenus.jl")
include("videos.jl")
include("calibrations.jl")
include("intervals.jl")
include("resfile.jl")
include("appinterface.jl")
include("pois.jl")
include("tests.jl")

# include("change2ms.jl")

# include("/home/yakir/dungProject/Databula/src/terminalmenus.jl")
# include("/home/yakir/dungProject/Databula/src/videos.jl")
# include("/home/yakir/dungProject/Databula/src/calibrations.jl")
# include("/home/yakir/dungProject/Databula/src/intervals.jl")

# JSON3.StructType(::Type{AbstractTimeLine}) = JSON3.AbstractType()
# JSON3.subtypekey(::Type{AbstractTimeLine}) = :type
# JSON3.subtypes(::Type{AbstractTimeLine}) = (wholevideo = WholeVideo, fragmentedvideo = FragmentedVideo, disjointvideo = DisjointVideo)

#=for y in (UUIDs, DungBase, Dates), x in names(y, all = true)
T = getproperty(y, x)
if T isa DataType
JSON3.StructType(::Type{T}) = JSON3.Struct()
elseif T isa UnionAll
JSON3.StructType(::Type{<:T}) = JSON3.Struct()
end
end=#

newvideoᵒ = Observable{AbstractTimeLine}(WholeVideo())
const videofile_menu = RadioMenu(["", " "])
empty!(videofile_menu.options)
videofile_menu.pagesize -= 2
const videofiles = VideoFile[]
const videos = AbstractTimeLine[]

h1 = on(newvideoᵒ) do v
    append!(videofile_menu.options, filenames(v))
    if videofile_menu.pagesize < 11
        videofile_menu.pagesize += length(files(v))
    end
    append!(videofiles, files(v))
    push!(videos, v)
end

vs = JLSO.load(joinpath(sourcefolder, "videos.jlso"))["data"]
for v in vs
    newvideoᵒ[] = v
end

function register_video()
    v = newvideo(videofile_menu.options)
    if v ≢ nothing
        newvideoᵒ[] = v
        JLSO.save(joinpath(sourcefolder, "videos.jlso"), videos)
    end
end

newcalibrationᵒ = Observable{Calibration}(Calibration())
const calibration_menu = MultiSelectMenu(["", " "])
empty!(calibration_menu.options)
calibration_menu.pagesize -= 2
const calibrations = Calibration[]

# ns2s(t::Nanosecond) = t/Nanosecond(1000000000)

_formatcalibration(x) = string("file: ", filenames(x), "; start: ", Time(0) + start(x))

h2 = on(newcalibrationᵒ) do c
    push!(calibration_menu.options, _formatcalibration(c))
    if isempty(calibration_menu.selected)
        push!(calibration_menu.selected, 1)
    end
    if calibration_menu.pagesize < 11
        calibration_menu.pagesize += 1
    end
    push!(calibrations, c)
end

cs = JLSO.load(joinpath(sourcefolder, "calibrations.jlso"))["data"]
for c in cs
    newcalibrationᵒ[] = c
end

#=function getnewkey(di)
local k
while haskey(di, (k=uuid1();)) end
k
end=#

function register_calibration()
    #=file = joinpath(sourcefolder, "calibrations.json")
    calibrations = open(file, "r") do i
    JSON3.read(i, Dict{String, Calibration})
    end=#
    boards = unique(c.board for c in values(calibrations))
    # k = getnewkey(calibrations)
    newcalibrationᵒ[] = newcalibration(boards)
    JLSO.save(joinpath(sourcefolder, "calibrations.jlso"), calibrations)
end

newpoiᵒ = Observable{Any}(" " => :k => POI{Calibration, Temporal}())
# newpoiᵒ = Observable{Pair{String, Pair{Symbol, POI{Calibration, Temporal}}}}(" " => :k => POI{Calibration, Temporal}())
h3 = on(newpoiᵒ) do file_type_poi
    pois = JLSO.load(joinpath(sourcefolder, "pois.jlso"))
    f, kp = file_type_poi
    pois[f] = kp
    JLSO.save(joinpath(sourcefolder, "pois.jlso"), pois)
end

function add_resfile()
    ids = newresfile()
    r = register_pois(ids)
end

# const COMMON_FACTORS = (:run, :experimenter, :date, :comment, :person)
# const FACTORS = (:displaced,:place,:pellet_manipulation,:id,:displace_location,:nest2feeder,:displace_direction,:azimuth,:species,:transfer,:nest_coverage)

# include("calibrations.jl")
# include("intervals.jl")
# include("experiments.jl")
# include("runs.jl")
# include("pois.jl")

end # module


text:
# Parser
# ======

mutable struct StreamReader
    tokenizer::Tokenizer
    # mutable state
    stack::Vector{Symbol}
    queue::Vector{Token}
end

function StreamReader(input::IO)
    return StreamReader(Tokenizer(input), Symbol[], Token[])
end

function linenumber(reader::StreamReader)
    linenum = reader.tokenizer.linenum
    for token in reader.queue
        if token.kind == :newline
            linenum -= 1
        end
    end
    return linenum
end

struct ParseError <: Exception
    msg::String
end

parse_error(msg, linenum) = throw(ParseError("$(msg) at line $(linenum)"))
unexpectedtoken(token, linenum) = parse_error("unexpected $(tokendesc(token))", token.kind == :newline ? linenum - 1 : linenum)

# Basically, this function parses line by line.
function parsetoken(reader::StreamReader)
    if !isempty(reader.queue)
        return popfirst!(reader.queue)
    end
    top() = isempty(reader.stack) ? :__empty__ : reader.stack[end]
    emit(token) = push!(reader.queue, token)
    tokenizer = reader.tokenizer
    stack = reader.stack
    token = readtoken(tokenizer, rhs=top() == :array)
    if token.kind == :whitespace
        return token
    elseif top() == :array
        if token.kind ∈ (:newline, :comment)
            return token
        elseif isatomicvalue(token)
            let token = peektoken(tokenizer, rhs=true)
                if token.kind == :whitespace
                    emit(readtoken(tokenizer, rhs=true))
                    token = peektoken(tokenizer, rhs=true)
                end
                if token.kind == :comma
                    emit(readtoken(tokenizer, rhs=true))
                elseif token.kind ∈ (:single_bracket_right, :newline, :comment)
                    # ok
                else
                    unexpectedtoken(token, tokenizer.linenum)
                end
            end
            return token
        elseif token.kind == :single_bracket_left
            # found a new array inside the current array (nested arrays)
            push!(stack, :array)
            emit(token)
            return TOKEN_INLINE_ARRAY_BEGIN
        elseif token.kind == :curly_brace_left
            # found an inline table inside the current array
            push!(stack, :table)
            emit(token)
            return TOKEN_INLINE_TABLE_BEGIN
        elseif token.kind == :single_bracket_right
            pop!(stack)
            emit(TOKEN_INLINE_ARRAY_END)
            if top() == :array
                let token = peektoken(tokenizer, rhs=true)
                    if token.kind == :whitespace
                        emit(readtoken(tokenizer, rhs=true))
                        token = peektoken(tokenizer, rhs=true)
                    end
                    if token.kind == :comma
                        emit(readtoken(tokenizer, rhs=true))
                    elseif token.kind ∈ (:single_bracket_right, :newline, :comment)
                        # ok
                    else
                        unexpectedtoken(token, tokenizer.linenum)
                    end
                end
            elseif top() == :table
                let token = peektoken(tokenizer)
                    if token.kind == :whitespace
                        emit(readtoken(tokenizer))
                        token = peektoken(tokenizer)
                    end
                    if token.kind == :comma
                        emit(readtoken(tokenizer))
                    elseif token.kind == :curly_brace_right
                        # ok
                    else
                        unexpectedtoken(token, tokenizer.linenum)
                    end
                end
            else @assert isempty(stack)
                parselineend(reader)
            end
            return token
        else
            unexpectedtoken(token, tokenizer.linenum)
        end
    elseif top() == :table
        if iskey(token)
            let token = readtoken(tokenizer)
                # read equal sign
                if token.kind == :whitespace
                    emit(token)
                    token = readtoken(tokenizer)
                end
                if token.kind != :equal
                    unexpectedtoken(token, tokenizer.linenum)
                end
                emit(token)
                # read value
                token = readtoken(tokenizer, rhs=true)
                if token.kind == :whitespace
                    emit(token)
                    token = readtoken(tokenizer, rhs=true)
                end
                if isatomicvalue(token)
                    emit(token)
                    # read comma or right curly brace
                    token = peektoken(tokenizer)
                    if token.kind == :whitespace
                        emit(readtoken(tokenizer))
                        token = peektoken(tokenizer)
                    end
                    if token.kind == :comma
                        emit(readtoken(tokenizer))
                    elseif token.kind == :curly_brace_right
                        # ok
                    else
                        unexpectedtoken(token, tokenizer.linenum)
                    end
                elseif token.kind == :single_bracket_left
                    push!(stack, :array)
                    emit(TOKEN_INLINE_ARRAY_BEGIN)
                    emit(token)
                elseif token.kind == :curly_brace_left
                    push!(stack, :table)
                    emit(TOKEN_INLINE_TABLE_BEGIN)
                    emit(token)
                else
                    unexpectedtoken(token, tokenizer.linenum)
                end
            end
            return token
        elseif token.kind == :curly_brace_right
            pop!(stack)
            emit(TOKEN_INLINE_TABLE_END)
            if top() == :array
                let token = peektoken(tokenizer, rhs=true)
                    if token.kind == :whitespace
                        emit(readtoken(tokenizer, rhs=true))
                        token = peektoken(tokenizer, rhs=true)
                    end
                    if token.kind == :comma
                        emit(readtoken(tokenizer, rhs=true))
                    elseif token.kind ∈ (:single_bracket_right, :newline, :comment)
                        # ok
                    else
                        unexpectedtoken(token, tokenizer.linenum)
                    end
                end
            elseif top() == :table
                let token = peektoken(tokenizer)
                    if token.kind == :whitespace
                        emit(readtoken(tokenizer))
                        token = peektoken(tokenizer)
                    end
                    if token.kind == :comma
                        emit(readtoken(tokenizer))
                    elseif token.kind == :curly_brace_right
                        # ok
                    else
                        unexpectedtoken(token, tokenizer.linenum)
                    end
                end
            else @assert isempty(stack)
                parselineend(reader)
            end
            return token
        else
            unexpectedtoken(token, tokenizer.linenum)
        end
    elseif iskey(token)
        let token = readtoken(tokenizer)
            if token.kind == :whitespace
                emit(token)
                token = readtoken(tokenizer)
            end
            if token.kind != :equal
                unexpectedtoken(token, tokenizer.linenum)
            end
            emit(token)
            token = readtoken(tokenizer, rhs=true)
            if token.kind == :whitespace
                emit(token)
                token = readtoken(tokenizer, rhs=true)
            end
            if isatomicvalue(token)
                emit(token)
                parselineend(reader)
            elseif token.kind == :single_bracket_left  # <key> = [
                push!(stack, :array)
                emit(TOKEN_INLINE_ARRAY_BEGIN)
                emit(token)
            elseif token.kind == :curly_brace_left  # <key> = {
                push!(stack, :table)
                emit(TOKEN_INLINE_TABLE_BEGIN)
                emit(token)
            else
                unexpectedtoken(token, tokenizer.linenum)
            end
        end
        return token
    elseif token.kind ∈ (:single_bracket_left, :double_brackets_left)
        endkind = token.kind == :single_bracket_left ? :single_bracket_right : :double_brackets_right
        emit(token)
        token = readtoken(tokenizer)
        if token.kind == :whitespace
            emit(token)
            token = readtoken(tokenizer)
        end
        if !iskey(token)
            unexpectedtoken(token, tokenizer.linenum)
        end
        emit(token)
        while (token = readtoken(tokenizer)).kind != endkind
            if token.kind == :whitespace
                emit(token)
                token = readtoken(tokenizer)
            end
            if token.kind == endkind
                break
            elseif token.kind != :dot
                unexpectedtoken(token, tokenizer.linenum)
            end
            emit(token)
            token = readtoken(tokenizer)
            if token.kind == :whitespace
                emit(token)
                token = readtoken(tokenizer)
            end
            if !iskey(token)
                unexpectedtoken(token, tokenizer.linenum)
            end
            emit(token)
        end
        emit(token)
        emit(endkind == :single_bracket_right ? TOKEN_TABLE_END : TOKEN_ARRAY_END)
        parselineend(reader)
        return endkind == :single_bracket_right ? TOKEN_TABLE_BEGIN : TOKEN_ARRAY_BEGIN
    elseif token.kind ∈ (:newline, :comment, :eof)
        return token
    else
        unexpectedtoken(token, tokenizer.linenum)
    end
end

function parselineend(reader::StreamReader)
    emit(token) = push!(reader.queue, token)
    tokenizer = reader.tokenizer
    token = readtoken(tokenizer)
    if token.kind == :whitespace
        emit(token)
        token = readtoken(tokenizer)
        if token.kind == :comment
            emit(token)
            token = readtoken(tokenizer)
        end
    elseif token.kind == :comment
        emit(token)
        token = readtoken(tokenizer)
    end
    if !iseol(token)
        unexpectedtoken(token, tokenizer.linenum)
    end
    emit(token)
    return token
end

# Infer the type of elements in the current array.
function arraytype(reader::StreamReader)
    t = Any
    token = parsetoken(reader)
    tokens = [token]
    while token.kind != :inline_array_end && t == Any
        if isatomicvalue(token)
            if issigned(token)
                t = Int
            elseif isunsigned(token)
                t = UInt
            elseif isfloat(token)
                t = Float64
            elseif isstring(token)
                t = String
            elseif isboolean(token)
                t = Bool
            else
                @assert(false)
            end
        elseif token.kind == :inline_array_begin
            t = Vector
        elseif token.kind == :inline_table_begin
            t = Dict{String}
        else
            # ignore
        end
        token = parsetoken(reader)
        push!(tokens, token)
    end
    prepend!(reader.queue, tokens)
    return t
end

parse(str::AbstractString) = parse(IOBuffer(str))
parsefile(filename::AbstractString) = open(parse, filename)

const Table = Dict{String,Any}

function parse(input::IO)
    reader = StreamReader(input)
    dupdef() = parse_error("found a duplicated definition", linenumber(reader))
    root = Table()
    tablekeys = Set{Vector{String}}()
    key = nothing
    node = root
    stack = []
    while (token = parsetoken(reader)).kind != :eof
        if iskey(token)
            key = keyname(token)
        elseif isatomicvalue(token)
            if node isa Array
                x = value(token)
                if !isempty(node) && typeof(node[1]) != typeof(x)
                    # This line number may be wrong.
                    parse_error("mixed array types", linenumber(reader))
                end
                push!(node, x)
            else
                haskey(node, key) && dupdef()
                node[key] = value(token)
            end
        elseif token.kind == :inline_array_begin
            push!(stack, node)
            if node isa Array
                push!(node, arraytype(reader)[])
                node = node[end]
            else
                haskey(node, key) && dupdef()
                node[key] = arraytype(reader)[]
                node = node[key]
            end
        elseif token.kind == :inline_array_end
            node = pop!(stack)
        elseif token.kind == :inline_table_begin
            push!(stack, node)
            if node isa Array
                push!(node, Table())
                node = node[end]
            else
                haskey(node, key) && dupdef()
                node = get!(node, key, Table())
            end
        elseif token.kind == :inline_table_end
            node = pop!(stack)
        elseif token.kind == :table_begin  # [foo.bar]
            keys = String[]
            while (token = parsetoken(reader)).kind != :table_end
                iskey(token) && push!(keys, keyname(token))
            end
            node = root
            for i in 1:length(keys)
                node = get!(node, keys[i], Table())
                if node isa Table
                    # ok
                elseif node isa Array
                    node = node[end]
                else
                    dupdef()
                end
            end
            keys ∈ tablekeys && dupdef()
            push!(tablekeys, keys)
            key = nothing
        elseif token.kind == :array_begin  # [[foo.bar]]
            keys = String[]
            while (token = parsetoken(reader)).kind != :array_end
                iskey(token) && push!(keys, keyname(token))
            end
            node = root
            for i in 1:length(keys)-1
                node = get!(node, keys[i], Table())
                if node isa Table
                    # ok
                elseif node isa Array
                    node = node[end]
                else
                    dupdef()
                end
            end
            node = push!(get!(node, keys[end], Table[]), Table())[end]
            key = nothing
        else
            # ignore
        end
    end
    return root
end


text:
using PCCA
using Test

@testset "PCCA.jl" begin
    # Write your tests here.
	x = rand(10, 10)
	x = x ./ sum(x, dims=2)
	c = pcca(x, 2)

	function randomstochasticmatrix(n, reversible=true)
		P = rand(n,n)
		if reversible
			P = (P + P') / 2
		end
		P ./= sum(P, dims=2)
	end

	@testset "Reversible = $rev" for rev in [true, false]
		P = [randomstochasticmatrix(3+mod(i,12), rev) for i in 1:10]
		@testset "Method $method" for method in [PCCA.BaseSolver, PCCA.ArnoldiSolver, PCCA.KrylovSolver]
			for P in P, n in 2:size(P,1)-1
				χ = pcca(P, n, solver=method(), optimize=true)
				a = PCCA.crispassignments(χ)

				#@show χ
				@test all(sum(χ, dims=2) .≈ 1)
			end
		end
	end

end


text:
using Multilane
using StatsBase
using MCTS
using JLD
# using Plots
# using StatPlots

solvers = Dict{String, Any}(
    "normal" => BehaviorSolver(Multilane.NORMAL, false, MersenneTwister(123)),
    "lane_seeking" => IDMLaneSeekingSolver(Multilane.NORMAL, MersenneTwister(123))
)

test = TestSet(lambda=10., solver_key="lane_seeking", behaviors="agents", N=500)

objects = gen_initials([test])

@show objects["param_table"]
objects["solvers"] = solvers

metrics=[MaxBrakeMetric(), NumBehaviorBrakesMetric(3,"aggressive",2.5)]
objects["metrics"] = metrics
# results = evaluate([test], objects, metrics=metrics, max_steps=200)
results_file_list = sbatch_spawn([test], objects,
                                 submit_command="sbatch",
                                 template_name="sherlock.sh",
                                 time_per_batch="10:00")

#=

results = gather_results(results_file_list)
stats = results["stats"]
@show names(results["stats"])

@show sum(results["stats"][:aggressive_nb_brakes])/sum(results["stats"][:nb_brakes])

# unicodeplots()
histogram(results["stats"][:max_brake], title="max brake histogram")
gui()

i = 20
problem, sim, policy = rerun(results, i);
f = write_tmp_gif(problem, sim)
run(`gifview $f`)

=#

# filename = string("results_", Dates.format(Dates.now(),"u_d_HH_MM"), ".jld")
# results["histories"] = nothing
# save(filename, results)
# println("results saved to $filename")


text:
import DelayEmbeddings: genembed, AbstractDataset, Dataset
export EmbeddingTE

function rc(x::Union{AbstractDataset, Vector{T}, Vector{Vector{T}}},
        dim::Union{Int, AbstractVector{Int}},
        τ::Union{Int, AbstractVector{Int}}, forward = false) where T <: Number

    if x isa Vector{T} where T <: Number
        # Fixed dimension, fixed lag
        if dim isa Int && τ isa Int
            pos = [1 for x in dim:-1:1]
            if τ > 0
                lags = [d*τ for d in dim:-1:1]
            elseif τ <= 0 && forward
                lags = [d*τ for d in dim:-1:1]
            elseif τ <= 0 && !forward
                lags = [d*τ for d in 0:(dim-1)]
            end
        end

        # Fixed dimension, multiple lags (number of lags must match dimension)
        if dim isa Int && length(τ) > 1
            length(τ) == dim || throw(ArgumentError("length(τ) must equal dim if multiple lags are specified manually (got length(τ) = $(length(τ)), dim=$(dim))"))
            length(unique(τ)) == length(τ) || throw(ArgumentError("Tried to specify reconstruction lags manually, but there are repeated lags. Provide unique lags."))
            pos = pos = [1 for x in dim:-1:1]
            return pos, τ
        end
    end

    # Multiple time series input
    if (x isa Vector{Vector{T}} where T <: Number) || (x isa AbstractDataset)
        if x isa AbstractDataset
            N = length(x)
        elseif x isa Vector{Vector{T}} where T <: Number
            N = size(x, 1)
        end

        if dim isa Int
            dim % N == 0 || throw(ArgumentError("If using multiple (`N` different) time series in a marginal, each time series is lagged `dim/N` times. Hence, `dim` must be a multiple of `N`."))

            lD = Int(dim / N) # "local" reconstruction dimension for this time series

            pos = Vector{Int}(undef, 0)

            for i = 1:N
                append!(pos, repeat([i], lD))
            end

            if τ isa Int
                if τ > 0
                    llags = [d*τ for d in lD:-1:1]
                elseif τ <= 0 && forward
                    llags = [d*τ for d in lD:-1:1]
                elseif τ <= 0 && !forward
                    llags = [d*τ for d in 0:(lD - 1)]
                end
                lags = repeat(llags, Int(dim/lD))
            elseif τ isa AbstractVector{Int}
                length(τ) == N || throw(ArgumentError("Tried using `N = $N` different time series in a marginal, but $(length(τ)) reconstruction lags were provided. The number of lags must equal `N`"))
                lags = Vector{Int}(undef, 0)
                for i = 1:length(x)
                    if τ[i] > 0
                        llags = [d*τ[i] for d in lD:-1:1]
                    elseif τ[i] <= 0 && forward
                        llags = [d*τ[i] for d in 0:(lD - 1)]
                    elseif τ[i] <= 0 && !forward
                        llags = [d*τ[i] for d in 0:(lD - 1)]
                    end
                    append!(lags, llags)
                end
            end
        end

        if dim isa AbstractVector{Int}
            length(dim) == N || throw(ArgumentError("There must be precisely one dimension specification per time series. Got $(length(dim)) specifications for $N time series."))
            pos = Vector{Int}(undef, 0)
            lags = Vector{Int}(undef, 0)

            for (i, lD) in enumerate(dim)
                append!(pos, repeat([i], lD))

                if τ isa Int
                    if τ > 0
                        llags = [d*τ for d in lD:-1:1]
                    elseif τ <= 0 && !forward
                        llags = [d*τ for d in lD:-1:1]
                    elseif τ <= 0 && forward
                        llags = [d*τ for d in 0:(lD - 1)]
                    end
                    append!(lags, llags)
                elseif τ isa AbstractVector{Int}
                    length(τ) == N || throw(ArgumentError("Tried using `N = $N` different time series in a marginal, but $(length(τ)) reconstruction lags were provided. The number of lags must equal `N`"))
                    if τ[i] > 0
                        llags = [d*τ[i] for d in lD:-1:1]
                    elseif τ[i] <= 0 && forward
                        llags = [d*τ[i] for d in lD:-1:1]
                    elseif τ[i] <= 0 && !forward
                        llags = [d*τ[i] for d in 0:(lD - 1)]
                    end
                    append!(lags, llags)
                end
            end
        end
    end
    return pos, lags
end

"""
    EmbeddingTE(; dS = 1, dT = 1, d𝒯 = 1, dC = 1, τS = -1, τT = -1, η𝒯 = 1, τC = -1)

Embedding parameters for transfer entropy analysis.

## Convention for generalized delay reconstruction

This struct contains instructions for transfer entropy computations using the following convention.
Let ``x(t)`` be time series for the source variable, ``y(t)`` be the time series for the target variable and
``z(t)`` the time series for any conditional variable. To compute transfer entropy, we need the
following marginals:

```math
\\begin{aligned}
\\mathcal{T}^{(d_{\\mathcal{T}})} &= \\{(y(t+\\eta^{d_{\\mathcal{T}}}), \\ldots, y(t+\\eta^2), y(t+\\eta^1) \\} \\\\
T^{(d_{T})} &= \\{ (y(t+\\tau^0_{T}), y(t+\\tau^1_{T}), y(t+\\tau^2_{T}), \\ldots, y(t + \\tau^{d_{T} - 1}_{T})) \\} \\\\
S^{(d_{S})} &= \\{ (x(t+\\tau^0_{S}), x(t+\\tau^1_{S}), x(t+\\tau^2_{S}), \\ldots, x(t + \\tau^{d_{S} - 1}_{S})) \\} \\\\
C^{(d_{C})} &= \\{ (z(t+\\tau^0_{C}), z(t+\\tau^1_{C}), z(t+\\tau^2_{C}), \\ldots, z(t + \\tau^{d_{C} - 1}_{C})) \\}
\\end{aligned}
```

Depending on the application, the delay reconstruction lags ``\\tau^k_{T} \\leq 0``, ``\\tau^k_{S} \\leq 0``, and ``\\tau^k_{C} \\leq 0``
may be equally spaced, or non-equally spaced. The predictions lags ``\\eta^k``may also be equally spaced
or non-equally spaced, but are always positive. For transfer entropy, convention dictates that at least one
``\\tau^k_{T}``, one ``\\tau^k_{S}`` and one ``\\tau^k_{C}`` equals zero. This way, the ``T``, ``S`` and ``C`` marginals
always contains present/past states,
while the ``\\mathcal T`` marginal contain future states relative to the other marginals.

Combined, we get the generalized delay reconstruction ``\\mathbb{E} = (\\mathcal{T}^{(d_{\\mathcal{T}})}, T^{(d_{T})}, S^{(d_{S})}, C^{(d_{C})})``. Transfer entropy is then computed as

```math
\\begin{aligned}
TE_{S \\rightarrow T | C} = \\int_{\\mathbb{E}} P(\\mathcal{T}, T, S, C) \\log_{b}{\\left(\\frac{P(\\mathcal{T} | T, S, C)}{P(\\mathcal{T} | T, C)}\\right)},
\\end{aligned}
```

or, if conditionals are not relevant,

```math
\\begin{aligned}
TE_{S \\rightarrow T} = \\int_{\\mathbb{E}} P(\\mathcal{T}, T, S) \\log_{b}{\\left(\\frac{P(\\mathcal{T} | T, S)}{P(\\mathcal{T} | T)}\\right)},
\\end{aligned}
```

Here,

- ``\\mathcal{T}`` denotes the ``d_{\\mathcal{T}}``-dimensional set of vectors furnishing the future states of ``T``,
- ``T`` denotes the ``d_{T}``-dimensional set of vectors furnishing the past and present states of ``T``,
- ``S`` denotes the ``d_{S}``-dimensional set of vectors furnishing the past and present of ``S``, and
- ``C`` denotes the ``d_{C}``-dimensional set of vectors furnishing the past and present of ``C``.

## Keyword arguments

### Specifying dimensions for generalized delay reconstructions of marginals

`dS`, `dT`, `d𝒯`, and `dC` are the dimensions of the ``S``, ``T``, ``\\mathcal{T}``,
and ``C`` marginals. The dimensions of each marginal can be specified manually by setting
either `dS`, `dT`, `d𝒯`, or `dC` to a *positive* integer number. Alternatively, the dimension
of each marginal can be optimised by setting either `dS`, `dT`, `d𝒯`, or `dC` to an
instance of [`OptimiseDim`](@ref)
(e.g. `EmbeddingTE(dT = OptimDim(method_delay = "ac_zero", method_dim = "f1nn")`).

### Specifying delays for generalized delay reconstructions of marginals

The corresponding embedding delay lags are given by `τS`, `τT` and `τC`. The delays
for each marginal can be specified manually by setting either `dS`, `dT`, `d𝒯`, or `dC`
to a *negative* integer number. The delay defaults for each marginal is -1 (but is set to zero
if the marginal is one-dimensional), and must always be negative. Alternatively, delays can
be estimated numerically by setting either `dS`, `dT`, `d𝒯`, and `dC`
to an instance of [`OptimiseDelay`](@ref) (e.g. `dS = OptimiseDelay(method_delay = "ac_zero")`).

The prediction lag `η` can be either positive or negative, but should not be zero.

In summary, one can provide

- A single delay ``\\tau``, in which case ``\\tau_{T} = \\{0, \\tau, 2\\tau, \\ldots, (d_{T}- 1)\\tau \\}``, or
- All the delays manually. If so, then the number of delays must match the dimension of the marginal).

For the prediction lag, one can provide

- A single delay ``\\eta_f``, in which case ``\\eta_{\\mathcal{T}} = \\{\\eta_f, 2\\eta_f, \\ldots, (d_{\\mathcal{T}} - 1)\\eta_f \\}``, or
- All the delays manually. If so, then the number of delays must equal ``d_{\\mathcal{T}}``, which is the dimension of the marginal).

!!! note
    If both the delay and the dimension for a given marginal is to be estimated numerically, make sure
    to use the same delay estimation method for both
    the [`OptimiseDelay`](@ref) and  [`OptimiseDim`](@ref) instances.

## Examples

Default parameters:

```jldoctest
using TransferEntropy
p = EmbeddingTE()

# output
EmbeddingTE(dS=1, dT=1, dC=1, d𝒯=1, τS=-1, τT=-1, τC=-1, η𝒯=1)
```
"""
@Base.kwdef struct EmbeddingTE
    dS::Union{Int, AbstractVector{Int}} = 1
    dT::Union{Int, AbstractVector{Int}} = 1
    d𝒯::Union{Int, AbstractVector{Int}} = 1
    dC::Union{Int, AbstractVector{Int}, Nothing} = 1
    τS::Union{Int, AbstractVector{Int}} = -1
    τT::Union{Int, AbstractVector{Int}} = -1
    η𝒯::Union{Int, AbstractVector{Int}} = 1
    τC::Union{Int, AbstractVector{Int}, Nothing} = -1

    function EmbeddingTE(
            dS::Union{Int, AbstractVector{Int}},
            dT::Union{Int, AbstractVector{Int}},
            d𝒯::Union{Int, AbstractVector{Int}},
            dC::Union{Int, AbstractVector{Int}},
            τS::Union{Int, AbstractVector{Int}},
            τT::Union{Int, AbstractVector{Int}},
            η𝒯::Union{Int, AbstractVector{Int}},
            τC::Union{Int, AbstractVector{Int}})

        if dS isa Int
            dS > 0 || throw(ArgumentError("dimension for marginal S must be a positive integer (got dS=$(dS))"))
        end
        if dT isa Int
            dT > 0 || throw(ArgumentError("dimension for marginal T must be a positive integer (got dT=$(dT))"))
        end
        if dC isa Int
            dC > 0 || throw(ArgumentError("dimension for marginal C must be a positive integer (got dC=$(dC))"))
        end
        if d𝒯 isa Int
            d𝒯 > 0 || throw(ArgumentError("dimension for marginal 𝒯 must be a positive integer (got d𝒯=$(d𝒯))"))
        end
        if τS isa Int
            τS < 0 || throw(ArgumentError("delay for marginal S must be a negative integer (got τS=$(τS))"))
        end
        if τT isa Int
            τT < 0 || throw(ArgumentError("delay for marginal T must be a negative integer (got τT=$(τT))"))
        end
        if τC isa Int
            τC < 0 || throw(ArgumentError("delay for marginal C must be a negative integer (got τC=$(τC))"))
        end

        if τS isa AbstractVector{Int} || τS isa AbstractUnitRange{Int64}
            all(τS .<= 0) || throw(ArgumentError("delays for marginal S must be <= 0 (got τS=$(τS))"))
        end

        if τT isa AbstractVector{Int} || τT isa AbstractUnitRange{Int64}
            all(τT .<= 0) || throw(ArgumentError("delays for marginal T must be <= 0 (got τT=$(τT))"))
        end

        if τC isa AbstractVector{Int} || τC isa AbstractUnitRange{Int64}
            all(τC .<= 0) || throw(ArgumentError("delays for marginal C must be <= 0 (got τC=$(τC))"))
        end

        new(dS, dT, d𝒯, dC, τS, τT, η𝒯, τC)
    end

end

function Base.show(io::IO, x::EmbeddingTE)
    s = "EmbeddingTE(dS=$(x.dS), dT=$(x.dT), dC=$(x.dC), d𝒯=$(x.d𝒯), τS=$(x.τS), τT=$(x.τT), τC=$(x.τC), η𝒯=$(x.η𝒯))"
    print(io, s)
end

function get_delay_reconstruction_params(source, target, p::EmbeddingTE)
    pos_𝒯, lags_𝒯 = rc(target, p.d𝒯, p.η𝒯, true)
    pos_T, lags_T = rc(target, p.dT, p.τT, false)
    pos_S, lags_S = rc(source, p.dS, p.τS, false)
    pos_C, lags_C = rc(source, p.dC, p.τC, false)

    js = ([pos_𝒯; pos_T; pos_S]...,)
    τs = ([lags_𝒯; lags_T; lags_S]...,)

    return τs, js
end

function get_delay_reconstruction_params(source, target, cond, p::EmbeddingTE)
    pos_𝒯, lags_𝒯 = rc(target, p.d𝒯, p.η𝒯, true)
    pos_T, lags_T = rc(target, p.dT, p.τT, false)
    pos_S, lags_S = rc(source, p.dS, p.τS, false)
    pos_C, lags_C = rc(cond, p.dC, p.τC, false)

    js = ([pos_𝒯; pos_T; pos_S; pos_C]...,)
    τs = ([lags_𝒯; lags_T; lags_S; pos_C]...,)

    return τs, js
end

"""
    te_embed(source::AbstractVector{T}, target::AbstractVector{T}, p::EmbeddingTE) → (points, vars, τs)
    te_embed(source::AbstractVector{T}, target::AbstractVector{T}, cond::AbstractVector{T}, p::EmbeddingTE) → (points, vars, τs)

Generalised delay reconstruction of `source` and `target` (and `cond` if provided)
for transfer entropy computation using embedding parameters provided by the [`EmbeddingTE`](@ref)
instance `p`.

Returns a tuple of the embedded `points`, `vars` (a [`TEVars`](@ref) instance that keeps track of which
variables of the embedding belong to which marginals of the reconstruction; indices are: source = 1,
target = 2, cond = 3), and a tuple `τs`, which stores the lags for each variable of the reconstruction.
"""
function te_embed(source::AbstractVector{T}, target::AbstractVector{T}, p::EmbeddingTE) where T

    #@show p.τS
    #if (p.τS isa Int && p.τS > 0) || (length(p.τS) > 1 && any(p.τS[p.τS .> 0]))
    #    @warn("Backwards lag τS should be negative. You might be getting nonsensical results!")
    #end

    # Get lags and posisions separately for each marginal
    pos_𝒯, lags_𝒯 = rc(target, p.d𝒯, p.η𝒯, true)
    pos_T, lags_T = rc(target, p.dT, p.τT, false)
    pos_S, lags_S = rc(source, p.dS, p.τS, false)

    # Add one to the index of the positions for the target (rc doesn't know it is in fact our second time series)
    # TODO: make sure this works when `source` and `target` are multiple time series
    pos_𝒯 .= pos_𝒯 .+ 1
    pos_T .= pos_T .+ 1

    js = ([pos_𝒯; pos_T; pos_S]...,)
    τs = ([lags_𝒯; lags_T; lags_S]...,)

    # TODO: This only works for single time series at the moment
    ts = Dataset(source, target)

    # The reconstructed points
    pts = genembed(ts, τs, js)
    d𝒯 = length(pos_𝒯)
    dT = length(pos_T)
    dS = length(pos_S)

    # Which columns/variables map to which marginals?
    vars = TEVars(
        𝒯  = 1:(d𝒯)           |> collect,
        T = 1+(d𝒯):dT+(d𝒯)     |> collect,
        S = 1+(dT+d𝒯):dS+(d𝒯+dT) |> collect)

    return pts, vars, τs, js
end

function te_embed(source::AbstractVector{T}, target::AbstractVector{T}, cond::AbstractVector{T}, p::EmbeddingTE) where T

    #@show p.τS
    #if (p.τS isa Int && p.τS > 0) || (length(p.τS) > 1 && any(p.τS[p.τS .> 0]))
    #    @warn("Backwards lag τS should be negative. You might be getting nonsensical results!")
    #end
    # Get lags and posisions separately for each marginal
    pos_𝒯, lags_𝒯 = rc(target, p.d𝒯, p.η𝒯,  true)
    pos_T, lags_T = rc(target, p.dT, p.τT, false)
    pos_S, lags_S = rc(source, p.dS, p.τS, false)
    pos_C, lags_C = rc(cond,   p.dC, p.τC, false)

    # Add one to the index of the positions for the target (rc doesn't know it is in fact our second time series)
    # TODO: make sure this works when `source` and `target` are multiple time series
    pos_𝒯 .= pos_𝒯 .+ 1
    pos_T .= pos_T .+ 1
    pos_C .= pos_C .+ 2

    js = ([pos_𝒯; pos_T; pos_S; pos_C]...,)
    τs = ([lags_𝒯; lags_T; lags_S; lags_C]...,)

    # TODO: This only works for single time series at the moment
    ts = Dataset(source, target, cond)

    # The reconstructed points
    pts = genembed(ts, τs, js)
    d𝒯 = length(pos_𝒯)
    dT = length(pos_T)
    dS = length(pos_S)
    dC = length(pos_C)

    # Which columns/variables map to which marginals?
    vars = TEVars(
        𝒯 = 1:(d𝒯)               |> collect,
        T = 1+(d𝒯):dT+(d𝒯)         |> collect,
        S = 1+(dT+d𝒯):dS+(d𝒯+dT)     |> collect,
        C = 1+(dT+d𝒯+dS):dC+(d𝒯+dT+dS) |> collect)

    return pts, vars, τs, js
end

"""
    TEVars(𝒯::Vector{Int}, T::Vector{Int}, S::Vector{Int})
    TEVars(𝒯::Vector{Int}, T::Vector{Int}, S::Vector{Int}, C::Vector{Int})
    TEVars(;𝒯 = Int[], T = Int[], S = Int[], C = Int[]) -> TEVars

Which axes of the state space correspond to the future of the target (`𝒯`),
the present/past of the target (`T`), the present/past of the source (`S`), and
the present/past of any conditioned variables (`C`)?  This information is used by
the transfer entropy estimators to ensure that marginal distributions are computed correctly.

Indices correspond to variables of the embedding, or, equivalently, colums of a `Dataset`.

- The three-argument constructor assumes there will be no conditional variables.
- The four-argument constructor assumes there will be conditional variables.

"""
struct TEVars
    𝒯::Vector{Int}
    T::Vector{Int}
    S::Vector{Int}
    C::Vector{Int}
end

"""
    TEVars(𝒯::Vector{Int},
            T::Vector{Int},
            S::Vector{Int})

Which axes of the state space correspond to the future of the target,
the present/past of the target, and the present/past of the source?
Indices correspond to variables of the embedding or colums of a `Dataset`.

This information is used by the transfer entropy estimators to ensure
the marginal distributions are computed correctly.

This three-argument constructor assumes there will be no conditional variables.
"""
TEVars(𝒯::Vector{Int}, T::Vector{Int}, S::Vector{Int}) = TEVars(𝒯, T, S, Int[])

TEVars(;𝒯::Vector{Int} = Int[],
	    T::Vector{Int} = Int[],
	    S::Vector{Int} = Int[],
		C::Vector{Int} = Int[]) =
	TEVars(𝒯, T, S, C)

function Base.show(io::IO, tv::TEVars)
    s = "$(typeof(tv))(𝒯 = $(tv.𝒯), T = $(tv.T), S = $(tv.S), C = $(tv.C))"
    print(io, s)
end


text:
module Vcov

using Combinatorics
using GroupedArrays
using LinearAlgebra
using StatsBase
using Tables

import Base: @propagate_inbounds, size, length, getindex, show

##############################################################################
##
## Mimimum RegressionModel used in Vcov
##
##############################################################################
struct VcovData{T, N} <: RegressionModel
    modelmatrix::Matrix{Float64} # X
    crossmodelmatrix::T          # X'X in the simplest case. Can be Matrix but preferably Factorization
    residuals::Array{Float64, N} # vector or matrix of residuals (matrix in the case of IV, residuals of Xendo on (Z, Xexo))
    dof_residual::Int
end
StatsBase.modelmatrix(x::VcovData) = x.modelmatrix
StatsBase.crossmodelmatrix(x::VcovData) = x.crossmodelmatrix
StatsBase.residuals(x::VcovData) = x.residuals
StatsBase.dof_residual(x::VcovData) = x.dof_residual

##############################################################################
##
## Any type used for standard errors must define the following methods:
##
##############################################################################
function completecases(table, ::CovarianceEstimator)
	Tables.istable(table) || throw(ArgumentError("completecases requires a table input"))
	return trues(length(Tables.rows(table)))
end
materialize(table, v::CovarianceEstimator) = v
S_hat(x::RegressionModel, ::CovarianceEstimator) = error("S_hat not defined for this type")
df_FStat(x::RegressionModel, ::CovarianceEstimator, hasintercept::Bool) = dof_residual(x) - hasintercept

include("utils.jl")
include("covarianceestimators/vcovsimple.jl")
include("covarianceestimators/vcovrobust.jl")
include("covarianceestimators/vcovcluster.jl")
include("ranktest.jl")

include("precompile.jl")
_precompile_()

end


text:
using GeometricIntegrators.CommonFunctions
using GeometricIntegrators.Solutions
using GeometricIntegrators.Solutions: createHDF5
using GeometricIntegrators.TestProblems.KuboOscillatorProblem
using HDF5: HDF5File
using Test

nd = 3
ns = 5
nt = 10
Δt = .1
dt = Float64

t0    = 0.
x0    = rand(2)
y0    = rand(2)
q0    = rand(1)
p0    = rand(1)

t1    = 1.
x1    = rand(2, ns)
y1    = rand(2, ns)
q1    = rand(1, ns)
p1    = rand(1, ns)

t2    = t1 + (t1-t0)
x2    = rand(2, ns)
y2    = rand(2, ns)
q2    = rand(1, ns)
p2    = rand(1, ns)

tx    = zero(x0)
ty    = zero(y0)
tq    = zero(q0)
tp    = zero(p0)

xs    = rand(2,nt)
ys    = rand(2,nt)
qs    = rand(1,nt)
ps    = rand(1,nt)

Xs    = rand(2,nt,ns)
Ys    = rand(2,nt,ns)
Qs    = rand(1,nt,ns)
Ps    = rand(1,nt,ns)

sde  = kubo_oscillator_sde_1()
psde = kubo_oscillator_psde_1()

x100  = rand(ndims(sde),  100)
q100  = rand(ndims(psde), 100)
p100  = rand(ndims(psde), 100)

h5file = "test.hdf5"

@testset "$(rpad("Wiener Process",80))" begin

    wp = WienerProcess(dt, 1, nt, 1,  Δt, :strong)
    @test wp == WienerProcess(Δt, wp.ΔW[:], wp.ΔZ[:], :strong)
    @test ndims(wp) == 2

    wp = WienerProcess(dt, 1, nt, 1,  Δt, :weak)
    @test wp == WienerProcess(Δt, wp.ΔW[:], wp.ΔZ[:], :weak)
    @test ndims(wp) == 2

    wp = WienerProcess(dt, nd, nt, 1,  Δt, :strong)
    @test wp == WienerProcess(Δt, wp.ΔW, wp.ΔZ, :strong)
    @test ndims(wp) == 2

    wp = WienerProcess(dt, nd, nt, 1,  Δt, :weak)
    @test wp == WienerProcess(Δt, wp.ΔW, wp.ΔZ, :weak)
    @test ndims(wp) == 2

    wp = WienerProcess(dt, nd, nt, ns, Δt, :strong)
    @test wp == WienerProcess(Δt, wp.ΔW, wp.ΔZ, :strong)
    @test ndims(wp) == 3

    wp = WienerProcess(dt, nd, nt, ns, Δt, :weak)
    @test wp == WienerProcess(Δt, wp.ΔW, wp.ΔZ, :weak)
    @test ndims(wp) == 3

end

@testset "$(rpad("SDE Solution",80))" begin
    asol = AtomicSolution(sde)

    # test constructors and general functionality
    sol = Solution(sde, Δt, nt)
    @test typeof(sol) <: SolutionSDE

    sol0 = Solution(similar(sde, x0), Δt, nt)
    @test typeof(sol0) <: SolutionSDE

    sol1 = Solution(similar(sde, x1), Δt, nt)
    @test typeof(sol1) <: SolutionSDE

    @test sol != sol0
    @test sol != sol1

    @test sol == SSolutionSDE(sde, Δt, sol.W.ΔW, sol.W.ΔZ, sol.ntime)

    @test nsave(sol) == 1
    @test ntime(sol) == nt
    @test timesteps(sol) == Δt .* collect(0:nt)

    # test initial conditions
    set_initial_conditions!(sol, t0, x0)
    get_initial_conditions!(sol, tx, 1)
    get_initial_conditions!(sol, asol, 1)

    @test tx == x0
    @test asol.t == t0
    @test asol.q == x0

    δt, δx = get_initial_conditions(sol, 1)
    @test δt == t0
    @test δx == x0

    set_initial_conditions!(sol1, similar(sde, t1, x2))
    for i in 1:ns
        get_initial_conditions!(sol1, tx, i)
        @test tx == x2[:,i]
    end

    # test set/get solution
    sol1 = Solution(similar(kubo_oscillator_sde_2(), x0), Δt, nt)
    sol2 = Solution(similar(kubo_oscillator_sde_2(), x0), Δt, nt)
    for i in 1:nt
        tx .= xs[:,i]
        asol.q .= xs[:,i]
        set_solution!(sol1, tx, i)
        set_solution!(sol2, asol, i)
    end
    @test sol1.q[:,1:nt] == xs
    @test sol2.q[:,1:nt] == xs

    sol1 = Solution(similar(sde, x1), Δt, nt)
    sol2 = Solution(similar(sde, x1), Δt, nt)
    for i in 1:nt
        for k in 1:ns
            tx .= Xs[:,i,k]
            asol.q .= Xs[:,i,k]
            set_solution!(sol1, tx, i, k)
            set_solution!(sol2, asol, i, k)
        end
    end
    @test sol1.q[:,1:nt,:] == Xs
    @test sol2.q[:,1:nt,:] == Xs

    # test hdf5 in- and output
    sol1 = Solution(kubo_oscillator_sde_2(x0), Δt, nt)
    create_hdf5!(sol1, h5file)
    write_to_hdf5(sol1)
    close(sol1)
    @test isfile(h5file)
    sol2 = SSolutionSDE(h5file)
    @test sol1   != sol2
    @test sol1.t == sol2.t
    @test sol1.q == sol2.q
    @test sol1.ntime == sol2.ntime
    @test sol1.nsave == sol2.nsave
    @test sol1.nd == sol2.nd
    @test sol1.nm == sol2.nm
    @test sol1.ns == sol2.ns
    @test sol1.K  == sol2.K
    @test conv(sol1) == conv(sol2)
    rm(h5file)

    sol1 = Solution(kubo_oscillator_sde_2(x1), Δt, nt)
    create_hdf5!(sol1, h5file)
    write_to_hdf5(sol1)
    close(sol1)
    @test isfile(h5file)
    sol2 = SSolutionSDE(h5file)
    @test sol1   != sol2
    @test sol1.t == sol2.t
    @test sol1.q == sol2.q
    @test sol1.ntime == sol2.ntime
    @test sol1.nsave == sol2.nsave
    @test sol1.nd == sol2.nd
    @test sol1.nm == sol2.nm
    @test sol1.ns == sol2.ns
    @test sol1.K  == sol2.K
    @test conv(sol1) == conv(sol2)
    rm(h5file)

    sde1 = kubo_oscillator_sde_2(x1)
    sol1 = Solution(sde1, Δt, nt)
    create_hdf5!(sol1, h5file; save_W=false)
    write_to_hdf5(sol1)
    close(sol1)
    @test isfile(h5file)
    sol2 = SSolutionSDE(h5file)
    @test sol2.W.ΔW == zeros(eltype(q1), 0, 0, 0)
    @test sol2.W.ΔZ == zeros(eltype(q1), 0, 0, 0)
    rm(h5file)

    # test nsave and nwrite parameters
    sol = Solution(sde, Δt, 20, nsave=2)
    @test sol.nt == 10

    sol = Solution(sde, Δt, 20, nsave=2, nwrite=10)
    @test sol.nt == 5

    # test reset
    sol = Solution(sde, Δt, nt)
    reset!(sol)
    @test sol.t[0]   == t1
    @test sol.t[end] == t2
    @test offset(sol) == nt

    sde1 = kubo_oscillator_sde_3()
    sol1 = Solution(sde1, Δt, 100, nwrite=10)
    create_hdf5!(sol1, h5file)
    ΔW = zeros(sde1.m, 100)
    ΔZ = zeros(sde1.m, 100)
    for i in 1:10
        ΔW[:,(i-1)*10+1:i*10] .= sol1.W.ΔW
        ΔZ[:,(i-1)*10+1:i*10] .= sol1.W.ΔZ
        for j in 1:10
            set_solution!(sol1, x100[:,(i-1)*10+j], j)
        end
        write_to_hdf5(sol1)
        reset!(sol1)
    end
    close(sol1)
    sol2 = SSolutionSDE(h5file)
    @test sol2.t.t ≈ Δt .* collect(0:100) atol=1E-14
    @test sol2.q.d == hcat(reshape(sde1.q₀, (ndims(sde1),1)), x100)
    @test sol2.W.ΔW == ΔW
    @test sol2.W.ΔZ == ΔZ
    @test sol2.ntime == 100
    @test sol2.nsave == 1
    rm(h5file)

    sol1 = Solution(sde1, Δt, 100, nsave=2, nwrite=10)
    create_hdf5!(sol1, h5file)
    ΔW = zeros(sde1.m, 100)
    ΔZ = zeros(sde1.m, 100)
    for i in 1:10
        ΔW[:,(i-1)*10+1:i*10] .= sol1.W.ΔW
        ΔZ[:,(i-1)*10+1:i*10] .= sol1.W.ΔZ
        for j in 1:10
            set_solution!(sol1, x100[:,(i-1)*10+j], j)
        end
        write_to_hdf5(sol1)
        reset!(sol1)
    end
    close(sol1)
    sol2 = SSolutionSDE(h5file)
    @test sol2.t.t ≈ Δt .* collect(0:2:100) atol=1E-14
    @test sol2.q.d == hcat(reshape(sde1.q₀, (ndims(sde1),1)), x100)[:,1:2:end]
    @test sol2.W.ΔW == ΔW
    @test sol2.W.ΔZ == ΔZ
    @test sol2.ntime == 100
    @test sol2.nsave == 2
    rm(h5file)
end

@testset "$(rpad("PSDE Solution",80))" begin
    asol = AtomicSolution(psde)

    # test constructors and general functionality
    sol = Solution(psde, Δt, nt)
    @test typeof(sol) <: SolutionPSDE

    sol0 = Solution(similar(psde, q0, p0), Δt, nt)
    @test typeof(sol0) <: SolutionPSDE

    sol1 = Solution(similar(psde, q1, p1), Δt, nt)
    @test typeof(sol1) <: SolutionPSDE

    @test sol != sol0
    @test sol != sol1

    @test sol == SSolutionPSDE(psde, Δt, sol.W.ΔW, sol.W.ΔZ, sol.ntime)

    @test nsave(sol) == 1
    @test ntime(sol) == nt
    @test timesteps(sol) == Δt .* collect(0:nt)

    # test initial conditions
    set_initial_conditions!(sol, t0, q0, p0)
    get_initial_conditions!(sol, tq, tp, 1)
    get_initial_conditions!(sol, asol, 1)

    @test tq == q0
    @test tp == p0
    @test asol.t == t0
    @test asol.q == q0
    @test asol.p == p0

    δt, δq, δp = get_initial_conditions(sol, 1)
    @test δt == t0
    @test δq == q0
    @test δp == p0

    set_initial_conditions!(sol1, similar(psde, t1, q2, p2))
    for i in 1:ns
        get_initial_conditions!(sol1, tq, tp, i)
        @test tq == q2[:,i]
        @test tp == p2[:,i]
    end

    # test set/get solution
    sol1 = Solution(similar(kubo_oscillator_psde_2(), q0, p0), Δt, nt)
    sol2 = Solution(similar(kubo_oscillator_psde_2(), q0, p0), Δt, nt)
    for i in 1:nt
        tq .= qs[:,i]
        tp .= ps[:,i]
        asol.q .= qs[:,i]
        asol.p .= ps[:,i]
        set_solution!(sol1, tq, tp, i)
        set_solution!(sol2, asol, i)
    end
    @test sol1.q[:,1:nt] == qs
    @test sol1.p[:,1:nt] == ps
    @test sol2.q[:,1:nt] == qs
    @test sol2.p[:,1:nt] == ps

    sol1 = Solution(similar(kubo_oscillator_psde_2(), q1, p1), Δt, nt)
    sol2 = Solution(similar(kubo_oscillator_psde_2(), q1, p1), Δt, nt)
    for i in 1:nt
        for k in 1:ns
            tq .= Qs[:,i,k]
            tp .= Ps[:,i,k]
            asol.q .= Qs[:,i,k]
            asol.p .= Ps[:,i,k]
            set_solution!(sol1, tq, tp, i, k)
            set_solution!(sol2, asol, i, k)
        end
    end
    @test sol1.q[:,1:nt,:] == Qs
    @test sol1.p[:,1:nt,:] == Ps
    @test sol2.q[:,1:nt,:] == Qs
    @test sol2.p[:,1:nt,:] == Ps

    # test hdf5 in- and output
    sol1 = Solution(kubo_oscillator_psde_2(q0, p0), Δt, nt)
    create_hdf5!(sol1, h5file)
    write_to_hdf5(sol1)
    close(sol1)
    @test isfile(h5file)
    sol2 = SSolutionPSDE(h5file)
    @test sol1   != sol2
    @test sol1.t == sol2.t
    @test sol1.q == sol2.q
    @test sol1.p == sol2.p
    @test sol1.ntime == sol2.ntime
    @test sol1.nsave == sol2.nsave
    @test sol1.nd == sol2.nd
    @test sol1.nm == sol2.nm
    @test sol1.ns == sol2.ns
    @test sol1.K  == sol2.K
    @test conv(sol1) == conv(sol2)
    rm(h5file)

    sol1 = Solution(kubo_oscillator_psde_2(q1, p1), Δt, nt)
    create_hdf5!(sol1, h5file)
    write_to_hdf5(sol1)
    close(sol1)
    @test isfile(h5file)
    sol2 = SSolutionPSDE(h5file)
    @test sol1   != sol2
    @test sol1.t == sol2.t
    @test sol1.q == sol2.q
    @test sol1.p == sol2.p
    @test sol1.ntime == sol2.ntime
    @test sol1.nsave == sol2.nsave
    @test sol1.nd == sol2.nd
    @test sol1.nm == sol2.nm
    @test sol1.ns == sol2.ns
    @test sol1.K  == sol2.K
    @test conv(sol1) == conv(sol2)
    rm(h5file)

    psde1 = kubo_oscillator_psde_2(q1, p1)
    sol1 = Solution(psde1, Δt, nt)
    create_hdf5!(sol1, h5file; save_W=false)
    write_to_hdf5(sol1)
    close(sol1)
    @test isfile(h5file)
    sol2 = SSolutionPSDE(h5file)
    @test sol2.W.ΔW == zeros(eltype(q1), 0, 0, 0)
    @test sol2.W.ΔZ == zeros(eltype(q1), 0, 0, 0)
    rm(h5file)

    # test nsave and nwrite parameters
    sol = Solution(psde, Δt, 20, nsave=2)
    @test sol.nt == 10

    sol = Solution(psde, Δt, 20, nsave=2, nwrite=10)
    @test sol.nt == 5

    # test reset
    sol = Solution(psde, Δt, nt)
    reset!(sol)
    @test sol.t[0]   == t1
    @test sol.t[end] == t2
    @test offset(sol) == nt

    psde1 = kubo_oscillator_psde_3()
    sol1 = Solution(psde1, Δt, 100, nwrite=10)
    create_hdf5!(sol1, h5file)
    ΔW = zeros(psde1.m, 100)
    ΔZ = zeros(psde1.m, 100)
    for i in 1:10
        ΔW[:,(i-1)*10+1:i*10] .= sol1.W.ΔW
        ΔZ[:,(i-1)*10+1:i*10] .= sol1.W.ΔZ
        for j in 1:10
            oj = (i-1)*10+j
            set_solution!(sol1, q100[:,oj], p100[:,oj], j)
        end
        write_to_hdf5(sol1)
        reset!(sol1)
    end
    close(sol1)
    sol2 = SSolutionPSDE(h5file)
    @test sol2.t.t ≈ Δt .* collect(0:100) atol=1E-14
    @test sol2.q.d == hcat(reshape(psde1.q₀, (ndims(psde1),1)), q100)
    @test sol2.p.d == hcat(reshape(psde1.p₀, (ndims(psde1),1)), p100)
    @test sol2.W.ΔW == ΔW
    @test sol2.W.ΔZ == ΔZ
    @test sol2.ntime == 100
    @test sol2.nsave == 1
    rm(h5file)

    sol1 = Solution(psde1, Δt, 100, nsave=2, nwrite=10)
    create_hdf5!(sol1, h5file)
    ΔW = zeros(psde1.m, 100)
    ΔZ = zeros(psde1.m, 100)
    for i in 1:10
        ΔW[:,(i-1)*10+1:i*10] .= sol1.W.ΔW
        ΔZ[:,(i-1)*10+1:i*10] .= sol1.W.ΔZ
        for j in 1:10
            oj = (i-1)*10+j
            set_solution!(sol1, q100[:,oj], p100[:,oj], j)
        end
        write_to_hdf5(sol1)
        reset!(sol1)
    end
    close(sol1)
    sol2 = SSolutionPSDE(h5file)
    @test sol2.t.t ≈ Δt .* collect(0:2:100) atol=1E-14
    @test sol2.q.d == hcat(reshape(psde1.q₀, (ndims(psde1),1)), q100)[:,1:2:end]
    @test sol2.p.d == hcat(reshape(psde1.p₀, (ndims(psde1),1)), p100)[:,1:2:end]
    @test sol2.W.ΔW == ΔW
    @test sol2.W.ΔZ == ΔZ
    @test sol2.ntime == 100
    @test sol2.nsave == 2
    rm(h5file)
end

@testset "$(rpad("SPSDE Solution",80))" begin
    spsde  = kubo_oscillator_spsde_1()
    sol = Solution(spsde, Δt, nt)
    @test typeof(sol) <: SolutionPSDE
end