text:
module SeismicVelocity

# This implements different methods to compute seismic velocities
#
# If you want to add a new method here, feel free to do so.
# Remember to also export the function name in GeoParams.jl (in addition to here)

using Parameters, LaTeXStrings, Unitful
using ..Units
using ..PhaseDiagrams
#using ..MaterialParameters
using GeoParams: AbstractMaterialParam, AbstractMaterialParamsStruct
import Base.show

abstract type AbstractSeismicVelocity <: AbstractMaterialParam end

export  ComputePwaveVelocity,  ComputeSwaveVelocity,    # calculation routines
        ComputePwaveVelocity!, ComputeSwaveVelocity!,   # in place calculation
        ConstantSeismicVelocity                         # constant

# Constant Velocity -------------------------------------------------------
"""
    ConstantSeismicVelocity(Vp=8.1 km/s, Vs=4.5km/s)

Set a constant seismic P and S-wave velocity:
```math
    V_p = cst \\
    V_s = cst
```
where ``V_p, V_s`` are the P-wave and S-wave velocities [``km/s``].
"""
@with_kw_noshow mutable struct ConstantSeismicVelocity <: AbstractSeismicVelocity
    equation::LaTeXString   =   L"v_p = cst \\ v_s = cst"
    Vp::GeoUnit             =  8.1km/s               # P-wave velocity
    Vs::GeoUnit             =  4.5km/s               # S-wave velocity
end

# Calculation routines
function ComputePwaveVelocity(P,T, s::ConstantSeismicVelocity)
    @unpack Vp   = s
    if length(T)>1
        return Value(Vp).*ones(size(T))
    else
        return NumValue(Vp)
    end
end

function ComputeSwaveVelocity(P,T, s::ConstantSeismicVelocity)
    @unpack Vs   = s
    if length(T)>1
        return Value(Vs).*ones(size(T))
    else
        return NumValue(Vs)
    end
end

function ComputePwaveVelocity!(Vp_array::AbstractArray{<:AbstractFloat,N},P::AbstractArray{<:AbstractFloat,N},T::AbstractArray{<:AbstractFloat,N}, s::ConstantSeismicVelocity) where N
    @unpack Vp   = s

    Vp_array .= NumValue(Vp)

    return nothing
end

function ComputeSwaveVelocity!(Vs_array::AbstractArray{<:AbstractFloat,N},P::AbstractArray{<:AbstractFloat,N},T::AbstractArray{<:AbstractFloat,N}, s::ConstantSeismicVelocity) where N
    @unpack Vs   = s

    Vs_array .= NumValue(Vs)

    return nothing
end

# Print info
function show(io::IO, g::ConstantSeismicVelocity)
    print(io, "Constant seismic velocity: Vp=$(Value(g.Vp)), Vs=$(Value(g.Vs))")
end
#-------------------------------------------------------------------------

#-------------------------------------------------------------------------
# Phase diagrams
"""
    ComputePwaveVelocity(P,T, s::PhaseDiagram_LookupTable)

Interpolates `Vp` as a function of `T,P`
"""
function ComputePwaveVelocity(P,T, s::PhaseDiagram_LookupTable)
    return s.Vp.(T,P)
end

"""
    ComputeSwaveVelocity(P,T, s::PhaseDiagram_LookupTable)

Interpolates `Vs` as a function of `T,P`
"""
function ComputeSwaveVelocity(P,T, s::PhaseDiagram_LookupTable)
    return s.Vs.(T,P)
end

"""
    ComputePwaveVelocity!(Vp_array::AbstractArray{<:AbstractFloat}, P::AbstractArray{<:AbstractFloat},T::AbstractArray{<:AbstractFloat}, s::PhaseDiagram_LookupTable)

In-place computation of P-wave velocity as a function of `T,P`, in case we are using a lookup table.
"""
function ComputePwaveVelocity!(Vp_array::AbstractArray{<:AbstractFloat}, P::AbstractArray{<:AbstractFloat},T::AbstractArray{<:AbstractFloat}, s::PhaseDiagram_LookupTable)
    Vp_array[:] = s.Vp.(T,P)
    return nothing
end

"""
    ComputeSwaveVelocity!(Vs_array::AbstractArray{<:AbstractFloat}, P::AbstractArray{<:AbstractFloat},T::AbstractArray{<:AbstractFloat}, s::PhaseDiagram_LookupTable)

In-place computation of S-wave velocity as a function of `T,P`, in case we are using a lookup table.
"""
function ComputeSwaveVelocity!(Vs_array::AbstractArray{<:AbstractFloat}, P::AbstractArray{<:AbstractFloat},T::AbstractArray{<:AbstractFloat}, s::PhaseDiagram_LookupTable)
    Vs_array[:] = s.Vs.(T,P)
    return nothing
end

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

"""
    ComputePwaveVelocity!(Vp_array::AbstractArray{<:AbstractFloat}, Phases::AbstractArray{<:Integer}, P::AbstractArray{<:AbstractFloat},T::AbstractArray{<:AbstractFloat}, MatParam::AbstractArray{<:AbstractMaterialParamsStruct})

In-place computation of P-wave velocity `Vp` for the whole domain and all phases, in case a vector with phase properties `MatParam` is provided, along with `P` and `T` arrays.
This assumes that the `Phase` of every point is specified as an Integer in the `Phases` array.

"""
function ComputePwaveVelocity!(Vp_array::AbstractArray{<:AbstractFloat, N}, Phases::AbstractArray{<:Integer, N}, P::AbstractArray{<:AbstractFloat, N},T::AbstractArray{<:AbstractFloat, N}, MatParam::AbstractArray{<:AbstractMaterialParamsStruct, 1}) where N

    for i = 1:length(MatParam)

        if !isnothing(MatParam[i].SeismicVelocity)
            # Create views into arrays (so we don't have to allocate)
            ind = Phases .== MatParam[i].Phase
            Vp_local    =   view(Vp_array, ind )
            P_local     =   view(P  , ind )
            T_local     =   view(T  , ind )

            ComputePwaveVelocity!(Vp_local, P_local, T_local, MatParam[i].SeismicVelocity[1] )
        end

    end

end

"""
    ComputeSwaveVelocity!(Vs_array::AbstractArray{<:AbstractFloat}, Phases::AbstractArray{<:Integer}, P::AbstractArray{<:AbstractFloat},T::AbstractArray{<:AbstractFloat}, MatParam::AbstractArray{<:AbstractMaterialParamsStruct})

In-place computation of S-wave velocity `Vp` for the whole domain and all phases, in case a vector with phase properties `MatParam` is provided, along with `P` and `T` arrays.
This assumes that the `Phase` of every point is specified as an Integer in the `Phases` array.

"""
function ComputeSwaveVelocity!(Vs_array::AbstractArray{<:AbstractFloat, N}, Phases::AbstractArray{<:Integer, N}, P::AbstractArray{<:AbstractFloat, N},T::AbstractArray{<:AbstractFloat, N}, MatParam::AbstractArray{<:AbstractMaterialParamsStruct, 1}) where N

    for i = 1:length(MatParam)

        if !isnothing(MatParam[i].SeismicVelocity)
            # Create views into arrays (so we don't have to allocate)
            ind = Phases .== MatParam[i].Phase
            Vs_local    =   view(Vs_array, ind )
            P_local     =   view(P  , ind )
            T_local     =   view(T  , ind )

            ComputeSwaveVelocity!(Vs_local, P_local, T_local, MatParam[i].SeismicVelocity[1] )
        end

    end

end

"""
    ComputePwaveVelocity!(Vp_array::AbstractArray{<:AbstractFloat,N}, PhaseRatios::AbstractArray{<:AbstractFloat, M}, P::AbstractArray{<:AbstractFloat,N},T::AbstractArray{<:AbstractFloat,N}, MatParam::AbstractArray{<:AbstractMaterialParamsStruct})

In-place computation of seismic P-wave velocity `Vp` for the whole domain and all phases, in case a vector with phase properties `MatParam` is provided, along with `P` and `T` arrays.
This assumes that the `PhaseRatio` of every point is specified as an Integer in the `PhaseRatios` array, which has one dimension more than the data arrays (and has a phase fraction between 0-1)

"""
function ComputePwaveVelocity!(Vp_array::AbstractArray{<:AbstractFloat, N}, PhaseRatios::AbstractArray{<:AbstractFloat, M}, P::AbstractArray{<:AbstractFloat, N},T::AbstractArray{<:AbstractFloat, N}, MatParam::AbstractArray{<:AbstractMaterialParamsStruct, 1}) where {N,M}

    if M!=(N+1)
        error("The PhaseRatios array should have one dimension more than the other arrays")
    end

    Vp_array .= 0.0;
    for i = 1:length(MatParam)

        Vp_local    = zeros(size(Vp_array))
        Fraction    = selectdim(PhaseRatios,M,i);
        if (maximum(Fraction)>0.0) & (!isnothing(MatParam[i].SeismicVelocity))

            ComputePwaveVelocity!(Vp_local, P, T, MatParam[i].SeismicVelocity[1] )

            Vp_array .= Vp_array .+ Vp_local.*Fraction
        end

    end

end

"""
    ComputeSwaveVelocity!(Vs_array::AbstractArray{<:AbstractFloat,N}, PhaseRatios::AbstractArray{<:AbstractFloat, M}, P::AbstractArray{<:AbstractFloat,N},T::AbstractArray{<:AbstractFloat,N}, MatParam::AbstractArray{<:AbstractMaterialParamsStruct})

In-place computation of seismic S-wave velocity `Vs` for the whole domain and all phases, in case a vector with phase properties `MatParam` is provided, along with `P` and `T` arrays.
This assumes that the `PhaseRatio` of every point is specified as an Integer in the `PhaseRatios` array, which has one dimension more than the data arrays (and has a phase fraction between 0-1)

"""
function ComputeSwaveVelocity!(Vs_array::AbstractArray{<:AbstractFloat, N}, PhaseRatios::AbstractArray{<:AbstractFloat, M}, P::AbstractArray{<:AbstractFloat, N},T::AbstractArray{<:AbstractFloat, N}, MatParam::AbstractArray{<:AbstractMaterialParamsStruct, 1}) where {N,M}

    if M!=(N+1)
        error("The PhaseRatios array should have one dimension more than the other arrays")
    end

    Vs_array .= 0.0;
    for i = 1:length(MatParam)

        Vs_local    = zeros(size(Vs_array))
        Fraction    = selectdim(PhaseRatios,M,i);
        if (maximum(Fraction)>0.0) & (!isnothing(MatParam[i].SeismicVelocity))

            ComputeSwaveVelocity!(Vs_local, P, T, MatParam[i].SeismicVelocity[1] )

            Vs_array .= Vs_array .+ Vs_local.*Fraction
        end

    end

end

end


text:
using DiffEqJump, DiffEqBase, OrdinaryDiffEq, StochasticDiffEq, Statistics
using Test
using StableRNGs
rng = StableRNG(12345)

rate = (u,p,t) -> 1.0*u[1]
affect! = function (integrator)
  integrator.u[1] = 1.0
end
jump1 = ConstantRateJump(rate,affect!)

prob = DiscreteProblem([10],(0.0,50.0))
prob_control = DiscreteProblem([10],(0.0,50.0))

jump_prob = JumpProblem(prob,Direct(),jump1; rng=rng)
jump_prob_control = JumpProblem(prob_control,Direct(),jump1; rng=rng)

coupling_map = [(1, 1)]
coupled_prob = SplitCoupledJumpProblem(jump_prob,jump_prob_control,Direct(),coupling_map; rng=rng)

@time sol =  solve(coupled_prob,FunctionMap())
@time solve(jump_prob,FunctionMap())
@test [s[1]-s[2] for s in sol.u] == zeros(length(sol.t)) # coupling two copies of the same process should give zero

rate = (u,p,t) -> 1.0
affect! = function (integrator)
  integrator.u[1] = 1.0
end
jump1 = ConstantRateJump(rate,affect!)
rate = (u,p,t) -> 2.0
jump2 = ConstantRateJump(rate,affect!)

f = function (du,u,p,t)
  du[1] = u[1]
end
g = function (du,u,p,t)
  du[1] = 0.1
end

# Jump ODE to jump ODE
prob = ODEProblem(f,[1.],(0.0,1.0))
prob_control = ODEProblem(f,[1.],(0.0,1.0))
jump_prob = JumpProblem(prob,Direct(),jump1; rng=rng)
jump_prob_control = JumpProblem(prob_control,Direct(),jump2; rng=rng)
coupled_prob = SplitCoupledJumpProblem(jump_prob,jump_prob_control,Direct(),coupling_map; rng=rng)
sol =  solve(coupled_prob,Tsit5())
@test mean([abs(s[1]-s[2]) for s in sol.u])<=5.

# Jump SDE to Jump SDE
prob = SDEProblem(f,g,[1.],(0.0,1.0))
prob_control = SDEProblem(f,g,[1.],(0.0,1.0))
jump_prob = JumpProblem(prob,Direct(),jump1; rng=rng)
jump_prob_control = JumpProblem(prob_control,Direct(),jump1; rng=rng)
coupled_prob = SplitCoupledJumpProblem(jump_prob,jump_prob_control,Direct(),coupling_map; rng=rng)
sol =  solve(coupled_prob,SRIW1())
@test mean([abs(s[1]-s[2]) for s in sol.u])<=5.

# Jump SDE to Jump ODE
prob = ODEProblem(f,[1.],(0.0,1.0))
prob_control = SDEProblem(f,g,[1.],(0.0,1.0))
jump_prob = JumpProblem(prob,Direct(),jump1; rng=rng)
jump_prob_control = JumpProblem(prob_control,Direct(),jump1; rng=rng)
coupled_prob = SplitCoupledJumpProblem(jump_prob,jump_prob_control,Direct(),coupling_map; rng=rng)
sol =  solve(coupled_prob,SRIW1())
@test mean([abs(s[1]-s[2]) for s in sol.u])<=5.

# Jump SDE to Discrete
rate = (u,p,t) -> 1.
affect! = function (integrator)
  integrator.u[1] += 1.
end
prob = DiscreteProblem([1.],(0.0,1.0))
prob_control = SDEProblem(f,g,[1.],(0.0,1.0))
jump_prob = JumpProblem(prob,Direct(),jump1; rng=rng)
jump_prob_control = JumpProblem(prob_control,Direct(),jump1; rng=rng)
coupled_prob = SplitCoupledJumpProblem(jump_prob,jump_prob_control,Direct(),coupling_map; rng=rng)
sol =  solve(coupled_prob,SRIW1())

# test mass action jumps coupled to ODE
# 0 -> A (stochasic) and A -> 0 (ODE)
rate        = [100.]
react_stoch = [Vector{Pair{Int,Int}}()]
net_stoch   = [[1 => 1]]
majumps     = MassActionJump(rate, react_stoch, net_stoch)
f = function (du,u,p,t)
  du[1] = -1.0*u[1]
end
odeprob     = ODEProblem(f,[10.0],(0.0,10.0))
jump_prob   = JumpProblem(odeprob, Direct(), majumps, save_positions=(false,false); rng=rng)
Nsims = 8000
Amean = 0.
for i in 1:Nsims
  global Amean
  local sol = solve(jump_prob,Tsit5(),saveat=10.)
  Amean += sol[1,end]
end
Amean /= Nsims
actmean = 100. + (10.0-100.0)*exp(-1.0*10.0)
#println(abs(Amean-actmean)/actmean)
@test abs(actmean - Amean) < .02 * actmean

ext: jl


text:
module SparseMatrixAssemblers

using Test
using Gridap.Arrays
using Gridap.TensorValues
using Gridap.ReferenceFEs
using Gridap.Geometry
using Gridap.Fields
using Gridap.Algebra
using SparseArrays
using SparseMatricesCSR
using Gridap.FESpaces
using Gridap.CellData
using Gridap.Algebra

domain =(0,1,0,1)
partition = (2,2)
model = CartesianDiscreteModel(domain,partition)

reffe = ReferenceFE(lagrangian,Float64,1)
V = FESpace(model,reffe,dirichlet_tags=[1,2,3,4,6,5])
U = V

b(x) = x[2]

v = get_fe_basis(V)
u = get_trial_fe_basis(U)

degree = 2
trian = get_triangulation(model)
quad = CellQuadrature(trian,degree)

btrian = BoundaryTriangulation(model)
bquad = CellQuadrature(btrian,degree)

b0trian = Triangulation(model,Int[])
b0quad = CellQuadrature(b0trian,degree)

cellmat = integrate(∇(v)⊙∇(u),quad)
cellvec = integrate(v⊙b,quad)
cellmatvec = pair_arrays(cellmat,cellvec)
rows = get_cell_dof_ids(V,trian)
cols = get_cell_dof_ids(U,trian)
cellmat_c = attach_constraints_cols(U,cellmat,trian)
cellmat_rc = attach_constraints_rows(V,cellmat_c,trian)
cellvec_r = attach_constraints_rows(V,cellvec,trian)
cellmatvec_c = attach_constraints_cols(U,cellmatvec,trian)
cellmatvec_rc = attach_constraints_rows(V,cellmatvec_c,trian)

bcellmat = integrate(v*u,bquad)
bcellvec = integrate(v*3,bquad)
bcellmatvec = pair_arrays(bcellmat,bcellvec)
brows = get_cell_dof_ids(V,btrian)
bcols = get_cell_dof_ids(U,btrian)
bcellmat_c = attach_constraints_cols(U,bcellmat,btrian)
bcellmat_rc = attach_constraints_rows(V,bcellmat_c,btrian)
bcellvec_r = attach_constraints_rows(V,bcellvec,btrian)
bcellmatvec_c = attach_constraints_cols(U,bcellmatvec,btrian)
bcellmatvec_rc = attach_constraints_rows(V,bcellmatvec_c,btrian)

b0cellmat = integrate(v*u,b0quad)
b0cellvec = integrate(v*3,b0quad)
b0cellmatvec = pair_arrays(b0cellmat,b0cellvec)
b0rows = get_cell_dof_ids(V,b0trian)
b0cols = get_cell_dof_ids(U,b0trian)
@test length(b0rows) == 0
b0cellmat_c = attach_constraints_cols(U,b0cellmat,b0trian)
b0cellmat_rc = attach_constraints_rows(V,b0cellmat_c,b0trian)
b0cellvec_r = attach_constraints_rows(V,b0cellvec,b0trian)
b0cellmatvec_c = attach_constraints_cols(U,b0cellmatvec,b0trian)
b0cellmatvec_rc = attach_constraints_rows(V,b0cellmatvec_c,b0trian)

term_to_cellmat = [cellmat_rc, bcellmat_rc, b0cellmat_rc]
term_to_cellvec = [cellvec, bcellvec, b0cellvec]
term_to_rows = [rows, brows, b0rows]
term_to_cols = [cols, bcols, b0cols]
term_to_cellmatvec = [ cellmatvec, bcellmatvec, b0cellmatvec ]

mtypes = [
  SparseMatrixCSC{Float64,Int},
  SparseMatrixCSR{0,Float64,Int},
  SparseMatrixCSR{1,Float64,Int},
  SymSparseMatrixCSR{0,Float64,Int},
  SymSparseMatrixCSR{1,Float64,Int}]

for T in mtypes

  matvecdata = ( term_to_cellmatvec , term_to_rows, term_to_cols)
  matdata = (term_to_cellmat,term_to_rows,term_to_cols)
  vecdata = (term_to_cellvec,term_to_rows)
  data = (matvecdata,matdata,vecdata)

  assem = SparseMatrixAssembler(T,Vector{Float64},U,V)
  test_sparse_matrix_assembler(assem,matdata,vecdata,data)

  strategy = GenericAssemblyStrategy(row->row,col->col,row->true,col->true)

  assem2 = SparseMatrixAssembler(T,Vector{Float64},U,V,strategy)
  test_sparse_matrix_assembler(assem2,matdata,vecdata,data)

  matdata = ([cellmat],[rows],[cols])
  vecdata = ([cellvec],[rows])

  mat = assemble_matrix(assem,matdata)
  vec = assemble_vector(assem,vecdata)

  x = mat \ vec

  assemble_matrix!(mat,assem,matdata)
  assemble_vector!(vec,assem,vecdata)

  x2 = mat \ vec
  @test x ≈ x2

  @test vec ≈ [0.0625, 0.125, 0.0625]
  @test mat[1, 1]  ≈  1.333333333333333
  @test mat[2, 1]  ≈ -0.33333333333333
  @test mat[1, 2]  ≈ -0.33333333333333
  @test mat[2, 2]  ≈ 2.666666666666666
  @test mat[3, 2]  ≈ -0.33333333333333
  @test mat[2, 3]  ≈ -0.33333333333333
  @test mat[3, 3]  ≈ 1.333333333333333

  data = (([cellmatvec],[rows],[cols]),([],[],[]),([],[]))
  mat, vec = allocate_matrix_and_vector(assem,data)
  assemble_matrix_and_vector!(mat,vec,assem,data)
  assemble_matrix_and_vector!(mat,vec,assem,data)

  @test vec ≈ [0.0625, 0.125, 0.0625]
  @test mat[1, 1]  ≈  1.333333333333333
  @test mat[2, 1]  ≈ -0.33333333333333

  x3 = mat \ vec
  @test x ≈ x3

  mat, vec = assemble_matrix_and_vector(assem,data)

  x4 = mat \ vec
  @test x ≈ x4

  @test vec ≈ [0.0625, 0.125, 0.0625]
  @test mat[1, 1]  ≈  1.333333333333333
  @test mat[2, 1]  ≈ -0.33333333333333

  mat, vec = assemble_matrix_and_vector(assem2,data)

  x4 = mat \ vec
  @test x ≈ x4

  @test vec ≈ [0.0625, 0.125, 0.0625]
  @test mat[1, 1]  ≈  1.333333333333333
  @test mat[2, 1]  ≈ -0.33333333333333

end

strian = SkeletonTriangulation(model)
squad = CellQuadrature(strian,degree)

scellmat = integrate(jump(v)*u.⁻,squad)
@test isa(scellmat[1],ArrayBlock)
@test scellmat[1][1,1] === nothing
@test scellmat[1][2,1] === nothing
@test scellmat[1][1,2] !== nothing
@test scellmat[1][2,2] !== nothing

scellvec = integrate(mean(v*3),squad)
@test isa(scellvec[1],ArrayBlock)
scellmatvec = pair_arrays(scellmat,scellvec)
srows = get_cell_dof_ids(V,strian)
scols = get_cell_dof_ids(U,strian)
zh = zero(V)
scellvals = get_cell_dof_values(zh,strian)
scellmatvec = attach_dirichlet(scellmatvec,scellvals)

matvecdata = ([scellmatvec],[srows],[scols])
matdata = ([scellmat],[srows],[scols])
vecdata = ([scellvec],[srows])
data = (matvecdata,matdata,vecdata)

assem = SparseMatrixAssembler(U,V)
test_sparse_matrix_assembler(assem,matdata,vecdata,data)

A = assemble_matrix(assem,matdata)
@test A == zeros(num_free_dofs(V),num_free_dofs(U))

assem = SparseMatrixAssembler(
  SparseMatrixBuilder(SparseMatrixCSC{Float64,Int},MinCPU()),U,V)
test_sparse_matrix_assembler(assem,matdata,vecdata,data)

assem = SparseMatrixAssembler(
  SparseMatrixBuilder(SparseMatrixCSC{Float64,Int},MinMemory()),U,V)
test_sparse_matrix_assembler(assem,matdata,vecdata,data)

end # module


text:
# hello.jl

println("Hello, World!");

max_stars_repo_name: neilpanchal/Julia-Algorithms

max_issues_repo_name: neilpanchal/Julia-Algorithms

max_forks_repo_name: neilpanchal/Julia-Algorithms


text:
@inline norm{T}(tv::TransposedVector{T}) = norm(transpose(tv))

"""
    norm(transposed_vector, [q = 2])

Takes the q-norm of a transposed vector, which is equivalent to the p-norm with
value `p = q/(1-q)`. They coincide at `p = q = 2`.

The difference in norm between a vector space and its dual arrises to preserve
the relationship between duality and the inner product, and the result is
consistent with the p-norm of `1 × n` matrix.
"""
@inline norm{T}(tv::TransposedVector{T}, q::Real) = q == Inf ? norm(transpose(tv), 1) : norm(transpose(tv), q/(1-q))


text:
using SafeTestsets

# @safetestset "BanditISTests" begin include("banditistests.jl") end


text:
@testset "Slater-Koster-Core" begin

@info("Slater Koster Core Tests...")
using ACEtb.SlaterKoster, Test, LinearAlgebra
import ACEtb.SlaterKoster.CodeGeneration
using ACEtb.SlaterKoster: SKH, sk2cart, cart2sk, sk2cart_num, cart2sk_num, allbonds, nbonds, index, max_locidx

SK = ACEtb.SlaterKoster
CG = ACEtb.SlaterKoster.CodeGeneration

# temporarily keep this to test against python code
using PyCall
_pysys = pyimport("sys")
push!(_pysys."path", @__DIR__()[1:end-4]*"src")
_codegen = pyimport("codegen")

alloc_block(H::SKH) = zeros(max_locidx(H::SKH), max_locidx(H::SKH))

@info("carttospher test")
for n = 1:5
   U = rand(3) .- 0.5
   U = U / norm(U)
   φ, θ = SK.carttospher(U...)
   r, α, β = _codegen.carttospher(U...)
   print((@test φ ≈ α), " ")
   print((@test θ ≈ β), " ")
   V = [cos(φ)*sin(θ), sin(φ)*sin(θ), cos(θ)]
   print((@test U ≈ V), " ")
end
println()

#@info("New implementation: sp")
#orbitals = [sko"s", sko"p"]
#bonds = [skb"ssσ", skb"spσ", skb"ppσ", skb"ppπ"]
#println(@test bonds == allbonds(orbitals))
## H = SKH(orbitals, bonds)   # TODO: need to fix this test somehow!!! (or drop it)
## println(@test H == SKH("sp"))
#H = skh"sp"
#for n = 1:5
#   V = rand(length(bonds))
#   U = rand(3) .- 0.5
#   U /= norm(U)
#   Hnew = sk2cart(H, U, V)
#   Hold = SK.OldSK.sk4!(U, V, zeros(4,4))
#   perm = [1,3,4,2]
#   println(@test Hnew ≈ Hold[perm, perm])
#end

#@info("New implementation: spd")
#orbitals = [sko"s", sko"p", sko"d"]
#H = SKH(orbitals) # generate bonds automatically
#println(@test H == SKH("spd"))
#H = SKH("spd")
#for n = 1:5
#   V = rand(nbonds(H))
#   U = rand(3) .- 0.5
#   U /= norm(U)
#   Hnew = sk2cart(H, U, V)
#   Hold = SK.OldSK.sk9!(U, V, zeros(9,9))
#   perm = [1,3,4,2,5,6,9,7,8]
#   println(@test Hnew ≈ Hold[perm, perm])
#end

#@info("New implementation: sp")
#orbitals = [sko"s", sko"p"]
#bonds = index([skb"ssσ", skb"spσ", skb"ppσ", skb"ppπ"])
#println(@test bonds == allbonds(orbitals))
## H = SKH(orbitals, bonds)   # TODO: need to fix this test somehow!!! (or drop it)
## println(@test H == SKH("sp"))
#H = skh"sp"
#for n = 1:5
#   V = rand(length(bonds))
#   U = rand(3) .- 0.5
#   U /= norm(U)
#   Hnew = sk2cart(H, U, V)
#   Hold = SK.OldSK.sk4!(U, V, zeros(4,4))
#   perm = [1,3,4,2]
#   println(@test Hnew ≈ Hold[perm, perm])
#end

@info("sk2cart to cart2sk : spd")
orbitals = [sko"s", sko"p", sko"d"]
H = SKH(orbitals) # generate bonds automatically
println(@test H == SKH("spd"))
H = SKH("spd")
#for n = 1:5
V = rand(nbonds(H))
U = rand(3) .- 0.5
U /= norm(U)
Hnew = sk2cart(H, U, V)
Vnew = cart2sk(H, U, Hnew)
println(@test V ≈ Vnew)
#end

@info("sk2cart to cart2sk : spspd")
orbitals = [sko"s", sko"p", sko"s", sko"p", sko"d"]
H = SKH(orbitals) # generate bonds automatically
println(@test H == SKH("spspd"))
H = SKH("spspd")
for n = 1:5
   V = rand(nbonds(H))
   U = rand(3) .- 0.5
   U /= norm(U)
   Hnew = sk2cart(H, U, V)
   Vnew = cart2sk(H, U, Hnew)
   println(@test V ≈ Vnew)
end

@info("sk2cart to cart2sk : sspp")
orbitals = [sko"s", sko"s", sko"p", sko"p"]
H = SKH(orbitals) # generate bonds automatically
println(@test H == SKH("sspp"))
H = SKH("spd")
V = rand(nbonds(H))
U = rand(3) .- 0.5
U /= norm(U)
Hnew = sk2cart(H, U, V)
Vnew = cart2sk(H, U, Hnew)
println(@test V ≈ Vnew)

#@info("sk2cart to cart2sk : sspp")
#orbitals = [sko"s", sko"s", sko"p", sko"p"]
#H = SKH(orbitals) # generate bonds automatically
#println(@test bonds == allbonds(orbitals))
#H = SKH("sspp")
#for n = 1:5
#   V = rand(nbonds(H))
#   U = rand(3) .- 0.5
#   U /= norm(U)
#   Hnew = sk2cart(H, U, V)
#   Vnew = cart2sk(H, U, Hnew)
#   Hnew2 = sk2cart(H, U, Vnew)
#   println(@test V ≈ Vnew)
#   println(@test Hnew ≈ Hnew2)
#end

#alloc_block(H::SKH) = zeros(max_locidx(H::SKH), max_locidx(H::SKH))

#@info("sk2cart to cart2sk : spspd")
#orbitals = [sko"s", sko"p", sko"s", sko"p", sko"d"]
#H = SKH(orbitals) # generate bonds automatically
#println(@test H == SKH("spspd"))
#H = SKH("spspd")
#for n = 1:5
#   HE = alloc_block(H)
#   for i=1:max_locidx(H)
#      for j=i:max_locidx(H)
#         HE[i,j] = rand(1)[1]
#         if i != j
#            HE[j,i] = HE[i,j]
#         end
#      end
#   end
#   U = rand(3) .- 0.5
#   U /= norm(U)
#   V = cart2sk(H, U, HE)
#   Hnew = sk2cart(H, U, V)
#   Vnew = cart2sk(H, U, Hnew)
#   println(norm(HE - Hnew, Inf))
#   println(@test V ≈ Vnew)
#   println(@test HE ≈ Hnew)
#end

@info("sk2cart vs sk2cart_num : spspd")
orbitals = [sko"s", sko"p", sko"s", sko"p", sko"d"]
H = SKH(orbitals) # generate bonds automatically
println(@test H == SKH("spspd"))
H = SKH("spspd")
#for n = 1:5
V = rand(nbonds(H))
U = rand(3) .- 0.5
U /= norm(U)
Hnew = sk2cart(H, U, V)
Hnew2 = sk2cart_num(H, U, V, FHIaims=true)
println(@test Hnew ≈ Hnew2)

#@info("cart2sk vs cart2sk_num : spspd")
#orbitals = [sko"s", sko"p", sko"s", sko"p", sko"d"]
#H = SKH(orbitals) # generate bonds automatically
#println(@test H == SKH("spspd"))
#H = SKH("spspd")
#HE = alloc_block(H)
#for i=1:max_locidx(H)
#   for j=i:max_locidx(H)
#      HE[i,j] = rand(1)[1]
#      if i != j
#         HE[j,i] = HE[i,j]
#      end
#   end
#end
#U = rand(3) .- 0.5
#U /= norm(U)
#V = cart2sk(H, U, HE)
#Vnew = cart2sk_num(H, U, HE)
#println(@test V ≈ Vnew)

#@info("cart2sk vs sk2cart : spspd")
#orbitals = [sko"s", sko"p", sko"s", sko"p", sko"d"]
#H = SKH(orbitals) # generate bonds automatically
#println(@test H == SKH("spspd"))
#H = SKH("spspd")
#HE = alloc_block(H)
#for i=1:max_locidx(H)
#   for j=i:max_locidx(H)
#      HE[i,j] = rand(1)[1]
#      if i != j
#         HE[j,i] = HE[i,j]
#      end
#   end
#end
#U = rand(3) .- 0.5
#U /= norm(U)
#V = cart2sk(H, U, HE)
#Hnew = sk2cart(H, U, V)
#println(@test HE ≈ Hnew)

end


text:
#path to data

#read velocity file
#filename = "C:\\Users\\Belal\\.julia\\dev\\My_iesh.jl\\data\\surface_Au111.dat"
#filename = "/home/razamaza/dev/My_iesh.jl/data/surface_Au111.dat"
filename = "C:\\Users\\braza2\\Documents\\GitHub\\My_iesh.jl\\data\\surface_Au111.dat"
f = open(filename)
data = readlines(f)
close(f)
global const Å = 1.0e-10
const global N = Int64(parse(Float64, data[1]))
const global g = parse(Float64, data[2]) * Å
const global a = sqrt(2) * g
cell_temp = zeros(Float64, 3)
cell_temp[1] = parse(Float64, data[3])
cell_temp[2] = parse(Float64, data[4])
cell_temp[3] = sqrt(6) * parse(Float64, data[2])* Float64(100)
const global cell = SVector{3, Float64}(copy(cell_temp)) * Å

x_au_temp = zeros(Float64, (Integer(N), 3))

for i in 5:length(data)
    data_strip = strip(data[i])
    data_split = split(data_strip)
    str_to_float = parse.(Float64, data_split)
    x_au_temp[i-4, :] = str_to_float
end
#x_no_temp = vcat(x_no, x_temp)
const global x_au0 = copy(x_au_temp)*Å

#filename = "C:\\Users\\Belal\\.julia\\dev\\My_iesh.jl\\data\\gold_velocities.dat" #units in meter per second
#filename = "/home/razamaza/dev/My_iesh.jl/data/gold_velocities.dat"
filename ="C:\\Users\\braza2\\Documents\\GitHub\\My_iesh.jl\\data\\gold_velocities.dat"
f = open(filename)
data = readlines(f)
close(f)

v_au_temp = zeros(Float64, (Integer(N), 3))

for i in 1:length(data)
    data_strip = strip(data[i])
    data_split = split(data_strip)
    str_to_float = parse.(Float64, data_split)
    v_au_temp[i, :] = str_to_float
end
#x_no_temp = vcat(x_no, x_temp)
const global v_au0 = copy(v_au_temp)


text:
---------------------------------------------
# the Rosenbrock-Lerner polynomial
N=60
b=15
@polyvar x[1:N]
f_R=10*sum((x[3:N-1].+x[2:N-2].-x[1:N-3].^2).^2)+sum((ntuple(i->1,N-3).-x[1:N-3].-x[4:N]).^2)
f_Q=sum([x[i]*i/j*x[j] for i=1:b,j=1:b])
f=f_R+f_Q
---------------------------------------------

---------------------------------------------
# the Broyden tridiagonal polynomial
l=5
p=20
n=l*p
@polyvar x[1:n]
f=((3-2*x[1])*x[1]-2*x[2]+1)^2
for i=2:n-1
    global f+=((3-2*x[i])*x[i]-x[i-1]-2*x[i+1]+1)^2
end
f+=((3-2*x[n])*x[n]-x[n-1]+1)^2
pop=[f]
for i=1:l
    push!(pop,1-sum(x[(i-1)*p+1:i*p].^2))
end
---------------------------------------------

---------------------------------------------
# the Chained singular polynomial
l=2
p=15
n=l*p
@polyvar x[1:n]
f=(x[1]+10*x[2])^2+5*(x[3]-x[4])^2+(x[2]-2*x[3])^4+10*(x[1]-10*x[4])^4
for i=3:2:n-3
    global f+=(x[i]+10*x[i+1])^2+5*(x[i+2]-x[i+3])^2+(x[i+1]-2*x[i+2])^4+10*(x[i]-10*x[i+3])^4
end
pop=[f]
for i=1:l
    push!(pop,1-sum(x[(i-1)*p+1:i*p].^2))
end
---------------------------------------------

---------------------------------------------
# the Broyden banded polynomial
n=400
@polyvar x[1:n]
f=x[1]
for i=1:n
    jset=max(1,i-5):min(n,i+1)
    jset=setdiff(jset,i)
    f0=sum([(1+x[jset[j]])*x[jset[j]] for j=1:length(jset)])
    global f+=(x[i]*(2+5*x[i]^2)+1-f0)^2
end
f-=x[1]
---------------------------------------------

---------------------------------------------
# the chained Wood polynomial
l=50
p=20
n=l*p
@polyvar x[1:n]
f=x[1]+1
for i=1:2:n-3
    global f+=100*(x[i+1]-x[i]^2)^2+(1-x[i])^2+90*(x[i+3]-x[i+2]^2)^2+(1-x[i+2])^2+10*(x[i+1]+x[i+3]-2)^2+0.1*(x[i+1]-x[i+3])^2
end
f-=x[1]
pop=[f]
for i=1:l
    push!(pop,1-sum(x[(i-1)*p+1:i*p].^2))
end
---------------------------------------------

---------------------------------------------
# the generalized Rosenbrock polynomial
l=50
p=20
n=l*p
@polyvar x[1:n]
f=x[1]+1
for i=2:n
    global f+=100*(x[i]-x[i-1]^2)^2+(1-x[i])^2
end
f-=x[1]
pop=[f]
for i=1:l
    push!(pop,1-sum(x[(i-1)*p+1:i*p].^2))
end
---------------------------------------------


text:
@info "Attempting to import PredictMD..."
@time import PredictMDFuture
@info "Successfully imported PredictMD"

using Test

import Aqua
import Literate
import Pkg
import TOML

const predictmd_directory = dirname(dirname(pathof(PredictMDFuture)))
const predictmd_projectfile = joinpath(predictmd_directory, "Project.toml")
const predictmd_project = TOML.parsefile(predictmd_projectfile)
const predictmd_direct_deps = predictmd_project["deps"]
const predictmd_project_compat = predictmd_project["compat"]

@info "" predictmd_directory

@time @testset "PredictMD.jl" begin
    include("aqua.jl")
    include("compat_entries.jl")

    include("package_tests.jl")
end


text:
# Top type in hierarchy of AST transformers.
abstract ASTTransformer

# ASTTransformers have methods to transform each type of Expr block. Their default
# methods do no transformation, they just leave the tree untouched.
[]

# Default transformer is identity, i.e. doesn't transform at all, but traverses into child trees.
function transform{AT <: ASTTransformer}(at::AT, ast::Expr)
	for i in 1:numsubtrees(ast)
		setchild(ast, i, transform(at, child(ast, i))
	end
	ast
end

function child(ast::Expr, i)
	if i == 1 && in(ast.head, [:call])
		ast

type DictDefinedASTTransformer <: ASTTransformer

end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

# Common update and digest functions which work across SHA1 and SHA2

# update! takes in variable-length data, buffering it into blocklen()-sized pieces,
# calling transform!() when necessary to update the internal hash state.
function update!(context::T, data::U, datalen=length(data)) where {T<:SHA_CTX, U<:AbstractBytes}
    # We need to do all our arithmetic in the proper bitwidth
    UIntXXX = typeof(context.bytecount)

    # Process as many complete blocks as possible
    0 ≤ datalen ≤ length(data) || throw(BoundsError(data, firstindex(data)+datalen-1))
    len = convert(UIntXXX, datalen)
    data_idx = convert(UIntXXX, firstindex(data)-1)
    usedspace = context.bytecount % blocklen(T)
    while len - data_idx + usedspace >= blocklen(T)
        # Fill up as much of the buffer as we can with the data given us
        copyto!(context.buffer, usedspace + 1, data, data_idx + 1, blocklen(T) - usedspace)

        transform!(context)
        context.bytecount += blocklen(T) - usedspace
        data_idx += blocklen(T) - usedspace
        usedspace = convert(UIntXXX, 0)
    end

    # There is less than a complete block left, but we need to save the leftovers into context.buffer:
    if len > data_idx
        copyto!(context.buffer, usedspace + 1, data, data_idx + 1, len - data_idx)
        context.bytecount += len - data_idx
    end
end

# Pad the remainder leaving space for the bitcount
function pad_remainder!(context::T) where T<:SHA_CTX
    usedspace = context.bytecount % blocklen(T)
    # If we have anything in the buffer still, pad and transform that data
    if usedspace > 0
        # Begin padding with a 1 bit:
        context.buffer[usedspace+1] = 0x80
        usedspace += 1

        # If we have room for the bitcount, then pad up to the short blocklen
        if usedspace <= short_blocklen(T)
            for i = 1:(short_blocklen(T) - usedspace)
                context.buffer[usedspace + i] = 0x0
            end
        else
            # Otherwise, pad out this entire block, transform it, then pad up to short blocklen
            for i = 1:(blocklen(T) - usedspace)
                context.buffer[usedspace + i] = 0x0
            end
            transform!(context)
            for i = 1:short_blocklen(T)
                context.buffer[i] = 0x0
            end
        end
    else
        # If we don't have anything in the buffer, pad an entire shortbuffer
        context.buffer[1] = 0x80
        for i = 2:short_blocklen(T)
            context.buffer[i] = 0x0
        end
    end
end

# Clear out any saved data in the buffer, append total bitlength, and return our precious hash!
# Note: SHA3_CTX has a more specialised method
function digest!(context::T) where T<:SHA_CTX
    pad_remainder!(context)
    # Store the length of the input data (in bits) at the end of the padding
    bitcount_idx = div(short_blocklen(T), sizeof(context.bytecount)) + 1
    pbuf = Ptr{typeof(context.bytecount)}(pointer(context.buffer))
    unsafe_store!(pbuf, bswap(context.bytecount * 8), bitcount_idx)

    # Final transform:
    transform!(context)

    # Return the digest
    return reinterpret(UInt8, bswap!(context.state))[1:digestlen(T)]
end

max_stars_repo_name: greimel/julia

max_issues_repo_name: greimel/julia

max_forks_repo_name: greimel/julia


text:
function ode(f, initcondition, timespan, alg=Tsit5())
    prob = ODEProblem(f,initcondition,timespan)
    return solve(prob,alg,reltol=1e-8,abstol=1e-8)
end

function plotode(sol)
    plot(sol,title="Solution to ODE")
end

function solve_plot_print(f, initcondition, timespan, alg=Tsit5(), isprint=true)
   sol =  ode(f, initcondition, timespan, alg)
   plotode(sol)

end

export ode,plotode, solve_plot_print

max_stars_repo_name: bkgsur/DifferentialEquationsWithJulia

max_issues_repo_name: bkgsur/DifferentialEquationsWithJulia

max_forks_repo_name: bkgsur/DifferentialEquationsWithJulia


text:
module Interpolation

    using DocStringExtensions
    using StaticArrays
    using LinearAlgebra

    using WaveProp.Utils
    using WaveProp.Geometry

    export
        barycentric_lagrange_matrix,
        barycentric_lagrange_weights,
        MonomialBasis,
        vandermond,
        lagrange_basis,
        PolynomialBasis,
        LagrangeBasis, gradLagrangeBasis,
        trace, grad

    include("lagrange.jl")
    include("monomial.jl")
    include("polynomialbasis.jl")

end

max_stars_repo_name: WaveProp/WaveProp

max_issues_repo_name: WaveProp/WaveProp

max_forks_repo_name: WaveProp/WaveProp


text:
#
# This file is part of the Actors.jl Julia package,
# MIT license, part of https://github.com/JuliaActors
#

include("delays.jl")

using Actors, Test, .Threads, .Delays
import Actors: spawn, info, diag, newLink

t1 = Ref{Task}()
t2 = Ref{Task}()
t3 = Ref{Task}()

sv = supervisor(taskref=t1)
a1 = diag(sv, :act)
@test t1[].state == :runnable
@test sv.mode == :supervisor
@test a1.mode == :supervisor
@test a1.bhv isa Actors.Supervisor
@test a1.bhv.option[:strategy] == :one_for_one
@test isempty(which_children(sv))

act2 = spawn(threadid, taskref=t2)
exec(act2, supervise, sv)
a2 = diag(act2, :act)
@test a2.conn[1] isa Actors.Super
@test a2.conn[1].lk == sv
@test a1.conn[1] isa Actors.Child
@test a1.conn[1].lk == act2
@test isnothing(a1.conn[1].start)
@test a1.conn[1].info.restart == :transient
@test isempty(a1.bhv.rtime)
@test length(which_children(sv)) == 1

println("Testing supervision with failures:")
# test actor restart
oldch2 = act2.chn
oldtsk = t2[]
# first failure
send(act2, "boom")
a2 = diag(changed(act2), :act)
@test t2[].state == :failed
@test t1[].state == :runnable
@test oldch2 != act2.chn
t2[] = diag(act2, :task)
@test oldtsk != t2[]
@test t2[].state == :runnable
@test length(a1.bhv.rtime) == 1
@test oldtsk == diag(sv, :err)[1]
@test a1.conn[1] isa Actors.Child
@test a1.conn[1].lk == act2
@test a1.bhv.childs[1].lk == act2
@test a2.conn[1].lk == sv
# second failure
send(act2, "boom")
changed(act2)
@test t1[].state == :runnable
@test t2[].state == :failed
t2[] = diag(act2, :task)
@test length(a1.bhv.rtime) == 2
# third failure
send(act2, "boom")
changed(act2)
@test t1[].state == :runnable
# fourth failure
send(act2, "boom")
sleep(2)
@test @delayed t1[].state == :done 2

# temporary actors
sv = supervisor(taskref=t1)
act2 = spawn(threadid, taskref=t2)
supervise(sv, act2, restart=:temporary)
a2 = diag(act2, :act)
a1 = diag(sv, :act)
@test t1[].state == :runnable
@test a1.conn[1].lk == act2
@test a1.conn[1].info.restart == :temporary
send(act2, "boom")
@test @delayed isempty(a1.bhv.childs)
@test isempty(a1.conn)
@test diag(sv, :err)[1] == t2[]
@test t1[].state == :runnable

# test task supervision
tvar = [0]
function ttask(tvar, delay, fail=false)
    sleep(delay)
    tvar[1] += 1
    fail && error("Task test error!")
end

println("supervisor restart!")
sv = supervisor(taskref=t1)
rt = start_task(()->ttask(tvar, 1), sv)
a1 = diag(sv, :act)
@test @delayed a1.conn[1] isa Actors.Child
@test a1.conn[1].lk == rt
@test a1.bhv.childs[1].lk == rt
@test @delayed isempty(a1.conn) 2
@test isempty(a1.bhv.childs)
@test rt[].state == :done
@test tvar[1] == 1
@test t1[].state == :runnable

# test task restart
rt = start_task(()->ttask(tvar, 0.2, true), sv)  # with failure
changed(rt)
@test @delayed t1[].state == :runnable
# 4 task errors
@test @delayed t1[].state == :done 3
@test @delayed tvar[1] == 5

# test restart strategies
function checkTasks(actors, tnr, equ)
    for i in eachindex(actors)
        t = diag(actors[i], :task)
        @test t.state == :runnable
        @test equ[i](tnr[i], convert(UInt, pointer_from_objref(t)))
    end
end
ctask = [Ref{Task}() for _ in 1:5]
sv = supervisor(:one_for_all, taskref=t1)
cact = [spawn(threadid, taskref=ctask[i]) for i in 1:5]
foreach(act->exec(act, supervise, sv), cact)
a1 = diag(sv, :act)
@test @delayed length(a1.bhv.childs) == 5
@test length(a1.conn) == 5
ptask = map(t->convert(UInt, pointer_from_objref(t[])), ctask)
send(cact[3], "boom")
changed(cact[3])
checkTasks(cact, ptask, fill(!=,5))
@test length(a1.bhv.childs) == 5
@test length(a1.conn) == 5

set_strategy(sv, :rest_for_one)
ptask = map(a->convert(UInt, pointer_from_objref(diag(a,:task))), cact)
send(cact[3], "boom")
changed(cact[3])
@test a1.bhv.option[:strategy] == :rest_for_one
checkTasks(cact, ptask, (==,==,!=,!=,!=))

# API
ch = Channel(0)
rt = start_task(()->take!(ch), sv, timeout=Inf)
@test @delayed length(which_children(sv)) == 6
@test length(which_children(sv, true)) == 6
unsupervise(sv, cact[5])
unsupervise(sv, rt)
put!(ch, 0)
@test @delayed length(which_children(sv)) == 4
@test Actors.diag(cact[5], :task).state == :runnable
act2 = start_actor(sv) do
    threadid()
end
t = start_task(sv, timeout=Inf) do
    take!(ch)
end
@test @delayed length(which_children(sv)) == 6
terminate_child(sv, act2)
@test @delayed length(which_children(sv)) == 5
@test info(act2) == :done
@test count_children(sv).all == 5

# supervisor shutdown
sv = supervisor(taskref=t1)
act2 = spawn(threadid, taskref=t2)
supervise(sv, act2)
sleep(0.5)
exit!(sv, :shutdown)
@test @delayed t1[].state == :done
@test t2[].state == :done


text:
using VoronoiCells
import VoronoiCells: LEFT, RIGHT, LOWER, UPPER
using VoronoiDelaunay
using Base.Test

# ------------------------------------------------------------
# Points inside the bounding box are untouched by clip

A = Point2D(LEFT+rand(), LOWER+rand())
B = Point2D(LEFT+rand(), LOWER+rand())
C, D = clip(A, B)

@test A == C
@test B == D

# ------------------------------------------------------------
# With one point on the boundary and the other outside (in the correct
# direction), clip should return two identical points

A = Point2D(LEFT, LOWER+rand())
B = Point2D(LEFT-0.5, LOWER+rand())
C, D = clip(A, B)
@test C == D

A = Point2D(LEFT+rand(), UPPER)
B = Point2D(LEFT+rand(), UPPER+0.5)
C, D = clip(A, B)
@test C == D

# ------------------------------------------------------------
# When clipping off part of a line segment the slope should not change

A = Point2D(LEFT+rand(), LOWER+rand())
B = Point2D(LEFT+rand(), UPPER+rand())
C, D = clip(A, B)

slope1 = (gety(A)-gety(B)) / (getx(A)-getx(B))
slope2 = (gety(C)-gety(D)) / (getx(C)-getx(D))

@test_approx_eq slope1 slope2


text:
# downsample.jl
# Copyright 2019-03-05, Jeff Fessler, University of Michigan

using Printf
using Test

"""
`y = downsample1(x, down; warn=true)`
downsample by factor m along first dimension by averaging

in
* x	[n1 (Nd)]
* down			integer downsampling factor

option
`warn`	warn if noninteger multiple; default true

out
* y	[n1/down (Nd)]

Copyright 2019-03-05, Jeff Fessler, University of Michigan
"""
function downsample1(x, down::Integer; warn::Bool=true)

	dim = size(x)
	dim1 = dim[1]
	x = reshape(x, dim1, :) # [n1 *Nd]
	m1 = Int(floor(dim1 / down))
	if m1 * down < dim1
		if warn
			@warn(@sprintf("truncating input size %d to %d", dim1, m1 * down))
		end
		x = x[1:(m1*down),:]
	end
	y = reshape(x, down, :)
	y = sum(y, dims=1) / down # mean
	y = reshape(y, m1, dim[2:end]...)
	return y
end

"""
downsample1_test()
"""
function downsample1_test()
	x = reshape(2:2:48, 4, 6)
	y = downsample1(x, 2)
	@test y == reshape(3:4:47, 2, 6)
	true
end

"""
`y = downsample2(x, down, warn=true)`

downsample by averaging by integer factors
in
* `x` [nx ny]
* `down` can be a scalar (same factor for both dimensions) or a 2-vector

option
* `warn`	warn if noninteger multiple; default true

out
`y`	[nx/down ny/down]
"""
function downsample2(x, down; warn::Bool=true)
	fun = (x, d) -> downsample1(x, d, warn=warn)
	y = fun(x, down[1])
	y = fun(y', down[end])'
	return y
end

"""
`downsample2_test()`
"""
function downsample2_test()
	x = reshape(1:24, 4, 6)
	y = downsample2(x, 2)
	@test y ==  [3.5 11.5 19.5; 5.5 13.5 21.5]
	true
end

"""
`y = downsample3(x, down, warn=true)`

downsample by averaging by integer factors
in
* `x` [nx ny nz]
* `down` can be a scalar (same factor for all dimensions) or a 3-vector

option
* `warn`	warn if noninteger multiple; default true

out
`y`	[nx/down ny/down nz/down]
"""
function downsample3(x, down; warn::Bool=true)

	if ndims(x) == 2
		@warn("2d case")
		if length(down) == 1; down = [down, down, 1]; end
		length(down) != 3 & throw(DimensionMismatch("$(length(down)) != 3"))&
		return downsample2(x, down[1:2])
	end

	ndims(x) != 3 && throw(DimensionMismatch("ndims(x)=$(ndims(x)) != 3"))

	if length(down) == 1
		down = [down,down,down]
	end

	length(down) != ndims(x) && throw(DimensionMismatch("length(down) != ndims(x)"))

	# downsample along each dimension
	y = downsample1(x, down[1])
	y = downsample1(permutedims(y, [2, 1, 3]), down[2])
	y = downsample1(permutedims(y, [3, 2, 1]), down[3]) # [3 1 2] order
	y = permutedims(y, [2, 3, 1])
	return y
end

function downsample3_test()
	x = [6, 5, 2]
	x = reshape(2*(1:prod(x)), x...)
	for down = 1:2
		y = downsample3(x, down)
		if down == 1
			@test y == x
		elseif down == 2
			@test y == reshape([39 63; 43 67; 47 71], 3,2,1) # squeeze
		end
	end
	true
end

#=
if has_aspire # todo
	filex = [test_dir filesep 'testx.fld']
	filey = [test_dir filesep 'testy.fld']
	fld_write(filex, x)
	delete(filey)
	com = ['op sample3 mean ' filey ' ' filex ' %d %d %d']
	com = sprintf(com, down, down, down)
	os_run(com)
	z = fld_read(filey)
	!isequal(y, z) && throw("aspire/matlab mismatch")
end
=#

"""
`downsample(:test)` self test
"""
function downsample(test::Symbol)
	test != :test && throw(ArgumentError("test $test"))
	downsample1_test()
	downsample2_test()
	downsample3_test()
	true
end


text:
# This file is a part of DataDeps.jl. License is MIT.

abstract type AbstractDataDep end

"""
    ManualDataDep(name, message)

A DataDep for if the installation needs to be handled manually.
This can be done via Pkg/git if you put the dependency into the packages repo's `/deps/data` directory.
More generally, message should give instructions on how to setup the data.
"""
struct ManualDataDep <: AbstractDataDep
    name::String
    message::String
end

struct DataDep{H, R, F, P} <: AbstractDataDep
    name::String
    remotepath::R
    hash::H
    fetch_method::F
    post_fetch_method::P
    extra_message::String
end

"""
```
DataDep(
    name::String,
    message::String,
    remote_path::Union{String,Vector{String}...},
    [checksum::Union{String,Vector{String}...},]; # Optional, if not provided will generate
    # keyword args (Optional):
    fetch_method=fetch_default # (remote_filepath, local_directory)->local_filepath
    post_fetch_method=identity # (local_filepath)->Any
)
```

#### Required Fields

 - *Name**: the name used to refer to this datadep, coresponds to a folder name where it will be stored
    - It can have spaces or any other character that is allowed in a Windows filestring (which is a strict subset of the restriction for unix filenames).
 - *Message*: A message displayed to the user for they are asked if they want to downloaded it.
    - This is normally used to give a link to the original source of the data, a paper to be cited etc.
 - *remote_path*: where to fetch the data from. Normally a string or strings) containing an URL
    - This is usually a string, or a vector of strings (or a vector of vector... see [Recursive Structure](Recursive Structure) below)

#### Optional Fields
 - *checksum* this is very flexible, it is used to check the files downloaded correctly
    - By far the most common use is to just provide a SHA256 sum as a hex-string for the files
    - If not provided, then a warning message with the  SHA256 sum is displayed. This is to help package devs work out the sum for their files without using an external tool.
    - If you want to use a different hashing algorithm, then you can provide a tuple `(hashfun, targethex)`
        - `hashfun` should be a function which takes an IOStream, and returns a `Vector{UInt8}`.
          - Such as any of the functions from [SHA.jl](https://github.com/staticfloat/SHA.jl), eg `sha3_384`, `sha1_512`
          - or `md5` from [MD5.jl](https://github.com/oxinabox/MD5.jl)
  - If you want to use a different hashing algorithm, but don't know the sum, you can provide just the `hashfun` and a warning message will be displayed, giving the correct tuple of `(hashfun, targethex)` that should be added to the registration block.

    - If you don't want to provide a checksum,  because your data can change pass in the type `Any` which will suppress the warning messages. (But see above warnings about "what if my data is dynamic")
    - Can take a vector of checksums, being one for each file, or a single checksum in which case the per file hashes are `xor`ed to get the target hash. (See [Recursive Structure](Recursive Structure) below)

 -  `fetch_method=fetch_default` a function to run to download the files.
    - Function should take 2 parameters (remotepath, local_directory), and must return a local filepath
    - It is responsible for determining what the local filename should be
    - Change this to change the transfer protocol, for example to use an auth'ed connection.
    - Default `fetch_default` which fully supports HTTP, and has fallbacks to support any type which overloads `Base.basename` and `Base.download` (see [`fetch_base`](@ref))
    - Can take a vector of methods, being one for each file, or a single method, in which case that method is used to download all of them. (See [Recursive Structure](Recursive Structure) below)
    - Very few people will need to override this if they are just downloading public HTTP files.

 - `post_fetch_method` a function to run after the files have download
    - Should take the local filepath as its first and only argument. Can return anything.
    - Default is to do nothing.
    - Can do what it wants from there, but most likes wants to extract the file into the data directory.
    - towards this end DataDeps includes a command: `unpack` which will extract an compressed folder, deleting the original.
    - It should be noted that it `post_fetch_method` runs from within the data directory
       - which means operations that just write to the current working directory (like `rm` or `mv` or ```run(`SOMECMD`))``` just work.
       - You can call `cwd()` to get the the data directory for your own functions. (Or `dirname(local_filepath)`)
    - Can take a vector of methods, being one for each file, or a single method, in which case that ame method is applied to all of the files. (See **Recursive Structure** in the README.md)
"""
function DataDep(name::String,
                 message::String,
                 remotepath, hash=nothing;
                 fetch_method=fetch_default,
                 post_fetch_method=identity)

    DataDep(name, remotepath, hash, fetch_method, post_fetch_method, message)
end


text:
# Common utilities

##### common types

abstract ClusteringResult

# generic functions

nclusters(R::ClusteringResult) = length(R.counts)
counts(R::ClusteringResult) = R.counts
assignments(R::ClusteringResult) = R.assignments

##### convert weight options

conv_weights{T}(::Type{T}, n::Int, w::Nothing) = nothing

function conv_weights{T}(::Type{T}, n::Int, w::Vector)
    length(w) == n || throw(DimensionMismatch("Incorrect length of weights."))
    convert(Vector{T}, w)::Vector{T}
end

##### convert display symbol to disp level

display_level(s::Symbol) =
    s == :none ? 0 :
    s == :final ? 1 :
    s == :iter ? 2 :
    error("Invalid value for the option 'display'.")

##### update minimum value

function updatemin!(r::AbstractArray, x::AbstractArray)
    n = length(r)
    length(x) == n || throw(DimensionMismatch("Inconsistent array lengths."))
    @inbounds for i = 1:n
        xi = x[i]
        if xi < r[i]
            r[i] = xi
        end
    end
    return r
end


text:
include("..\\libs\\FunctionsRobot.jl")

function moveAroundAndPut!(r::Robot)#, side::HorizonSide)
    for side in instances(HorizonSide)
        putmarker!(r)
        marksLine!(r, side)
    end
end

function perimetrWithPartitions!(r::Robot)
    sides = moveAndReturnDirections!(r)
    moveAroundAndPut!(r)
    moveToStartplace!(r)
    moveToBeginplace!(r, sides)
end


text:
bitstring(x::Posit8) = bitstring(reinterpret(UInt8,x))
bitstring(x::Posit16) = bitstring(reinterpret(UInt16,x))
bitstring(x::Posit32) = bitstring(reinterpret(UInt32,x))

bitstring(x::Posit8_1) = bitstring(UInt8(reinterpret(UInt32,x) >> 24))
bitstring(x::Posit8_2) = bitstring(UInt8(reinterpret(UInt32,x) >> 24))

bitstring(x::Posit16_1) = bitstring(UInt16(reinterpret(UInt32,x) >> 16))
bitstring(x::Posit16_2) = bitstring(UInt16(reinterpret(UInt32,x) >> 16))

# cut off the string as UInt24 does not exist
bitstring(x::Posit24_1) = bitstring(reinterpret(UInt32,x))[1:24]
bitstring(x::Posit24_2) = bitstring(reinterpret(UInt32,x))[1:24]

"""Bitstring split into sign, regime, exponent and fraction bits.
# Example
julia> bitstring(Posit16(123)," ")
"0 11110 0 111011000"
"""
function bitstring(x::AbstractPosit,mode::Symbol)
    if mode == :split
        # number of exponent bits
        ne = if x isa PositEx1 1 elseif x isa PositEx2 2 else 0 end
        s = bitstring(x)
        n = n_regimebits(s)

        # concatenate the string
        if n >= length(s)-2     # only exponent bits (possibly including terminating bit)
            return "$(s[1]) $(s[2:end])"
        elseif n > length(s)-ne-3
            return "$(s[1]) $(s[2:2+n]) $(s[2+n+1:end])"
        else
            if ne == 0
                return "$(s[1]) $(s[2:2+n]) $(s[2+n+1:end])"
            else
                return "$(s[1]) $(s[2:2+n]) $(s[2+n+1:2+n+ne]) $(s[2+n+ne+1:end])"
            end
        end
    else
        return bitstring(x)
    end
end

"""Number of regime bits for a posit bitstring, excluding the terminating bit.
# Example
julia> n_regimebits("0111100111011000")
4"""
function n_regimebits(s::String)
    n = 2   # skip the sign bit
    slength = length(s)
    while n < slength && s[n+1] == s[n]
        n += 1
    end
    n-1     # subtract the sign bit
end

function Base.show(io::IO, x::AbstractPosit)
    if isnan(x)
        print(io, "NaR")
    else
		io2 = IOBuffer()
        print(io2,Float32(x))
        f = String(take!(io2))
        print(io,string(typeof(x))*"("*f*")")
    end
end


text:
export  NeuralNetwork2D, getStress

mutable struct NeuralNetwork2D
    ρ::Float64 # density
    # hardening parameter, yield function = f - (σY + Kα)
    σ0::Array{Float64} # stress at last time step
    σ0_::Array{Float64} # σ0 to be updated in `commitHistory`
    ε0::Array{Float64}
    ε0_::Array{Float64}
    α0::Union{Missing, Array{Float64}}
    α0_::Union{Missing, Array{Float64}}
    nn::Function
end

function NeuralNetwork2D(prop::Dict{String, Any})
    ρ = prop["rho"];
    nn = prop["nn"]
    α0 = missing; α0_ = missing
    if haskey(prop, "n_internal")
        n_internal = prop["n_internal"]
        @assert isa(n_internal, Integer) && n_internal>0
        α0 = zeros(n_internal); α0_ = zeros(n_internal)
        # @show α0
    end
    σ0 = zeros(3); σ0_ = zeros(3); ε0 = zeros(3); ε0_ = zeros(3);
    NeuralNetwork2D(ρ, σ0, σ0_,ε0,ε0_, α0, α0_, nn)
end

function getStress(self::NeuralNetwork2D,  strain::Array{Float64},  Dstrain::Array{Float64}, Δt::Float64 = 0.0)
    # #@show "***", strain, Dstrain
    ε = strain
    ε0 = Dstrain
    σ0 = self.σ0
    if ismissing(self.α0_)
        σ, dΔσdΔε = self.nn(ε, ε0, σ0, Δt)
        self.σ0_ = σ
        self.ε0_ = ε
        return σ, dΔσdΔε
    else
        # @show self.α0
        σ, dΔσdΔε, α = self.nn(ε, ε0, σ0, Δt, self.α0)
        self.σ0_ = σ
        self.ε0_ = ε
        self.α0_ = α
        return σ, dΔσdΔε
    end
end

function getTangent(self::NeuralNetwork2D)
    error("Not implemented")
end

function commitHistory(self::NeuralNetwork2D)
    self.σ0 = self.σ0_
    self.ε0 = self.ε0_
    self.α0 = self.α0_
end


text:
using Gadfly, Colors, BlackBoxOptimizationBenchmarking
using Distributed, Statistics
##

run_bench = true

nprocs()< 2 && addprocs(2)

## Setup

bbob_dir = dirname(pathof(BlackBoxOptimizationBenchmarking))

@everywhere begin

    using BlackBoxOptimizationBenchmarking
    const BBOB = BlackBoxOptimizationBenchmarking
    bbob_dir = dirname(pathof(BBOB))
    include(joinpath(bbob_dir,"../scripts/optimizers_interface.jl"))

    dimensions = [3 6 12]
    Ntrials = 1
    Δf = 1e-6
    run_lengths = round.(Int,range(20, stop=3_000, length=15))
    #run_lengths = round.(Int,linspace(20,20_000,15))
    funcs = BBOB.list_functions()
end

c = opt -> Chain(opt,NelderMead(),0.9)

optimizers = [
    c(SAMIN(;verbosity=0)),
    NelderMead(),
    OptimRestart(NelderMead()),
    c(SimulatedAnnealing()),
    c(BlackBoxOptimMethod(:adaptive_de_rand_1_bin_radiuslimited)),
    c(BlackBoxOptimMethod(:adaptive_de_rand_1_bin)),
    c(NLoptOptimMethod(:GN_ISRES)),
    c(NLoptOptimMethod(:GN_ESCH)),
    c(NLoptOptimMethod(:GD_STOGO)),
    c(BlackBoxOptimMethod(:xnes)),
    c(BlackBoxOptimMethod(:generating_set_search)),
    c(BlackBoxOptimMethod(:de_rand_2_bin)),
    #BlackBoxOptimMethod(:resampling_memetic_search),#poor performance
    PyMinimize("Nelder-Mead"),
    #PyCMA(),
]

opt_strings = map(string,optimizers)

## run benchmark

if run_bench

    results = pmap(opt->BBOB.benchmark(opt, funcs, run_lengths, Ntrials, dimensions, Δf), optimizers)
    error("done")
    @info("Done, saving data and making plots.")

    # save output
    outdir = joinpath(bbob_dir,"..","data")

    writedlm(joinpath(outdir,"mean_succ.txt"),mean_succ)
    writedlm(joinpath(outdir,"mean_dist.txt"),mean_dist)
    writedlm(joinpath(outdir,"mean_fmin.txt"),mean_fmin)
    writedlm(joinpath(outdir,"runtime.txt"),runtime)
else
    # load output
    outdir = joinpath(bbob_dir,"..","data")

    s = (length(optimizers),length(funcs),length(run_lengths),length(dimensions))
    mean_succ = reshape( readdlm(joinpath(outdir,"mean_succ.txt")), s)
    mean_dist = reshape( readdlm(joinpath(outdir,"mean_dist.txt")), s)
    mean_fmin = reshape( readdlm(joinpath(outdir,"mean_fmin.txt")), s)
    runtime = readdlm(joinpath(outdir,"runtime.txt"))
end

## make plots

outdir = joinpath(bbob_dir,"..","data","plots")
cols = Colors.distinguishable_colors(size(mean_succ,1)+1,colorant"white")[2:end]
line_style = i -> mod(i,2)==0 ? [:solid] : [:dash]

#all together

idx = sortperm([mean(mean_succ[i,:,end,:]) for i=1:length(optimizers)],rev=true)

p = plot(
    [layer(
        x=collect(run_lengths),y=mean(mean_succ[i,:,:,:],dims=(1,3)),Geom.line,
        Theme(default_color=cols[i], line_style=line_style(i), line_width=2pt)
    ) for i in 1:size(mean_succ,1)]...,
    Coord.cartesian(ymax=1),
    Guide.title("All functions"), Guide.xlabel("Run Length"), Guide.ylabel("Success rate"),
    Guide.manual_color_key("", opt_strings[idx], cols[idx])
)

#

#draw(PDF(joinpath(outdir,"mean_succ.pdf"),22cm,14cm),p)
#draw(PNG(joinpath(outdir,"mean_succ.png"),dpi=150,16cm,10cm),p)
draw(SVG(joinpath(outdir,"mean_succ.svg"),16cm,10cm),p)

# per dimension
for (k,D) in zip(1:length(dimensions), dimensions)

    idx = sortperm([mean(mean_succ[i,:,end,k]) for i=1:length(optimizers)],rev=true)

    p = plot(
        [layer(
            x=collect(run_lengths),y=mean(mean_succ[i,:,:,k],dims=1),Geom.line,
            Theme(default_color=cols[i], line_style=line_style(i), line_width=2pt)
        ) for i in 1:size(mean_succ,1)]...,
        Coord.cartesian(ymax=1),
        Guide.title("All functions, D: $D"), Guide.xlabel("Run Length"), Guide.ylabel("Success rate"),
        Guide.manual_color_key("", opt_strings[idx], cols[idx])
    )
    #draw(PDF(joinpath(outdir,"per_dimension","mean_succ_$(D).pdf"),22cm,14cm),p)
    #draw(PNG(joinpath(outdir,"per_dimension","mean_succ_$(D).png"),dpi=150,16cm,10cm),p)
    draw(SVG(joinpath(outdir,"per_dimension","mean_succ_$(D).svg"),16cm,10cm),p)

end

## f min

idx = sortperm([exp(mean(log.(mean_fmin[i,:,end,:].+eps()))) for i=1:length(optimizers)],rev=true)

# all together
p = plot(
    [layer(
        x=collect(run_lengths),y=exp.(mean(log.(mean_fmin[i,:,:,:].+eps()),dims=(1,3))),Geom.line,
        Theme(default_color=cols[i],line_width=2pt)
    )
    for i in 1:size(mean_succ,1)]...,
    layer(yintercept=[Δf],Geom.hline,Theme(default_color=colorant"gray")),
    Guide.title("All functions"), Guide.xlabel("Run Length"), Guide.ylabel("fmin"),
    Guide.manual_color_key("", opt_strings[idx], cols[idx]),
#    Scale.x_log10,
    Scale.y_log10,
)
#draw(PDF(joinpath(outdir,"mean_fmin.pdf"),22cm,14cm),p)
#draw(PNG(joinpath(outdir,"mean_fmin.png"),dpi=150,16cm,10cm),p)
draw(SVG(joinpath(outdir,"mean_fmin.svg"),16cm,10cm),p)

#per dimension
for (k,D) in zip(1:length(dimensions), dimensions)

    idx = sortperm([exp.(mean(log.(mean_fmin[i,:,end,k].+eps()))) for i=1:length(optimizers)],rev=true)

    p = plot(
        [layer(
            x=collect(run_lengths)/D,y=exp.(mean(log.(mean_fmin[i,:,:,k].+eps()),dims=1)),Geom.line,
            Theme(default_color=cols[i],line_width=2pt)
        ) for i in 1:size(mean_succ,1)]...,
        layer(yintercept=[Δf],Geom.hline,Theme(default_color=colorant"gray")),
        Guide.title("All functions"), Guide.xlabel("Run Length / D"), Guide.ylabel("fmin"),
        Guide.manual_color_key("", opt_strings[idx], cols[idx]),
    #    Scale.x_log10,
        Scale.y_log10,
    )

    #draw(PDF(joinpath(outdir,"per_dimension","mean_fmin_$(D).pdf"),20cm,12cm),p)
    #draw(PNG(joinpath(outdir,"per_dimension","mean_fmin_$(D).png"),dpi=150,16cm,10cm),p)
    draw(SVG(joinpath(outdir,"per_dimension","mean_fmin_$(D).svg"),16cm,10cm),p)
end

## distance to minimizer

for (k,D) in zip(1:length(dimensions), dimensions)

    p = plot(
        [layer(x=collect(run_lengths)/D,y=median(mean_dist[i,:,:,k],dims=1),Geom.line,Theme(default_color=cols[i])) for i in 1:size(mean_succ,1)]...,
        Guide.title("All functions"), Guide.xlabel("Run Length / D"), Guide.ylabel("Distance to xmin"),
        Guide.manual_color_key("", opt_strings, cols),
    #    Scale.x_log10,
        Scale.y_log10,
    )

#    #draw(PDF(joinpath(outdir,"per_dimension","mean_dist_$(D).pdf"),20cm,12cm),p)
    draw(SVG(joinpath(outdir,"per_dimension","mean_dist_$(D).svg"),16cm,10cm),p)
end

## runtime

idx = sortperm(mean(runtime,dims=2)[:],rev=true)

p = plot(
    layer(
        y=opt_strings[idx], x=(mean(runtime,dims=2)/Base.minimum(mean(runtime,dims=2)))[idx],
        color = cols[idx],
        Geom.bar(orientation=:horizontal)
    ),
    Guide.title("All functions"), Guide.ylabel(""),
    Guide.xlabel("Relative Run Time",orientation=:horizontal),
    Scale.x_log10,
#    Scale.y_log10,
)

#draw(PDF(joinpath(outdir,"runtime.pdf"),20cm,12cm),p)
#draw(PNG(joinpath(outdir,"runtime.png"),dpi=150,16cm,10cm),p)
draw(SVG(joinpath(outdir,"runtime.svg"),16cm,10cm),p)
nothing

# convert svg to png (Gadfly's png export is broken at the moment)
# for f in *.svg; do
#     svgexport "$f" "${f%.svg}.png" 1.5x
# done

##


text:
module LowLevelParticleFilters

export KalmanFilter, UnscentedKalmanFilter, ParticleFilter, AuxiliaryParticleFilter, AdvancedParticleFilter, SigmaFilter, PFstate, index, state, covariance, num_particles, effective_particles, weights, expweights, particles, particletype, smooth, sample_measurement, simulate, loglik, log_likelihood_fun, forward_trajectory, mean_trajectory, mode_trajectory, update!, predict!, correct!, reset!, metropolis, shouldresample, TupleProduct

export densityplot, debugplot, commandplot, trajectorydensity, dimensiondensity

using StatsBase, Parameters, Lazy, Random, LinearAlgebra, Printf, LoopVectorization
import PDMats # To extend some methods on static arrays
using StaticArrays
using Distributions
using RecipesBase

abstract type ResamplingStrategy end
struct ResampleSystematic <: ResamplingStrategy end

abstract type AbstractFilter end

include("PFtypes.jl")
include("kalman.jl")
include("filtering.jl")
include("resample.jl")
include("utils.jl")
include("smoothing.jl")
include("plotting.jl")

index(f::AbstractFilter) = f.t[]

using Requires
function __init__()
    @require Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" include("plots.jl")
end

end # module


text:
const changeheight!(tree, height) = setheight!(tree, height)
const changetarget!(tree, branch, new_target) = settarget!(tree, branch, new_target)
const changesource!(tree, branch, new_source) = setsource!(tree, branch, new_source)
const changelength!(tree, branch, new_length) = setlength!(tree, branch, new_length)

"""
setheight!(tree::Tree,
           height::Float64)

Sets the height of a `Tree`
"""
function setheight!(tree::Tree,
                    height::Float64)
    tree.height = height
    return tree
end

"""
addnode!(tree::Tree)

Adds a single `Node` to a `Tree`
"""
function addnode!(tree::Tree)
  # Find highest node index
  if length(tree.nodes) == 0
    max_node = 0
  else
    max_node = maximum(keys(tree.nodes))
  end
  # Add the new node
  tree.nodes[max_node+1] = Node()
  # Return the updated tree
  return tree
end

"""
addnodes!(tree::Tree,
          nodes::Int64)

Adds multiple `Node`s to a `Tree`
"""
function addnodes!(tree::Tree,
                   nodes::Int64)
  # Find highest node index
  if length(tree.nodes) == 0
    max_node = 0
  else
    max_node = maximum(keys(tree.nodes))
  end
  # Add the new nodes
  for i = 1:nodes
    tree.nodes[max_node+i] = Node()
  end
  # Return the updated tree
  return tree
end

"""
deletenode!(tree::Tree,
            nodes::Int64)

Deletes a specified `Node` from a `Tree`
"""
function deletenode!(tree::Tree,
                     node::Int64)
  # Check if node exists
  if !haskey(tree.nodes, node)
    error("Node does not exist")
  # Check if any branches are connected to node
  elseif length(tree.nodes[node].in) > 0 || length(tree.nodes[node].out) > 0
    error("There are currently branches connected to this node")
  end
  # Remove node itself
  delete!(tree.nodes, node)
  return tree
end

"""
addbranch!(tree::Tree,
           source::Int64,
           target::Int64,
           branch_length::Float64)

Adds a `Branch` of specified length to a `Tree`
"""
function addbranch!(tree::Tree,
                    source::Int64,
                    target::Int64,
                    branch_length::Float64)
  # Error checking
  if target == source
    error("Branch must connect unique nodes")
  elseif !haskey(tree.nodes, source)
    error("Source node does not exist")
  elseif !haskey(tree.nodes, target)
    error("Target node does not exist")
  elseif length(tree.nodes[target].in) == 1
    error("The in degree of the target node is > 1")
  end
  # Find highest branch index
  if length(tree.branches) == 0
    max_branch = 0
  else
    max_branch = maximum(keys(tree.branches))
  end
  # Add the new branch
  tree.branches[max_branch+1] = Branch(source, target, branch_length)
  # Update the associated source and target nodes
  push!(tree.nodes[source].out, max_branch+1)
  push!(tree.nodes[target].in, max_branch+1)
  # Return updated tree
  return tree
end

"""
branch!(tree::Tree,
        source::Int64,
        branch_length::Float64)

Adds a `Branch` of a specified length, and a `Node` to a `Tree`
"""
function branch!(tree::Tree,
                 source::Int64,
                 branch_length::Float64)
  # Error checking
  if !haskey(tree.nodes, source)
    error("Source node does not exist")
  end
  # Find highest node index
  if length(tree.nodes) == 0
    max_node = 0
  else
    max_node = maximum(keys(tree.nodes))
  end
  # Add the new node
  tree.nodes[max_node+1] = Node()
  # Find highest branch index
  if length(tree.branches) == 0
    max_branch = 0
  else
    max_branch = maximum(keys(tree.branches))
  end  # Add the new branch
  tree.branches[max_branch+1] = Branch(source, max_node+1, branch_length)
  # Update the associated source and target nodes
  push!(tree.nodes[source].out, max_branch+1)
  push!(tree.nodes[max_node+1].in, max_branch+1)
  # Return updated tree
  return tree
end

"""
setsource!(tree::Tree,
              branch::Int64,
              new_source::Int64)

Set the source `Node` of a `Branch`
"""
function setsource!(tree::Tree,
                    branch::Int64,
                    new_source::Int64)
  # Error checking
  if !haskey(tree.branches, branch)
    error("Branch does not exist")
  elseif !haskey(tree.nodes, new_source)
    error("New source node does not exist")
  elseif new_source == tree.branches[branch].target
    error("Branch must connect unique nodes")
  end
  # Identify old source node
  old_source = tree.branches[branch].source
  # Remove branch reference from old source node
  splice!(tree.nodes[old_source].out, findfirst(tree.nodes[old_source].out .== branch))
  # Update new source node with branch reference
  push!(tree.nodes[new_source].out, branch)
  # Update branch with new source node reference
  tree.branches[branch].source = new_source
  # Return updated tree
  return tree
end

"""
settarget!(tree::Tree,
              branch::Int64,
              new_target::Int64)

Set the source `Node` of a `Branch`
"""
function settarget!(tree::Tree,
                    branch::Int64,
                    new_target::Int64)
  if !haskey(tree.branches, branch)
    error("Branch does not exist")
  elseif !haskey(tree.nodes, new_target)
    error("New target node does not exist")
  elseif length(tree.nodes[new_target].in) != 0
    error("New target node has an in degree > 1")
  end
  # Identify old target node
  old_target = tree.branches[branch].target
  # Remove branch reference from old target node
  splice!(tree.nodes[old_target].in, findfirst(tree.nodes[old_target].in .== branch))
  # Add branch reference to new target node
  push!(tree.nodes[new_target].in, branch)
  # Update branch with new target node reference
  tree.branches[branch].target = new_target
  # Return updated tree
  return tree
end

"""
setlength!(tree::Tree,
              branch::Int64,
              new_length::Float64)

Set the length of a `Branch`
"""
function setlength!(tree::Tree,
                    branch::Int64,
                    new_length::Float64)
    if !haskey(tree.branches, branch)
        error("Branch does not exist")
    elseif new_length < 0.0
        error("New branch length is < 0.0")
    end
    tree.branches[branch].length = new_length
end

"""
deletebranch!(tree::Tree,
              branch::Int64)

Delete a `Branch` from a `Tree`
"""
function deletebranch!(tree::Tree,
                       branch::Int64)
  if !haskey(tree.branches, branch)
    error("Branch does not exist")
  end
  # Identify old source node
  old_source = tree.branches[branch].source
  # Identify old target node
  old_target = tree.branches[branch].target
  # Remove branch reference from its source node
  splice!(tree.nodes[old_source].out, findfirst(tree.nodes[old_source].out .== branch))
  # Remove branch reference from its target node
  splice!(tree.nodes[old_target].in, findfirst(tree.nodes[old_target].in .== branch))
  # Remove branch itself
  delete!(tree.branches, branch)
  # Return the updated tree
  return tree
end


text:
"""
    struct ImpurityScale{T, N, S, AT} <: AbstractArray{T, N}

Impurity scalar field of the simulation.
It is kinda an alpha map for the impurity density of semiconductors.
It can takes values between `0` and `1`:
 * `1`: The impurity density has its full value. For grid points in depleted region of the semiconductor.
 * `]0,1[`: The impurity density is scaled down but not zero. For grid points at the edge of the depleted region which are partially depleted.
 * `0`: The impurity density is set to `0`. For grid points in undepleted regions of the semiconductor.

## Parametric types
* `T`: Element type of `data`.
* `N`: Dimension of the `grid` and `data` array.
* `S`: Coordinate system (`Cartesian` or `Cylindrical`).
* `AT`: Axes type.

## Fields
* `data::Array{T, N}`: Array containing the values of the impurity scale at the discrete points of the `grid`.
* `grid::Grid{T, N, S, AT}`: [`Grid`](@ref) defining the discrete points for which the impurity scale is determined.
"""
struct ImpurityScale{T, N, S, AT} <: AbstractArray{T, N}
    data::Array{T, N}
    grid::Grid{T, N, S, AT}
end

@inline size(impscale::ImpurityScale{T, N, S}) where {T, N, S} = size(impscale.data)
@inline length(impscale::ImpurityScale{T, N, S}) where {T, N, S} = length(impscale.data)
@inline getindex(impscale::ImpurityScale{T, N, S}, I::Vararg{Int, N}) where {T, N, S} = getindex(impscale.data, I...)
@inline getindex(impscale::ImpurityScale{T, N, S}, i::Int) where {T, N, S} = getindex(impscale.data, i)
@inline getindex(impscale::ImpurityScale{T, N, S}, s::Symbol) where {T, N, S} = getindex(impscale.grid, s)

function NamedTuple(impscale::ImpurityScale{T, 3}) where {T}
    return (
        grid = NamedTuple(impscale.grid),
        values = impscale.data,
    )
end
Base.convert(T::Type{NamedTuple}, x::ImpurityScale) = T(x)

function ImpurityScale(nt::NamedTuple)
    grid = Grid(nt.grid)
    T = typeof(ustrip(nt.values[1]))
    S = get_coordinate_system(grid)
    N = get_number_of_dimensions(grid)
    ImpurityScale{T, N, S, typeof(grid.axes)}(nt.values, grid)
end
Base.convert(T::Type{ImpurityScale}, x::NamedTuple) = T(x)

"""
    get_ticks_at_positions_of_edge_of_depleted_volumes(impscale::ImpurityScale{T, 3})

The impurity scale field is analyzed in order to find and return ticks where the gradient is strong,
which is the case at the surface of the depleted volume of a semiconductor.
"""
function get_ticks_at_positions_of_edge_of_depleted_volumes(impscale::ImpurityScale{T, 3}) where T
    dims = (1, 2, 3)
    mps::NTuple{3, Vector{T}} = broadcast(idim -> StatsBase.midpoints(impscale.grid[idim]), dims)
    extended = length.(mps) .> 0
    max_diffs = broadcast(idim -> extended[idim] ? map(i->maximum(abs.(selectdim(impscale.data, idim, i+1) .- selectdim(impscale.data, idim, i))), eachindex(mps[idim]))::Vector{T} : T[], dims);
    inds = broadcast(idim -> extended[idim] ? findall(Δ -> 0.5 < Δ < 1, max_diffs[idim]) : Int[], dims) # 0.5 seems to be a good limit as we don't want to add too many ticks.
    # As the impurity scale only has values between `0` and `1` the change (`max_diffs`) can also only between `0` and `1`.
    # We want to exclude a change of `1` as this happens at the surface of a semiconductor where the detector is undepleted
    # and we only want to generate ticks in the semiconductor at the position between undepleted and depleted regions.
    return broadcast(idim -> extended[idim] ? mps[idim][inds[idim]] : T[], dims)
end


text:
"""
    targets2bfr(redis, message, outdir, telinfo_file)

Subscriber worker function for Redis `target-selector:targets` channel.  This
function requires extra parameters beyond what the Redis `subscribe()` callback
supports.  It is expected that `subscribe()` will be passed a helper function
that will call this function with suitable parameters.

Messages published to the `target-selector:targets` channel have the format:

    targets:<telescope_name>:<subarray_name>:<timestamp>

This entire message is also the name of a Redis key that contains a JSON encoded
list of targets.  For each subscribed message, it performs these steps:

1. Parses targets `message` to extract the various components.
2. Gets a GuppiRaw.Header object by calling `subarray2grh()`
3. Retrieves the targets from Redis using key `targets:<telescope_name>`
3. Creates a BeamformerRecipe object
4. Outputs that object to an HDF5 file in `outdir` using the filename
   `<telescope_name>-<subarray_name>-<timestamp>.bfr5`.
"""
function targets2bfr(redis::Redis.RedisConnectionBase, message, outdir,
    telinfo_file::AbstractString=joinpath(ENV["HOME"], "telinfo.yml"))

    subarray_name = split(message*":::", ':')[3]
    grh = subarray2grh(redis, subarray_name)
    targets2bfr(grh, redis, message, outdir; telinfo_file=telinfo_file)
end

function targets2bfr(grh::GuppiRaw.Header, redis, message, outdir,
    telinfo_file::AbstractString=joinpath(ENV["HOME"], "telinfo.yml"))

    if isempty(grh)
        @warn "no GUPPI RAW metadata found for subarray $subarray_name" _module=nothing _file=nothing
    else
        _, telescope_name, subarray_name, timestamp = split(message, ':')
        beam_names, beam_positions = gettargets(redis, message)
        bfr = BeamformerRecipe(grh, beam_names, beam_positions; redis=redis, telinfo_file=telinfo_file)
        mkpath(outdir)
        h5name = joinpath(outdir, "$telescope_name-$subarray_name-$timestamp.bfr5")
        if ispath(h5name)
            @warn "refusing to overwrite existing path" _module=nothing _file=nothing
        else
            to_hdf5(h5name, bfr)
        end
    end
end

"""
Retrieves a JSON encoded targets list from `redis` using `key`, parses it, and
returns `beam_names` and `beam_positions` Arrays suitable for passing to
`BeamformerRecipe`.
"""
function gettargets(redis, key)
    targets_json = get(redis, key)
    targets = YAML.load(targets_json)
    beam_names = getindex.(targets, "source_id")
    beam_positions = permutedims(getindex.(targets, ["ra" "dec"]))
    beam_names, deg2rad.(beam_positions)
end

# This is a long running function.  It never returns until interrupted by
# Ctrl-C.  It is intended to be invoked from the command line by:
#
# julia --project=DIR -e 'using BluseBeamformerRecipes; run_targets_listener()' OUTDIR [TELINFO_FILE]
function run_targets_listener(args...)
    if isempty(args)
        args = ARGS
        if isempty(args)
            error("no directory given")
        end
    end
    outdir = args[1]
    telinfo_file = get(args, 2, joinpath(ENV["HOME"], "telinfo.yml"))

    redishost = get(ENV, "REDISHOST", "redishost")
    redis = Redis.RedisConnection(host=redishost)
    sub = open_subscription(redis)

    subscribe(sub, "target-selector:targets", message->begin
        @info "got message" message
        try
            targets2bfr(redis, message, outdir, telinfo_file)
        catch e
            @error e _module=nothing _file=nothing
        end
    end)
    @info "subscribed to target-selector:targets"

    Base.exit_on_sigint(false)
    try
        while(true)
            sleep(1)
        end
    catch e
        e isa InterruptException || rethrow()
    finally
        disconnect(sub)
        disconnect(redis)
    end

    @info "done"
end


text:
import Pkg

pkgs = ["NNlib"]

Pkg.test(pkgs; coverage = true)


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 = "librdkafka_build"
version = v"0.1.0"

# Collection of sources required to build librdkafka_build
sources = [
    "https://github.com/edenhill/librdkafka/archive/v1.0.0.tar.gz" =>
    "b00a0d9f0e8c7ceb67b93b4ee67f3c68279a843a15bf4a6742eb64897519aa09",

]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd librdkafka-1.0.0/
./configure
make
cp src/librdkafka*.so $prefix/ || true
cp src/librdkafka*.so.1 $prefix/ || true
cp src/librdkafka*.dll $prefix/ || true
cp src/librdkafka*.dylib $prefix || true
exit

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:i686, libc=:glibc),
    Linux(:x86_64, libc=:glibc),
    Linux(:aarch64, libc=:glibc),
    Linux(:powerpc64le, libc=:glibc),
    Linux(:i686, libc=:musl),
    Linux(:x86_64, libc=:musl),
    Linux(:aarch64, libc=:musl),
    Linux(:armv7l, libc=:musl, call_abi=:eabihf),
    # MacOS(:x86_64),
    FreeBSD(:x86_64)
]

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

# Dependencies that must be installed before this package can be built
dependencies = [

]

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


text:
function sanitize_component!(line::Union{MonitoredLine, Line}, sys::System)
    sanitize_angle_limits!(line)
end

function validate_component_with_system(line::Union{MonitoredLine, Line}, sys::System)
    is_valid = true
    if !check_endpoint_voltages(line)
        is_valid = false
    elseif !calculate_thermal_limits!(line, get_base_power(sys))
        is_valid = false
    end
    return is_valid
end

function sanitize_angle_limits!(line)
    max_limit = pi / 2
    min_limit = -pi / 2

    orderedlimits(line.angle_limits, "Angles")

    if (line.angle_limits.max / max_limit > 3) ||
       (-1 * line.angle_limits.min / max_limit > 3)
        @warn "The angle limits provided is larger than 3π/2 radians.\n " *
              "PowerSystems inferred the data provided in degrees and will transform it to radians" maxlog =
            PS_MAX_LOG

        if line.angle_limits.max / max_limit >= 0.99
            line.angle_limits = (
                min = line.angle_limits.min,
                max = min(line.angle_limits.max * (π / 180), max_limit),
            )
        else
            line.angle_limits =
                (min = line.angle_limits.min, max = min(line.angle_limits.max, max_limit))
        end

        if (-1 * line.angle_limits.min / max_limit > 0.99)
            line.angle_limits = (
                min = max(line.angle_limits.min * (π / 180), min_limit),
                max = line.angle_limits.max,
            )
        else
            line.angle_limits =
                (min = max(line.angle_limits.min, min_limit), max = line.angle_limits.max)
        end
    else
        if line.angle_limits.max >= max_limit && line.angle_limits.min <= min_limit
            line.angle_limits = (min = min_limit, max = max_limit)
        elseif line.angle_limits.max >= max_limit && line.angle_limits.min >= min_limit
            line.angle_limits = (min = line.angle_limits.min, max = max_limit)
        elseif line.angle_limits.max <= max_limit && line.angle_limits.min <= min_limit
            line.angle_limits = (min = min_limit, max = line.angle_limits.max)
        elseif line.angle_limits.max == 0.0 && line.angle_limits.min == 0.0
            line.angle_limits = (min = min_limit, max = max_limit)
        end
    end
end

function linerate_calculation(l::Union{Line, MonitoredLine})
    theta_max = max(abs(l.angle_limits.min), abs(l.angle_limits.max))
    g = l.r / (l.r^2 + l.x^2)
    b = -l.x / (l.r^2 + l.x^2)
    y_mag = sqrt(g^2 + b^2)
    fr_vmax = isnothing(l.arc.from.voltage_limits) ? nothing : l.arc.from.voltage_limits.max
    to_vmax = isnothing(l.arc.to.voltage_limits) ? nothing : l.arc.to.voltage_limits.max

    if isa(fr_vmax, Nothing) || isa(to_vmax, Nothing)
        fr_vmax = 1.0
        to_vmax = 0.9
        diff_angle = abs(l.arc.from.angle - l.arc.to.angle)
        new_rate = y_mag * fr_vmax * to_vmax * cos(theta_max)
    else
        m_vmax = max(fr_vmax, to_vmax)
        c_max = sqrt(fr_vmax^2 + to_vmax^2 - 2 * fr_vmax * to_vmax * cos(theta_max))
        new_rate = y_mag * m_vmax * c_max
    end

    return new_rate
end

function calculate_thermal_limits!(branch, basemva::Float64)
    is_valid = true
    if get_rate(branch) < 0.0
        @error "PowerSystems does not support negative line rates"
        is_valid = false

    elseif get_rate(branch) == 0.0
        @warn "Data for line rating is not provided, PowerSystems will infer a rate from line parameters" maxlog =
            PS_MAX_LOG
        if get_angle_limits(branch) == get_angle_limits(Line(nothing))
            branch.rate =
                min(calculate_sil(branch, basemva), linerate_calculation(branch)) / basemva
        else
            branch.rate = linerate_calculation(branch) / basemva
        end

    elseif get_rate(branch) > linerate_calculation(branch)
        mult = get_rate(branch) / linerate_calculation(branch)
        if mult > 50
            @warn "Data for line rating is $(mult) times larger than the base MVA for the system" maxlog =
                PS_MAX_LOG
        end
    end

    check_SIL(branch, basemva)

    return is_valid
end

const SIL_STANDARDS = Dict( #from https://neos-guide.org/sites/default/files/line_flow_approximation.pdf
    69.0 => (min = 12.0, max = 13.0),
    138.0 => (min = 47.0, max = 52.0),
    230.0 => (min = 134.0, max = 145.0),
    345.0 => (min = 325.0, max = 425.0),
    500.0 => (min = 850.0, max = 1075.0),
    765.0 => (min = 2200.0, max = 2300.0),
)

# calculation from https://neos-guide.org/sites/default/files/line_flow_approximation.pdf
function calculate_sil(line, basemva::Float64)
    arc = get_arc(line)
    #Assumess voltage at both ends of the arc is the same
    vrated = get_base_voltage(get_to(arc))
    zbase = vrated^2 / basemva
    l = get_x(line) / (2 * pi * 60) * zbase
    r = get_r(line) * zbase
    c = sum(get_b(line)) / (2 * pi * 60 * zbase)
    # A line with no C doesn't allow calculation of SIL
    isapprox(c, 0.0) && return Inf
    zc = sqrt((r + im * 2 * pi * 60 * l) / (im * 2 * pi * 60 * c))
    sil = vrated^2 / abs(zc)
    return sil
end

function check_SIL(line, basemva::Float64)
    arc = get_arc(line)
    vrated = get_base_voltage(get_to(arc))
    SIL_levels = collect(keys(SIL_STANDARDS))
    rate = get_rate(line)
    closestV = findmin(abs.(SIL_levels .- vrated))
    closestSIL = SIL_STANDARDS[SIL_levels[closestV[2]]]

    # Assuming that the rate is in pu
    if (rate >= closestSIL.max / basemva)
        # rate outside of expected SIL range
        sil = calculate_sil(line, basemva)
        @warn "Rate $(round(rate*basemva; digits=2)) MW for $(line.name) is larger than the max expected in the range of $(closestSIL)." maxlog =
            PS_MAX_LOG
    end
    return
end

function check_endpoint_voltages(line)
    is_valid = true
    arc = get_arc(line)
    if get_base_voltage(get_from(arc)) != get_base_voltage(get_to(arc))
        is_valid = false
        @error "Voltage endpoints of $(line.name) are different, cannot create Line"
    end
    return is_valid
end


text:
function to_non_normal(l::Vector{FieldsTower}, G::GAP.GapObj, deg::Int)
  for x in l
    assure_automorphisms(x)
    assure_isomorphism(x, G)
  end
  lC = GAP.Globals.ConjugacyClassesSubgroups(G)
  ind = 0
  for i = 1:length(lC)
    r = GAP.Globals.Representative(lC[i])
    if GAP.Globals.Size(r) == divexact(degree(l[1].field), deg) && GAP.Globals.Size(GAP.Globals.Core(G, r)) == 1
      ind = i
      break
    end
  end
  if iszero(ind)
    error("Representation not possible")
  end
  rep = GAP.Globals.Representative(lC[ind])
  ffields = Vector{AnticNumberField}(undef, length(l))
  for i = 1:length(ffields)
    ffields[i] = fixed_field(l[i], rep)
  end
  return ffields
end

function fixed_field(x::FieldsTower, H::GAP.GapObj)
  gH = GAP.Globals.SmallGeneratingSet(H)

  auts = NfToNfMor[]
  found = 0
  D = x.isomorphism
  autsx = automorphisms(number_field(x))
  i = 0
  while length(gH) != found
    i += 1
    if D[autsx[i]] in gH
      push!(auts, autsx[i])
      found += 1
    end
  end
  return fixed_field(number_field(x), auts)[1]
end

function _find_discriminant_bound(n, i, disc)
  Gt = GAP.Globals.TransitiveGroup(n, i)
  @assert GAP.Globals.IsSolvable(Gt)
  id = GAP.Globals.IdGroup(Gt)
  G1 = GAP.Globals.SmallGroup(id)
  lC = GAP.Globals.ConjugacyClassesSubgroups(G1)
  ind = 0
  for k = 1:length(lC)
    H = GAP.Globals.Representative(lC[k])
    if GAP.Globals.Index(G1, H) != n
      continue
    end
    core = GAP.Globals.Core(G1, H)
    if !isone(GAP.Globals.Size(core))
      continue
    end
    mp = GAP.Globals.FactorCosetAction(G1, H)
    idT = GAP.Globals.TransitiveIdentification(GAP.Globals.Image(mp))
    if idT == i
      ind = i
      break
    end
  end
  cg = lC[ind]
  H = GAP.Globals.Representative(cg)
  conjs = GAP.Globals.Elements(cg)
  #I check if it is a Frobenius group!
  isfrobenius = true
  for i = 1:length(conjs)
    for j = i+1:length(conjs)
      if GAP.Globals.Size(GAP.Globals.Intersection(conjs[i], conjs[j])) != 1
        isfrobenius = false
        break
      end
    end
    if !isfrobenius
      break
    end
  end
  if isfrobenius
    @show "Frobenius!"
    #In this case, we can find a better bound for the closure!
    m = divexact(id[1], n)
    bdisc = disc^(2*m*n-m)
    return root(bdisc, n-1)
  end
  j = 1
  for i = 2:length(conjs)
    H = GAP.Globals.Intersection(H, conjs[i])
    j += 1
    if GAP.Globals.Size(H) == 1
      break
    end
  end
  return disc^(j*divexact(id[1], n))
end

function fields_transitive_group(n::Int, i::Int, disc::fmpz)
  Gt = GAP.Globals.TransitiveGroup(n, i)
  @assert GAP.Globals.IsSolvable(Gt)
  id = GAP.Globals.IdGroup(Gt)
  G1 = GAP.Globals.SmallGroup(id)
  @show disc_NC = _find_discriminant_bound(n, i, disc)
  lf = fields(id[1], id[2], disc_NC)
  ln = to_non_normal(lf, G1, n)
  indices = Int[]
  for i = 1:length(ln)
    if abs(discriminant(maximal_order(ln[i]))) <= disc
      push!(indices, i)
    end
  end
  return ln[indices]
end


text:
module OceanBasins

using DataDeps
using DelimitedFiles
using StaticArrays
using PolygonOps

function readdata(datadepname)
    registerfile(datadepname)
    datadeppath = @datadep_str datadepname
    filepath = joinpath(datadeppath, "Limits_of_oceans_and_seas.tab")
    return readdlm(filepath, '\t'; header=true, skipstart=15)
end

struct OceanOrSea{T}
    name::String
    polygon::T
end
name(o::OceanOrSea) = o.name
polygon(o::OceanOrSea) = o.polygon

# lat, lon because that's how the polygons are made
struct Point2D{T} <: FieldVector{2,T}
    lat::T
    lon::T
end

function oceanpolygons(datadepname="Oceans_and_seas")
    cells, headers = readdata(datadepname)
    names = unique(string.(cells[:,1]))
    out = OceanOrSea[]
    for name in names
        iocn = findall(isequal(name), cells[:,1])
        poly = [Point2D(cells[j,2], fixmeridian(cells[j,3])) for j in iocn]
        push!(out, OceanOrSea(name, poly))
    end
    return out
end

# Shortlist of oceans/seas of interest for me
# Don't hesitate to add your own in PRs
const SHORTLIST = [60,61,85,120,146,147,100]
west_pacific() = 60
east_pacific() = 61
atlantic() = 85
indian() = 120
west_arctic() = 146
east_arctic() = 147
mediterranean() = 100
# that have their functions built-in
ispacific(P, oceans) = P ∈ oceans[west_pacific()] || P ∈ oceans[east_pacific()]
isatlantic(P, oceans) = P ∈ oceans[atlantic()]
isindian(P, oceans) = P ∈ oceans[indian()]
isarctic(P, oceans) = P ∈ oceans[west_arctic()] || P ∈ oceans[east_arctic()] || P.lat > 70
ismediterranean(P, oceans) = P ∈ oceans[mediterranean()]
isantarctic(P, oceans) = P.lat ≤ -40
isatlantic2(P, oceans) = isatlantic(P, oceans) && !isantarctic(P, oceans)
ispacific2(P, oceans) = ispacific(P, oceans) && !isantarctic(P, oceans)
isindian2(P, oceans) = isindian(P, oceans) && !isantarctic(P, oceans)
convert_65S_to_90S(P::Point2D{T}) where T = Point2D(P.lat < -65 ? T(-65) : P.lat, P.lon)
isSOatlantic(P, oceans) = isatlantic(convert_65S_to_90S(P), oceans) && isantarctic(P, oceans)
isSOpacific(P, oceans) = ispacific(convert_65S_to_90S(P), oceans) && isantarctic(P, oceans)
isSOindian(P, oceans) = isindian(convert_65S_to_90S(P), oceans) && isantarctic(P, oceans)
isindonesian(P, oceans) = foldl(|, P ∈ oceans[i] for i in 11:37) # Indo throughflow (TODO fix name?)
ismediterranean2(P, oceans) = ismediterranean(P, oceans) || ((30 ≤ P.lat ≤ 50) && !isatlantic(P, oceans) && !ispacific(P, oceans))
isNPAC(P, oceans) = ispacific(P, oceans) && P.lat > 0
isNATL(P, oceans) = isatlantic(P, oceans) && P.lat > 0
isSPAC(P, oceans) = ispacific(P, oceans) && !isantarctic(P, oceans) && P.lat ≤ 0
isSATL(P, oceans) = isatlantic(P, oceans) && !isantarctic(P, oceans) && P.lat ≤ 0

for ocn in (:pacific, :atlantic, :indian, :arctic, :mediterranean, :antarctic,
            :pacific2, :atlantic2, :indian2, :indonesian, :SOatlantic, :SOpacific,
            :SOindian, :mediterranean2, :NPAC, :SPAC, :NATL, :SATL)
    f = Symbol(:is, ocn)
    @eval begin
        """
            $($f)(NT::NamedTuple, oceans)

        Returns `$($f)(NT.lat, NT.lon, oceans)`.
        """
        $f(NT::NamedTuple, oceans) = $f(NT.lat, NT.lon, oceans)
        """
            $($f)(Ps::Vector, oceans)

        Returns the bit vector `[$($f)(P, oceans) for P in Ps]`.
        """
        $f(Ps::Vector, oceans) = [$f(P, oceans) for P in Ps]
        """
            $($f)(lat, lon, oceans)

        Returns `true` if the `(lat,lon)` coordinate is the polygon.
        """
        $f(lat, lon, oceans) = $f(Point2D(lat,lon), oceans)
        """
            $($f)(lat::Vector, lon::Vector, oceans)

        Returns `$($f)(Point2D.(lat,lon), oceans)`.
        """
        $f(lat::Vector, lon::Vector, oceans) = $f(Point2D.(lat,lon), oceans)
        export $f
    end
end

# Convert ±179.99899 to ±180
function fixmeridian(lon)
    if abs(lon) == 179.99899
        return 181.0sign(lon)
    elseif abs(lon) == 0.001
        return 0.0
    else
        return lon
    end
end

whichoceans(P, oceans) = [name(ocn) for ocn in oceans if P ∈ ocn]
whichoceans(lat, lon, oceans) = whichoceans(Point2D(lat,lon), oceans)

whichshortlistoceans(P, oceans) = [name(ocn) for ocn in oceans[SHORTLIST] if P ∈ ocn]
whichshortlistoceans(lat, lon, oceans) = whichshortlistoceans(Point2D(lat,lon), oceans)

function whichshortlistocean(lat, lon, oceans)
    list = whichshortlistoceans(lat, lon, oceans)
    if isempty(list)
        return "None"
    elseif (length(list) == 1)
        return shorten_name(list[1])
    else
        return shorten_name.(list)
    end
end

export OceanOrSea, oceanpolygons, whichshortlistoceans, whichshortlistocean, whichoceans

Base.in(P, ocn::OceanOrSea) = inpolygon(center(P), polygon(ocn), in=true, on=true, out=false)
center(P) = [P[1], mod(P[2] + 180, 360) - 180]

function shorten_name(name)
    name = replace(name, " Ocean" => "")
    name = replace(name, "South " => "S")
    name = replace(name, "North " => "N")
    name = replace(name, ", eastern part" => "e")
    name = replace(name, ", western part" => "w")
    name = replace(name, "Atlantic" => "ATL")
    name = replace(name, "Pacific" => "PAC")
    name = replace(name, "Indian" => "IND")
    name = replace(name, "Arctic" => "ARC")
    return name
end

function registerfile(datadepname)
    register(DataDep(
        datadepname,
        """
        Dataset: Limits of oceans and seas in digitized, machine readable form
        Website: https://doi.pangaea.de/10.1594/PANGAEA.777975
        Author: IHO International Hydrographic Organization, Rainer Sieger
        Date of Publication: March 20, 2012
        License: https://creativecommons.org/licenses/by/3.0/legalcode

        This dataset is the digitized version of the printed "Limits of Oceans and Seas" (IHO, 1953).
        The report describes the boundaries of 148 oceans and seas. The given positions were typed into a spreadsheet and were completed to a right ordered polygon by hand using Google Earth. The dataset consists of five columns.
        The first column contains the name of the ocean or sea, column 2 and 3 contain the position of each polygon point (latitude, longitude). Column 4 contains an area index, column 5 a polygon index. These two indices structurize the dataset.
        The other version link provides a Google Earth layer in KML format (zipped).
        This layer consists of all polygons and names of oceans or seas at area centroid position given in dataset doi:10.1594/PANGAEA.777976.

        Please cite this dataset: International Hydrographic Organization, I. H. O., & Sieger, R. (2012). <i>Limits of oceans and seas in digitized, machine readable form</i> [Data set]. PANGAEA - Data Publisher for Earth and Environmental Science. https://doi.org/10.1594/PANGAEA.777975
        """,
        ["https://doi.pangaea.de/10.1594/PANGAEA.777975?format=textfile"],
        "29c1ec5a1df3e575ef53f0a29de8e3351189666e63e8c2827265aed2025a204b"
    ))
end

end # module


text:
module LoziMap
export lozi_map

using Parameters: @with_kw, @unpack
using ..ExampleBase: LEDemo, DiscreteExample

@with_kw struct LoziMapParam
    a::Float64 = 1.4
    b::Float64 = 0.3
end

@inline function phase_dynamics!(u_next, u, p, t)
    @unpack a, b = p
    u_next[1] = 1 + u[2] - a * abs(u[1])
    u_next[2] = b * u[1]
end

@inline @views function tangent_dynamics!(u_next, u, p, t)
    @unpack a, b = p
    phase_dynamics!(u_next[:, 1], u[:, 1], p, t)
    u_next[1, 2:end] .= u[2, 2:end] .- a * sign(u[1, 1]) .* u[1, 2:end]
    u_next[2, 2:end] .= b .* u[1, 2:end]
end

"""
Return a [`LEDemo`](@ref) for the Lozi map.
"""
function lozi_map(;
        u0=[0.1, 0.1],
        t_renorm=10,
        t_attr=100000,
        atol=0, rtol=0.2,
        # terminator_options = [:atol => atol, :rtol => rtol],  # TODO: use it
        terminator_options = [],
        kwargs...)
    LEDemo(DiscreteExample(
        "Lozi map",
        phase_dynamics!, u0, t_renorm, LoziMapParam(),
        tangent_dynamics!,
        t_attr,
        Float64[],   # known_exponents
        atol, rtol,
        terminator_options,
    ); kwargs...)
end

end