text:
#needed within OrdinaryDiffEq somehow:
import Base: zero
zero(tmp::Array) = zero.(tmp)

#needed to avoid allocations:
Flo=Union{Float32,Float64}
mydt(tim::Flo,𝑇::Array{Float32,1})=(tim-𝑇[1])/(𝑇[2]-𝑇[1])
mydt(tim::Flo,𝑇::Array{Float64,1})=(tim-𝑇[1])/(𝑇[2]-𝑇[1])

#needed to avoid allocations:
fSize(f::Array{Tuple{Int,Int}},i::Int) = (f[i][1],f[i][2])
fview(f::Array{Array{Float32,2},1},i::Int) = view(f[i],:,:)
fview(f::Array{Array{Float64,2},1},i::Int) = view(f[i],:,:)
fview(f::Array{Array{Float32,2},2},i::Int,j::Int) = view(f[i,j],:,:)
fview(f::Array{Array{Float64,2},2},i::Int,j::Int) = view(f[i,j],:,:)

"""
    dxyz_dt!(du,u,p::𝐹_MeshArray3D,tim)

Interpolate velocity from gridded fields (3D; with halos) to position `u`
(`x,y,z,fIndex`) to compute the derivative of position v time  `du_dt`.

```jldoctest
using IndividualDisplacements, Statistics
p=dirname(pathof(IndividualDisplacements))
include(joinpath(p,"../examples/worldwide/three_dimensional_ocean.jl"))
ref=[211. 34. -70.]
prod(isapprox.([mean(𝐼.🔴.lon) mean(𝐼.🔴.lat) mean(𝐼.🔴.z)],ref,atol=50.0))

# output

true
```
"""
function dxyz_dt!(du::Array{T,1},u::Array{T,1},𝑃::𝐹_MeshArray3D,tim) where T
    dt=mydt(tim,𝑃.𝑇)
    g=𝑃.u0.grid
    #
    while location_is_out(u,g)
        𝑃.update_location!(u)
    end

    x=u[1]
    y=u[2]
    z=u[3]
    fIndex = Int(u[4])
    (nx,ny)=fSize(g.fSize,fIndex)
    nz=size(𝑃.u0,2)
    #
    (dx,dy,dz)=(x - floor(x),y - floor(y),z - floor(z))
    i_c = Int32(floor(x)) + 2
    j_c = Int32(floor(y)) + 2
    k_c = Int32(floor(z)) + 1
    #
    (i_w,i_e)=(i_c,i_c+1)
    (j_s,j_n)=(j_c,j_c+1)
    (k_l,k_r)=(k_c,k_c+1)
    #
    k_c=min(k_c,nz)
    k_l=min(k_l,nz)
    k_r=min(k_r,nz)
    k_c=max(k_c,1)
    k_l=max(k_l,1)
    k_r=max(k_r,1)
    #
    onemdt=(1.0-dt)
    onemdx=(1.0-dx)
    onemdy=(1.0-dy)
    onemdz=(1.0-dz)
    #interpolate u to position and time
    u0=fview(𝑃.u0.f,fIndex,k_c)
    u1=fview(𝑃.u1.f,fIndex,k_c)
    du[1]=onemdx*onemdt*u0[i_w,j_c]+
    dx*onemdt*u0[i_e,j_c]+
    onemdx*dt*u1[i_w,j_c]+
    dx*dt*u1[i_e,j_c]
    #interpolate v to position and time
    v0=fview(𝑃.v0.f,fIndex,k_c)
    v1=fview(𝑃.v1.f,fIndex,k_c)
    du[2]=onemdy*onemdt*v0[i_c,j_s]+
    dy*onemdt*v0[i_c,j_n]+
    onemdy*dt*v1[i_c,j_s]+
    dy*dt*v1[i_c,j_n]
    #interpolate w to position and time
    du[3]=onemdz*onemdt*𝑃.w0.f[fIndex,k_l][i_c,j_c]+
    dz*onemdt*𝑃.w0.f[fIndex,k_r][i_c,j_c]+
    onemdz*dt*𝑃.w1.f[fIndex,k_l][i_c,j_c]+
    dz*dt*𝑃.w1.f[fIndex,k_r][i_c,j_c]
    #leave face index unchanged
    du[4]=0.0
    #
    return du
end

function dxyz_dt!(du::Array{T,2},u::Array{T,2},𝑃::𝐹_MeshArray3D,tim) where T
    [dxyz_dt!(du[i],u[i],𝑃,tim) for i=1:size(u,2)]
end

"""
    dxy_dt!(du,u,p::𝐹_MeshArray2D,tim)

Interpolate velocity from gridded fields (2D; with halos) to position `u`
(`x,y,fIndex`) to compute the derivative of position v time  `du_dt`.

```jldoctest
using IndividualDisplacements, Statistics
p=dirname(pathof(IndividualDisplacements))
include(joinpath(p,"../examples/flow_fields.jl"));

u,v,ϕ=random_flow_field()
#𝐹=𝐹_Array2D(u,u,v,v,[0.,10.])
𝐹=convert_to_FlowFields(u,v,10.0)

np,nq=size(u)
x=np*(0.4 .+ 0.2*rand(100))
y=nq*(0.4 .+ 0.2*rand(100))

a=ones(size(x))
𝐼=Individuals(𝐹,x,y,a)
∫!(𝐼)

ref=size(u) ./2
prod(isapprox.([mean(𝐼.🔴.x) mean(𝐼.🔴.y)],ref,atol=10.0))

# output

true
```

```jldoctest
using IndividualDisplacements, Statistics
p=dirname(pathof(IndividualDisplacements))
include(joinpath(p,"../examples/worldwide/global_ocean_circulation.jl"))
ref=[78. 88.]
prod(isapprox.([mean(𝐼.🔴.x) mean(𝐼.🔴.y)],ref,atol=10.0))

# output

true
```
"""
function dxy_dt!(du::Array{T,1},u::Array{T,1},𝑃::𝐹_MeshArray2D,tim) where T
    #compute positions in index units
    dt=mydt(tim,𝑃.𝑇)
    g=𝑃.u0.grid
    #
    while location_is_out(u,g)
        𝑃.update_location!(u)
    end

    x=u[1]
    y=u[2]
    fIndex = Int(u[3])
    (nx,ny)=fSize(g.fSize,fIndex)
    #
    (dx,dy)=(x - floor(x),y - floor(y))
    i_c = Int32(floor(x)) + 2
    j_c = Int32(floor(y)) + 2
    #
    (i_w,i_e)=(i_c,i_c+1)
    (j_s,j_n)=(j_c,j_c+1)
    #interpolate u to position and time
    u0=fview(𝑃.u0.f,fIndex)
    u1=fview(𝑃.u1.f,fIndex)
    du[1]=(1.0-dx)*(1.0-dt)*u0[i_w,j_c]+
    dx*(1.0-dt)*u0[i_e,j_c]+
    (1.0-dx)*dt*u1[i_w,j_c]+
    dx*dt*u1[i_e,j_c]
    #interpolate v to position and time
    v0=fview(𝑃.v0.f,fIndex)
    v1=fview(𝑃.v1.f,fIndex)
    du[2]=(1.0-dy)*(1.0-dt)*v0[i_c,j_s]+
    dy*(1.0-dt)*v0[i_c,j_n]+
    (1.0-dy)*dt*v1[i_c,j_s]+
    dy*dt*v1[i_c,j_n]
    #leave face index unchanged
    du[3]=0.0
    #
    return du
end

function dxy_dt!(du::Array{T,2},u::Array{T,2},𝑃::𝐹_MeshArray2D,tim) where T
    [dxy_dt!(du[i],u[i],𝑃,tim) for i=1:size(u,2)]
end

"""
    dxyz_dt(du,u,𝑃::𝐹_Array3D,tim)

Interpolate velocity from gridded fields (3D; NO halos) to position `u`
(`x,y,z`) to compute the derivative of position v time  `du_dt`.

```jldoctest
using IndividualDisplacements
p=dirname(pathof(IndividualDisplacements))
include(joinpath(p,"../examples/basics/solid_body_rotation.jl"))
ref=[7.767441577479032 9.513402495574852 0.7065855989421701]
prod(isapprox.(𝐼.📌',ref,atol=1.0))

# output

true
```
"""
function dxyz_dt(du::Array{T,1},u::Array{T,1},𝑃::𝐹_Array3D,tim) where T
    #compute positions in index units
    dt=mydt(tim,𝑃.𝑇)
    #
    (nx,ny,nz) = size(𝑃.u0)
    x=mod(u[1],nx)
    y=mod(u[2],ny)
    z=mod(u[3],nz)
    #
    dx=x - floor(x)
    dy=y - floor(y)
    dz=z - floor(z)
    #
    i_c = Int32(floor(x)) + 1
    j_c = Int32(floor(y)) + 1
    k_c = Int32(floor(z)) + 1
    #
    i_c==0 ? i_c=nx : nothing
    j_c==0 ? j_c=ny : nothing
    #k_c==0 ? k_c=nz : nothing
    #
    (i_w,i_e)=(i_c,i_c+1)
    x>=nx-1 ? (i_w,i_e)=(nx,1) : nothing
    (j_s,j_n)=(j_c,j_c+1)
    y>=ny-1 ? (j_s,j_n)=(ny,1) : nothing
    (k_l,k_r)=(k_c,k_c+1)
    #z>=nz-1 ? (k_l,k_r)=(nz,1) : nothing

    #interpolate u to position and time
    du[1]=(1.0-dx)*(1.0-dt)*𝑃.u0[i_w,j_c,k_c]+
    dx*(1.0-dt)*𝑃.u0[i_e,j_c,k_c]+
    (1.0-dx)*dt*𝑃.u1[i_w,j_c,k_c]+
    dx*dt*𝑃.u1[i_e,j_c,k_c]
    #interpolate v to position and time
    du[2]=(1.0-dy)*(1.0-dt)*𝑃.v0[i_c,j_s,k_c]+
    dy*(1.0-dt)*𝑃.v0[i_c,j_n,k_c]+
    (1.0-dy)*dt*𝑃.v1[i_c,j_s,k_c]+
    dy*dt*𝑃.v1[i_c,j_n,k_c]
    #interpolate w to position and time
    du[3]=(1.0-dz)*(1.0-dt)*𝑃.w0[i_c,j_c,k_l]+
    dz*(1.0-dt)*𝑃.w0[i_c,j_c,k_r]+
    (1.0-dz)*dt*𝑃.w1[i_c,j_c,k_l]+
    dz*dt*𝑃.w1[i_c,j_c,k_r]
    #
    return du
end

function dxyz_dt(du::Array{T,2},u::Array{T,2},𝑃::𝐹_Array3D,tim) where T
    [dxyz_dt(du[i],u[i],𝑃,tim) for i=1:size(u,2)]
end

"""
    dxy_dt(du,u,𝑃::𝐹_Array2D,tim)

Interpolate velocity from gridded fields (2D; NO halos) to position `u`
(`x,y`) to compute the derivative of position v time  `du_dt`.

```jldoctest
using IndividualDisplacements, Statistics
p=dirname(pathof(IndividualDisplacements))
include(joinpath(p,"../examples/basics/particle_cloud.jl"))
ref=[29.381183342468674  19.890831699436823]
prod(isapprox.([mean(𝐼.🔴.x) mean(𝐼.🔴.y)],ref,atol=1.0))

# output

true
```
"""
function dxy_dt(du::Array{T,1},u::Array{T,1},𝑃::𝐹_Array2D,tim) where T
    dt=mydt(tim,𝑃.𝑇)
    #
    (nx,ny) = size(𝑃.u0)
    x=mod(u[1],nx)
    y=mod(u[2],ny)
    #
    dx=x - floor(x)
    dy=y - floor(y)
    #
    i_c = Int32(floor(x)) + 1
    j_c = Int32(floor(y)) + 1
    #
    i_c==0 ? i_c=nx : nothing
    j_c==0 ? j_c=ny : nothing
    #
    (i_w,i_e)=(i_c,i_c+1)
    x>=nx-1 ? (i_w,i_e)=(nx,1) : nothing
    (j_s,j_n)=(j_c,j_c+1)
    y>=ny-1 ? (j_s,j_n)=(ny,1) : nothing
    #interpolate u to position and time
    du[1]=(1.0-dx)*(1.0-dt)*𝑃.u0[i_w,j_c]+
    dx*(1.0-dt)*𝑃.u0[i_e,j_c]+
    (1.0-dx)*dt*𝑃.u1[i_w,j_c]+
    dx*dt*𝑃.u1[i_e,j_c]
    #interpolate v to position and time
    du[2]=(1.0-dy)*(1.0-dt)*𝑃.v0[i_c,j_s]+
    dy*(1.0-dt)*𝑃.v0[i_c,j_n]+
    (1.0-dy)*dt*𝑃.v1[i_c,j_s]+
    dy*dt*𝑃.v1[i_c,j_n]
    #
    return du
end

function dxy_dt(du::Array{T,2},u::Array{T,2},𝑃::𝐹_Array2D,tim) where T
    [dxy_dt(du[i],u[i],𝑃,tim) for i=1:size(u,2)]
end

"""
    dxy_dt_CyclicArray(du,u,𝑃::NamedTuple,tim)

Nearest neighbor (?) velocity from gridded fields (2D; NO halos but
not needed when CyclicArrays is used to extend valid indice ranges).

_notes:_ spatial interpolation & temporal interpolation are lacking

```jldoctest
using IndividualDisplacements, Statistics
p=dirname(pathof(IndividualDisplacements))
include(joinpath(p,"../examples/example_CyclicArray.jl"))
(x,y)=cyclicarray_example()
ref=[270. 243.]
prod(isapprox.([mean(x) mean(y)],ref,atol=1.0))

# output

true
```
"""
function dxy_dt_CyclicArray(du::Array{T,2},u::Array{T,2},𝑃::NamedTuple,tim) where T
    np=size(du,2)
    xi,yi=(u[1,:],u[2,:])

    #not needed:
    #xi=floor.(𝑃.xg[1,Int.(sign.(xi).*floor.(abs.(xi)))])+rem.(xi,1)
    #yi=floor.(𝑃.yg[Int.(sign.(yi).*floor.(abs.(yi))),1])+rem.(yi,1)

    i=Int.(floor.(𝑃.xg[1,Int.(floor.(xi))]))
    j=Int.(floor.(𝑃.yg[Int.(floor.(yi)),1]))
    du[1,:]=[𝑃.u[i[ii],j[ii]] for ii in 1:np]
    du[2,:]=[𝑃.v[i[ii],j[ii]] for ii in 1:np]
    return du
end

"""
    dict_to_nt(tmp::Dict)

Attempt to convert from Dict to NamedTuple
"""
dict_to_nt(tmp::Dict) = tmp=(; zip(Symbol.(keys(tmp)), values(tmp))...)

dxy_dt!(du,u,𝑃::Dict,tim) = dxy_dt!(du,u,dict_to_nt(𝑃),tim)
dxy_dt(du,u,𝑃::Dict,tim) = dxy_dt(du,u,dict_to_nt(𝑃),tim)
dxyz_dt(du,u,𝑃::Dict,tim) = dxyz_dt(du,u,dict_to_nt(𝑃),tim)

"""
    dxy_dt_replay(du,u,p::DataFrame,t)

Interpolate velocity from MITgcm float_trajectories output and return
position increment `du`.

```jldoctest
using IndividualDisplacements, Statistics
p=dirname(pathof(IndividualDisplacements))
include(joinpath(p,"../examples/basics/detailed_look.jl"))
prod(isapprox.(sol[:,end],ref[:,end],atol=1.0))

# output

true
```
"""
function dxy_dt_replay(du,u,p::DataFrame,t)
    tt=t/3600.0
    tt0=Int32(floor(tt))
    w=tt-tt0
    du[1]=(1.0-w)*p[tt0+1,:uVel]+w*p[tt0+2,:uVel]
    du[2]=(1.0-w)*p[tt0+1,:vVel]+w*p[tt0+2,:vVel]
end


text:
# This file contains function to convert events dataframe to SPADL representation.

include("types.jl")
include("tags.jl")
include("actions.jl")
include("positions.jl")
include("fixing_events.jl")

"""
    function convert_to_spadl()

Convert Wyscout events to SPADL actions
"""
function convert_to_spadl(event_df::Vector{WyscoutEvent})::Vector{RegularSPADL}
    # Create the vector of spadl actions
    # spadl_df = [RegularSPADL() for i in eachindex(event_df)]
    # tags_df = get_tags(event_df)

    # spadl_df = make_new_positions(spadl_df, tags_df, event_df)

    return action_df
end

ext: jl


text:
function extract_kwargs(kwargs, keys)
    return (k => v for (k, v) in pairs(kwargs) if k in keys)
end

extract_fields(x) = Dict(String(name) => getfield(x, name) for name in fieldnames(typeof(x)))

# Trailing directory separator matters for `basename`, so let's strip it to
# avoid confusion
strip_path_separator(path::AbstractString) = isdirpath(path) ? dirname(path) : path

function valid_dl_path(path, platform)
    try
        parse_dl_name_version(path, os(platform))
        return true
    catch
        return false
    end
end

# If we don't have `stdlib_version` from Pkg, recreate it ourselves
if !isdefined(Pkg.Types, :stdlib_version)
    function stdlib_version(uuid::Base.UUID, julia_version::Union{VersionNumber,Nothing})::Union{VersionNumber,Nothing}
        last_stdlibs = Pkg.Types.get_last_stdlibs(julia_version)
        if !(uuid in keys(last_stdlibs))
            return nothing
        end
        return last_stdlibs[uuid][2]
    end
else
    const stdlib_version = Pkg.Types.stdlib_version
end

if isdefined(Pkg, :Registry) && isdefined(Pkg.Registry, :registry_info)
    const registry_info = Pkg.Registry.registry_info
elseif isdefined(Pkg, :RegistryHandling) && isdefined(Pkg.RegistryHandling, :registry_info)
    const registry_info = Pkg.RegistryHandling.registry_info
end


text:
# Bogumił Kamiński, 2019-03-25

using Statistics, Distributions

function v_asian_sample1(T, r, K, σ, X₀, m)
    X = X₀
    hatx = zero(X)
    Δ = T / m
    for i in 1:m
        X *= exp((r-σ^2/2)*Δ + σ*√Δ*randn())
        hatx += X
    end
    exp(-r*T)*max(hatx/m - K, 0)
end

function v_asian_sample2(T, r, K, σ, X₀, m)
    Δ = T / m
    X = rand(Normal((r-σ^2/2)*Δ, σ*√Δ), m)
    exp(-r*T)*max(mean(exp.(cumsum(X))) * X₀ - K, 0)
end

function v_asian(T, r, K, σ, X₀, m, n, sampler)
    mean(sampler(T, r, K, σ, X₀, m) for i in 1:n)
end

for f in [v_asian_sample1, v_asian_sample2]
    @time println(v_asian(1.0, 0.05, 55.0, 0.3, 50.0, 10_000, 10_000, f))
end

max_stars_repo_path: 2. Asianoption/asianoption.jl

max_issues_repo_path: 2. Asianoption/asianoption.jl

max_forks_repo_path: 2. Asianoption/asianoption.jl


text:
using ArrayInterfaceCore
using ArrayInterfaceCore: zeromatrix
import ArrayInterfaceCore: has_sparsestruct, findstructralnz, fast_scalar_indexing, lu_instance,
        parent_type, zeromatrix
using Base: setindex
using LinearAlgebra
using Random
using SparseArrays
using Test

using Aqua
Aqua.test_all(ArrayInterfaceCore)

@test zeromatrix(rand(4,4,4)) == zeros(4*4*4,4*4*4)

@testset "matrix colors" begin
    @test ArrayInterfaceCore.fast_matrix_colors(1) == false
    @test ArrayInterfaceCore.fast_matrix_colors(Diagonal{Int,Vector{Int}})

    @test ArrayInterfaceCore.matrix_colors(Diagonal([1,2,3,4])) == [1, 1, 1, 1]
    @test ArrayInterfaceCore.matrix_colors(Bidiagonal([1,2,3,4], [7,8,9], :U)) == [1, 2, 1, 2]
    @test ArrayInterfaceCore.matrix_colors(Tridiagonal([1,2,3],[1,2,3,4],[4,5,6])) == [1, 2, 3, 1]
    @test ArrayInterfaceCore.matrix_colors(SymTridiagonal([1,2,3,4],[5,6,7])) == [1, 2, 3, 1]
    @test ArrayInterfaceCore.matrix_colors(rand(4,4)) == Base.OneTo(4)
end

@testset "parent_type" begin
    x = ones(4, 4)
    @test parent_type(view(x, 1:2, 1:2)) <: typeof(x)
    @test parent_type(reshape(x, 2, :)) <: typeof(x)
    @test parent_type(transpose(x)) <: typeof(x)
    @test parent_type(Symmetric(x)) <: typeof(x)
    @test parent_type(UpperTriangular(x)) <: typeof(x)
    @test parent_type(PermutedDimsArray(x, (2,1))) <: typeof(x)
    @test parent_type(Base.Slice(1:10)) <: UnitRange{Int}
    @test parent_type(Diagonal{Int,Vector{Int}}) <: Vector{Int}
    @test parent_type(UpperTriangular{Int,Matrix{Int}}) <: Matrix{Int}
    @test parent_type(LowerTriangular{Int,Matrix{Int}}) <: Matrix{Int}
end

@testset "buffer" begin
    @test ArrayInterfaceCore.buffer(sparse([1,2,3],[1,2,3],[1,2,3])) == [1, 2, 3]
    @test ArrayInterfaceCore.buffer(sparsevec([1, 2, 0, 0, 3, 0])) == [1, 2, 3]
    @test ArrayInterfaceCore.buffer(Diagonal([1,2,3])) == [1, 2, 3]
end

@test ArrayInterfaceCore.can_avx(ArrayInterfaceCore.can_avx) == false

@testset "lu_instance" begin
    A = randn(5, 5)
    @test lu_instance(A) isa typeof(lu(A))
    A = sprand(50, 50, 0.5)
    @test lu_instance(A) isa typeof(lu(A))
    @test lu_instance(1) === 1
end

@testset "ismutable" begin
    @test ArrayInterfaceCore.ismutable(rand(3))
    @test ArrayInterfaceCore.ismutable((0.1,1.0)) == false
    @test ArrayInterfaceCore.ismutable(Base.ImmutableDict{Symbol,Int64}) == false
    @test ArrayInterfaceCore.ismutable((;x=1)) == false
    @test ArrayInterfaceCore.ismutable(UnitRange{Int}) == false
    @test ArrayInterfaceCore.ismutable(Dict{Any,Any})
    @test ArrayInterfaceCore.ismutable(spzeros(1, 1))
    @test ArrayInterfaceCore.ismutable(spzeros(1))
end

@testset "can_change_size" begin
    @test ArrayInterfaceCore.can_change_size([1])
    @test ArrayInterfaceCore.can_change_size(Vector{Int})
    @test ArrayInterfaceCore.can_change_size(Dict{Symbol,Any})
    @test !ArrayInterfaceCore.can_change_size(Base.ImmutableDict{Symbol,Int64})
    @test !ArrayInterfaceCore.can_change_size(Tuple{})
end

@testset "can_setindex" begin
    @test !@inferred(ArrayInterfaceCore.can_setindex(1:2))
    @test @inferred(ArrayInterfaceCore.can_setindex(Vector{Int}))
    @test !@inferred(ArrayInterfaceCore.can_setindex(UnitRange{Int}))
    @test !@inferred(ArrayInterfaceCore.can_setindex(Base.ImmutableDict{Int,Int}))
    @test !@inferred(ArrayInterfaceCore.can_setindex(Tuple{}))
    @test !@inferred(ArrayInterfaceCore.can_setindex(NamedTuple{(),Tuple{}}))
    @test @inferred(ArrayInterfaceCore.can_setindex(Dict{Int,Int}))
end

@testset "fast_scalar_indexing" begin
    @test !fast_scalar_indexing(qr(rand(10, 10)).Q)
    if VERSION >= v"1.7"
        @test !fast_scalar_indexing(qr(rand(10, 10), ColumnNorm()).Q)
    else
        @test !fast_scalar_indexing(qr(rand(10, 10), Val(true)).Q)
    end
    @test !fast_scalar_indexing(lq(rand(10, 10)).Q)
    @test fast_scalar_indexing(Nothing)  # test default
end

@testset "isstructured" begin
    Sp=sparse([1,2,3],[1,2,3],[1,2,3])
    @test has_sparsestruct(Sp)
    rowind,colind=findstructralnz(Sp)
    @test [Tridiagonal([1,2,3],[1,2,3,4],[4,5,6])[rowind[i],colind[i]] for i in 1:length(rowind)]==[1,2,3]
end

@testset "restructure" begin
    x = rand(Float32, 2, 2)
    y = rand(4)
    yr = ArrayInterfaceCore.restructure(x, y)
    @test yr isa Matrix{Float64}
    @test size(yr) == (2,2)
    @test vec(yr) == vec(y)

    @testset "views" begin
        x = @view rand(4)[1:2]
        y = rand(2)
        yr = ArrayInterfaceCore.restructure(x, y)
        @test yr isa Vector{Float64}
        @test size(yr) == (2,)
        @test yr == y

        x = @view rand(4,4)[1:2,1:2]
        y = rand(2,2)
        yr = ArrayInterfaceCore.restructure(x, y)
        @test yr isa Matrix{Float64}
        @test size(yr) == (2,2)
        @test yr == y

        x = @view rand(4,4)[1]
        y = @view rand(2,2)[1]
        yr = ArrayInterfaceCore.restructure(x, y)
        @test yr isa Array{Float64,0}
        @test size(yr) == ()
        @test yr == y
    end
end

@testset "isstructured" begin
    @test !@inferred(ArrayInterfaceCore.isstructured(Matrix{Int}))  # default
    @test @inferred(ArrayInterfaceCore.isstructured(Hermitian{Complex{Int64}, Matrix{Complex{Int64}}}))
    @test @inferred(ArrayInterfaceCore.isstructured(Symmetric{Int,Matrix{Int}}))
    @test @inferred(ArrayInterfaceCore.isstructured(LowerTriangular{Int,Matrix{Int}}))
    @test @inferred(ArrayInterfaceCore.isstructured(UpperTriangular{Int,Matrix{Int}}))
    @test @inferred(ArrayInterfaceCore.isstructured(Diagonal{Int64, Vector{Int64}}))
    @test @inferred(ArrayInterfaceCore.isstructured(Bidiagonal{Int64, Vector{Int64}}))
    @test @inferred(ArrayInterfaceCore.isstructured(Tridiagonal([1,2,3],[1,2,3,4],[4,5,6])))
    @test @inferred(ArrayInterfaceCore.isstructured(SymTridiagonal{Int64, Vector{Int64}}))
end

@testset "ArrayInterfaceCore.issingular" begin
    for T in [Float64, ComplexF64]
        R = randn(MersenneTwister(2), T, 5, 5)
        S = Symmetric(R)
        L = UpperTriangular(R)
        U = LowerTriangular(R)
        @test all(!ArrayInterfaceCore.issingular, [R, S, L, U, U'])
        R[:, 2] .= 0
        @test all(ArrayInterfaceCore.issingular, [R, L, U, U'])
        @test !ArrayInterfaceCore.issingular(S)
        R[2, :] .= 0
        @test ArrayInterfaceCore.issingular(S)
        @test all(!ArrayInterfaceCore.issingular, [UnitLowerTriangular(R), UnitUpperTriangular(R), UnitUpperTriangular(R)'])
    end
    @test !@inferred(ArrayInterfaceCore.issingular(Diagonal([1,2,3,4])))
    @test @inferred(ArrayInterfaceCore.issingular(UniformScaling(0)))
    @test !@inferred(ArrayInterfaceCore.issingular(Bidiagonal([1,2,3,4], [7,8,9], :U)))
    @test !@inferred(ArrayInterfaceCore.issingular(SymTridiagonal([1,2,3,4],[5,6,7])))
    @test !@inferred(ArrayInterfaceCore.issingular(Tridiagonal([1,2,3],[1,2,3,4],[4,5,6])))
end

@testset "setindex" begin
    @testset "$(typeof(x))" for x in [
        zeros(3),
        falses(3),
        spzeros(3),
    ]
        y = setindex(x, true, 1)
        @test iszero(x)  # x is not mutated
        @test y[1] == true
        @test iszero(x[CartesianIndices(size(x)) .== [CartesianIndex(1)]])

        y2 = setindex(x, one.(x), :)
        @test iszero(x)
        @test all(isone, y2)
    end

    @testset "$(typeof(x))" for x in [
        zeros(3, 3),
        falses(3, 3),
        spzeros(3, 3),
    ]
        y = setindex(x, true, 1, 1)
        @test iszero(x)  # x is not mutated
        @test y[1, 1] == true
        @test iszero(x[CartesianIndices(size(x)) .== [CartesianIndex(1, 1)]])

        y2 = setindex(x, one.(x), :, :)
        @test iszero(x)
        @test all(isone, y2)
    end

    @testset "$(typeof(x))" for x in [
        zeros(3, 3, 3),
        falses(3, 3, 3),
    ]
        y = setindex(x, true, 1, 1, 1)
        @test iszero(x)  # x is not mutated
        @test y[1, 1, 1] == true
        @test iszero(x[CartesianIndices(size(x)) .== [CartesianIndex(1, 1, 1)]])

        y2 = setindex(x, one.(x), :, :, :)
        @test iszero(x)
        @test all(isone, y2)
    end
end

@testset "Sparsity Structure" begin
    D=Diagonal([1,2,3,4])
    @test has_sparsestruct(D)
    rowind,colind=findstructralnz(D)
    @test [D[rowind[i],colind[i]] for i in 1:length(rowind)]==[1,2,3,4]
    @test length(rowind)==4
    @test length(rowind)==length(colind)

    Bu = Bidiagonal([1,2,3,4], [7,8,9], :U)
    @test has_sparsestruct(Bu)
    rowind,colind=findstructralnz(Bu)
    @test [Bu[rowind[i],colind[i]] for i in 1:length(rowind)]==[1,7,2,8,3,9,4]
    Bl = Bidiagonal([1,2,3,4], [7,8,9], :L)
    @test has_sparsestruct(Bl)
    rowind,colind=findstructralnz(Bl)
    @test [Bl[rowind[i],colind[i]] for i in 1:length(rowind)]==[1,7,2,8,3,9,4]

    Tri=Tridiagonal([1,2,3],[1,2,3,4],[4,5,6])
    @test has_sparsestruct(Tri)
    rowind,colind=findstructralnz(Tri)
    @test [Tri[rowind[i],colind[i]] for i in 1:length(rowind)]==[1,2,3,4,4,5,6,1,2,3]

    STri=SymTridiagonal([1,2,3,4],[5,6,7])
    @test has_sparsestruct(STri)
    rowind,colind=findstructralnz(STri)
    @test [STri[rowind[i],colind[i]] for i in 1:length(rowind)]==[1,2,3,4,5,6,7,5,6,7]
end

@testset "known values" begin
    CI = CartesianIndices((2, 2))

    @test isnothing(@inferred(ArrayInterfaceCore.known_first(typeof(1:4))))
    @test isone(@inferred(ArrayInterfaceCore.known_first(Base.OneTo(4))))
    @test isone(@inferred(ArrayInterfaceCore.known_first(typeof(Base.OneTo(4)))))
    @test @inferred(ArrayInterfaceCore.known_first(typeof(CI))) == CartesianIndex(1, 1)
    @test @inferred(ArrayInterfaceCore.known_first(typeof(CI))) == CartesianIndex(1, 1)

    @test isnothing(@inferred(ArrayInterfaceCore.known_last(1:4)))
    @test isnothing(@inferred(ArrayInterfaceCore.known_last(typeof(1:4))))
    @test @inferred(ArrayInterfaceCore.known_last(typeof(CI))) === nothing

    @test isnothing(@inferred(ArrayInterfaceCore.known_step(typeof(1:0.2:4))))
    @test isone(@inferred(ArrayInterfaceCore.known_step(1:4)))
    @test isone(@inferred(ArrayInterfaceCore.known_step(typeof(1:4))))
    @test isone(@inferred(ArrayInterfaceCore.known_step(typeof(Base.Slice(1:4)))))
    @test isone(@inferred(ArrayInterfaceCore.known_step(typeof(view(1:4, 1:2)))))
end

@testset "ndims_index" begin
    @test @inferred(ArrayInterfaceCore.ndims_index(CartesianIndices(()))) == 1
    @test @inferred(ArrayInterfaceCore.ndims_index(trues(2, 2))) == 2
    @test @inferred(ArrayInterfaceCore.ndims_index(CartesianIndex(2,2))) == 2
    @test @inferred(ArrayInterfaceCore.ndims_index(1)) == 1
end

@testset "ndims_shape" begin
    @test @inferred(ArrayInterfaceCore.ndims_shape(1)) === 0
    @test @inferred(ArrayInterfaceCore.ndims_shape(:)) === 1
    @test @inferred(ArrayInterfaceCore.ndims_shape(CartesianIndex(1, 2))) === (0, 0)
    @test @inferred(ArrayInterfaceCore.ndims_shape(CartesianIndices((2,2)))) === (1, 1)
    @test @inferred(ArrayInterfaceCore.ndims_shape([1 1])) === 2
end

@testset "indices_do_not_alias" begin
  @test ArrayInterfaceCore.instances_do_not_alias(Float64)
  @test !ArrayInterfaceCore.instances_do_not_alias(Matrix{Float64})
  @test ArrayInterfaceCore.indices_do_not_alias(Matrix{Float64})
  @test !ArrayInterfaceCore.indices_do_not_alias(BitMatrix)
  @test !ArrayInterfaceCore.indices_do_not_alias(Matrix{Matrix{Float64}})
  @test ArrayInterfaceCore.indices_do_not_alias(Adjoint{Float64,Matrix{Float64}})
  @test ArrayInterfaceCore.indices_do_not_alias(Transpose{Float64,Matrix{Float64}})
  @test ArrayInterfaceCore.indices_do_not_alias(typeof(view(rand(4,4)', 2:3, 1:2)))
  @test ArrayInterfaceCore.indices_do_not_alias(typeof(view(rand(4,4,4), CartesianIndex(1,2), 2:3)))
  @test ArrayInterfaceCore.indices_do_not_alias(typeof(view(rand(4,4)', 1:2, 2)))
  @test !ArrayInterfaceCore.indices_do_not_alias(typeof(view(rand(7),ones(Int,7))))
  @test !ArrayInterfaceCore.indices_do_not_alias(Adjoint{Matrix{Float64},Matrix{Matrix{Float64}}})
  @test !ArrayInterfaceCore.indices_do_not_alias(Transpose{Matrix{Float64},Matrix{Matrix{Float64}}})
  @test !ArrayInterfaceCore.indices_do_not_alias(typeof(view(fill(rand(4,4),4,4)', 2:3, 1:2)))
  @test !ArrayInterfaceCore.indices_do_not_alias(typeof(view(rand(4,4)', StepRangeLen(1,0,5), 1:2)))
end


text:
# fisher.jl
# Fisher's exact test
#
# Copyright (C) 2013   Simon Kornblith
#
# Permission is hereby granted, free of charge, to any person obtaining
# a copy of this software and associated documentation files (the
# "Software"), to deal in the Software without restriction, including
# without limitation the rights to use, copy, modify, merge, publish,
# distribute, sublicense, and/or sell copies of the Software, and to
# permit persons to whom the Software is furnished to do so, subject to
# the following conditions:
#
# The above copyright notice and this permission notice shall be
# included in all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
# NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
# LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
# OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
# WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

export FisherExactTest

"""
    FisherExactTest(a::Integer, b::Integer, c::Integer, d::Integer)

Perform Fisher's exact test of the null hypothesis that the success probabilities ``a/c``
and ``b/d`` are equal, that is the odds ratio ``(a/c) / (b/d)`` is one, against the
alternative hypothesis that they are not equal.

See [`pvalue(::FisherExactTest)`](@ref) and [`confint(::FisherExactTest)`](@ref) for details
about the computation of the default p-value and confidence interval, respectively.

The contingency table is structured as:

| -  | X1 | X2 |
|:--:|:--:|:--:|
|*Y1*| a  | b  |
|*Y2*| c  | d  |

!!! note
    The `show` function output contains the conditional maximum likelihood estimate of the
    odds ratio rather than the sample odds ratio; it maximizes the likelihood given by
    Fisher's non-central hypergeometric distribution.

Implements: [`pvalue(::FisherExactTest)`](@ref), [`confint(::FisherExactTest)`](@ref)

# References

  * Fay, M.P., Supplementary material to "Confidence intervals that match Fisher’s exact or
    Blaker’s exact tests". Biostatistics, Volume 11, Issue 2, 1 April 2010, Pages 373–374,
    [link](https://doi.org/10.1093/biostatistics/kxp050)
"""
struct FisherExactTest <: HypothesisTest
    # Format:
    # X1  X2
    # Y1  a  b
    # Y2  c  d
    a::Int
    b::Int
    c::Int
    d::Int

    # conditional maximum likehood estimate of odd ratio
    ω::Float64

    function FisherExactTest(a::Int, b::Int, c::Int, d::Int)
        ω = cond_mle_odds_ratio(a, b, c, d)
        new(a, b, c, d, ω)
    end
end

testname(::FisherExactTest) = "Fisher's exact test"
population_param_of_interest(x::FisherExactTest) = ("Odds ratio", 1.0, x.ω) # parameter of interest: name, value under h0, point estimate
default_tail(test::FisherExactTest) = :both

# The sizing argument to print_matrix was removed during the 0.5 dev period
function _print_matrix(io::IO, X::AbstractVecOrMat, pre::AbstractString)
    Base.print_matrix(io, X, pre)
end

function show_params(io::IO, x::FisherExactTest, ident="")
    println(io, ident, "contingency table:")
    _print_matrix(io, [x.a x.b; x.c x.d], repeat(ident, 2))
    println(io)
end

# DOC: for tail=:both there exist multiple ``method``s for computing a pvalue and the corresponding ci.
"""
    pvalue(x::FisherExactTest; tail = :both, method = :central)

Compute the p-value for a given Fisher exact test.

The one-sided p-values are based on Fisher's non-central hypergeometric distribution
``f_ω(i)`` with odds ratio ``ω``:
```math
    \\begin{align*}
        p_ω^{(\\text{left})} &=\\sum_{i ≤ a} f_ω(i)\\\\
        p_ω^{(\\text{right})} &=\\sum_{i ≥ a} f_ω(i)
    \\end{align*}
```
For `tail = :both`, possible values for `method` are:

  - `:central` (default): Central interval, i.e. the p-value is two times the minimum of the
    one-sided p-values.
  - `:minlike`: Minimum likelihood interval, i.e. the p-value is computed by summing all
    tables with the same marginals that are equally or less probable:
    ```math
        p_ω = \\sum_{f_ω(i)≤ f_ω(a)} f_ω(i)
    ```

# References

  * Gibbons, J.D., Pratt, J.W., P-values: Interpretation and Methodology, American
    Statistican, 29(1):20-25, 1975.
  * Fay, M.P., Supplementary material to "Confidence intervals that match Fisher’s exact or
    Blaker’s exact tests". Biostatistics, Volume 11, Issue 2, 1 April 2010, Pages 373–374,
    [link](https://doi.org/10.1093/biostatistics/kxp050)
"""
function pvalue(x::FisherExactTest; tail=:both, method=:central)
    if tail == :both && method != :central
        if method == :minlike
            p = pvalue_both_minlike(x)
        else
            throw(ArgumentError("method=$(method) is not implemented yet"))
        end
    else
        p = pvalue(Hypergeometric(x.a + x.b, x.c + x.d, x.a + x.c), x.a, tail=tail)
    end
    p = max(min(p, 1.0), 0.0)

    return p
end

function pvalue_both_minlike(x::FisherExactTest, ω::Float64=1.0)
    a, b, c, d = reorder(x.a, x.b, x.c, x.d)

    dist = FisherNoncentralHypergeometric(a+b, c+d, a+c, ω)

    p = pdf(dist, a)
    v = nextfloat(p)
    if a != 0
        p += cdf(dist, a-1)
    end

    # Add p-values of all tables in other tail equally or less probable
    for i = a+c:-1:a+1
        curp = pdf(dist, i)
        if curp > v
            break
        end
        p += curp
    end
    p
end

# confidence interval by inversion of p-value
"""
    confint(x::FisherExactTest, alpha::Float64=0.05; tail=:both, method=:central)

Compute a confidence interval with coverage 1 - `alpha`. One-sided intervals are based on
Fisher's non-central hypergeometric distribution. For `tail = :both`, the only
`method` implemented yet is the central interval (`:central`).

!!! note
    Since the p-value is not necessarily unimodal, the corresponding confidence region might
    not be an interval.

# References

  * Gibbons, J.D, Pratt, J.W. P-values: Interpretation and Methodology, American
    Statistican, 29(1):20-25, 1975.
  * Fay, M.P., Supplementary material to "Confidence intervals that match Fisher’s exact or
    Blaker’s exact tests". Biostatistics, Volume 11, Issue 2, 1 April 2010, Pages 373–374,
    [link](https://doi.org/10.1093/biostatistics/kxp050)
"""
function StatsBase.confint(x::FisherExactTest, alpha::Float64=0.05; tail=:both, method=:central)
    check_alpha(alpha)
    dist(ω) = FisherNoncentralHypergeometric(x.a+x.b, x.c+x.d, x.a+x.c, ω)
    obj(ω) = pvalue(dist(ω), x.a, tail=tail) - alpha

    if tail == :left # upper bound
        if (x.a == maximum(dist(1.0)))
            (0.0, Inf)
        else
            (0.0, fzero(obj, find_brackets(obj)...))
        end
    elseif tail == :right # lower bound
        if (x.a == minimum(dist(1.0)))
            (0.0, Inf)
        else
            (fzero(obj, find_brackets(obj)...), Inf)
        end
    elseif tail == :both
        if method == :central
            (StatsBase.confint(x, alpha/2; tail=:right)[1], StatsBase.confint(x, alpha/2; tail=:left)[2])
        else
            throw(ArgumentError("method=$(method) is not implemented yet"))
        end
    else
        throw(ArgumentError("tail=$(tail) is invalid"))
    end
end

## helpers

function reorder(a,b,c,d)
    if a + c > b + d
        a, b, c, d = b, a, d, c
    end
    if a/c > b/d
        a, b, c, d = c, d, a, b
    end
    (a, b, c, d)
end

# find values x_lower and x_upper, s.t. f(x_lower) < 0 and f(x_upper) > 0 or vice versa
function find_brackets(f::Function, x_init::Float64=1.0)
    f_init = f(x_init)

    if f_init > f(x_init + 1.0)
        find_brackets(x -> -f(x), x_init)
    else
        x_upper = x_lower = x_init
        if f_init > 0.0
            while x_lower > eps(0.0) && f(x_lower) > 0.0
                x_lower /= 2
            end
        else
            while f(x_upper) < 0.0
                x_upper *= 2
            end
        end
        (x_lower, x_upper)
    end
end

# find odds ratio ω that maximizes the (conditional) likelihood;
# since the mode and mean of Fisher's Noncentral Hypergeometric distribution
# coincide, this is equivalent to find ω, s.t., mean(dist(ω)) = a
function cond_mle_odds_ratio(a::Int, b::Int, c::Int, d::Int)
    dist(ω) = FisherNoncentralHypergeometric(a+b, c+d, a+c, ω)

    if (a == minimum(dist(1.0)))
        0.0
    elseif a == maximum(dist(1.0))
        Inf
    else
        obj(ω) = mean(dist(ω))-a
        fzero(obj, find_brackets(obj)...)
    end
end


text:
using Test

import MathOptInterface
const MOI = MathOptInterface

import ComplexOptInterface
const COI = ComplexOptInterface

"""
    projection_test(optimizer, config)

Test computation of the projection of a hermitian matrix to the cone
of hermitian positive semidefinite matrices.

[ 1    -1+im]   [ 1-im  ?]   [√3    ]   [ 1+im  -1+√3]
[-1-im -1   ] = [-1+√3  ?] * [   -√3] * [  ?      ?  ] / (6 - 2√3)
So it's projection to the Hermitian PSD cone is
                 [1-im]  [ 1+im  1-√3]
√3 / (6 - 2√3) * [1-√3]
which is
                 [1-im]  [ 1+im  1-√3]
(1 + √3) / 4   * [1-√3]
which is
[(1+√3)/2  -1/2+im/2]
[-1/2+im/2 (-1+√3)/2]
"""
function projection_test(optimizer, config)
    atol = config.atol
    rtol = config.rtol

    MOI.empty!(optimizer)
    set = COI.HermitianPositiveSemidefiniteConeTriangle(2)
    x, cx = MOI.add_constrained_variables(optimizer, set)
    x11 = x[1:3]
    x12 = x[4]
    t = MOI.add_variable(optimizer)
    ft = MOI.SingleVariable(t)
    MOI.add_constraint(optimizer, MOI.Utilities.operate(vcat, Float64, ft, x[1] - 1.0, √2 * (x[2] + 1.0), x[3] + 1.0, √2 * (x[4] - 1.0)),
                       MOI.SecondOrderCone(5))
    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    MOI.set(optimizer, MOI.ObjectiveFunction{typeof(ft)}(), ft)
    MOI.optimize!(optimizer)
    primal = [(1 + √3) / 2, -1/2, (-1 + √3)/2, 1/2]
    dual = [(3 - √3) / 6, 0.2886751198, (3 + √3) / 6, -0.2886751197]
    @test MOI.Utilities.set_dot(primal, dual, set) ≈ 0 atol=atol rtol=rtol
    @test MOI.get(optimizer, MOI.VariablePrimal(), x) ≈ primal atol=atol rtol=rtol
    @test MOI.get(optimizer, MOI.ConstraintPrimal(), cx) ≈ primal atol=atol rtol=rtol
    @test MOI.get(optimizer, MOI.ConstraintDual(), cx) ≈ dual atol=atol rtol=rtol
end

function hermitian_psd_test(optimizer, config)
    atol = config.atol
    rtol = config.rtol

    MOI.empty!(optimizer)
    set = COI.HermitianPositiveSemidefiniteConeTriangle(3)
    x, cx = MOI.add_constrained_variables(optimizer, set)
    MOI.add_constraints(optimizer, x, MOI.EqualTo.([1.0, 0.0, 1.0, 0.0, 0.0, 0.0, -1.0, 0.0, 0.0]))
    MOI.optimize!(optimizer)
    primal = [1.0, 0.0, 1.0, 0.0, 0.0, 0.0, -1.0, 0.0, 0.0]
    @test MOI.get(optimizer, MOI.VariablePrimal(), x) ≈ primal atol=atol rtol=rtol
    @test MOI.get(optimizer, MOI.ConstraintPrimal(), cx) ≈ primal atol=atol rtol=rtol
    @test MOI.get(optimizer, MOI.ConstraintDual(), cx) ≈ zeros(9) atol=atol rtol=rtol
end

function zero_1_test(optimizer, config)
    atol = config.atol
    rtol = config.rtol

    MOI.empty!(optimizer)
    x, cx = MOI.add_constrained_variables(optimizer, MOI.Nonnegatives(2))
    func = (1.0 + 2.0im) * x[1] + (1.0 - 1.0im) * x[2] + (-1.0 + -1.0im)
    c = MOI.add_constraint(
        optimizer,
        MOI.Utilities.operate(vcat, Complex{Float64}, func),
        MOI.Zeros(1)
    )
    MOI.optimize!(optimizer)
    @test MOI.get(optimizer, MOI.TerminationStatus()) == MOI.OPTIMAL
    @test MOI.get(optimizer, MOI.VariablePrimal(), x) ≈ [2/3, 1/3] atol=atol rtol=rtol
    @test MOI.get(optimizer, MOI.ConstraintPrimal(), cx) ≈ [2/3, 1/3] atol=atol rtol=rtol
    @test MOI.get(optimizer, MOI.ConstraintDual(), cx) ≈ zeros(2) atol=atol rtol=rtol
    @test MOI.get(optimizer, MOI.ConstraintPrimal(), c) ≈ [0.0 + 0.0im] atol=atol rtol=rtol
    @test MOI.get(optimizer, MOI.ConstraintDual(), c) ≈ [0.0 + 0.0im] atol=atol rtol=rtol
end

function zero_2_test(optimizer, config)
    atol = config.atol
    rtol = config.rtol

    MOI.empty!(optimizer)
    x, cx = MOI.add_constrained_variables(optimizer, MOI.Nonnegatives(1))
    func = (1.0 + 0.0im) * x[1] + 1.0im * x[1] - 2.0im - (1.0 + 0.0im) * x[1]
    c = MOI.add_constraint(
        optimizer,
        MOI.Utilities.operate(vcat, Complex{Float64}, func),
        MOI.Zeros(1)
    )
    MOI.optimize!(optimizer)
    @test MOI.get(optimizer, MOI.TerminationStatus()) == MOI.OPTIMAL
    @test MOI.get(optimizer, MOI.VariablePrimal(), x) ≈ [2.0] atol=atol rtol=rtol
    @test MOI.get(optimizer, MOI.ConstraintPrimal(), cx) ≈ [2.0] atol=atol rtol=rtol
    @test MOI.get(optimizer, MOI.ConstraintDual(), cx) ≈ zeros(1) atol=atol rtol=rtol
    @test MOI.get(optimizer, MOI.ConstraintPrimal(), c) ≈ [0.0 + 0.0im] atol=atol rtol=rtol
    @test MOI.get(optimizer, MOI.ConstraintDual(), c) ≈ [0.0 + 0.0im] atol=atol rtol=rtol
end

import CSDP
@testset "CSDP" begin
    config = MOI.Test.Config(atol=1e-4, rtol=1e-4)
    bridged = MOI.instantiate(CSDP.Optimizer, with_bridge_type=Float64)
    MOI.Bridges.add_bridge(bridged, COI.Bridges.Variable.HermitianToSymmetricPSDBridge{Float64})
    projection_test(bridged, config)
    hermitian_psd_test(bridged, config)
    bridged = MOI.Bridges.LazyBridgeOptimizer(MOI.Utilities.CachingOptimizer(MOI.Utilities.Model{Float64}(), CSDP.Optimizer()))
    MOI.Bridges.add_bridge(bridged, MOI.Bridges.Constraint.ScalarizeBridge{Float64})
    MOI.Bridges.add_bridge(bridged, COI.Bridges.Constraint.SplitZeroBridge{Float64})
    zero_1_test(bridged, config)
    cis = MOI.get(bridged.model, MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}}())
    @test length(cis) == 2
    zero_2_test(bridged, config)
    cis = MOI.get(bridged.model, MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}}())
    @test length(cis) == 1
end


text:
# This file is a part of LegendHDF5IO.jl, licensed under the MIT License (MIT

function to_table(x::AbstractVector{<:RDWaveform})
    TypedTables.Table(
        t0 = first.(x.time),
        dt = step.(x.time),
        values = x.value
    )
end

function _dtt02range(dt::RealQuantity, t0::RealQuantity, values::AbstractVector)
    # TODO: Handle different units for dt and t0
    t_idxs = Int32(0):Int32(size(values,1) - 1)
    t0 .+ t_idxs .* dt
end

function from_table(tbl, ::Type{<:AbstractVector{<:RDWaveform}})
    StructArray{RDWaveform}((
        _dtt02range.(tbl.dt, tbl.t0, tbl.values),
        tbl.values
    ))
end

function LegendDataTypes.writedata(
    output::HDF5.DataFile, name::AbstractString,
    x::AbstractVector{<:RDWaveform},
    fulldatatype::DataType = typeof(x)
) where {T}
    @assert fulldatatype == typeof(x)
    writedata(output, name, to_table(x))
end

function LegendDataTypes.readdata(
    input::HDF5.DataFile, name::AbstractString,
    AT::Type{<:AbstractVector{<:RDWaveform}}
)
    tbl = readdata(input, name, TypedTables.Table{<:NamedTuple{(:t0, :dt, :values)}})
    from_table(tbl, AbstractVector{<:RDWaveform})
end


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

name = "Tk"
version = v"8.6.11" # current version number is actually 8.6.11.1

# Collection of sources required to build Tk
sources = [
    ArchiveSource("https://downloads.sourceforge.net/sourceforge/tcl/tk$(version).1-src.tar.gz",
                  "006CAB171BEECA6A968B6D617588538176F27BE232A2B334A0E96173E89909BE"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
if [[ "${target}" == *-mingw* ]]; then
    cd $WORKSPACE/srcdir/tk*/win/
else
    cd $WORKSPACE/srcdir/tk*/unix/
fi

export CFLAGS="-I${prefix}/include ${CFLAGS}"

FLAGS=(--enable-threads --disable-rpath)
if [[ "${target}" == x86_64-* ]]; then
    FLAGS+=(--enable-64bit)
fi
if [[ "${target}" == *-apple-* ]] || [[ "${target}" == -*mingw* ]]; then
    FLAGS+=(--with-x=no)
fi
if [[ "${target}" == *-apple-* ]]; then
    FLAGS+=(--enable-aqua=yes)

    # The following patch replaces the hard-coded path of Cocoa framework
    # with the actual path on our system.
    atomic_patch -p1 "${WORKSPACE}/srcdir/patches/apple_cocoa_configure.patch"
fi
if [[ "${target}" == *mingw* ]]; then
    # `windres` invocations don't get the proper tk include path; just hack it in
    atomic_patch -p2 "${WORKSPACE}/srcdir/patches/win_tk_rc_include.patch"
fi

./configure --prefix=${prefix} --host=${target} "${FLAGS[@]}"
make -j${nproc}
make install
make install-private-headers

# Install license file
install_license $WORKSPACE/srcdir/tk*/license.terms
"""

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

# The products that we will ensure are always built
products = [
    LibraryProduct(["libtk8.6", "libtk8", "tk86"], :libtk),
    ExecutableProduct(["wish8.6", "wish86"], :wish),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    BuildDependency("Xorg_xorgproto_jll"),
    Dependency("Tcl_jll"),
    Dependency("Xorg_libXft_jll"),
]

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

max_stars_repo_name: blenessy/Yggdrasil

max_issues_repo_name: blenessy/Yggdrasil

max_forks_repo_name: blenessy/Yggdrasil


text:
"""
    scaled_squared_distance(log_ell::AbstractArray{Float64, 1},
        x::AbstractArray{Float64, 2}, z::AbstractArray{Float64, 2})

Compute the scaled squared distance between `x` and `z`:
```math
r_{ij} = \\sum_{k=1}^d\\frac{(x_{ik}-z_{jk})^2}{l_k^2}
```
The gradient of this function with respect to
length scale hyperparameter(s) is returned by [`scaled_squared_distance_grad`](@ref).

# Arguments
- `x, z`: Input data, reshaped into 2-d arrays.
  `x` must have dimensions ``n \\times d``; `z` must have dimensions ``m \\times d``.
- `log_ell`: logarithm of length scale(s). Can either be an array of size one (isotropic),
  or an array of size `d` (ARD)

# Return
An ``n \\times m`` matrix of scaled squared distances
"""
function scaled_squared_distance(log_ell::AbstractArray{Float64, 1},
        x1::AbstractArray{Float64, 2}, x2::AbstractArray{Float64, 2})
    ell = exp.(-log_ell .* 2)
    if x1 === x2
        mu = mean(x1, dims=1)
    else
        n1 = size(x1, 1)
        n2 = size(x2, 1)
        mu = n1 / (n1 + n2) * mean(x1, dims=1) + n2 / (n1 + n2) * mean(x2, dims=1)
    end
    x1 = x1 .- mu               # subtract mean for stability, in case x1 or x2 are huge
    ax1 = x1 .* ell'            # a * x1
    sax1 = sum(x1 .* ax1, dims=2)    # a * x1 ^ 2
    if x1 === x2                # a shortcut if we need K(x, x) for log likelihood
        ax2 = ax1
        sax2 = sax1
    else                        # repeat computations for x2
        x2 = x2 .- mu
        ax2 = x2 .* ell'
        sax2 = sum(x2 .* ax2, dims=2)
    end
    D2 = sax1 .- 2 * x1 * ax2' .+ sax2'   # a((b-c)^2) = ab^2 - 2abc + ac^2
end

"""
    scaled_squared_distance_grad(log_ell::AbstractArray{Float64, 1},
        x::AbstractArray{Float64, 2}, z::AbstractArray{Float64, 2}, R::AbstractArray{Float64, 2})

Return the gradient of the [`scaled_squared_distance`](@ref) function with respect to
logarigthms of length scales, based on the provided direction matrix.

# Arguments
- `x, z`: Input data, reshaped into 2-d arrays.
  `x` must have dimensions ``n \\times d``; `z` must have dimensions ``m \\times d``.
- `log_ell`: logarithm of length scale(s). Can either be an array of size one (isotropic),
  or an array of size `d` (ARD)
- `R` the direction matrix, ``n \\times m``. This can be used to compute the gradient
  of a function that depends on [`scaled_squared_distance`](@ref) via the chain rule.

# Return
A vector of size `length(log_ell)`, whose ``k``'th element is equal to
```math
\\text{tr}(R \\frac{\\partial K}{\\partial l_k})
```
"""
function scaled_squared_distance_grad(log_ell::AbstractArray{Float64, 1},
        x1::AbstractArray{Float64, 2}, x2::AbstractArray{Float64, 2}, R::AbstractArray{Float64, 2})
    n1 = size(x1, 1)
    n2 = size(x2, 2)
    if x1 === x2                                        # a shortcut if x1 and x2 are the same - handy for log likelihood optimisation
        x1 = x1 .- mean(x1, dims=1)                          # subtract mean for stability if x1 is huge
        x12_2 = (sum(R, dims=2)' + sum(R, dims=1)) * (x1 .^ 2)    # R * x1 ^ 2 + R * x2 ^ 2
        x1x2 = sum((R * x1) .* x1, dims=1)                   # R * x1 * x2
        rk = x12_2' - 2 * x1x2'                         # R * (x1 - x2) ^ 2 = R * x1 ^ 2 - 2 * R * x1 * x2 + R * x2 ^ 2
    else
        mu = n1 / (n1 + n2) * mean(x1, dims=1) + n2 / (n1 + n2) * mean(x2, dims=1)
        x1 = x1 .- mu                                   # subtract mean for stability if x1 or x2 is huge
        x2 = x2 .- mu
        x1_2 = sum(R, dims=2)' * (x1 .^ 2)                   # R * x1 ^ 2
        x2_2 = sum(R, dims=1) * (x2 .^ 2)                    # R * x2 ^ 2
        x1x2 = sum((R * x2) .* x1, dims=1)                   # R * x1 * x2
        rk = x1_2' - 2 * x1x2' + x2_2'                  # R * (x1 - x2) ^ 2 = R * x1 ^ 2 - 2 * R * x1 * x2 + R * x2 ^ 2
    end
    ell = exp.(-2 * log_ell)
    return reshape(-2 * ell .* rk, length(ell))
end


text:
using FastCap

cl = 500E-6 # side length of the chip
mt = 0.1E-6 # metal thickness
st = 50E-6 # substrate thickness

cw = 10E-6 # center cconuctor width
gg = 6E-6 # ground gap width

il = 249E-6 # transmon island length
iw = 10E-6 # transmon island width (sans digits)
ig = 3E-6 # transmon island gap to ground
ii = 6E-6 # transmon island to island gap

dl = 14E-6 # digit length
dw = 6E-6 # digit width
ds = 12E-6 # digit spacing

ks = 11.9 # substrate relative permitivitty

# Compute some distances
tl = il+2ig # Transmon length, including ground gap
tw = 2(ig+iw)+ii+dl # Transmon width, including ground gap
s1 = (cl-tl)/2 # length of smaller ground section on transmon side of tl
tf = (cl+tl)/2 # y distance to far side of transmon (up to ground plane)
cn = (cl-cw)/2 # y distance to near side of center conductor
cf = (cl+cw)/2 # y distance to far side of center conductor

m = Model(recombineall=true)

# Center conductor
cen = add_rectangle!(m,0,cn,cl,cn+cw,lc=gg,transfinite=[0,(8,"Bump",0.1),0,(8,"Bump",0.1)])
add_physicalsurf!(m,"center_bot",cen)
ext = extrude!(m,cen,(0,0,mt),transfinite=(6,"Bump",0.1))
add_physicalsurf!(m,"center_top",ext)

# Island
isl = add_rectangle!(m,s1+ig,cf+ig,s1+ig+il,cf+ig+tw-2ig,lc=gg,transfinite=[0,(8,"Bump",0.1),0,(8,"Bump",0.1)])
add_physicalsurf!(m,"island_bot",isl)
ext = extrude!(m,isl,(0,0,mt),transfinite=(6,"Bump",0.1))
add_physicalsurf!(m,"island_top",ext)

# Substrate around island and centerline, including holes for island and center cond
xs = [0-gg, cl+gg,cl+gg,s1+tl,s1+tl,s1,   s1,   0-gg]
ys = [cn-gg,cn-gg,cf+gg,cf+gg,cf+tw,cf+tw,cf+gg,cf+gg]
sub = add_polygon!(m,xs,ys,holes=[cen.lineloops[1],isl.lineloops[1]],lc=gg)

# Ground plane around substrate
gnd = add_rectangle!(m,-10gg,-10gg,cl+10gg,cl+10gg,holes=sub.lineloops[1],lc=4gg)
add_physicalsurf!(m,"ground_bot",gnd)
ext = extrude!(m,gnd,(0,0,mt),transfinite=(6,"Bump",0.1))
add_physicalsurf!(m,"ground_top",ext)

# Make temporary surface for substrate
tmp = FastCap.add_surface!(m,gnd.lineloops[1])
ext = extrude!(m,tmp,(0,0,-st),transfinite=(6,"Progression",2))
deleteat!(m.surfaces,length(m.surfaces)-5)
add_physicalsurf!(m,"substrate",[sub;ext])

# Mesh
open(f->println(f,FastCap.print_geo_file(m)),"examples\\transmon2\\transmon2.geo","w")
readstring(`gmsh examples\\transmon2\\transmon2.geo -2 -algo front2d -o examples\\transmon2\\transmon2.msh`)
spawn(`gmsh examples\\transmon2\\transmon2.geo examples\\transmon2\\transmon2.msh`)
msh = FastCap.read_msh("examples\\transmon2\\transmon2.msh")
@time shapes = FastCap.get_shapes(msh)

geoms = [
    ConductorGroup("center",[Conductor("center_bot",ks,shapes[1]),Conductor("center_top",1.,shapes[2])]),
    ConductorGroup("island",[Conductor("island_bot",ks,shapes[3]),Conductor("island_top",1.,shapes[4])]),
    ConductorGroup("gnd",[Conductor("ground_bot",ks,shapes[5]),Conductor("ground_top",1.,shapes[6])]),
    Dielectric("substrate",1.,ks,(cl/2,cl/2,-st/2),shapes[7])
    ]

FastCap.writefilestodir(geoms,"examples\\transmon2")

#res = fastcapwr(geoms)
res = fastcap2(geoms)

max_stars_repo_name: jebej/FastCap

max_issues_repo_name: jebej/FastCap

max_forks_repo_name: jebej/FastCap


text:
"""
    Circumsphere(simplex)

Compute the circumsphere of a simplex (a matrix whose columns represent the vertices of the
simplex). Returns a column vector where the first entry is the radius of the circumsphere
and the remaining entries represent the centroid.

"""
function Circumsphere(simplex)
    n::Int = size(simplex, 1) # Dimension of the space the simplex lives in
    centroid = simplex * ones(n + 1, 1) / (n + 1)
    centroidmatrix = repeat(centroid, 1, n + 1)

    radius::Float64 = sqrt(maximum(ones(1, n) * ((simplex - centroidmatrix).^2)))

    return [radius; centroid]
end


text:
using DataFrames
using Distributions
using StatsBase
using ATools

"""
	true_int(b1::Hit, b2::Hit, df1::DataFrame, df2::DataFrame)

Return the true positions and interaction energies for
each hemisphers's interaction.
"""
function true_int(b1::Hit, df1::DataFrame, df2::DataFrame)
	# distance between baricenter and true
	d1 = dxyz([b1.x, b1.y, b1.z], [df1.x[1], df1.y[1], df1.z[1]])
	d2 = dxyz([b1.x, b1.y, b1.z], [df2.x[1], df2.y[1], df2.z[1]])

	if d2 < d1
		xt2 = [df1.x[1], df1.y[1], df1.z[1]]
		Et2 = df1.pre_KE[1]
		xt1 = [df2.x[1], df2.y[1], df2.z[1]]
		Et1 = df2.pre_KE[1]
	else
		xt1 = [df1.x[1], df1.y[1], df1.z[1]]
		Et1 = df1.pre_KE[1]
		xt2 = [df2.x[1], df2.y[1], df2.z[1]]
		Et2 = df2.pre_KE[1]
	end
	return xt1, Et1, xt2, Et2
end

"""
	check_valid_hits(hitdf::DataFrame)

Check the validity of the event hitmaps for continued processing
"""
function check_valid_hits(hitdf::Union{DataFrame, Nothing}, event::Int64)
	if hitdf === nothing
		@warn "Warning, hidtf evaluates to nothing for event = $event"
		return false
	end

	if nrow(hitdf) < 2
		@warn "Warning, hidtf is <2 for event = $event"
		return false
	end
	return true
end

"""
	check_valid_charge(h1::DataFrame, h2::DataFrame, dc::DetConf)

Check if the hemispheres have detected charge above threshold
"""
function check_valid_charge(h1::DataFrame, h2::DataFrame, dc::DetConf)
	q1 = sum(h1.q)
	q2 = sum(h2.q)
	@info " total charge: q1 = $(q1), q2 = $(q2)"
	if q1 < dc.qmin || q1 > dc.qmax
		@info "Warning, q1 is $(q1) for event"
		return false
	end

	if q2 < dc.qmin || q2 > dc.qmax
		@info "Warning, q2 is $(q2) for event"
		return false
	end
	return true
end

"""
	recovent(event       ::Integer,
	dc          		 ::DetConf,
	df1         		 ::DataFrame,
	df2         		 ::DataFrame,
	primaries   		 ::SubDataFrame,
	sensor_xyz  		 ::DataFrame,
	waveform    		 ::SubDataFrame,
	lor_algo    		 ::Function)

Return a dictionary with the variables characterising the event.
"""
function recovent(event     ::Integer,
		  dc        ::DetConf,
		  df1       ::DataFrame,
		  df2       ::DataFrame,
		  primaries ::SubDataFrame,
		  sensor_xyz::DataFrame,
		  waveform  ::SubDataFrame,
		  lor_algo  ::Function)
	# hit dataframe
	hitdf = recohits(event,sensor_xyz, waveform, dc.ecut, dc.pde, dc.sigma_tof)

	if !check_valid_hits(hitdf, event)
		return nothing
	end

	# reconstruct (x,y) : barycenter
	b1, b2, hq1df, hq2df = lor_algo(hitdf)
	@info " barycenters" b1 b2

	if !check_valid_charge(hq1df, hq2df, dc)
		return nothing
	end

	# find true position (and correlate with barycenter)
	xt1, Et1, xt2, Et2 = true_int(b1, df1, df2)
	@info " True position in hemisphere 1" xt1
	@info " True position in hemisphere 2" xt2

	## Weights and phi positions
	int1_weights = FrequencyWeights(hq1df.q)
	phi1_values  = fphi(hq1df)
	int2_weights = FrequencyWeights(hq2df.q)
	phi2_values  = fphi(hq2df)

	true_rad1 = rxy(xt1[1], xt1[2])
	true_rad2 = rxy(xt2[1], xt2[2])
	# New (x,y) positions estimated from r1, r2
	xyz1 = radial_correction(b1, true_rad1)
	xyz2 = radial_correction(b2, true_rad2)

	ntof1 = min(dc.ntof, nrow(hq1df))
	ntof2 = min(dc.ntof, nrow(hq2df))

	# sort reco times in ascending order
	t1s = sort(hq1df.trmin)
	t2s = sort(hq2df.trmin)

	tmin1 = minimum(hq1df.tmin)
	tmin2 = minimum(hq2df.tmin)
	ht1   = select_by_column_value(hq1df, "tmin", tmin1)
	ht2   = select_by_column_value(hq2df, "tmin", tmin2)

	@info " New (x,y,z) positions estimated from r1, r2 & r1q, r2q"
	@info " from r1:  x1 = $(xyz1[1]), y1=$(xyz1[2]), z1=$(xyz1[3])"
	@info " from r2:  x2 = $(xyz2[1]), y1=$(xyz2[2]), z1=$(xyz2[3])"

	@info " hit dataframe: size = $size(hitdf)"
	@debug first(hitdf, 5)

	reduced_event = ATools.EventParameters(
	event_id=event,
	phot1 = df1.process_id[1] == 1,
	phot2 = df2.process_id[1] == 1,
	# Primary particle origins
	xs = primaries.x[1],
	ys = primaries.y[1],
	zs = primaries.z[1],
	ux = primaries.vx[1],
	uy = primaries.vy[1],
	uz = primaries.vz[1],
	# barycenter on SiPM plane
	xr1    = b1.x,
	yr1    = b1.y,
	zr1    = b1.z,
	nsipm1 = nrow(hq1df),
	xr2    = b2.x,
	yr2    = b2.y,
	zr2    = b2.z,
	nsipm2 = nrow(hq1df),
	# total charge
	q1 = sum(hq1df.q),
	q2 = sum(hq2df.q),
	# Gamma energy at interaction
	E1 = Et1,
	E2 = Et2,
	# True interaction positions
	xt1 = xt1[1],
	yt1 = xt1[2],
	zt1 = xt1[3],
	t1  = tmin1,
	tr1 = minimum(hq1df.trmin),
	xt2 = xt2[1],
	yt2 = xt2[2],
	zt2 = xt2[3],
	t2  = tmin2,
	tr2 = minimum(hq2df.trmin),
	# find r1 and r2 (from True info)
	r1 = true_rad1,
	r2 = true_rad2,
	# Compute phistd and zstd of detected light
	phistd1  = std(phi1_values, int1_weights, corrected=true),
	zstd1     = std(hq1df.z   , int1_weights, corrected=true),
	widz1     = maximum(hq1df.z) - minimum(hq1df.z),
	widphi1   = maximum(phi1_values) - minimum(phi1_values),
	corrzphi1 = cor(hcat(hq1df.z, phi1_values), int1_weights)[1,2],
	phistd2   = std(phi2_values, int2_weights, corrected=true),
	zstd2     = std(hq2df.z   , int2_weights, corrected=true),
	widz2     = maximum(hq2df.z) - minimum(hq2df.z),
	widphi2   = maximum(phi2_values) - minimum(phi2_values),
	corrzphi2 = cor(hcat(hq2df.z, phi2_values), int2_weights)[1,2],
	# XYZ given true radius and barycentre
	x1 = xyz1[1],
	y1 = xyz1[2],
	z1 = xyz1[3],
	x2 = xyz2[1],
	y2 = xyz2[2],
	z2 = xyz2[3],
	# Average arrival time of first ntof SiPMs
	ta1 = mean(t1s[1:ntof1]),
	ta2 = mean(t2s[1:ntof2]),
	# Positions of SiPM with tmin
	xb1 = ht1.x[1],
	yb1 = ht1.y[1],
	zb1 = ht1.z[1],
	xb2 = ht2.x[1],
	yb2 = ht2.y[1],
	zb2 = ht2.z[1])

	return hq1df, hq2df, reduced_event
end

"""
	nema_analysis!(     data_vec    ::Vector{ATools.EventParameters},
						event       ::Integer,
						dc          ::DetConf,
						df1         ::DataFrame,
						df2         ::DataFrame,
						primaries   ::SubDataFrame,
						sensor_xyz  ::DataFrame,
						waveform    ::SubDataFrame,
						lor_algo    ::Function)

Fill the DataFrame for nema analysis
"""
function nema_dict!(data_vec      ::Vector{ATools.EventParameters}   ,
		    event     ::Integer     ,
		    dc        ::DetConf     ,
		    df1       ::DataFrame   ,
		    df2       ::DataFrame   ,
		    primaries ::SubDataFrame,
		    sensor_xyz::DataFrame   ,
		    waveform  ::SubDataFrame,
		    lor_algo  ::Function    )

	result = recovent(event, dc, df1, df2, primaries, sensor_xyz, waveform, lor_algo)

	if result !== nothing
		push!(data_vec, result[3])
	end
end

"""
	function nemareco(files    ::Vector{String},
					  dconf    ::DetConf,
		              file_i   ::Integer=1,
					  file_l   ::Integer=1,
					  phot     ::Bool=true
					  lor_algo ::Function=lor_maxq)

Return the evtdf DataFrame

   nsipm1  => Number of sipms in hemisphere 1 (nsipm2 for hemisphere 2)
   q1      => Total charge in hemisphere
   E1      => Energy of the gamma at interaction.
   r1      => True radius, gamma interaction point,
   phistd  =>  Phi standard deviations (weigthed by charge)
   zstd    => Z standard deviations (weigthed by charge)

   xs, ys, zs    => Position of the pair of gammas (in target),
   ux, uy, uz    => Direction vectors of gammas

   xt1, yt1, zt1 => True position of gammas,
   x1,   y1, z1  => Reco position of gammas (from barycenter and r1),
   xr1, yr1, zr  => Baricenter
   xb1, yb1, zb1 => Position of the sipm giving time stamp

   t1            => Time stamp of first photon,
   tr1           => Time stamp of first photon, after pdf and smearing
   ta1           => Average of time stamps, first five photons

"""
function nemareco(files    ::Vector{String},
				  dconf    ::DetConf,
	              file_i   ::Integer=1,
				  file_l   ::Integer=1,
				  lor_algo ::Function=lor_maxq)

	output_vector = ATools.EventParameters[]
	total_events  = 0

	for file in files[file_i:file_l]               # loop on files
		println("reading file = ", file)
		pdf = read_abc(file)            # read file

		# count the total number of events generated.
		# This, like many other sections assumes primary generation of gammas.
		total_events += nrow(pdf.primaries)

		dfs = primary_in_lxe(pdf.vertices)       # primary photons in LXe

		## We are interested in events with two primary photons in LXe
		grp_dfs = filter(x -> any(x.track_id .== 1) && any(x.track_id .== 2),
							groupby(dfs, :event_id))
		primaries = groupby(pdf.primaries, :event_id)
		waveforms = groupby(pdf.waveform , :event_id)

		for (event, vdf) in pairs(grp_dfs)
			df1 = vdf[vdf.track_id .== 1, :]
			df2 = vdf[vdf.track_id .== 2, :]

			try
				wvf = waveforms[values(event)]
				nema_dict!(output_vector, event.event_id, dconf, df1, df2,
					   primaries[values(event)], pdf.sensor_xyz,
					   wvf, lor_algo)
			catch err
				if err isa KeyError
					continue
				else
					rethrow(err)
				end
			end
    	end
	end
	return total_events, output_vector
end

max_stars_repo_name: andLaing/NReco

max_issues_repo_name: andLaing/NReco

max_forks_repo_name: andLaing/NReco


text:
include("computation/derivatives.jl")
using Plots

f = x-> x[1]*x[2] + x[2]^3 #derivative
df = x -> [x[2];  x[1]+3*x[2]^2] #Jacobian, pass in as column vector
d2f = x ->  [0 1; 1 6*x[2]] #Hessian

(h, A) = checkDerivative(f, df, [1.0;3.0], d2f, 3, 6)
polyfit(h, A, 1)
plt1 = Plots.plot(h, A, title = "Finite Difference Derivative Checker", xlabel = "log of h", ylabel = "log of error",fontfamily=font(48, "Courier") , reuse = false)
#plot(polyfit(h, A, 1), reuse = true)
println("derivative of p(z0|z)")
println(polyfit(h, A, 1))
display(plt1)

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
@doc """
      Fourth-Order Runge-Kutta Method
      Supported by:
      -PhD Kelvyn B. Sánchez kelvyn.baruc@gmail.com
 """ ->

function rk4(f, tspan, y0; h = 1)
tint = floor((tspan[end] - tspan[1])/h) + 1
      t = range(tspan[1],tspan[end],length=Int(tint))

#t = linspace(tspan[1], tspan[end],tint)
y = zeros(Int(tint), length(y0))
y[1,:] = y0

      try
           error = f(t[1],y0)

            for i=1:1:(length(t) - 1)
                  ymid = zeros(1,Int(length(y0)))
                  k1 = zeros(1,Int(length(y0)))
                  k2 = zeros(1,Int(length(y0)))
                  k3 = zeros(1,Int(length(y0)))
                  k4 = zeros(1,Int(length(y0)))

                  for j=1:length(y0)
                        k1[1,j] = f(t[i],y[i,:])[j]
                        ymid[1,j] = y[i,j] + k1[1,j]*(h/2)
                  end

                  for j=1:length(y0)
                        k2[1,j] = f(t[i] + (1/2)*h,ymid)[j]
                        ymid[1,j] = y[i,j] + k2[1,j]*(h/2)
                  end

                  for j=1:length(y0)
                        k3[1,j] = f(t[i] + (1/2)*h,ymid)[j]
                        ymid[1,j] = y[i,j] + k3[1,j]*h
                  end

                  for j=1:length(y0)
                        k4[1,j] = f(t[i] + h,ymid)[j]
                  end

                  for j=1:length(y0)
                        y[i+1,j] = y[i,j] + (1/6)*(k1[1,j] + 2*(k2[1,j] + k3[1,j]) + k4[1,j])*h
                  end
            end

            return t, y
      catch
           throw(ErrorException("The number of initial values does not match the number of equations"))
      end

 end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:5] <= 1)

                     @variable(m, -25.780717722472463 <= q <= 30.482165922896684)

                     add_NL_constraint(m, :(sigmoid(0.6681236842141707 + 0.4190368949234773*sigmoid(-0.3928529255603399 + 0.13822049866887154*sigmoid(-0.9290908758950072 + -0.9885689195552865*$(x[1]) + 0.968971534656105*$(x[2]) + 0.0899903458574629*$(x[3]) + 0.3621960275577596*$(x[4]) + 0.6428275846795262*$(x[5])) + -0.4562953764921889*sigmoid(-0.6068769875537994 + 0.23936570066126128*$(x[1]) + -0.8813597146688212*$(x[2]) + 0.7513604373502818*$(x[3]) + -0.8801053294405703*$(x[4]) + -0.09321218509736573*$(x[5])) + -0.9204155811731316*sigmoid(-0.4931215033246312 + -0.7619598530221103*$(x[1]) + -0.2076077765621065*$(x[2]) + -0.31174640065939396*$(x[3]) + -0.2011575506947998*$(x[4]) + -0.9004282460591928*$(x[5]))) + 0.7912756220587216*sigmoid(-0.5063192893501953 + 0.4673927080857103*sigmoid(-0.9290908758950072 + -0.9885689195552865*$(x[1]) + 0.968971534656105*$(x[2]) + 0.0899903458574629*$(x[3]) + 0.3621960275577596*$(x[4]) + 0.6428275846795262*$(x[5])) + -0.2643988862479607*sigmoid(-0.6068769875537994 + 0.23936570066126128*$(x[1]) + -0.8813597146688212*$(x[2]) + 0.7513604373502818*$(x[3]) + -0.8801053294405703*$(x[4]) + -0.09321218509736573*$(x[5])) + -0.8391655243989531*sigmoid(-0.4931215033246312 + -0.7619598530221103*$(x[1]) + -0.2076077765621065*$(x[2]) + -0.31174640065939396*$(x[3]) + -0.2011575506947998*$(x[4]) + -0.9004282460591928*$(x[5]))) + 0.030565575094480213*sigmoid(-0.34808094533003686 + 0.7360164287507445*sigmoid(-0.9290908758950072 + -0.9885689195552865*$(x[1]) + 0.968971534656105*$(x[2]) + 0.0899903458574629*$(x[3]) + 0.3621960275577596*$(x[4]) + 0.6428275846795262*$(x[5])) + -0.8890315695148892*sigmoid(-0.6068769875537994 + 0.23936570066126128*$(x[1]) + -0.8813597146688212*$(x[2]) + 0.7513604373502818*$(x[3]) + -0.8801053294405703*$(x[4]) + -0.09321218509736573*$(x[5])) + -0.7245715776404307*sigmoid(-0.4931215033246312 + -0.7619598530221103*$(x[1]) + -0.2076077765621065*$(x[2]) + -0.31174640065939396*$(x[3]) + -0.2011575506947998*$(x[4]) + -0.9004282460591928*$(x[5])))) + sigmoid(-0.037313120688143186 + -0.8220585173358868*sigmoid(-0.3928529255603399 + 0.13822049866887154*sigmoid(-0.9290908758950072 + -0.9885689195552865*$(x[1]) + 0.968971534656105*$(x[2]) + 0.0899903458574629*$(x[3]) + 0.3621960275577596*$(x[4]) + 0.6428275846795262*$(x[5])) + -0.4562953764921889*sigmoid(-0.6068769875537994 + 0.23936570066126128*$(x[1]) + -0.8813597146688212*$(x[2]) + 0.7513604373502818*$(x[3]) + -0.8801053294405703*$(x[4]) + -0.09321218509736573*$(x[5])) + -0.9204155811731316*sigmoid(-0.4931215033246312 + -0.7619598530221103*$(x[1]) + -0.2076077765621065*$(x[2]) + -0.31174640065939396*$(x[3]) + -0.2011575506947998*$(x[4]) + -0.9004282460591928*$(x[5]))) + 0.9374338955858046*sigmoid(-0.5063192893501953 + 0.4673927080857103*sigmoid(-0.9290908758950072 + -0.9885689195552865*$(x[1]) + 0.968971534656105*$(x[2]) + 0.0899903458574629*$(x[3]) + 0.3621960275577596*$(x[4]) + 0.6428275846795262*$(x[5])) + -0.2643988862479607*sigmoid(-0.6068769875537994 + 0.23936570066126128*$(x[1]) + -0.8813597146688212*$(x[2]) + 0.7513604373502818*$(x[3]) + -0.8801053294405703*$(x[4]) + -0.09321218509736573*$(x[5])) + -0.8391655243989531*sigmoid(-0.4931215033246312 + -0.7619598530221103*$(x[1]) + -0.2076077765621065*$(x[2]) + -0.31174640065939396*$(x[3]) + -0.2011575506947998*$(x[4]) + -0.9004282460591928*$(x[5]))) + -0.1479923276913464*sigmoid(-0.34808094533003686 + 0.7360164287507445*sigmoid(-0.9290908758950072 + -0.9885689195552865*$(x[1]) + 0.968971534656105*$(x[2]) + 0.0899903458574629*$(x[3]) + 0.3621960275577596*$(x[4]) + 0.6428275846795262*$(x[5])) + -0.8890315695148892*sigmoid(-0.6068769875537994 + 0.23936570066126128*$(x[1]) + -0.8813597146688212*$(x[2]) + 0.7513604373502818*$(x[3]) + -0.8801053294405703*$(x[4]) + -0.09321218509736573*$(x[5])) + -0.7245715776404307*sigmoid(-0.4931215033246312 + -0.7619598530221103*$(x[1]) + -0.2076077765621065*$(x[2]) + -0.31174640065939396*$(x[3]) + -0.2011575506947998*$(x[4]) + -0.9004282460591928*$(x[5])))) + sigmoid(0.5136304015740518 + -0.41595218402127365*sigmoid(-0.3928529255603399 + 0.13822049866887154*sigmoid(-0.9290908758950072 + -0.9885689195552865*$(x[1]) + 0.968971534656105*$(x[2]) + 0.0899903458574629*$(x[3]) + 0.3621960275577596*$(x[4]) + 0.6428275846795262*$(x[5])) + -0.4562953764921889*sigmoid(-0.6068769875537994 + 0.23936570066126128*$(x[1]) + -0.8813597146688212*$(x[2]) + 0.7513604373502818*$(x[3]) + -0.8801053294405703*$(x[4]) + -0.09321218509736573*$(x[5])) + -0.9204155811731316*sigmoid(-0.4931215033246312 + -0.7619598530221103*$(x[1]) + -0.2076077765621065*$(x[2]) + -0.31174640065939396*$(x[3]) + -0.2011575506947998*$(x[4]) + -0.9004282460591928*$(x[5]))) + -0.01955723268679055*sigmoid(-0.5063192893501953 + 0.4673927080857103*sigmoid(-0.9290908758950072 + -0.9885689195552865*$(x[1]) + 0.968971534656105*$(x[2]) + 0.0899903458574629*$(x[3]) + 0.3621960275577596*$(x[4]) + 0.6428275846795262*$(x[5])) + -0.2643988862479607*sigmoid(-0.6068769875537994 + 0.23936570066126128*$(x[1]) + -0.8813597146688212*$(x[2]) + 0.7513604373502818*$(x[3]) + -0.8801053294405703*$(x[4]) + -0.09321218509736573*$(x[5])) + -0.8391655243989531*sigmoid(-0.4931215033246312 + -0.7619598530221103*$(x[1]) + -0.2076077765621065*$(x[2]) + -0.31174640065939396*$(x[3]) + -0.2011575506947998*$(x[4]) + -0.9004282460591928*$(x[5]))) + 0.13292433203840748*sigmoid(-0.34808094533003686 + 0.7360164287507445*sigmoid(-0.9290908758950072 + -0.9885689195552865*$(x[1]) + 0.968971534656105*$(x[2]) + 0.0899903458574629*$(x[3]) + 0.3621960275577596*$(x[4]) + 0.6428275846795262*$(x[5])) + -0.8890315695148892*sigmoid(-0.6068769875537994 + 0.23936570066126128*$(x[1]) + -0.8813597146688212*$(x[2]) + 0.7513604373502818*$(x[3]) + -0.8801053294405703*$(x[4]) + -0.09321218509736573*$(x[5])) + -0.7245715776404307*sigmoid(-0.4931215033246312 + -0.7619598530221103*$(x[1]) + -0.2076077765621065*$(x[2]) + -0.31174640065939396*$(x[3]) + -0.2011575506947998*$(x[4]) + -0.9004282460591928*$(x[5])))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
import HarmonicBalance.load

# load the previously-saved result
current_path = @__DIR__
@test load(current_path * "/parametron_result.jld2") isa HarmonicBalance.Result


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

using JuliaFEM, Test

# Example, 2d nonlinear elasticity with surface load (mesh is read from file)

meshfile = joinpath("test_elasticity_2d_nonlinear_with_surface_load", "BLOCK_1elem.med")
mesh = aster_read_mesh(meshfile)

# define field problem
block = Problem(Elasticity, "BLOCK", 2)
block.properties.formulation = :plane_stress
block.properties.finite_strain = true
block.properties.geometric_stiffness = true

# Add volume elements
block_elements = create_elements(mesh, "BLOCK")
update!(block_elements, "youngs modulus", 288.0)
update!(block_elements, "poissons ratio", 1/3)
update!(block_elements, "displacement load 2", 576.0)
add_elements!(block, block_elements)

# Add surface elements
traction_elements = create_elements(mesh, "TOP")
update!(traction_elements, "displacement traction force 2", 288.0)
add_elements!(block, traction_elements)

# Create boundary problem, add boundary conditions:
bc_sym = Problem(Dirichlet, "symmetry bc", 2, "displacement")
bc_elements_left = create_elements(mesh, "LEFT")
bc_elements_bottom = create_elements(mesh, "BOTTOM")
update!(bc_elements_left, "displacement 1", 0.0)
update!(bc_elements_bottom, "displacement 2", 0.0)
add_elements!(bc_sym, bc_elements_left)
add_elements!(bc_sym, bc_elements_bottom)

# Create analysis, add problems to analysis and run analysis:
analysis = Analysis(Nonlinear)
add_problems!(analysis, block, bc_sym)
run!(analysis)

# Results from Code Aster:
u3_expected = [-4.92316106779943E-01, 7.96321884292103E-01]
eps_zz = -3.71128811855451E-01
eps_expected = [-3.71128532282463E-01, 1.11338615599337E+00, 0.0]
sig_expected = [ 3.36174888827909E-05, 2.23478729403118E+03, 0.0]

u3 = block("displacement", 0.0)[3]
@debug("displacement", u3)
@test isapprox(u3, u3_expected, atol=1.0e-5)


text: # used as a dummy


text:
mutable struct Channel{T}
    value::T
    isset::Bool
    moretakers::Bool
    waiters::Vector{OneWayCondition}
    lock::ReentrantLock
end

function _tryput!(channel, value)
    sender = nothing
    receiver = nothing
    lock(channel.lock)
    try
        if channel.isset
            @assert !channel.moretakers
            todelete = Int[]
            for (i, w) in pairs(channel.waiters)
                if w.task === current_task()
                    if Atomic.cas!(w.state, OWC_EMPTY, OWC_WAITING) === OWC_EMPTY
                        sender = w
                        break
                    end
                elseif w.state[] === OWC_EMPTY
                    if Atomic.cas!(w.state, OWC_EMPTY, OWC_CLOSED) === OWC_EMPTY
                        push!(todelete, i)
                    end
                end
            end
            deleteat!(channel.waiters, todelete)
            if sender === nothing
                sender = OneWayCondition()
                sender.task = current_task()
                push!(channel.waiters, sender)
            end
        else
            if channel.moretakers
                todelete = Int[]
                for (i, w) in pairs(channel.waiters)
                    if Atomic.cas!(w.state, OWC_WAITING, OWC_WAITING) === OWC_EMPTY
                        sender = w
                        break
                    end
                end
            else
                channel.isset = true
                channel.value = value
            end
        end
    finally
        unlock(channel.lock)
    end
end

function _trytake!(channel)
    lock(channel.lock)
    try
        if channel.isset
            channel.isset = false
            return Some(channel.value)
        end
        for waiter in channel.waiters
        end
    finally
        unlock(channel.lock)
    end
end


text:
using CairoMakie
using Dagger

T    = zeros(62, 62)
node = Observable(T)
fig  = heatmap(node, colorrange = (0.0, 2.0))

record(fig, "output.mkv", 1:1) do t
    node[] = T
end

fetch(Dagger.@spawn single=1 1)

max_stars_repo_name: JuliaLabs/DaggerDiffusion

max_issues_repo_name: JuliaLabs/DaggerDiffusion

max_forks_repo_name: JuliaLabs/DaggerDiffusion


text:
include("main.jl")

function runProcess(config, val, out, loadTime)
    f, dist = val
    params = first(values(dist))
    r = try
        main(config, dist, true, loadTime = loadTime, resultFolder=out)
    catch
        (NaN, NaN, NaN, false, NaN, NaN, NaN, NaN, NaN, NaN, NaN, "")
    end
    return (f, params, r)
end


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

using JuliaFEM
using JuliaFEM.Preprocess
using JuliaFEM.Testing

@testset "test 2d nonlinear residual" begin
    X = Dict{Int64, Vector{Float64}}(
        1 => [0.0, 0.0],
        2 => [1.0, 0.0],
        3 => [1.0, 1.0],
        4 => [0.0, 1.0])
    u = Dict{Int64, Vector{Float64}}(
        1 => [0.1, 0.2],
        2 => [0.3, 0.4],
        3 => [0.5, 0.6],
        4 => [0.7, 0.8])
    T = Dict{Int64, Vector{Float64}}(
        3 => [0.0, 288.0],
        4 => [0.0, 288.0])
    element = Element(Quad4, [1, 2, 3, 4])
    update!(element, "geometry", X)
    update!(element, "displacement", u)
    update!(element, "youngs modulus", 288.0)
    update!(element, "poissons ratio", 1/3)
    traction = Element(Seg2, [3, 4])
    update!(traction, "geometry", X)
    update!(traction, "displacement", u)
    update!(traction, "displacement traction force", T)

    # field problem
    block = Problem(Elasticity, "BLOCK", 2)
    block.properties.formulation = :plane_stress
    block.properties.finite_strain = true
    block.properties.geometric_stiffness = true
    push!(block, element)
    #push!(block, traction)
    assemble!(block, 0.0)
    Km = full(block.assembly.K)
    Kg = full(block.assembly.Kg)
    K = Km + Kg
    f = full(block.assembly.f)

    K_expected = [
      401.76   200.88  -123.12    -5.76  -191.52  -117.36  -87.12   -77.76
      200.88   473.76    -5.76    28.08  -117.36  -205.92  -77.76  -295.92
     -123.12    -5.76   197.28    -2.16   -12.24   -25.92  -61.92    33.84
       -5.76    28.08    -2.16   298.08   -25.92  -163.44   33.84  -162.72
     -191.52  -117.36   -12.24   -25.92   240.48   120.24  -36.72    23.04
     -117.36  -205.92   -25.92  -163.44   120.24   312.48   23.04    56.88
      -87.12   -77.76   -61.92    33.84   -36.72    23.04  185.76    20.88
      -77.76  -295.92    33.84  -162.72    23.04    56.88   20.88   401.76]
#   f_expected = [142.272, 214.272, -13.824, 58.176, -75.456, 19.584, -52.992, -4.032]
    f_expected = [142.272, 214.272, -13.824, 58.176, -75.456, -124.416, -52.992, -148.032]
    @test isapprox(K, K_expected)
    @test isapprox(f, f_expected)
end


text:
using BlackBoxOptim
using BlackBoxOptim: num_func_evals # We will update progress based on num func evals so far
using ProgressMeter

# Case 1: Optimization for a given number of function evaluations
MaxFuncEvals = 5000000

# Minimum update interval 0.5 seconds
const Prog = Progress(MaxFuncEvals, 0.5, "Optimizing...")

function callback_progress_stepper(optcontroller)
    global Prog
    ProgressMeter.update!(Prog, num_func_evals(optcontroller))
end

# Let's try this on 100-dimensional rosenbrock
function rosenbrock(x)
    return( sum( 100*( x[2:end] .- x[1:end-1].^2 ).^2 .+ ( x[1:end-1] .- 1 ).^2 ) )
end

println("Starting optimization!")
res = bboptimize(rosenbrock;
        SearchRange=(-5.0,5.0), NumDimensions = 100,
        PopulationSize=10, MaxFuncEvals=MaxFuncEvals,
        CallbackFunction = callback_progress_stepper,
        CallbackInterval = 0.0,
        TraceMode = :silent # So that output from BBO doesn't interfere with progress bar...
);

# We can adapt this to other use cases (such as max time etc) by keeping track of start time
# and then updating based on delta time etc.


text:
export graph_denman_beavers

"""
        (graph,cref)=graph_denman_beavers(k;T=ComplexF64,cref_mode=0)

Creates the graph corresponding to the Denman–Beavers iteration for the matrix
square root, using k iterations. The kwarg `cref_mode` specifies if references
to all X and Y (`cref_mode=0`) should be stored or just Y (`cref_mode=-1`) or
just X (`cref_mode=1`).

**Reference**

1. E. Denman and A. Beavers. "The matrix sign function and computations in
   systems". Applied Mathematics and Computation, 2(1), 1976. DOI:
   [10.1016/0096-3003(76)90020-5](https://doi.org/10.1016/0096-3003(76)90020-5)
"""
function graph_denman_beavers(k; T = Float64, cref_mode = 0)
    graph = Compgraph(T)

    cref = []
    for j = 0:k
        Xj = Symbol("X$j")
        Yj = Symbol("Y$j")

        if (j == 0)
            # X0=A
            # Y0=I
            Xj = :A
            Yj = :I
        end

        # Yjinv=Yj^{-1}
        Yinvj = Symbol("Yinv$j")
        if (j > 0)
            add_ldiv!(graph, Yinvj, Yj, :I)
        else
            Yinvj = :I # First step Yinv0=inv(I)=I;
        end

        # Xjinv=Xj^{-1}
        add_ldiv!(graph, Symbol("Xinv$j"), Xj, :I)

        # X_{j+1}=0.5*Xj+0.5*inv(Yj)
        add_lincomb!(graph, Symbol("X" * string(j + 1)), 0.5, Xj, 0.5, Yinvj)
        if (cref_mode >= 0)
            push!(cref, (Symbol("X" * string(j + 1)), 1))
            push!(cref, (Symbol("X" * string(j + 1)), 2))
        end

        # Y_{j+1}=0.5*Yj+0.5*inv(Xj)
        add_lincomb!(
            graph,
            Symbol("Y" * string(j + 1)),
            0.5,
            Yj,
            0.5,
            Symbol("Xinv$j"),
        )
        if (cref_mode <= 0)
            push!(cref, (Symbol("Y" * string(j + 1)), 1))
            push!(cref, (Symbol("Y" * string(j + 1)), 2))
        end

        add_output!(graph, Symbol("X" * string(j + 1)))
    end
    # Not needed to get the output
    del_node!(graph, Symbol("Y" * string(k + 1)))
    setdiff!(
        cref,
        [(Symbol("Y" * string(k + 1)), 1), (Symbol("Y" * string(k + 1)), 2)],
    )

    del_node!(graph, Symbol("Xinv" * string(k)))
    #
    return (graph, cref)
end


text:
# Glue contains wrapper functions to allow Mass2 steps to call out to other functions with
# convenient syntax. I hope to move functions out of here to better locations. Eventually this file may
# not need to exist.
import TESOptimalFiltering: filter5lag, calculate_filter, autocorrelation
using LsqFit

function compute_whitenoise_filter(pulse, selection_good)
	sumpulse = zeros(Int64, length(pulse[1]))
	n=0
	for i=1:length(selection_good) # take the mean of the pulses in a way that avoids int overflows
		if selection_good[i]
			sumpulse+=pulse[i]
			n+=1
		end
	end
	meanpulse = sumpulse/n
	filter = meanpulse-mean(meanpulse) # 0 average
	normalization = (maximum(filter)-minimum(filter))./dot(filter, filter)
	filter*normalization
end
function filter1lag(pulse, filter)
	out = Array(Float64, length(pulse))
	for i=1:length(pulse)
		out[i] = dot(pulse[i], filter)
	end
	out
end

function selectfromcriteria(x...) # placeholder, kinda ugly to use and probalby a bit slow
	iseven(length(x)) || error("x must be indicator,criteria,indicator,criteria...")
	out = trues(length(x[1]))
	for i = 1:2:length(x)
		low, high = x[i+1]
		out &= low .< x[i] .< high
	end
	out
end

function calc_dc(indicator, uncorrected, selection_good)
	σ_smooth = 1.0
	drift_correct(indicator[:][selection_good], uncorrected[:][selection_good], σ_smooth)
end

function two_exponential_pulse_for_fittng(x,p)
	# check for bad values of rise_point and fall_points
	# return vector record containing points at the quiesence value
	if p[1]<0 || p[2]<0 || p[1]==p[2]
		return zeros(Float64,length(x))
	end
	two_exponential_pulses(length(x), p[1], p[2],p[3], (p[4],), (p[5],))
end
function fit_pulse_two_exponential(pulse, pretrig_nsamples, frametime)
	pguess = zeros(Float64,5)
	pguess[1] = length(pulse)/10 # rise points
	pguess[2] = 3*pguess[1] # fall points
	pguess[3] = pulse[1] # quiescent_value
	pguess[4] = pretrig_nsamples # arrival time (really should be based on ljh)
	pguess[5] = maximum(pulse)-pulse[1] # amplitude
	fit = curve_fit(two_exponential_pulse_for_fittng, collect(1:length(pulse)), pulse, pguess)
	rise_tau_s, fall_tau_s = fit.param[1]*frametime, fit.param[2]*frametime
end

"selectfromcriteria(x...)
selectfromcriteria(indicator1, criteria1, indicator2, criteria2,...) takes 1 or more indicator,criteria inputs and
returns a bitvector that is true if all indictors pass criteria[1] .< indicator .< criteria[2]"
function selectfromcriteria(x...) # placeholder, kinda ugly to use and probalby a bit slow
	iseven(length(x)) || error("x must be indicator,criteria,indicator,criteria...")
	out = trues(length(x[1]))
	for i = 1:2:length(x)
		low, high = x[i+1]
		out &= low .< x[i] .< high
	end
	out
end

"readcontinuous(ljh,maxnpoints=50_000_000)
read from an LJH file, return pulses concatonated together to form a continuour vector
reads either the whole file, or until maxnpoints"
function readcontinuous(ljh,maxnpoints=50_000_000)
	nsamp = record_nsamples(ljh)
	maxnpoints = nsamp*div(maxnpoints, nsamp)
	npoints = min(maxnpoints, nsamp*length(ljh))
	data = zeros(UInt16,npoints)
	for (i,r) in enumerate(ljh[1:div(npoints, nsamp)])
		data[1+(i-1)*nsamp:i*nsamp] = r.data
	end
	data
end

"compute_noise_autocorr(ljh_filename, nsamp)
computer noise autorcorrelation of length namp from an ljh file of name ljh_filename"
function compute_noise_autocorr(ljh_filename, nsamp)
	data = readcontinuous(LJHGroup(ljh_filename))
	autocorrelation(data, nsamp, 1000) # final argument is max_excursion
	# this needs a revamp to reject fewer noise records and also return delta_f
end

"compute_average_pulse(pulse, selection_good)
compute the average (mean) the good pulses in pulse. good determined by bitvector selection_good"
function compute_average_pulse(pulse, selection_good)
	sumpulse = zeros(Int64, length(pulse[1]))
	n=0
	for i=1:length(selection_good) # take the mean of the pulses in a way that avoids int overflows
		if selection_good[i]
			sumpulse+=pulse[i]
			n+=1
		end
	end
	meanpulse = sumpulse/n
end

"compute_filter(average_pulse, noise_autocorr, f_3db, dt)\n
compute filter given average_pulse and noise_autorr of same length, f_3db in hz, dt in seconds
return `filter`, `vdv`"
function compute_filter(average_pulse, noise_autocorr, f_3db, dt)
	filter, variance = calculate_filter(average_pulse, noise_autocorr, f_3db, dt)
	filter*=maximum(average_pulse)-minimum(average_pulse) # like mass, normalize so you get pulseheight similar to raw units
	vdv = 1/sqrt(variance)/2.355
	filter,vdv
end

function estimate_pretrig_rms_and_postpeak_deriv_criteria(fname, pretrig_nsamples)
	# using the noise file we will try to estimate values to use for actual cuts on pulses
	# pretrig_rms and post_peak deriv should be fairly well estiamated by the noise file
	ljh = LJHGroup(fname)

	traces=[r.data for r in ljh[1:end]]
	pretrig_mean, pretrig_rms, pulse_average, pulse_rms, rise_time, postpeak_deriv,
		peak_index, peak_value, min_value=compute_summary(traces,pretrig_nsamples,frametime(ljh))

	# the distribution of pretrig_rms should follow a chisq distributions if the points are independent (IID) and have gaussian noise
	# but the chisq distribution for 256 degrees of freedom looks a lot like a gaussian, so we're just going to stick with the
	# gaussian approximation for now, plus a backup measure that says at most cut 1 % of pulses
	n_std = 10
	pretrig_rms_max = max(median(pretrig_rms)+n_std*std(pretrig_rms), StatsBase.percentile(pretrig_rms,99))
	pretrig_rms_min = max(0, median(pretrig_rms)-n_std*std(pretrig_rms))

	postpeak_deriv_max = max(median(postpeak_deriv)+n_std*std(postpeak_deriv), StatsBase.percentile(postpeak_deriv,99))
	postpeak_deriv_min = median(postpeak_deriv)-n_std*std(postpeak_deriv)
	Float64[pretrig_rms_min, pretrig_rms_max], Float64[postpeak_deriv_min, postpeak_deriv_max]
end

function estimate_peak_index_criteria(peak_index)
	mad = mean(abs(peak_index-median(peak_index)))
	med = median(peak_index)
	Float64[med-10*mad, med+10*mad]
end


text:
using Test

struct Point
    x::Float64
    y::Float64
end

# valid directions
@enum Directions e se sw w nw ne
directions = ["e", "se", "sw", "w", "nw", "ne"]
dmap = Dict(
    "e" => e,
    "ne" => ne,
    "nw" => nw,
    "se" => se,
    "sw" => sw,
    "w" => w,
)
offsets = Dict(
    e => Point(1, 0),
    ne => Point(0.5, 0.5),
    se => Point(0.5, -0.5),
    w => Point(-1, 0),
    nw => Point(-0.5, 0.5),
    sw => Point(-0.5, -0.5),
)

# determine if a series of moves are equivalent

function get_neighbors(p::Point)
    return map((o) -> Point(p.x + o.x, p.y + o.y), collect(values(offsets)))
end

function move(path::String)
    x = 0
    y = 0
    ds = Nothing
    i = 1
    while i <= length(path)
        # get next direction
        if (length(path) >= i + 1) && in(path[i:i + 1], directions)
            ds = path[i:i + 1]
            i += 2
        else
            ds = string(path[i])
            i += 1
        end

        # move
        o = offsets[dmap[ds]]
        x += o.x
        y += o.y
    end

    return Point(x, y)
end

function run(fname)
    f = open(fname, "r")
    lines = readlines(f)
    close(f)

    flipped = Dict()
    for l in lines
        p = move(l)
        if ! haskey(flipped, p)
            flipped[p] = 0
        end
        flipped[p] += 1
    end

    pts = map(x -> x[1], filter(x -> isodd(x[2]), collect(flipped)))
    # @show pts
    result = length(pts)
    println("[part 1] Ran $(fname) and got result = $(result)")
    return result, pts
end

function day(pts)
    min_y  = minimum((p) -> p.y, pts)
    max_y = maximum((p) -> p.y, pts)
    min_x  = minimum((p) -> p.x, pts)
    max_x = maximum((p) -> p.x, pts)

    # @show min_x, min_y, max_x, max_y
    before = pts
    after = Set()

    y = min_y - 1
    while y <= (max_y + 1)
        x = min_x - 1
        while x <= (max_x + 1)
            cur = Point(x, y)
            ns = get_neighbors(cur)
            num_black_tiles_adj = length(filter(n -> in(n, before), ns))
            if in(cur, before)
                # Any black tile with zero or more than 2 black tiles immediately
                # adjacent to it is flipped to white.
                if 1 <= num_black_tiles_adj <= 2
                    push!(after, cur)
                end
            else
                # Any white tile with exactly 2 black tiles immediately adjacent to
                # it is flipped to black.
                if num_black_tiles_adj == 2
                    push!(after, cur)
                end
            end

            x += 0.5
        end
        y += 0.5
    end

    return after
end

function run2(fname)
    _, pts = run(fname)

    black_tiles = Set()
    for p in pts
        push!(black_tiles, p)
    end

    println("Day = 0 $(length(black_tiles))")
    for d in 1:100
        black_tiles = day(black_tiles)
        println("Day = $(d) $(length(black_tiles))")
    end

    result = length(black_tiles)
    println("[part 2] Ran $(fname) and got result = $(result)")
    return result
end

@test move("esew") == Point(+0.5, -0.5)
@test move("nwwswee") == Point(0, 0)

@test run("24ex.txt")[1] == 10

run("24.txt")

# TODO: run2
# - run() gives us initial state
# - then, do conway style simulation
run2("24ex.txt")

run2("24.txt")


text:
abstract type AbstractMatchingPursuit{T} <: Update{T} end

# TODO: lazy evaluation for OMP via approximate submodularity? early dropping?
# TODO: rename Ar for δ to be consistent with generalized matching pursuits

############################# Matching Pursuit #################################
# should MP, OMP, have a k parameter?
# SP needs one, but for MP could be handled in outer loop
# WARNING: requires A to have unit norm columns
struct MatchingPursuit{T, AT<:AbstractMatOrFac{T}, B<:AbstractVector{T}} <: AbstractMatchingPursuit{T}
    A::AT
    b::B
    # temporary storage
    r::B # residual
    Ar::B # inner products between measurement matrix and residual
    # a approximation # can be updated efficiently
end
const MP = MatchingPursuit
function MP(A::AbstractMatOrFac, b::AbstractVector)
    n, m = size(A)
    T = eltype(A)
    r, Ar = zeros(T, n), zeros(T, m)
    MP(A, b, r, Ar)
end

function update!(P::MP, x::AbstractVector = spzeros(size(P.A, 2)))
    residual!(P, x) # IDEA: could just update approximation directly for mp
    i = argmaxinner!(P)
    x[i] += dot(@view(P.A[:, i]), P.r) # add non-zero index to x
    return x
end

# calculates k-sparse approximation to Ax = b via matching pursuit
function mp(A::AbstractMatOrFac, b::AbstractVector, k::Int, x = spzeros(size(A, 2)))
    P = MP(A, b)
    for i in 1:k
        update!(P, x)
    end
    return x
end

###################### Orthogonal Matching Pursuit #############################
# could extend: preconditioning, non-negativity constraint
struct OrthogonalMatchingPursuit{T, AT<:AbstractMatOrFac{T}, B<:AbstractVector{T},
                            V<:AbstractVector{T}, FT} <: AbstractMatchingPursuit{T}
    A::AT
    b::B
    # temporary storage
    r::V # residual
    Ar::V # inner products between measurement matrix and residual
    AiQR::FT # updatable QR factorization
end
const OMP = OrthogonalMatchingPursuit
function OMP(A::AbstractMatOrFac, b::AbstractVector, k::Integer = size(A, 1))
    n, m = size(A)
    T = eltype(A)
    r, Ar = zeros(T, n), zeros(T, m)
    AiQR = UpdatableQR(T, n, k) # initializing empty qr factorization with maximum rank k
    OMP(A, b, r, Ar, AiQR)
end

function update!(P::OMP, x::AbstractVector = spzeros(size(P.A, 2)))
    nnz(x) < size(P.A, 1) || return x
    residual!(P, x)
    i = argmaxinner!(P)
    ∉(i, x.nzind) || return x
    addindex!(x, P, i)
    ldiv!!(x.nzval, P.AiQR, P.b, P.r)
    return x
end

# approximately solves Ax = b with error tolerance ε it at most k steps
function omp(A::AbstractMatOrFac, b::AbstractVector, ε::Real, k::Int = size(A, 1))
    ε ≥ 0 || throw("ε = $ε has to be non-negative")
    P = OMP(A, b, k)
    x = spzeros(size(A, 2))
    for i in 1:k
        update!(P, x)
        norm(residual!(P, x)) ≥ ε || break
    end
    return x
end
# calculates k-sparse approximation to Ax = b via orthogonal matching pursuit
function omp(A::AbstractMatOrFac, b::AbstractVector, k::Int)
    omp(A, b, eps(eltype(A)), k)
end

function omp(A::AbstractMatOrFac, b::AbstractVector;
    max_residual = eps(eltype(A)), sparsity = size(A, 2))
    omp(A, b, max_residual, sparsity)
end

################## Generalized Orthogonal Matching Pursuit ######################
# see https://arxiv.org/pdf/1111.6664.pdf
struct GeneralizedOrthogonalMatchingPursuit{T, AT<:AbstractMatOrFac{T}, B<:AbstractVector{T},
                            V<:AbstractVector{T}, FT} <: AbstractMatchingPursuit{T}
    A::AT
    b::B
    l::Int # how many atoms to add in each iteration
    # temporary storage
    r::V # residual
    Ar::V # inner products between measurement matrix and residual
    AiQR::FT # updatable QR factorization
end
const GOMP = GeneralizedOrthogonalMatchingPursuit
# k is maximum number of non-zeros we expect to find, only important
# if an n x n is too big to allocate densely
function GOMP(A::AbstractMatOrFac, b::AbstractVector, l::Int, k::Integer = size(A, 1))
    n, m = size(A)
    T = eltype(A)
    r, Ar = zeros(T, n), zeros(T, m)
    AiQR = UpdatableQR(T, n, k) # initializing empty qr factorization with maximum rank k
    GOMP(A, b, l, r, Ar, AiQR)
end

function update!(P::GOMP, x::AbstractVector = spzeros(size(P.A, 2)), l = P.l)
    nnz(x) < size(P.A, 1) || return x
    residual!(P, x)
    i = argmaxinner!(P, l) # returns indices of l largest inner products
    addindex!(x, P, i)
    ldiv!!(x.nzval, P.AiQR, P.b, P.r)
    return x
end

# approximately solves Ax = b with error tolerance ε with maximum sparsity k
function gomp(A::AbstractMatOrFac, b::AbstractVector, l::Int, ε::Real, k::Int = size(A, 1))
    ε ≥ 0 || throw("ε = $ε has to be non-negative")
    P = GOMP(A, b, l)
    x = spzeros(size(A, 2))
    for i in 1:(k ÷ l)
        update!(P, x, l)
        norm(residual!(P, x)) ≥ ε || break
    end
    rem = mod(k, l)
    if rem > 0
        update!(P, x, rem)
    end
    return x
end
# calculates k-sparse approximation to Ax = b via orthogonal matching pursuit
function gomp(A::AbstractMatOrFac, b::AbstractVector, l::Int, k::Int)
    gomp(A, b, l, eps(eltype(A)), k)
end

function gomp(A::AbstractMatOrFac, b::AbstractVector, l::Int;
              max_residual = eps(eltype(A)), sparsity = size(A, 2))
    gomp(A, b, l, max_residual, sparsity)
end

####################### Matching Pursuit Helpers ###############################
# calculates residual of
function residual!(P::AbstractMatchingPursuit, x::AbstractVector)
    residual!(P.r, P.A, x, P.b)
end
function residual!(P, x::AbstractVector)
    residual!(P.r, P.A, x, P.b)
end
function residual!(r::AbstractVector, A::AbstractMatOrFac, x::AbstractVector, b::AbstractVector)
    copyto!(r, b)
    mul!(r, A, x, -1, 1)
end

# adds non-zero at i, adds ith column of A to AiQR, and solves resulting ls-problem
function addindex!(x::SparseVector, P::AbstractMatchingPursuit, i)
    addindex!(x, P.AiQR, @view(P.A[:, i]), i)
    return x, P
end

# overwrites r with intermediate result and y with final result, leaves b intact
function ldiv!!(y::AbstractVector, AiQR::UpdatableQR, b::AbstractVector, r::AbstractVector)
    @. r = b # reuse P.r as temporary storage to execute the least-squares solve with P.AiQR
    # println(y)
    # println(ldiv!(AiQR, r))
    @. r = b # reuse P.r as temporary storage to execute the least-squares solve with P.AiQR
    y .= ldiv!(AiQR, r)
end

########################### acquisition criteria ###############################
# WARNING: need to compute residual first
# returns index of atom with largest dot product with residual
@inline function argmaxinner!(P::AbstractMatchingPursuit) # 0 alloc
    mul!(P.Ar, P.A', P.r)
    @. P.Ar = abs(P.Ar)
    argmax(P.Ar)
end

# returns indices of k atoms with largest inner products with residual
# IDEA: use partialsortperm(P.ix, P.Ar, 1:k, rev = true)
@inline function argmaxinner!(P::AbstractMatchingPursuit, k::Int)
    mul!(P.Ar, P.A', P.r)
    @. P.Ar = abs(P.Ar)
    partialsortperm(P.Ar, 1:k, rev = true)
end

function random_acquisition!(P::AbstractMatchingPursuit, x::SparseVector, k::Int)
    ind = sample(1:length(x), k, replace = false) # or could initialize with max inner products
    sort!(ind)
    @. x[ind] = NaN
    for i in ind
        add_column!(P.AiQR, @view(P.A[:, i]))
    end
    ldiv!(x.nzval, P.AiQR, P.b)
    return x
end

# using oblivious algorithm to initialize P
function oblivious_acquisition!(P::AbstractMatchingPursuit, x::SparseVector, k::Int = P.k)
    residual!(P, x)
    ind = argmaxinner!(P, k)
    sort!(ind)
    @. x[ind] = NaN
    for i in ind
        add_column!(P.AiQR, @view P.A[:, i])
    end
    ldiv!(x.nzval, P.AiQR, P.b)
end

################################################################################
factorize!(P::AbstractMatchingPursuit, x::SparseVector) = factorize!(P, x.nzind)
# IDEA: could have non-qr-based solution
@inline function factorize!(P::AbstractMatchingPursuit, nzind::Vector{<:Int})
    n = size(P.A, 1)
    k = length(nzind)
    Ai = reshape(@view(P.Ai[1:n*k]), n, k)
    @. Ai = @view(P.A[:, nzind])
    return qr!(Ai) # this still allocates a little memory
end


text:
module HetaSimulator

  using DiffEqBase: isempty
  using Reexport
  # heta compiler support
  using NodeJS
  # diffeq-related pkgs
  using LabelledArrays
  using DataStructures
  using DiffEqBase
  @reexport using OrdinaryDiffEq
  @reexport using DiffEqBase.EnsembleAnalysis
  using Sundials
  # fitting

  using NLopt
  # utils
  @reexport using DataFrames
  @reexport using Distributions
  using LinearAlgebra: pinv, diag
  using Distributed
  using RecursiveArrayTools: vecvec_to_mat, VectorOfArray
  using ProgressMeter
  # measurements
  using CSV
  using XLSX
  #plots
  using RecipesBase

  const HetaSimulatorDir = dirname(Base.@__DIR__)

  include("types.jl")
  include("heta_connect.jl")
  include("utils.jl")
  include("events.jl")
  include("measurements.jl")
  include("ode_problem.jl")
  include("scenario.jl")
  include("parameters.jl")
  include("simulate.jl")
  include("saving.jl")
  include("solution_interface.jl")
  include("plots.jl")
  include("loss.jl")
  include("fit.jl")
  include("monte_carlo.jl")
  include("import_platform.jl")
  include("gsa.jl")
  include("save_as_heta.jl")

  export heta_update, heta_update_dev, heta_build, load_platform, load_jlplatform, load_jlmodel
  export Platform, Model, Scenario, Params
  export read_scenarios, add_scenarios!
  export read_measurements, add_measurements!, measurements_as_table
  export read_parameters
  export models, scenarios, constants, records, switchers, events, parameters, events_active, events_save, observables  # variables, dynamic, static
  export measurements, tspan, saveat
  export CVODE_BDF, CVODE_Adams
  export optim, obj
  export sim, mc
  export fit, loss
  export HetaSimulatorDir
  export update
  export times, vals, status, status_summary
  export save_results, read_mcvecs
  export gsa, pearson, partial, standard
  export save_as_heta
end