text:
# FiniteRange gives the nonzero entries in a row/column
struct FiniteRange end

getindex(A::AbstractMatrix,::Type{FiniteRange},j::Integer) = A[1:colstop(A,j),j]
getindex(A::AbstractMatrix,k::Integer,::Type{FiniteRange}) = A[k,1:rowstop(A,k)]

const ⤓ = FiniteRange

mutable struct RaggedMatrix{T} <: AbstractMatrix{T}
    data::Vector{T} # a Vector of non-zero entries
    cols::Vector{Int} # a Vector specifying the first index of each column
    m::Int #Number of rows
    function RaggedMatrix{T}(data::Vector{T}, cols::Vector{Int}, m::Int) where T
        # make sure the cols are monitonically increasing
        @assert 1==cols[1]
        for j=1:length(cols)-1
            @assert cols[j] ≤ cols[j+1]
        end
        @assert cols[end] == length(data)+1

        # make sure we have less entries than the size of the matrix
        @assert length(data) ≤ m*(length(cols)-1)

        new{T}(data,cols,m)
    end
end

RaggedMatrix(dat::Vector,cols::Vector{Int},m::Int) =
    RaggedMatrix{eltype(dat)}(dat,cols,m)

RaggedMatrix{T}(::UndefInitializer, m::Int, colns::AbstractVector{Int}) where {T} =
    RaggedMatrix(Vector{T}(undef, sum(colns)),Int[1;1 .+ cumsum(colns)],m)

Base.size(A::RaggedMatrix) = (A.m,length(A.cols)-1)

colstart(A::RaggedMatrix,j::Integer) = 1
colstop(A::RaggedMatrix,j::Integer) = min(A.cols[j+1]-A.cols[j],size(A,1))

Base.IndexStyle(::Type{RM}) where {RM<:RaggedMatrix} = IndexCartesian()

function getindex(A::RaggedMatrix,k::Int,j::Int)
    if k>size(A,1) || k < 1 || j>size(A,2) || j < 1
        throw(BoundsError(A,(k,j)))
    end

    if A.cols[j]+k-1 < A.cols[j+1]
        A.data[A.cols[j]+k-1]
    else
        zero(eltype(A))
    end
end

function Base.setindex!(A::RaggedMatrix,v,k::Int,j::Int)
    if k>size(A,1) || k < 1 || j>size(A,2) || j < 1
        throw(BoundsError(A,(k,j)))
    end

    if A.cols[j]+k-1 < A.cols[j+1]
        A.data[A.cols[j]+k-1]=v
    elseif v ≠ 0
        throw(BoundsError(A,(k,j)))
    end
    v
end

convert(::Type{RaggedMatrix{T}}, R::RaggedMatrix{T}) where T = R

convert(::Type{RaggedMatrix{T}}, R::RaggedMatrix) where T =
    RaggedMatrix{T}(Vector{T}(R.data), R.cols, R.m)

function convert(::Type{Matrix{T}}, A::RaggedMatrix) where T
    ret = zeros(T,size(A,1),size(A,2))
    for j=1:size(A,2)
        ret[1:colstop(A,j),j] = view(A,1:colstop(A,j),j)
    end
    ret
end

convert(::Type{Matrix}, A::RaggedMatrix) = Matrix{eltype(A)}(A)

function convert(::Type{RaggedMatrix{T}}, B::BandedMatrix) where T
    l = bandwidth(B,1)
    ret = RaggedMatrix(Zeros{T}(size(B)), Int[colstop(B,j) for j=1:size(B,2)])
    for j=1:size(B,2),k=colrange(B,j)
        ret[k,j] = B[k,j]
    end
    ret
end

convert(::Type{RaggedMatrix}, B::BandedMatrix) = RaggedMatrix{eltype(B)}(B)

function convert(::Type{RaggedMatrix{T}}, B::AbstractMatrix) where T
    ret = RaggedMatrix(Zeros{T}(size(B)), Int[colstop(B,j) for j=1:size(B,2)])
    for j=1:size(B,2), k=colrange(B,j)
        ret[k,j] = B[k,j]
    end
    ret
end

convert(::Type{RaggedMatrix}, B::AbstractMatrix) = RaggedMatrix{eltype(B)}(B)

RaggedMatrix(B::AbstractMatrix) = convert(RaggedMatrix, B)
RaggedMatrix{T}(B::AbstractMatrix) where T = convert(RaggedMatrix{T}, B)

Base.similar(B::RaggedMatrix,::Type{T}) where {T} = RaggedMatrix(Vector{T}(length(B.data)),copy(B.cols),B.m)

for (op,bop) in ((:(Base.rand), :rrand),)
    @eval begin
        $bop(::Type{T}, m::Int, colns::AbstractVector{Int}) where {T} =
            RaggedMatrix($op(T,sum(colns)),[1; (1 .+ cumsum(colns))],m)
        $bop(m::Int, colns::AbstractVector{Int}) = $bop(Float64, m, colns)
    end
end

function RaggedMatrix{T}(Z::Zeros, colns::AbstractVector{Int}) where {T}
    if size(Z,2) ≠ length(colns)
        throw(DimensionMismatch())
    end
    RaggedMatrix(zeros(T,sum(colns)), [1; (1 .+cumsum(colns))], size(Z,1))
end

function RaggedMatrix{T}(A::AbstractMatrix, colns::AbstractVector{Int}) where T
    ret = RaggedMatrix{T}(undef, size(A,1), colns)
    @inbounds for j = 1:length(colns), k = 1:colns[j]
        ret[k,j] = A[k,j]
    end
    ret
end

RaggedMatrix(A::AbstractMatrix, colns::AbstractVector{Int}) = RaggedMatrix{eltype(A)}(A, colns)

## BLAS

function mul!(y::Vector, A::RaggedMatrix, b::Vector)
    m=size(A,2)

    if m ≠ length(b) || size(A,1) ≠ length(y)
        throw(BoundsError())
    end
    T=eltype(y)
    fill!(y,zero(T))
    for j=1:m
        kr=A.cols[j]:A.cols[j+1]-1
        BLAS.axpy!(b[j],view(A.data,kr),view(y,1:length(kr)))
    end
    y
end

function BLAS.axpy!(a, X::RaggedMatrix, Y::RaggedMatrix)
    if size(X) ≠ size(Y)
        throw(BoundsError())
    end

    if X.cols == Y.cols
        BLAS.axpy!(a,X.data,Y.data)
    else
        for j = 1:size(X,2)
            Xn = colstop(X,j)
            Yn = colstop(Y,j)
            if Xn > Yn  # check zeros otherwise
                for k = Yn+1:Xn
                    @assert iszero(X[k,j])
                end
            end
            cs = min(Xn,Yn)
            BLAS.axpy!(a,view(X.data,X.cols[j]:X.cols[j]+cs-1),
                         view(Y.data,Y.cols[j]:Y.cols[j]+cs-1))
        end
    end
    Y
end

colstop(X::SubArray{T,2,RaggedMatrix{T},Tuple{UnitRange{Int},UnitRange{Int}}},
     j::Integer) where {T} = min(colstop(parent(X),j + first(parentindices(X)[2])-1) -
                                            first(parentindices(X)[1]) + 1,
                            size(X,1))

function BLAS.axpy!(a,X::RaggedMatrix,
                    Y::SubArray{T,2,RaggedMatrix{T},Tuple{UnitRange{Int},UnitRange{Int}}}) where T
    if size(X) ≠ size(Y)
        throw(BoundsError())
    end

    P = parent(Y)
    ksh = first(parentindices(Y)[1]) - 1  # how much to shift
    jsh = first(parentindices(Y)[2]) - 1  # how much to shift

    for j=1:size(X,2)
        cx=colstop(X,j)
        cy=colstop(Y,j)
        if cx > cy
            for k=cy+1:cx
                if X[k,j] ≠ 0
                    throw(BoundsError("Trying to add a non-zero to a zero."))
                end
            end
            kr = X.cols[j]:X.cols[j]+cy-1
        else
            kr = X.cols[j]:X.cols[j+1]-1
        end

        BLAS.axpy!(a,view(X.data,kr),
                    view(P.data,(P.cols[j + jsh] + ksh-1) .+ (1:length(kr))))
    end

    Y
end

function *(A::RaggedMatrix,B::RaggedMatrix)
    cols = zeros(Int,size(B,2))
    T = promote_type(eltype(A),eltype(B))
    for j=1:size(B,2),k=1:colstop(B,j)
        cols[j] = max(cols[j],colstop(A,k))
    end

    unsafe_mul!(RaggedMatrix{T}(undef, size(A,1), cols), A, B)
end

function unsafe_mul!(Y::RaggedMatrix,A::RaggedMatrix,B::RaggedMatrix)
    fill!(Y.data,0)

    for j=1:size(B,2),k=1:colstop(B,j)
        BLAS.axpy!(B[k,j],view(A,1:colstop(A,k),k),view(Y.data,Y.cols[j] .- 1 .+ (1:colstop(A,k))))
    end

    Y
end

function mul!(Y::RaggedMatrix,A::RaggedMatrix,B::RaggedMatrix)
    for j=1:size(B,2)
        col = 0
        for k=1:colstop(B,j)
            col = max(col,colstop(A,k))
        end

        if col > colstop(Y,j)
            throw(BoundsError())
        end
    end

    unsafe_mul!(Y,A,B)
end


text:
function part1(input, number_of_moves = "100")
    order = parse.(Int, split(readchomp(input), ""))
    state, current = order_to_state(order)
    for _ in 1:parse(Int, number_of_moves)
        current = move!(state, current)
    end
    order = state_to_order(state, 1)
    return string(order[2:end]...)
end

function part2(input, number_of_moves = "10000000")
    order = parse.(Int, split(readchomp(input), ""))
    append!(order, (length(order) + 1):1_000_000)
    @assert length(order) == 1_000_000
    state, current = order_to_state(order)
    for _ in 1:parse(Int, number_of_moves)
        current = move!(state, current)
    end
    order = state_to_order(state, 1)
    return order[2] * order[3]
end

function order_to_state(order)
    state = zeros(Int, length(order))
    for i = 2:length(order)
        state[order[i - 1]] = order[i]
    end
    state[order[end]] = order[1]
    return state, order[1]
end

function state_to_order(state, first)
    order = zeros(Int, length(state))
    i = first
    for j = 1:length(state)
        order[j] = i
        i = state[i]
    end
    return order
end

function move!(state, current)
    removed1 = state[current]
    removed2 = state[removed1]
    removed3 = state[removed2]
    next_current = state[removed3]
    state[current] = next_current
    destination = current - 1
    while true
        if destination == 0
            destination = length(state)
        end
        if destination == removed1 || destination == removed2 || destination == removed3
            destination -= 1
            continue
        end
        break
    end
    state[removed3] = state[destination]
    state[destination] = removed1
    return next_current
end

ext: jl


text:
module EarthWrapper

using JavaCall

const AbsoluteDate = @jimport org.orekit.time.AbsoluteDate
const BodyShape = @jimport org.orekit.bodies.BodyShape
const Class = @jimport java.lang.Class
const Collection = @jimport java.util.Collection
const DataProvidersManager = @jimport org.orekit.data.DataProvidersManager
const EarthITU453AtmosphereRefraction = @jimport org.orekit.models.earth.EarthITU453AtmosphereRefraction
const EarthShape = @jimport org.orekit.models.earth.EarthShape
const EarthStandardAtmosphereRefraction = @jimport org.orekit.models.earth.EarthStandardAtmosphereRefraction
const Ellipse = @jimport org.orekit.bodies.Ellipse
const Ellipsoid = @jimport org.orekit.bodies.Ellipsoid
const FieldAbsoluteDate = @jimport org.orekit.time.FieldAbsoluteDate
const FieldGeodeticPoint = @jimport org.orekit.bodies.FieldGeodeticPoint
const FieldLine = @jimport org.hipparchus.geometry.euclidean.threed.FieldLine
const FieldVector3D = @jimport org.hipparchus.geometry.euclidean.threed.FieldVector3D
const Frame = @jimport org.orekit.frames.Frame
const GeoMagneticElements = @jimport org.orekit.models.earth.GeoMagneticElements
const GeoMagneticField = @jimport org.orekit.models.earth.GeoMagneticField
const GeoMagneticFieldFactory = @jimport org.orekit.models.earth.GeoMagneticFieldFactory
const GeoMagneticFieldFactory_FieldModel = @jimport org.orekit.models.earth.GeoMagneticFieldFactory$FieldModel
const GeoMagneticFields = @jimport org.orekit.models.earth.GeoMagneticFields
const GeoMagneticModelLoader = @jimport org.orekit.models.earth.GeoMagneticModelLoader
const GeodeticPoint = @jimport org.orekit.bodies.GeodeticPoint
const Geoid = @jimport org.orekit.models.earth.Geoid
const InputStream = @jimport java.io.InputStream
const JString = @jimport java.lang.String
const LazyLoadedGeoMagneticFields = @jimport org.orekit.models.earth.LazyLoadedGeoMagneticFields
const Line = @jimport org.hipparchus.geometry.euclidean.threed.Line
const NormalizedSphericalHarmonicsProvider = @jimport org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider
const Object = @jimport java.lang.Object
const OneAxisEllipsoid = @jimport org.orekit.bodies.OneAxisEllipsoid
const PVCoordinates = @jimport org.orekit.utils.PVCoordinates
const ReferenceEllipsoid = @jimport org.orekit.models.earth.ReferenceEllipsoid
const TimeStampedPVCoordinates = @jimport org.orekit.utils.TimeStampedPVCoordinates
const Vector3D = @jimport org.hipparchus.geometry.euclidean.threed.Vector3D

export EarthITU453AtmosphereRefraction
export EarthShape
export EarthStandardAtmosphereRefraction
export GeoMagneticElements
export GeoMagneticField
export GeoMagneticFieldFactory
export GeoMagneticFields
export GeoMagneticModelLoader
export Geoid
export LazyLoadedGeoMagneticFields
export ReferenceEllipsoid
export calculate_field
export equals
export get_a
export get_b
export get_body_frame
export get_c
export get_c2n0
export get_cartesian_intersection_point
export get_class
export get_decimal_year
export get_declination
export get_ellipsoid
export get_epoch
export get_equatorial_radius
export get_field
export get_field_vector
export get_flattening
export get_frame
export get_geo_magnetic_fields
export get_gm
export get_grs80
export get_horizontal_intensity
export get_iers2003
export get_iers2010
export get_iers96
export get_igrf
export get_inclination
export get_intersection_point
export get_model_name
export get_models
export get_normal_gravity
export get_plane_section
export get_polar_radius
export get_pressure
export get_refraction
export get_spin
export get_temperature
export get_theta0
export get_theta_min
export get_total_intensity
export get_undulation
export get_wgs84
export get_wmm
export hash_code
export is_inside
export load_data
export notify
export notify_all
export point_on_limb
export project_to_ground
export set_angular_threshold
export set_pressure
export set_temperature
export still_accepts_data
export supports_time_transform
export to_string
export transform
export transform_model
export valid_from
export valid_to
export wait

include("AtmosphereWrapper/AtmosphereWrapper.jl")
include("DisplacementWrapper/DisplacementWrapper.jl")
include("IonosphereWrapper/IonosphereWrapper.jl")
include("TessellationWrapper/TessellationWrapper.jl")
include("TroposphereWrapper/TroposphereWrapper.jl")
include("WeatherWrapper/WeatherWrapper.jl")
include("earth_it_u453_atmosphere_refraction.jl")
include("earth_shape.jl")
include("earth_standard_atmosphere_refraction.jl")
include("geo_magnetic_elements.jl")
include("geo_magnetic_field.jl")
include("geo_magnetic_field_factory.jl")
include("geo_magnetic_fields.jl")
include("geo_magnetic_model_loader.jl")
include("geoid.jl")
include("lazy_loaded_geo_magnetic_fields.jl")
include("reference_ellipsoid.jl")

end


text:
using Diagonalizations, LinearAlgebra, PosDefManifold, Test

const err=1e-6

# method (1) real
t, n, k=50, 10, 10
A=randn(n, n) # mixing matrix in model x=As
Xset = [genDataMatrix(t, n) for i = 1:k]
Xfixed=randn(t, n)./1
for i=1:length(Xset) Xset[i]+=Xfixed end
Cset = ℍVector([ℍ((Xset[s]'*Xset[s])/t) for s=1:k])
aC=ajd(Cset; algorithm=:OJoB, simple=true)
aC2=ajd(Cset; algorithm=:NoJoB, simple=true)
aC3=ajd(Cset; algorithm=:LogLike, simple=true)
aC4=ajd(Cset; algorithm=:LogLikeR, simple=true)
aC5=ajd(Cset; algorithm=:JADE, simple=true)
aC6=ajd(Cset; algorithm=:JADEmax, simple=true)
aC7=ajd(Cset; algorithm=:GAJD, simple=true)
aC8=ajd(Cset; algorithm=:QNLogLike, simple=true)

# a=ajd(Cset; algorithm=:GAJD2, simple=true, verbose=true)

# method (2) real
aX=ajd(Xset; algorithm=:OJoB, simple=true)
aX2=ajd(Xset; algorithm=:NoJoB, simple=true)
aX3=ajd(Xset; algorithm=:LogLike, simple=true)
aX4=ajd(Xset; algorithm=:LogLikeR, simple=true)
aX5=ajd(Xset; algorithm=:JADE, simple=true)
aX6=ajd(Xset; algorithm=:JADEmax, simple=true)
aX7=ajd(Xset; algorithm=:GAJD, simple=true)
aX8=ajd(Xset; algorithm=:QNLogLike, simple=true)

@test aX≈aC
@test aX2≈aC2
@test aX3≈aC3
@test aX4≈aC4
@test aX5≈aC5
@test aX6≈aC6
@test aX7≈aC7
@test aX8≈aC8

# method (1) complex
t, n, k=50, 10, 10
Ac=randn(ComplexF64, n, n) # mixing matrix in model x=As
Xcset = [genDataMatrix(ComplexF64, t, n) for i = 1:k]
Xcfixed=randn(ComplexF64, t, n)./1
for i=1:length(Xcset) Xcset[i]+=Xcfixed end
Ccset = ℍVector([ℍ((Xcset[s]'*Xcset[s])/t) for s=1:k])
aCc=ajd(Ccset; algorithm=:OJoB, simple=true)
aCc2=ajd(Ccset; algorithm=:NoJoB, simple=true)
aCc3=ajd(Ccset; algorithm=:LogLike, simple=true)
aCc4=ajd(Ccset; algorithm=:JADE, simple=true)
aCc5=ajd(Ccset; algorithm=:JADEmax, simple=true)

# method (2) complex
aXc=ajd(Xcset; algorithm=:OJoB, simple=true)
aXc2=ajd(Xcset; algorithm=:NoJoB, simple=true)
aXc3=ajd(Xcset; algorithm=:LogLike, simple=true)
aXc4=ajd(Xcset; algorithm=:JADE, simple=true)
aXc5=ajd(Xcset; algorithm=:JADEmax, simple=true)

@test aXc≈aCc
@test aXc2≈aCc2
@test aXc3≈aCc3
@test aXc4≈aCc4
@test aXc5≈aCc5

# create 20 REAL random commuting matrices
# they all have the same eigenvectors
Cset2=PosDefManifold.randP(3, 20; eigvalsSNR=Inf, commuting=true)
# estimate the approximate joint diagonalizer (AJD)
a=ajd(Cset2; algorithm=:OJoB)
# the orthogonal AJD must be equivalent to the eigenvector matrix
# of any of the matrices in Cset
@test norm([spForm(a.F'*eigvecs(C)) for C ∈ Cset2])/20 < err
# do the same for JADE algorithm
a=ajd(Cset2; algorithm=:JADE)
@test norm([spForm(a.F'*eigvecs(C)) for C ∈ Cset2])/20 < err
# do the same for JADEmax algorithm
a=ajd(Cset2; algorithm=:JADEmax)
@test norm([spForm(a.F'*eigvecs(C)) for C ∈ Cset2])/20 < err

# generate positive definite matrices with model A*D_κ*D, where
# A is the mixing matrix and D_κ, for all κ=1:k, are diagonal matrices.
# The estimated AJD matrix must be the inverse of A
# and all transformed matrices bust be diagonal
n, k=3, 10
Dest=PosDefManifold.randΛ(eigvalsSNR=10, n, k)
# make the problem identifiable
for i=1:k Dest[k][1, 1]*=i/(k/2) end
for i=1:k Dest[k][3, 3]/=i/(k/2) end
A=randn(n, n) # non-singular mixing matrix
Cset3=Vector{Hermitian}([Hermitian(A*D*A') for D ∈ Dest])
a=ajd(Cset3; algorithm=:NoJoB, eVarC=n)
@test spForm(a.F'*A)<√err
@test mean(nonD(a.F'*Cset3[i]*a.F) for i=1:k)<err
a=ajd(Cset3; algorithm=:LogLike, eVarC=n)
@test spForm(a.F'*A)<√err
@test mean(nonD(a.F'*Cset3[i]*a.F) for i=1:k)<err
a=ajd(Cset3; algorithm=:LogLikeR, eVarC=n)
@test spForm(a.F'*A)<√err
@test mean(nonD(a.F'*Cset3[i]*a.F) for i=1:k)<err
a=ajd(Cset3; algorithm=:GAJD, eVarC=n)
@test spForm(a.F'*A)<√err
@test mean(nonD(a.F'*Cset3[i]*a.F) for i=1:k)<err
a=ajd(Cset3; algorithm=:QNLogLike, eVarC=n)
@test spForm(a.F'*A)<√err
@test mean(nonD(a.F'*Cset3[i]*a.F) for i=1:k)<err

# Do the same thing for orthogonal diagonalizers:
# now A will be orthogonal
O=randU(n) # orthogonal mixing matrix
Cset4=Vector{Hermitian}([Hermitian(O*D*O') for D ∈ Dest])
a=ajd(Cset4; algorithm=:OJoB, eVarC=n)
@test spForm(a.F'*O)<√err
@test mean(nonD(a.F'*Cset4[i]*a.F) for i=1:k)<err
a=ajd(Cset4; algorithm=:JADE, eVarC=n)
@test spForm(a.F'*O)<√err
@test mean(nonD(a.F'*Cset4[i]*a.F) for i=1:k)<√err
a=ajd(Cset4; algorithm=:JADEmax, eVarC=n)
@test spForm(a.F'*O)<√err
@test mean(nonD(a.F'*Cset4[i]*a.F) for i=1:k)<√err

# repeat the test adding noise; now the model is no more exactly identifiable
for k=1:length(Cset3) Cset3[k]+=randP(n)/1000 end
a=ajd(Cset3; algorithm=:NoJoB, eVarC=n)
@test spForm(a.F'*A)<err^(1/6)
@test mean(nonD(a.F'*Cset3[i]*a.F) for i=1:k)<√err
a=ajd(Cset3; algorithm=:LogLike, eVarC=n)
@test spForm(a.F'*A)<err^(1/6)
@test mean(nonD(a.F'*Cset3[i]*a.F) for i=1:k)<√err
a=ajd(Cset3; algorithm=:LogLikeR, eVarC=n)
@test spForm(a.F'*A)<err^(1/6)
@test mean(nonD(a.F'*Cset3[i]*a.F) for i=1:k)<√err
a=ajd(Cset3; algorithm=:GAJD, eVarC=n)
@test spForm(a.F'*A)<err^(1/6)
@test mean(nonD(a.F'*Cset3[i]*a.F) for i=1:k)<√err
a=ajd(Cset3; algorithm=:QNLogLike, eVarC=n)
@test spForm(a.F'*A)<err^(1/6)
@test mean(nonD(a.F'*Cset3[i]*a.F) for i=1:k)<√err

# the same thing for orthogonal diagonalizers
for k=1:length(Cset4) Cset4[k]+=randP(n)/1000 end
a=ajd(Cset4; algorithm=:OJoB, eVarC=n)
@test spForm(a.F'*O)<err^(1/6)
@test mean(nonD(a.F'*Cset4[i]*a.F) for i=1:k)<√err
a=ajd(Cset4; algorithm=:JADE, eVarC=n)
@test spForm(a.F'*O)<err^(1/6)
@test mean(nonD(a.F'*Cset4[i]*a.F) for i=1:k)<√err
a=ajd(Cset4; algorithm=:JADEmax, eVarC=n)
@test spForm(a.F'*O)<err^(1/6)
@test mean(nonD(a.F'*Cset4[i]*a.F) for i=1:k)<√err

# create 20 COMPLEX random commuting matrices
# they all have the same eigenvectors
Ccset2=PosDefManifold.randP(ComplexF64, 3, 20; eigvalsSNR=Inf, commuting=true)
# estimate the approximate joint diagonalizer (AJD)
ac=ajd(Ccset2; algorithm=:OJoB)
# he AJD must be equivalent to the eigenvector matrix of any of the matrices in Cset
# just a sanity check as rounding errors appears for complex data
@test norm([spForm(ac.F'*eigvecs(C)) for C ∈ Ccset2])/20<√err
# do the same for JADE algorithm
ac=ajd(Ccset2; algorithm=:JADE)
@test norm([spForm(ac.F'*eigvecs(C)) for C ∈ Ccset2])/20<√err
# do the same for JADEmax algorithm
ac=ajd(Ccset2; algorithm=:JADEmax)
@test norm([spForm(ac.F'*eigvecs(C)) for C ∈ Ccset2])/3<√err

# the same thing using the NoJoB and LogLike algorithms. Require less precision
# as the NoJoB solution is not constrained in the orthogonal group
ac=ajd(Ccset2; algorithm=:NoJoB)
@test norm([spForm(ac.F'*eigvecs(C)) for C ∈ Ccset2])/20<√err
ac=ajd(Ccset2; algorithm=:LogLike)
@test norm([spForm(ac.F'*eigvecs(C)) for C ∈ Ccset2])/20<√err

# REAL data:
# normalize the trace of input matrices,
# give them weights according to the `nonDiagonality` function
# apply pre-whitening and limit the explained variance both
# at the pre-whitening level and at the level of final vector selection
Cset=PosDefManifold.randP(20, 80; eigvalsSNR=10, SNR=10, commuting=false)

a=ajd(Cset; algorithm=:OJoB, trace1=true, w=nonD, preWhite=true, eVarC=4, eVar=0.99)
a=ajd(Cset; algorithm=:NoJoB, trace1=true, w=nonD, preWhite=true, eVarC=4, eVar=0.99)
a=ajd(Cset; algorithm=:LogLike, w=nonD, preWhite=true, eVarC=4, eVar=0.99)
a=ajd(Cset; algorithm=:LogLikeR, w=nonD, preWhite=true, eVarC=4, eVar=0.99)
a=ajd(Cset; algorithm=:JADE, w=nonD, preWhite=true, eVarC=4, eVar=0.99)
a=ajd(Cset; algorithm=:JADEmax, w=nonD, preWhite=true, eVarC=4, eVar=0.99)
a=ajd(Cset; algorithm=:GAJD, w=nonD, preWhite=true, eVarC=4, eVar=0.99)
a=ajd(Cset; algorithm=:QNLogLike, w=nonD, preWhite=true, eVarC=4, eVar=0.99)

# AJD for plots below
a=ajd(Cset; algorithm=:QNLogLike, verbose=true, preWhite=true)

using Plots
# plot the original covariance matrices
# and their transformed counterpart
CMax=maximum(maximum(abs.(C)) for C ∈ Cset);
 h1 = heatmap(Cset[1], clim=(-CMax, CMax), title="C1", yflip=true, c=:bluesreds);
 h2 = heatmap(Cset[2], clim=(-CMax, CMax), title="C2", yflip=true, c=:bluesreds);
 h3 = heatmap(Cset[3], clim=(-CMax, CMax), title="C3", yflip=true, c=:bluesreds);
 h4 = heatmap(Cset[4], clim=(-CMax, CMax), title="C4", yflip=true, c=:bluesreds);
 📈=plot(h1, h2, h3, h4, size=(700,400))
# savefig(📈, homedir()*"\\Documents\\Code\\julia\\Diagonalizations\\docs\\src\\assets\\FigAJD1.png")

Dset=[a.F'*C*a.F for C ∈ Cset];
 DMax=maximum(maximum(abs.(D)) for D ∈ Dset);
 h5 = heatmap(Dset[1], clim=(-DMax, DMax), title="F'*C1*F", yflip=true, c=:bluesreds);
 h6 = heatmap(Dset[2], clim=(-DMax, DMax), title="F'*C2*F", yflip=true, c=:bluesreds);
 h7 = heatmap(Dset[3], clim=(-DMax, DMax), title="F'*C3*F", yflip=true, c=:bluesreds);
 h8 = heatmap(Dset[4], clim=(-DMax, DMax), title="F'*C4*F", yflip=true, c=:bluesreds);
 📉=plot(h5, h6, h7, h8, size=(700,400))
# savefig(📉, homedir()*"\\Documents\\Code\\julia\\Diagonalizations\\docs\\src\\assets\\FigAJD2.png")

# COMPLEX data:
# normalize the trace of input matrices,
# give them weights according to the `nonDiagonality` function
# apply pre-whitening and limit the explained variance both
# at the pre-whitening level and at the level of final vector selection
Ccset=PosDefManifold.randP(3, 20; eigvalsSNR=10, SNR=2, commuting=false)

ac=ajd(Ccset; trace1=true, w=nonD, preWhite=true,
       algorithm=:OJoB, eVarC=8, eVar=0.99)
ac=ajd(Ccset; eVarC=8, eVar=0.99)
ac=ajd(Ccset; algorithm=:LogLike, eVarC=8, eVar=0.99)
ac=ajd(Ccset; algorithm=:JADE, eVarC=8, eVar=0.99)
ac=ajd(Ccset; algorithm=:JADEmax, eVarC=8, eVar=0.99)


text:
module SetRounding
# Extracted from Julia base. License is MIT: https://julialang.org/license

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

import Base.Rounding:
    setrounding

const Float32_64 = Union{Float32, Float64}

let fenv_consts = Vector{Cint}(undef, 9)
    ccall(:jl_get_fenv_consts, Nothing, (Ptr{Cint},), fenv_consts)
    global const JL_FE_INEXACT = fenv_consts[1]
    global const JL_FE_UNDERFLOW = fenv_consts[2]
    global const JL_FE_OVERFLOW = fenv_consts[3]
    global const JL_FE_DIVBYZERO = fenv_consts[4]
    global const JL_FE_INVALID = fenv_consts[5]

    global const JL_FE_TONEAREST = fenv_consts[6]
    global const JL_FE_UPWARD = fenv_consts[7]
    global const JL_FE_DOWNWARD = fenv_consts[8]
    global const JL_FE_TOWARDZERO = fenv_consts[9]
end

to_fenv(::RoundingMode{:Nearest}) = JL_FE_TONEAREST
to_fenv(::RoundingMode{:ToZero}) = JL_FE_TOWARDZERO
to_fenv(::RoundingMode{:Up}) = JL_FE_UPWARD
to_fenv(::RoundingMode{:Down}) = JL_FE_DOWNWARD

function from_fenv(r::Integer)
    if r == JL_FE_TONEAREST
        return RoundNearest
    elseif r == JL_FE_DOWNWARD
        return RoundDown
    elseif r == JL_FE_UPWARD
        return RoundUp
    elseif r == JL_FE_TOWARDZERO
        return RoundToZero
    else
        throw(ArgumentError("invalid rounding mode code: $r"))
    end
end

setrounding_raw(::Type{<:Float32_64}, i::Integer) = ccall(:fesetround, Int32, (Int32,), i)

@noinline setrounding(::Type{T}, r::RoundingMode) where {T<:Float32_64} = setrounding_raw(T, to_fenv(r))

end #module


text:
module TestChaining

using Test, Random
using DataFrames
using DataFramesMeta
using Lazy, Chain
using Statistics
using Random

Random.seed!(1)
n = 100
df = DataFrame(a = rand(1:3, n),
               b = ["a","b","c","d"][rand(1:4, n)],
               x = rand(n))

x = @where(df, :a .> 2)
x = @transform(x, y = 10 * :x)
x = @by(x, :b, meanX = mean(:x), meanY = mean(:y))
x = @orderby(x, :b, -:meanX)
x = @select(x, var = :b, :meanX, :meanY)

x_as = @as _x_ begin
    df
    @where(_x_, :a .> 2)
    @transform(_x_, y = 10 * :x)
    @by(_x_, :b, meanX = mean(:x), meanY = mean(:y))
    @orderby(_x_, :b, -:meanX)
    @select(_x_, var = :b, :meanX, :meanY)
end

x_thread = @> begin
   df
   @where(:a .> 2)
   @transform(y = 10 * :x)
   @by(:b, meanX = mean(:x), meanY = mean(:y))
   @orderby(:b, -:meanX)
   @select(var = :b, :meanX, :meanY)
end

x_chain = @chain df begin
   @where(:a .> 2)
   @transform(y = 10 * :x)
   @by(:b, meanX = mean(:x), meanY = mean(:y))
   @orderby(:b, -:meanX)
   @select(var = :b, :meanX, :meanY)
end

@test x == x_as
@test x == x_thread

end # module


text:
"""
    AbstractDiffBackend

An abstract type for diff backends. See [`FiniteDifferencesBackend`](@ref) for
an example.
"""
abstract type AbstractDiffBackend end

struct NoneDiffBackend <: AbstractDiffBackend end

"""
    _derivative(f, t[, backend::AbstractDiffBackend])

Compute the derivative of a callable `f` at time `t` computed using the given `backend`,
an object of type [`Manifolds.AbstractDiffBackend`](@ref). If the backend is not explicitly
specified, it is obtained using the function [`Manifolds.diff_backend`](@ref).

This function calculates plain Euclidean derivatives, for Riemannian differentiation see
for example [`differential`](@ref Manifolds.differential(::Manifold, ::Any, ::Real, ::AbstractRiemannianDiffBackend)).

!!! note

    Not specifying the backend explicitly will usually result in a type instability
    and decreased performance.
"""
function _derivative end

function _derivative!(f, X, t, backend::AbstractDiffBackend)
    return copyto!(X, _derivative(f, t, backend))
end

"""
    _gradient(f, p[, backend::AbstractDiffBackend])

Compute the gradient of a callable `f` at point `p` computed using the given `backend`,
an object of type [`AbstractDiffBackend`](@ref). If the backend is not explicitly
specified, it is obtained using the function [`diff_backend`](@ref).

This function calculates plain Euclidean gradients, for Riemannian gradient calculation see
for example [`gradient`](@ref Manifolds.gradient(::Manifold, ::Any, ::Any, ::AbstractRiemannianDiffBackend)).

!!! note

    Not specifying the backend explicitly will usually result in a type instability
    and decreased performance.
"""
function _gradient end

function _gradient!(f, X, p, backend::AbstractDiffBackend)
    return copyto!(X, _gradient(f, p, backend))
end

"""
    CurrentDiffBackend(backend::AbstractDiffBackend)

A mutable struct for storing the current differentiation backend in a global
constant [`_current_diff_backend`](@ref).

# See also

[`AbstractDiffBackend`](@ref), [`diff_backend`](@ref), [`diff_backend!`](@ref)
"""
mutable struct CurrentDiffBackend
    backend::AbstractDiffBackend
end

"""
    _current_diff_backend

The instance of [`Manifolds.CurrentDiffBackend`](@ref) that stores the globally default
differentiation backend.
"""
const _current_diff_backend = CurrentDiffBackend(NoneDiffBackend())

"""
    _diff_backends

A vector of valid [`Manifolds.AbstractDiffBackend`](@ref).
"""
const _diff_backends = AbstractDiffBackend[]

"""
    diff_backend() -> AbstractDiffBackend

Get the current differentiation backend.
"""
diff_backend() = _current_diff_backend.backend

"""
    diff_backend!(backend::AbstractDiffBackend)

Set current backend for differentiation to `backend`.
"""
function diff_backend!(backend::AbstractDiffBackend)
    _current_diff_backend.backend = backend
    return backend
end

"""
    diff_backends() -> Vector{AbstractDiffBackend}

Get vector of currently valid differentiation backends.
"""
diff_backends() = _diff_backends

_derivative(f, t) = _derivative(f, t, diff_backend())

_derivative!(f, X, t) = _derivative!(f, X, t, diff_backend())

_gradient(f, p) = _gradient(f, p, diff_backend())

_gradient!(f, X, p) = _gradient!(f, X, p, diff_backend())

# Finite differences

"""
    FiniteDifferencesBackend(method::FiniteDifferenceMethod = central_fdm(5, 1))

Differentiation backend based on the FiniteDifferences package.
"""
struct FiniteDifferencesBackend{TM<:FiniteDifferenceMethod} <: AbstractDiffBackend
    method::TM
end

function FiniteDifferencesBackend()
    return FiniteDifferencesBackend(central_fdm(5, 1))
end

push!(_diff_backends, FiniteDifferencesBackend())

diff_backend!(_diff_backends[end])

function _derivative(f, t, backend::FiniteDifferencesBackend)
    return backend.method(f, t)
end

function _gradient(f, p, backend::FiniteDifferencesBackend)
    return FiniteDifferences.grad(backend.method, f, p)[1]
end

function _jacobian(f, p, backend::FiniteDifferencesBackend)
    return FiniteDifferences.jacobian(backend.method, f, p)[1]
end


text:
struct AMGXException <: Exception
    e::String
end

function error_string(err_code)
    buf = zeros(UInt8, 1024)
    GC.@preserve buf begin
        cstr = Cstring(pointer(buf))
        v = API.AMGX_get_error_string(err_code, cstr, length(buf))
        if v == API.AMGX_RC_OK
            return unsafe_string(cstr)
        else
            return "error occured when getting error string for error code $(err_code)"
        end
    end
end

macro checked(expr)
    return quote
        v = $(esc(expr))::API.AMGX_RC
        if v != API.AMGX_RC_OK
            throw(AMGXException(error_string(v)))
        end
        nothing
    end
end


text:
# This provides various transformations (GeoData <=> Cartesian; UTMData <=> Cartesian)
#

export ProjectCartData

"""
    d_cart = ProjectCartData(d_cart::CartData, d::GeoData, p::ProjectionPoint)

Projects all datafields from the GeoData struct `d` to the CartData struct `d_cart`, around the projection point `p`.
`d_cart` *must* be an orthogonal cartesian grid (deformed doesn't work; use `Convert2CartData(d, proj)`, where `proj` is a projection point in that case).

# Note:
- If `d_cart` and `d` are horizontal surfaces (3rd dimension has size==1), it also interpolates the depth coordinate.

"""
function ProjectCartData(d_cart::CartData, d::GeoData, p::ProjectionPoint)
    Data_UTM    = Convert2UTMzone(d_cart, p)
    Data_lonlat = convert(GeoData,Data_UTM)

    # Check whether the data sets have the same sign. If not, we may have to shift one by 360 degrees
    min_lon_cart, max_lon_cart = minimum(Data_lonlat.lon.val), maximum(Data_lonlat.lon.val)
    min_lon, max_lon = minimum(d.lon.val), maximum(d.lon.val)

    if (sign(min_lon)!=sign(min_lon_cart)) &&  (sign(max_lon)!=sign(max_lon_cart))
        # the longitude data has a different sign. This can happen if one of them is "West" (and has negative values), whereas the other has
        if (min_lon_cart<0)
            Data_lonlat = GeoData(Data_lonlat.lon.val .+ 360, Data_lonlat.lat.val, ustrip.(Data_lonlat.depth.val), Data_lonlat.fields)
        end
    end

    if size(Data_lonlat.lon.val,3)==1
        z_new, fields_new = InterpolateDataFields2D(d,Data_lonlat.lon.val, Data_lonlat.lat.val)

        # Create new struct
        d_cart = CartData(d_cart.x.val,d_cart.y.val,z_new,fields_new)

    else
        d_data = InterpolateDataFields(d, Data_lonlat.lon.val, Data_lonlat.lat.val, Data_lonlat.depth.val)
        d_cart = CartData(d_cart.x.val,d_cart.y.val,d_cart.z.val,d_data.fields)

    end

    return d_cart
end

"""
    d_cart = ProjectCartData(d_cart::CartData, d::UTMData, p::ProjectionPoint)

Projects all datafields from the UTMData struct `d` to the CartData struct `d_cart`, around the projection point `p`.
    `d_cart` *must* be an orthogonal cartesian grid (deformed doesn't work; use `Convert2CartData(d, proj)`, where `proj` is a projection point in that case).

    # Note:
    - If `d_cart` and `d` are horizontal surfaces (3rd dimension has size==1), it also interpolates the depth coordinate.

"""
function ProjectCartData(d_cart::CartData, d::UTMData, p::ProjectionPoint)
    Data_UTM    = Convert2UTMzone(d_cart, p)

    if size(Data_UTM.EW.val,3)==1
        z_new, fields_new = InterpolateDataFields2D(d,Data_UTM.EW.val, Data_UTM.NS.val)

        # Create new struct
        d_cart = CartData(d_cart.x.val,d_cart.y.val,z_new,fields_new)

    else
        d_data = InterpolateDataFields(d, Data_UTM.EW.val, Data_UTM.NS.val, Data_UTM.depth.val)
        d_cart = CartData(d_cart.x.val,d_cart.y.val,d_cart.z.val,d_data.fields)

    end

    return d_cart
end


text:
@doc raw"""
    IntrinsicCoregionMOKernel(; kernel::Kernel, B::AbstractMatrix)

Kernel associated with the intrinsic coregionalization model.

# Definition

For inputs ``x, x'`` and output dimensions ``p, p'``, the kernel is defined as[^ARL]
```math
k\big((x, p), (x', p'); B, \tilde{k}\big) = B_{p, p'} \tilde{k}\big(x, x'\big),
```
where ``B`` is a positive semidefinite matrix of size ``m \times m``, with ``m`` being the
number of outputs, and ``\tilde{k}`` is a scalar-valued kernel shared by the latent
processes.

[^ARL]: M. Álvarez, L. Rosasco, & N. Lawrence (2012). [Kernels for Vector-Valued Functions: a Review](https://arxiv.org/pdf/1106.6251.pdf).
"""
struct IntrinsicCoregionMOKernel{K<:Kernel,T<:AbstractMatrix} <: MOKernel
    kernel::K
    B::T

    function IntrinsicCoregionMOKernel{K,T}(kernel::K, B::T) where {K,T}
        @check_args(
            IntrinsicCoregionMOKernel,
            B,
            eigmin(B) >= 0,
            "B has to be positive semi-definite"
        )
        return new{K,T}(kernel, B)
    end
end

function IntrinsicCoregionMOKernel(; kernel::Kernel, B::AbstractMatrix)
    return IntrinsicCoregionMOKernel{typeof(kernel),typeof(B)}(kernel, B)
end

function IntrinsicCoregionMOKernel(kernel::Kernel, B::AbstractMatrix)
    return IntrinsicCoregionMOKernel{typeof(kernel),typeof(B)}(kernel, B)
end

function (k::IntrinsicCoregionMOKernel)((x, px)::Tuple{Any,Int}, (y, py)::Tuple{Any,Int})
    return k.B[px, py] * k.kernel(x, y)
end

function _mo_output_covariance(k::IntrinsicCoregionMOKernel, out_dim)
    @assert size(k.B) == (out_dim, out_dim)
    return k.B
end

function kernelmatrix(
    k::IntrinsicCoregionMOKernel, x::MOI, y::MOI
) where {MOI<:IsotopicMOInputsUnion}
    x.out_dim == y.out_dim ||
        throw(DimensionMismatch("`x` and `y` must have the same `out_dim`"))
    Kfeatures = kernelmatrix(k.kernel, x.x, y.x)
    Koutputs = _mo_output_covariance(k, x.out_dim)
    return _kernelmatrix_kron_helper(MOI, Kfeatures, Koutputs)
end

if VERSION >= v"1.6"
    function kernelmatrix!(
        K::AbstractMatrix, k::IntrinsicCoregionMOKernel, x::MOI, y::MOI
    ) where {MOI<:IsotopicMOInputsUnion}
        x.out_dim == y.out_dim ||
            throw(DimensionMismatch("`x` and `y` must have the same `out_dim`"))
        Kfeatures = kernelmatrix(k.kernel, x.x, y.x)
        Koutputs = _mo_output_covariance(k, x.out_dim)
        return _kernelmatrix_kron_helper!(K, MOI, Kfeatures, Koutputs)
    end
end

function Base.show(io::IO, k::IntrinsicCoregionMOKernel)
    return print(
        io, "Intrinsic Coregion Kernel: ", k.kernel, " with ", size(k.B, 1), " outputs"
    )
end


text:
using JLD2, Interpolations
@load joinpath(@__DIR__,"../src/exactcore.jld2") ψi

##
using VortexDistributions


text:
    @testset "Histogram Equalisation" begin

    for T in (Gray{N0f8}, Gray{N0f16}, Gray{Float32}, Gray{Float64})
        #=
        Create an image that spans a narrow graylevel range. Then quantize
        the 256 bins down to 32 and determine how many bins have non-zero
        counts.
        =#

        img = Gray{Float32}.([i/255.0 for i = 64:128, j = 1:10])
        img = T.(img)
        _, counts_before = build_histogram(img, 32, minval = 0, maxval = 1)
        nonzero_before = sum(counts_before .!= 0)

        #=
        Equalize the image histogram. Then quantize the 256 bins down to 32
        and verify that all 32 bins have non-zero counts. This will confirm
        that the dynamic range of the original image has been increased.
        =#
        imgeq = adjust_histogram(img, Equalization(nbins = 256, minval = 0, maxval = 1))
        edges, counts_after = build_histogram(imgeq, 32, minval = 0, maxval = 1)
        nonzero_after = sum(counts_after .!= 0)
        @test nonzero_before < nonzero_after
        @test nonzero_after == 32
    end

    for T in (RGB{N0f8}, RGB{N0f16}, RGB{Float32}, RGB{Float64})
        #=
        Create a color image that spans a narrow graylevel range.  Then
        quantize the 256 bins down to 32 and determine how many bins have
        non-zero counts.
        =#

        imgg = Gray{Float32}.([i/255.0 for i = 64:128, j = 1:10])
        img = colorview(RGB,imgg,imgg,imgg)
        img = T.(img)
        _, counts_before = build_histogram(img, 32, minval = 0, maxval = 1)
        nonzero_before = sum(counts_before .!= 0)

        #=
        Equalize the histogram. Then quantize the 256 bins down to 32 and
        verify that all 32 bins have non-zero counts. This will confirm that
        the dynamic range of the original image has been increased.
        =#
        imgeq = adjust_histogram(img,Equalization(nbins = 256, minval = 0, maxval = 1))
        edges, counts_after = build_histogram(imgeq, 32, minval = 0, maxval = 1)
        nonzero_after = sum(counts_after .!= 0)
        @test nonzero_before < nonzero_after
        @test nonzero_after == 32
    end

    # Verify that the minimum and maximum values of the equalised image match the
    # specified minimum and maximum values, i.e. that the intensities of the equalised
    # image are in the interval [minvalue, maxvalue].
    imgeq = adjust_histogram(collect(0:1:255), Equalization(nbins = 256, minval = 64, maxval = 128))
    @test all(imgeq[1:65] .== 64)
    @test all(imgeq[128+1:end] .== 128)

    imgeq = adjust_histogram(collect(0:1/255:1), Equalization(nbins = 256, minval = 64/255, maxval = 128/255))
    @test all(imgeq[1:65] .== 64/255)
    @test all(imgeq[128+1:end] .== 128/255)
end


text:
using MeasureBase
using SpecialFunctions: besseli
using Random
using KeywordCalls
import Statistics
import Distributions as Dists

export Rice, CPVonMises, CPNormal

@parameterized Rice(ν, σ) ≃ Lebesgue(ℝ₊)

@kwstruct Rice(ν, σ)

function MeasureBase.logdensity(d::Rice{(:ν, :σ)}, x)
    li0 = log(besselix(0, x*d.ν/d.σ^2)) + abs(x*d.ν/d.σ^2)
    log(x/d.σ^2) + li0 - (x^2 + d.ν^2)/(2*d.σ^2)
end

function Base.rand(rng::AbstractRNG, T::Type, d::Rice{(:ν, :σ)})
    x1 = randn(rng, T)*d.σ + d.ν
    x2 = randn(rng, T)*d.σ
    return hypot(x1,x2)
end

#### Statistics
L12(x) = exp(x/2)*( (1-x)*besseli(0,-x/2) - x*besseli(1,-x/2) )
Statistics.mean(d::Rice{(:ν,:σ)}) = d.σ*sqrt(π/2)*L12(-0.5*(d.ν/(d.σ))^2)
Statistics.var(d::Rice{(:ν, :σ)}) = 2*d.σ^2 + d.ν^2 - π*d.σ^2/2*L12(-0.5*(d.ν/(d.σ))^2)^2

@parameterized CPVonMises(μ, κ)

@kwstruct CPVonMises(μ, κ)
@kwstruct CPVonMises(μ, σ)
const log2π = log(2π)
function MeasureBase.logdensity(d::CPVonMises{(:μ, :κ)}, x)
    T = eltype(d.μ)
    sum = zero(T)
    @inbounds @simd for i = eachindex(x)
        sum += d.κ[i]*(cos(x[i]-d.μ[i])-1)
    end
    return sum
end

function MeasureBase.logdensity(d::CPVonMises{(:μ, :σ)},x)
    sum = zero(eltype(d.μ))
    @inbounds @simd for i = eachindex(x)
        sum += (cos(x[i]-d.μ[i])-1)/d.σ[i]^2
    end
    return sum
end

function Base.rand(rng::AbstractRNG, T::Type, d::CPVonMises{(:μ, :κ)})
    d = Dists.VonMises.(d.μ, d.κ)
    return rand.(Ref(rng), Ref(T), d)
end

function Base.rand(rng::AbstractRNG, T::Type, d::CPVonMises{(:μ, :σ)})
    d = @. Dists.VonMises(d.μ, 1/d.σ^2)
    return rand.(Ref(rng), Ref(T), d)
end

@parameterized CPNormal(μ, σ)

Base.minimum(::CPNormal) = -Inf
Base.maximum(::CPNormal) = Inf

const log2π = log(2π)
function MeasureBase.logdensity(dist::CPNormal{(:μ, :σ)}, x::Real)
    μ,σ = dist.μ, dist.σ
    s,c = sincos(x)
    dθ = atan(s, c)
    return  -(abs2(dθ/σ) + log2π)/2
end

function Base.rand(rng::AbstractRNG, T::Type, d::CPNormal)
    return d.μ + d.σ*randn(rng, T)
end


text:
using CSV
using DataFrames

cd("/Users/kiristern/Documents/GitHub/SDM/data/prey_genus/")

# Import and read all saved GBIF csv files from directory
file_list = readdir()
#remove .DS_Store
file_list = file_list[2:end]
#create an empty array of dimension Any
files = Any[]
#fill array with data from downloaded csv files
for i in 1:length(file_list)
    #see which data file is causing the error problems
    global j = i
    #push! into files, info from file_list
    push!(files, CSV.read(file_list[i], delim="\t"))
end
#create one dataframe
df_prey = vcat(files...)

# Select specific columns only
df = df_prey[:, [:genus, :species, :infraspecificEpithet, :taxonRank, :decimalLatitude, :decimalLongitude, :year]]
# create new column "new_sp" for subspecies naming
df.new_sp = copy(df.species)
for i in 1:length(df)
    if df.taxonRank[i] == "SUBSPECIES"
        df.new_sp[i] = string(df.species[i], "_", df.infraspecificEpithet[i])
    end
end

#save csv file
CSV.write("prey_gen_onedf.csv", df)

max_stars_repo_name: kiristern/SDM

max_issues_repo_name: kiristern/SDM

max_forks_repo_name: kiristern/SDM


text:
@doc raw"""
    Euclidean{T<:Tuple,𝔽} <: Manifold

Euclidean vector space.

# Constructor

    Euclidean(n)

Generate the $n$-dimensional vector space $ℝ^n$.

    Euclidean(n₁,n₂,...,nᵢ; field=ℝ)
    𝔽^(n₁,n₂,...,nᵢ) = Euclidean(n₁,n₂,...,nᵢ; field=𝔽)

Generate the vector space of $k = n_1 \cdot n_2 \cdot … \cdot n_i$ values, i.e. the
manifold $𝔽^{n_1, n_2, …, n_i}$ whose
elements are interpreted as $n_1 × n_2 × … × n_i$ arrays.
For $i=2$ we obtain a matrix space.
The default `field=ℝ` can also be set to `field=ℂ`.
The dimension of this space is $k \dim_ℝ 𝔽$, where $\dim_ℝ 𝔽$ is the
[`real_dimension`](@ref) of the field $𝔽$.
"""
struct Euclidean{N,𝔽} <: Manifold where {N<:Tuple, 𝔽<:AbstractNumbers} end

function Euclidean(n::Vararg{Int,I}; field::AbstractNumbers = ℝ) where {I}
    return Euclidean{Tuple{n...},field}()
end

^(𝔽::AbstractNumbers, n) = Euclidean(n...; field = 𝔽)

"""
    EuclideanMetric <: RiemannianMetric

A general type for any manifold that employs the Euclidean Metric, for example
the [`Euclidean`](@ref) manifold itself, or the [`Sphere`](@ref), where every
tangent space (as a plane in the embedding) uses this metric (in the embedding).

Since the metric is independent of the field type, this metric is also used for
the Hermitian metrics, i.e. metrics that are analogous to the `EuclideanMetric`
but where the field type of the manifold is `ℂ`.

This metric is the default metric for example for the [`Euclidean`](@ref) manifold.
"""
struct EuclideanMetric <: RiemannianMetric end

^(M::Euclidean, n::Int) = ^(M, (n,))
function ^(::Euclidean{T,𝔽}, n::NTuple{N,Int}) where {T,𝔽,N}
    return Euclidean{Tuple{T.parameters...,n...},𝔽}()
end

function allocation_promotion_function(
    M::Euclidean{<:Tuple,ℂ},
    ::Union{typeof(get_vector),typeof(get_coordinates)},
    args::Tuple,
)
    return complex
end

function check_manifold_point(M::Euclidean{N,𝔽}, p) where {N,𝔽}
    if (𝔽 === ℝ) && !(eltype(p) <: Real)
        return DomainError(
            eltype(p),
            "The matrix $(p) is not a real-valued matrix, so it does not lie on $(M).",
        )
    end
    if (𝔽 === ℂ) && !(eltype(p) <: Real) && !(eltype(p) <: Complex)
        return DomainError(
            eltype(p),
            "The matrix $(p) is neither a real- nor complex-valued matrix, so it does not lie on $(M).",
        )
    end
    if size(p) != representation_size(M)
        return DomainError(
            size(p),
            "The matrix $(p) does not lie on $(M), since its dimensions ($(size(p))) are wrong (expected: $(representation_size(M))).",
        )
    end
end

function check_tangent_vector(M::Euclidean{N,𝔽}, p, X; check_base_point = true, kwargs...) where {N,𝔽}
    if check_base_point
        mpe = check_manifold_point(M, p; kwargs...)
        mpe === nothing || return mpe
    end
    if (𝔽 === ℝ) && !(eltype(X) <: Real)
        return DomainError(
            eltype(X),
            "The matrix $(X) is not a real-valued matrix, so it can not be a tangent vector to $(p) on $(M).",
        )
    end
    if (𝔽 === ℂ) && !(eltype(X) <: Real) && !(eltype(X) <: Complex)
        return DomainError(
            eltype(X),
            "The matrix $(X) is neither a real- nor complex-valued matrix, so it can not be a tangent vector to $(p) on $(M).",
        )
    end
    if size(X) != representation_size(M)
        return DomainError(
            size(X),
            "The matrix $(X) does not lie in the tangent space of $(p) on $(M), since its dimensions $(size(X)) are wrong  (expected: $(representation_size(M))).",
        )
    end
end

det_local_metric(M::MetricManifold{<:Manifold,EuclideanMetric}, p) = one(eltype(p))

"""
    distance(M::Euclidean, p, q)

Compute the Euclidean distance between two points on the [`Euclidean`](@ref)
manifold `M`, i.e. for vectors it's just the norm of the difference, for matrices
and higher order arrays, the matrix and ternsor Frobenius norm, respectively.
"""
distance(::Euclidean, p, q) = norm(p .- q)

@doc raw"""
    exp(M::Euclidean, p, X)

Compute the exponential map on the [`Euclidean`](@ref) manifold `M` from `p` in direction
`X`, which in this case is just
````math
\exp_p X = p + X.
````
"""
exp(::Euclidean, ::Any...)

exp!(M::Euclidean, q, p, X) = (q .= p .+ X)

"""
    flat(M::Euclidean, p, X)

Transform a tangent vector `X` into a cotangent. Since they can directly be identified in the
[`Euclidean`](@ref) case, this yields just the identity for a tangent vector `w` in the
tangent space of `p` on `M`.
"""
flat(::Euclidean, ::Any...)

flat!(M::Euclidean, ξ::CoTFVector, p, X::TFVector) = copyto!(ξ, X)

function get_basis(M::Euclidean{<:Tuple,ℝ}, p, B::DefaultOrthonormalBasis)
    vecs = [_euclidean_basis_vector(p, i) for i in eachindex(p)]
    return CachedBasis(B,vecs)
end
function get_basis(M::Euclidean{<:Tuple,ℂ}, p, B::DefaultOrthonormalBasis)
    vecs = [_euclidean_basis_vector(p, i) for i in eachindex(p)]
    return CachedBasis(B,[vecs; im * vecs])
end
function get_basis(M::Euclidean, p, B::DiagonalizingOrthonormalBasis)
    vecs = get_vectors(M, p, get_basis(M, p, DefaultOrthonormalBasis()))
    eigenvalues = zeros(real(eltype(p)), manifold_dimension(M))
    return CachedBasis(B, DiagonalizingBasisData(B.frame_direction, eigenvalues, vecs))
end

function get_coordinates!(M::Euclidean{<:Tuple,ℝ}, Y, p, X, B::DefaultOrDiagonalizingBasis)
    S = representation_size(M)
    PS = prod(S)
    copyto!(Y, reshape(X, PS))
    return Y
end
function get_coordinates!(M::Euclidean{<:Tuple,ℂ}, Y, p, X, B::DefaultOrDiagonalizingBasis)
    S = representation_size(M)
    PS = prod(S)
    Y .= [reshape(real(X), PS)..., reshape(imag(X), PS)...]
    return Y
end

function get_vector!(M::Euclidean{<:Tuple,ℝ}, Y, p, X, B::DefaultOrDiagonalizingBasis)
    S = representation_size(M)
    Y .= reshape(X, S)
    return Y
end
function get_vector!(M::Euclidean{<:Tuple,ℂ}, Y, p, X, B::DefaultOrDiagonalizingBasis)
    S = representation_size(M)
    N = div(length(X), 2)
    Y .= reshape(X[1:N] + im * X[N+1:end], S)
    return Y
end

@doc raw"""
    injectivity_radius(M::Euclidean)

Return the injectivity radius on the [`Euclidean`](@ref) `M`, which is $∞$.
"""
injectivity_radius(::Euclidean) = Inf

@doc raw"""
    inner(M::Euclidean, p, X, Y)

Compute the inner product on the [`Euclidean`](@ref) `M`, which is just
the inner product on the real-valued or complex valued vector space
of arrays (or tensors) of size $n_1 × n_2  ×  …  × n_i$, i.e.

````math
g_p(X,Y) = \sum_{k ∈ I} \overline{X}_{k} Y_{k},
````
where $I$ is the set of vectors $k ∈ ℕ^i$, such that for all
$1 ≤ j ≤ i$ it holds $1 ≤ k_j ≤ n_j$.

For the special case of $i ≤ 2$, i.e. matrices and vectors, this simplifies to
````math
g_p(X,Y) = X^{\mathrm{H}}Y,
````
where $\cdot^{\mathrm{H}}$ denotes the Hermitian, i.e. complex conjugate transposed.
"""
inner(::Euclidean, ::Any...)
@inline inner(::Euclidean, p, X, Y) = dot(X, Y)
@inline inner(::MetricManifold{<:Manifold,EuclideanMetric}, p, X, Y) = dot(X, Y)

inverse_local_metric(M::MetricManifold{<:Manifold,EuclideanMetric}, p) = local_metric(M, p)

default_metric_dispatch(::Euclidean, ::EuclideanMetric) = Val(true)

function local_metric(::MetricManifold{<:Manifold,EuclideanMetric}, p)
    return Diagonal(ones(SVector{size(p, 1),eltype(p)}))
end

@doc raw"""
    log(M::Euclidean, p, q)

Compute the logarithmic map on the [`Euclidean`](@ref) `M` from `p` to `q`,
which in this case is just
````math
\log_p q = q-p.
````
"""
log(::Euclidean, ::Any...)

log!(M::Euclidean, X, p, q) = (X .= q .- p)

log_local_metric_density(M::MetricManifold{<:Manifold,EuclideanMetric}, p) = zero(eltype(p))

@generated _product_of_dimensions(::Euclidean{N}) where {N} = prod(N.parameters)

"""
    manifold_dimension(M::Euclidean)

Return the manifold dimension of the [`Euclidean`](@ref) `M`, i.e.
the product of all array dimensions and the [`real_dimension`](@ref) of the
underlying number system.
"""
function manifold_dimension(M::Euclidean{N,𝔽}) where {N,𝔽}
    return _product_of_dimensions(M) * real_dimension(𝔽)
end

mean(::Euclidean{Tuple{1}}, x::AbstractVector{<:Number}; kwargs...) = mean(x)
function mean(
    ::Euclidean{Tuple{1}},
    x::AbstractVector{<:Number},
    w::AbstractWeights;
    kwargs...,
)
    return mean(x, w)
end
mean(::Euclidean, x::AbstractVector; kwargs...) = mean(x)

function mean!(M::Euclidean, p, x::AbstractVector, w::AbstractVector; kwargs...)
    return mean!(M, p, x, w, GeodesicInterpolation(); kwargs...)
end

function mean_and_var(::Euclidean{Tuple{1}}, x::AbstractVector{<:Number}; kwargs...)
    m, v = mean_and_var(x; kwargs...)
    return m, sum(v)
end
function mean_and_var(
    ::Euclidean{Tuple{1}},
    x::AbstractVector{<:Number},
    w::AbstractWeights;
    corrected = false,
    kwargs...,
)
    m, v = mean_and_var(x, w; corrected = corrected, kwargs...)
    return m, sum(v)
end
function mean_and_var(M::Euclidean, x::AbstractVector, w::AbstractWeights; kwargs...)
    return mean_and_var(M, x, w, GeodesicInterpolation(); kwargs...)
end

median(::Euclidean{Tuple{1}}, x::AbstractVector{<:Number}; kwargs...) = median(x)
function median(
    ::Euclidean{Tuple{1}},
    x::AbstractVector{<:Number},
    w::AbstractWeights;
    kwargs...,
)
    return median(x, w)
end

function median!(::Euclidean{Tuple{1}}, p, x::AbstractVector; kwargs...)
    return copyto!(p, [median(vcat(x...))])
end
function median!(::Euclidean{Tuple{1}}, p, x::AbstractVector, w::AbstractWeights; kwargs...)
    return copyto!(p, [median(vcat(x...), w)])
end

@doc raw"""
    norm(M::Euclidean, p, X)

Compute the norm of a tangent vector `X` at `p` on the [`Euclidean`](@ref)
`M`, i.e. since every tangent space can be identified with `M` itself
in this case, just the (Frobenius) norm of `X`.
"""
norm(::Euclidean, p, X) = norm(X)
norm(::MetricManifold{<:Manifold,EuclideanMetric}, p, X) = norm(X)

"""
    normal_tvector_distribution(M::Euclidean, p, σ)

Normal distribution in ambient space with standard deviation `σ`
projected to tangent space at `p`.
"""
function normal_tvector_distribution(M::Euclidean{Tuple{N}}, p, σ) where {N}
    d = Distributions.MvNormal(zero(p), σ)
    return ProjectedFVectorDistribution(TangentBundleFibers(M), p, d, project_vector!, p)
end

@doc raw"""
    project(M::Euclidean, p)

Project an arbitrary point `p` onto the [`Euclidean`](@ref) manifold `M`, which
is of course just the identity map.
"""
project(::Euclidean, ::Any)

project!(M::Euclidean, q, p) = copyto!(q, p)

"""
    project(M::Euclidean, p, X)

Project an arbitrary vector `X` into the tangent space of a point `p` on the
[`Euclidean`](@ref) `M`, which is just the identity, since any tangent
space of `M` can be identified with all of `M`.
"""
project(::Euclidean, ::Any, ::Any)

project!(M::Euclidean, Y, p, X) = copyto!(Y, X)

"""
    projected_distribution(M::Euclidean, d, [p])

Wrap the standard distribution `d` into a manifold-valued distribution. Generated
points will be of similar type to `p`. By default, the type is not changed.
"""
function projected_distribution(M::Euclidean, d, p)
    return ProjectedPointDistribution(M, d, project!, p)
end
function projected_distribution(M::Euclidean, d)
    return ProjectedPointDistribution(M, d, project!, rand(d))
end

"""
    representation_size(M::Euclidean)

Return the array dimensions required to represent an element on the
[`Euclidean`](@ref) `M`, i.e. the vector of all array dimensions.
"""
@generated representation_size(::Euclidean{N}) where {N} = size_to_tuple(N)

"""
    sharp(M::Euclidean, p, ξ)

Transform the cotangent vector `ξ` at `p` on the [`Euclidean`](@ref) `M` to a tangent vector `X`.
Since cotangent and tangent vectors can directly be identified in the [`Euclidean`](@ref)
case, this yields just the identity.
"""
sharp(::Euclidean, ::Any...)

sharp!(M::Euclidean, X::TFVector, p, ξ::CoTFVector) = copyto!(X, ξ)

function show(io::IO, ::Euclidean{N,𝔽}) where {N,𝔽}
    print(io, "Euclidean($(join(N.parameters, ", ")); field = $(𝔽))")
end

"""
    vector_transport_to(M::Euclidean, p, X, q, ::ParallelTransport)

Parallely transport the vector `X` from the tangent space at `p` to the tangent space at `q`
on the [`Euclidean`](@ref) `M`, which simplifies to the identity.
"""
vector_transport_to(::Euclidean, ::Any, ::Any, ::Any, ::ParallelTransport)

vector_transport_to!(M::Euclidean, Y, p, X, q, ::ParallelTransport) = copyto!(Y, X)

var(::Euclidean, x::AbstractVector; kwargs...) = sum(var(x; kwargs...))
function var(::Euclidean, x::AbstractVector{T}, m::T; kwargs...) where {T}
    return sum(var(x; mean = m, kwargs...))
end

"""
    zero_tangent_vector(M::Euclidean, x)

Return the zero vector in the tangent space of `x` on the [`Euclidean`](@ref)
`M`, which here is just a zero filled array the same size as `x`.
"""
zero_tangent_vector(::Euclidean, ::Any...)

zero_tangent_vector!(M::Euclidean, v, x) = fill!(v, 0)


text:
using YaoCompiler
using YaoCompiler.Compiler: optimize
using ZXCalculus
using YaoArrayRegister
using Test

@device function test_cir()
    5 => H
    @ctrl 4 5 => X
    5 => shift($(7 / 4 * π))
    @ctrl 1 5 => X
    5 => T
    @ctrl 4 5 => X
    5 => shift($(7 / 4 * π))
    @ctrl 1 5 => X
    4 => T
    5 => T
    @ctrl 1 4 => X
    4 => shift($(7 / 4 * π))
    1 => T
    @ctrl 1 4 => X
    @ctrl 4 5 => X
    5 => shift($(7 / 4 * π))
    @ctrl 3 5 => X
    5 => T
    @ctrl 4 5 => X
    5 => shift($(7 / 4 * π))
    @ctrl 3 5 => X
    4 => T
    5 => T
    @ctrl 3 4 => X
    4 => shift($(7 / 4 * π))
    5 => H
    3 => T
    @ctrl 3 4 => X
    @ctrl 4 5 => X
    5 => H
    @ctrl 3 5 => X
    5 => shift($(7 / 4 * π))
    @ctrl 2 5 => X
    5 => T
    @ctrl 3 5 => X
    5 => shift($(7 / 4 * π))
    @ctrl 2 5 => X
    3 => T
    5 => T
    @ctrl 2 3 => X
    3 => shift($(7 / 4 * π))
    5 => H
    2 => T
    @ctrl 2 3 => X
    @ctrl 3 5 => X
    5 => H
    @ctrl 2 5 => X
    5 => shift($(7 / 4 * π))
    @ctrl 1 5 => X
    5 => T
    @ctrl 2 5 => X
    5 => shift($(7 / 4 * π))
    @ctrl 1 5 => X
    2 => T
    5 => T
    @ctrl 1 2 => X
    2 => shift($(7 / 4 * π))
    5 => H
    1 => T
    @ctrl 1 2 => X
    @ctrl 2 5 => X
    @ctrl 1 5 => X
end
cir = test_cir()
r = rand_state(5)
r_pt = copy(r)
r_pt_cl = copy(r)
r |> cir

@device optimizer = [:zx_teleport] function teleport_cir()
    5 => H
    @ctrl 4 5 => X
    5 => shift($(7 / 4 * π))
    @ctrl 1 5 => X
    5 => T
    @ctrl 4 5 => X
    5 => shift($(7 / 4 * π))
    @ctrl 1 5 => X
    4 => T
    5 => T
    @ctrl 1 4 => X
    4 => shift($(7 / 4 * π))
    1 => T
    @ctrl 1 4 => X
    @ctrl 4 5 => X
    5 => shift($(7 / 4 * π))
    @ctrl 3 5 => X
    5 => T
    @ctrl 4 5 => X
    5 => shift($(7 / 4 * π))
    @ctrl 3 5 => X
    4 => T
    5 => T
    @ctrl 3 4 => X
    4 => shift($(7 / 4 * π))
    5 => H
    3 => T
    @ctrl 3 4 => X
    @ctrl 4 5 => X
    5 => H
    @ctrl 3 5 => X
    5 => shift($(7 / 4 * π))
    @ctrl 2 5 => X
    5 => T
    @ctrl 3 5 => X
    5 => shift($(7 / 4 * π))
    @ctrl 2 5 => X
    3 => T
    5 => T
    @ctrl 2 3 => X
    3 => shift($(7 / 4 * π))
    5 => H
    2 => T
    @ctrl 2 3 => X
    @ctrl 3 5 => X
    5 => H
    @ctrl 2 5 => X
    5 => shift($(7 / 4 * π))
    @ctrl 1 5 => X
    5 => T
    @ctrl 2 5 => X
    5 => shift($(7 / 4 * π))
    @ctrl 1 5 => X
    2 => T
    5 => T
    @ctrl 1 2 => X
    2 => shift($(7 / 4 * π))
    5 => H
    1 => T
    @ctrl 1 2 => X
    @ctrl 2 5 => X
    @ctrl 1 5 => X
end
tp_cir = teleport_cir()
r_pt |> tp_cir
@test fidelity(r, r_pt) ≈ 1

@device optimizer = [:zx_clifford, :zx_teleport] function clifford_teleport_cir()
    5 => H
    @ctrl 4 5 => X
    5 => shift($(7 / 4 * π))
    @ctrl 1 5 => X
    5 => T
    @ctrl 4 5 => X
    5 => shift($(7 / 4 * π))
    @ctrl 1 5 => X
    4 => T
    5 => T
    @ctrl 1 4 => X
    4 => shift($(7 / 4 * π))
    1 => T
    @ctrl 1 4 => X
    @ctrl 4 5 => X
    5 => shift($(7 / 4 * π))
    @ctrl 3 5 => X
    5 => T
    @ctrl 4 5 => X
    5 => shift($(7 / 4 * π))
    @ctrl 3 5 => X
    4 => T
    5 => T
    @ctrl 3 4 => X
    4 => shift($(7 / 4 * π))
    5 => H
    3 => T
    @ctrl 3 4 => X
    @ctrl 4 5 => X
    5 => H
    @ctrl 3 5 => X
    5 => shift($(7 / 4 * π))
    @ctrl 2 5 => X
    5 => T
    @ctrl 3 5 => X
    5 => shift($(7 / 4 * π))
    @ctrl 2 5 => X
    3 => T
    5 => T
    @ctrl 2 3 => X
    3 => shift($(7 / 4 * π))
    5 => H
    2 => T
    @ctrl 2 3 => X
    @ctrl 3 5 => X
    5 => H
    @ctrl 2 5 => X
    5 => shift($(7 / 4 * π))
    @ctrl 1 5 => X
    5 => T
    @ctrl 2 5 => X
    5 => shift($(7 / 4 * π))
    @ctrl 1 5 => X
    2 => T
    5 => T
    @ctrl 1 2 => X
    2 => shift($(7 / 4 * π))
    5 => H
    1 => T
    @ctrl 1 2 => X
    @ctrl 2 5 => X
    @ctrl 1 5 => X
end
cl_tp_cir = clifford_teleport_cir()
r_pt_cl |> cl_tp_cir
@test fidelity(r, r_pt_cl) ≈ 1

code = @code_yao test_cir()
zxd = ZXDiagram(code)
code2 = optimize(code, [:zx_teleport])
zxd2 = ZXDiagram(code2)
@test tcount(zxd2) == 8


text:
using Documenter, TensorCore

makedocs(;
    modules=[TensorCore],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/JuliaMath/TensorCore.jl/blob/{commit}{path}#L{line}",
    sitename="TensorCore.jl",
    authors="Tim Holy <tim.holy@gmail.com>",
)

deploydocs(;
    repo="github.com/JuliaMath/TensorCore.jl",
    push_preview=true
)


text:
# list of srids for known datums (allows converting LLA / ECEF to SRID points)
# TODO: Check if these are actually correct!

# by default we don't know
SRID{T}(::Type{T}) = error("No known SRID for type $(T).\nPlease overload the SRID constructor to return an SRID for the $(T) point type:\nSRID(::Type{$(T)}) = SRID{..., ...}")
SRID{T}(::T) = SRID(T)

# build the list
known_srids = [

############
# WGS84
############

LLA{WGS84}      :EPSG       4326;    # EPSG code for lon lat wgs84 (GPS).
ECEF{WGS84}     :EPSG       4978;    # EPSG code for lon ecef wgs84 (GPS).

##############
# NAD27
##############

LLA{NAD27}      :EPSG       4267;

##############
# ED50
##############

LLA{ED50}       :EPSG       4230;

##############
# OSGB36
##############

LLA{OSGB36}     :EPSG       4277;

##############
# GDA94
##############

# Australia
LLA{GDA94}     :EPSG        4283;

# New Zealand
LLA{NZGD2000}  :EPSG        4167;

##############
# ETRS89
##############

# Europia
LLA{ETRS89}    :EPSG        4258;

##############
# NAD83
##############

# Americania
LLA{NAD83}     :EPSG        4269

]

# and code gen them up
for i = 1:size(known_srids,1)

    q = quote

        # constructor from type
        auth = symbol($(string(known_srids[i,2])))  # there must be an easier way to insert a raw symbol
        SRID(::Type{$(known_srids[i,1])}) = SRID{auth, $(known_srids[i,3])}

        # constructor from a point
        SRID(::$(known_srids[i,1])) = SRID($(known_srids[i,1]))

        # add a get_srid method
        get_srid(::Type{$(known_srids[i,1])}) = SRID($(known_srids[i,1]))
        get_srid(::$(known_srids[i,1])) = SRID($(known_srids[i,1]))

        # be able to grab the datum as well
        get_datum(::Type{SRID{auth, $(known_srids[i,3])}}) = get_datum($(known_srids[i,1]))
        get_datum(::SRID{auth, $(known_srids[i,3])}) = get_datum($(known_srids[i,1]))

    end
    eval(q)
end

# defaults
# get_srid(x) = UnknownSRID
# SRID(x) = UnknownSRID


text:
# MIT License
# Copyright (c) 2017: Xavier Gandibleux, Anthony Przybylski, Gauthier Soleilhac, and contributors.
function solveWeightedSum(pb::problem, λ1::Int, λ2::Int, lb::Int=0)
    mono_p = mono_problem(pb, λ1, λ2, false)
    mono_s = solve_mono(mono_p,lb)
    return solution(mono_s)
end

#Calcul de X_SE1m par dichotomie
function solveRecursion(xr::solution, xs::solution, X::Vector{solution})

    λ1 = obj_2(xr) - obj_2(xs)
    λ2 = obj_1(xs) - obj_1(xr)

    lb = λ1*xr.obj_1 + λ2*xr.obj_2

    x = solveWeightedSum(xr.pb, λ1, λ2, lb)

    if λ1*x.obj_1 + λ2*x.obj_2 > lb
        push!(X, x)
        solveRecursion(xr, x, X)
        solveRecursion(x, xs, X)
    end
    nothing
end

function first_phase(pb::problem)

    xr = solveWeightedSum(pb, 0, 1)
    xs = solveWeightedSum(pb, 1, 0)
    X = [xr,xs]

    solveRecursion(xr, xs, X)

    sort!(X, by = obj, lt = isless, alg=QuickSort)

    del_inds = Int[]
    for i = 1:length(X)-1
        X[i] >= X[i+1] && push!(del_inds, i+1)
        X[i] <= X[i+1] && push!(del_inds, i)
    end
    deleteat!(X, del_inds)
    return unique(obj, X)
end

max_stars_repo_name: vOptSolver/vOptSpecific

max_issues_repo_name: vOptSolver/vOptSpecific

max_forks_repo_name: vOptSolver/vOptSpecific


text:
"""
Collection of logic for interfacing with HTTP APIs
"""
module APIClients

using Chain
using Dates

using Stonks: DataClientError
using Stonks.Parsers
using Stonks.Models:
  AbstractStonksRecord,
  AssetPrice,
  AssetInfo,
  ExchangeRate,
  IncomeStatement,
  BalanceSheet,
  CashflowStatement,
  Earnings

export APIClient, APIResource, AlphavantageJSONClient, YahooClient

"""
Abstract type to be subtyped by any type responsible for data retrieval, like HTTP API or web socket.
"""
abstract type AbstractDataClient end

"""
    APIResource{T<:AbstractStonksRecord}

Stores data required to make requests to an API resource.

### Fields
- `url::String`: the url of the API resource, excluding query parameters
- `query_params::Dict{String, String}`: parameters used in the request
- `parser::AbstractContentParser`: a subtype of `AbstractContentParser` implementing `parse_content`
- `headers::Dict{String,String}`: HTTP headers
- `symbol_key::String`: indicates the identifier for the symbol
- `max_batch_size::Integer`: the maximum number of symbols allowed in a single request
- `max_retries::Integer`: how many times to retry a failed request
- `rank_order::Integer`: if multiple `APIResource` capable of handling `T`, the one with highest values is preffered.

### Constructors
```julia
APIResource{T}(;
  url::String,
  parser::AbstractContentParser,
  headers::Dict{String, String} = Dict(),
  query_params::Dict{String, String} = Dict(),
  symbol_key::String = max_batch_size > 1 ? "symbols" : "symbol",
  max_batch_size::String = 1,
  max_retries::Integer = 0,
  rank_order::Integer = 1,
)
```
"""
mutable struct APIResource{T<:AbstractStonksRecord}
  url::String
  query_params::Dict{String,String}
  parser::Parsers.AbstractContentParser
  headers::Dict{String,String}
  symbol_key::String
  max_batch_size::Int8
  max_retries::Int
  rank_order::Int8
end
# constructor from kwargs
function APIResource{T}(;
  url,
  query_params=Dict(),
  parser=Parsers.JSONParser(()),
  headers=Dict(),
  symbol_key="",
  max_batch_size=1,
  max_retries=0,
  rank_order=1,
) where {T<:AbstractStonksRecord}
  symbol_key = isempty(symbol_key) && max_batch_size > 1 ? "symbols" : "symbol"
  return APIResource{T}(
    url, query_params, parser, headers, symbol_key, max_batch_size, max_retries, rank_order
  )
end
get_type_param(::APIResource{T}) where {T} = T

"""
Groups collection of `APIResource`. Can hold resources from the same or different API.

### Fields
- `resources::Dict{String,APIResource}`
- `[url::String]`: for descriptive purposes. should be ommited if you mix in resources from different APIs.

### Constructors
```julia
APIClient(resources::Dict{String, APIResource}, url::String="")
```
"""
mutable struct APIClient <: AbstractDataClient
  resources::Dict{String,APIResource}
  url::String
  function APIClient(resources, url="")
    return new(resources, url)
  end
end
get_supported_types(client::APIClient) = [get_type_param(v) for (k, v) in client.resources]

function get_resource(
  client::Union{APIClient,Nothing}, ::Type{T}
) where {T<:AbstractStonksRecord}
  if client !== nothing
    return first([r for (k, r) in client.resources if typeof(r) == APIResource{T}])
  else
    @chain begin
      build_clients_from_env()
      get_resource_from_clients(_, T)
    end
  end
end

"""
    build_clients_from_env() -> Union{Vector{APIClient}, DataClientError}

Utility function which will will try to create instances of `APIClient` from environment variables.
"""
function build_clients_from_env()::Union{Vector{APIClient},Exception}
  clients = APIClient[]
  function try_build_from_keys(keys::Vector{String}, clientBuilder::Function)
    for key in keys
      if haskey(ENV, key)
        !isempty(ENV[key]) && push!(clients, clientBuilder(ENV[key]))
        break
      end
    end
  end
  yahoo_keys = ["YAHOO_APIKEY", "YAHOO_TOKEN", "YAHOOFINANCE_APIKEY", "YAHOOFINANCE_TOKEN"]
  alphavantage_keys = ["ALPHAVANTAGE_APIKEY", "ALPHAVANTAGE_TOKEN"]
  try_build_from_keys(yahoo_keys, YahooClient)
  try_build_from_keys(alphavantage_keys, AlphavantageJSONClient)
  all_keys = k -> join(k, ", ")(vcat(yahoo_keys, alphavantage_keys))
  isempty(clients) && return DataClientError(
    "Could not build API client from environment variables. Please set one of the following: $all_keys",
  )
  return clients
end

"""
    get_resource_from_clients(clients::Vector{APIClient}, ::Type{T}) -> APIResource

Given a list of clients, returns the `APIResource` capable of handling `T`.
If multiple resources are capable:
  - the one with the highest priority (`rank_order`) is returned.
  - if there is a tie based on `rank_order`, one will be chosen at random.
"""
function get_resource_from_clients(
  clients::Union{Vector{APIClient},Exception}, ::Type{T}
) where {T<:AbstractStonksRecord}
  typeof(clients) <: Exception && return clients
  capable_clients = [client for client in clients if T in get_supported_types(client)]
  isempty(capable_clients) &&
    return DataClientError("No client capable of handling type $T")
  length(capable_clients) == 1 && return get_resource(first(capable_clients), T)
  #@debug "Found $(length(capable_clients)) capable clients for handling $T. Will chose them based on provided ranking"
  resources = [get_resource(client, T) for client in capable_clients]
  min_rank = minimum(map(r -> r.rank_order, resources))
  res = filter(r -> r.rank_order == min_rank, resources)
  length(res) == 1 && return first(res)
  #@debug "Found multiple resources with the same rank: $min_rank. WIll chose one at random"
  return res[rand(1:length(res))]
end

"""
    YahooClient(api_key::String) -> APIClient

Utility function for creating a client for accessing yahoofinance API.
Contains the following `resources`:
  - info => `APIResource{AssetInfo}`
  - price => `APIResource{AssetPrice}`
  - exchange => `APIResource{ExchangeRate}`
  - income_statement => `APIResource{IncomeStatement}`
  - balance_sheet => `APIResource{BalanceSheet}`
  - cashflow_statement => `APIResource{CashflowStatement}`
  - earnings => `APIResource{Earnings}`
"""
function YahooClient(api_key::String)::APIClient
  url = "https://yfapi.net"
  headers = Dict("accept" => "application/json", "X-API-KEY" => api_key)
  price = APIResource{AssetPrice}(;
    url="$url/v8/finance/spark",
    query_params=Dict(
      "interval" => "{interval}", "range" => "{range}", "symbols" => "{symbols}"
    ),
    parser=Parsers.YahooPriceParser,
    headers=headers,
    max_batch_size=10,
    max_retries=1,
  )
  info = APIResource{AssetInfo}(;
    url="$url/v11/finance/quoteSummary/{symbol}",
    query_params=Dict("modules" => "assetProfile,quoteType,price"),
    parser=Parsers.YahooInfoParser,
    headers=headers,
    max_retries=1,
  )
  exchange = APIResource{ExchangeRate}(;
    url="$url/v8/finance/spark",
    query_params=Dict(
      "interval" => "{interval}", "range" => "{range}", "symbols" => "{symbols}"
    ),
    parser=Parsers.YahooExchangeRateParser,
    headers=headers,
    max_batch_size=10,
    max_retries=1,
  )
  income_statement = APIResource{IncomeStatement}(;
    url="$url/v11/finance/quoteSummary/{symbol}",
    query_params=Dict(
      "modules" => "price,incomeStatementHistory,incomeStatementHistoryQuarterly"
    ),
    parser=Parsers.YahooIncomeStatementParser,
    headers=headers,
    max_retries=1,
    rank_order=2,
  )
  balance_sheet = APIResource{BalanceSheet}(;
    url="$url/v11/finance/quoteSummary/{symbol}",
    query_params=Dict(
      "modules" => "price,balanceSheetHistory,balanceSheetHistoryQuarterly"
    ),
    parser=Parsers.YahooBalanceSheetParser,
    headers=headers,
    max_retries=1,
    rank_order=2,
  )
  cashflow_statement = APIResource{CashflowStatement}(;
    url="$url/v11/finance/quoteSummary/{symbol}",
    query_params=Dict(
      "modules" => "price,cashflowStatementHistory,cashflowStatementHistoryQuarterly"
    ),
    parser=Parsers.YahooCashflowParser,
    headers=headers,
    max_retries=1,
    rank_order=2,
  )
  earnings = APIResource{Earnings}(;
    url="$url/v11/finance/quoteSummary/{symbol}",
    query_params=Dict("modules" => "price,earningsHistory"),
    parser=Parsers.YahooEarningsParser,
    headers=headers,
    max_retries=1,
    rank_order=2,
  )
  resources = Dict(
    "price" => price,
    "info" => info,
    "exchange" => exchange,
    "income_statement" => income_statement,
    "balance_sheet" => balance_sheet,
    "cashflow_statement" => cashflow_statement,
    "earnings" => earnings,
  )
  return APIClient(resources, url)
end

"""
    AlphavantageJSONClient(api_key::String) -> APIClient

Utility function for creating a client for accessing alphavantage API.
Contains the following `resources`:
  - info => `APIResource{AssetInfo}`
  - price => `APIResource{AssetPrice}`
  - exchange => `APIResource{ExchangeRate}`
  - income_statement => `APIResource{IncomeStatement}`
  - balance_sheet => `APIResource{BalanceSheet}`
  - cashflow_statement => `APIResource{CashflowStatement}`
  - earnings => `APIResource{Earnings}`
"""
function AlphavantageJSONClient(api_key::String)::APIClient
  url = "https://www.alphavantage.co"
  headers = Dict("accept" => "application/json")
  price = APIResource{AssetPrice}(;
    url="$url/query",
    query_params=Dict(
      "function" => "TIME_SERIES_DAILY", "symbol" => "{symbol}", "apikey" => api_key
    ),
    parser=Parsers.AlphavantagePriceParser,
    headers=headers,
    max_retries=1,
    rank_order=2,
  )
  info = APIResource{AssetInfo}(;
    url="$url/query",
    query_params=Dict(
      "function" => "OVERVIEW", "symbol" => "{symbol}", "apikey" => api_key
    ),
    parser=Parsers.AlphavantageInfoParser,
    headers=headers,
    max_retries=1,
    rank_order=2,
  )
  exchange = APIResource{ExchangeRate}(;
    url="$url/query",
    query_params=Dict(
      "function" => "FX_DAILY",
      "from_symbol" => "{base}",
      "to_symbol" => "{target}",
      "apikey" => api_key,
    ),
    parser=Parsers.AlphavantageExchangeRateParser,
    headers=headers,
    max_retries=1,
    rank_order=2,
  )
  income_statement = APIResource{IncomeStatement}(;
    url="$url/query",
    query_params=Dict(
      "function" => "INCOME_STATEMENT", "symbol" => "{symbol}", "apikey" => api_key
    ),
    parser=Parsers.AlphavantageIncomeStatementParser,
    headers=headers,
    max_retries=1,
    rank_order=1,
  )
  balance_sheet = APIResource{BalanceSheet}(;
    url="$url/query",
    query_params=Dict(
      "function" => "BALANCE_SHEET", "symbol" => "{symbol}", "apikey" => api_key
    ),
    parser=Parsers.AlphavantageBalanceSheetParser,
    headers=headers,
    max_retries=1,
    rank_order=1,
  )
  cashflow_statement = APIResource{CashflowStatement}(;
    url="$url/query",
    query_params=Dict(
      "function" => "CASH_FLOW", "symbol" => "{symbol}", "apikey" => api_key
    ),
    parser=Parsers.AlphavantageCashflowParser,
    headers=headers,
    max_retries=1,
    rank_order=1,
  )
  earnings = APIResource{Earnings}(;
    url="$url/query",
    query_params=Dict(
      "function" => "EARNINGS", "symbol" => "{symbol}", "apikey" => api_key
    ),
    parser=Parsers.AlphavantageEarningsParser,
    headers=headers,
    max_retries=1,
    rank_order=1,
  )
  resources = Dict(
    "price" => price,
    "info" => info,
    "exchange" => exchange,
    "income_statement" => income_statement,
    "balance_sheet" => balance_sheet,
    "cashflow_statement" => cashflow_statement,
    "earnings" => earnings,
  )
  return APIClient(resources, url)
end

end


text:
module TestComposites

using Test
using MLJBase
using ..Models
using CategoricalArrays
import Random.seed!
seed!(1234)

@load KNNRegressor

N = 50
Xin = (a=rand(N), b=rand(N), c=rand(N))
yin = rand(N)

train, test = partition(eachindex(yin), 0.7);
Xtrain = MLJBase.selectrows(Xin, train)
ytrain = yin[train]

ridge_model = FooBarRegressor(lambda=0.1)
selector_model = FeatureSelector()

@testset "first test of hand-exported network" begin
    composite = SimpleDeterministicCompositeModel(model=ridge_model,
                                                  transformer=selector_model)

    fitresult, cache, report = MLJBase.fit(composite, 3, Xtrain, ytrain)

    # to check internals:
    ridge = MLJBase.machines(fitresult)[1]
    selector = MLJBase.machines(fitresult)[2]
    ridge_old = deepcopy(ridge)
    selector_old = deepcopy(selector)

    # this should trigger no retraining:
    fitresult, cache, report =
        @test_logs(
            (:info, r"^Not"),
            (:info, r"^Not"),
            MLJBase.update(composite, 2, fitresult, cache, Xtrain, ytrain))
    @test ridge.fitresult == ridge_old.fitresult
    @test selector.fitresult == selector_old.fitresult

    # this should trigger update of selector and training of ridge:
    selector_model.features = [:a, :b]
    fitresult, cache, report =
        @test_logs(
            (:info, r"^Updating"),
            (:info, r"^Training"),
            MLJBase.update(composite, 2, fitresult, cache, Xtrain, ytrain))
    @test ridge.fitresult != ridge_old.fitresult
    @test selector.fitresult != selector_old.fitresult
    ridge_old = deepcopy(ridge)
    selector_old = deepcopy(selector)

    # this should trigger updating of ridge only:
    ridge_model.lambda = 1.0
    fitresult, cache, report =
        @test_logs(
            (:info, r"^Not"),
            (:info, r"^Updating"),
            MLJBase.update(composite, 2, fitresult, cache, Xtrain, ytrain))
    @test ridge.fitresult != ridge_old.fitresult
    @test selector.fitresult == selector_old.fitresult

    predict(composite, fitresult, MLJBase.selectrows(Xin, test))

    Xs = source(Xtrain)
    ys = source(ytrain, kind=:target)

    mach = machine(composite, Xs, ys)
    yhat = predict(mach, Xs)
    fit!(yhat, verbosity=3)
    composite.transformer.features = [:b, :c]
    fit!(yhat, verbosity=3)
    fit!(yhat, rows=1:20, verbosity=3)
    yhat(MLJBase.selectrows(Xin, test))

end

mutable struct WrappedRidge <: DeterministicNetwork
    ridge
end

@testset "second test of hand-exported network" begin

    function MLJBase.fit(model::WrappedRidge, verbosity::Integer, X, y)
        Xs = source(X)
        ys = source(y, kind=:target)

        stand = Standardizer()
        standM = machine(stand, Xs)
        W = transform(standM, Xs)

        boxcox = UnivariateBoxCoxTransformer()
        boxcoxM = machine(boxcox, ys)
        z = transform(boxcoxM, ys)

        ridgeM = machine(model.ridge, W, z)
        zhat = predict(ridgeM, W)
        yhat = inverse_transform(boxcoxM, zhat)

        fit!(yhat)
        return fitresults(yhat)
    end

    MLJBase.input_scitype(::Type{<:WrappedRidge}) = Table(Continuous)
    MLJBase.target_scitype(::Type{<:WrappedRidge}) = AbstractVector{<:Continuous}

    ridge = FooBarRegressor(lambda=0.1)
    model_ = WrappedRidge(ridge)
    mach = machine(model_, Xin, yin)
    fit!(mach)
    yhat=predict(mach, Xin)
    ridge.lambda = 1.0
    fit!(mach)
    @test predict(mach, Xin) != yhat

end

@load DecisionTreeRegressor
@load DecisionTreeClassifier

## FROM_NETWORK_PREPROCESS

# supervised:
Xs = source(nothing)
ys = source(nothing, kind=:target)
z  = log(ys)
stand = UnivariateStandardizer()
standM = machine(stand, z)
u = transform(standM, z)
hot = OneHotEncoder()
hotM = machine(hot, Xs)
W = transform(hotM, Xs)
knn = KNNRegressor()
knnM = machine(knn, W, u)
oak = DecisionTreeRegressor()
oakM = machine(oak, W, u)
uhat = 0.5*(predict(knnM, W) + predict(oakM, W))
zhat = inverse_transform(standM, uhat)
yhat = exp(zhat)

ex = Meta.parse("Composite(knn_rgs=knn, one_hot_enc=hot) <= yhat")
modeltype_ex, fieldname_exs, model_exs, N_ex, kind, trait_dic =
    MLJBase.from_network_preprocess(TestComposites, ex)
@test modeltype_ex == :Composite
@test fieldname_exs == [:knn_rgs, :one_hot_enc]
@test model_exs == [:knn, :hot]
@test N_ex == :yhat
@test kind == :DeterministicNetwork
@test !(haskey(trait_dic, :supports_weights)) ||
    !trait_dic[:supports_weights]

# supervised with sample weights:
ws = source(kind=:weights)
knnM = machine(knn, W, u, ws)
uhat = 0.5*(predict(knnM, W) + predict(oakM, W))
zhat = inverse_transform(standM, uhat)
yhat = exp(zhat)
ex = Meta.parse("Composite(knn_rgs=knn, one_hot_enc=hot) <= yhat")
modeltype_ex, fieldname_exs, model_exs, N_ex, kind, trait_dic =
    MLJBase.from_network_preprocess(TestComposites, ex)
@test trait_dic[:supports_weights]

# unsupervised:
ex = Meta.parse("Composite(one_hot_enc=hot) <= W")
modeltype_ex, fieldname_exs, model_exs, N_ex, kind, trait_dic =
    MLJBase.from_network_preprocess(TestComposites, ex)
@test modeltype_ex == :Composite
@test fieldname_exs == [:one_hot_enc,]
@test model_exs == [:hot,]
@test N_ex == :W
@test kind == :UnsupervisedNetwork
@test !(haskey(trait_dic, :supports_weights))

# second supervised test:
fea = FeatureSelector()
feaM = machine(fea, Xs)
G = transform(feaM, Xs)
hotM = machine(hot, G)
H = transform(hotM, G)
elm = DecisionTreeClassifier()
elmM = machine(elm, H, ys)
yhat = predict(elmM, H)

ex = Meta.parse("Composite(selector=fea,one_hot=hot,tree=elm) <= yhat")
modeltype_ex, fieldname_exs, model_exs, N_ex, kind, trait_dic =
    MLJBase.from_network_preprocess(TestComposites,
                                ex, :(prediction_type=:probabilistic))
@test modeltype_ex == :Composite
@test fieldname_exs == [:selector, :one_hot, :tree]
@test model_exs == [:fea, :hot, :elm]
@test N_ex == :yhat
@test kind == :ProbabilisticNetwork
@test !(haskey(trait_dic, :supports_weights)) ||
    !trait_dic[:supports_weights]

ex = Meta.parse("45")
@test_throws(ArgumentError,
             MLJBase.from_network_preprocess(TestComposites,
                                         ex, :(prediction_type=:probabilistic)))

ex = Meta.parse("Composite(elm=elm) << yhat")
@test_throws(ArgumentError,
             MLJBase.from_network_preprocess(TestComposites,
                                         ex, :(prediction_type=:probabilistic)))

ex = Meta.parse("45")
@test_throws(ArgumentError,
             MLJBase.from_network_preprocess(TestComposites,
                                         ex, :(prediction_type=:probabilistic)))

ex = Meta.parse("45 <= yhat")
@test_throws(ArgumentError,
             MLJBase.from_network_preprocess(TestComposites,
                                         ex, :(prediction_type=:probabilistic)))
ex = Meta.parse("Comp(elm=45) <= yhat")
@test_throws(ArgumentError,
             MLJBase.from_network_preprocess(TestComposites,
                                         ex, :(prediction_type=:probabilistic)))

ex = Meta.parse("Comp(elm=>elm) <= yhat")
@test_throws(ArgumentError,
             MLJBase.from_network_preprocess(TestComposites,
                                         ex, :(prediction_type=:probabilistic)))

ex = Meta.parse("Comp(34=elm) <= yhat")
@test_throws(ArgumentError,
             MLJBase.from_network_preprocess(TestComposites,
                                         ex, :(prediction_type=:probabilistic)))

ex = Meta.parse("Comp(elm=elm) <= 45")
@test_throws(ArgumentError,
             MLJBase.from_network_preprocess(TestComposites,
                                         ex, :(prediction_type=:probabilistic)))

z = vcat(ys, ys)
ex = Meta.parse("Comp() <= z")
@test_throws(ArgumentError,
             MLJBase.from_network_preprocess(TestComposites,ex))

X2s = source(nothing)
# z = @test_logs (:warn, r"^A node ref") vcat(Xs, X2s)
z = vcat(Xs, X2s)
ex = Meta.parse("Comp() <= z")
@test_throws(ArgumentError,
             MLJBase.from_network_preprocess(TestComposites, ex))

y2s = source(nothing, kind=:target)
# z = @test_logs (:warn, r"^A node ref") vcat(ys, y2s, Xs)
z = vcat(ys, y2s, Xs)
@test_throws Exception z(Xs())

ex = Meta.parse("Comp() <= z")
@test_throws(ArgumentError,
             MLJBase.from_network_preprocess(TestComposites, ex))

ex = Meta.parse("Composite(one_hot_enc=hot) <= W")
@test_throws(ArgumentError,
             MLJBase.from_network_preprocess(TestComposites,
                                         ex, :(prediction_type=:probabilistic)))

## TEST MACRO-EXPORTED SUPERVISED NETWORK
# (CANNOT WRAP IN @testset)

x1 = map(n -> mod(n,3), rand(UInt8, 100)) |> categorical;
x2 = randn(100);
X = (x1=x1, x2=x2);
y = x2.^2;

Xs = source(X)
ys = source(y, kind=:target)
z = log(ys)
stand = UnivariateStandardizer()
standM = machine(stand, z)
u = transform(standM, z)
hot = OneHotEncoder()
hotM = machine(hot, Xs)
W = transform(hotM, Xs)
knn = KNNRegressor()
knnM = machine(knn, W, u)
oak = DecisionTreeRegressor()
oakM = machine(oak, W, u)
uhat = 0.5*(predict(knnM, W) + predict(oakM, W))
zhat = inverse_transform(standM, uhat)
yhat = exp(zhat)

# test that state changes after fit:
@test sum(MLJBase.state(yhat) |> MLJBase.flat_values) == 0
fit!(yhat)
@test sum(MLJBase.state(W) |> MLJBase.flat_values) == 1

# test nested reporting:
r = MLJBase.report(yhat)
@test r isa NamedTuple
@test length(r.reports) == 4
@test r.reports[1] == NamedTuple()

hot2 = deepcopy(hot)
knn2 = deepcopy(knn)
ys2 = source(nothing, kind=:target)

# duplicate a network:
yhat2 = replace(yhat, hot=>hot2, knn=>knn2,
                ys=>source(ys.data, kind=:target);
                empty_unspecified_sources=true)
@test isempty(sources(yhat2, kind=:input)[1])
yhat2 = @test_logs((:warn, r"No replacement"),
                   replace(yhat, hot=>hot2, knn=>knn2,
                           ys=>source(ys.data, kind=:target)))
@test !isempty(sources(yhat2, kind=:input)[1])

@test_logs((:info, r"^Train.*OneHot"),
           (:info, r"^Spawn"),
           (:info, r"^Train.*Univ"),
           (:info, r"^Train.*KNN"),
           (:info, r"^Train.*Dec"), fit!(yhat2))
@test length(MLJBase.machines(yhat)) == length(MLJBase.machines(yhat2))
@test models(yhat) == models(yhat2)
@test sources(yhat) == sources(yhat2)
@test MLJBase.tree(yhat) == MLJBase.tree(yhat2)
@test yhat() ≈ yhat2()

# this change should trigger retraining of all machines except the
# univariate standardizer:
hot2.drop_last = true
@test_logs((:info, r"^Updating.*OneHot"),
           (:info, r"^Spawn"),
           (:info, r"^Not.*Univ"),
           (:info, r"^Train.*KNN"),
           (:info, r"^Train.*Dec"), fit!(yhat2))

# export a supervised network:
model_ = @from_network Composite(knn_rgs=knn, one_hot_enc=hot) <= yhat

mach = machine(model_, X, y)
@test_logs((:info, r"^Train.*Composite"),
           (:info, r"^Train.*OneHot"),
           (:info, r"^Spawn"),
           (:info, r"^Train.*Univ"),
           (:info, r"^Train.*KNN"),
           (:info, r"^Train.*Dec"), fit!(mach))
model_.knn_rgs.K = 55
@test_logs((:info, r"^Updat.*Composite"),
           (:info, r"^Not.*OneHot"),
           (:info, r"^Not.*Univ"),
           (:info, r"^Updat.*KNN"),
           (:info, r"^Not.*Dec"), fit!(mach))

@test MLJBase.tree(mach.fitresult).arg1.arg1.arg1.arg1.model.K == 55

# check data anomynity:
@test all(x->(x===nothing), [s.data for s in sources(mach.fitresult)])

## TEST MACRO-EXPORTED UNSUPERVISED NETWORK
# (CANNOT WRAP IN @testset)

multistand = Standardizer()
multistandM = machine(multistand, W)
W2 = transform(multistandM, W)
model_ = @from_network Transf(one_hot=hot) <= W2
mach = machine(model_, X)
@test_logs((:info, r"^Training.*Transf"),
               (:info, r"^Train.*OneHot"),
               (:info, r"^Spawn"),
               (:info, r"Train.*Stand"), fit!(mach))
model_.one_hot.drop_last=true
@test_logs((:info, r"^Updating.*Transf"),
               (:info, r"^Updating.*OneHot"),
               (:info, r"^Spawn"),
               (:info, r"Train.*Stand"), fit!(mach))
@test(fitted_params(mach).fitted_params[1] isa
          NamedTuple{(:mean_and_std_given_feature,)})

# check data anomynity:
@test all(x->(x===nothing), [s.data for s in sources(mach.fitresult)])

transform(mach)

## TEST MACRO-EXPORTED SUPERVISED NETWORK WITH SAMPLE WEIGHTS

seed!(1234)
N = 100
X = (x = rand(3N), );
y = categorical(rand("abc", 3N));
# define class weights :a, :b, :c in ration 2:4:1
w = map(y) do η
    if η == 'a'
        return 2
    elseif η == 'b'
        return 4
    else
        return 1
    end
end;
Xs = source(X)
ys = source(y, kind=:target)
ws = source(w, kind=:weights)

standM = machine(Standardizer(), Xs)
W = transform(standM, Xs)

rgs = ConstantClassifier() # supports weights
rgsM = machine(rgs, W, ys, ws)
yhat = predict(rgsM, W)

fit!(yhat)
fit!(yhat, rows=1:div(N,2))
yhat(rows=1:div(N,2));

composite = @from_network Composite3(regressor=rgs) <= yhat

@test MLJBase.supports_weights(composite)
mach = fit!(machine(composite, X, y))
predict(mach, rows=1:div(N,2))[1]
posterior = predict(mach, rows=1:div(N,2))[1]

# "posterior" is roughly uniform:
@test abs(pdf(posterior, 'b')/(pdf(posterior, 'a'))  - 1) < 0.15
@test abs(pdf(posterior, 'b')/(pdf(posterior, 'c'))  - 1) < 0.15

# now add weights:
mach = fit!(machine(composite, X, y, w), rows=1:div(N,2))
posterior = predict(mach, rows=1:div(N,2))[1]

# "posterior" is skewed appropriately in weighted case:
@test abs(pdf(posterior, 'b')/(2*pdf(posterior, 'a'))  - 1) < 0.15
@test abs(pdf(posterior, 'b')/(4*pdf(posterior, 'c'))  - 1) < 0.15

end
true


text:
module WelcomeToJulia

end # module


text:
@testset "testReactCode." begin
    params = ones(15) * 0.5

    @testset "Can successfully assemble the parameters." begin
        TAMode.runTAM(tps, params, 10.0)
    end

    function gradFunc(x)
        return dot(TAMode.runTAM([10.0], x, 1.0), TAMode.pY)
    end

    @testset "Check that we can diff though the solver." begin
        g = similar(params)

        ForwardDiff.gradient!(g, gradFunc, params)
        @time ForwardDiff.gradient!(g, gradFunc, params)

        @test(length(g) == length(params))
    end

    @testset "Make sure code upholds mass conservation." begin
        tt = TAMode.param(params)
        firstV = TAMode.getAutocrine(tt)

        tt.kDeg = 0
        tt.TAMs[1].expression = 0
        tt.TAMs[2].expression = 0
        tt.TAMs[3].expression = 0
        tt.gasCur *= 1000

        secondV = TAMode.runTAMinit([1000000.0], tt, firstV)

        @test dot(firstV, TAMode.total) - dot(secondV, TAMode.total) < 0.0001
    end

    @testset "Test that receptor and ligand amounts match expectations." begin
        tt = TAMode.param(params)

        tt.gasCur = 0.0
        tt.TAMs[1].expression = 10.0
        tt.TAMs[2].expression = 10.0
        tt.TAMs[3].expression = 10.0
        tt.kRec = 0.0
        tt.internalize = 10.0

        outt = TAMode.getAutocrine(tt)

        # Expect no ligand
        @test outt[end] ≈ 0.0

        for i = 1:3
            @test isapprox(dot(outt, TAMode.recpSpecific[i] .* TAMode.total .* TAMode.surface), 1.0, rtol = 1.0e-5)
        end
    end

    @testset "Check for detailed balance." begin
        # PMID: 16698778
        rr = TAMode.param(params)
        autoC = TAMode.getAutocrine(rr)

        rr.TAMs[1].expression = 0.0
        rr.TAMs[2].expression = 0.0
        rr.TAMs[3].expression = 0.0
        rr.kDeg = 0.0

        uLong = TAMode.runTAMinit([1000000.0], rr, autoC)

        # Get the Jacobian matrix
        du = zero(uLong)
        J = ForwardDiff.jacobian((y, x) -> TAMode.TAMreact(y, x, rr, 0.0), du, uLong)
        GK = J * diagm(vec(uLong))

        @test norm(GK - transpose(GK)) < 1.0e-5
    end

    @testset "Swapping twice leads to the same rates." begin
        rr = TAMode.param(params)
        rrB = TAMode.swapIgs(TAMode.swapIgs(rr))

        for ii = 1:3
            @test all(rr.TAMs[ii].xRev .≈ rrB.TAMs[ii].xRev)
            @test all(rr.TAMs[ii].binding .≈ rrB.TAMs[ii].binding)
            @test all(rr.hetR[ii].xRev .≈ rrB.hetR[ii].xRev)
        end
    end

    @testset "Swapping Ig1 with Ig2 doesn't change anything." begin
        rr = TAMode.param(params)
        rrSwap = TAMode.swapIgs(rr)

        dataDiff = TAMode.runTAMinit(tps, rr, zeros(55)) .- TAMode.runTAMinit(tps, rrSwap, zeros(55))
        @test all(aboutZero.(dataDiff * TAMode.pY))
        @test all(aboutZero.(dataDiff * TAMode.total))
        @test all(aboutZero.(dataDiff * TAMode.surface))

        for ii = 1:3
            @test all(aboutZero.(dataDiff * TAMode.recpSpecific[ii]))
        end

        autoDiff = TAMode.getAutocrine(rr) .- TAMode.getAutocrine(rrSwap)
        @test aboutZero(dot(autoDiff, TAMode.pY))
        @test aboutZero(dot(autoDiff, TAMode.total))
        @test aboutZero(dot(autoDiff, TAMode.surface))

        for ii = 1:3
            @test aboutZero(dot(autoDiff, TAMode.recpSpecific[ii]))
        end

        dataDiff = TAMode.runTAM(tps, rr, 1.0) .- TAMode.runTAM(tps, rrSwap, 1.0)
        @test all(aboutZero.(dataDiff * TAMode.pY))
        @test all(aboutZero.(dataDiff * TAMode.total))
        @test all(aboutZero.(dataDiff * TAMode.surface))
        @test all(aboutZero.(dataDiff * TAMode.boundLig))

        for ii = 1:3
            @test all(aboutZero.(dataDiff * TAMode.recpSpecific[ii]))
        end
    end

    @testset "Test that if none of the ligand is expressed, we don't end up seeing any." begin
        rr = TAMode.param(params)
        rr.gasCur = 0.0

        data = TAMode.runTAM(tps, rr, 0.0)

        @test all(data * TAMode.pY .≈ 0)
        @test all(data * TAMode.boundLig .≈ 0)
    end

    @testset "Check if there's no receptor that we don't see any." begin
        for i = 1:3
            rr = TAMode.param(params)
            rr.TAMs[i].expression = 0.0

            data = TAMode.runTAM(tps, rr, 1.0)

            @test all(data * TAMode.recpSpecific[i] .≈ 0.0)
        end
    end

    @testset "Make sure that TAM surface don't explode at long time in reaction code." begin
        tt = TAMode.param(params)
        firstSurf = TAMode.getAutocrine(tt)

        tt.TAMs[1].expression = 0.0
        tt.TAMs[2].expression = 0.0
        tt.TAMs[3].expression = 0.0
        tt.kRec = 0.0
        tt.internalize = 0.0
        tt.pYinternalize = 0.0
        tt.gasCur *= 10.0

        secondSurf = TAMode.runTAMinit([100.0], tt, firstSurf)
        reduce = TAMode.surface .* TAMode.total

        @test aboutZero(dot(vec(firstSurf) - vec(secondSurf), reduce))
    end

    @testset "Ensure that system reaches equilibrium." begin
        rr = TAMode.param(params)
        du = zeros(55)
        ddnorm = TAMode.TAMreact(du, TAMode.getAutocrine(rr), rr, 0.0)
        @test norm(du) < 1.0e-4
    end
end


text:
#script for the 2D Linear elasticity problem
#plate with a hole under tension at infinity benchmark
include("../src/Solver2D.jl")
include("../src/SplinePlotting.jl")

#Material properties
Emod = 1e5
nu = 0.3
Cmat = Emod/(1-nu^2)*[1  nu  0; nu  1  0; 0  0  (1-nu)/2]
material = Material2D(Emod, nu, Cmat)

#define the domain geometry
R = 1.
L = 4.
center = [0., 0.]
sang = π/2
eang = π

bottom_crv =  nrbcirc(R, center, sang, eang)
bottom_crv = nrbreverse(bottom_crv)
right_seg = nrbline([0, R], [0, L])
left_seg = nrbline([-R, 0], [-L, 0])
up_seg = nrbline([-L, 0], [-L, L], [0, L])
nrb = nrbcoons(bottom_crv, up_seg, left_seg, right_seg)

#degree elevate and refine
degP = [3, 3]
num_subdiv = [10, 20]
nrb = nrbdegelev(nrb, degP.+1-nrb.order)
_, _, new_knots = kntrefine(nrb.knots, num_subdiv, degP, degP.-1)
nrb = nrbkntins(nrb, new_knots)
# Plots.surface(reuse = false)
# plt = nrbctrlplot(nrb)
# display(plt)

#symmetry Dirichlet B.C., u_y = 0 for y=0 and u_x=0 for x=0
u_bound_dir_symy0(x,y) = [undef, 0]
u_bound_dir_symx0(x,y) = [0, undef]
# traction at infinity
tx = 10

#define the exact solution
function exact_disp(x,y)
    cart2pol(x, y) = hypot(x, y), atan(y, x)
    r, th = cart2pol(x, y)
    ux = (1+nu)/Emod*tx*(1/(1+nu)*r*cos(th)+2*R^2/((1+nu)*r)*cos(th)+R^2/(2*r)*cos(3*th)-R^4/(2*r^3)*cos(3*th))
    uy = (1+nu)/Emod*tx*(-nu/(1+nu)*r*sin(th)-(1-nu)*R^2/((1+nu)*r)*sin(th)+
            R^2/(2*r)*sin(3*th)-R^4/(2*r^3)*sin(3*th))
    return [ux, uy]
end

function exact_stress(x,y)
    cart2pol(x, y) = hypot(x, y), atan(y, x)
    r, th = cart2pol(x, y)
    sigma_rr = tx/2*(1-R^2/r^2)+tx/2*(1-4*R^2/r^2+3*R^4/r^4)*cos(2*th)
    sigma_tt = tx/2*(1+R^2/r^2)-tx/2*(1+3*R^4/r^4)*cos(2*th)
    sigma_rt = -tx/2*(1+2*R^2/r^2-3*R^4/r^4)*sin(2*th)

    A = [cos(th)^2 sin(th)^2 2*sin(th)*cos(th); sin(th)^2 cos(th)^2 -2*sin(th)*cos(th);
        -sin(th)*cos(th) sin(th)*cos(th) cos(th)^2-sin(th)^2]
    stress = A\[sigma_rr;sigma_tt;sigma_rt]
    return stress
end
"""
Define the Neumann boundary condition derived from the exact stresses
"""
function u_bound_neu(x,y)
    tol_eq=1e-10
    ex_stress = exact_stress(x,y)
    if abs(x+L)<tol_eq
        #on the boundary x=-L the outer normal is (-1, 0)
        nx = -1
        ny = 0
    elseif abs(y-L)<tol_eq
        #on the boundary y=L, the outer normal is (0, 1)
        nx = 0
        ny = 1
    else
        error("Point not on the Neumann boundary")
    end
    tx = nx*ex_stress[1]+ny*ex_stress[3]
    ty = nx*ex_stress[3]+ny*ex_stress[2]
    return tx, ty
end

bound_left = Boundary2D("Dirichlet", "Left", u_bound_dir_symy0)
bound_right = Boundary2D("Dirichlet", "Right", u_bound_dir_symx0)
bound_up = Boundary2D("Neumann", "Up", u_bound_neu)
bound_all = [bound_left, bound_right, bound_up]

#generate the NURBS mesh
IGAmesh = genMesh(nrb)

gauss_rule = [genGaussLegendre(degP[1]+1), genGaussLegendre(degP[2]+1)]
stiff = assemble_stiff_elast2D(IGAmesh, Cmat, gauss_rule)
#Assume no body forces
rhs = zeros(2*IGAmesh.numBasis)
bcdof_all, elem_all = classifyBoundary2D(IGAmesh)
lhs, rhs = applyBCnurbsElast(IGAmesh, bound_all, stiff, rhs, bcdof_all, elem_all, gauss_rule)
println("Solving linear system with $(length(rhs)) DOFs")
@time sol0 = lhs\rhs
@time plotSolElast(IGAmesh, Cmat, sol0, "PlateHoleC0")
@time plotSolErrorElast(IGAmesh, Cmat, sol0, exact_disp, exact_stress, "PlateHoleC0Err")
@time relL2Err, relH1Err = compErrorNormElast(IGAmesh, Cmat, sol0, exact_disp, exact_stress, gauss_rule)
println("Relative L2-norm error is $relL2Err")
println("Relative energy-norm error is $relH1Err")


text:
using Test
using BenchmarkTools, LinearAlgebra, StaticArrays
using Plots; pyplot()
using SolidStateDetectors

T = Float64

np = SSD.CartesianPoint{T}(0, 0, 0)
pt = SSD.CartesianPoint{T}(0.64, 0.15, 0.54)

ex = SSD.CartesianVector{T}(1.0, 0.0, 0.0)
ey = SSD.CartesianVector{T}(0.0, 1.0, 0.0)
ez = SSD.CartesianVector{T}(0.0, 0.0, 1.0)

rect = SSD.RectangularCuboid{T}( np, ex, ey, ez )
@test in(np, rect)
@test in(np + ex, rect)
@test in(np + ey, rect)
@test in(np + ez, rect)
@test in(pt, rect)
@test !in(pt + ex, rect)
@test !in(pt + ey, rect)
@test !in(pt + ez, rect)

cyl = SSD.Cylinder{T}( np, ez, T(0.4) )

l1 = SSD.Line( SSD.CartesianPoint{T}(1, 0, 0), SSD.CartesianVector{T}(0, 1, 0) )
l2 = SSD.Line( SSD.CartesianPoint{T}(0, 1, 0), SSD.CartesianVector{T}(1, 0, 0) )
r1 = SSD.Ray( np, ex )
r2 = SSD.Ray( np, ey )
edge = SSD.LineSegment(np, SSD.CartesianPoint{T}(0, 1, 0))
edge1 = SSD.LineSegment(SSD.CartesianPoint{T}(0, 1, 0),  SSD.CartesianPoint{T}(1, 0, 0))
edge2 = SSD.LineSegment(SSD.CartesianPoint{T}(-10, 0, 0),  SSD.CartesianPoint{T}(1, 1, 0))

csg_diff = rect - cyl
@test !in(np, csg_diff)
@test in(np + ey, csg_diff)


text:
#useful functions that do no directly call C code

function getAllVertCoords()
numV = num_entities[1]

vertCoords = zeros(3, numV)  # storage for all vertex coordinates
coords_tmp = zeros(3,1)  # temporary storage for vetex coordinates

for i=1:numV
  apf.getVertCoords(coords_tmp, 3, 1)
  vertCoords[:, i] = coords_tmp
  incrementVertIt()
end

resetVertIt()

return vertCoords
end

function makeSBPMatrix()
# make matrix needed for summation by parts

numEl = num_entities[4]  # number of elements
numPerEl = downward_counts[4, 1]  # nodes per element

coords_tmp = zeros(3, numPerEl)
sbpMatrix = zeros(3, numPerEl, numEl)

for i = 1:numEl
    apf.getElCoords( coords_tmp, 3,  numPerEl)
    sbpMatrix[:, :, i] = coords_tmp
    incrementElIt()
end

resetElIt()
return sbpMatrix
end

function incrementIt(dim::Int)
# increment an iterator of a given dimension
if (dim == 0)
  incrementVertIt()
elseif (dim==1)
  incrementEdgeIt()
elseif (dim == 2)
  incrementElIt()
elseif (dim == 3)
  incrementElIt()
end

end  # end function


text:
push!(LOAD_PATH,"../src/")
using Documenter, Triangle

makedocs(
	format = Documenter.HTML(),
	sitename = "Triangle.jl",
	pages = [
		"Home" => "index.md"
	]
)