text:
#=
Utilities for combining and manipulating TS objects with their indexes
=#

import Base: hcat, vcat, merge

# function partner(x::TS, y::TS)
#     yy = !overlaps(x.index, y.index) .* NaN
#     yy[!isnan(yy),:] = y.values
#     return ts(yy, x.index, y.fields)
# end

"""
    ojoin(x::TS,y::TS)::TS

Outer join two TS objects by index.

Equivalent to `x` OUTER JOIN `y` ON `x.index` = `y.index`.

...
# Arguments
- `x::TS`: Left side of the join.
- `y::TS`: Right side of the join.
...
"""
function ojoin(x::TS, y::TS)::TS
    if isempty(x) && !isempty(y)
        return y
    elseif isempty(y) && !isempty(x)
        return x
    elseif isempty(x) && isempty(y)
        return ts()
    end
    idx = union(x.index, y.index)
    xna = setdiff(idx, x.index)
    yna = setdiff(idx, y.index)
    xi = sortperm(unique([x.index; xna]))
    yi = sortperm(unique([y.index; yna]))
    xvals = [x.values; fill(NaN, (length(xna), size(x,2)))][xi,:]
    yvals = [y.values; fill(NaN, (length(yna), size(y,2)))][yi,:]
    return ts([xvals yvals], sort(idx), [x.fields; y.fields])
end

"""
    ijoin(x::TS,y::TS)::TS

Inner join two TS objects by index.

Equivalent to `x` INNER JOIN `y` on `x.index` = `y.index`.

...
# Arguments
- `x::TS`: Left side of the join.
- `y::TS`: Right side of the join.
...
"""
function ijoin(x::TS, y::TS)::TS
    if isempty(x) && !isempty(y)
        return y
    elseif isempty(y) && !isempty(x)
        return x
    elseif isempty(x) && isempty(y)
        return ts()
    end
    idx = intersect(x.index, y.index)
    return ts([x[idx].values y[idx].values], idx, [x.fields; y.fields])
end

"""
    ljoin(x::TS, y::TS)::TS

Left join two TS objects by index.

Equivalent to `x` LEFT JOIN `y` ON `x.index` = `y.index`.

...
# Arguments
- `x::TS`: Left side of the join.
- `y::TS`: Right side of the join.
...
"""
function ljoin(x::TS, y::TS)::TS
    return [x y[intersect(x.index, y.index)]]
end

"""
    rjoin(x::TS, y::TS)::TS

Right join two TS objects by index.

Equivalent to `x` RIGHT JOIN `y` ON `x.index` = `y.index`.

...
# Arguments
- `x::TS`: Left side of the join.
- `y::TS`: Right side of the join.
...
"""
function rjoin(x::TS, y::TS)::TS
    return [x[intersect(x.index, y.index)] y]
end

hcat(x::TS, y::TS)::TS = ojoin(x, y)
hcat(x::TS)::TS = x
function hcat(series::TS...)
    out = series[1]
    @inbounds for j = 2:length(series)
        out = [out series[j]]
    end
    return out
end
function vcat(x::TS, y::TS)
    @assert size(x,2) == size(y,2) "Dimension mismatch: Number of columns must be equal."
    return TS([x.values;y.values], [x.index;y.index], x.fields)
end
function vcat(series::TS...)
    out = series[1]
    @inbounds for j = 2:length(series)
        out = vcat(out, series[j])
    end
    return out
end

"""
    merge(x::TS,y::TS;join::Char='o')::TS

Merge two time series objects together by index with an optionally specified join type parameter.

...
# Arguments
- `x::TS`: Left side of the join.
- `y::TS`: Right side of the join.
Optional args:
- `join::Char='o'::TS`: Specifies the logic used to perform the merge, and may take on the values 'o' (outer join), 'i' (inner join), 'l' (left join), or 'r' (right join). Defaults to outer join, whose result is the same as `hcat(x, y)` or `[x y]`.
...
"""
function merge(x::TS, y::TS; join::Char='o')::TS
    @assert join == 'o' || join == 'i' || join == 'l' || join == 'r' "`join` must be 'o', 'i', 'l', or 'r'."
    if join == 'o'
        return ojoin(x, y)
    elseif join == 'i'
        return ijoin(x, y)
    elseif join == 'l'
        return ljoin(x, y)
    elseif join == 'r'
        return rjoin(x, y)
    end
end

#===============================================================================
                COMBINING/MERGING WITH OTHER TYPES
===============================================================================#
hcat(x::TS, y::AbstractArray)::TS = ojoin(x, ts(y, x.index))
hcat(y::AbstractArray, x::TS)::TS = ojoin(ts(y, x.index), x)
hcat(x::TS, y::Number)::TS = ojoin(x, ts(fill(y,size(x,1)), x.index))
hcat(y::Number, x::TS)::TS = ojoin(ts(fill(y,size(x,1)), x.index), x)

function getmaxtype(arrs)::Type
    result::Type = Any
    @inbounds for k in 2:length(arrs)
        result = eltype(promote(arrs[k-1][1], arrs[k][1]))
    end
    return result
end

function hcat(series::TS, arrs::AbstractArray...)::TS
    n = size(series,1)
    cols = map(arr->size(arr,2), arrs)
    rows = map(arr->size(arr,1), arrs)
    @assert all(rows.==n) "All arrays must have same number of rows as TS object."
    out = zeros(getmaxtype(arrs), (n, sum(cols)))
    first_col, last_col = 0, 0
    @inbounds for j in 1:length(arrs)
        last_col = first_col + cols[j]
        first_col = last_col - first_col
        out[:,first_col:last_col] = arrs[j]
    end
    return [series out]
end

function hcat(series::TS, nums::Number...)::TS
    n = size(series,1)
    k = length(nums)
    T::Type = getmaxtype(nums)
    out = zeros(T, (n,k))
    @inbounds for j in 1:k
        out[:,j] = ones(T,n) * nums[j]
    end
    return [series out]
end


text:
"""
Copyright Marcus Greiff 2020

Tests of the discontinuous attitude FSF controller on SU(2)
"""

# Include modules and recompile if running test standalone
if !isdefined(Main, :runAllTests)
    using Test, AerialVehicleControl, LinearAlgebra
    recompile();
end

# Testing fixture
function fixture()
    J  = rand(3,3); J  = J + J'; J  = J + I*maximum(abs.(eigvals(J)));
    kR, kc, kw = 1.0, 1.0, 1.0;
    return J, kR, kc, kw
end

@verbose(1, "Testing FSF utilities...")
@testset "Test FSF utilities" begin
    @testset "SU2 - Feasibility check - Non-symmetric inertia" begin
        # Fixture
        J, kX, kc, kw = fixture()
        J[1,2], J[2,1] = 2.0, 1.0;
        C      = con_state_qw_fsf_t(J, kX, kc, kw);
        C_copy = deepcopy(C)

        # Test
        status = ccall((:assert_attitude_FSF_SU2, AerialVehicleControl.CONT_LIB_PATH),
            Cint,
            (Ref{con_state_qw_fsf_t},),
            C)

        @test isapprox(status, 0);
        allFieldsEqual(C, C_copy, true);
    end

    @testset "SU2 - Feasibility check - Non-PSD inertia" begin
        # Fixture
        J, kX, kc, kw = fixture()
        J  = J - I*1.01*minimum(abs.(eigvals(J)));
        C      = con_state_qw_fsf_t(J, kX, kc, kw);
        C_copy = deepcopy(C)

        # Test
        status = ccall((:assert_attitude_FSF_SU2, AerialVehicleControl.CONT_LIB_PATH),
            Cint,
            (Ref{con_state_qw_fsf_t},),
            C)

        @test isapprox(status, 0);
        allFieldsEqual(C, C_copy, true);
    end
    @testset "SU2 - Feasibility check - Controller parameters - kc large" begin
        # Fixture
        J, kX, _, kw = fixture()
        kc     = 1.01 * maximal_feasible_kc_SU2(kX, kw, J);
        C      = con_state_qw_fsf_t(J, kX, kc, kw);
        C_copy = deepcopy(C)

        # Test
        status = ccall((:assert_attitude_FSF_SU2, AerialVehicleControl.CONT_LIB_PATH),
            Cint,
            (Ref{con_state_qw_fsf_t},),
            C)

        # Assert
        @test isapprox(status, 0);
        allFieldsEqual(C, C_copy, true);
    end
    @testset "SU2 - Feasibility check - Controller parameters - kc small" begin
        # Fixture
        J, kX, _, kw = fixture()
        kc     = 0.99 * maximal_feasible_kc_SU2(kX, kw, J);
        C      = con_state_qw_fsf_t(J, kX, kc, kw);
        C_copy = deepcopy(C)

        # Test
        status = ccall((:assert_attitude_FSF_SU2, AerialVehicleControl.CONT_LIB_PATH),
            Cint,
            (Ref{con_state_qw_fsf_t},),
            C)

        # Assert
        @test isapprox(status, 1);
        allFieldsEqual(C, C_copy, true);
    end
    @testset "SO3 - Feasibility check - Non-symmetric inertia" begin
        # Fixture
        J, kR, kc, kw = fixture()
        J[1,2], J[2,1] = 2.0, 1.0;
        C      = con_state_qw_fsf_t(J, kR, kc, kw);
        C_copy = deepcopy(C)

        # Test
        status = ccall((:assert_attitude_FSF_SO3, AerialVehicleControl.CONT_LIB_PATH),
            Cint,
            (Ref{con_state_qw_fsf_t},),
            C)

        @test isapprox(status, 0);
        allFieldsEqual(C, C_copy, true);
    end

    @testset "SO3 - Feasibility check - Non-PSD inertia" begin
        # Fixture
        J, kR, kc, kw = fixture()
        J  = J - I*1.01*minimum(abs.(eigvals(J)));
        C      = con_state_qw_fsf_t(J, kR, kc, kw);
        C_copy = deepcopy(C)

        # Test
        status = ccall((:assert_attitude_FSF_SO3, AerialVehicleControl.CONT_LIB_PATH),
            Cint,
            (Ref{con_state_qw_fsf_t},),
            C)

        @test isapprox(status, 0);
        allFieldsEqual(C, C_copy, true);
    end
    @testset "SO3 - Feasibility check - Controller parameters - kc large" begin
        # Fixture
        J, kR, _, kw = fixture()
        kc     = 1.01 * maximal_feasible_kc_SO3(kR, kw, J);
        C      = con_state_qw_fsf_t(J, kR, kc, kw);
        C_copy = deepcopy(C)

        # Test
        status = ccall((:assert_attitude_FSF_SO3, AerialVehicleControl.CONT_LIB_PATH),
            Cint,
            (Ref{con_state_qw_fsf_t},),
            C)

        # Assert
        @test isapprox(status, 0);
        allFieldsEqual(C, C_copy, true);
    end
    @testset "SO3 - Feasibility check - Controller parameters - kc small" begin
        # Fixture
        J, kR, _, kw = fixture()
        kc     = 0.99 * maximal_feasible_kc_SO3(kR, kw, J);
        C      = con_state_qw_fsf_t(J, kR, kc, kw);
        C_copy = deepcopy(C)

        # Test
        status = ccall((:assert_attitude_FSF_SO3, AerialVehicleControl.CONT_LIB_PATH),
            Cint,
            (Ref{con_state_qw_fsf_t},),
            C)

        # Assert
        @test isapprox(status, 1);
        allFieldsEqual(C, C_copy, true);
    end
end

ext: jl


text:
module SynapseTables
using DataFrames
using Query
using ..Utils.BoundingBoxes
using OffsetArrays
using ProgressMeter

export SynapseTable
const SynapseTable = DataFrame

function preprocess(self::SynapseTable, voxelSize::Tuple)
    @assert length(voxelSize) == 3
    voxelSize = map(Float32, voxelSize)

    println("drop missing values...")
	@time DataFrames.dropmissing!(self)

    # this is not working due to a DataFrames bug
    # https://github.com/JuliaData/DataFrames.jl/issues/1862
    # we need to add true to make it work as old ways...
	@showprogress 1 "transform datatype to Float32..." for (key, value) in DataFrames.eachcol(self, true)
		#if key!=:presyn_wt && key!=:postsyn_wt
		#	self[key] = round.(Int, value)
        #end
        self[!, key] = Vector{Float32}(value)
    end

    println("remove self connections...")
    nrow1 = DataFrames.nrow(self)
    self = @from i in self begin
        @where i.presyn_segid != i.postsyn_segid
        @select i
        @collect DataFrame
    end
    nrow2 = DataFrames.nrow(self)
    println("removed $(nrow1-nrow2) self connections in $nrow1 synaptic connections.")

    try
        DataFrames.rename!(self, [  :centroid_x => :psd_x,
                                    :centroid_y => :psd_y,
                                    :centroid_z => :psd_z])
    catch err
        @info("did not find centroid column in the synapse table file.")
    end

    println("transform to physical coordinate...")
    self[!, :BBOX_bx]      = (self[!, :BBOX_bx]   .- one(Float32)) .* voxelSize[1]
    self[!, :BBOX_by]      = (self[!, :BBOX_by]   .- one(Float32)) .* voxelSize[2]
    self[!, :BBOX_bz]      = (self[!, :BBOX_bz]   .- one(Float32)) .* voxelSize[3]

    self[!, :BBOX_ex]      = (self[!, :BBOX_ex]   .- one(Float32)) .* voxelSize[1]
    self[!, :BBOX_ey]      = (self[!, :BBOX_ey]   .- one(Float32)) .* voxelSize[2]
    self[!, :BBOX_ez]      = (self[!, :BBOX_ez]   .- one(Float32)) .* voxelSize[3]

    self[!, :psd_x]        = (self[!, :psd_x]     .-one(Float32)) .* voxelSize[1]
    self[!, :psd_y]        = (self[!, :psd_y]     .-one(Float32)) .* voxelSize[2]
    self[!, :psd_z]        = (self[!, :psd_z]     .-one(Float32)) .* voxelSize[3]

    self[!, :presyn_x]     = (self[!, :presyn_x]  .- one(Float32)) .* voxelSize[1]
    self[!, :presyn_y]     = (self[!, :presyn_y]  .- one(Float32)) .* voxelSize[2]
    self[!, :presyn_z]     = (self[!, :presyn_z]  .- one(Float32)) .* voxelSize[3]

    self[!, :postsyn_x]    = (self[!, :postsyn_x] .- one(Float32)) .* voxelSize[1]
    self[!, :postsyn_y]    = (self[!, :postsyn_y] .- one(Float32)) .* voxelSize[2]
    self[!, :postsyn_z]    = (self[!, :postsyn_z] .- one(Float32)) .* voxelSize[3]
	return self
end

function postprocess(self::SynapseTable, voxelSize::Union{Vector, Tuple})
    @assert length(voxelSize) == 3
    voxelSize = map(Float32, voxelSize)

    println("transform to voxel coordinate...")
    self[!, :BBOX_bx]      = (self[!, :BBOX_bx]./voxelSize[1]) .+ one(Float32)
    self[!, :BBOX_by]      = (self[!, :BBOX_by]./voxelSize[2]) .+ one(Float32)
    self[!, :BBOX_bz]      = (self[!, :BBOX_bz]./voxelSize[3]) .+ one(Float32)

    self[!, :BBOX_ex]      = (self[!, :BBOX_ex]./voxelSize[1]) .+ one(Float32)
    self[!, :BBOX_ey]      = (self[!, :BBOX_ey]./voxelSize[2]) .+ one(Float32)
    self[!, :BBOX_ez]      = (self[!, :BBOX_ez]./voxelSize[3]) .+ one(Float32)

    self[!, :psd_x]        = (self[!, :psd_x]./voxelSize[1]) .+ one(Float32)
    self[!, :psd_y]        = (self[!, :psd_y]./voxelSize[2]) .+ one(Float32)
    self[!, :psd_z]        = (self[!, :psd_z]./voxelSize[3]) .+ one(Float32)

    self[!, :presyn_x]     = (self[!, :presyn_x]./voxelSize[1]) .+ one(Float32)
    self[!, :presyn_y]     = (self[!, :presyn_y]./voxelSize[2]) .+ one(Float32)
    self[!, :presyn_z]     = (self[!, :presyn_z]./voxelSize[3]) .+ one(Float32)

    self[!, :postsyn_x]     = (self[!, :postsyn_x]./voxelSize[1]) .+ one(Float32)
    self[!, :postsyn_y]     = (self[!, :postsyn_y]./voxelSize[2]) .+ one(Float32)
    self[!, :postsyn_z]     = (self[!, :postsyn_z]./voxelSize[3]) .+ one(Float32)

    return self
end

function get_coordinate_array(self::SynapseTable, prefix::String)
	coordinateNames = get_coordinate_names(prefix)
	hcat(	self[!, coordinateNames[1]],
			self[!, coordinateNames[2]],
			self[!, coordinateNames[3]])
end

function get_coordinate(self::DataFrameRow, prefix::String)
	coordinateNames = get_coordinate_names(prefix)
    map(x->self[x], coordinateNames)
end

function get_pre_synapses_of_a_neuron(self::SynapseTable, neuronId::Int)
    ret =  @from i in self begin
        @where i.presyn_segid == neuronId
        @select i
        @collect DataFrame
    end
    ret
end

function get_post_synapses_of_a_neuron(self::SynapseTable, neuronId::Int)
    ret =  @from i in self begin
        @where i.postsyn_segid == neuronId
        @select i
        @collect DataFrame
    end
    ret
end

function get_synapses_of_a_neuron(self::SynapseTable, neuronId::Int)
	ret =  @from i in self begin
        @where i.postsyn_segid == neuronId || i.presyn_segid == neuronId
        @select i
        @collect DataFrame
    end
    ret
end

@inline function get_coordinate_names( prefix::String )
    map(x->Symbol(prefix*x), ("x", "y", "z"))
end

function initialize_mask(self::SynapseTable, voxelSize::NTuple{3,Int};
                coordinatePrefixList = ["presyn_", "postsyn_"], T::DataType = Bool)
	@assert !isempty(self)
    range = mapreduce(x->BoundingBox(self,x), union, coordinatePrefixList) |> BoundingBoxes.get_unit_range
    range = map((r,s)->fld(r.start, s):cld(r.stop, s), range, voxelSize)
    mask = OffsetArray{T}(undef, range...)
    fill!(mask, zero(T))
    mask
end

function get_mask(self::SynapseTable;
                  voxelSize::Union{Vector, Tuple} = (1000, 1000, 1000) #=nm=#,
                coordinatePrefixList::Vector{String} = ["presyn_", "postsyn_"],
                mask::OffsetArray{T,3,Array{T,3}} = initialize_mask(self, voxelSize;
                    coordinatePrefixList = coordinatePrefixList, T=Float32)) where T
    @assert !isempty(self)
    @assert length(voxelSize) == 3
    voxelSize = map(Float32, voxelSize)
	for coordinatePrefix in coordinatePrefixList
		coordinateNames = SynapseTables.get_coordinate_names( coordinatePrefix )
		coordinateList = map((c,s)->self[!, c]/s, coordinateNames, voxelSize)
        coordinateList = map(x->round.(Int, x), coordinateList)
		for i in 1:length(coordinateList[1])
			mask[   coordinateList[1][i],
				    coordinateList[2][i],
					coordinateList[3][i] ] = one(eltype(mask))
		end
	end
    mask
end

@inline function BoundingBox(self::SynapseTable, prefix::String)
    coordinateNames = get_coordinate_names( prefix )
    BoundingBox(self, coordinateNames)
end

@inline function BoundingBox(self::SynapseTable, coordinateNames::NTuple{3,Symbol})
    if isempty(self)
        @warn("empty synapse table!")
        return nothing
    else
        minCorner = map(x->minimum(self[!, x]), coordinateNames)
        maxCorner = map(x->maximum(self[!, x]), coordinateNames)
        return BoundingBoxes.BoundingBox(minCorner, maxCorner)
    end
end

end # module


text:
@testset "Exceptions" begin
  @test_throws DomainError SimpleBirthDeathProcess.log1mexp(0.01)
  @test_throws DomainError SimpleBirthDeathProcess.logexpm1(-0.01)
end

@testset "Type stability (log1mexp)" begin
  @testset for F = [Float32, Float64, BigFloat]
    @test isa(SimpleBirthDeathProcess.log1mexp(F(-1)), F)
    @test isa(SimpleBirthDeathProcess.log1mexp(F(-0.5)), F)
    @test isa(SimpleBirthDeathProcess.log1mexp(F(0)), F)
  end
end

@testset "Type stability (logexpm1)" begin
  @testset for F = [Float32, Float64, BigFloat]
    @test isa(SimpleBirthDeathProcess.log1pexp(F(40)), F)
    @test isa(SimpleBirthDeathProcess.log1pexp(F(20)), F)
    @test isa(SimpleBirthDeathProcess.log1pexp(F(10)), F)
  end
end

@testset "Type stability (log1pexp)" begin
  @testset for F = [Float32, Float64, BigFloat]
    @test isa(SimpleBirthDeathProcess.log1pexp(F(-40)), F)
    @test isa(SimpleBirthDeathProcess.log1pexp(F(-10)), F)
    @test isa(SimpleBirthDeathProcess.log1pexp(F(10)), F)
    @test isa(SimpleBirthDeathProcess.log1pexp(F(40)), F)
  end
end

@testset "Accuracy (log1mexp)" begin
  @test SimpleBirthDeathProcess.log1mexp(-10.0) ≈ -0.45400960370489209504e-4
  @test SimpleBirthDeathProcess.log1mexp(-1.0) ≈ -0.45867514538708189103
  @test SimpleBirthDeathProcess.log1mexp(-0.7) ≈ -0.68634100280838510968
  @test SimpleBirthDeathProcess.log1mexp(-0.69) ≈ -0.69630429714405665212
  @test SimpleBirthDeathProcess.log1mexp(-0.01) ≈ -4.6101660193248969177
end

@testset "Accuracy (logexpm1)" begin
  @test SimpleBirthDeathProcess.logexpm1(0.5) ≈ -0.432752129567188571894641000155
  @test SimpleBirthDeathProcess.logexpm1(1.0) ≈ 0.541324854612918108978356354931
  @test SimpleBirthDeathProcess.logexpm1(10.0) ≈ 9.99995459903962951079049555364
  @test SimpleBirthDeathProcess.logexpm1(19.0) ≈ 18.9999999943972035467670684411
  @test SimpleBirthDeathProcess.logexpm1(20.0) ≈ 19.9999999979388463754372650415
  @test SimpleBirthDeathProcess.logexpm1(30.0) ≈ 29.9999999999999064237703115939
  @test SimpleBirthDeathProcess.logexpm1(37.0) ≈ 36.9999999999999999146695237426
  @test SimpleBirthDeathProcess.logexpm1(40.0) ≈ 39.9999999999999999957516457447
  @test SimpleBirthDeathProcess.logexpm1(50.0) ≈ 49.9999999999999999999998071250
end

@testset "Accuracy (log1pexp)" begin
  @test SimpleBirthDeathProcess.log1pexp(-40.0) ≈ 4.24835425528999999097574306078e-18
  @test SimpleBirthDeathProcess.log1pexp(-37.0) ≈ 8.53304762574399963593549108392e-17
  @test SimpleBirthDeathProcess.log1pexp(-10.0) ≈ 0.453988992168646467694878337464e-4
  @test SimpleBirthDeathProcess.log1pexp(0.0) ≈ 0.693147180559945309417232121458
  @test SimpleBirthDeathProcess.log1pexp(10.0) ≈ 10.0000453988992168646467694878
  @test SimpleBirthDeathProcess.log1pexp(18.0) ≈ 18.0000000152299796287364881015
  @test SimpleBirthDeathProcess.log1pexp(20.0) ≈ 20.0000000020611536203143807032
  @test SimpleBirthDeathProcess.log1pexp(30.0) ≈ 30.0000000000000935762296883974
  @test SimpleBirthDeathProcess.log1pexp(40.0) ≈ 40.0000000000000000042483542553
end


text:
# ------------------------------------
# Tests adapted from Git2Go Library
# ------------------------------------
context("test lookup tree") do
    test_path = joinpath(pwd(), "testrepo")
    repo = create_test_repo(test_path)
    try
        @test isdir(workdir(repo))
        idx = GitIndex(repo)
        @test isa(idx, GitIndex)
        push!(idx, "README")
        tree_id = write_tree!(idx)
        @test isa(tree_id, Oid)
        @test isequal(hex(tree_id), "b7119b11e8ef7a1a5a34d3ac87f5b075228ac81e")
        tree1 = lookup(GitTree, repo, tree_id)
        tree2 = lookup_tree(repo, tree_id)

        @test tree1 == tree2
        tree = tree1
        @test isa(tree, GitTree)
        @test tree.ptr != C_NULL
        @test isa(Oid(tree), Oid)
        @test isa(hex(tree), ASCIIString)
    finally
        close(repo)
        LibGit2.free!(repo)
        Base.gc()
        rm(test_path, recursive=true)
    end
end

# -----------------------------------------
# Tests adapted from Ruby's Rugged Library
# -----------------------------------------

function new_idx_entry()
    now = time()
    return GitIndexEntry("new_path",
                      Oid("d385f264afb75a56a5bec74243be9b367ba4ca08"),
                      now,
                      now,
                      1000,
                      234881027,
                      88888,
                      33188,
                      502,
                      502,
                      false,
                      3)
end

with_test_index("test index size") do test_index, path
    @test length(test_index) == 2
end

with_test_index("test empty index") do test_index, path
    clear!(test_index)
    @test length(test_index) == 0
end

with_test_index("test remove entries") do test_index, path
    remove!(test_index, "new.txt")
    @test length(test_index) == 1
end

with_test_index("test rmeove dir") do test_index, path
    removedir!(test_index, "does-not-exist")
    @test length(test_index) == 2

    removedir!(test_index, "", 2)
    @test length(test_index) == 2

    removedir!(test_index, "")
    @test length(test_index) == 0
end

with_test_index("test get entry data") do test_index, path
    entry = test_index[1]
    @test "README" == entry.path
    @test Oid("1385f264afb75a56a5bec74243be9b367ba4ca08") == Oid(entry)
    @test 1273360380 == int(entry.mtime)
    @test 1273360380 == int(entry.ctime)
    @test 4 == entry.file_size

    @test 234881026 == entry.dev
    @test 6674088 == entry.ino
    @test 33188 == entry.mode
    @test 501 == entry.uid
    @test 0 == entry.gid
    @test false == entry.valid
    @test 0 == entry.stage

    entry = test_index[2]
    @test "new.txt" == entry.path
    @test Oid("fa49b077972391ad58037050f2a75f74e3671e92") == Oid(entry)
end

with_test_index("test iterate entries") do test_index, path
    c = sort(collect(test_index), by=x->Oid(x))
    @test join(map(x->x.path, c), ":") == "README:new.txt"
end

with_test_index("test update entries") do test_index, path
    now = int(time())
    entry = test_index[1]
    Oid(entry) = Oid("12ea3153a78002a988bb92f4123e7e831fd1138a")
    entry.mtime = now
    entry.ctime = now
    entry.file_size = 1000
    entry.dev = 234881027
    entry.ino = 88888
    entry.mode = 33188
    entry.uid = 502
    entry.gid = 502
    entry.stage = 3
    push!(test_index, entry)
    newentry = test_index[entry.path, 3]
    @test isa(newentry, GitIndexEntry)
    #TODO: rounding in times causes this test to fail
    #@test newentry == entry
end

with_test_index("test add new entries") do test_index, path
    push!(test_index, new_idx_entry())
    @test length(test_index) == 3
    c = sort(collect(test_index), by=x->Oid(x))
    @test join(map(x -> x.path, c), ":") == "README:new_path:new.txt"
end

with_test_index("test can write index") do test_index, path
    tmp_path, tmp_io = mktemp()
    try
        write(tmp_io, readall(path))
        close(tmp_io)
        test_index1 = GitIndex(tmp_path)

        entry = new_idx_entry()
        push!(test_index1, entry)
        entry.path = "else.txt"
        push!(test_index1, entry)
        write!(test_index1)

        test_index2 = GitIndex(tmp_path)
        c = sort(collect(test_index2), by=x->Oid(x))
        @test join(map(x -> x.path, c), ":") == "README:else.txt:new_path:new.txt"
        @test length(test_index2) == 4
    finally
        close(tmp_io)
        rm(tmp_path, recursive=true)
    end
end

context("test adding a path") do
    tmp_path = mktempdir()
    test_repo = init_repo(tmp_path, bare=false)
    try
        test_index1 = GitIndex(test_repo)
        fh = open(joinpath(tmp_path, "test.txt"), "w")
        write(fh, "test content")
        close(fh)
        push!(test_index1, "test.txt")
        write!(test_index1)

        test_index2 = GitIndex(joinpath(tmp_path, ".git", "index"))
        @test test_index2[1].path == "test.txt"
    finally
        close(test_repo)
        LibGit2.free!(test_repo)
        Base.gc()
        rm(tmp_path, recursive=true)
    end
end

context("test reloading index") do
    tmp_path = mktempdir()
    test_repo = init_repo(tmp_path, bare=false)
    try
        index = GitIndex(test_repo)
        fh = open(joinpath(tmp_path, "test.txt"), "w")
        write(fh, "test content")
        close(fh)
        push!(index, "test.txt")
        write!(index)

        rindex = GitIndex(joinpath(tmp_path, ".git", "index"))
        entry = rindex["test.txt"]
        @test entry.stage == 0

        push!(rindex, new_idx_entry())
        write!(rindex)

        @test length(index) == 1
        reload!(index)
        @test length(index) == 2

        entry = index["new_path", 3]
        @test entry.mode == 33188
    finally
        close(test_repo)
        LibGit2.free!(test_repo)
        Base.gc()
        rm(tmp_path, recursive=true)
    end
end

with_tmp_repo_access("test idempotent read write") do test_repo, path
    head_id = lookup_ref(test_repo, "HEAD") |> resolve |> target
    tree = GitTree(lookup(test_repo, head_id))
    index = GitIndex(test_repo)
    read_tree!(index, tree)

    index_tree_id = write_tree!(index)
    index_tree = test_repo[index_tree_id]
    @test Oid(index_tree) == Oid(tree)
end

with_tmp_repo_access("test build tree from index") do test_repo, path
    head_id = lookup_ref(test_repo,
                "refs/remotes/origin/packed") |> resolve |> target
    tree = GitTree(test_repo[head_id])
    index = GitIndex(test_repo)
    read_tree!(index, tree)
    remove!(index, "second.txt")

    new_tree_id = write_tree!(index)
    @test head_id != new_tree_id
    @test test_repo[new_tree_id]["second.txt"] == nothing
end

# --------------------
# test index all
# --------------------

# test add all lifecycle
sandboxed_test("testrepo.git") do _, path
    test_repo = init_repo(joinpath(path, "add-all"))
    cd(workdir(test_repo)) do
        open("file.foo", "w") do fh
            write(fh, "a file")
        end
        open("file.bar", "w") do fh
            write(fh, "another file")
        end
        open("file.zzz", "w") do fh
            write(fh, "yet another one")
        end
        open("other.zzz", "w") do fh
            write(fh, "yet another one")
        end
        open("more.zzz", "w") do fh
            write(fh, "yet another one")
        end
        open(".gitignore", "w") do fh
            write(fh, "*.foo\n")
        end

        index = GitIndex(test_repo)

        add_all!(index, "file.*")
        @test index["file.foo"]  == nothing
        @test index["file.bar"]  != nothing
        @test index["file.zzz"]  != nothing
        @test index["other.zzz"] == nothing
        @test index["more.zzz"]  == nothing

        add_all!(index, "*.zzz")
        @test index["file.foo"]  == nothing
        @test index["file.bar"]  != nothing
        @test index["file.zzz"]  != nothing
        @test index["other.zzz"] != nothing
        @test index["more.zzz"]  != nothing
    end
    close(test_repo)
    LibGit2.free!(test_repo)
end

# test update all
sandboxed_test("testrepo.git") do _, path
    test_repo = init_repo(joinpath(path, "add-all"))
    cd(workdir(test_repo)) do
        open("file.foo", "w") do fh
            write(fh, "a file")
        end
        open("file.bar", "w") do fh
            write(fh, "another file")
        end
        open("file.zzz", "w") do fh
            write(fh, "yet another one")
        end
        open("other.zzz", "w") do fh
            write(fh, "yet another one")
        end
        open("more.zzz", "w") do fh
            write(fh, "yet another one")
        end
        open(".gitignore", "w") do fh
            write(fh, "*.foo\n")
        end

        index = GitIndex(test_repo)
        add_all!(index, "file.*")

        open("file.bar", "w") do fh
            write(fh, "new content for file")
        end
        update_all!(index, "file.*")

        @test index["file.bar"] != nothing
        @test test_repo[Oid(index["file.bar"])] |> bytestring == "new content for file"

        @test index["other.zzz"] == nothing
        @test index["more.zzz"]  == nothing

        rm("file.bar")
        update_all!(index)

        @test index["file.bar"] == nothing
    end
    close(test_repo)
    LibGit2.free!(test_repo)
end

# test remove all
sandboxed_test("testrepo.git") do _, path
    test_repo = init_repo(joinpath(path, "add-all"))
    cd(workdir(test_repo)) do
        open("file.foo", "w") do fh
            write(fh, "a file")
        end
        open("file.bar", "w") do fh
            write(fh, "another file")
        end
        open("file.zzz", "w") do fh
            write(fh, "yet another one")
        end
        open("other.zzz", "w") do fh
            write(fh, "yet another one")
        end
        open("more.zzz", "w") do fh
            write(fh, "yet another one")
        end
        open(".gitignore", "w") do fh
            write(fh, "*.foo\n")
        end

        index = GitIndex(test_repo)

        add_all!(index, "file.*")
        remove_all!(index, "*.zzz")

        @test index["file.bar"] != nothing
        @test index["file.zzz"] == nothing
    end
    close(test_repo)
    LibGit2.free!(test_repo)
end


text:
@testset "create_surfaces_cylinder.jl" begin
    # Test create cylinders
    setup = CylinderSetup(
        name = "test",
        ncell = 10,
        rmin = 1.0,
        rmax = 5.0,
        dmin = 0.2,
        dmax = 0.3,
        bend = 0.0,
        twist = 0.0,
        height = 3.0,
        include_in = true,
        in_ratio = 0.6,
        ecs_shape = :convex_hull,
        ecs_ratio = 0.5,
        refinement = 0.5,
        D_in = 0.002 * I(3),
        D_out = 0.002 * I(3),
        D_ecs = 0.002 * I(3),
        T₂_in = Inf,
        T₂_out = Inf,
        T₂_ecs = Inf,
        ρ_in = 1.0,
        ρ_out = 1.0,
        ρ_ecs = 1.0,
        κ_in_out = 1e-5,
        κ_out_ecs = 1e-5,
        κ_in = 0.0,
        κ_out = 0.0,
        κ_ecs = 0.0,
    )

    cells = create_cells(setup)
    surfaces = create_surfaces_cylinder(cells, setup)

    ncompartment = 2 * setup.ncell + 1

    @test size(surfaces.points)[1] == 3
    @test size(surfaces.facets)[1] == 3
    @test length(surfaces.facetmarkers) == size(surfaces.facets)[2]
    @test size(surfaces.regions) == (3, ncompartment)
end


text:
#=
multivariate_filter.jl

    Multivariate Kalman filtering routines for a linear Gaussian State Space
	model

@author: Quint Wiersma <q.wiersma@vu.nl>

@date: 2020/07/28
=#

# Structs
struct MultivariateFilter{Ta, TP, Tv, TF, TK, Tm} <: KalmanFilter
	a::Ta	    # filtered state
	P::TP	    # filtered state variance
	v::Tv  	    # forecast error
	F::TF	    # forecast error variance
	K::TK	    # Kalman gain
    tmp_pn::Tm  # buffer
    tmp_p::Tm   # buffer
    tmp_n::Tm   # buffer
end
# Constructor
function MultivariateFilter(n::Integer, p::Integer, T_len::Integer, T::Type)
    # filter output
    a= Matrix{T}(undef, p, T_len)
    P= Array{T,3}(undef, p, p, T_len)
    v= Matrix{T}(undef, n, T_len)
    F= Array{T,3}(undef, n, n, T_len)
    K= Array{T,3}(undef, p, n, T_len)

    # buffers
    tmp_pn= Matrix{T}(undef, p, n)
    tmp_p= Matrix{T}(undef, p, p)
    tmp_n= Matrix{T}(undef, n, n)

    return MultivariateFilter(a, P, v, F, K, tmp_pn, tmp_p, tmp_n)
end

struct WoodburyFilter{Ta, TP, Tv, TFi, TK, Tm} <: KalmanFilter
	a::Ta	    # filtered state
	P::TP	    # filtered state variance
	v::Tv	    # forecast error
	Fi::TFi	    # inverse forecast error variance
	K::TK	    # Kalman gain
    tmp_np::Tm  # buffer
    tmp_pn::Tm  # buffer
    tmp_p::Tm   # buffer
end
# Constructor
function WoodburyFilter(n::Integer, p::Integer, T_len::Integer, T::Type)
    # filter output
    a= Matrix{T}(undef, p, T_len)
    P= Array{T,3}(undef, p, p, T_len)
    v= Matrix{T}(undef, n, T_len)
    Fi= Array{T,3}(undef, n, n, T_len)
    K= Array{T,3}(undef, p, n, T_len)

    # buffers
    tmp_np= Matrix{T}(undef, n, p)
    tmp_pn= Matrix{T}(undef, p, n)
    tmp_p= Matrix{T}(undef, p, p)

    return WoodburyFilter(a, P, v, Fi, K, tmp_np, tmp_pn, tmp_p)
end

"""
	error!(v_t, y, Z, a_t, d)

Compute forecast error ``v`` at time ``t``, storing the result in `v_t`.

#### Arguments
  - `y::AbstractVector`		: data (n x 1)
  - `Z::AbstractMatrix`		: system matrix ``Z`` (n x p)
  - `a_t::AbstractVector`	: predicted state (p x 1)
  - `d::AbstractVector`		: mean adjustment (n x 1)

#### Returns
  - `v_t::AbstractVector`	: forecast error (n x 1)
"""
function error!(v_t::AbstractVector, y::AbstractVector, Z::AbstractMatrix,
				a_t::AbstractVector, d::AbstractVector)
	# yₜ - d
	v_t.= y .- d
    # vₜ = yₜ - Z×aₜ - d
    mul!(v_t, Z, a_t, -1., 1.)

	return nothing
end

"""
	error_var!(F_t, P_t, Z, H, tmp)

Compute forecast error variance ``F`` at time ``t``, storing the result in
`F_t`.

#### Arguments
  - `P_t::AbstractMatrix`	: predicted state variance (p x p)
  - `Z::AbstractMatrix`		: system matrix ``Z`` (n x p)
  - `H::AbstractMatrix`		: system matrix ``H`` (n x n)
  - `tmp::AbstractMatrix`	: tmp storage array (p x n)

#### Returns
  - `F_t::AbstractMatrix`	: forecast error variance (n x n)
"""
function error_var!(F_t::AbstractMatrix, P_t::AbstractMatrix, Z::AbstractMatrix,
					H::AbstractMatrix, tmp::AbstractMatrix)
    # PₜｘZ'
	mul!(tmp, P_t, transpose(Z))
    # ZｘPₜｘZ'
	mul!(F_t, Z, tmp)
	# Fₜ = ZｘPₜｘZ' + H
	F_t.+= H

	return nothing
end

function error_var!(F_t::AbstractMatrix, P_t::AbstractMatrix, Z::AbstractMatrix,
					H::Diagonal, tmp::AbstractMatrix)
    # PₜｘZ'
	mul!(tmp, P_t, transpose(Z))
    # ZｘPₜｘZ'
	mul!(F_t, Z, tmp)
	# Fₜ = ZｘPₜｘZ' + H
	@inbounds @fastmath for i in axes(H,1)
		F_t[i,i]+= H.diag[i]
	end

	return nothing
end

"""
	error_prec!(Fi_t, P_t, Z, H, tmp)

Compute forecast error precision ``F⁻¹`` at time ``t`` using Woodbury's
identity, storing the result in `Fi_t`.

#### Arguments
  - `P_t::AbstractMatrix`	: predicted state variance (p x p)
  - `Z::AbstractMatrix`		: system matrix ``Z`` (n x p)
  - `Hi::AbstractMatrix`	: inverse system matrix ``H`` (n x n)
  - `tmp_np::AbstractMatrix`: tmp storage array (n x p)
  - `tmp_pn::AbstractMatrix`: tmp storage array (p x n)
  - `tmp_p::AbstractMatrix`	: tmp storage array (p x p)

#### Returns
  - `Fi_t::AbstractMatrix`	: forecast error precision (n x n)
"""
function error_prec!(Fi_t::AbstractMatrix, P_t::AbstractMatrix, Z::AbstractMatrix,
					Hi::AbstractMatrix, tmp_np::AbstractMatrix, tmp_pn::AbstractMatrix,
					tmp_p::AbstractMatrix)
	# Pₜ⁻¹
	copyto!(tmp_p, P_t)
	LinearAlgebra.inv!(cholesky!(Hermitian(tmp_p)))
	# H⁻¹×Z
	mul!(tmp_np, Hi, Z)
	# Pₜ⁻¹ + Z'×H⁻¹×Z
	mul!(tmp_p, transpose(Z), tmp_np, 1., 1.)
	# (Pₜ⁻¹ + Z'×H⁻¹×Z)⁻¹
	LinearAlgebra.inv!(cholesky!(Hermitian(tmp_p)))
	# (Pₜ⁻¹ + Z′×H⁻¹×Z)⁻¹×Z′×H⁻¹
	mul!(tmp_pn, tmp_p, transpose(tmp_np))
	# -H⁻¹×Z×(Pₜ⁻¹ + Z'×H⁻¹×Z)⁻¹×Z'×H⁻¹
	mul!(Fi_t, tmp_np, tmp_pn, -1., .0)
	# Fₜ⁻¹ = H⁻¹ - H⁻¹×Z×(Pₜ⁻¹ + Z'×H⁻¹×Z)⁻¹×Z'×H⁻¹
	Fi_t.+= Hi

	return nothing
end

function error_prec!(Fi_t::AbstractMatrix, P_t::AbstractMatrix, Z::AbstractMatrix,
					Hi::Diagonal, tmp_np::AbstractMatrix, tmp_pn::AbstractMatrix,
					tmp_p::AbstractMatrix)
	# Pₜ⁻¹
	copyto!(tmp_p, P_t)
	C= cholesky!(Hermitian(tmp_p))		# pointer(C.factors) = pointer(tmp_p)
	LinearAlgebra.inv!(C)
	# H⁻¹×Z
	mul!(tmp_np, Hi, Z)
	# Pₜ⁻¹ + Z'×H⁻¹×Z
	mul!(Pi, transpose(Z), tmp_np, 1., 1.)
	# (Pₜ⁻¹ + Z'×H⁻¹×Z)⁻¹
	C= cholesky!(Hermitian(tmp_p))		# pointer(C.factors) = pointer(tmp_p)
	LinearAlgebra.inv!(C)
	# (Pₜ⁻¹ + Z′×H⁻¹×Z)⁻¹×Z′×H⁻¹
	mul!(tmp_pn, tmp_p, transpose(tmp_np))
	# -H⁻¹×Z×(Pₜ⁻¹ + Z'×H⁻¹×Z)⁻¹×Z'×H⁻¹
	mul!(Fi_t, tmp_np, tmp_pn, -1., .0)
	# Fₜ⁻¹ = H⁻¹ - H⁻¹×Z×(Pₜ⁻¹ + Z'×H⁻¹×Z)⁻¹×Z'×H⁻¹
	@inbounds @fastmath for i in axes(Hi,1)
		Fi_t[i,i]+= Hi.diag[i]
	end

	return nothing
end

"""
	gain!(K_t, P_t, Fi_t, Z, T, tmp)

Compute Kalman gain ``K`` at time ``t`` using the inverse of ``F``, storing
the result in `K_t`.

#### Arguments
  - `P_t::AbstractMatrix`	: predicted state variance (p x p)
  - `Fi_t::AbstractMatrix`	: inverse of ``Fₜ`` (n x n)
  - `Z::AbstractMatrix`		: system matrix ``Z`` (n x p)
  - `T::AbstractMatrix`		: system matrix ``T`` (p x p)
  - `tmp::AbstractMatrix`	: tmp storage array (p x n)

#### Returns
  - `K_t::AbstractMatrix`	: Kalman gain (p x n)
"""
function gain!(K_t::AbstractMatrix, P_t::AbstractMatrix, Fi_t::AbstractMatrix,
				Z::AbstractMatrix, T::AbstractMatrix, tmp::AbstractMatrix)
    # PₜｘZ'
	mul!(K_t, P_t, transpose(Z))
    # PₜｘZ'ｘFₜ⁻¹
	mul!(tmp, K_t, Fi_t)
	# Kₜ = TｘPₜｘZ'ｘFₜ⁻¹
	mul!(K_t, T, tmp)

	return nothing
end

"""
	gain!(K_t, P_t, fac, Z, T, tmp)

Compute Kalman gain ``K`` at time ``t`` using the factorization of ``F``, storing
the result in `K_t`.
"""
function gain!(K_t::AbstractMatrix, P_t::AbstractMatrix, fac::Factorization,
				Z::AbstractMatrix, T::AbstractMatrix, tmp::AbstractMatrix)
    # PₜｘZ'
	mul!(tmp, P_t, transpose(Z))
    # PₜｘZ'ｘFₜ⁻¹
	rdiv!(tmp, fac)
	# Kₜ = TｘPₜｘZ'ｘFₜ⁻¹
	mul!(K_t, T, tmp)

	return nothing
end

"""
	predict_state!(a_p, a_t, K_t, v_t, T, c)

Predict states ``a`` at time ``t+1``, storing the result in `a_p`.

#### Arguments
  - `a_t::AbstractVector`	: predicted state at time ``t`` (p x 1)
  - `K_t::AbstractMatrix`	: Kalman gain (p x n)
  - `v_t::AbstractVector`	: forecast error (n x 1)
  - `T::AbstractMatrix`		: system matrix ``T`` (p x p)
  - `c::AbstractVector`		: mean adjustment (p x 1)

#### Returns
  - `a_p::AbstractVector`	: predicted state at time ``t+1`` (p x 1)
"""
function predict_state!(a_p::AbstractVector, a_t::AbstractVector, K_t::AbstractMatrix,
						v_t::AbstractVector, T::AbstractMatrix, c::AbstractVector)
    # Kₜｘvₜ
    mul!(a_p, K_t, v_t)
    # aₜ₊₁ = Tｘaₜ + Kₜｘvₜ
	mul!(a_p, T, a_t, 1., 1.)
	# aₜ₊₁ = Tｘaₜ + Kₜｘvₜ + c
	a_p.+= c

	return nothing
end

"""
	predict_state_var!(P_p, P_t, K_t, Z, T, Q, tmp)

Predict states variance ``P`` at time ``t+1``, storing the result in `P_p`.

#### Arguments
  - `P_t::AbstractMatrix`	: predicted states variance at time ``t`` (p x p)
  - `K_t::AbstractMatrix`	: Kalman gain (p x n)
  - `Z::AbstractMatrix`		: system matrix ``Z`` (n x p)
  - `T::AbstractMatrix`		: system matrix ``T`` (p x p)
  - `Q::AbstractMatrix`		: system matrix ``Q`` (p x p)
  - `tmp::AbstractMatrix`	: tmp storage array (p x p)

#### Returns
  - `P_p::AbstractMatrix`	: predicted states variance at time ``t+1`` (p x p)
"""
function predict_state_var!(P_p::AbstractMatrix, P_t::AbstractMatrix, K_t::AbstractMatrix,
						Z::AbstractMatrix, T::AbstractMatrix, Q::AbstractMatrix,
						tmp::AbstractMatrix)
    # -KₜｘZ
	mul!(P_p, K_t, Z, -1., .0)
	# T - KₜｘZ
	P_p.+= T
    # Pₜｘ(T - KₜｘZ)'
	mul!(tmp, P_t, transpose(P_p))
	# TｘPₜｘ(T - KₜｘZ)'
	mul!(P_p, T, tmp)
	# TｘPₜｘ(T - KₜｘZ)' + Q
	P_p.+= Q

	return nothing
end

function predict_state_var!(P_p::AbstractMatrix, P_t::AbstractMatrix, K_t::AbstractMatrix,
						Z::AbstractMatrix, T::Diagonal, Q::AbstractMatrix,
						tmp::AbstractMatrix)
    # -KₜｘZ
	mul!(P_p, K_t, Z, -1., .0)
	# T - KₜｘZ
	@inbounds @fastmath for i in axes(T,1)
		P_p[i,i]+= T.diag[i]
	end
    # Pₜｘ(T - KₜｘZ)'
	mul!(tmp, P_t, transpose(P_p))
	# TｘPₜｘ(T - KₜｘZ)'
	mul!(P_p, T, tmp)
	# TｘPₜｘ(T - KₜｘZ)' + Q
	P_p.+= Q

	return nothing
end

function predict_state_var!(P_p::AbstractMatrix, P_t::AbstractMatrix, K_t::AbstractMatrix,
						Z::AbstractMatrix, T::AbstractMatrix, Q::Diagonal,
						tmp::AbstractMatrix)
    # -KₜｘZ
	mul!(P_p, K_t, Z, -1., .0)
	# T - KₜｘZ
	P_p.+= T
    # Pₜｘ(T - KₜｘZ)'
	mul!(tmp, P_t, transpose(P_p))
	# TｘPₜｘ(T - KₜｘZ)'
	mul!(P_p, T, tmp)
	# TｘPₜｘ(T - KₜｘZ)' + Q
	@inbounds @fastmath for i in axes(Q,1)
		P_p[i,i]+= Q.diag[i]
	end

	return nothing
end

function predict_state_var!(P_p::AbstractMatrix, P_t::AbstractMatrix, K_t::AbstractMatrix,
						Z::AbstractMatrix, T::Diagonal, Q::Diagonal,
						tmp::AbstractMatrix)
    # -KₜｘZ
	mul!(P_p, K_t, Z, -1., .0)
	# T - KₜｘZ
	@inbounds @fastmath for i in axes(T,1)
		P_p[i,i]+= T.diag[i]
	end
    # Pₜｘ(T - KₜｘZ)'
	mul!(tmp, P_t, transpose(P_p))
	# TｘPₜｘ(T - KₜｘZ)'
	mul!(P_p, T, tmp)
	# TｘPₜｘ(T - KₜｘZ)' + Q
	@inbounds @fastmath for i in axes(Q,1)
		P_p[i,i]+= Q.diag[i]
	end

	return nothing
end

"""
	predict_state_var!(P_p, P_t, T, Q, tmp)

Predict states variance ``P`` at time ``t+1`` when observation ``y`` at time
``t`` is missing, i.e. NaN, storing the result in `P_p`.

#### Arguments
  - `P_t::AbstractMatrix`	: predicted states variance at time ``t`` (p x p)
  - `T::AbstractMatrix`		: system matrix ``T`` (p x p)
  - `Q::AbstractMatrix`		: system matrix ``Q`` (p x p)
  - `tmp::AbstractMatrix`	: tmp storage array (p x p)

#### Returns
  - `P_p::AbstractMatrix`	: predicted states variance at time ``t+1`` (p x p)
"""
function predict_state_var!(P_p::AbstractMatrix, P_t::AbstractMatrix,
                            T::AbstractMatrix, Q::AbstractMatrix,
						    tmp::AbstractMatrix)
    # PₜｘT'
    mul!(tmp, P_t, transpose(T))
	# TｘPₜｘT'
	mul!(P_p, T, tmp)
	# TｘPₜｘT' + Q
	P_p.+= Q

	return nothing
end

function predict_state_var!(P_p::AbstractMatrix, P_t::AbstractMatrix, T::Diagonal,
                            Q::AbstractMatrix, tmp::AbstractMatrix)
	# TｘPₜｘT' + Q
	P_p.= T.diag .* P_t .* transpose(T.diag) .+ Q

	return nothing
end

function predict_state_var!(P_p::AbstractMatrix, P_t::AbstractMatrix,
                            T::AbstractMatrix, Q::Diagonal, tmp::AbstractMatrix)
    # PₜｘT'
    mul!(tmp, P_t, transpose(T))
	# TｘPₜｘT'
	mul!(P_p, T, tmp)
	# TｘPₜｘT' + Q
	@inbounds @fastmath for i in axes(Q,1)
		P_p[i,i]+= Q.diag[i]
	end

	return nothing
end

function predict_state_var!(P_p::AbstractMatrix, P_t::AbstractMatrix, T::Diagonal,
                            Q::Diagonal, tmp::AbstractMatrix)
    # TｘPₜｘT'
	P_p.= T.diag .* P_t .* transpose(T.diag)
    # TｘPₜｘT' + Q
    @inbounds @fastmath for i in axes(Q,1)
		P_p[i,i]+= Q.diag[i]
	end

	return nothing
end

"""
    update_filter!(filter, y, Z, T, d, c, H, Q, t, T_len)

Update Kalman filter components at time `t`, storing the results in `filter`.

#### Arguments
  - `y::AbstractVector`     : data
  - `Z::AbstractMatrix`     : system matrix ``Z``
  - `T::AbstractMatrix`     : system matrix ``T``
  - `d::AbstractVector`     : mean adjustment observation
  - `c::AbstractVector`     : mean adjustment state
  - `H::AbstractMatrix`     : system matrix ``H``
  - `Q::AbstractMatrix`     : system matrix ``Q``
  - `t::Integer`            : time
  - `T_len::Integer`        : number of obs

#### Returns
  - `filter::KalmanFilter`  : Kalman filter components
"""
function update_filter!(filter::MultivariateFilter,
                        y::AbstractVector,
                        Z::AbstractMatrix, T::AbstractMatrix,
                        d::AbstractVector, c::AbstractVector,
                        H::AbstractMatrix, Q::AbstractMatrix,
                        t::Integer, T_len::Integer)
    # Store views
    a= view(filter.a,:,t)
    P= view(filter.P,:,:,t)
    v= view(filter.v,:,t)
    F= view(filter.F,:,:,t)
    K= view(filter.K,:,:,t)

    # Forecast error
    error!(v, y, Z, a, d)

    # Forecast error variance
    error_var!(F, P, Z, H, filter.tmp_pn)

    # Cholesky factorization of Fₜ
    copyto!(filter.tmp_n, F)
    fac= cholesky!(Hermitian(filter.tmp_n))

    # Kalman gain
    gain!(K, P, fac, Z, T, filter.tmp_pn)

    if t < T_len
        # Store views
        a_p= view(filter.a,:,t+1)
        P_p= view(filter.P,:,:,t+1)

        # Predict states
        predict_state!(a_p, a, K, v, T, c)

        # Predict states variance
        predict_state_var!(P_p, P, K, Z, T, Q, filter.tmp_p)
    end

    return nothing
end

function update_filter!(filter::WoodburyFilter,
                        y::AbstractVector,
                        Z::AbstractMatrix, T::AbstractMatrix,
                        d::AbstractVector, c::AbstractVector,
                        Hi::AbstractMatrix, Q::AbstractMatrix,
                        t::Integer, T_len::Integer)
    # Store views
    a= view(filter.a,:,t)
    P= view(filter.P,:,:,t)
    v= view(filter.v,:,t)
    Fi= view(filter.Fi,:,:,t)
    K= view(filter.K,:,:,t)

    # Forecast error
    error!(v, y, Z, a, d)

    # Forecast error precision
    error_prec!(Fi, P, Z, Hi, filter.tmp_np, filter.tmp_pn, filter.tmp_p)

    # Kalman gain
    gain!(K, P, Fi, Z, T, filter.tmp_pn)

    if t < T_len
        # Store views
        a_p= view(filter.a,:,t+1)
        P_p= view(filter.P,:,:,t+1)

        # Predict states
        predict_state!(a_p, a, K, v, T, c)

        # Predict states variance
        predict_state_var!(P_p, P, K, Z, T, Q, filter.tmp_p)
    end

    return nothing
end

"""
    update_filter!(filter, T, d, c, H, Q, t, T_len)

Update Kalman filter components at time `t` when observation ``y_t`` is missing,
i.e. NaN, storing the results in `filter`.

#### Arguments
  - `T::AbstractMatrix`     : system matrix ``T``
  - `d::AbstractVector`     : mean adjustment observation
  - `c::AbstractVector`     : mean adjustment state
  - `H::AbstractMatrix`     : system matrix ``H``
  - `Q::AbstractMatrix`     : system matrix ``Q``
  - `t::Integer`            : time
  - `T_len::Integer`        : number of obs

#### Returns
  - `filter::KalmanFilter`  : Kalman filter components
"""
function update_filter!(filter::MultivariateFilter,
                        Z::AbstractMatrix, T::AbstractMatrix,
                        c::AbstractVector,
                        H::AbstractMatrix, Q::AbstractMatrix,
                        t::Integer, T_len::Integer)
    # Store views
    a= view(filter.a,:,t)
    P= view(filter.P,:,:,t)
    v= view(filter.v,:,t)
    F= view(filter.F,:,:,t)
    K= view(filter.K,:,:,t)

    # Forecast error
    v.= NaN

    # Forecast error variance
    error_var!(F, P, Z, H, filter.tmp_pn)

    # Kalman gain
    K.= zero(eltype(K))

    if t < T_len
        # Store views
        a_p= view(filter.a,:,t+1)
        P_p= view(filter.P,:,:,t+1)

        # Predict states
        a_p.= c
        mul!(a_p, T, a, 1., 1.)

        # Predict states variance
        predict_state_var!(P_p, P, T, Q, filter.tmp_p)
    end

    return nothing
end

function update_filter!(filter::WoodburyFilter,
                        Z::AbstractMatrix, T::AbstractMatrix,
                        c::AbstractVector,
                        Hi::AbstractMatrix, Q::AbstractMatrix,
                        t::Integer, T_len::Integer)
    # Store views
    a= view(filter.a,:,t)
    P= view(filter.P,:,:,t)
    v= view(filter.v,:,t)
    Fi= view(filter.F,:,:,t)
    K= view(filter.K,:,:,t)

    # Forecast error
    v.= NaN

    # Forecast error precision
    error_prec!(Fi, P, Z, Hi, filter.tmp_np, filter.tmp_pn, filter.tmp_p)

    # Kalman gain
    K.= zero(eltype(K))

    if t < T_len
        # Store views
        a_p= view(filter.a,:,t+1)
        P_p= view(filter.P,:,:,t+1)

        # Predict states
        a_p.= c
        mul!(a_p, T, a, 1., 1.)

        # Predict states variance
        predict_state_var!(P_p, P, T, Q, filter.tmp_p)
    end

    return nothing
end

"""
	kalman_filter!(filter, sys)

Compute predicted states ``a`` and forecast errors ``v`` with corresponding
variances ``P`` and ``F`` and Kalman gain ``K`` for a linear Gaussian State
Space model with system matrices `sys` using the Kalman filter, storing the
results in `filter`.

#### Arguments
  - `sys::StateSpaceSystem`	: state space system matrices

#### Returns
  - `filter::MultivariateFilter`: Kalman filter output
"""
function kalman_filter!(filter::MultivariateFilter, sys::LinearTimeInvariant)
	# Get dims
	T_len= size(sys.y,2)

	# Initialize filter
	filter.a[:,1]= sys.a1
	filter.P[:,:,1]= sys.P1

	# Filter
	@inbounds for t in 1:T_len
        # data
        y_t= view(sys.y, :, t)

        # update filter
        if any(isnan, y_t)
            update_filter!(filter, sys.Z, sys.T, sys.c, sys.H, sys.Q, t, T_len)
        else
            update_filter!(filter, y_t, sys.Z, sys.T, sys.d, sys.c, sys.H, sys.Q, t, T_len)
        end
	end

	return nothing
end

function kalman_filter!(filter::MultivariateFilter, sys::LinearTimeVariant)
	# Get dims
	T_len= size(sys.y,2)

	# Initialize filter
	filter.a[:,1]= sys.a1
	filter.P[:,:,1]= sys.P1

	# Filter
	@inbounds for t in 1:T_len
        # data
        y_t= view(sys.y, :, t)

        # means
        d_t= view(sys.d, :, t)
        c_t= view(sys.c, :, t)

        # update filter
        if any(isnan, y_t)
            update_filter!(filter, sys.Z[t], sys.T[t], c_t, sys.H[t], sys.Q[t], t, T_len)
        else
            update_filter!(filter, y_t, sys.Z[t], sys.T[t], d_t, c_t, sys.H[t], sys.Q[t], t, T_len)
        end
	end

	return nothing
end

"""
	kalman_filter!(filter, sys)

Compute predicted states ``a`` and forecast errors ``v`` with corresponding
variance ``P`` and precision ``F⁻¹`` and Kalman gain ``K`` for a linear Gaussian
State Space model with system matrices `sys` using the Kalman filter based on
Woodbury's Identity, storing the results in `filter`.

Woodbury's Identity allows direct computation of the inverse variance
(precision) ``F⁻¹``.

#### Arguments
  - `sys::StateSpaceSystem`	: state space system matrices

#### Returns
  - `filter::WoodburyFilter`: Kalman filter output
"""
function kalman_filter!(filter::WoodburyFilter, sys::LinearTimeInvariant)
	# Get dims
	T_len= size(sys.y,2)

	# Inverse of H
	Hi= copy(sys.H)
	if sys.H isa Diagonal
		Hi.diag.= inv.(sys.H.diag)
	else
		LinearAlgebra.inv!(cholesky!(Hi))
	end

	# Initialize filter
	filter.a[:,1]= sys.a1
	filter.P[:,:,1]= sys.P1

	# Filter
	@inbounds for t in 1:T_len
        # data
        y_t= view(sys.y, :, t)

        # update filter
        if any(isnan, y_t)
            update_filter!(filter, sys.Z, sys.T, sys.c, Hi, sys.Q, t, T_len)
        else
            update_filter!(filter, y_t, sys.Z, sys.T, sys.d, sys.c, Hi, sys.Q, t, T_len)
        end
	end

	return nothing
end

function kalman_filter!(filter::WoodburyFilter, sys::LinearTimeVariant)
	# Get dims
	T_len= size(sys.y,2)

    # Initialize temp. containers
	Hi= similar(sys.H[1])

	# Initialize filter
	filter.a[:,1]= sys.a1
	filter.P[:,:,1]= sys.P1

	# Filter
	@inbounds for t in 1:T_len
        # data
        y_t= view(sys.y, :, t)

        # means
        d_t= view(sys.d, :, t)
        c_t= view(sys.c, :, t)

        # Inverse of H
        Hi.= sys.H[t]
        if sys.H[t] isa Diagonal
            Hi.diag.= inv.(sys.H[t].diag)
        else
            LinearAlgebra.inv!(cholesky!(Hi))
        end

        # update filter
        if any(isnan, y_t)
            update_filter!(filter, sys.Z[t], sys.T[t], c_t, Hi, sys.Q[t], t, T_len)
        else
            update_filter!(filter, y_t, sys.Z[t], sys.T[t], d_t, c_t, Hi, sys.Q[t], t, T_len)
        end
	end

	return nothing
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    g1 ∩ g2

Return the intersection of two geometries `g1` and `g2`.
"""
Base.intersect(g1::Geometry, g2::Geometry) = get(intersecttype(g1, g2))

"""
    intersecttype(g1, g2)

Return the intersection type of two geometries `g1` and `g2`.
"""
function intersecttype end

"""
    Intersection

An intersection type (e.g. crossing line segments, overlapping boxes).
"""
abstract type Intersection end

Base.get(I::Intersection) = I.value

# ------------------------------
# SEGMENT-SEGMENT INTERSECTIONS
# ------------------------------

struct CrossingSegments{P<:Point} <: Intersection
  value::P
end

struct MidTouchingSegments{P<:Point} <: Intersection
  value::P
end

struct CornerTouchingSegments{P<:Point} <: Intersection
  value::P
end

struct OverlappingSegments{S<:Segment} <: Intersection
  value::S
end

struct NonIntersectingSegments <: Intersection end

Base.get(::NonIntersectingSegments) = nothing

# ----------------------
# BOX-BOX INTERSECTIONS
# ----------------------

struct OverlappingBoxes{B<:Box} <: Intersection
  value::B
end

struct FaceTouchingBoxes{B<:Box} <: Intersection
  value::B
end

struct CornerTouchingBoxes{P<:Point} <: Intersection
  value::P
end

struct NonIntersectingBoxes <: Intersection end

Base.get(::NonIntersectingBoxes) = nothing

# ----------------
# IMPLEMENTATIONS
# ----------------

include("intersections/segments.jl")
include("intersections/boxes.jl")


text:
@testset "test shortest path" begin
    errcode, path = dubins_shortest_path(zeros(3), [1.0, 0.0, 0.0], 1.0)
    @test errcode == EDUBOK
end

@testset "test invalid ρ" begin
    errcode, path = dubins_shortest_path(zeros(3), [1.0, 0.0, 0.0], -1.0)
    @test errcode == EDUBBADRHO
    @test path == nothing
end

@testset "test no path" begin
    errcode, path = dubins_path(zeros(3), [10.0, 0.0, 0.0], 1.0, LRL)
    @test errcode == EDUBNOPATH
    @test path == nothing
end

@testset "test path length" begin
    errcode, path = dubins_shortest_path(zeros(3), [4.0, 0.0, 0.0], 1.0)
    @test errcode == EDUBOK
    path_length = dubins_path_length(path)
    @test isapprox(path_length, 4., atol=1e-3)
end

@testset "test simple path" begin
    errcode, path = dubins_path(zeros(3), [1.0, 0.0, 0.0], 1.0, LSL)
    @test errcode == EDUBOK
end

@testset "test segment lengths" begin
    errcode, path = dubins_path(zeros(3), [4.0, 0.0, 0.0], 1.0, LSL)
    @test errcode == EDUBOK
    @test dubins_segment_length_normalized(path, 0) == Inf
    @test dubins_segment_length_normalized(path, 1) == 0.0
    @test dubins_segment_length_normalized(path, 2) == 4.0
    @test dubins_segment_length_normalized(path, 3) == 0.0
    @test dubins_segment_length_normalized(path, 4) == Inf
end

@testset "test sample" begin
    errcode, path = dubins_path(zeros(3), [4.0, 0.0, 0.], 1.0, LSL)
    @test errcode == EDUBOK

    errcode, qsamp = dubins_path_sample(path, 0.0)
    @test errcode == EDUBOK
    @test qsamp == zeros(3)

    errcode, qsamp = dubins_path_sample(path, 4.0)
    @test errcode == EDUBOK
    @test qsamp == [4.0, 0.0, 0.0]
end

@testset "test sample out of bounds" begin
    errcode, path = dubins_path(zeros(3), [4.0, 0.0, 0.0], 1.0, LSL)
    @test errcode == EDUBOK

    errcode, qsamp = dubins_path_sample(path, -1.0)
    @test errcode == EDUBPARAM
    @test qsamp == nothing

    errcode, qsamp = dubins_path_sample(path, 5.0)
    @test errcode == EDUBPARAM
    @test qsamp == nothing
end

@testset "test sample many LSL" begin
    errcode, path = dubins_path([0.0, 0.0, π/2], [4.0, 0.0, -π/2], 1.0, LSL)
    @test errcode == EDUBOK

    errcode, configurations = dubins_path_sample_many(path, 1.0)
    @test errcode == 0
    @test isapprox(configurations[1], [0.0, 0.0, π/2], atol=1e-8)
end

@testset "test sample many RLR" begin
    errcode, path = dubins_path([0.0, 0.0, π/2], [4.0, 0.0, -π/2], 1.0, RLR)
    @test errcode == EDUBOK

    errcode, configurations = dubins_path_sample_many(path, 1.0)
    @test errcode == 0
    @test isapprox(configurations[1], [0.0, 0.0, π/2], atol=1e-8)
end

@testset "test path type" begin
    for i in 0:5
        path_type::DubinsPathType = (DubinsPathType)(i)
        errcode, path = dubins_path(zeros(3), [1.0, 0.0, 0.0], 1.0, path_type)
        (errcode == EDUBOK) && (@test dubins_path_type(path) == path_type)
    end
end

@testset "test end point" begin
    errcode, path = dubins_path(zeros(3), [4.0, 0.0, 0.0], 1.0, LSL)
    @test errcode == EDUBOK

    errcode, qsamp = dubins_path_endpoint(path)
    @test isapprox(qsamp, [4.0, 0.0, 0.0], atol=1e-8)
end

@testset "test extract sub-path" begin
    errcode, path = dubins_path(zeros(3), [4.0, 0.0, 0.0], 1.0, LSL)
    @test errcode == EDUBOK

    errcode, subpath = dubins_extract_subpath(path, 2.0)
    @test errcode == EDUBOK

    errcode, qsamp = dubins_path_endpoint(subpath)
    @test isapprox(qsamp, [2.0, 0.0, 0.0], atol=1e-8)
end

@testset "test extract invalid sub-path" begin
    errcode, path = dubins_path(zeros(3), [4.0, 0.0, 0.0], 1.0, LSL)
    @test errcode == EDUBOK

    errcode, subpath = dubins_extract_subpath(path, 8.0)
    @test errcode == EDUBPARAM
    @test subpath == nothing
end


text:
function potnucl!(ptnucl::Bool, nr, r, zn, vn )
    SMALL = 1e-10
    if abs(zn) <= SMALL
        vn[:] .= 0.0
        return
    end
    if ptnucl
        # nucleus is taken to be a point particle
        for ir in 1:nr
            vn[ir] = zn/r[ir]
        end
    else
        # approximate nuclear radius
        rn = radnucl(zn)
        t1 = zn/(2.0*rn^3)
        t2 = 3.0*rn^2
        for ir in 1:nr
          if r[ir] < rn
            vn[ir] = t1*(t2 - r[ir]^2)
          else
            vn[ir] = zn/r[ir]
          end
        end
    end
    return
end


text:
module Constructors

export classconstructors, ConstructorDescriptor

using JavaCall.JNI
using JavaCall.Conversions
using JavaCall.Core

using JavaCall.Reflection: Classes

#=
Struct to hold information about java constructors
used to generate constructors that call the jni
This is not a replacement for java.lang.reflect.Constructor
as it only stores essential information

paramtypes: List of ClassDescriptor with the parameter
            types accepted by the constructor
=#
struct ConstructorDescriptor
    paramtypes::Vector{Classes.ClassDescriptor}
end

function Base.show(io::IO, c::ConstructorDescriptor)
    print(io, "ConstructorDescriptor{paramtypes: ", c.paramtypes, "}");
end

function Base.:(==)(x::ConstructorDescriptor, y::ConstructorDescriptor)
    x.paramtypes == y.paramtypes
end

function descriptorfromconstructor(constructor::jobject)
    paramtypes = callinstancemethod(
        constructor,
        :getParameterTypes,
        Vector{Symbol("java.lang.Class")},
        [])
    ConstructorDescriptor(map(
        Classes.descriptorfromclass,
        convert_to_vector(Vector{jclass}, paramtypes)
    ))
end

classconstructors(classname::Symbol) = classconstructors(Classes.findclass(classname))

function classconstructors(classdescriptor::Classes.ClassDescriptor)
    array = convert_to_vector(Vector{jobject}, callinstancemethod(
        classdescriptor.jniclass,
        :getConstructors,
        Vector{Symbol("java.lang.reflect.Constructor")},
        []))
    map(descriptorfromconstructor, array)
end

end


text:
facts("[nl_write] Temp file handling") do
    # Turn on debug mode so files persist
    old_debug = AmplNLWriter.debug
    AmplNLWriter.setdebug(true)

    filename = "test"
    filepath = joinpath(AmplNLWriter.solverdata_dir, filename)
    probfile = "$filepath.nl"
    solfile = "$filepath.sol"
    AmplNLWriter.clean_solverdata()

    context("all temp files deleted successfully") do
        @fact length(readdir(AmplNLWriter.solverdata_dir)) --> 1
    end

    solver = BonminNLSolver(filename=filename)
    m = Model(solver=solver)
    @variable(m, x >= 0)
    @objective(m, Min, x)
    solve(m)

    context("temp files present after solve in debug mode") do
        @fact length(readdir(AmplNLWriter.solverdata_dir)) --> 3
    end
    context("temp files used custom name") do
        @fact isfile(probfile) --> true
        @fact isfile(solfile) --> true
    end

    # Reset debug mode and clean up
    AmplNLWriter.setdebug(old_debug)
    AmplNLWriter.clean_solverdata()
end


text:
# BioGenerics.jl
# ==============
#
# Core types and methods common to many packages in the BioJulia ecosystem.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/BioGenerics/blob/master/LICENSE.md

module BioGenerics

include("methods.jl")
include("Exceptions.jl")
include("Automa.jl")
include("IO.jl")
#include("Mem.jl")
#include("ReaderHelper.jl")
#include("RecordHelper.jl")
include("Testing.jl")

end # module BioGenerics

max_forks_repo_name: BioJulia/BioGenerics


text:
using Dates: Millisecond
using TimeToLive: TTL, Node
using Test

const p = Millisecond(250)

@testset "TimeToLive.jl" begin
    @testset "Constructors" begin
        t = TTL(p)
        @test t.ttl == p
        @test t.d isa Dict{Any, Node{Any}}
        t = TTL{Int, String}(p)
        @test t.d isa Dict{Int, Node{String}}
    end

    @testset "Basic expiry" begin
        t = TTL(p)
        t[0] = "!"
        @test get(t, 0, nothing) == "!"
        sleep(2p)
        @test isempty(t)
    end

    @testset "Refreshing expiry" begin
        t = TTL(p)
        t[0] = "!"
        sleep(p/2)
        touch(t, 0)
        sleep(p)
        @test get(t, 0, nothing) == "!"
        sleep(2p)
        @test isempty(t)
    end

    @testset "Iteration" begin
        t = TTL(p)
        t[1] = t[2] = t[3] = t[4] = t[5] = "!"
        count = 0
        for pair in t
            @test pair isa Pair{Int, String}
            count += 1
        end
        @test count == 5
    end

    @testset "Refresh on access" begin
        t = TTL(p; refresh_on_access=true)
        t[0] = "!"
        sleep(p/2)
        t[0]
        sleep(p)
        @test get(t, 0, nothing) == "!"
    end

    @testset "Disabled refresh on access" begin
        t = TTL(p; refresh_on_access=false)
        t[0] = "!"
        sleep(p/2)
        t[0]
        sleep(p)
        @test get(t, 0, nothing) === nothing
    end

    @testset "Troublesome Base methods" begin
        t = TTL{Int, Int}(p)
        t[1] = 2
        t[2] = 3
        t[3] = 4
        t[4] = 5
        t[5] = 6

        @test all(v -> v isa Int, values(t))
        pair = pop!(t)
        @test pair isa Pair{Int, Int}
        t[pair.first] = pair.second
        @test pop!(t, 4) == 5
        @test !haskey(t, 4)
        @test length(filter!(p -> p.second > 2, t)) == 3
        @test !haskey(t, 1)
    end
end


text:
using Distributions
using IntervalArithmetic
using Plots
using SlicedNormals

n = 500

# Step 1
θ = rand(Normal(π / 2, 1.3), n)
r = 3 .+ rand(Uniform(0, 0.2), n) .* (θ .- π / 2)

δ1 = r .* cos.(θ)
δ2 = r .* sin.(θ)

idx = δ1 .< 0
δ2[idx] = δ2[idx] * -1

δ = [δ1 δ2]

# Step 2
θ = rand(Normal(π / 2, 1.3), n)
r = 3 .+ rand(Uniform(0, 0.2), n) .* (θ .- π / 2)

δ1 = r .* cos.(θ)
δ2 = r .* sin.(θ)

idx = δ1 .< 0
δ2[idx] = δ2[idx] * -1

δ = vcat(δ, [δ1 δ2] .* -1)

# Fit Sliced Normal Distribution
d = 5
b = 1000

# Use baseline fit
sn, lh = fit_baseline(δ, d)

println("Baseline likelihood: $lh")

samples = rand(sn, 500)

p = scatter(
    δ[:, 1], δ[:, 2]; aspect_ratio=:equal, lims=[-4, 4], xlab="δ1", ylab="δ2", label="data"
)
scatter!(p, samples[:, 1], samples[:, 2]; label="samples")

display(p)

# Use baseline fit and then scale P
sn, lh = fit_scaling(δ, d)

println("Scaling likelihood: $lh")

samples = rand(sn, 500)

p = scatter(
    δ[:, 1], δ[:, 2]; aspect_ratio=:equal, lims=[-4, 4], xlab="δ1", ylab="δ2", label="data"
)
scatter!(p, samples[:, 1], samples[:, 2]; label="samples")

display(p)

# Plot density
xs = range(-4, 4; length=200)
ys = range(-4, 4; length=200)

contour!(xs, ys, (x, y) -> SlicedNormals.pdf(sn, [x, y], false))


text:
using LinearAlgebra: Matrix
using Random, LinearAlgebra

f(X) = 3sin(X[1]/3.0) + 2cos(X[2]/2.0)

τ(n) = 1/√(10n)

mutable struct CMAES{T}
    dim::Int
    λ::Int
    μ::Int

    centroid::Array{T, 1}
    c_m::T

    weights::Array{T, 1}
    μ_eff::T

    σ::T
    p_σ::Array{T, 1}
    c_σ::T
    d_σ::T

    C::Array{T, 2}
    p_c::Array{T, 1}
    c_c::T
    c_1::T
    c_μ::T
end

# A の各列を引数としてFの最も高い1~n位のインデックスを求める
function best_n(A::Array, n::Int64, F::Function)
    score = zeros(size(A)[2])
    for i in 1:size(A)[2]
        score[i] = F(A[:, i])
    end

    return sortperm(score, rev=true)[1:n]
end

function test_commaES(;μ::Int64, ρ::Int64, λ::Int64, F::Function, t=100)
    if λ > μ
        throw(DomainError(λ, "λ must be larger than μ"))
    end
    # 親群の初期化 (x, y, σ) x 10
    P = 2 * rand(Float64, (3, μ)) .- 1

    for n in 1:t
        # 子供を作る
        perm = best_n(P[1:2, :], ρ, F)
        Pρ = copy(P[:, perm])

        # 組み替え & 突然変異
        C = zeros(3, λ)
        ave_σ = sum(Pρ[3, :])/ρ
        C[3, :] = ave_σ * exp.(τ(n) .* randn(λ))

        ave_y = (sum(Pρ[1, :])./ρ, sum(Pρ[2, :])./ρ)
        for i in 1:λ
            C[1:2, i] = ave_y .+ C[3, i] .* randn(2)
        end

        # 次の親を選ぶ
        perm = best_n(C[1:2, :], μ, F)
        P = copy(C[:, perm])
    end

    return P
end

function test_plusES(;μ::Int64, ρ::Int64, λ::Int64, F::Function, t=100)
    # 親群の初期化 (x, y, σ) x 10
    P = 100 * randn(Float64, (3, μ))

    for n in 1:t
        # 子供を作る
        perm = best_n(P[1:2, :], ρ, F)
        Pρ = copy(P[:, perm])

        # 組み替え & 突然変異
        C = zeros(3, λ)
        ave_σ = sum(Pρ[3, :])/ρ
        C[3, :] = ave_σ * exp.(τ(n) .* randn(λ))

        ave_y = (sum(Pρ[1, :])./ρ, sum(Pρ[2, :])./ρ)
        for i in 1:λ
            C[1:2, i] = ave_y .+ C[3, i] .* randn(2)
        end

        # 次の親を選ぶ
        candidate = hcat(C, P)
        perm = best_n(candidate[1:2, :], μ, F)
        P = copy(candidate[:, perm])
    end

    return P
end

# size(X) == (λ, dim)
function get_fitness(X::Matrix{Float64}, f::Function)
    ToReturn = zeros(Float64, size(X)[2])

    for i in 1:size(X)[2]
        ToReturn[i] = f(X[:, i])
    end

    return ToReturn
end

# CMAESパラメータの初期化
function init_CMAES(center::Array{T, 1}, σ::T, λ::Int) where T <: Real
    ## 初期化フェーズ
    # 次元数
    dim = length(center)

    # 世代の総個体数 λ と エリート数 μ
    λ = λ > 0 ? λ : round(Int, 4 + 3 * log(dim))
    μ = λ÷2

    # 正規分布の中心と学習率
    centroid = T.(center)
    c_m = one(T)

    # 順位重み係数
    weights = [log(0.5(λ+1)) - log(i) for i in 1:μ]
    weights = T.(weights ./ sum(weights))
    μ_eff = 1 / sum(weights.^2)

    # ステップサイズ
    p_σ = zeros(T, dim)
    c_σ = (μ_eff + 2) / (dim + μ_eff + 5)
    d_σ = 1 + 2 * max(0, sqrt((μ_eff - 1)/(dim + 1)) - 1) + c_σ

    # 共分散行列 進化パス p と rank-μ, rank-one更新の学習率c
    C = Matrix{T}(I, dim, dim)
    p_c = zeros(T, dim)
    c_c = (4 + μ_eff / dim) / (dim + 4 + 2 * μ_eff/dim)
    c_1 = 2 / ((dim + 1.3f0)^2 + μ_eff)
    c_μ = min(1 - c_1, 2 * (μ_eff - 2 + 1/μ_eff) / ((dim+2)^2 + μ_eff))

    ToReturn = CMAES(dim, λ, μ, centroid, c_m, weights, μ_eff, σ, p_σ, c_σ, d_σ,
                C, p_c, c_c, c_1, c_μ)
    return ToReturn
end

# 個体の生成
function samplePopulation(self::CMAES{T}; rng=MersenneTwister(123)) where T <: Real
    ## 探索フェーズ
    # z ~ N(0, I) なる個体を λ 個生成
    Z = randn(rng, (self.dim, self.λ))

    # C を固有値分解
    Ei = eigen(self.C)
    diagD = sqrt(abs.(Diagonal(Ei.values)))
    B = Ei.vectors
    BD = B * diagD

    # y ~ N(0, C)
    Y = BD * Z
    # X ~ N(μ, σC)
    X = self.centroid .+ self.σ * Y

    return X
end

# 個体および行列 C の更新
function update!(self::CMAES{T}, X, fitnesses, gen) where T <: Real
    ### 更新フェーズ
    """ 正規分布パラメータの更新
    X : 個体群, shape == (λ, μ)
    fitnesses : 適合度
    gen : 現代の世代数
    """

    ## 1. Selection and recombination
    old_centroid = self.centroid
    old_σ = self.σ

    # fitnesses が上位 μ までのインデックスを抽出
    elite_indices = sortperm(fitnesses)[1:self.μ]

    X_elite = X[:, elite_indices]
    Y_elite = (X_elite .- old_centroid) ./ old_σ

    X_w = (X_elite * self.weights)
    Y_w = (Y_elite * self.weights)

    # 正規分布中心の更新
    self.centroid = (1 - self.c_m) * old_centroid .+ self.c_m * X_w

    ## 2. Step-size control
    Ei = eigen(self.C)
    diagD = sqrt(abs.(Diagonal(Ei.values)))
    B = Ei.vectors
    inv_diagD = inv(diagD)

    # Note. B*Z == C_ * Y
    C_ = B * inv_diagD * B'

    new_p_σ = (1 - self.c_σ) * self.p_σ
    new_p_σ += sqrt(self.c_σ * (2 - self.c_σ) * self.μ_eff) * C_ * Y_w
    self.p_σ = new_p_σ

    E_normal = sqrt(self.dim) * (1 - 1/(4*self.dim) + 1/(21 * self.dim ^ 2)) # 定数パラメータ
    self.σ = self.σ * exp((self.c_σ / self.d_σ) * (sqrt(sum(self.p_σ .^ 2))/E_normal - 1))

    if isinf(self.σ)
        self.σ = old_σ
    end

    # 3. Covariance matrix adaptation (CMA)
    # Note. h_σ(heaviside 関数) はステップサイズ σ が大きいときに C の更新を中断させるのに使う
    left = sqrt(sum(self.p_σ.^2)) / sqrt(1 - (1 - self.c_σ) ^ (2 * (gen + 1)))
    right = (1.4 + 2 / (self.dim + 1)) * E_normal
    hσ = left < right ? 1 : 0
    d_hσ = (1 - hσ) * self.c_c * (2 - self.c_c)

    # p_c の更新
    new_p_c = (1 - self.c_c) * self.p_c
    new_p_c += hσ * sqrt(self.c_c * (2 - self.c_c) * self.μ_eff) * Y_w
    self.p_c = new_p_c

    # C の更新
    new_C = (1 + self.c_1 * d_hσ - self.c_1 - self.c_μ) * self.C
    new_C += self.c_1 * [i * j for i in self.p_c, j in self.p_c]

    # 愚直な実装(スマートな実装はdeapのcma.pyを参照)
    wyy = zeros(T, (self.dim, self.dim))
    for i in 1:self.μ
        y_i = Y_elite[i]
        wyy .+= self.weights[i] * [y1 * y2 for y1 in y_i, y2 in y_i]
    end
    new_C += self.c_μ * wyy

    self.C = new_C
end

max_stars_repo_name: toshiki-kosen/SolvingTangramPuzzle

max_issues_repo_name: toshiki-kosen/SolvingTangramPuzzle

max_forks_repo_name: toshiki-kosen/SolvingTangramPuzzle


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMBase.jl/blob/master/LICENSE

using FEMBase
using FEMBase: get_problems
using Test

import FEMBase: run!

struct MyAnalysis1 <: AbstractAnalysis end

struct MyAnalysis2 <: AbstractAnalysis end

mutable struct MyProblem <: FieldProblem
    value :: Bool
end

MyProblem() = MyProblem(false)

function run!(analysis::Analysis{MyAnalysis2})
    for problem in get_problems(analysis)
        problem.properties.value = true
    end
end

@testset "test creating new analysis" begin
    analysis1 = Analysis(MyAnalysis1, "my test analysis 1")
    analysis2 = Analysis(MyAnalysis2, "my test analysis 2")
    problem = Problem(MyProblem, "my test problem", 1)
    add_problems!(analysis1, [problem])
    add_problems!(analysis2, [problem])
    run!(analysis1)
    @test problem.properties.value == false
    run!(analysis2)
    @test problem.properties.value == true
end

## Result writers

# Writing results of some analysis is done, in abstract level, by using
# results writers. They must be subtypes of `AbstractResultsWriter`:

mutable struct TestResultsWriter <: AbstractResultsWriter
end

# Results writer is added to analysis using `add_results_writer!`:

writer = TestResultsWriter()
analysis1 = Analysis(MyAnalysis1)
add_results_writer!(analysis1, writer)
@test first(get_results_writers(analysis1)) === writer

# Results are written with function `write_results!`-function, taking analysis
# as input argument. Function applies all results writers (there can be many!)
# to the analysis in the order they are attached to analysis using
# `add_results_writer!`. Internally the function runs
# `write_results!(analysis, results_writer)`, which must be implemented for each
# analysis / results_writer pair.

write_results!(analysis1)

# This of course doesn't do anything useful but inform that writing the results
# of MyAnalysis1 using TestResultsWriter is not implemented as the actual
# implementation is missing. In pseudo-level, the actual implementation would
# be something like:

function FEMBase.write_results!(analysis::MyAnalysis1, writer::TestResultsWriter)
    # fid = open(writer.filename)
    # for problem in get_problems(analysis)
    #     field = problem("displacement", analysis.time)
    #     write(fid, field)
    # end
    return nothing!
end

# Lastly, it is not mandatory to have a results writer in the analysis at all
# or the results writer can write directly to stdout or do some other fancy
# stuff. If the results writer is not defined, user gets only info message
# about missing results writers.

empty!(analysis1.results_writers)
write_results!(analysis1)


text:
const Callback = Union{Function, Nothing}

# Generate code for a global callback
macro callback(ex)
	ref = gensym()
	esc(quote
		const $ref = Ref{Callback}(nothing)
		$(callbackcode(extractargs(ex)..., :($ref[])))
	end)
end

# Number of callback functions per window
const _window_callback_num = Ref(0)

# Generate code for a window-specific callback
macro windowcallback(ex)
	_window_callback_num[] += 1
	i = _window_callback_num[]
	ref = :(callbacks(window)[$i])
	esc(callbackcode(extractargs(ex)..., ref, [:(window::Window)]))
end

# Generate expression with functions for [un]setting a callback
function callbackcode(
	name,
	callback_params,  # Signature of the C-compatible wrapper function
	callback_args,    # How the wrapper passes arguments to the callback
	callback_ref,     # Persisted reference to callback function
	setter_params=[], # Initial parameter(s) to setter (e.g. window handle)
)
	# Construct function names
	setter = Symbol("Set", name, "Callback")          # SetFooCallback
	libsetter = Expr(:quote, Symbol("glfw", setter))  # glfwSetFooCallback
	wrapper = Symbol('_', name, "CallbackWrapper")    # _FooCallbackWrapper

	# Separate names and types of parameters
	callback_param_types = Expr(:tuple, map(paramtype, callback_params)...)
	setter_param_names = map(paramname, setter_params)
	setter_param_types = map(paramtype, setter_params)

	quote
		# Callback wrapper that can be passed to `cfunction`
		$wrapper($(callback_params...)) = ($callback_ref($(callback_args...)); return nothing)

		# Set the callback function
		function $setter($(setter_param_names...), callback::Function)
			old_callback = $callback_ref
			$callback_ref = callback  # Prevent callback function from being garbage-collected
			cfunptr = @cfunction($wrapper, Cvoid, $callback_param_types)
			ccall(($libsetter, libglfw), Ptr{Cvoid}, ($(setter_param_types...), Ptr{Cvoid}), $(setter_param_names...), cfunptr)
			return old_callback
		end

		# Unset the callback function
		function $setter($(setter_param_names...), ::Nothing)
			ccall(($libsetter, libglfw), Ptr{Cvoid}, ($(setter_param_types...), Ptr{Cvoid}), $(setter_param_names...), C_NULL)
			old_callback = $callback_ref
			$callback_ref = nothing  # Allow former callback function to be garbage-collected
			return old_callback
		end
	end
end

function extractargs(ex)
	before, after = arrowsplit(ex)
	name = string(before.args[1])
	params = before.args[2:end]
	args = isa(after, Expr) ? after.args : map(paramname, params)
	name, params, args
end

function arrowsplit(ex)
	if ex.head == :->
		ex.args[1], ex.args[2].args[2]
	else
		ex, nothing
	end
end

paramname(param_ex) = param_ex.args[1]
paramtype(param_ex) = param_ex.args[2]


text:
module Example

using Glimmer, Glimmer.FlexUI
using SparseArrays
using LinearAlgebra
using Makie
import WGLMakie, JSServe

println("Start [$(splitext(basename(@__FILE__))[1])]")

#---------------------------------------------------------------
# function to convert a WGLMakie scene to HTML
#---------------------------------------------------------------
function renderWGLMakieHTML(obj)
    io = IOBuffer()
    show(io, MIME"text/html"(), JSServe.Page(exportable=true, offline=true))
    app = JSServe.App() do
        return JSServe.DOM.div(
            obj,
        )
    end
    show(io, MIME"text/html"(), app)
    page_html = String(take!(io))
end

#---------------------------------------------------------------
# rendering the two figures
#---------------------------------------------------------------
WGLMakie.activate!()
Makie.inline!(true)

fig1 = WGLMakie.surface(rand(4, 4))

x, y = collect(-8:0.5:8), collect(-8:0.5:8)
z = [sinc(√(X^2 + Y^2) / π) for X ∈ x, Y ∈ y]

fig2 = wireframe(x, y, z, axis=(type=Axis3,), color=:black)

#---------------------------------------------------------------
# define the application and some basic properties such as title and initial window size
#---------------------------------------------------------------
app = App()
prop!(app, :title, "Glimmer Example - Makie with WGLMakie backend")
prop!(app, :winInitWidth, 1200)
prop!(app, :winInitHeight, 800)

#---------------------------------------------------------------
# Define Controls
#---------------------------------------------------------------
ui = VContainer(
    Card(
        title="This is an example of a WGLMakie scene",
        subtitle="WGLMakie allow us to embed an interactive Makie figure into any webpage. Note: Not every figure is interactive (see thge second one), but according to the developers, more and more will becode so in the future.",
        content=VContainer(
            Card(
                title="Interactive Figure (use mouse to rotate)",
                content=VContainer(
                    RawHTML(                                # Container for figure 1
                        html=renderWGLMakieHTML(fig1),
                        useFrame=true,
                        style="width: 100%; height: 650px;"
                    ),
                ),
            ),
            Card(
                title="Non-Interactive Figure",
                content=VContainer(
                    RawHTML(
                        html=renderWGLMakieHTML(fig2),      # Container for figure 2
                        useFrame=true,
                        style="width: 100%; height: 650px;"
                    ),
                ),
            ),
        ),
    ),

    Glimmer.exampleSourceAsCard(@__FILE__),     # add the source code of the example as the last control
)
# set the controls for the application
controls!(app, ui)

#---------------------------------------------------------------
# the render function
#---------------------------------------------------------------
function render()
end
renderFunction!(app, render)

#---------------------------------------------------------------
# Run the application
#---------------------------------------------------------------
run(app)

println("End [$(splitext(basename(@__FILE__))[1])]")

end # module


text:
using StaticCompiler
using StaticTools

function times_table(argc::Int, argv::Ptr{Ptr{UInt8}})
    argc == 3 || return printf(stderrp(), c"Incorrect number of command-line arguments\n")
    rows = parse(Int64, argv, 2)            # First command-line argument
    cols = parse(Int64, argv, 3)            # Second command-line argument

    M = MallocArray{Int64}(undef, rows, cols)
    @inbounds for i=1:rows
        for j=1:cols
           M[i,j] = i*j
        end
    end
    # Print to stdout
    printf(M)
    # Also print to file
    fp = fopen(c"table.tsv",c"w")
    printf(fp, M)
    fclose(fp)
    # Clean up matrix
    free(M)
end

# Attempt to compile
path = compile_executable(times_table, (Int64, Ptr{Ptr{UInt8}}), "./")


text:
export SoftMax

"""
`SoftMax([domainType=Float64::Type,] dim_in::Tuple)`

Creates the softmax non-linear operator with input dimensions `dim_in`.
```math
\\sigma(\\mathbf{x}) = \\frac{e^{\\mathbf{x} }}{ \\sum e^{\\mathbf{x} } }
```

"""
struct SoftMax{T,N} <: NonLinearOperator
	dim::NTuple{N,Int}
	buf::Array{T,N}
end

function SoftMax(DomainType::Type, DomainDim::NTuple{N,Int}) where {N}
	SoftMax{DomainType,N}(DomainDim,zeros(DomainType,DomainDim))
end

SoftMax(DomainDim::NTuple{N,Int}) where {N} = SoftMax(Float64,DomainDim)

function mul!(y::AbstractArray{T,N}, L::SoftMax{T,N}, x::AbstractArray{T,N}) where {T,N}
	y .= exp.(x.-maximum(x))
	y ./= sum(y)
end

function mul!(y::AbstractArray,
              J::AdjointOperator{Jacobian{A,TT}},
              b::AbstractArray) where {T, N, A<: SoftMax{T,N}, TT <: AbstractArray{T,N} }
    L = J.A
	fill!(y,zero(T))
	L.A.buf .= exp.(L.x.-maximum(L.x))
	L.A.buf ./= sum(L.A.buf)
	for i in eachindex(y)
		y[i] = -L.A.buf[i]*dot(L.A.buf,b)
		y[i] += L.A.buf[i]*b[i]
	end
	return y
end

fun_name(L::SoftMax) = "σ"

size(L::SoftMax) = (L.dim, L.dim)

domainType(L::SoftMax{T,N}) where {T,N} = T
codomainType(L::SoftMax{T,N}) where {T,N} = T


text:
const INTERVAL_REGEX = r"""
    ^
    (\[|\()                      # Left bound type (1)
    (?|                          # Left bound value (2)
        "((?:[^"]|(?<=\\)")*?)"|
        ((?:[^,.]|\.(?!\.))*?)
    )
    (?:,|\ ?\.\.)\ ?             # Delimiter
    (?|                          # Right bound value (3)
        "((?:[^"]|(?<=\\)")*?)"|
        ((?:[^,.]|\.(?!\.))*?)
    )
    (\]|\))                      # Right bound type (4)
    $
    """x

"""
    parse(::Type{Interval{T}}, str; element_parser=parse) -> Interval{T}

Parse a string of the form `<left-type><left-value><delim><right-value><right-type>`
(e.g. `[1 .. 2)`) as an `Interval{T}`. The format above is interpreted as:

- `left-type`: Must be either "[" or "(" which indicates if the left-endpoint of the
  interval is either `Closed` or `Open`.

- `left-value`: Specifies the value of the left-endpoint which will be parsed as the type
  `T`. If the value string has a length of zero then the left-endpoint will be specified as
  `Unbounded`. If the value string contains the delimiter (see below) then you may
  double-quote the value string to avoid any ambiguity.

- `delim`: Must be either ".." or "," which indicates the delimiter separating the
  left/right endpoint values.

- `right-value`: Specifies the value of the right-endpoint. See `left-value` for more
  details.

- `right-type`: Must be either "]" or ")" which indicates if the right-endpoint of the
   interval is either `Closed` or `Open`.

The `element_parser` keyword allows a custom parser to be used when parsing the left/right
values. The function is expected to take two arguments: `Type{T}` and `AbstractString`.
This is useful for supplying additional arguments/keywords, alternative parser functions, or
for types that do not define `parse` (e.g. `String`).
"""
function Base.parse(::Type{Interval{T}}, str::AbstractString; element_parser=parse) where T
    m = match(INTERVAL_REGEX, str)

    if m === nothing
        throw(ArgumentError("Unable to parse \"$str\" as an interval"))
    end

    if isempty(m[2])
        L = Unbounded
        left = nothing
    else
        L = m[1] == "[" ? Closed : Open
        left = element_parser(T, m[2])
    end

    if isempty(m[3])
        R = Unbounded
        right = nothing
    else
        R = m[4] == "]" ? Closed : Open
        right = element_parser(T, m[3])
    end

    return Interval{T,L,R}(left, right)
end


text:
# A class of Matern isotropic functions including the Matrern 3/2 and 5/2, where d= 3 or 5. Also the exponential function, where d=1
abstract type MaternIso <: Isotropic{Euclidean} end
abstract type MaternARD <: StationaryARD{WeightedEuclidean} end

@inline function dKij_dθp(mat::MaternARD, X::MatF64, i::Int, j::Int, p::Int, dim::Int)
    r=distij(metric(mat),X,i,j,dim)
    if p <= dim
        wdiffp=dist2ijk(metric(mat),X,i,j,p)
        return wdiffp == 0.0 ? 0.0 : dk_dll(mat,r,wdiffp)
    elseif p==dim+1
        return dk_dlσ(mat, r)
    else
        return NaN
    end
end
@inline function dKij_dθp(mat::MaternARD, X::MatF64, data::StationaryARDData, i::Int, j::Int, p::Int, dim::Int)
    return dKij_dθp(mat,X,i,j,p,dim)
end

@inline function dk_dθp(mat::MaternIso, r::Float64, p::Int)
    if p==1
        return r == 0.0 ? 0.0 : dk_dll(mat, r)
    elseif p==2
        return dk_dlσ(mat, r)
    else
        return NaN
    end
end

include("mat12_iso.jl")
include("mat12_ard.jl")

include("mat32_iso.jl")
include("mat32_ard.jl")

include("mat52_iso.jl")
include("mat52_ard.jl")

"""
    Matern(ν::Float64, ll::Union{Float64,Vector{Float64}}, lσ::Float64)

Create Matérn kernel of type `ν` (i.e. `ν = 1/2`, `ν = 3/2`, or `ν = 5/2`) with length scale
`exp.(ll)` and signal standard deviation `exp(σ)`.

See also [`Mat12Iso`](@ref), [`Mat12Ard`](@ref), [`Mat32Iso`](@ref), [`Mat32Ard`](@ref),
[`Mat52Iso`](@ref), and [`Mat52Ard`](@ref).
"""
function Matern(ν::Float64, ll::Float64, lσ::Float64)
    if ν==1/2
        kern = Mat12Iso(ll, lσ)
    elseif ν==3/2
        kern = Mat32Iso(ll, lσ)
    elseif ν==5/2
        kern = Mat52Iso(ll, lσ)
    else throw(ArgumentError("Only Matern 1/2, 3/2 and 5/2 are implementable"))
    end
    return kern
end

function Matern(ν::Float64, ll::Vector{Float64}, lσ::Float64)
    if ν==1/2
        kern = Mat12Ard(ll, lσ)
    elseif ν==3/2
        kern = Mat32Ard(ll, lσ)
    elseif ν==5/2
        kern = Mat52Ard(ll, lσ)
    else throw(ArgumentError("Only Matern 1/2, 3/2 and 5/2 are implementable"))
    end
    return kern
end

@deprecate Mat Matern


text:
## Aetherlang core systems for the interpreter
include("AetherlangTypes.jl")
include("AetherlangExceptions.jl")

max_stars_repo_name: Aetherlang/aetherlang

max_issues_repo_name: gojakuch/aetherlang

max_forks_repo_name: gojakuch/aetherlang


text:
########## action ############
POMDPs.actions(DP::DrivePOMDP) = [:giveup, :takeover]
POMDPs.n_actions(DP::DrivePOMDP) = 2
function POMDPs.actionindex(DP::DrivePOMDP, a::Symbol)
    if a == :giveup
        return 1
    elseif a == :takeover
        return 2
    end
    error("Invalid action: $a")
end

function Acclimit!(DP::DrivePOMDP, acc::Float64) # will be moved inside function Actiontransform() later
    acc = min(acc, DP.Aset.max)
    acc = max(acc, DP.Aset.min)
end

function ToVref(DP::DrivePOMDP, Car::CarSt)
    step = UInt16(Car.s/DP.Δs)+1
    Vref = DP.Routes[DP.SsInit.Other.r].Vref[step] # the real route of other car is known to other car itself
    acc = IDM(Vego=Car.v, Vfront=Vref, Vref=Vref, Snet=Inf, T=1.0, Amax=DP.Aset.max, Bdec=DP.Aset.comfort)
    acc = Acclimit!(DP, acc)
end

# calculate the acceleration of the ego car using IDM based on stopline
function AccCalculate(DP::DrivePOMDP, Ego::CarSt, Aego::Symbol)
    acc_ego = 0.0
    Vref = DP.Routes[Ego.r].Vref[min(UInt16(floor(Ego.s/DP.Δs)+1), 21)]
    if Ego.s < DP.Stopline # strategy before stop line
        if Aego == :giveup
            acc_ego = IDM(Vego=Ego.v, Vfront=0.0, Vref=Vref, Snet=DP.Stopline-Ego.s, T=0.0, Amax=DP.Aset.max, Bdec=DP.Aset.comfort, Smin=0.0)
        elseif Aego == :takeover
            acc_ego = IDM(Vego=Ego.v, Vfront=Vref, Vref=Vref, Snet=Inf, T=1.0, Amax=DP.Aset.max, Bdec=DP.Aset.comfort, Smin=0.0)
        else
            error("Transform failed: Invalid Action $Aego !")
        end
    else # strategy after stop line
        if Aego == :giveup # brake with min a
            if Ego.v > 0
                acc_ego = DP.Aset.min
            else
                acc_ego == 0.0
            end
        elseif Aego == :takeover
            acc_ego = IDM(Vego=Ego.v, Vfront=Vref, Vref=Vref, Snet=Inf, T=1.0, Amax=DP.Aset.max, Bdec=DP.Aset.comfort, Smin=DP.Smin)
        else
            error("Transform failed: Invalid Action $Aego !")
        end
    end
    return acc_ego
end

# calculate the acceleration of the ego car using IDM based on Other and Ego cars; simulate ACC
function AccCalculate(DP::DrivePOMDP, Ego::CarSt, Other::CarSt, Aego::Symbol, Δs::Float64)
    acc_ego = 0.0
    Vref = DP.Routes[Ego.r].Vref[min(UInt16(floor(Ego.s/DP.Δs)+1), 21)]
    if Aego == :giveup # if giveup, ego car follows other car
        acc_ego = IDM(Vego=Ego.v, Vfront=min(Vref, Other.v), Vref=Vref, Snet=Δs, T=0.9, Amax=DP.Aset.max, Bdec=DP.Aset.comfort, Smin=DP.Smin)
    elseif Aego == :takeover # if take over, ego car drives as on free road.
        acc_ego = IDM(Vego=Ego.v, Vfront=Vref, Vref=Vref, Snet=Inf, T=1.0, Amax=DP.Aset.max, Bdec=DP.Aset.comfort, Smin=DP.Smin)
    else
        error("Transform failed: Invalid Action $Aego !")
    end
    return acc_ego
end


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: jl:light,ipynb
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.5.1
#   kernelspec:
#     display_name: Julia 1.4.2
#     language: julia
#     name: julia-1.4
# ---

# +
"""
Tutorial: A reference implementation of configuration interactions singles.
"""

__authors__   = ["D. Menendez", "Boyi Zhang", "Adam S. Abbott"]
__credits__   = ["D. Menendez", "Boyi Zhang", "Adam S. Abbott", "Justin M. Turney"]

__copyright_amp__ = "(c) 2014-2020, The Psi4Julia Developers"
__license__   = "BSD-3-Clause"
__date__      = "2020-08-03"
# -

# # Configuration Interaction Singles (CIS)

# ## I. Theoretical Overview

# In this tutorial, we will implement the configuration interaction singles method in the spin orbital notation. The groundwork for working in the spin orbital notation has been laid out in "Introduction to the Spin Orbital Formulation of Post-HF methods" [tutorial](../08_CEPA0_and_CCD/8a_Intro_to_spin_orbital_postHF.ipynb). It is highly recommended to work through that introduction before starting this tutorial.

# ### Configuration Interaction (CI)
#
# The configuration interaction wavefunction is constructed as a linear combination of the reference determinants and all singly, doubly, ... n-tuple excited determinants where n is the number of electrons in a given system:
#
# \begin{equation}
# \Psi_\mathrm{CI} = (1 + \hat{C_1} + \hat{C_2} + ...\hat{C_n)}\Phi
# \end{equation}
#
# Here, $\hat{C_n}$ is the n configuration excitation operator.
#
# In Full CI, all possible excitations are included in the wavefunction expansion. In truncated CI methods, only a subset of excitations are included.

# ## CIS
#
# In CIS, only single excitations from the occupied (indices i,j,k...) to the virtual (indices a,b,c...) orbitals are included. As a result, CIS gives transition energies to an excited state.
#
# Assuming we are using canonical Hartree-Fock spin orbitals($\{\mathrm{\psi_p}\}$) with orbital energies $\{\epsilon_p\}$, we can build a shifted CIS Hamiltonian matrix:
#
# \begin{equation}
# \tilde{\textbf{H}} = \textbf{H} - E_0 \textbf{I} = [\langle \Phi_P | \hat{H_e} - E_0|\Phi_Q \rangle],\,
# \Phi_P \in {\Phi_i^a}
# \end{equation}
#
# where $E_0$ is the ground state Hartree-Fock state energy given by $\langle \Phi | \hat{H_e}|\Phi \rangle$.
#
# The matrix elements of this shifted CIS Hamiltonian matrix can be evaluated using Slater's rules to give:
#
# \begin{equation}
# \langle \Phi_i^a | \hat{H_e} - E_0|\Phi_j^b \rangle = (\epsilon_a - \epsilon_i)\delta_{ij} \delta_{ab}
# + \langle aj || ib \rangle
# \end{equation}
#
# This then becomes a standard eigenvalue equation from which we can solve for the excitation energies and the wavefunction expansion coefficients:
#
# \begin{equation}
# \tilde{\textbf{H}} \textbf{c}_K = \Delta E_K\textbf{c}_K, \,\Delta E_K = E_K - E_0
# \end{equation}
#
#

# ## II. Implementation

# As with previous tutorials, let's begin by importing Psi4, NumPy, and TensorOperations and setting memory and output file options.

# +
# ==> Import Psi4, NumPy, & TensorOperations <==
using PyCall: pyimport
psi4 = pyimport("psi4")
np   = pyimport("numpy") # used only to cast to Psi4 arrays
using TensorOperations: @tensor

# ==> Set Basic Psi4 Options <==

# Memory specifications
psi4.set_memory(Int(2e9))
numpy_memory = 2

# Output options
psi4.core.set_output_file("output.dat", false)
# -

# We now define the molecule and set Psi4 options:

# +
mol = psi4.geometry("""
0 1
O
H 1 1.1
H 1 1.1 2 104
symmetry c1
""")

psi4.set_options(Dict("basis"         => "sto-3g",
                      "scf_type"      => "pk",
                      "reference"     => "rhf",
                      "mp2_type"      => "conv",
                      "e_convergence" => 1e-8,
                      "d_convergence" => 1e-8))
# -

# We use Psi4 to compute the RHF energy and wavefunction and store them in variables `scf_e` and `scf_wfn`. We also check the memory requirements for computation:

# +
# Get the SCF wavefunction & energies
scf_e, scf_wfn = psi4.energy("scf", return_wfn=true)

# Check memory requirements
nmo = scf_wfn.nmo()
I_size = nmo^4 * 8e-9
println("\nSize of the ERI tensor will be $I_size GB.")
memory_footprint = I_size * 1.5
if I_size > numpy_memory
    psi4.core.clean()
    throw(OutOfMemoryError("Estimated memory utilization ($memory_footprint GB) exceeds " *
                           "allotted memory limit of $numpy_memory GB."))
end
# -

# We first obtain orbital information from our wavefunction. We also create an instance of MintsHelper to help build our molecular integrals:

# +
# Create instance of MintsHelper class
mints = psi4.core.MintsHelper(scf_wfn.basisset())

# Get basis and orbital information
nbf = mints.nbf()          # Number of basis functions
nalpha = scf_wfn.nalpha()  # Number of alpha electrons
nbeta = scf_wfn.nbeta()    # Number of beta electrons
nocc = nalpha + nbeta      # Total number of electrons
nso = 2nbf                 # Total number of spin orbitals
nvirt = nso - nocc         # Number of virtual orbitals
# -

# We now build our 2-electron integral, a 4D tensor, in the spin orbital formulation. We also convert it into physicist's notation and antisymmetrize for easier manipulation of the tensor later on.

# +
"""
Spin blocks 2-electron integrals
Using np.kron, we project I and I tranpose into the space of the 2x2 ide
The result is our 2-electron integral tensor in spin orbital notation
"""
function spin_block_tei(I)
    identity = [ 1.0 0.0; 0.0 1.0]
    I = np.kron(identity, I)
    np.kron(identity, permutedims(I, reverse(1:4)))
end

I = np.asarray(mints.ao_eri())
I_spinblock = spin_block_tei(I)

# Convert chemist's notation to physicist's notation, and antisymmetrize
# (pq | rs) ---> <pr | qs>
# <pr||qs> = <pr | qs> - <pr | sq>
gao = permutedims(I_spinblock, (1, 3, 2, 4)) - permutedims(I_spinblock, (1, 3, 4, 2));
# -

# We get the orbital energies from alpha and beta electrons and append them together. We spin-block the coefficients obtained from the reference wavefunction and convert them into NumPy arrays. There is a set corresponding to coefficients from alpha electrons and a set of coefficients from beta electrons. We then sort them according to the order of the orbital energies using argsort():

# +
# Get orbital energies, cast into NumPy array, and extend eigenvalues
eps_a = np.asarray(scf_wfn.epsilon_a())
eps_b = np.asarray(scf_wfn.epsilon_b())
eps = vcat(eps_a, eps_b)

# Get coefficients, block, and sort
Ca = np.asarray(scf_wfn.Ca())
Cb = np.asarray(scf_wfn.Cb())
C = [Ca zero(Ca); zero(Cb) Cb]; # direct sum

# Sort the columns of C according to the order of orbital energies
C = C[:, sortperm(eps)]

# Sort orbital energies
sort!(eps);
# -

# We now transform the 2-electron integral from the AO basis into the MO basis using the coefficients:
#

# Transform gao, which is the spin-blocked 4d array of physicist's notation,
# antisymmetric two-electron integrals, into the MO basis using MO coefficients
gmo = @tensor begin
   gmo[P,Q,R,S] := gao[p,Q,R,S] * C[p,P]
   gmo[p,Q,R,S] := gmo[p,q,R,S] * C[q,Q]
   gmo[p,q,R,S] := gmo[p,q,r,S] * C[r,R]
   gmo[p,q,r,S] := gmo[p,q,r,s] * C[s,S]
end
nothing

# Now that we have our integrals, coefficents, and orbital energies set up in with spin orbitals, we can start our CIS procedure. We first start by initializing the shifted Hamiltonion matrix $\tilde{\textbf{H}}$ (`HCIS`). Let's think about the size of $\tilde{\textbf{H}}$. We need all possible single excitations from the occupied to virtual orbitals. This is given by the number of occupied orbitals times the number of virtual orbitals  (`nocc * nvirt`).
#
# The size of our matrix is thus `nocc * nvirt` by `nocc * nvirt`.

# Initialize CIS matrix.
# The dimensions are the number of possible single excitations
HCIS = zeros(nocc * nvirt, nocc * nvirt);

# Next, we want to build all possible excitations from occupied to virtual orbitals. We create two for-loops that will loop over the number of occupied orbitals and number of virtual orbitals, respectively, and store the combination of occupied and virtual indices as a tuple `(i, a)`. We put all tuples in a list called `excitations`.

# Build the possible excitations, collect indices into a list
excitations = []
for i in 1:nocc, a in nocc+1:nso
    push!(excitations,(i, a))
end

# Now we can evaluate the matrix elements of the shifted CIS Hamiltonian matrix using the equation given above. For each element, there are several layers of indexing that we have to consider.
# First, there are the indices of the element itself, which gives the position of the element in the matrix. Indices `p` and `q` are used:
#
# `HCIS[p, q]`
#
# Second, there are two sets of excitations from occupied to virtual orbitals corresponding to the bra and ket of each matrix element. For these, we will take advantage of the `excitations` list that we build with the list of all possible excitations. We will use indices `i` and `a` to denote the excitation in the bra (`left_excitation`) and `j` and `b` to denote the excitation in the ket (`right_excitation`).
#
# To manage these indices, we will use the `enumerate` function.
#
# Note that a Kronecker delta $\delta_{pq}$ can be represented as `p == q`.
#

# Form matrix elements of shifted CIS Hamiltonian
for (p, left_excitation) in enumerate(excitations)
    i, a = left_excitation
    for (q, right_excitation) in enumerate(excitations)
        j, b = right_excitation
        HCIS[p, q] = (eps[a] - eps[i]) * (i == j) * (a == b) + gmo[a, j, i, b]
    end
end

# We now use the composed function `eigen ∘ Hermitian` (for hermitian matrices) to diagonalize the shifted CIS Hamiltonian. This will give us the excitation energies (`ECIS`). These eigenvalues correspond to the CIS total energies for various excited states. The columns of matrix `CCIS` give us the coefficients which describe the relative contribution of each singly excited determinant to the excitation energy.
#

# +
# Diagonalize the shifted CIS Hamiltonian
using LinearAlgebra: eigen, Hermitian

ECIS, CCIS = (eigen ∘ Hermitian)(HCIS) ;
# -

# For a given excitation energy, each coefficent in the linear combination of excitations represents the amount that a particular excitation contributes to the overall excitation energy. The percentage contribution of each coefficient can be calculated by squaring the coefficent and multiplying by 100.

# Percentage contributions for each state vector
percent_contrib = @. round(CCIS^2 * 100);

# In addition to excitation energies, we want to print the excitations that contribute 10% or more to the overall energy, as well as their percent contribution.
#
# Note that `printfmt` allows us to print different sections to the same line without a line break.

# Print detailed information on significant excitations
using Formatting: printfmt
println("CIS:")
for state in eachindex(ECIS)
    # Print state, energy
    printfmt("State {1:3d} Energy (Eh) {2:10.7f} ", state, ECIS[state])
    for (idx, excitation) in enumerate(excitations)
        if percent_contrib[idx, state] > 10
            i, a = excitation
            # Print percentage contribution and the excitation
            printfmt("{1:4d}% {2:2d} -> {3:2d} ", percent_contrib[idx, state], i, a)
        end
    end
    printfmt("\n")
end

# ## References
# 1. Background paper:
#  >"Toward a systematic molecular orbital theory for excited states"
# [[Foresman:1992:96](http://pubs.acs.org/doi/abs/10.1021/j100180a030)] J. B. Foresman, M. Head-Gordon, J. A. Pople, M. J. Frisch, *J. Phys. Chem.* **96**, 135 (1992).
#
#
# 2. Algorithms from:
# 	> [[CCQC:CIS](https://github.com/CCQC/summer-program/tree/master/7)] CCQC Summer Program, "CIS" accessed with https://github.com/CCQC/summer-program/tree/master/7.
#