text:
# MIT License

# Copyright (c) Microsoft Corporation.

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:

# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.

# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE

"""
    WrapperSurface{T,S<:Surface{T}} <: Surface{T}

A generic surface type which serves as a basis for extension of [`Surface`](@ref)s for custom [`OpticalInterface`](@ref) subclasses.
Essentially just forwards all `Surface` and `ParametricSurface` methods to a field of the `WrapperSurface` named `surface`.
Also provides a generic implementation of [`surfaceintersection`](@ref) which tests for an intersection with the underlying surface and returns either an [`EmptyInterval`](@ref) or a half space (never a closed interval).
"""
abstract type WrapperSurface{T,S<:Surface{T}} <: Surface{T} end
export WrapperSurface

# forward all methods, can be specialised if not applicable for subclass
interface(r::WrapperSurface{T}) where {T<:Real} = r.interface
centroid(r::WrapperSurface{T,S}) where {T<:Real,S<:Surface{T}} = centroid(r.surface::S)
normal(r::WrapperSurface{T,S}) where {T<:Real,S<:Surface{T}} = normal(r.surface::S)
normal(r::WrapperSurface{T}, ::T, ::T) where {T<:Real} = normal(r)
point(r::WrapperSurface{T,S}, u::T, v::T) where {T<:Real,S<:Surface{T}} = point(r.surface::S, u, v)
uv(r::WrapperSurface{T,S}, x::T, y::T, z::T) where {T<:Real,S<:Surface{T}} = uv(r, SVector{3,T}(x, y, z))
uv(r::WrapperSurface{T,S}, p::SVector{3,T}) where {T<:Real,S<:Surface{T}} = uv(r.surface::S, p)
uvrange(::Type{WrapperSurface{T,S}}) where {T<:Real,S<:Surface{T}} = uvrange(S)
uvrange(::WrapperSurface{T,S}) where {T<:Real,S<:Surface{T}} = uvrange(S)
onsurface(a::WrapperSurface{T,S}, x::T, y::T, z::T) where {T<:Real,S<:Surface{T}} = onsurface(a.surface::S, x, y, z)
onsurface(a::WrapperSurface{T,S}, p::SVector{3,T}) where {T<:Real,S<:Surface{T}} = onsurface(a.surface::S, p)
inside(a::WrapperSurface{T,S}, x::T, y::T, z::T) where {T<:Real,S<:Surface{T}} = inside(a.surface::S, x, y, z)
inside(a::WrapperSurface{T,S}, p::SVector{3,T}) where {T<:Real,S<:Surface{T}} = inside(a.surface::S, p)
partials(a::WrapperSurface{T,S}, u::T, v::T) where {T<:Real,S<:Surface{T}} = partials(a.surface::S, u, v)
makemesh(l::WrapperSurface{T,S}, subdivisions::Int = 20) where {T<:Real,S<:Surface{T}} = makemesh(l.surface::S, subdivisions)

function surfaceintersection(l::WrapperSurface{T,S}, r::AbstractRay{T,3}) where {T<:Real,S<:Surface{T}}
    itvl = surfaceintersection(l.surface::S, r)
    if itvl isa EmptyInterval{T}
        return EmptyInterval(T)
    elseif itvl isa DisjointUnion{T}
        throw(ErrorException("Grating/HOE must use a simple (planar) surface for now"))
    else
        intsct = halfspaceintersection(itvl)
        u, v = uv(intsct)
        intsct = Intersection(α(intsct), point(intsct), normal(intsct), u, v, interface(l), flippednormal = flippednormal(intsct))
        if dot(normal(intsct), direction(r)) > zero(T)
            return rayorigininterval(intsct)
        else
            return positivehalfspace(intsct)
        end
    end
end

######################################################################

"""
    ThinGratingSurface{T,S} <: WrapperSurface{T,S}

Surface type for use with [`ThinGratingInterface`](@ref).

```julia
ThinGratingSurface(surface::Surface{T}, interface::ThinGratingInterface{T})
```
"""
struct ThinGratingSurface{T,S} <: WrapperSurface{T,S}
    surface::S
    interface::ThinGratingInterface{T}
    ThinGratingSurface(surface::S, interface::ThinGratingInterface{T}) where {T<:Real,S<:Surface{T}} = new{T,S}(surface, interface)
end
export ThinGratingSurface

interface(r::ThinGratingSurface{T}) where {T<:Real} = r.interface

function processintersection(opticalinterface::ThinGratingInterface{T}, point::SVector{N,T}, normal::SVector{N,T}, incidentray::OpticalRay{T,N}, temperature::T, pressure::T, test::Bool, firstray::Bool = false) where {T<:Real,N}
    # we want the surface to work if hit from either side, so we just reverse the normal if it is hit on the back side
    if dot(direction(incidentray), normal) > zero(T)
        normal = -normal
    end

    mᵢ, mₜ = mᵢandmₜ(outsidematerialid(opticalinterface), insidematerialid(opticalinterface), normal, incidentray)
    nᵢ = one(T)
    nₜ = one(T)
    α = zero(T)
    if !isair(mᵢ)
        mat = glassforid(mᵢ)::OpticSim.GlassCat.Glass
        nᵢ = index(mat, λ, temperature = temperature, pressure = pressure)::T
        α = absorption(mat, λ, temperature = temperature, pressure = pressure)::T
    end
    if !isair(mₜ)
        nₜ = index(glassforid(mₜ)::OpticSim.GlassCat.Glass, λ, temperature = temperature, pressure = pressure)::T
    end

    incident_pow = power(incidentray)
    # compute updated power based on absorption coefficient of material using Beer's law
    internal_trans = one(T)
    if α > zero(T)
        internal_trans = exp(-α * geometricpathlength)
    end

    # optical distance from ray origin to point of intersection in mm. Compensate for the fact that the ray has been slightly shortened.
    geometricpathlength = norm(point - origin(incidentray)) + (firstray ? zero(T) : T(RAY_OFFSET))
    thisraypathlength = nᵢ * geometricpathlength
    raypathlength = pathlength(incidentray) + thisraypathlength

    λ = wavelength(incidentray)
    r = !test * rand()
    if opticalinterface.period < λ / nₜ
        # TODO!! what to do in this case?
        raypower = transmission(opticalinterface, 0) * internal_trans
        if r >= raypower
            return nothing
        end
        raydirection = direction(incidentray)
    else
        # do the diffraction
        order = rand((opticalinterface.minorder):(opticalinterface.maxorder)) # TODO not the most sensible way to sample this...
        refl = reflectance(opticalinterface, order) * internal_trans
        trans = transmission(opticalinterface, order) * internal_trans
        if r < refl
            raypower = refl * incident_pow
            normal = -normal
        elseif r < trans
            raypower = trans * incident_pow
        else
            return nothing
        end
        incident_mag = 2π * nᵢ / λ
        grating_dir = order * opticalinterface.vector * (2π / opticalinterface.period)
        u = cross(normal, ((direction(incidentray) * incident_mag) + grating_dir))
        output_mag = 2π * nₜ / λ
        q = output_mag^2 - sum(u .^ 2)
        if q < zero(T)
            @warn "Order $order not valid for this configuration, skipping this ray" maxlog = 1
            return nothing
        else
            raydirection = normalize(normal * sqrt(q) - cross(normal, u))
        end
    end
    return raydirection, raypower, raypathlength
end

######################################################################

"""
    HologramSurface{T,S} <: WrapperSurface{T,S}

Surface type for use with [`HologramInterface`](@ref).

```julia
HologramSurface(surface::Surface{T}, interface::HologramInterface{T})
```
"""
struct HologramSurface{T,S} <: WrapperSurface{T,S}
    surface::S
    interface::HologramInterface{T}
    HologramSurface(surface::S, interface::HologramInterface{T}) where {T<:Real,S<:Surface{T}} = new{T,S}(surface, interface)
end
export HologramSurface

"""
    MultiHologramSurface{T,S} <: WrapperSurface{T,S}

Surface type for use with [`MultiHologramInterface`](@ref).

```julia
MultiHologramSurface(surface::Surface{T}, interface::MultiHologramInterface{T})
```
"""
struct MultiHologramSurface{T,S} <: WrapperSurface{T,S}
    surface::S
    interface::MultiHologramInterface{T}
    MultiHologramSurface(surface::S, interface::MultiHologramInterface{T}) where {T<:Real,S<:Surface{T}} = new{T,S}(surface, interface)
end
export MultiHologramSurface

function processintersection(opticalinterface::HologramInterface{T}, point::SVector{N,T}, normal::SVector{N,T}, incidentray::OpticalRay{T,N}, temperature::T, pressure::T, test::Bool, firstray::Bool = false) where {T<:Real,N}
    hitback = dot(direction(incidentray), normal) > zero(T)
    # we want the surface to work if hit from either side, so we just reverse the normal and interfaces if it is hit on the back side
    if hitback
        frontinterface = FresnelInterface{T}(opticalinterface.substratematerial, opticalinterface.aftermaterial)
        backinterface = FresnelInterface{T}(opticalinterface.substratematerial, opticalinterface.beforematerial)
        normal = -normal
    else
        frontinterface = FresnelInterface{T}(opticalinterface.substratematerial, opticalinterface.beforematerial)
        backinterface = FresnelInterface{T}(opticalinterface.substratematerial, opticalinterface.aftermaterial)
    end

    # refract the incoming ray going from beforematerial to substrate material, can have Fresnel/TI reflections
    tmp = processintersection(frontinterface, point, normal, incidentray, temperature, pressure, test, firstray)
    if tmp === nothing
        return nothing
    end
    input_ray, input_power, input_opl = tmp
    if dot(input_ray, normal) > zero(T)
        # just reflected of the surface through fresenel reflection or TIR, so return without diffraction
        return input_ray, input_power, input_opl
    end

    λ = wavelength(incidentray)
    mat = glassforid(opticalinterface.substratematerial)::OpticSim.GlassCat.Glass
    # get the index of the playback ray in the substrate
    nₛ = index(mat, λ, temperature = temperature, pressure = pressure)::T
    # get the index of the recording ray in the substrate
    nₛrec = index(mat, opticalinterface.recordingλ, temperature = temperature, pressure = pressure)::T

    # get the index of the recording ray in the material of the signal beam
    if !isair(opticalinterface.signalrecordingmaterial)
        signalbeammaterial = glassforid(opticalinterface.signalrecordingmaterial)::OpticSim.GlassCat.Glass
        nsig = index(signalbeammaterial, opticalinterface.recordingλ, temperature = temperature, pressure = pressure)::T
    else
        nsig = one(T)
    end
    # get the signal magnitue in the substrate
    mag_Ksig = 2π * nₛ / opticalinterface.recordingλ
    # and direction
    if opticalinterface.signalbeamstate === CollimatedBeam
        Ksig_raw = opticalinterface.signalpointordir
    elseif opticalinterface.signalbeamstate === ConvergingBeam
        Ksig_raw = normalize(opticalinterface.signalpointordir - point)
    elseif opticalinterface.signalbeamstate === DivergingBeam
        Ksig_raw = normalize(point - opticalinterface.signalpointordir)
    else
        throw(ErrorException("Invalide beam state"))
    end
    # refract it for the substrate interface
    sigonside = dot(Ksig_raw, normal) < zero(T)
    refr = refractedray(nsig, nₛrec, sigonside ? -normal : normal, Ksig_raw) # TODO not sure normal flip is needed
    if refr === nothing
        return nothing
    end
    Ksig_refracted = mag_Ksig * refr

    # get the index of the recording ray in the material of the reference beam
    if !isair(opticalinterface.referencerecordingmaterial)
        referencebeammaterial = glassforid(opticalinterface.referencerecordingmaterial)::OpticSim.GlassCat.Glass
        nref = index(referencebeammaterial, opticalinterface.recordingλ, temperature = temperature, pressure = pressure)::T
    else
        nref = one(T)
    end
    # get the reference magnitue in the substrate
    mag_Kref = 2π * nₛ / opticalinterface.recordingλ
    if opticalinterface.referencebeamstate === CollimatedBeam
        Kref_raw = opticalinterface.referencepointordir
    elseif opticalinterface.referencebeamstate === ConvergingBeam
        Kref_raw = normalize(opticalinterface.referencepointordir - point)
    elseif opticalinterface.referencebeamstate === DivergingBeam
        Kref_raw = normalize(point - opticalinterface.referencepointordir)
    else
        throw(ErrorException("Invalide beam state"))
    end
    # refract it for the substrate interface
    refonside = dot(Kref_raw, normal) < zero(T)
    refr = refractedray(nref, nₛrec, refonside ? -normal : normal, Kref_raw) # TODO not sure normal flip is needed
    if refr === nothing
        return nothing
    end
    Kref_refracted = mag_Kref * refr

    # check whether we should be doing reflection or transmission
    isreflection = sigonside != refonside
    if isreflection
        # if the vectors are opposite in relation to the normal then the HOE should work in reflection
        dnormal = normal
    else
        # otherwise it is transmission
        dnormal = -normal
    end

    # calulcate the grating vector
    grating_vec = (Ksig_refracted - Kref_refracted)
    # choose random order for which Kogelnik is valid
    if !opticalinterface.include0order
        order = 1
    else
        order = rand(0:1)
    end

    # calculate the magnitudes of input and output
    incident_mag = 2π * nₛ / λ
    incident_ray = incident_mag * input_ray
    # do the diffraction
    u = cross(dnormal, (incident_ray + order * grating_vec))
    output_mag = 2π * nₛ / λ
    q = output_mag^2 - sum(u .^ 2)
    if q < zero(T)
        @warn "Order $order not valid for this configuration, skipping this ray" maxlog = 1
        return nothing
    else
        diffracted_out = dnormal * sqrt(q) - cross(dnormal, u)
    end
    diffracted_out = normalize(diffracted_out)

    # calculate the efficiency in this case
    d = opticalinterface.thickness # in microns as λ in microns
    Δn = opticalinterface.RImodulation
    σ = incident_ray + grating_vec
    cₒ = dot(σ, -normal) / incident_mag
    cᵣ = dot(incident_ray, -normal) / incident_mag
    υ = (norm(incident_ray)^2 - norm(σ)^2) / (2 * incident_mag) # Kick from Kogelnik eq. 17
    ξ = (υ * d) / (2 * cₒ) # Kogelnik eq. 42
    ν = (π * Δn * d) / (λ * sqrt(complex(cₒ * cᵣ))) # Kogelnik eq. 42
    refl = dot(diffracted_out, normal) > zero(T)
    if refl # reflected
        # modulate the power according to Kogelnik - Lossless Dielectric Grating (Reflection)
        ν = 1im * ν # Kogelnik eq. 55
        ξ = -ξ # Kogelnik eq. 55
        η = 1 / (1 + (1 - (ξ / ν)^2) / (sinh(sqrt(ν^2 - ξ^2))^2)) # Kogelnik eq. 57, Kick eq. 6
    else # refracted
        # modulate the power according to Kogelnik - Lossless Dielectric Grating (Transmission)
        η = sin(sqrt(ν^2 + ξ^2))^2 / (1 + (ξ / ν)^2) # Kogelnik eq. 43, Kick eq. 4
    end

    # kill the ray based on efficiency for monte carlo
    if !(imag(η) == zero(T) && real(η) > zero(T))
        @warn "Invalid η: $η"
        return nothing
    else
        if order === 0
            k = (1 - real(η))
            if !test * rand() >= k
                return nothing
            end
            diffracted_power = input_power * k
        else
            k = real(η)
            if !test * rand() >= k
                return nothing
            end
            diffracted_power = input_power * k
        end
    end
    outray = OpticalRay(point, diffracted_out, diffracted_power, λ)
    # refract the output ray for the substrate/aftermaterial interface
    # for refraction back interface is the other way around so reverse normal
    tmp = processintersection(refl ? frontinterface : backinterface, point, refl ? normal : -normal, outray, temperature, pressure, test, true)
    if tmp === nothing
        return nothing
    end
    outpur_dir, output_pow, _ = tmp
    # TODO nothing happening to OPL currently
    return outpur_dir, output_pow, input_opl
end

function processintersection(opticalinterface::MultiHologramInterface{T}, point::SVector{N,T}, normal::SVector{N,T}, incidentray::OpticalRay{T,N}, temperature::T, pressure::T, test::Bool, firstray::Bool = false) where {T<:Real,N}
    # minη = typemax(T)
    # r = rand(T)
    # Ση = zero(T)
    # for i in 1:(opticalinterface.numinterfaces)
    #     int = interface(opticalinterface, i)::HologramInterface{T}
    #     tmp = processintersection(int, point, normal, incidentray, temperature, pressure, scatter, firstray)
    #     if tmp !== nothing
    #         dir, pow, opl = tmp
    #         η = pow / power(incidentray)
    #         Ση += η / opticalinterface.numinterfaces
    #         # @assert Ση <= 1.0
    #         if r < Ση
    #             # don't modulate the power because we sample the efficiencies proportionally
    #             return dir, power(incidentray), opl
    #         end
    #     end
    # end
    # return nothing
    # when test = true behavior is totally invalid so that the test is deteministic...
    i = test ? 1 : rand(1:(opticalinterface.numinterfaces)) # TODO definitely not the best way to sample this...
    int = interface(opticalinterface, i)::HologramInterface{T}
    return processintersection(int, point, normal, incidentray, temperature, pressure, test, firstray)
end


text:
import CompEcon
using PyPlot
using PyCall

reload("CompEcon")
## DEMAPP01 Approximating functions on R

# This m-file illustrates how to use CompEcon Toolbox routines to construct
# and operate with an approximant for a function defined on an interval of
# the real line.

# In particular, we construct an approximant for f(x)=exp(-x) on the
# interval [-1,1].  The function used in this illustration posseses a
# closed-form, which will allow us to measure approximation error precisely.
# Of course, in practical applications, the function to be approximated
# will not possess a known closed-form.

# In order to carry out the exercise, one must first code the function to
# be approximated at arbitrary points.  The required code is presented at
# the end of this m-file (see below). Let's begin:

# Preliminary tasks
function main()
        f(x) = exp(-x)

    # Set the endpoints of approximation interval:
    a =  -1                            # left endpoint
    b =   1                            # right endpoint

    # Choose an approximation scheme. In this case, let us use an order 10
    # Chebychev approximation scheme:
    n = 10                             # order of approximation
    basis = CompEcon.fundefn(:cheb, n, a, b)      # define basis

    # Compute the basis coefficients c.  There are various way to do this:
    # One may use funfitf:
    c = CompEcon.funfitf(basis, f)

    # ... or one may compute the standard approximation nodes x and
    # corresponding function values y and use funfitxy:
    x = CompEcon.funnode(basis)[1]
    y = f(x)
    c = CompEcon.funfitxy(basis, x, y)

    # ... or one compute the standard approximation nodes x, corresponding
    # function values y, and the interpolation matrix phi, and solve the
    # interpolation equation directly using the backslash operator:
    x = CompEcon.funnode(basis)[1]
    y = f(x)
    phi = CompEcon.funbase(basis)
    c = phi\y

    # Having computed the basis coefficients, one may now evaluate the
    # approximant at any point x using funeval:
    x = [0.0]
    y = CompEcon.funeval(c, basis, x)[1]
    println("The approximate value of exp(-x) at x=0 is $y")
    println("The ''exact'' value of exp(-x) at x=0 is 1")

    # ... one may also evaluate the approximant's first and second derivatives
    # at x:
    d1 = CompEcon.funeval(c, basis, x, 1)[1][1]
    d2 = CompEcon.funeval(c, basis, x, 2)[1][1]
    println("The approximate first derivative of exp(-x) at x=0 is $d1")
    println("The ''exact'' first derivative of exp(-x) at x=0 is -1")
    println("The approximate second derivative of exp(-x) at x=0 is $d2")
    println("The ''exact'' second derivative of exp(-x) at x=0 is 1")

    # ... and one may even evaluate the approximant's definite integral between
    # the left endpoint a and x:
    int = CompEcon.funeval(c, basis, x, -1)[1]
    println("The approximate integral of exp(-x) between x=-1 and x=0 is $int")
    println("The `exact` integral of exp(-x) between x ∈[-1, 0] is $(exp(1)-1)")

    # One may evaluate the accuracy of the Chebychev polynomial approximant by
    # computing the approximation error on a highly refined grid of points:
    ngrid = 5001                       # number of grid nodes
    xgrid = range(a, stop=b, length=ngrid)        # generate refined grid for plotting

    function plot_approx(f::Function, basistype::Symbol, n, a, b, xgrid, k=0)
        basis = basistype == :spli ? CompEcon.fundefn(:spli, n, a, b, k) :
                basistype == :cheb ? CompEcon.fundefn(:cheb, n, a, b) :
                error("Not doing this")

        # compute approximation
        c = CompEcon.funfitf(basis, f)
        yapp = CompEcon.funeval(c, basis, xgrid)[1]

        # compute exact
        yact = f(xgrid)

        # plot error
        fig, ax = plt.subplots()
        ax[:plot](xgrid,yapp-yact)
        ax[:plot](xgrid, zeros(ngrid), "k--", linewidth=2)
        ax[:set_xlabel]("x")
        ax[:set_ylabel]("Error")
        nm = basistype == :spli ? "Degree $k spline" : "Degree $n Chebyshev"
        ax[:set_title]("$nm Approximation Error for exp(-x)")
    end

    plot_approx(f, :cheb, n, a, b, xgrid)
    plot_approx(f, :spli, 21, a, b, xgrid, 3)
    # plot_approx(f, :cheb, 31, a, b, xgrid, 1)
end

main()

ext: jl


text:
abstract type Dimension end
abstract type _1D <: Dimension end
abstract type _2D <: Dimension end
abstract type _3D <: Dimension end

function checkDim(D::Type{<:Dimension})
    if D === _1D
        return 1
    elseif D === _2D
        return 2
    elseif D === _3D
        return 3
    end
    return 0
end

mutable struct FVMPDEGrid
    X
    Y
    Z
    nx
    ny
    nz
    dX
    dY
    dZ
    dX_C
    dY_C
    dZ_C
    U
    Uboundries
    dimension::Type{<:Dimension}
    indices
    function FVMPDEGrid(X,Y,Z,nx,ny,nz,dX,dY,dZ,dX_C,dY_C,dZ_C,U,Uboundries,dimension,indices)
        return new(X,Y,Z,nx,ny,nz,dX,dY,dZ,dX_C,dY_C,dZ_C,U,Uboundries,dimension,indices)
    end
end

function FVMPDEGrid(_X::AbstractVector, _Y::AbstractVector, _Z::AbstractVector)
    X = _X
    Y = _Y
    Z = _Z
    Xb = [_X[1];(_X[1:end-1].+_X[2:end])./2;_X[end]]
    Yb = [_Y[1];(_Y[1:end-1].+_Y[2:end])./2;_Y[end]]
    Zb = [_Z[1];(_Z[1:end-1].+_Z[2:end])./2;_Z[end]]

    nx = length(X)
    ny = length(Y)
    nz = length(Z)
    dX = _X[2:end] .- _X[1:end-1]
    dY = _Y[2:end] .- _Y[1:end-1]
    dZ = _Z[2:end] .- _Z[1:end-1]
    dX_C = Xb[2:end] .- Xb[1:end-1]
    dY_C = Yb[2:end] .- Yb[1:end-1]
    dZ_C = Zb[2:end] .- Zb[1:end-1]
    #Uboundries = zeros(length(_X),length(_Y))
    dimension = _3D
    indices = CartesianIndices((Base.OneTo(nx),Base.OneTo(ny),Base.OneTo(nz)))
    return FVMPDEGrid(X,Y,Z,nx,ny,nz,dX,dY,dZ,dX_C,dY_C,dZ_C,nothing,nothing,dimension,indices)
end

function FVMPDEGrid(_X::AbstractVector, _Y::AbstractVector)
    X = _X
    Y = _Y
    Xb = [_X[1];(_X[1:end-1].+_X[2:end])./2;_X[end]]
    Yb = [_Y[1];(_Y[1:end-1].+_Y[2:end])./2;_Y[end]]
    nx = length(X)
    ny = length(Y)
    dX = _X[2:end] .- _X[1:end-1]
    dY = _Y[2:end] .- _Y[1:end-1]
    dX_C = Xb[2:end] .- Xb[1:end-1]
    dY_C = Yb[2:end] .- Yb[1:end-1]
    #Uboundries = zeros(length(_X),length(_Y))
    dimension = _2D
    indices = CartesianIndices((Base.OneTo(nx),Base.OneTo(ny)))
    return FVMPDEGrid(X,Y,[],nx,ny,1,dX,dY,nothing,dX_C,dY_C,nothing,nothing,nothing,dimension,indices)
end

function FVMPDEGrid(_X::AbstractVector)
    X = _X
    Xb = [_X[1];(_X[1:end-1].+_X[2:end])./2;_X[end]]
    nx = length(X)
    dX = _X[2:end] .- _X[1:end-1]
    dX_C = Xb[2:end] .- Xb[1:end-1]
    #Uboundries = zeros(length(_X))
    dimension = _1D
    indices = CartesianIndices((Base.OneTo(nx)))
    return FVMPDEGrid(X,[],[],nx,1,1,dX,nothing,nothing,dX_C,nothing,nothing,nothing,nothing,dimension,indices)
end


text:
using Revise
using CUDA
using StaticArrays
using StatsBase: sample
using BenchmarkTools

function hist_cpu!(hist, δ, idx)
    Threads.@threads for j in 1:size(idx, 2)
        for i in 1:size(idx, 1)
            @inbounds hist[idx[i,j], j] += δ[i,j]
        end
    end
    return
end

# base kernel
function kernel1!(h, x, id)
    i = threadIdx().x + (blockIdx().x - 1) * blockDim().x
    j = threadIdx().y + (blockIdx().y - 1) * blockDim().y
    if i <= size(id, 1) && j <= size(id, 2)
        @inbounds k = Base._to_linear_index(h, id[i,j], j)
        # @inbounds k = id[i,j] + 32 * (j-1)
        @inbounds CUDA.atomic_add!(pointer(h, k), x[i,j])
    end
    return
end

# base approach - block built along the cols first, the rows (limit collisions)
function hist_gpu1!(h::AbstractMatrix{T}, x::AbstractMatrix{T}, id::AbstractMatrix{S}; MAX_THREADS=256) where {T,S}
    # thread_i = min(MAX_THREADS, size(id, 1))
    # thread_j = min(MAX_THREADS ÷ thread_i, size(id, 2))
    thread_j = min(MAX_THREADS, size(id, 2))
    thread_i = min(MAX_THREADS ÷ thread_j, size(id, 1))
    threads = (thread_i, thread_j)
    blocks = ceil.(Int, (size(id, 1), size(id, 2)) .÷ threads)
    # println("threads:", threads)
    # println("blocks:", blocks)
    @cuda blocks = blocks threads = threads kernel1!(h, x, id)
    return
end

nbins = 32
ncol = 100
items = Int32(1e6)
hist = zeros(Float32, nbins, ncol)
δ = rand(Float32, items, ncol)
idx = UInt8.(rand(1:nbins, items, ncol))
# idx = Int64.(rand(1:nbins, items, ncol))

hist_gpu = CuArray(hist)
δ_gpu = CuArray(δ)
idx_gpu = CuArray(idx)

hist .- Array(hist_gpu)
sum(hist) - sum(Array(hist_gpu))

@CUDA.time hist_gpu1!(hist_gpu, δ_gpu, idx_gpu, MAX_THREADS=1024)
@time hist_cpu!(hist, δ, idx)
@btime hist_cpu!($hist, $δ, $idx)
@btime CUDA.@sync hist_gpu1!($hist_gpu, $δ_gpu, $idx_gpu, MAX_THREADS=1024)
# test on view
@CUDA.time hist_gpu1!(hist_gpu, view(δ_gpu, 1:items ÷ 2, 1:ncol ÷ 2), view(idx_gpu, 1:items ÷ 2, 1:ncol ÷ 2), MAX_THREADS=1024)

size(δ_gpu)
size(view(δ_gpu, 1:items ÷ 2, 1:ncol ÷ 2))

##############################################################
## Build histogram from a subsample idx
# base kernel
function kernel2!(h::CuDeviceMatrix{T}, x::CuDeviceMatrix{T}, id, 𝑖, 𝑗) where {T <: AbstractFloat}
    i = threadIdx().x + (blockIdx().x - 1) * blockDim().x
    j = threadIdx().y + (blockIdx().y - 1) * blockDim().y
    if i <= length(𝑖) && j <= length(𝑗)
        @inbounds k = Base._to_linear_index(h, id[𝑖[i], 𝑗[j]], 𝑗[j])
        @inbounds CUDA.atomic_add!(pointer(h, k), x[𝑖[i], 𝑗[j]])
    end
    return
end

# base approach - block built along the cols first, the rows (limit collisions)
function hist_gpu2!(h::CuMatrix{T}, x::CuMatrix{T}, id, 𝑖, 𝑗; MAX_THREADS=1024) where {T <: AbstractFloat}
    thread_j = min(MAX_THREADS, length(𝑗))
    thread_i = min(MAX_THREADS ÷ thread_j, length(𝑖))
    threads = (thread_i, thread_j)
    blocks = ceil.(Int, (length(𝑖), length(𝑗)) ./ threads)
    # println("threads:", threads)
    # println("blocks:", blocks)
    CUDA.@sync begin
        @cuda blocks = blocks threads = threads kernel2!(h, x, id, 𝑖, 𝑗)
    end
    return
end

hist = zeros(Float32, nbins, ncol)
δ = rand(Float32, items, ncol)
idx = rand(1:nbins, items, ncol)
𝑖 = sample(1:items, items ÷ 2, replace=false, ordered=true)
𝑗 = sample(1:ncol, ncol ÷ 2, replace=false, ordered=true)
hist_gpu = CuArray(hist)
δ_gpu = CuArray(δ)
idx_gpu = CuArray(idx)
𝑖_gpu = CuArray(𝑖)
𝑗_gpu = CuArray(𝑗)

@CUDA.time hist_gpu2!(hist_gpu, δ_gpu, idx_gpu, 𝑖_gpu, 𝑗_gpu, MAX_THREADS=1024)
@btime hist_gpu2!($hist_gpu, $δ_gpu, $idx_gpu, 𝑖_gpu, 𝑗_gpu, MAX_THREADS=1024)

#############################################
# test for SVector - basic test - success!
function kernel!(x, y)
    i = threadIdx().x + (blockIdx().x - 1) * blockDim().x
    if i <= length(x)
        # @inbounds x[i] += y[i]
        k = Base._to_linear_index(x, i)
        CUDA.atomic_add!(pointer(x, k), y[i])
    end
    return
end

# base approach - block built along the cols first, the rows (limit collisions)
function hist_gpu!(x, y; MAX_THREADS=1024)
    thread_i = min(MAX_THREADS, length(x))
    threads = (thread_i)
    blocks = ceil.(Int, length(x) .÷ threads)
    CUDA.@sync begin
        @cuda blocks = blocks threads = threads kernel!(x, y)
    end
    return
end

x = rand(SVector{2,Float32}, Int(1e7))
y = rand(SVector{2,Float32}, Int(1e7))
x = rand(Float32, Int(1e7))
y = rand(Float32, Int(1e7))

x_gpu = CuArray(x)
y_gpu = CuArray(y)

@CuArrays.time hist_gpu!(x_gpu, y_gpu)
@btime hist_gpu!($x_gpu, $y_gpu)

#############################################
# test for SVector - real test
function kernelS2!(h, x, id)
    i = threadIdx().x + (blockIdx().x - 1) * blockDim().x
    j = threadIdx().y + (blockIdx().y - 1) * blockDim().y
    if i <= size(id, 1) && j <= size(id, 2)
        @inbounds k = Base._to_linear_index(h, id[i,j], j)
        # @inbounds k = id[i,j] + 32 * (j-1)
        # @inbounds CUDAnative.atomic_add!(pointer(h, k), x[i,j])
        # h[id[i,j],j] += x[i,j]
    end
    return
end

# base approach - block built along the cols first, the rows (limit collisions)
function hist_gpuS2!(h, x, id; MAX_THREADS=256) where {T}
    thread_j = min(MAX_THREADS, size(id, 2))
    thread_i = min(MAX_THREADS ÷ thread_j, size(id, 1))
    threads = (thread_i, thread_j)
    blocks = ceil.(Int, (size(id, 1), size(id, 2)) .÷ threads)
    println("threads:", threads)
    println("blocks:", blocks)
    CUDA.@sync begin
        @cuda blocks = blocks threads = threads kernelS2!(h, x, id)
    end
    return
end

hist = zeros(SVector{2,Float32}, nbins, ncol)
δ = rand(SVector{2,Float32}, items, ncol)
idx = rand(1:nbins, items, ncol)
hist_gpu = CuArray(hist)
δ_gpu = CuArray(δ)
idx_gpu = CuArray(idx)

@CuArrays.time hist_gpuS2!(hist_gpu, δ_gpu, idx_gpu)
@btime hist_gpuS2!($hist_gpu, $δ_gpu, $idx_gpu)

##############################################################
## Build histogram from a subsample idx
# accumulate all gradient single pass
function kernel3!(h, x, id, 𝑖, 𝑗)
    i = threadIdx().x + (blockIdx().x - 1) * blockDim().x
    j = threadIdx().y + (blockIdx().y - 1) * blockDim().y
    if i <= length(𝑖) && j <= length(𝑗)
        @inbounds k = Base._to_linear_index(h, id[𝑖[i], 𝑗[j]], 𝑗[j], 1)
        @inbounds CUDAnative.atomic_add!(pointer(h, k), x[𝑖[i], 𝑗[j], 1])
        @inbounds k = Base._to_linear_index(h, id[𝑖[i], 𝑗[j]], 𝑗[j], 2)
        @inbounds CUDAnative.atomic_add!(pointer(h, k), x[𝑖[i], 𝑗[j], 2])
    end
    return
end

# base approach - block built along the cols first, the rows (limit collisions)
function hist_gpu3!(h, x, id, 𝑖, 𝑗; MAX_THREADS=1024)
    thread_j = min(MAX_THREADS, length(𝑗))
    thread_i = min(MAX_THREADS ÷ thread_j, length(𝑖))
    threads = (thread_i, thread_j)
    blocks = ceil.(Int, (length(𝑖), length(𝑗)) ./ threads)
    # println("threads:", threads)
    # println("blocks:", blocks)
    CuArrays.@sync begin
        @cuda blocks = blocks threads = threads kernel3!(h, x, id, 𝑖, 𝑗)
    end
    return
end

hist = zeros(Float32, nbins, ncol, 2)
δ = rand(Float32, items, ncol, 2)
idx = rand(1:nbins, items, ncol)
𝑖 = sample(1:items, items ÷ 2, replace=false, ordered=true)
𝑗 = sample(1:ncol, ncol ÷ 2, replace=false, ordered=true)
hist_gpu = CuArray(hist)
δ_gpu = CuArray(δ)
idx_gpu = CuArray(idx)
𝑖_gpu = CuArray(𝑖)
𝑗_gpu = CuArray(𝑗)

@CuArrays.time hist_gpu3!(hist_gpu, δ_gpu, idx_gpu, 𝑖_gpu, 𝑗_gpu, MAX_THREADS=1024)
@btime hist_gpu3!($hist_gpu, $δ_gpu, $idx_gpu, 𝑖_gpu, 𝑗_gpu, MAX_THREADS=1024)

# accumulate in shared memory histograms
function kernel2!(h::CuDeviceMatrix{T}, x::CuDeviceMatrix{T}, id, nbins) where {T <: AbstractFloat}
    tid = threadIdx().x
    i = threadIdx().x + (blockIdx().x - 1) * blockDim().x
    j = threadIdx().y + (blockIdx().y - 1) * blockDim().y

    # shared memory on block of length nbins
    # To Do: nbins cannot be passed as argument - dynamic shared memory generate kernel through macro
    # shared = CUDAnative.@cuStaticSharedMem(T, 32)
    # fill!(shared, 0)
    # sync_threads()

    # accumulate in per block histogram
    # Why is the atomic add on shared much longer than atomic on h in global mem?
    if i <= size(id, 1) && j <= size(id, 2)
        # should be the legit way to go - 70ms
        # @inbounds CUDAnative.atomic_add!(pointer(shared, id[i,j]), x[i,j])

        # unsafe (collisions) on shared mem: 3.0ms
        # @inbounds shared[id[i,j]] = x[i,j]

        # unsafe (collisions) add on global memory - 3.6ms
        # @inbounds h[id[i,j],j] += x[i,j]

        # atomic add on global hist - 3.2ms
        @inbounds k = id[i,j] + nbins * (j - 1)
        @inbounds CUDA.atomic_add!(pointer(h, k), x[i,j])
    end
    # sync_threads()

    # if blockIdx().x == 1
    #     if tid <= nbins
    #         CUDA.atomic_add!(pointer(h,tid), shared[tid])
    #     end
    # end
    return
end

# shared memory -
function hist_gpu2!(h::CuMatrix{T}, x::CuMatrix{T}, id::CuMatrix{Int}, nbins; MAX_THREADS=256) where {T <: AbstractFloat}
    # thread_i = min(MAX_THREADS, size(id, 1))
    # thread_j = min(MAX_THREADS ÷ thread_i, size(id, 2))
    thread_j = min(MAX_THREADS, size(id, 2))
    thread_i = min(MAX_THREADS ÷ thread_j, size(id, 1))
    threads = (thread_i, thread_j)
    blocks = ceil.(Int, (size(id, 1), size(id, 2)) ./ threads)
    CUDA.@sync begin
        @cuda blocks = blocks threads = threads kernel2!(h, x, id, nbins)
    end
    return h
end

@CuArrays.time hist_gpu2!(hist_gpu, δ_gpu, idx_gpu, 32, MAX_THREADS=1024)
@btime hist_gpu2!($hist_gpu, $δ_gpu, $idx_gpu, 32, MAX_THREADS=1024)
@device_code_warntype hist_gpu2!(hist_gpu, δ_gpu, idx_gpu, 32, MAX_THREADS=1024)

######################################
# Appoach 1
######################################
# GPU - apply along the features axis
function kernel!(h::CuDeviceMatrix{T}, x::CuDeviceVector{T}, id, 𝑖, 𝑗) where {T <: AbstractFloat}
    i = threadIdx().x + (blockIdx().x - 1) * blockDim().x
    j = threadIdx().y + (blockIdx().y - 1) * blockDim().y
    if i <= length(𝑖) && j <= length(𝑗)
        @inbounds k = Base._to_linear_index(h, id[𝑖[i], 𝑗[j]], 𝑗[j])
        @inbounds CUDA.atomic_add!(pointer(h, k), x[𝑖[i]])
    end
    return
end

# base approach - block built along the cols first, the rows (limit collisions)
function hist_gpu!(h::CuMatrix{T}, x::CuVector{T}, id, 𝑖, 𝑗; MAX_THREADS=1024) where {T <: AbstractFloat}
    thread_j = min(MAX_THREADS, length(𝑗))
    thread_i = min(MAX_THREADS ÷ thread_j, length(𝑖))
    threads = (thread_i, thread_j)
    blocks = ceil.(Int, (length(𝑖), length(𝑗)) ./ threads)
    @cuda blocks = blocks threads = threads kernel!(h, x, id, 𝑖, 𝑗)
    return
end

hist = zeros(Float32, nbins, ncol)
δ = rand(Float32, items)
idx = rand(1:nbins, items, ncol)
𝑖 = sample(1:items, items ÷ 2, replace=false, ordered=true)
𝑗 = sample(1:ncol, ncol ÷ 2, replace=false, ordered=true)
hist_gpu = CuArray(hist)
δ_gpu = CuArray(δ)
idx_gpu = CuArray(idx)
𝑖_gpu = CuArray(𝑖)
𝑗_gpu = CuArray(𝑗)

@CUDA.time hist_gpu!(hist_gpu, δ_gpu, idx_gpu, 𝑖_gpu, 𝑗_gpu, MAX_THREADS=1024)
@btime hist_gpu!($hist_gpu, $δ_gpu, $idx_gpu, 𝑖_gpu, 𝑗_gpu, MAX_THREADS=1024)

######################################
# Appoach 2 - Loop for assigning command grad to appropriate bin per column
# Idea: exploit the fact that there's a single grad per row: take that grad and add it to each column bin
######################################
# GPU - apply along the features axis
function kernel!(h::CuDeviceMatrix{T}, x::CuDeviceVector{T}, id, 𝑖, 𝑗) where {T <: AbstractFloat}
    i = threadIdx().x + (blockIdx().x - 1) * blockDim().x
    j = threadIdx().y + (blockIdx().y - 1) * blockDim().y
    if i <= length(𝑖) && j <= length(𝑗)
        @inbounds k = Base._to_linear_index(h, id[𝑖[i], 𝑗[j]], 𝑗[j])
        @inbounds CUDAnative.atomic_add!(pointer(h, k), x[𝑖[i]])
    end
    return
end

# base approach - block built along the cols first, the rows (limit collisions)
function hist_gpu!(h::CuMatrix{T}, x::CuVector{T}, id::CuMatrix{UInt8}, 𝑖, 𝑗; MAX_THREADS=1024) where {T <: AbstractFloat}
    thread_j = min(MAX_THREADS, length(𝑗))
    thread_i = min(MAX_THREADS ÷ thread_j, length(𝑖))
    threads = (thread_i, thread_j)
    blocks = ceil.(Int, (length(𝑖), length(𝑗)) ./ threads)
    @cuda blocks = blocks threads = threads kernel!(h, x, id, 𝑖, 𝑗)
    return
end

using CUDA
using BenchmarkTools
N1 = Int(2^12)
x1 = rand(Float32, N1);
x2 = rand(Float32, N1);
x1g = CuArray(x1);
x2g = CuArray(x2);

function dot_atomic!(x, y, z)
    idx = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    if idx <= length(x)
        CUDA.atomic_add!(pointer(z, 1), x[idx] * y[idx])
    end
    return nothing
end

function bench_dot_atomic!(x, y, z, threads)
    numblocks = ceil(Int, N1 / threads)
    @cuda threads = threads blocks = numblocks dot_atomic!(x, y, z)
end

threads = 512
numblocks = ceil(Int, N1 / threads)
z0 = CUDA.zeros(1)
# @cuda threads=gthreads blocks=numblocks dot_atomic!(x1g, x2g, z0)
@btime CUDA.@sync bench_dot_atomic!($x1g, $x2g, $z0, threads)
#  17.323 ms (50 allocations: 1.67 KiB)

function dot_share!(x::CuDeviceVector{T}, y::CuDeviceVector{T}, z::CuDeviceVector{T}) where {T}

    tid = threadIdx().x
    idx = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    shared = CUDA.@cuStaticSharedMem(T, 64)
    fill!(shared, 0)
    sync_threads()

    if idx <= length(x)
        @inbounds shared[tid] = x[idx] * y[idx]
    end
    sync_threads()

    i = blockDim().x ÷ 2
    while i > 0
        if tid <= i
            @inbounds shared[tid] += shared[tid + i] # valid non atomic operation
            # CUDA.atomic_add!(pointer(shared, tid), shared[tid+i]) # invalid op - results in error
        end
        sync_threads()
        i ÷= 2
    end
    if tid == 1
        CUDA.atomic_add!(pointer(z, 1), shared[1])
    end
    return nothing
end

function bench_dot_share!(x, y, z, threads, numblocks)
    CUDA.@sync @cuda threads = threads blocks = numblocks dot_share!(x, y, z)
    return z
end

function wrap_share(x, y, threads)
    numblocks = ceil(Int, N1 / threads)
    z = CUDA.zeros(1)
    x = bench_dot_share!(x, y, z, threads, numblocks)
    return x
end

threads = 64
numblocks = ceil(Int, N1 / threads)
z = CUDA.zeros(1)
@cuda threads = threads blocks = numblocks dot_share!(x1g, x2g, z)
@time CUDA.@sync wrap_share(x1g, x2g, threads)
@btime CUDA.@sync wrap_share($x1g, $x2g, threads)
x = CUDA.@sync wrap_share(x1g, x2g, threads)
x1g' * x2g
@btime x1g' * x2g
@btime x1' * x2

function dot_share2!(x::CuDeviceVector{T}, z::CuDeviceVector{T}) where {T}

    tid = threadIdx().x
    bid = blockIdx().x
    idx = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    shared = CUDA.@cuStaticSharedMem(T, 128)
    fill!(shared, 0f0)
    sync_threads()

    # if idx <= length(x)
    #     shared[tid] = x[idx] * y[idx]
    # end
    # sync_threads()

    i = blockDim().x ÷ 2
    while i > 0
        # if tid <= i
        if tid == 1
            # shared[tid] += shared[tid + i] # valid non atomic operation
            CUDA.atomic_add!(pointer(shared, tid), shared[tid + 1]) # invalid op - results in error
        end
        sync_threads()
        i ÷= 2
    end
    if tid == 1
        z[bid] += shared[1]
    end
    return nothing
end

function bench_dot_share2!(x, z, threads, numblocks)
    CUDA.@sync @cuda threads = threads blocks = numblocks dot_share2!(x, z)
    return sum(z)
end

function wrap_share2(x, threads)
    numblocks = ceil(Int, N1 / threads)
    z = CUDA.zeros(numblocks)
    # x = bench_dot_share2!(x, z, threads, numblocks)
    CUDA.@sync @cuda threads = threads blocks = numblocks dot_share2!(x, z)
    return(x)
end

threads = 128
numblocks = ceil(Int, N1 / threads)
z = CUDA.zeros(numblocks)
sum(z)
x1g' * x2g
CUDA.@sync @cuda threads = threads blocks = numblocks dot_share2!(x1g, x2g, z)
@btime CUDA.@sync wrap_share2($x1g, $x2g, threads)
x = CUDA.@sync wrap_share2(x1g, x2g, threads)

@btime x1g' * x2g

using CUDA
function kernel(x)
    shared = @cuStaticSharedMem(Float32, 2)
    fill!(shared, 1f0)
    sync_threads()
    # @atomic shared[threadIdx().x] += 0f0
    tid = threadIdx().x
    CUDA.atomic_add!(pointer(shared, tid), shared[tid + 1])
    CUDA.atomic_add!(pointer(x, 1), shared[1])
    return
end

x = CUDA.zeros(1)
@cuda kernel(x)
synchronize()


text:
@rule typeof(-)(:out, Marginalisation) (m_in1::Any, m_in2::Any, meta::Any) = begin
    @call_rule typeof(+)(:in1, Marginalisation) (m_out = m_in1, m_in2 = m_in2, meta = meta)
end


text:
# ucsdet.jl - Wrapper for ICU (International Components for Unicode) library

# Some content of the documentation strings was derived from the ICU header files ucsdet.h
# (Those portions copyright (C) 1996-2015, International Business Machines Corporation and others)

"""
   Charset Detection API

   This API provides a facility for detecting the charset or encoding of character data in an unknown text format.
   The input data can be from an array of bytes.

   Character set detection is at best an imprecise operation.  The detection process will attempt to identify the
   charset that best matches the characteristics of the byte data, but the process is partly statistical in nature,
   and the results can not be guaranteed to always be correct.

   For best accuracy in charset detection, the input data should be primarily in a single language, and a minimum of
   a few hundred bytes worth of plain text in the language are needed.  The detection process will attempt to
   ignore html or xml style markup that could otherwise obscure the content.
"""
module ucsdet end

macro libucsdet(s)  ; _libicu(s, iculibi18n, "ucsdet_")   ; end

export UCharsetDetector

const _empty_vec8 = UInt8[]

mutable struct UCharsetDetector
    p::Ptr{Cvoid}
    s

    function UCharsetDetector()
        err = Ref{UErrorCode}(0)
        p = ccall(@libucsdet(open), Ptr{Cvoid}, (Ptr{UErrorCode},), err)
        SUCCESS(err[]) || error("ICU: $(err[]), could not open charset detector")
        self = new(p, _empty_vec8)
        finalizer(self, close)
        self
    end
end
close(c::UCharsetDetector) =
    c.p == C_NULL || (ccall(@libucsdet(close), Cvoid, (Ptr{Cvoid},), c.p); c.p = C_NULL)

"""
   Set the input byte data whose charset is to detected.

   Ownership of the input text byte array remains with the caller.
   The input string must not be altered or deleted until the charset detector is either closed or reset
   to refer to different input text.

   Arguments:
   csd   the charset detector to be used.
   s     the input text of unknown encoding.
"""
function set!(csd::UCharsetDetector, s, p::Ptr{UInt8}, len)
    csd.s = s
    err = Ref{UErrorCode}(0)
    ccall(@libucsdet(setText), Cvoid,
          (Ptr{Cvoid}, Ptr{UInt8}, Int32, Ptr{UErrorCode}), csd.p, p, len, err)
    SUCCESS(err[]) || error("ICU: $(err[]), could not set charset detector text")
    nothing
end
set!(csd::UCharsetDetector, s::Str) =
    @preserve s set!(csd, s, reinterpret(Ptr{UInt8}, pointer(s)), ncodeunits(s))
set!(csd::UCharsetDetector, s::Vector{UInt8}) =
    @preserve s set!(csd, s, pointer(s), length(s))
set!(csd::UCharsetDetector, s::CodeUnits) = set!(csd, Vector{UInt8}(s))
set!(csd::UCharsetDetector, s::AbstractString) = set!(csd, codeunits(s))

"""
   Set the declared encoding for charset detection.
   The declared encoding of an input text is an encoding obtained by the user from an http header or xml
   declaration or similar source that can be provided as an additional hint to the charset detector.

   How and whether the declared encoding will be used during the detection process is TBD.

   Arguments:
   csd      the charset detector to be used.
   encoding encoding for the current data obtained from a header or declaration or other source outside
            of the byte data itself.
"""
function set_declared_encoding!(csd::UCharsetDetector, s::Ptr{UInt8}, len)
    err = Ref{UErrorCode}(0)
    ccall(@libucsdet(setDeclaredEncoding), Cvoid,
          (Ptr{Cvoid},Ptr{UInt8},Int32,Ptr{UErrorCode}), csd.p, s, len, err)
    SUCCESS(err[]) || error("ICU: $(err[]), could not set charset detector declared encoding")
    nothing
end
set_declared_encoding!(csd::UCharsetDetector, s::T) where {T<:Str} =
    @preserve s set_declared_encoding(csd, reinterpret(Ptr{UInt8}, pointer(s)), sizeof(s))
set_declared_encoding!(csd::UCharsetDetector, s::AbstractString) =
    @preserve s set_declared_encoding(csd, pointer(s), sizeof(s))
set_declared_encoding!(csd::UCharsetDetector, s::Vector{T}) where {T<:Union{UInt8,UInt16,UInt32}} =
    @preserve s set_declared_encoding(csd, pointer(s), sizeof(s))

"""
   Opaque structure representing a match that was identified from a charset detection operation.
"""
mutable struct UCharsetMatch
    p::Ptr{Cvoid}
    # This is needed to keep the UCharsetDetector from being GCed until after all dependent UCharsetMatch's are GCed
    ucd::UCharsetDetector
end

"""
   Return the charset that best matches the supplied input data.

   Note though, that because the detection only looks at the start of the input data, there is a possibility
   that the returned charset will fail to handle the full set of input data.

   The returned UCharsetMatch object is owned by the UCharsetDetector.
   It will remain valid until the detector input is reset, or until the detector is closed.

   The function will fail if

no charset appears to match the data.

no input text has been provided

   Arguments:
   csd      the charset detector to be used.
   status   any error conditions are reported back in this variable.

   Returns: a UCharsetMatch representing the best matching charset, or nothing if no charset matches the byte data.
"""
function detect(csd::UCharsetDetector)
    err = Ref{UErrorCode}(0)
    p = ccall(@libucsdet(detect), Ptr{Cvoid}, (Ptr{Cvoid},Ptr{UErrorCode}), csd.p, err)
    SUCCESS(err[]) || error("ICU: $(err[]), could not detect encoding")
    p != C_NULL ? UCharsetMatch(p, csd) : nothing
end

"""
   Find all charset matches that appear to be consistent with the input, returning an array of results.

   The results are ordered with the best quality match first.

   Because the detection only looks at a limited amount of the input byte data, some of the returned charsets may fail
   to handle the all of input data.

   The returned UCharsetMatch objects are owned by the UCharsetDetector.
   They will remain valid until the detector is closed or modified

   Return an error if

no charsets appear to match the input data.

no input text has been provided

   Arguments:
   ucsd          the charset detector to be used.
   matches       Found pointer to a variable that will be set to the number of charsets identified that
                 are consistent with athe input data.  Output only.
   status        any error conditions are reported back in this variable.

   Returns:      A pointer to an array of pointers to UCharSetMatch objects.
                 This array, and the UCharSetMatch instances to which it refers, are owned by the UCharsetDetector,
                 and will remain valid until the detector is closed or modified.
"""
function detectall(csd::UCharsetDetector)
    err = Ref{UErrorCode}(0)
    n = Ref{Int32}(0)
    p = ccall(@libucsdet(detectAll), Ptr{Ptr{Cvoid}},
              (Ptr{Cvoid},Ptr{Int32},Ptr{UErrorCode}), csd.p, n, err)
    SUCCESS(err[]) || error("ICU: $(err[]), could not detect encoding")
    n[] > 0 || return Vector{UCharsetMatch}()

    [UCharsetMatch(x) for x in pointer_to_array(p, int(n[]))]
end

"""
   Get the name of the charset represented by a UCharsetMatch.

   The storage for the returned name string is owned by the UCharsetMatch, and will remain valid
   while the UCharsetMatch is valid.

   The name returned is suitable for use with the ICU conversion APIs.

   Arguments:
   ucsm     The charset match object.
   status   Any error conditions are reported back in this variable.

   Returns: The name of the matching charset.
"""
function get_name(csmatch::UCharsetMatch)
    csmatch.p != C_NULL || throw(UndefRefError())
    err = Ref{UErrorCode}(0)
    name = ccall(@libucsdet(getName), Ptr{UInt8}, (Ptr{Cvoid},Ptr{UErrorCode}), csmatch.p, err)
    SUCCESS(err[]) || error("ICU: $(err[]), could not get name of matching encoding")
    unsafe_string(name)
end

"""
   Get a confidence number for the quality of the match of the byte data with the charset.
   Confidence numbers range from zero to 100, with 100 representing complete confidence and zero
   representing no confidence.

   The confidence values are somewhat arbitrary.  They define an an ordering within the results for any
   single detection operation but are not generally comparable between the results for different input.

   A confidence value of ten does have a general meaning - it is used for charsets that can represent the
   input data, but for which there is no other indication that suggests that the charset is the correct one.
   Pure 7 bit ASCII data, for example, is compatible with a great many charsets, most of which will appear as
   possible matches with a confidence of 10.

   Arguments:
   ucsm     The charset match object.
   status   Any error conditions are reported back in this variable.

   Returns: A confidence number for the charset match.
"""
function get_confidence(csmatch::UCharsetMatch)
    csmatch.p != C_NULL || throw(UndefRefError())
    err = Ref{UErrorCode}(0)
    confidence = ccall(@libucsdet(getConfidence), Int32,
                       (Ptr{Cvoid},Ptr{UErrorCode}), csmatch.p, err)
    SUCCESS(err[]) || error("ICU: $(err[]), could not get confidence")
    Int(confidence)
end

"""
   Get the RFC 3066 code for the language of the input data.

   The Charset Detection service is intended primarily for detecting charsets, not language.
   For some, but not all, charsets, a language is identified as a byproduct of the detection process,
   and that is what is returned by this function.

   CAUTION:
      1.  Language information is not available for input data encoded in all charsets. In particular,
          no language is identified for UTF-8 input data.
      2.  Closely related languages may sometimes be confused.

   If more accurate language detection is required, a linguistic analysis package should be used.

   The storage for the returned name string is owned by the UCharsetMatch, and will remain valid
   while the UCharsetMatch is valid.

   Arguments:
   ucsm    The charset match object.
   status  Any error conditions are reported back in this variable.

   Returns:
   The RFC 3066 code for the language of the input data, or an empty string if the language could not be determined.
"""
function get_language(csmatch::UCharsetMatch)
    csmatch.p != C_NULL || throw(UndefRefError())
    err = Ref{UErrorCode}(0)
    p = ccall(@libucsdet(getLanguage), Ptr{UInt8}, (Ptr{Cvoid},Ptr{UErrorCode}), csmatch.p, err)
    SUCCESS(err[]) || error("ICU: $(err[]), could not get language")
    unsafe_string(p)
end

"""
   Get the entire input text as a UChar string, placing it into a caller-supplied buffer.  A terminating
   NUL character will be appended to the buffer if space is available.

   The number of UChars in the output string, not including the terminating NUL, is returned.

   If the supplied buffer is smaller than required to hold the output, the contents of the buffer are undefined.
   The full output string length (in UChars) is returned as always, and can be used to allocate a buffer of the
   correct size.

   Arguments:
   ucsm     The charset match object.
   buf      A UChar buffer to be filled with the converted text data.
   cap      The capacity of the buffer in UChars.
   status   Any error conditions are reported back in this variable.

   Returns: The number of UChars in the output string.
"""
function get_uchars(csmatch::UCharsetMatch)
    csmatch.p != C_NULL || throw(UndefRefError())
    err = Ref{UErrorCode}(0)
    len = _get_uchars(csmatch, C_NULL, 0, err)
    dest = _allocate(UInt16, len+1)
    len = _get_uchars(csmatch, dest, len+1, err)
    SUCCESS(err[]) || error("ICU: $(err[]), could not get string from UCharsetMatch object")
    Str(UTF16CSE, dest[1:len])
end

_get_uchars(p::Ptr{Cvoid}, buf, siz, err) =
    ccall(@libucsdet(getUChars), Int32,
          (Ptr{Cvoid},Ptr{UInt8},Int32,Ptr{UErrorCode}), p, buf, siz, err)

"""
   Get an iterator over the set of all detectable charsets -
   over the charsets that are known to the charset detection service.

   The returned UEnumeration provides access to the names of the charsets.

   The state of the Charset detector that is passed in does not affect the result of this function,
   but requiring a valid, open charset detector as a parameter insures that the charset detection
   service has been safely initialized and that the required detection data is available.

   Note: Multiple different charset encodings in a same family may use a single shared name
   in this implementation. For example, this method returns an array including "ISO-8859-1" (ISO Latin 1),
   but not including "windows-1252" (Windows Latin 1). However, actual detection result could be "windows-1252"
   when the input data matches Latin 1 code points with any points only available in "windows-1252".

   Arguments:
   ucsd      Charset detector.

   Returns:  iterator providing access to the detectable charset names.
"""
function get_all_detectable_charsets(csd::UCharsetDetector)
    csd.p != C_NULL || throw(UndefRefError())
    err = Ref{UErrorCode}(0)
    p = ccall(@libucsdet(getAllDetectableCharsets), Ptr{Cvoid},
              (Ptr{Cvoid}, Ptr{UErrorCode}), csd.p, err)
    SUCCESS(err[]) || error("ICU: $(err[]), could not get all detected charsets")
    p
end

"""
    Test whether input filtering is enabled for this charset detector.

    Input filtering removes text that appears to be HTML or xml markup from the input
    before applying the code page detection heuristics.

    Arguments:
    ucsd   The charset detector to check.

    Returns: TRUE if filtering is enabled.
"""
input_filter_enabled(csd::UCharsetDetector) =
    ccall(@libucsdet(isInputFilterEnabled), UBool, (Ptr{Cvoid},), csd.p)

"""
   Enable filtering of input text. If filtering is enabled, text within angle brackets ("<" and ">") will be removed
   before detection, which will remove most HTML or xml markup.

   Arguments:
   ucsd    the charset detector to be modified.
   filter  true to enable input text filtering.

   Returns: The previous setting.
"""
enable_input_filter(csd::UCharsetDetector, filter::Bool=true) =
    ccall(@libucsdet(enableInputFilter), UBool, (Ptr{Cvoid},Bool), csd.p, filter)


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

__precompile__(true)

"""
Arnoldi and Lanczos iteration for computing eigenvalues
"""
module IterativeEigenSolvers

using Base.LinAlg: BlasFloat, BlasInt, SVD, checksquare

export eigs, svds

include("arpack.jl")

using .ARPACK

## eigs
"""
    eigs(A; nev=6, ncv=max(20,2*nev+1), which=:LM, tol=0.0, maxiter=300, sigma=nothing, ritzvec=true, v0=zeros((0,))) -> (d,[v,],nconv,niter,nmult,resid)

Computes eigenvalues `d` of `A` using implicitly restarted Lanczos or Arnoldi iterations for real symmetric or
general nonsymmetric matrices respectively.

The following keyword arguments are supported:

* `nev`: Number of eigenvalues
* `ncv`: Number of Krylov vectors used in the computation; should satisfy `nev+1 <= ncv <= n`
  for real symmetric problems and `nev+2 <= ncv <= n` for other problems, where `n` is the
  size of the input matrix `A`. The default is `ncv = max(20,2*nev+1)`. Note that these
  restrictions limit the input matrix `A` to be of dimension at least 2.
* `which`: type of eigenvalues to compute. See the note below.

| `which` | type of eigenvalues                                                                                |
|:--------|:--------------------------------------------------------------------------------------------------------------------------|
| `:LM`   | eigenvalues of largest magnitude (default)                                                                                |
| `:SM`   | eigenvalues of smallest magnitude                                                                                |
| `:LR`   | eigenvalues of largest real part                                                                                |
| `:SR`   | eigenvalues of smallest real part                                                                                |
| `:LI`   | eigenvalues of largest imaginary part (nonsymmetric or complex `A` only)                                                  |
| `:SI`   | eigenvalues of smallest imaginary part (nonsymmetric or complex `A` only)                                                 |
| `:BE`   | compute half of the eigenvalues from each end of the spectrum, biased in favor of the high end. (real symmetric `A` only) |

* `tol`: parameter defining the relative tolerance for convergence of Ritz values (eigenvalue estimates).
     A Ritz value ``θ`` is considered converged when its associated residual
     is less than or equal to the product of `tol` and ``max(ɛ^{2/3}, |θ|)``,
     where `ɛ = eps(real(eltype(A)))/2` is LAPACK's machine epsilon.
     The residual associated with ``θ`` and its corresponding Ritz vector ``v``
     is defined as the norm ``||Av - vθ||``.
     The specified value of `tol` should be positive; otherwise, it is ignored
     and ``ɛ`` is used instead.
     Default: ``ɛ``.

* `maxiter`: Maximum number of iterations (default = 300)
* `sigma`: Specifies the level shift used in inverse iteration. If `nothing` (default),
  defaults to ordinary (forward) iterations. Otherwise, find eigenvalues close to `sigma`
  using shift and invert iterations.
* `ritzvec`: Returns the Ritz vectors `v` (eigenvectors) if `true`
* `v0`: starting vector from which to start the iterations

`eigs` returns the `nev` requested eigenvalues in `d`, the corresponding Ritz vectors `v`
(only if `ritzvec=true`), the number of converged eigenvalues `nconv`, the number of
iterations `niter` and the number of matrix vector multiplications `nmult`, as well as the
final residual vector `resid`.

# Examples
```jldoctest
julia> A = Diagonal(1:4);

julia> λ, ϕ = eigs(A, nev = 2);

julia> λ
2-element Array{Float64,1}:
 4.0
 3.0
```

!!! note
    The `sigma` and `which` keywords interact: the description of eigenvalues
    searched for by `which` do *not* necessarily refer to the eigenvalues of
    `A`, but rather the linear operator constructed by the specification of the
    iteration mode implied by `sigma`.

    | `sigma`         | iteration mode                   | `which` refers to eigenvalues of |
    |:----------------|:---------------------------------|:---------------------------------|
    | `nothing`       | ordinary (forward)               | ``A``                            |
    | real or complex | inverse with level shift `sigma` | ``(A - \\sigma I )^{-1}``        |

!!! note
    Although `tol` has a default value, the best choice depends strongly on the
    matrix `A`. We recommend that users _always_ specify a value for `tol`
    which suits their specific needs.

    For details of how the errors in the computed eigenvalues are estimated, see:

    * B. N. Parlett, "The Symmetric Eigenvalue Problem", SIAM: Philadelphia, 2/e
      (1998), Ch. 13.2, "Accessing Accuracy in Lanczos Problems", pp. 290-292 ff.
    * R. B. Lehoucq and D. C. Sorensen, "Deflation Techniques for an Implicitly
      Restarted Arnoldi Iteration", SIAM Journal on Matrix Analysis and
      Applications (1996), 17(4), 789–821.  doi:10.1137/S0895479895281484
"""
eigs(A; kwargs...) = eigs(A, I; kwargs...)
eigs(A::AbstractMatrix{<:BlasFloat}, ::UniformScaling; kwargs...) = _eigs(A, I; kwargs...)

eigs(A::AbstractMatrix{T}, B::AbstractMatrix{T}; kwargs...) where {T<:BlasFloat} = _eigs(A, B; kwargs...)
eigs(A::AbstractMatrix{BigFloat}, B::AbstractMatrix...; kwargs...) = throw(MethodError(eigs, Any[A,B,kwargs...]))
eigs(A::AbstractMatrix{BigFloat}, B::UniformScaling; kwargs...) = throw(MethodError(eigs, Any[A,B,kwargs...]))
function eigs(A::AbstractMatrix{T}, ::UniformScaling; kwargs...) where T
    Tnew = typeof(zero(T)/sqrt(one(T)))
    eigs(convert(AbstractMatrix{Tnew}, A), I; kwargs...)
end
function eigs(A::AbstractMatrix, B::AbstractMatrix; kwargs...)
    T = promote_type(eltype(A), eltype(B))
    Tnew = typeof(zero(T)/sqrt(one(T)))
    eigs(convert(AbstractMatrix{Tnew}, A), convert(AbstractMatrix{Tnew}, B); kwargs...)
end
"""
    eigs(A, B; nev=6, ncv=max(20,2*nev+1), which=:LM, tol=0.0, maxiter=300, sigma=nothing, ritzvec=true, v0=zeros((0,))) -> (d,[v,],nconv,niter,nmult,resid)

Computes generalized eigenvalues `d` of `A` and `B` using implicitly restarted Lanczos or Arnoldi iterations for
real symmetric or general nonsymmetric matrices respectively.

The following keyword arguments are supported:

* `nev`: Number of eigenvalues
* `ncv`: Number of Krylov vectors used in the computation; should satisfy `nev+1 <= ncv <= n`
  for real symmetric problems and `nev+2 <= ncv <= n` for other problems, where `n` is the
  size of the input matrices `A` and `B`. The default is `ncv = max(20,2*nev+1)`. Note that
  these restrictions limit the input matrix `A` to be of dimension at least 2.
* `which`: type of eigenvalues to compute. See the note below.

| `which` | type of eigenvalues                                                                                |
|:--------|:--------------------------------------------------------------------------------------------------------------------------|
| `:LM`   | eigenvalues of largest magnitude (default)                                                                                |
| `:SM`   | eigenvalues of smallest magnitude                                                                                |
| `:LR`   | eigenvalues of largest real part                                                                                |
| `:SR`   | eigenvalues of smallest real part                                                                                |
| `:LI`   | eigenvalues of largest imaginary part (nonsymmetric or complex `A` only)                                                  |
| `:SI`   | eigenvalues of smallest imaginary part (nonsymmetric or complex `A` only)                                                 |
| `:BE`   | compute half of the eigenvalues from each end of the spectrum, biased in favor of the high end. (real symmetric `A` only) |

* `tol`: relative tolerance used in the convergence criterion for eigenvalues, similar to
     `tol` in the [`eigs(A)`](@ref) method for the ordinary eigenvalue
     problem, but effectively for the eigenvalues of ``B^{-1} A`` instead of ``A``.
     See the documentation for the ordinary eigenvalue problem in
     [`eigs(A)`](@ref) and the accompanying note about `tol`.
* `maxiter`: Maximum number of iterations (default = 300)
* `sigma`: Specifies the level shift used in inverse iteration. If `nothing` (default),
  defaults to ordinary (forward) iterations. Otherwise, find eigenvalues close to `sigma`
  using shift and invert iterations.
* `ritzvec`: Returns the Ritz vectors `v` (eigenvectors) if `true`
* `v0`: starting vector from which to start the iterations

`eigs` returns the `nev` requested eigenvalues in `d`, the corresponding Ritz vectors `v`
(only if `ritzvec=true`), the number of converged eigenvalues `nconv`, the number of
iterations `niter` and the number of matrix vector multiplications `nmult`, as well as the
final residual vector `resid`.

# Examples
```jldoctest
julia> A = sparse(1.0I, 4, 4); B = Diagonal(1:4);

julia> λ, ϕ = eigs(A, B, nev = 2);

julia> λ
2-element Array{Float64,1}:
 1.0
 0.4999999999999999
```

!!! note
    The `sigma` and `which` keywords interact: the description of eigenvalues searched for by
    `which` do *not* necessarily refer to the eigenvalue problem ``Av = Bv\\lambda``, but rather
    the linear operator constructed by the specification of the iteration mode implied by `sigma`.

    | `sigma`         | iteration mode                   | `which` refers to the problem      |
    |:----------------|:---------------------------------|:-----------------------------------|
    | `nothing`       | ordinary (forward)               | ``Av = Bv\\lambda``                |
    | real or complex | inverse with level shift `sigma` | ``(A - \\sigma B )^{-1}B = v\\nu`` |
"""
eigs(A, B; kwargs...) = _eigs(A, B; kwargs...)
function _eigs(A, B;
               nev::Integer=6, ncv::Integer=max(20,2*nev+1), which=:LM,
               tol=0.0, maxiter::Integer=300, sigma=nothing, v0::Vector=zeros(eltype(A),(0,)),
               ritzvec::Bool=true)
    n = checksquare(A)

    T = eltype(A)
    iscmplx = T <: Complex
    isgeneral = B !== I
    sym = !iscmplx && issymmetric(A) && issymmetric(B)
    nevmax = sym ? n-1 : n-2
    if nevmax <= 0
        throw(ArgumentError("input matrix A is too small. Use eigfact instead."))
    end
    if nev > nevmax
        warn("Adjusting nev from $nev to $nevmax")
        nev = nevmax
    end
    if nev <= 0
        throw(ArgumentError("requested number of eigenvalues (nev) must be ≥ 1, got $nev"))
    end
    ncvmin = nev + (sym ? 1 : 2)
    if ncv < ncvmin
        warn("Adjusting ncv from $ncv to $ncvmin")
        ncv = ncvmin
    end
    ncv = BlasInt(min(ncv, n))
    bmat = isgeneral ? "G" : "I"
    isshift = sigma !== nothing

    if isa(which,AbstractString)
        warn("Use symbols instead of strings for specifying which eigenvalues to compute")
        which=Symbol(which)
    end
    if (which != :LM && which != :SM && which != :LR && which != :SR &&
        which != :LI && which != :SI && which != :BE)
        throw(ArgumentError("which must be :LM, :SM, :LR, :SR, :LI, :SI, or :BE, got $(repr(which))"))
    end
    if which == :BE && !sym
        throw(ArgumentError("which=:BE only possible for real symmetric problem"))
    end
    isshift && which == :SM && warn("use of :SM in shift-and-invert mode is not recommended, use :LM to find eigenvalues closest to sigma")

    if which==:SM && !isshift # transform into shift-and-invert method with sigma = 0
        isshift=true
        sigma=zero(T)
        which=:LM
    end

    if sigma !== nothing && !iscmplx && isa(sigma,Complex)
        throw(ArgumentError("complex shifts for real problems are not yet supported"))
    end
    sigma = isshift ? convert(T,sigma) : zero(T)

    if !isempty(v0)
        if length(v0) != n
            throw(DimensionMismatch())
        end
        if eltype(v0) != T
            throw(ArgumentError("starting vector must have element type $T, got $(eltype(v0))"))
        end
    end

    whichstr = "LM"
    if which == :BE
        whichstr = "BE"
    end
    if which == :LR
        whichstr = (!sym ? "LR" : "LA")
    end
    if which == :SR
        whichstr = (!sym ? "SR" : "SA")
    end
    if which == :LI
        if !sym
            whichstr = "LI"
        else
            throw(ArgumentError("largest imaginary is meaningless for symmetric eigenvalue problems"))
        end
    end
    if which == :SI
        if !sym
            whichstr = "SI"
        else
            throw(ArgumentError("smallest imaginary is meaningless for symmetric eigenvalue problems"))
        end
    end

    # Refer to ex-*.doc files in ARPACK/DOCUMENTS for calling sequence
    matvecA!(y, x) = A_mul_B!(y, A, x)
    if !isgeneral           # Standard problem
        matvecB = x -> x
        if !isshift         #    Regular mode
            mode       = 1
            solveSI = x->x
        else                #    Shift-invert mode
            mode       = 3
            F = factorize(A - UniformScaling(sigma))
            solveSI = x -> F \ x
        end
    else                    # Generalized eigenproblem
        matvecB = x -> B * x
        if !isshift         #    Regular inverse mode
            mode       = 2
            F = factorize(B)
            solveSI = x -> F \ x
        else                #    Shift-invert mode
            mode       = 3
            F = factorize(A - sigma*B)
            solveSI = x -> F \ x
        end
    end

    # Compute the Ritz values and Ritz vectors
    (resid, v, ldv, iparam, ipntr, workd, workl, lworkl, rwork, TOL) =
       ARPACK.aupd_wrapper(T, matvecA!, matvecB, solveSI, n, sym, iscmplx, bmat, nev, ncv, whichstr, tol, maxiter, mode, v0)

    # Postprocessing to get eigenvalues and eigenvectors
    output = ARPACK.eupd_wrapper(T, n, sym, iscmplx, bmat, nev, whichstr, ritzvec, TOL,
                                 resid, ncv, v, ldv, sigma, iparam, ipntr, workd, workl, lworkl, rwork)

    # Issue 10495, 10701: Check that all eigenvalues are converged
    nev = length(output[1])
    nconv = output[ritzvec ? 3 : 2]
    nev ≤ nconv || warn("not all wanted Ritz pairs converged. Requested: $nev, converged: $nconv")

    return output
end

## svds
struct SVDAugmented{T,S} <: AbstractArray{T, 2}
    X::S
    SVDAugmented{T,S}(X::AbstractMatrix) where {T,S} = new(X)
end

function SVDAugmented(A::AbstractMatrix{T}) where T
    Tnew = typeof(zero(T)/sqrt(one(T)))
    Anew = convert(AbstractMatrix{Tnew}, A)
    SVDAugmented{Tnew,typeof(Anew)}(Anew)
end

function Base.A_mul_B!(y::StridedVector{T}, A::SVDAugmented{T}, x::StridedVector{T}) where T
    m, mn = size(A.X, 1), length(x)
    A_mul_B!( view(y, 1:m), A.X, view(x, m + 1:mn)) # left singular vector
    Ac_mul_B!(view(y, m + 1:mn), A.X, view(x, 1:m)) # right singular vector
    return y
end
Base.size(A::SVDAugmented)  = ((+)(size(A.X)...), (+)(size(A.X)...))
Base.ishermitian(A::SVDAugmented) = true

struct AtA_or_AAt{T,S} <: AbstractArray{T, 2}
    A::S
    buffer::Vector{T}
end

function AtA_or_AAt(A::AbstractMatrix{T}) where T
    Tnew = typeof(zero(T)/sqrt(one(T)))
    Anew = convert(AbstractMatrix{Tnew}, A)
    AtA_or_AAt{Tnew,typeof(Anew)}(Anew, Vector{Tnew}(uninitialized, max(size(A)...)))
end

function Base.A_mul_B!(y::StridedVector{T}, A::AtA_or_AAt{T}, x::StridedVector{T}) where T
    if size(A.A, 1) >= size(A.A, 2)
        A_mul_B!(A.buffer, A.A, x)
        return Ac_mul_B!(y, A.A, A.buffer)
    else
        Ac_mul_B!(A.buffer, A.A, x)
        return A_mul_B!(y, A.A, A.buffer)
    end
end
Base.size(A::AtA_or_AAt) = ntuple(i -> min(size(A.A)...), Val(2))
Base.ishermitian(s::AtA_or_AAt) = true

svds(A::AbstractMatrix{<:BlasFloat}; kwargs...) = _svds(A; kwargs...)
svds(A::AbstractMatrix{BigFloat}; kwargs...) = throw(MethodError(svds, Any[A, kwargs...]))
function svds(A::AbstractMatrix{T}; kwargs...) where T
    Tnew = typeof(zero(T)/sqrt(one(T)))
    svds(convert(AbstractMatrix{Tnew}, A); kwargs...)
end

"""
    svds(A; nsv=6, ritzvec=true, tol=0.0, maxiter=1000, ncv=2*nsv, v0=zeros((0,))) -> (SVD([left_sv,] s, [right_sv,]), nconv, niter, nmult, resid)

Computes the largest singular values `s` of `A` using implicitly restarted Lanczos
iterations derived from [`eigs`](@ref).

**Inputs**

* `A`: Linear operator whose singular values are desired. `A` may be represented as a
  subtype of `AbstractArray`, e.g., a sparse matrix, or any other type supporting the four
  methods `size(A)`, `eltype(A)`, `A * vector`, and `A' * vector`.
* `nsv`: Number of singular values. Default: 6.
* `ritzvec`: If `true`, return the left and right singular vectors `left_sv` and `right_sv`.
   If `false`, omit the singular vectors. Default: `true`.
* `tol`: tolerance, see [`eigs`](@ref).
* `maxiter`: Maximum number of iterations, see [`eigs`](@ref). Default: 1000.
* `ncv`: Maximum size of the Krylov subspace, see [`eigs`](@ref) (there called `nev`). Default: `2*nsv`.
* `v0`: Initial guess for the first Krylov vector. It may have length `min(size(A)...)`, or 0.

**Outputs**

* `svd`: An `SVD` object containing the left singular vectors, the requested values, and the
  right singular vectors. If `ritzvec = false`, the left and right singular vectors will be
  empty.
* `nconv`: Number of converged singular values.
* `niter`: Number of iterations.
* `nmult`: Number of matrix--vector products used.
* `resid`: Final residual vector.

# Examples
```jldoctest
julia> A = Diagonal(1:4);

julia> s = svds(A, nsv = 2)[1];

julia> s[:S]
2-element Array{Float64,1}:
 4.0
 2.9999999999999996
```

!!! note "Implementation"
    `svds(A)` is formally equivalent to calling [`eigs`](@ref) to perform implicitly restarted
    Lanczos tridiagonalization on the Hermitian matrix ``A^\\prime A`` or ``AA^\\prime`` such
    that the size is smallest.
"""
svds(A; kwargs...) = _svds(A; kwargs...)
function _svds(X; nsv::Int = 6, ritzvec::Bool = true, tol::Float64 = 0.0, maxiter::Int = 1000, ncv::Int = 2*nsv, v0::Vector=zeros(eltype(X),(0,)))
    if nsv < 1
        throw(ArgumentError("number of singular values (nsv) must be ≥ 1, got $nsv"))
    end
    if nsv > minimum(size(X))
        throw(ArgumentError("number of singular values (nsv) must be ≤ $(minimum(size(X))), got $nsv"))
    end
    m, n = size(X)
    otype = eltype(X)
    if length(v0) ∉ [0,n]
        throw(DimensionMismatch("length of v0, the guess for the starting right Krylov vector, must be 0, or $n, got $(length(v0))"))
    end
    ex    = eigs(AtA_or_AAt(X), I; which = :LM, ritzvec = ritzvec, nev = nsv, tol = tol, maxiter = maxiter, v0=v0)
    # ind   = [1:2:ncv;]
    # sval  = abs.(ex[1][ind])

    svals = sqrt.(real.(ex[1]))
    if ritzvec
        # calculating singular vectors
        # left_sv  = sqrt(2) * ex[2][ 1:size(X,1),     ind ] .* sign.(ex[1][ind]')
        if size(X, 1) >= size(X, 2)
            V = ex[2]
            U = qr(scale!(X*V, inv.(svals)))[1]
        else
            U = ex[2]
            V = qr(scale!(X'U, inv.(svals)))[1]
        end

        # right_sv = sqrt(2) * ex[2][ size(X,1)+1:end, ind ]
        return (SVD(U, svals, V'), ex[3], ex[4], ex[5], ex[6])
    else
        #The sort is necessary to work around #10329
        return (SVD(zeros(eltype(svals), n, 0),
                    svals,
                    zeros(eltype(svals), 0, m)),
                    ex[2], ex[3], ex[4], ex[5])
    end
end

end # module

max_stars_repo_name: ihnorton/julia

max_issues_repo_name: ihnorton/julia

max_forks_repo_name: ihnorton/julia


text:
module EnsembleSampler

using Distributed
using NPZ
using Logging

# This is a direct Julia port of [emcee](http://dan.iel.fm/emcee/current/), the ensemble sampler by Dan Foreman-Mackey et al.

export Sampler, run_mcmc, run_schedule, reset_mcmc, write_samples, emcee_chain

# There is a type, called the sampler.
mutable struct Sampler
    nwalkers::Int
    ndim::Int
    lnprobfn::Function
    a::Float64
    chain # Array of shape (nwalkers, iterations, dim) in python.
    # For us, this is (ndim, iterations, nwalkers) in Julia, since we will be filling in
    # an entire row of ndim at once for a single iteration, making it the fastest changing dimension.
    lnprob # storing the values of the lnprob. # lnprob in python is (nwalkers, niterations)
    # in julia, it's the same
    iterations::Int # How many iterations has the ensemble taken
end

function Sampler(nwalkers::Int, ndim::Int, lnprobfn::Function, test::Bool = false)
    a = 2.0

    if !test
        @assert ndim < nwalkers "I don't believe you have more parameters than walkers! Try flipping them."
    end

    chain = Array{Float64}(undef, (ndim, 0, nwalkers))
    lnprob = Array{Float64}(undef, (nwalkers, 0))
    iterations = 0

    sampler = Sampler(nwalkers, ndim, lnprobfn, a, chain, lnprob, iterations)
end

# Transpose the chain order from Julia column-major to Python row-major
function emcee_chain(sampler::Sampler)
    # ndim, niter, nwalkers = size(sampler.chain)
    # return reshape(sampler.chain, (nwalkers, niter, ndim))
    return permutedims(sampler.chain, [3, 2, 1])
end

# Clear the samples after burn in
function reset_mcmc(sampler::Sampler)
    # self.naccepted = np.zeros(self.k)
    sampler.chain = Array{Float64}(sampler.ndim, 0, sampler.nwalkers)
    sampler.lnprob = Array{Float64}(sampler.nwalkers, 0)
    sampler.iterations = 0
end

function sample(sampler::Sampler, p0, lnprob0 = nothing, iterations = 1)

    p = p0

    # What is the index to divide the number of walkers in half
    halfk = floor(Int, sampler.nwalkers / 2)

    # If the initial log-probabilities were not provided, calculate them
    # now.
    # This step is calculating the lnprob for the total number of walkers.
    lnprob = lnprob0
    if lnprob == nothing
        lnprob = get_lnprob(sampler, p)
    end

    # Check to make sure that the probability function didn't return
    # ``np.nan``.
    if any(isnan.(lnprob))
        println("The initial lnprob was NaN.")
        throw(DomainError())
    end

    # Store the initial size of the stored chain.
    i0 = size(sampler.chain)[2]
    # resize the chain by adding new rows to the array.

    # (ndim, iterations, nwalkers)
    sampler.chain = cat(sampler.chain, Array{Float64}(undef, sampler.ndim, iterations, sampler.nwalkers), dims = 2)

    # (niterations, nwalkers)
    sampler.lnprob = cat(sampler.lnprob, Array{Float64}(undef, sampler.nwalkers, iterations), dims = 2)

    for i = 1:iterations
        sampler.iterations += 1
        @info "Iteration " sampler.iterations
        # Loop over the two ensembles, calculating the proposed positions.

        # Slices for the first and second halves
        first = 1:halfk
        second = (halfk + 1):sampler.nwalkers

        # Do this for both halves, alternating
        for (S0, S1) in [(first, second), (second, first)]
            q, newlnp, acc = propose_stretch(sampler, p[:, S0], p[:, S1], lnprob[S0])

            if any(acc)
                # Update the positions and log probabilities

                # Note that the indexing needed to be done differently than emcee, since double
                # indexing in Julia does not address the same memory in the array.
                lnprob[S0[acc]] = newlnp[acc]
                p[:, S0[acc]] = q[:, acc]

            end

            ind = i0 + i
            sampler.chain[:, ind, :] = p
            sampler.lnprob[:, ind] = lnprob

        end

    end
    # Return the current position
    return p
end

# Propose a new position for one sub-ensemble given the positions of
# another.
# :param p0:
#  The positions from which to jump.
# :param p1:
#  The positions of the other ensemble.
# :param lnprob0:
#  The log-probabilities at ``p0``.
# This method returns:
# * ``q`` - The new proposed positions for the walkers in ``ensemble``.
# * ``newlnprob`` - The vector of log-probabilities at the positions
# given by ``q``.
# * ``accept`` - A vector of type ``bool`` indicating whether or not
# the proposed position for each walker should be accepted..

function propose_stretch(sampler::Sampler, p0, p1, lnprob0)
    # In python, these are (nwalkers, ndim)
    # In Julia, these are (ndim, nwalkers)

    s = p0
    Ns = size(s)[2]
    c = p1
    Nc = size(c)[2]

    # Ns and Nc refer to the number of walkers in the subset and complimentary set, respectively

    # Generate the vectors of random numbers that will produce the
    # proposal.

    # self._random.rand(Ns) provides a 1D array of values in the range [0,1.) of size Ns
    zz = ((sampler.a - 1.0) .* rand(Ns) .+ 1.0).^2.0 ./ sampler.a

    # An array of random integers the size of the subset, designed to slice into the complimentary sample.
    rint = rand(1:Nc, Ns)

    # Index into the complementary sample, and advance based upon the scale factor.
    # c[:, rint] is a 2D array (ndim, Ns)
    q = c[:, rint] - (reshape(zz, (1, Ns)) .*  (c[:, rint] - s))

    # This get_lnprob is only calculating the lnprob for the proposed positions, which is Ns, and is probably half the total number of walkers.
    newlnprob = get_lnprob(sampler, q)

    # Decide whether or not the proposals should be accepted.
    lnpdiff = (sampler.ndim - 1.0) .* log.(zz) .+ newlnprob .- lnprob0
    accept = (lnpdiff .> log.(rand(Ns)))

    return q, newlnprob, accept

end

# Calculate the vector of log-probability for the walkers.
# :param pos: (optional)
#     The position vector in parameter space where the probability
#     should be calculated. This defaults to the current position
#     unless a different one is provided.

#     (ndim, nwalkers)
# This method returns:
# * ``lnprob`` - A vector of log-probabilities with one entry for each
#   walker in this sub-ensemble.
function get_lnprob(sampler::Sampler, pos)

    p = pos

    # Check that the parameters are in physical ranges.
    if any(isinf.(p))
        println("At least one parameter value was infinite.")
        throw(DomainError())
    end
    if any(isnan.(p))
        println("At least one parameter value was NaN.")
        throw(DomainError())
    end

    # What is the shape of p here in this subset?
    nwalkers_sub = size(pos)[2]
    lst = [p[:,i] for i = 1:nwalkers_sub]

    # In Python, it seems like each row corresponds to a different walker.
    # For Julia, we really want each column to the parameters corresponding to a different walker.
    result = pmap(sampler.lnprobfn, lst)

    # The array may contain one or two RemoteExceptions, so let's check to see what caused them.
    for (res, par) in zip(result, lst)
        if typeof(res) == RemoteException
            println("RemoteException found for parameters ", par)
            println(res)
            throw(DomainError)
        end
    end

    # If we've made it to here, everything's ok.
    lnprob = convert(Array{Float64,1}, result)

    # Check for lnprob returning NaN.
    if any(isnan.(lnprob))
        # Print some debugging stuff.
        println("NaN value of lnprob for parameters: ")
        for pars in p[isnan.(lnprob)]
            println(pars)
        end
        # Finally raise exception.
        throw(DomainError())
    end

    return lnprob

end

# Advance the Sampler for N iterations, staring from pos0.
function run_mcmc(sampler::Sampler, pos0, N::Int)
    pos = sample(sampler, pos0, nothing, N)
    return pos
end

# Write out the samples
function write_samples(sampler::Sampler, outdir = "")

    # fchain = flatchain(sampler)

    npzwrite(joinpath(outdir, "chain.npy"), emcee_chain(sampler))
    # npzwrite(outdir * "flatchain.npy", fchain)
    npzwrite(joinpath(outdir, "lnprob.npy"), sampler.lnprob)

    # Needs to be reshaped to remove singleton dimension
    npzwrite(joinpath(outdir, "pos0.npy"), reshape(sampler.chain[:, end, :], (sampler.ndim, sampler.nwalkers)))

end

# Run the sampler on a periodic save schedule, to prevent losses on long-running calculations
# run N iterations for each loop
function run_schedule(sampler::Sampler, pos0, N::Int, loops::Int, outdir, func = nothing)
    for i = 1:loops
        pos0 = sample(sampler, pos0, nothing, N)
        @info "Finished loop $i of $loops"
        write_samples(sampler, outdir)
        if func != nothing
            func(sampler, outdir)
        end

    end
    return pos0
end

end # module


text:
export BioviewError

type BioviewError <: Exception
    message::AbstractString
end

Base.show(io::IO, e::BioviewError) = show(io, "Bioview error: $(e.message)")


text:
module OdinSon

using PyPlot
# PyPlot exports so only OdinSon is needed for an using import
#TODO: thin this out to only be useful plotting features
export acorr,annotate,arrow,autoscale,autumn,axes,axhline,axhspan,axis,axvline,axvspan
export bar,barbs,barh,bone,box,boxplot,broken_barh,cla,clabel,clf,clim,cohere,colorbar
export colors,contour,contourf,cool,copper,csd,delaxes,disconnect,draw,errorbar,eventplot
export figaspect,figimage,figlegend,figtext,figure,fill_between,fill_betweenx,findobj
export flag,gca,gcf,gci,get_current_fig_manager,get_figlabels,get_fignums,get_plot_commands
export ginput,gray,grid,hexbin,hist2D,hlines,hold,hot,hsv,imread,imsave,imshow,ioff,ion
export ishold,jet,legend,locator_params,loglog,margins,matshow,minorticks_off,minorticks_on
export over,pause,pcolor,pcolormesh,pie,pink,plot,plot_date,plotfile,polar,prism,psd
export quiver,quiverkey,rc,rc_context,rcdefaults,rgrids,savefig,sca,scatter,sci
export semilogx,semilogy,set_cmap,setp,show,specgram,spectral,spring,spy,stackplot
export stem,step,streamplot,subplot,subplot2grid,subplot_tool,subplots,subplots_adjust
export summer,suptitle,table,text,thetagrids,tick_params,ticklabel_format,tight_layout
export title,tricontour,tricontourf,tripcolor,triplot,twinx,twiny,vlines
export winter,xkcd,xlabel,xlim,xscale,xticks,ylabel,ylim,yscale,yticks
# 3D exports
export plot3D
# some nice utilities from PyPlot
export @L_str

# core types
include("coretypes.jl")
export Point

# Colors
import ColorTypes: HSV, RGB, HSL
export HSV, RGB, HSL
include("color_data/crayon.jl")
include("color_data/xkcd.jl")
include("color_data/svg.jl")
include("color.jl")
export PyObject # my overload so that RGB types are turned to python tuples
export crayons, xkcd_rgb, svg_rgb
export @NC_str, @HC_str
# Seaborn like exports TODO: make these names more julian
export set_palette, color_cycle

include("util.jl")
export despine!
# setup the default styling
_set_mpl_kws()

# Styling
include("style.jl")
export Style

# Containers
include("containers.jl")
export Canvas, Viewport

# Basic SVG/Mathematica like language
include("shapes.jl")
export PolyLine, Polygon, Circle

# basic method for all graphic types
export render

end # module

max_stars_repo_name: gabrielgellner/OdinSon

max_issues_repo_name: gabrielgellner/OdinSon

max_forks_repo_name: gabrielgellner/OdinSon


text: include("/opt/bcbi/predictmd/bin/startup.jl");


text:
# Eventually, we will hopefully support multiarch MachO files
@io struct MachOFatArch
    cputype::UInt32
    cpusubtype::UInt32
    offset::UInt32
    size::UInt32
    align::UInt32
end

struct MachOFatHeader{H <: ObjectHandle} <: MachOHeader{H}
    archs::Vector{MachOFatArch}
end