text:
############################ Gradient Algebra ##################################
################################### Sum ########################################
# allocates space for gradient kernel evaluation but does not evaluate
# the separation from evaluation is useful for ValueGradientKernel
function gradient_kernel(k::Sum, x, y, ::GenericInput)
    H = [gradient_kernel(h, x, y, input_trait(h)) for h in k.args]
    LazyMatrixSum(H)
end

# NOTE: K should have been previously allocated with gradient_kernel(k, x, y)
function gradient_kernel!(K::LazyMatrixSum, k::Sum, x::AbstractVector, y::AbstractVector, ::GenericInput)
    for i in eachindex(k.args)
        h = k.args[i]
        K.args[i] = gradient_kernel!(K.args[i], h, x, y, input_trait(h))
    end
    return K
end

function gramian(k::GradientKernel{<:Any, <:Sum}, x::AbstractVector, y::AbstractVector)
    gramian(k, x, y, input_trait(k))
end

# this is more efficient for global structure if the constituents are structured (saves allocations)
function gramian(k::GradientKernel{<:Any, <:Sum}, x::AbstractVector, y::AbstractVector, ::GenericInput)
    LazyMatrixSum((gramian(GradientKernel(ki), x, y) for ki in k.k.args)...)
end

# if all constituents have the same structue, use lazy representation
function gramian(k::GradientKernel{<:Any, <:Sum}, x::AbstractVector, y::AbstractVector, ::InputTrait)
    G = Gramian(k, x, y)
    BlockFactorization(G, isstrided = true)
end

################################ Product #######################################
# for product kernel with generic input
function allocate_gradient_kernel(k::Product, x, y, ::GenericInput)
    d, r = length(x), length(k.args)
    H = (gradient_kernel(h, x, y, input_trait(h)) for h in k.args)
    T = typeof(k(x, y))
    A = LazyMatrixSum(
                (LazyMatrixProduct([Diagonal(zeros(T, d)), h]) for h in H)...
                )
    U = zeros(T, (d, r)) # storage for Jacobian
    V = zeros(T, (d, r))
    C = Woodbury((-I)(r), ones(r), ones(r)')
    W = Woodbury(A, U, C, V')
end
function gradient_kernel(k::Product, x, y, ::GenericInput)
    W = allocate_gradient_kernel(k, x, y, GenericInput())
    gradient_kernel!(W, k, x, y, GenericInput())
end

function gradient_kernel!(W::Woodbury, k::Product, x::AbstractVector, y::AbstractVector, ::GenericInput = input_trait(k))
    r = length(k.args)
    A = W.A # this is a LazyMatrixSum of LazyMatrixProducts
    prod_k_xy = k(x, y)
    if !iszero(prod_k_xy) # main case
        for i in 1:r # parallelize this?
            h, H = k.args[i], A.args[i]
            ui, vi = @views W.U[:, i], W.V[i, :]
            product_gradient_kernel_helper!(h, H, prod_k_xy, ui, vi, x, y)
        end
    else # requires less efficient O(dr²) special case if prod_k_xy is zero, for now, throw error
        throw(DomainError("GradientKernel of Product not supported for products that are zero, since it would require a O(dr²) special case."))
    end
    return W
end

@inline function product_gradient_kernel_helper!(ki, Ki, prod_k_xy, ui, vi, x, y, IT::InputTrait = input_trait(ki))
    ki_xy = ki(x, y)
    D = Ki.args[1]
    @. D.diag = prod_k_xy / ki_xy
    Ki.args[2] = gradient_kernel!(Ki.args[2], ki, x, y, IT)
    value_gradient_covariance!(ui, vi, ki, x, y, IT)
    @. ui *= prod_k_xy / ki_xy
    @. vi /= ki_xy
    return Ki
end

############################# Separable Product ################################
# for product kernel with generic input
function allocate_gradient_kernel(k::SeparableProduct, x::AbstractVector{<:Number},
                                  y::AbstractVector{<:Number}, ::GenericInput)
    d = length(x)
    H = (allocate_gradient_kernel(h, x, y, input_trait(h)) for h in k.args)
    T = typeof(k(x, y))
    A = LazyMatrixProduct(Diagonal(zeros(T, d)), Diagonal(zeros(T, d)))
    U = Diagonal(zeros(T, d))
    V = Diagonal(zeros(T, d))
    r = length(k.args)
    d == r || throw(DimensionMismatch("d = $d ≠ $r = r where r is number of product kernel constituents"))
    C = Woodbury((-1I)(r), ones(r), ones(r)')
    Woodbury(A, U, C, V)
end
function gradient_kernel(k::SeparableProduct, x, y, ::GenericInput)
    W = allocate_gradient_kernel(k, x, y, GenericInput())
    gradient_kernel!(W, k, x, y, GenericInput())
end

function gradient_kernel!(W::Woodbury, k::SeparableProduct, x::AbstractVector, y::AbstractVector, ::GenericInput = input_trait(k))
    A = W.A # this is a LazyMatrixProducts of Diagonals
    prod_k_xy = k(x, y)
    D, H = A.args # first is scaling matrix by leave_one_out_products, second is diagonal of derivative kernels
    for (i, ki) in enumerate(k.args)
        xi, yi = x[i], y[i]
        kixy = ki(xi, yi)
        # D[i, i] = prod_k_xy / kixy
        D[i, i] = ki(xi, yi)
        W.U[i, i] = ForwardDiff.derivative(z->ki(z, yi), xi)
        W.V[i, i] = ForwardDiff.derivative(z->ki(xi, z), yi)
        W.U[i, i] *= prod_k_xy / kixy
        W.V[i, i] /= kixy
        H[i, i] = DerivativeKernel(ki)(xi, yi)
    end
    leave_one_out_products!(D.diag)
    return W
end

############################# Separable Sum ####################################
# IDEA: implement block separable with x::AbstractVecOfVec
function allocate_gradient_kernel(k::SeparableSum, x::AbstractVector{<:Number},
                                  y::AbstractVector{<:Number}, ::GenericInput)
    f, h, d = k.f, k.k, length(x)
    H = allocate_gradient_kernel(h, x, y, input_trait(h))
    D = Diagonal(d)
end

function gradient_kernel!(D::Diagonal, k::SeparableSum, x::AbstractVector{<:Number},
                          y::AbstractVector{<:Number}, ::GenericInput)
    for (i, ki) in enumerate(k.args)
        D[i, i] = DerivativeKernel(ki)(x[i], y[i])
    end
    return D
end

############################## Input Transformations ###########################
# can be most efficiently represented by factoring out the Jacobian w.r.t. input warping
function gramian(G::GradientKernel{<:Real, <:Warped},  x::AbstractVector, y::AbstractVector)
    W = G.k
    U(x) = BlockFactorization(Diagonal([ForwardDiff.jacobian(W.u, xi) for xi in x]))
    k = GradientKernel(W.k)
    LazyMatrixProduct(U(x)', gramian(k, x, y), U(y))
end

function gramian(G::GradientKernel{<:Real, <:ScaledInputKernel},  x::AbstractVector, y::AbstractVector)
    n, m = length(x), length(y)
    S = G.k
    Ux = kronecker(I(n), S.U)
    Uy = n == m ? Ux : kronecker(I(m), S.U)
    k = GradientKernel(S.k)
    LazyMatrixProduct(Ux', gramian(k, x, y), Uy)
end

# I don't think this needs a special case, since we can take care of it in
# function gramian(G::GradientKernel{<:Real, <:Lengthscale}, x::AbstractVector, y::AbstractVector)
#     n, m = length(x), length(y)
#     L = G.k
#     Ux = Diagonal(fill(L.l, d*n)) # IDEA: Fill for lazy uniform array
#     Uy = n == m ? Ux : Diagonal(fill(L.l, d*m))
#     k = GradientKernel(L.k)
#     LazyMatrixProduct(Ux', gramian(k, x, y), Uy)
# end

############################### VerticalRescaling ##############################
# gradient element can be expressed with WoodburyIdentity and LazyMatrixProduct
function allocate_gradient_kernel(k::VerticalRescaling, x, y, ::GenericInput = GenericInput())
    f, h, d = k.f, k.k, length(x)
    H = allocate_gradient_kernel(h, x, y, input_trait(h))
    A = LazyMatrixProduct(Diagonal(fill(f(x), d)), H, Diagonal(fill(f(y), d)))
    U = zeros(d, 2)
    V = zeros(d, 2)
    C = zeros(2, 2)
    Woodbury(A, U, C, V')
end

function gradient_kernel!(W::Woodbury, k::VerticalRescaling, x, y, ::GenericInput = GenericInput())
    f, h, A = k.f, k.k, W.A
    fx, fy = f(x), f(y)
    @. A.args[1].diag = fx
    H = A.args[2] # LazyMatrixProduct: first and third are the diagonal scaling matrices, second is the gradient_kernel_matrix of h
    @. A.args[3].diag = fy
    A.args[2] = gradient_kernel!(H, h, x, y, input_trait(h))
    ForwardDiff.gradient!(@view(W.U[:, 1]), f, x)
    ForwardDiff.gradient!(@view(W.U[:, 2]), z->h(z, y), x)
    ForwardDiff.gradient!(@view(W.V[1, :]), f, y)
    ForwardDiff.gradient!(@view(W.V[2, :]), z->h(x, z), y)
    W.C[1, 1] = h(x, y)
    W.C[1, 2] = fy
    W.C[2, 1] = fx
    return W
end

############################ Scalar Chain Rule #################################
# generic implementation of scalar chain rule, does not require input kernel to have a basic input type
# gradient element can be expressed with WoodburyIdentity and LazyMatrixProduct
function allocate_gradient_kernel(k::Chained, x, y, ::GenericInput)
    f, h, d = k.f, k.k, length(x)
    H = allocate_gradient_kernel(h, x, y, input_trait(h))
    A = LazyMatrixProduct(Diagonal(fill(f(h(x, y)), d)), H)
    U = zeros(d, 1)
    V = zeros(d, 1)
    C = zeros(1, 1)
    Woodbury(A, U, C, V')
end

function gradient_kernel!(W::Woodbury, k::Chained, x, y, ::GenericInput)
    f, h, A = k.f, k.k, W.A
    f1, f2 = derivative_laplacian(f, h(x, y))
    @. A.args[1].diag = f1
    H = A.args[2] # LazyMatrixProduct: first argument is diagonal scaling, second is the gradient_kernel_matrix of h
    H.A.args[2] = gradient_kernel!(H, h, x, y, input_trait(h))
    ForwardDiff.gradient!(@view(W.U[:]), z->h(z, y), x)
    ForwardDiff.gradient!(@view(W.V[1, :]), z->h(x, z), y)
    @. W.C = f2
    return W
end


text:
@testset "xkcd" begin
	xkcd = load_xkcd()
	@test length(xkcd) == 949
	@test xkcd["purple"] == RGB24(0.494,0.118,0.612)
end

@testset "Resene" begin
	resene = load_resene()
	@test length(resene) == 1383
	@test resene["Cold Purple"] == reinterpret(RGB24, 0xaba0d9)
end

@testset "NBS" begin
	data = load_nbs()
	@test length(data) == 267
	@test data["moderatepurple"] == reinterpret(RGB24, 0x86608e)
end

@testset "X11" begin
	data = load_x11()
	@test length(data) > 200
	@test data["purple"] == reinterpret(RGB24, 0xa020f0)
end

@testset "CSS3" begin
	data = load_css3()
	@test length(data) == 16
	@test data["purple"] == reinterpret(RGB24, 0x800080)
end

@testset "Crayola" begin
	data = load_crayola()
	@test length(data) == 195 #Note, there seems to be a duplicate
	@test data["Maximum Purple"] == reinterpret(RGB24, 0x733380)
end

@testset "WinsorNewton" begin
	data = load_winsor_newton()
	@test length(data) == 600
	@test data["purple lake light"] == reinterpret(RGB24, 0x784D63)
end

ext: jl


text:
using OMEinsumContractionOrders, OMEinsum
using Test, Random
using OMEinsum: getixs, getiy, isleaf

function have_identity(ne::NestedEinsum)
    if isleaf(ne)
        return false
    elseif length(getixs(ne.eins)) == 1 && getixs(ne.eins)[1] == getiy(ne.eins)
        return true
    else
        return any(have_identity, ne.args)
    end
end

@testset "simplify vectors" begin
    tn = ein"ab,bc,cd,de,a,b,c,d,f,f,e->ab"
    xs = [randn(fill(2, length(ix))...) for ix in OMEinsum.getixs(tn)]
    simplifier, tn2 = merge_vectors(tn)
    @test tn2 == ein"ab,bc,cd,de,f->ab"
    @test tn(xs...) ≈ tn2(simplifier(xs...)...)
    tn3 = optimize_greedy(tn2, uniformsize(tn2, 2))
    tn4 = embed_simplifier(tn3, simplifier)
    @test !have_identity(tn4)
    @test tn(xs...) ≈ tn4(xs...)

    tn = OMEinsum.DynamicEinCode(tn)
    xs = [randn(fill(2, length(ix))...) for ix in OMEinsum.getixs(tn)]
    simplifier, tn2 = merge_vectors(tn)
    @test tn2 == OMEinsum.DynamicEinCode(ein"ab,bc,cd,de,f->ab")
    @test tn(xs...) ≈ tn2(simplifier(xs...)...)
    tn3 = optimize_greedy(tn2, uniformsize(tn2, 2))
    tn4 = embed_simplifier(tn3, simplifier)
    @test !have_identity(tn4)
    @test tn(xs...) ≈ tn4(xs...)
end

@testset "simplify greedy" begin
    tn = ein"ab,bc,cd,de,a,b,c,d,f,f,e->ab"
    size_dict = uniformsize(tn, 2)
    xs = [randn(fill(2, length(ix))...) for ix in OMEinsum.getixs(tn)]
    simplifier, tn2 = merge_greedy(tn, size_dict)
    @test tn2 == ein"ab,->ab"
    @test tn(xs...) ≈ tn2(simplifier(xs...)...)
    tn3 = optimize_greedy(tn2, uniformsize(tn2, 2))
    tn4 = embed_simplifier(tn3, simplifier)
    @test !have_identity(tn4)
    @test tn(xs...) ≈ tn4(xs...)

    tn = OMEinsum.DynamicEinCode(tn)
    xs = [randn(fill(2, length(ix))...) for ix in OMEinsum.getixs(tn)]
    simplifier, tn2 = merge_greedy(tn, size_dict)
    @test tn2 == OMEinsum.DynamicEinCode(ein"ab,->ab")
    @test tn(xs...) ≈ tn2(simplifier(xs...)...)
    tn3 = optimize_greedy(tn2, uniformsize(tn2, 2))
    tn4 = embed_simplifier(tn3, simplifier)
    @test !have_identity(tn4)
    @test tn(xs...) ≈ tn4(xs...)
end

@testset "optimize permute" begin
    code = ein"lcij,lcjk->kcli"
    @test OMEinsumContractionOrders.optimize_output_permute(getixsv(code), getiyv(code)) == collect("iklc")
    code = ein"lcij,lcjk,c->kcli"
    @test OMEinsumContractionOrders.optimize_output_permute(getixsv(code), getiyv(code)) == collect("kcli")
    code = ein"lcij,lcjk->"
    @test OMEinsumContractionOrders.optimize_output_permute(getixsv(code), getiyv(code)) == Char[]
end

@testset "optimize permute" begin
    code = ein"(lcij,lcjk),lcik->"
    res = NestedEinsum((NestedEinsum((NestedEinsum{StaticEinCode}(1),NestedEinsum{StaticEinCode}(2)), ein"lcij,lcjk->iklc"), NestedEinsum{StaticEinCode}(3)), ein"iklc,lcik->")
    @test optimize_permute(code) == res
    code = SlicedEinsum(Slicing(['c']), ein"(lcij,lcjk),lcik->")
    @test optimize_permute(code) == SlicedEinsum(Slicing(['c']), res)
end


text:
module GNSSSignals

    using Core: toInt16
    using DocStringExtensions
    using FixedPointNumbers
    using Statistics
    using Unitful: Frequency, Hz

    using CUDA
    const use_gpu = Ref(false)

    import Base.show

    export
        AbstractGNSS,
        GPSL1,
        GPSL5,
        GalileoE1B,
        BOCcos,
        gen_code!,
        gen_code,
        get_codes,
        get_code_length,
        get_secondary_code_length,
        get_center_frequency,
        get_code_frequency,
        get_code,
        get_code_unsafe,
        get_data_frequency,
        get_code_center_frequency_ratio,
        get_subcarrier_frequency,
        get_code_spectrum,
        get_system_string,
        min_bits_for_code_length

    abstract type AbstractGNSS{C} end
    abstract type AbstractGNSSBOCcos{C, M, N} <: AbstractGNSS{C} end

    Base.Broadcast.broadcastable(system::AbstractGNSS) = Ref(system)

    Base.show(io::IO, x::AbstractGNSS) = print("$(typeof(x))()")

    """
    $(SIGNATURES)

    `GNSSSignals.jl` checks if there is a working installation of CUDA on the system and informs the user
    to activate GPU acceleration if they wish to do so.
    """
    function __init__()
        if CUDA.functional()
            @info "Found a working CUDA installation. To activate GPU acceleration set use_gpu = Val(true), e.g. GPSL1(use_gpu = Val(true))."
        end
    end

    """
    $(SIGNATURES)

    Reads Int8 encoded codes from a file with filename `filename` (including the path). The
    code length must be provided by `code_length` and the number of PRNs by `num_prns`.
    # Examples
    ```julia-repl
    julia> read_in_codes("/data/gpsl1codes.bin", 32, 1023)
    ```
    """
    function read_in_codes(type, filename, num_prns, code_length)
        open(filename) do file_stream
            read!(file_stream, Array{type}(undef, code_length, num_prns))
        end
    end

    include("gps_l1.jl")
    include("gps_l5.jl")
    include("galileo_e1b.jl")
    include("boc.jl")
    include("common.jl")
end


text:
type LaplacianModel{C <: CoordinateSystem} <: GP
	k::Kernel
	lik::Likelihood
	coordinates::C
	jmax::Integer
	js::Array{Int64, 1}
	X::Array{Float64, 2}
	y::Array{Float64, 1}
	y_raw::Array{Float64, 1}
	stale::Bool
	Phi::Array{Float64, 2}
	lambdas::Array{Float64, 1}
	PhiTPhi::Array{Float64, 2}
	PhiTy::Array{Float64, 1}
	yTy::Float64
	Z::Array{Float64, 2}
	Zchol::Base.LinAlg.Cholesky{Float64}
	mu::Float64
	std::Float64
end

LaplacianModel{C <: CoordinateSystem}(k::Kernel, lik::Likelihood, c::C, jmax::Integer) = LaplacianModel(k, lik, c, jmax, Array(Int64, (0, )),
		Array(Float64, (0, 0)), Array(Float64, (0, )), Array(Float64, (0, )),
		true,
		Array(Float64, (0, 0)), Array(Float64, (0, )),
		Array(Float64, (0, 0)), Array(Float64, (0,)), 0.0,
		Array(Float64, (0, 0)), cholfact(eye(1)),
		0.0, 1.0)

function observe!{T <: FloatingPoint}(m::LaplacianModel, X::Array{T}, y::Array{T})
	if size(m.X) == (0,0)
		m.X = X
		m.y_raw = y
	else
		m.X = vcat(m.X, X)
		m.y_raw = vcat(m.y_raw, y)
	end
	m.mu = mean(m.y_raw)
	m.std = std(m.y_raw)
	m.y = (m.y_raw .- m.mu) / (m.std + eps())
	refresh_cache!(m)
	m.stale = true
end

function predict(m::LaplacianModel, Xp::Array{Float64, 2})
	if m.stale
		refresh_cache!(m)
	end
	sigma2 = exp(2m.lik.log_sigma)
	Ssqrtλ = Float64[ S( m.k, √(λ) ) for λ in m.lambdas ]
	Λinv = eye(diagm(Ssqrtλ)) / diagm(Ssqrtλ)
	size(m.js, 1)
	phi_star = Float64[φ(m.coordinates, Xp[r, :][:], m.js[jr]) for r = 1:size(Xp, 1), jr = 1:size(m.js, 1)]
	Z = (m.PhiTPhi + sigma2 * Λinv)
	size(Z), size(phi_star), size(m.PhiTy)
	mu_p = vec(phi_star * (Z \ m.PhiTy)) * m.std .+ m.mu
	sigma_p = m.std^2 * vec(sigma2 * sum(phi_star' .* (Z \ phi_star'), 1))
	mu_p, sigma_p
end

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

function refresh_cache!(m::LaplacianModel)
	m.js = get_js(m)
	m.lambdas = [λ(m.coordinates,  m.js[jr]) for jr = 1:size(m.js, 1)]

	m.Phi = [φ(m.coordinates, m.X[r, :][:], m.js[jr]) for r = 1:size(m.X, 1), jr = 1:size(m.js, 1)]

	m.PhiTPhi = m.Phi'm.Phi
	m.PhiTy = m.Phi'*m.y
	m.yTy = dot(m.y,m.y)
	m.stale = false
end

function get_js{C <: Cartesian}(m::LaplacianModel{C})
	Int64[int(mod(floor(i/(m.jmax^(d-1))) , m.jmax)) for i = 0:((m.jmax^m.coordinates.dimensions)-1), d = 1:m.coordinates.dimensions] .+ 1
end

function get_js{C <: OneD}(m::LaplacianModel{C})
	[1:m.jmax]
end


text:
struct VariableSet{T <: AbstractVariable}
    vars::Vector{T}
    nmes
end

function VariableSet(vars::AbstractVariable...)
    nmes = Dict([Symbol(i) => i for i ∈ eachindex(vars)])
    VariableSet([vars...], nmes)
end

function VariableSet(vars::Pair{Symbol,}...)
    nmes = Dict([first(vars[i]) => i for i ∈ eachindex(vars)])
    vars = [last(vars[i]) for i ∈ eachindex(vars)]
    VariableSet(vars, nmes)
end

function Base.names(vs)
    nm = collect(keys(vs.nmes))
    idx = collect(values(vs.nmes))
    nn = copy(nm)
    for i ∈ eachindex(idx)
        nn[idx[i]] = nm[i]
    end
    tuple(nn...)
end

Base.getindex(vs::VariableSet, i::Int) = vs.vars[i]

Base.getindex(vs::VariableSet, n) = vs.vars[vs.nmes[n]]

Base.setindex!(vs::VariableSet, w, i::Int) = (vs.vars[i] = w)

Base.setindex!(vs::VariableSet, w, n) = (vs.vars[vs.nmes[n]] = w)

Base.length(vs::VariableSet) = length(vs.vars)

Base.firstindex(vs::VariableSet) = 1

Base.lastindex(vs::VariableSet) = length(vs)

Base.iterate(vs::VariableSet) = iterate(vs, 1)

function Base.iterate(vs::VariableSet, st)
    st > length(vs) && return nothing
    (vs[st], st + 1)
end

function Base.push!(vs::VariableSet, var::Pair)
    push!(vs.vars, last(var))
    n = length(vs)
    push!(vs.nmes, first(var) => n)
    vs
end

function Base.push!(vs::VariableSet, var::AbstractVariable)
    push!(vs.vars, var)
    n = length(vs)
    push!(vs.nmes, Symbol(n) => n)
    vs
end

function Base.push!(vs::VariableSet, vars::AbstractVariable...)
    for v ∈ vars
        push!(vs.vars, v)
    end
    vs
end


text:
struct ParameterGrid2D
    r_range::StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}
    y₀_range::StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}
    z₀_range::StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}
end
ParameterGrid2D(Ly::Number, Lz::Number, Ncourse::Int) = begin
    ParameterGrid2D(range(0,Ly/2,length=Ncourse),
                    range(-Ly/2,Ly/2,length=Ncourse),
                    range(-Lz,0,length=Ncourse))
end
function grid_search(pgrid::ParameterGrid2D, bitmap_input, guess_cache, ynodes, znodes)
    i_loc, j_loc, k_loc = 1, 1, 1

    val_min = circ_error(bitmap_input, guess_cache, ynodes, znodes,
            pgrid.r_range[1], pgrid.y₀_range[1], pgrid.z₀_range[1])
    for (i, r) in enumerate(pgrid.r_range)
        for (j, y) in enumerate(pgrid.y₀_range)
            for (k, z) in enumerate(pgrid.z₀_range)
                err = circ_error(bitmap_input, guess_cache, ynodes, znodes, r, y, z)
                if err < val_min
                    val_min = err
                    i_loc, j_loc, k_loc = i, j, k
                end
            end
        end
    end
    return (i_loc, j_loc, k_loc), val_min
end
function coarse_fine_search(pgrid::ParameterGrid2D, bitmap_input, guess_cache, ynodes, znodes)
    # Coarse search
    res_course = grid_search(pgrid, bitmap_input, guess_cache, ynodes, znodes)

    # Fine search
    to_minim(x) = circ_error(bitmap_input, guess_cache, ynodes, znodes, x[1], x[2], x[3])
    x0 = [pgrid.r_range[res_course[1][1]],
          pgrid.y₀_range[res_course[1][2]],
          pgrid.z₀_range[res_course[1][3]]]
    res = optimize(to_minim, x0)

    return res.minimizer, res.minimum, x0
end


text:
function nnls(A::AbstractMatrix, b::AbstractVector, ops::Options)
    # Set up the preconditioner as the inverse diagonal of the Hessian
    a = sum(A.^2, 1)
    @defaults ops precondprep=(P, x, l, u, mu)->precondprepnnls(P, x, mu, a)
    ops = copy(ops)
    @set_options ops precondprep=precondprep
    # Create the initial guess (an interior point)
    T = promote_type(eltype(A), eltype(b))
    x = fill(one(T), size(A, 2))
    # Set up constraints
    l = zeros(eltype(x), length(x))
    u = fill(inf(eltype(x)), length(x))
    # Perform the optimization
    func = (g,x) -> nnlsobjective(g, x, A, b)
    fminbox(func, x, l, u, ops)
end
nnls(A::AbstractMatrix, b::AbstractVector) = nnls(A, b, Options())

function nnlsobjective(g, x::AbstractVector, A::AbstractMatrix, b::AbstractVector)
    d = A*x - b
    val = sum(d.^2)/2
    if !(g === nothing)
        At_mul_B(g, A, d)
    end
    val
end

function precondprepnnls(P, x, mu, a)
    for i = 1:length(x)
        P[i] = 1/(mu/x[i]^2 + a[i])
    end
end

export nnls


text:
#   Unit linearAlgebra.jl, part of PosDefManifold Package for julia language
#   v 0.0.3 - last update 16th of May 2019
#
#   MIT License
#   Copyright (c) 2019, Marco Congedo, CNRS, Grenobe, France:
#   https://sites.google.com/site/marcocongedo/home
#
#   DESCRIPTION
#   This Unit implements Linear Algebra methods and procedures
#   useful in relation to the Riemannian Geometry on the manifold
#   of Symmetric Positive Definite (SPD) or Hermitian matrices.
#
#   CONTENT
#   0. Utilities
#   1. Matrix normalizations
#   2. Boolean functions of matrices
#   3. Scalar functions of matrices
#   4. Diagonal functions of matrices
#   5. Unitary functions of matrices
#   6. Matrix function of matrices
#   7. Spectral decompositions of positive matrices
#   8. Decompositions involving triangular matrices
# __________________________________________________________________

#  ------------------------
## 0. Utilities
#  ------------------------
"""
    function typeofMatrix(array::Union{ AnyMatrix,
                                        AnyMatrixVector,
                                        AnyMatrixVector₂ })

 **alias**: `typeofMat`

 Return the type of a matrix, either `Hermitian`,
 `Diagonal`, `LowerTriangular`, or `Matrix`.
 Argument `array` may be a matrix of one of these types, but also one of the following:

 ℍVector, ℍVector₂, 𝔻Vector, 𝔻Vector₂, 𝕃Vector, 𝕃Vector₂, 𝕄Vector, 𝕄Vector₂.

 Those are [Array of Matrices types](@ref).
 See also [aliases](@ref) for the symbols `ℍ`, `𝔻`, `𝕃` and `𝕄`.

 Note that this function is different from Julia function
 [typeof](https://docs.julialang.org/en/v1/base/base/#Core.typeof),
 which returns the concrete type (see example below), thus
 cannot be used for [typecasting matrices](@ref).

 ## Examples
    using LinearAlgebra, PosDefManifold
    P=randP(3) # generate a 3x3 Hermitian matrix
    typeofMatrix(P) # returns `Hermitian`
    typeof(P) # returns `Hermitian{Float64,Array{Float64,2}}`
    # typecast P as a `Matrix` M
    M=Matrix(P)
    # typecast M as a matrix of the same type as P and write the result in A
    A=typeofMatrix(P)(M)

    Pset=randP(3, 4) # generate a set of 4 3x3 Hermitian matrix
    # Pset is an ℍVector type
    typeofMatrix(Pset) # again returns `Hermitian`

"""
typeofMatrix(H::Union{ℍ, ℍVector, ℍVector₂}) = ℍ
typeofMatrix(D::Union{𝔻, 𝔻Vector, 𝔻Vector₂}) = 𝔻
typeofMatrix(L::Union{𝕃, 𝕃Vector, 𝕃Vector₂}) = 𝕃
typeofMatrix(M::Union{𝕄, 𝕄Vector, 𝕄Vector₂}) = 𝕄
typeofMat=typeofMatrix

"""
    (1) function dim(X::AnyMatrix, [d])
    (2) function dim(vector::AnyMatrixVector, [d])
    (3) function dim(vector₂::AnyMatrixVector₂, [d])

 (1) ``X`` is a real or complex `Matrix`, `Diagonal`,
 `LowerTriangular` or `Hermitian` matrix.
 Return a 2-tuple containing the dimensions of ``X``,
 which is two times the same dimension for all possible types of ``X``
 with the exception of the `Matrix` type, which can be rectangular.
 Optionally you can specify a dimension (1 or 2)
 to get just the length of that dimension.

 (2) `vector` is an 𝕄Vector, 𝔻Vector, 𝕃Vector or ℍVector type
 (see [AnyMatrixVector type](@ref)).
 Return a 3-tuple containing the number of matrices it holds
 (dimension 1) and their dimensions (dimension 2 and 3).
 Optionally you can specify a dimension (1, 2, or 3)
 to get just the length of that dimension.

 (3) `vector₂` is an 𝕄Vector₂, 𝔻Vector₂, 𝕃Vector₂ or ℍVector₂ type
 (see [AnyMatrixVector type](@ref)).
 Return a 4-tuple containing
 - the number of vectors of matrices it holds (dimension 1),
 - a vector holding the number of matrices in each vector of matrices (dimensions 2),
 - the two dimensions of the matrices (dimension 3 and 4).
 Optionally you can specify a dimension (1, 2, 3 or 4)
 to get just the length of that dimension.

 `vector` and `vector₂` are [Array of Matrices types](@ref).
 See also [aliases](@ref) for the symbols `ℍ`, `𝔻`, `𝕃` and `𝕄`.

!!! note "nota bene"
    If you specify a dimension and this is out of the valid range,
    the function returns zero.

    Both the `vector`(2) and the `vector₂`(3) object are meant to hold
    matrices living in the same manifold, therefore it is assumed
    that all matrices they holds are of the same dimension.
    The dimensions of the matrices are retrived from
    - the first matrix in `vector`(2),
    - the first matrix in the first vector of `vector₂`(3).

 This function replaces Julia [size](https://docs.julialang.org/en/v1/base/arrays/#Base.size)
 function, which cannot be used to retrive dimension for matrix vectors.
 It is not possible to overload the `size` function for matrix vectors
 since this causes problems to other Julia functions.

 ## Examples
    using LinearAlgebra, PosDefManifold
    # (1)
    M=randn(3, 4) # generate a 3x4 `Matrix`
    dim(M) # returns (3, 4)
    dim(M, 1) # returns 3
    dim(M, 2) # returns 4
    dim(M, 3) # out of range: returns 0

    # (2)
    Pset=randP(3, 4) # generate an ℍVector holding 4 3x3 Hermitian matrices
    dim(Pset) # returns (4, 3, 3)
    dim(Pset, 1) # returns 4
    dim(Pset, 2) # returns 3
    dim(Pset, 3) # returns 3

    # (3)
    # Generate a set of 4 random 3x3 SPD matrices
    Pset=randP(3, 4)
    # Generate a set of 40 random 4x4 SPD matrices
    Qset=randP(3, 40)
    A=ℍVector₂([Pset, Qset])
    dim(A) # return (2, [4, 40], 3, 3)
    dim(A, 1) # return 2
    dim(A, 2) # return [4, 40]
    dim(A, 2)[1] # return 4
    dim(A, 3) # return 3
    dim(A, 4) # return 3
    dim(A, 5) # out of range: return 0

    # note: to create an ℍVector₂ object holding k ℍVector objects use
    sets=ℍVector₂(undef, k) # and then fill them
"""
dim(X::AnyMatrix, d::Int) = 1<=d<=2 ? size(X, d) : 0
dim(X::AnyMatrix) = size(X)

function dim(vector::AnyMatrixVector, d::Int) # change to bold X
    if      d==1    return length(vector)
    elseif  2<=d<=3 return size(vector[1], d-1)
    elseif          return 0
    end
end
dim(vector::AnyMatrixVector) = (length(vector), size(vector[1], 1), size(vector[1], 2))

function dim(vector₂::AnyMatrixVector₂, d::Int)
    if      d==1    return length(vector₂)
    elseif  d==2    return collect(length(vec) for vec in vector₂)
    elseif  3<=d<=4 return size(vector₂[1][1], d-2)
    elseif          return 0
    end
end
dim(vector₂::AnyMatrixVector₂) =
    (length(vector₂), collect(length(vec) for vec in vector₂), size(vector₂[1][1], 1), size(vector₂[1][1], 2))

#  ------------------------
## 1. Matrix Normalizations
#  ------------------------

"""
    function det1(X::AnyMatrix; <tol::Real=0>)

 Return the argument matrix ``X`` normalized so as to have *unit determinant*.
 For square positive definite matrices this is the best approximant
 from the set of matrices in the [special linear group](https://bit.ly/2W5jDZ6) -
 see Bhatia and Jain (2014)[🎓].

 ``X`` can be a real or complex `Diagonal`, `LowerTriangular`,
 `Matrix`, or `Hermitian` matrix. (see [AnyMatrix type](@ref))

 If the determinant is not greater to `tol` (which defalts to zero)
 a warning is printed and ``X`` is returned.

!!! note "Nota Bene"
    This function is meant for positive definite matrices.
    Julia may throws an error while computing
    the determinant if the matrix is defective.

 **See** [Julia det function](https://bit.ly/2Y4MnTF).

 **See also**: [`tr1`](@ref).

 ## Examples
    using LinearAlgebra, PosDefManifold
    P=randP(5) # generate a random real positive definite matrix 5x5
    Q=det1(P)
    det(Q) # must be 1
    # using a tolerance
    Q=det1(P; tol=1e-12)
"""
function det1(X::AnyMatrix; tol::Real=0)
    tol>=0 ? tolerance=tol : tolerance = 0
    Det = det(X)
    if Det>tolerance X/Det^(1/size(X, 1)) else @warn det1Msg Det tolerance; X end
end
det1Msg="function det1 in LinearAlgebra.jl of PosDefMaifold package: the determinant of the input matrix is not greater than the tolerance."

"""
    tr1(X::AnyMatrix)

 Return the argument matrix ``X`` normalized so as to have *unit trace*.

 ``X`` can be a real or complex `Diagonal`, `LowerTriangular`,
 `Matrix` or `Hermitian` matrix (see [AnyMatrix type](@ref)).

 If the trace is not greater to `tol`
 (which defalts to zero) a warning is printed and ``X`` is returned.

 **See**: [Julia trace function](https://bit.ly/2HoOLiM).

 **See also**: [`tr`](@ref), [`det1`](@ref).

 ## Examples
    using LinearAlgebra, PosDefManifold
    P=randP(5) # generate a random real positive definite matrix 5x5
    Q=tr1(P)
    tr(Q)  # must be 1
    # using a tolerance
    Q=tr1(P; tol=1e-12)

"""
function tr1(X::AnyMatrix; tol::Real=0)
    tol>=0 ? tolerance=tol : tolerance = 0
    trace = tr(X)
    if trace>tolerance X/trace else @warn tr1Msg trace tolerance; X end
end
tr1Msg="function tr1 in LinearAlgebra.jl of PosDefMaifold package: the trace of the input matrix is not greater than the tolerance."

"""
    (1) normalizeCol!(X::𝕄{T}, j::Int)
    (2) normalizeCol!(X::𝕄{T}, j::Int, by::Number)
    (3) normalizeCol!(X::𝕄{T}, range::UnitRange)
    (4) normalizeCol!(X::𝕄{T}, range::UnitRange, by::Number)
                     for all above: where T<:RealOrComplex

 Given a `Matrix` type ``X`` comprised of real or complex elements,
 - (1) normalize the ``j^{th}`` column to unit norm
 - (2) divide the elements of the ``j^{th}`` column by number ``by``
 - (3) normalize the columns in ``range`` to unit norm
 - (4) divide the elements of columns in ``range``  by number ``by``.

 ``by`` is a number of abstract supertype [Number](https://bit.ly/2JwXjGr).
 It should be an integer, real or complex number.
 For efficiency, it should be of the same type as the elements of ``X``.

 ``range`` is a [UnitRange](https://bit.ly/2HSfK5J) type.

 Methods (1) and (3) call the
 [BLAS.nrm2](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.BLAS.nrm2)
 routine for computing the norm of concerned columns.
 See [Threads](@ref).

!!! note "Nota Bene"
    Julia does not allow normalizing the columns of `Hermitian` matrices.
    If you want to call this function for an `Hermitian` matrix see [typecasting matrices](@ref).

 **See** [norm](https://bit.ly/2TaAkR0) and also [randn](https://bit.ly/2I1Vgrg)
 for the example below.

 **See also**: [`colNorm`](@ref), [`colProd`](@ref).

 ## Examples
    using PosDefManifold
    X=randn(10, 20)
    normalizeCol!(X, 2)                  # (1) normalize columns 2
    normalizeCol!(X, 2, 10.0)            # (2) divide columns 2 by 10.0
    normalizeCol!(X, 2:4)                # (3) normalize columns 2 to 4
    X=randn(ComplexF64, 10, 20)
    normalizeCol!(X, 3)                  # (1) normalize columns 3
    normalizeCol!(X, 3:6, (2.0 + 0.5im)) # (4) divide columns 3 to 5 by (2.0 + 0.5im)

"""
function normalizeCol!(X::𝕄{T}, j::Int) where T<:RealOrComplex
    w=colNorm(X, j)
    for i=1:size(X, 1) @inbounds X[i, j]/=w end
end

normalizeCol!(X::𝕄{T}, j::Int, by::Number) where T<:RealOrComplex =
             for i=1:size(X, 1) @inbounds X[i, j]/=by end

function normalizeCol!(X::𝕄{T}, range::UnitRange) where T<:RealOrComplex
    l=range[1]-1
    w=[colNorm(X, j) for j in range]
    for j in range normalizeCol!(X, j, w[j-l]) end
end

normalizeCol!(X::𝕄{T}, range::UnitRange, by::Number) where T<:RealOrComplex =
             for j in range normalizeCol!(X, j, by) end

#  -------------------------------
## 2. Boolean Functions of Matrices
#  -------------------------------

"""
```
    (1) ispos(   λ::Vector{T};
                <tol::Real=0, rev=true, 🔔=true, msg="">) where T<:Real

    (2) ispos(   Λ::𝔻{T};
                <tol::Real=0, rev=true, 🔔=true, msg="">) where T<:Real
```

 Return ``true`` if all numbers in (1) real vector ``λ`` or in (2) real `Diagonal`
 matrix ``Λ`` are not inferior to ``tol``, otherwise return ``false``. This is used,
 for example, in spectral functions to check that all eigenvalues are positive.

!!! note "Nota Bene"
    ``tol`` defaults to the square root of `Base.eps` of the type of ``λ`` (1)
     or ``Λ`` (2). This corresponds to requiring positivity beyond about half of
     the significant digits.

 The following are *<optional keyword arguments>*:

 - If ``rev=true`` the (1) elements in ``λ`` or (2) the diagonal elements in ``Λ`` will be chacked in reverse order.
 This is done for allowing a very fast check when the elements
 are sorted and it is known from where is best to start checking.

 If the result is ``false``:
 - if ``🔔=true`` a bell character will be printed. In most systems this will ring a bell on the computer.
 - if string ``msg`` is provided, a warning will print ``msg`` followed by:
 "at position *pos*", where *pos* is the position where the
 first non-positive element has been found.

```
 ## Examples
 using PosDefManifold
 a=[1, 0, 2, 8]
 ispos(a, msg="non-positive element found")

 # it will print:
 # ┌ Warning: non-positive element found at position 2
 # └ @ [here julie will point to the line of code issuing the warning]
```
 """
function ispos(λ::Vector{T}; tol::Real=0, rev=true, 🔔=true, msg="") where T<:Real
    tol==0 ? tolerance = √eps(T) : tolerance = tol
    rev ? iterations = (length(λ):-1:1) : iterations=(1:length(λ))
    for i in iterations
        if λ[i]<tolerance
            🔔 && print('\a') # print('\a') sounds a bell
            length(msg)>0 && @warn("function ispos(linearAlgebra.jl) "*msg* " at position $i")
            return false; break
        end
    end
    return true
end

ispos(Λ::Diagonal{T}; tol::Real=0, rev=true, 🔔=true, msg="") where T<:Real =
      ispos( diag(Λ); tol=tol, rev=rev, 🔔=🔔, msg=msg)

#  -------------------------------
## 3. Scalar functions of matrices
#  -------------------------------

"""
    (1) colProd(X::Union{𝕄{T}, ℍ{T}}, j::Int, l::Int)
    (2) colProd(X::Union{𝕄{T}, ℍ{T}}, Y::Union{𝕄{T}, ℍ{T}}, j::Int, l::Int)
               for all above: where T<:RealOrComplex

 (1) Given a real or complex `Matrix` or `Hermitian` matrix ``X``,
 return the dot product of the ``j^{th}`` and ``l^{th}`` columns, defined as,

 ``\\sum_{i=1}^{r} \\big(x_{ij}^*x_{il}\\big), ``

 where ``r`` is the number of rows of ``X`` and ``^*`` denotes complex
 conjugate (nothing if the matrix is real).

 (2) Given real or complex `Matrix` or `Hermitian` matrices ``X`` and ``Y``,
 return the dot product of the ``j^{th}`` column of ``X`` and the ``l^{th}`` column
 of ``Y``, defined as,

 ``\\sum_{i=1}^{r} \\big(x_{ij}^*y_{il}\\big), ``

 where ``r`` is the number of rows of ``X`` and of ``Y`` and ``^*`` is as above.

!!! note "Nota Bene"
    ``X`` and of ``Y`` may have a different number of columns, but must have
    the same number of rows.

 Arguments ``j`` and ``l`` must be positive integers in range
 - (1) `j,l in 1:size(X, 2)`,
 - (2) `j in 1:size(X, 2), l in 1:size(Y, 2)`.

 **See also**: [`normalizeCol!`](@ref), [`colNorm`](@ref).

 ## Examples
    using PosDefManifold
    X=randn(10, 20)
    p=colProd(X, 1, 3)
    Y=randn(10, 30)
    q=colProd(X, Y, 2, 25)

"""
colProd(X::Union{𝕄{T}, ℍ{T}}, j::Int, l::Int) where T<:RealOrComplex =
        𝚺(conj(x1)*x2 for (x1, x2) in zip(X[:, j], X[:, l]))

colProd(X::Union{𝕄{T}, ℍ{T}}, Y::Union{𝕄{T}, ℍ{T}}, j::Int, l::Int) where T<:RealOrComplex =
        𝚺(conj(x1)*x2 for (x1, x2) in zip(X[:, j], Y[:, l]))

"""
    colNorm(X::Union{𝕄{T}, ℍ{T}}, j::Int) where T<:RealOrComplex

 Given a real or complex `Matrix` or `Hermitian` matrix ``X``,
 return the Euclidean norm of its ``j^{th}`` column.

 This function calls the
 [BLAS.nrm2](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.BLAS.nrm2)
 routine. See [Threads](@ref).

 **See also**: [`normalizeCol!`](@ref), [`colProd`](@ref), [`sumOfSqr`](@ref).

 ## Examples
    using PosDefManifold
    X=randn(10, 20)
    normOfSecondColumn=colNorm(X, 2)

"""
colNorm(X::Union{𝕄{T}, ℍ{T}}, j::Int) where T<:RealOrComplex =
        BLAS.nrm2(size(X, 1), X[:, j], 1)

"""
    (1) sumOfSqr(A::Array)
    (2) sumOfSqr(H::ℍ{T})
    (3) sumOfSqr(L::𝕃{T})
    (4) sumOfSqr(D::𝔻{T})
    (5) sumOfSqr(X::Union{𝕄{T}, ℍ{T}}, j::Int)
    (6) sumOfSqr(X::Union{𝕄{T}, ℍ{T}}, range::UnitRange)
                 for (1)-(6) above: where T<:RealOrComplex

**alias**: `ss`

 Return
 - (1) the sum of squares of the elements in an array ``A`` of any dimensions.
 - (2) as in (1), but for an `Hermitian` matrix ``H``, using only the lower triangular part.
 - (3) as in (1), but for a `LowerTriangular` matrix ``L``.
 - (4) as in (1), but for a `Diagonal` matrix ``D`` (sum of squares of diagonal elements).
 - (5) the sum of square of the ``j^{th}`` column of a `Matrix` or `Hermitian` ``X``.
 - (6) the sum of square of the columns of a `Matrix` or `Hermitian` ``X`` in a given range.

 All methods support real and complex matrices.

 Only method (1) works for arrays of any dimensions.

 Methods (1)-(4) return the square of the [Frobenius norm](https://bit.ly/2Fi10eH).

 For method (5), ``j`` is a positive integer in range `1:size(X, 1)`.

 For method (6), ``range`` is a [UnitRange type](https://bit.ly/2HDoFbk).

 **See also**: [`colNorm`](@ref), [`sumOfSqrDiag`](@ref), [`sumOfSqrTril`](@ref).

 ## Examples
    using PosDefManifold
    X=randn(10, 20)
    sum2=sumOfSqr(X)        # (1) sum of squares of all elements
    sum2=sumOfSqr(X, 1)     # (2) sum of squares of elements in column 1
    sum2=sumOfSqr(X, 2:4)   # (3) sum of squares of elements in column 2 to 4

"""
sumOfSqr(A::Array) = 𝚺(abs2(a) for a in A)

function sumOfSqr(H::ℍ{T}) where T<:RealOrComplex
    r=size(H, 1)
    s=real(T)(0)
    for j=1:size(H, 2)-1
        @inbounds s+=abs2(H[j, j])
        for i=j+1:r @inbounds s+=2*abs2(H[i, j]) end
    end
    @inbounds s+=abs2(H[r, r])
    return s
end

function sumOfSqr(L::𝕃{T}) where T<:RealOrComplex
    s=real(T)(0)
    for j=1:size(L, 2), i=j:size(L, 1)
        @inbounds s+=abs2(L[i, j])
    end
    return s
end

sumOfSqr(D::𝔻{T}) where T<:RealOrComplex = sumOfSqrDiag(D)

sumOfSqr(X::Union{𝕄{T}, ℍ{T}}, j::Int) where T<:RealOrComplex = 𝚺(abs2.(X[:, j]))

sumOfSqr(X::Union{𝕄{T}, ℍ{T}}, range::UnitRange) where T<:RealOrComplex =
         𝚺(sumOfSqr(X, j) for j in range)

ss=sumOfSqr

"""
    sumOfSqrDiag(X::AnyMatrix)

 **alias**: `ssd`

 Sum of squares of the diagonal elements in real or complex `Matrix`,
 `Diagonal`, `Hermitian` or `LowerTriangular` matrix ``X``.
 If ``X`` is rectangular (which can be only if it is of the `Matrix` type),
 the main diagonal is considered.

 **See** [AnyMatrix type](@ref)

 **See also**: [`sumOfSqr`](@ref), [`sumOfSqrTril`](@ref).

 ## Examples
    using LinearAlgebra, PosDefManifold
    X=randn(10, 20)
    sumDiag2=sumOfSqrDiag(X) # (1)
    sumDiag2=sumOfSqrDiag(𝔻(X)) # (2) 𝔻=LinearAlgebra.Diagonal

"""
sumOfSqrDiag(X::𝕄{T}) where T<:RealOrComplex =
    𝚺(abs2(X[i, i]) for i=1:minimum(size(X)))

sumOfSqrDiag(X::Union{𝔻{T}, ℍ{T}, 𝕃{T}}) where T<:RealOrComplex =
    𝚺(abs2(X[i, i]) for i=1:size(X, 1))

ssd=sumOfSqrDiag

"""
    sumOfSqrTril(X::AnyMatrix, k::Int=0)

**alias**: `sst`

 Given a real or complex `Matrix`, `Diagonal`, `Hermitian` or
 `LowerTriangular` matrix ``X`` (see [AnyMatrix type](@ref)),
 return the sum of squares of the elements
 in its lower triangle up to the ``k^{th}`` underdiagonal.

 `Matrix` ``X`` may be rectangular.

 ``k`` must be in range
 - `1-size(X, 1):c-1` for ``X`` `Matrix`, `Diagonal` or `Hermitian`,
 - `1-size(X, 1):0` for ``X`` `LowerTriangular`.

 For ``X`` `Diagonal` the result is
 - ``0`` if ``k<0``,
 - the sum of the squares of the diagonal elements otherwise.

 See julia [tril(M, k::Integer)](https://bit.ly/2Tbx8o7) function
 for numbering of diagonals.

 **See also**: [`sumOfSqr`](@ref), [`sumOfSqrDiag`](@ref).

 ## Examples
    using PosDefManifold
    A=[4. 3.; 2. 5.; 1. 2.]
    #3×2 Array{Float64,2}:
    # 4.0  3.0
    # 2.0  5.0
    # 1.0  2.0

    s=sumOfSqrTril(A, -1)
    # 9.0 = 1²+2²+2²

    s=sumOfSqrTril(A, 0)
    # 50.0 = 1²+2²+2²+4²+5²

"""
function sumOfSqrTril(X::AnyMatrix, k::Int=0)
    (r, c) = size(X)
    X isa 𝕃 ? range = (1-r:0) : range = (1-r:c-1)
    if k in range
        s=eltype(X)(0)
        for j=1:c, i=max(j-k, 1):r @inbounds s+=abs2(X[i, j]) end
        return s
    else
        @error "in LinearAmgebraInP.sumOfSqrTRil function: argument k is out of bounds"
    end
end
sst=sumOfSqrTril

"""
    (1) tr(P::ℍ{T}, Q::ℍ{T})
    (2) tr(P::ℍ{T}, M::𝕄{T})
    (3) tr(D::𝔻{T}, H::Union{ℍ{T}, 𝕄{T}})
    (4) tr(H::Union{ℍ{T}, 𝕄{T}}, D::𝔻{T})
          for all above: where T<:RealOrComplex

 Given (1) two `Hermitian` positive definite matrix ``P`` and ``Q``,
 return the trace of the product ``PQ``.
 This is real even if ``P`` and ``Q`` are complex.

 ``P`` must always be flagged as `Hermitian`. See [typecasting matrices](@ref).

 In (2) ``Q`` is a `Matrix` object,
 in which case return
 - a real trace if the product ``PQ`` is real or if it has all positive real eigenvalues.
 - a complex trace if the product ``PQ`` is not real and has complex eigenvalues.

 Methods (3) and (4) return the trace of the product ``DH`` or ``HD``,
 where ``D`` is a `Diagonal` matrix and ``H`` an ``Hermitian``
 or ``Matrix`` object. The result is of the same type as the input matrices.

 For all methods all arguments must be of the same type.

 ## Math
 Let ``P`` and ``Q`` be `Hermitian` matrices, using the properties of the trace
 (e.g., the cyclic property and the similarity invariance) you can use this
 function to fast compute the trace of several expressions. For example:

 ``\\textrm{tr}(PQ)=\\textrm{tr}(P^{1/2}QP^{1/2})``

 and

 ``\\textrm{tr}(PQP)=\\textrm{tr}(P^{2}Q)`` (see example below).

 **See**: [trace](https://bit.ly/2HoOLiM).

 **See also**: [`DiagOfProd`](@ref), [`tr1`](@ref).

 ## Examples
    using PosDefManifold
    P=randP(ComplexF64, 5) # generate a random complex positive definite matrix 5x5
    Q=randP(ComplexF64, 5) # generate a random complex positive definite matrix 5x5
    tr(P, Q) ≈ tr(P*Q) ? println(" ⭐ ") : println(" ⛔ ")
    tr(P, Q) ≈ tr(sqrt(P)*Q*sqrt(P)) ? println(" ⭐ ") : println(" ⛔ ")
    tr(sqr(P), Q) ≈ tr(P*Q*P) ? println(" ⭐ ") : println(" ⛔ ")

"""
tr(P::ℍ{T}, Q::ℍ{T}) where T<:RealOrComplex = real(tr(DiagOfProd(P, Q)))

function tr(P::ℍ{T}, M::𝕄{T}) where T<:RealOrComplex
    λ = [colProd(P, M, i, i) for i=1:size(P, 2)]
    OK=true
    for l in λ
        if imag(l) ≉  0
            OK=false
            break
        end
    end
    if OK return real(𝚺(λ)) else return 𝚺(λ) end
end

function tr(D::𝔻{T}, H::Union{ℍ{T}, 𝕄{T}}) where T<:RealOrComplex
    s=T(0)
    for i=1:size(D, 1) @inbounds s += D[i, i] * H[i, i] end
    return s
end

tr(H::Union{ℍ{T}, 𝕄{T}}, D::𝔻{T}) where T<:RealOrComplex = tr(D, H)

"""
    (1) quadraticForm(v::Vector{T}, P::ℍ{T}) where T<:Real
    (2) quadraticForm(v::Vector{T}, L::𝕃{T}) where T<:Real
    (3) quadraticForm(v::Vector{T}, X::𝕄{T}, forceLower::Bool=false) where T<:Real
    (4) quadraticForm(v::Vector{S}, X::Union{𝕄{S}, ℍ{S}}) where S<:Complex

 **alias**: `qf`

 (1) Given a real vector ``v`` and a real `Hermitian` matrix ``P``,
 compute the quadratic form

 ``v^TPv``,

 where the superscript *T* denotes transpose,
 using only the lower triangular part of ``P``.

 (2) As in (1), given a real vector ``v``
 and the `LowerTriangular` view ``L`` of a real symmetric matrix.

 (3) As in (1), given a real vector ``v``
 and a real generic `Matrix` ``M``, if `forceLower=true`. If `forceLower=false`,
 the product ``v^TMv`` is evaluated using the whole matrix ``M``.

 (4) Quadratic form ``v^HPv``, where superscript *H* denotes complex conjugate
 and transpose, for a complex vector `v` and  `Matrix` or `Hermitian` matrix.
 The whole matrix is then used.

 ## Math

 For ``v`` and ``X`` real and ``X`` symmetric, the quadratic form is

 ``\\sum_i(v_i^2x_{ii})+\\sum_{i>j}(2v_iv_jx_{ij})``.

 This is used in (1), (2) and (3).

 ## Examples
    using PosDefManifold
    P=randP(5) # generate a random real positive definite matrix 5x5
    v=randn(5)
    q1=quadraticForm(v, P) # or q1=quad(v, P)
    # obtain a lower Triangular view of P
    L=LowerTriangular(Matrix(P)) # or L=𝕃(Matrix(P))
    q2=quadraticForm(v, L)
    q1 ≈ q2 ? println(" ⭐ ") : println(" ⛔ ")

"""
quadraticForm(v::Vector{T}, P::ℍ{T}) where T<:Real = qf(v, 𝕃(P))

quadraticForm(v::Vector{T}, X::𝕄{T}, forceLower::Bool=false) where T<: Real = forceLower==true ? qf(v, 𝕃(X)) : v'*X*v

quadraticForm(v::Vector{T}, X::Union{𝕄{T}, ℍ{T}}) where T<:Complex = v'*X*v

function quadraticForm(v::Vector{T}, L::𝕃{T}) where T<:Real
    r=length(v)
    s=T(0)
    for j=1:r-1
        @inbounds s+=(v[j]^2 * L[j, j])
        for i=j+1:r @inbounds s+=2*v[i]*v[j]*L[i, j]
        end
    end
    @inbounds s+=(v[r]^2 * L[r, r])
    return s
end
qf=quadraticForm

"""
    fidelity(P::ℍ{T}, Q::ℍ{T}) where T<:RealOrComplex

 Given two positive definte `Hermitian` matrices ``P`` and ``Q``,
 return their *fidelity*:

  ``tr\\big(P^{1/2}QP^{1/2}\\big)^{1/2}.``

  This is used in quantum physics and is related to the
 [Wasserstein](@ref) metric. See for example Bhatia, Jain and Lim (2019b)[🎓](@ref).

 ## Examples
    using PosDefManifold
    P=randP(5);
    Q=randP(5);
    f=fidelity(P, Q)

"""
function fidelity(P::ℍ{T}, Q::ℍ{T}) where T<:RealOrComplex
    A = √(P)
    return tr(√ℍ(A*Q*A'))
end

#  ---------------------------------
## 4. Diagonal functions of matrices
#  ---------------------------------
"""
    fDiag(func::Function, X::AnyMatrix, k::Int=0)

 **alias**: `𝑓𝔻`

 Applies function `func` element-wise to the elements of the ``k^{th}``
 diagonal of real or complex `Diagonal`, `LowerTriangular`, `Matrix`
 or `Hermitian` matrix ``X`` and return a diagonal matrix with these elements.
 ``X`` must be square in all cases, but for the 𝕄=`Matrix` type argument,
 in which case it may be of dimension *r⋅c*, with *r ≠ c*.

 See julia [tril(M, k::Integer)](https://bit.ly/2Tbx8o7) function
 for numbering of diagonals.

 Bt default the main diagonal is considered.
 - If ``X`` is `Diagonal`, ``k`` must be zero (main diagonal).
 - If ``X`` is `LowerTriangular`, ``k`` cannot be positive.

 Note that if ``X`` is rectangular the dimension of the result depends
 on the size of ``X`` and on the chosen diagonal.
 For example,
 - *r ≠ c* and ``k``=0 (main diagonal), the result will be of dimension min*(r,c)*⋅*min(r,c)*,
 - ``X`` *3⋅4* and ``k=-1``, the result will be *2⋅2*,
 - ``X`` *3⋅4* and ``k=1``, the result will be *3⋅3*, etc.

!!! note "Nota Bene"
    The function `func` must support the `func.` syntax and therefore
    must be able to apply element-wise to the elements of the chosen diagonal
    (this includes anonymous functions). If the input matrix is complex, the function `func`
    must be able to support complex arguments.

 **See also**: [`DiagOfProd`](@ref), [`tr`](@ref).

 ## Examples
    using PosDefManifold
    P=randP(5) # use P=randP(ComplexF64, 5) for generating an Hermitian matrix
    D=fDiag(inv, P, -1)   # diagonal matrix with the inverse of the first sub-diagonal of P
    (Λ, U) = evd(P)       # Λ holds the eigenvalues of P, see evd
    Δ=fDiag(log, Λ)       # diagonal matrix with the log of the eigenvalues
    Δ=fDiag(x->x^2, Λ)    # using an anonymous function for the square of the eigenvalues
"""
fDiag(func::Function, X::𝔻{T}, k::Int=0) where T<:RealOrComplex = func.(X)

function fDiag(func::Function, X::𝕃{T}, k::Int=0)  where T<:RealOrComplex
 if k>0 @error("in function fDiag (linearAlgebra.jl): k argument cannot be positive.")
 else return 𝔻(func.(diag(X, k)))
 end
end

fDiag(func::Function, X::Union{𝕄{T}, ℍ{T}}, k::Int=0)  where T<:RealOrComplex =
      𝔻(func.(diag(X, k)))

𝑓𝔻=fDiag

"""
    DiagOfProd(P::ℍ{T}, Q::ℍ{T}) where T<:RealOrComplex

 **alias**: `dop`

 Return the `Diagonal` matrix holding the diagonal of the product ``PQ``
 of two `Hermitian` matrices `P` and `Q`. Only the diagoanl part
 of the product is computed.

 **See also**: [`tr`](@ref), [`fDiag`](@ref).

 ## Examples
    using PosDefManifold, LinearAlgebra
    P, Q=randP(5), randP(5)
    DiagOfProd(P, Q)≈Diagonal(P*Q) ? println("⭐ ") : println("⛔ ")
"""
DiagOfProd(P::ℍ{T}, Q::ℍ{T}) where T<:RealOrComplex =
           𝔻([colProd(P, Q, i, i) for i=1:size(P, 1)])

dop=DiagOfProd

#  -------------------------------
## 5. Unitary functions of matrices
#  -------------------------------
"""
    mgs(X::𝕄{T}, numCol::Int=0) where T<:RealOrComplex

 Modified (stabilized) [Gram-Schmidt orthogonalization](https://bit.ly/2YE6zvy)
 of the columns of square or tall matrix ``X``, which can be comprised of real
 or complex elements.
 The orthogonalized ``X`` is returned by the function. ``X`` is not changed.

 All columns are orthogonalized by default. If instead argument `numCol` is provided,
 then only the first `numCol` columns of ``X`` are orthogonalized.
 In this case only the firt `numCol` columns will be returned.

 ## Examples
    using LinearAlgebra, PosDefManifold
    X=randn(10, 10);
    U=mgs(X)        # result is 10⋅10
    U=mgs(X, 3)     # result is 10⋅3
    U'*U ≈ I ? println(" ⭐ ") : println(" ⛔ ")
    # julia undertands also:
    U'U ≈ I ? println(" ⭐ ") : println(" ⛔ ")

"""
function mgs(X::𝕄{T}, numCol::Int=0) where T<:RealOrComplex
    r=size(X, 1)
    if numCol != 0 && numCol in 2:size(X, 2) c=numCol else c=size(X, 2) end
    U=Matrix{T}(undef, r, c)
    U[:, 1] = X[:, 1]/colNorm(X, 1);
    for i in 2:c
        U[:, i]=X[:, i]
        for j in 1:i-1
            s = colProd(U, j, i) / sumOfSqr(U, j)
            U[:, i] -= s * U[:, j]
        end
        normalizeCol!(U, i)
    end
    return U
end # mgs function

#  ------------------------------
## 6. Matrix function of matrices
#  ------------------------------

#  -----------------------------------------------
## 7. Spectral decompositions of positive matrices
#  -----------------------------------------------

"""
    evd(S::ℍ{T}) where T<:RealOrComplex

 Given a positive semi-definite matrix ``S``,
 returns a 2-tuple ``(Λ, U)``, where ``U`` is the matrix holding in columns
 the eigenvectors and ``Λ`` is the matrix holding the eigenvalues on the diagonal.
 This is the output of Julia
 [eigen](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.eigen)
 function in ``UΛU'=S`` form.

 As for the `eigen` function, the eigenvalues and associated
 eigenvectors are sorted by increasing values of eigenvalues.

 ``S`` may be real or complex and must be flagged by Julia as `Hermitian`.
 See [typecasting matrices](@ref).

 **See also**: [`spectralFunctions`](@ref).

 ## Examples
    using PosDefManifold
    A=randn(3, 3);
    S=ℍ(A+A');
    Λ, U=evd(S); # which is equivalent to (Λ, U)=evd(P)
    (U*Λ*U') ≈ S ? println(" ⭐ ") : println(" ⛔ ")
    # => UΛU'=S, UΛ=SU, ΛU'=U'S
"""
function evd(S::ℍ{T}) where T<:RealOrComplex # return tuple (Λ, U)
    F = eigen(S)
    return  𝔻(F.values), F.vectors # 𝔻=LinearAlgebra.Diagonal
end

"""
    (1) spectralFunctions(P::ℍ{T}, func) where T<:RealOrComplex
    (2) spectralFunctions(D::𝔻{S}, func) where S<:Real

 (1) This is the *mother function* for all spectral functions of eigenvalues implemented
 in this library, which are:
 - `pow`     (power),
 - `isqrt`   (inverse square root).

 The function `sqr` (square) does not use it, as it can be obtained more
 efficiently by simple multiplication.

 You can use this function if you need another spectral function of eigenvalues
 besides those and those already implemented in the standard package `LinearAlgebra`.
 In general, you won't call it directly.

 `func` is the function that will be applied on the eigenvalues.

 ``P`` must be flagged as Hermitian. See [typecasting matrices](@ref).
 It must be a positive definite or positive semi-definite matrix,
 depending on `func`.

 A special method is provided for real `Diagonal` matrices (2).

!!! note "Nota Bene"
    The function `func` must support the `func.` syntax and therefore
    must be able to apply element-wise to the eigenvalues
    (those include anonymous functions).

 ### Maths

 The definition of spectral functions for a positive definite matrix ``P``
 is at it follows:

 ``f\\big(P\\big)=Uf\\big(Λ\\big)U^H,``

 where ``U`` is the matrix holding in columns the eigenvectors of ``P``,
 ``Λ`` is the matrix holding on diagonal its eigenvalues and ``f`` is
 a function applying element-wise to the eigenvalues.

 **See also**: [`evd`](@ref).

 ## Examples
    using LinearAlgebra, PosDefManifold
    n=5
    P=randP(n) # P=randP(ComplexF64, 5) to generate an Hermitian complex matrix
    noise=0.1;
    Q=spectralFunctions(P, x->x+noise) # add white noise to the eigenvalues
    tr(Q)-tr(P) ≈ noise*n ? println(" ⭐ ") : println(" ⛔ ")

"""
function spectralFunctions(P::ℍ{T}, func::Function) where T<:RealOrComplex
    F = eigen(P)
    ispos(F.values, msg="function "*string(func)*": at least one eigenvalue is smaller than the default tolerance")
    # optimize by computing only the upper trinagular part
    return ℍ(F.vectors * 𝔻(func.(F.values)) * F.vectors')
end

spectralFunctions(D::𝔻{T}, func::Function) where T<:Real = func.(D)

"""
    (1) pow(P::ℍ{T}, args...) where T<:RealOrComplex
    (2) pow(D::𝔻{S}, args...) where S<:Real

 (1) Given a positive semi-definite `Hermitian` matrix ``P``, return the power
 ``P^{r_1}, P^{r_2},...``
 for any number of exponents ``r_1, r_2,...``.
 It returns a tuple comprising as many elements as arguments passed after ``P``.

 ``P`` must be flagged as `Hermitian`. See [typecasting matrices](@ref).

 ``arg1, arg2,...`` are real numbers.

 A special method is provided for real `Diagonal` matrices (2).

 **See also**: [`invsqrt`](@ref).

 ## Examples
    using LinearAlgebra, PosDefManifold
    P=randP(5);     # use P=randP(ComplexF64, 5) for generating an Hermitian matrix
    Q=pow(P, 0.5);            # =>  QQ=P
    Q, W=pow(P, 0.5, -0.5);
    W*P*W ≈ I ? println(" ⭐ ") : println(" ⛔ ")
    Q*Q ≈ P ? println(" ⭐ ") : println(" ⛔ ")
    R, S=pow(P, 0.3, 0.7);
    R*S ≈ P ? println(" ⭐ ") : println(" ⛔ ")

"""
pow(P::ℍ{T}, p) where T<:RealOrComplex = spectralFunctions(P, x->x^p) # one argument

pow(D::𝔻{T}, p)  where T<:Real = spectralFunctions(D, x->x^p) # one argument

function pow(P::ℍ{T}, args...) where T<:RealOrComplex # several arguments
    (Λ, U) = evd(P)
    ispos(Λ, msg="function pow: at least one eigenvalue is smaller than the default tolerance")
    # optimize by computing only the upper trinagular part
    return  (ℍ(U * Λ^p * U') for p in args)
end

function pow(D::𝔻{T}, args...) where T<:Real  # several arguments
    ispos(D, msg="function pow: at least one eigenvalue is smaller than the default tolerance")
    return  (D^p for p in args)
end

"""
    (1) invsqrt(P{T}::ℍ) where T<:RealOrComplex
    (2) invsqrt(D{S}::𝔻) where S<:Real

 Given a positive definite `Hermitian` matrix ``P``,
 compute the inverse of the principal
 square root ``P^{-1/2}``.

 ``P`` must be flagged as Hermitian. See [typecasting matrices](@ref).

 A special method is provided for real `Diagonal` matrices (2).

 ### Maths

 The principal square root of a positive definite matrix ``P`` is the only
 symmetric (if ``P`` is real) or Hermitian (if ``P`` is complex) square root.
 Its inverse ``P^{-1/2}`` is also named the **whitening** or **sphering**
 matrix since``P^{-1/2}PP^{-1/2}=I``.

 **See**: [typecasting matrices](@ref).

 **See also**: [`pow`](@ref).

 ## Examples
    using LinearAlgebra, PosDefManifold
    P=randP(ComplexF64, 5);
    Q=invsqrt(P);
    Q*P*Q ≈ I ? println(" ⭐ ") : println(" ⛔ ")

"""
invsqrt(P::ℍ{T}) where T<:RealOrComplex = spectralFunctions(P, x->1/sqrt(x))

invsqrt(D::𝔻{T}) where T<:Real = spectralFunctions(D, x->1/sqrt(x))

"""
    (1) sqr(P::ℍ{T}) where T<:RealOrComplex
    (2) sqr(X::Union{𝕄{T}, 𝕃{T}, 𝔻{S}}) where T<:RealOrComplex where S<:Real

 (1) Given a positive semi-definite `Hermitian` matrix ``P``,
 compute its square ``P^{2}``.

 ``P`` must be flagged as Hermitian. See [typecasting matrices](@ref).

 A method is provided also for generic matrices of the `Matrix` type,
 `LowerTriangular` matrices and real `Diagonal` matrices (2). The output
 is of the same type as the input.

 **See also**: [`pow`](@ref).

 ## Examples
    using PosDefManifold
    P=randP(5);
    P²=sqr(P);  # =>  P²=PP
    sqrt(P²)≈ P ? println(" ⭐ ") : println(" ⛔ ")

"""
sqr(P::ℍ{T}) where T<:RealOrComplex = ℍ(P*P)

sqr(X::Union{𝕄{T}, 𝕃{T}, 𝔻{S}}) where T<:RealOrComplex where S<:Real = X*X

"""
    powerIterations(H::Union{ℍ{T}, 𝕄{T}}, q::Int;
    <evalues=false, tol::Real=0, maxiter=300, ⍰=false>)  where T<:RealOrComplex

    powerIterations(L::𝕃{S}, q::Int;
    <evalues=false, tol::Real=0, maxiter=300, ⍰=false)> where S<:Real

 **alias**: `powIter`

 (1) Compute the ``q`` eigenvectors associated to the ``q`` largest (real) eigenvalues
 of real or complex `Hermitian` or `Matrix` ``H`` using the
 [power iterations](https://bit.ly/2JSo0pb) +
 [Gram-Schmidt orthogonalization](https://bit.ly/2YE6zvy) as suggested by Strang.
 The eigenvectors are returned with the same type as the elements of ``H``.

 ``H`` must have real eigenvalues, that is, it must be a symmetric matrix if it is real
 or an Hermitian matrix if it is complex.

 (2) as in (1), but using only the `LowerTriangular` view ``L`` of a matrix.
 This option is available only for real matrices (see below).

 The following are *<optional keyword arguments>*:
 - ``tol`` is the tolerance for the convergence of the power method (see below),
 - ``maxiter`` is the maximum number of iterations allowed for the power method,
 - if ``⍰=true``, the convergence of all iterations will be printed,
 - if ``evalues=true``, return the 4-tuple ``(Λ, U, iterations, covergence)``,
 - if ``evalues=false`` return the 3-tuple ``(U, iterations, covergence)``.

!!! note "Nota Bene"
    Differently from the [`evd`](@ref) function, the eigenvectors and
    eigenvalues are sorted by decreasing order of eigenvalues.

    If ``H`` is `Hermitian` and real, only its lower triangular part is used
    for computing the power iterations, like in (2). In this case the
    [BLAS.symm](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.BLAS.symm)
    routine is used.
    Otherwise the [BLAS.gemm](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.BLAS.gemm)
    routine is used. See [Threads](@ref).

    ``tol`` defaults to 100 times the square root of `Base.eps` of the type
    of ``H``. This corresponds to requiring the relative convergence criterion
    over two successive iterations to vanish for about half the significant
    digits minus 2.

**See also**: [`mgs`](@ref).

 ## Examples
    using LinearAlgebra, PosDefManifold
    # Generate an Hermitian (complex) matrix
    H=randP(ComplexF64, 10);
    # 3 eigenvectors and eigenvalues
    U, iterations, convergence=powIter(H, 3, ⍰=true)
    # all eigenvectors
    Λ, U, iterations, convergence=powIter(H, size(H, 2), evalues=true, ⍰=true);
    U'*U ≈ I && U*Λ*U'≈H ? println(" ⭐ ") : println(" ⛔ ")

    # passing a `Matrix` object
    Λ, U, iterations, convergence=powIter(Matrix(H), 3, evalues=true)

    # passing a `LowerTriangular` object (must be a real matrix in this case)
    L=𝕃(randP(10))
    Λ, U, iterations, convergence=powIter(L, 3, evalues=true)

"""
function powerIterations(H::𝕄{T}, q::Int;
  evalues=false, tol::Real=0, maxiter=300, ⍰=false) where T<:RealOrComplex

    (n, sqrtn, type) = size(H, 1), √(size(H, 1)), eltype(H)
    tol==0 ? tolerance = √eps(real(type))*1e2 : tolerance = tol
    msg1="Power Iterations reached a saddle point at:"
    msg2="Power Iterations reached the max number of iterations at:"
    U=randn(type, n, q) # initialization
    normalizeCol!(U, 1:q)
    💡=similar(U) # 💡 is the poweriteration matrix
    (iter, conv, oldconv) = 1, 0., maxpos
    ⍰ && @info("Running Power Iterations...")
    while true
        # power iteration U-<H*U of the q vectors of U and their Gram-Schmidt Orth.
        type<:Real ? 💡=mgs(BLAS.symm('L', 'L', H, U)) : 💡=mgs(BLAS.gemm('N', 'N', H, U))
        conv = √norm((💡)' * U - I) / sqrtn # relative difference to identity
        (saddlePoint = conv ≈ oldconv)  && @info(msg1, iter)
        (overRun     = iter == maxiter) && @warn(msg2, iter)
        #diverged    = conv > oldconv && @warn(msg3, iter)
        ⍰ && println("iteration: ", iter, "; convergence: ", conv)
        if conv<=tolerance || saddlePoint==true || overRun==true
            break
        else U = 💡 end
        oldconv=conv
        iter += 1
    end # while
    if evalues == false
        return (💡, iter, conv)
    else
        type<:Real ? d=[qf(💡[:, i], H, true) for i=1:q] : d=[qf(💡[:, i], H) for i =1:q]
        return (𝔻(real(d)), 💡, iter, conv)
    end
end

powerIterations(H::ℍ{T}, q::Int;
    evalues=false, tol::Real=0, maxiter=300, ⍰=false) where T<:RealOrComplex =
    powerIterations(Matrix(H), q; evalues=evalues, tol=tol, maxiter=maxiter, ⍰=⍰)

powerIterations(L::𝕃{T}, q::Int;
        evalues=false, tol::Real=0, maxiter=300, ⍰=false) where T<:Real =
    powerIterations(𝕄(L), q; evalues=evalues, tol=tol, maxiter=maxiter, ⍰=⍰)

powIter=powerIterations

#  -----------------------------------------------
## 8. Decompositions involving triangular matrices
#  -----------------------------------------------
"""
    (1) choL(P::ℍ{T}) where T<:RealOrComplex
    (2) choL(D::𝔻{S}) where S<:Real

 (1) Given a real or complex positive definite `Hermitian` matrix ``P``,
 return the *Cholesky lower triangular factor* ``L``
 such that ``LL^H=P``. To obtain ``L^H`` or both ``L`` and ``L^H``, use instead
 julia function [cholesky](https://bit.ly/2u9Hw5P).

 On output, ``L`` is of type [`LowerTriangular`](https://bit.ly/2U511f3).

 (2) For a real `Diagonal` matrix ``D``, return ``D^{1/2}``.

 ## Examples
    using PosDefManifold
    P=randP(5);
    L=choL(P);
    L*L'≈ P ? println(" ⭐ ") : println(" ⛔ ")

"""
function choL(P::ℍ{T}) where T<:RealOrComplex
    choP = cholesky(P)
    return choP.L
end

choL(D::𝔻{T}) where T<:Real = √D


text:
using SimplexQuad

using Test

function integrate(f, X::AbstractMatrix, W::AbstractVector)
    np = size(X, 1)
    @assert size(W, 1) == np
    return sum(W[i] * f(X[i, :]) for i in 1:np)
end

const Types = [Float32, Float64]
@testset "Integrate polynomials T=$T D=$D N=$N" for T in Types,
D in 1:5,
N in 1:(11 - D)

    X, W = simplexquad(T, N, D)

    # Check returned types
    X::AbstractMatrix{T}
    W::AbstractVector{T}
    # Check returned sizes
    np = N^D
    @test size(X, 1) == np
    @test size(X, 2) == D
    @test size(W, 1) == size(X, 1)
    # Ensure all points lie strictly inside the simplex
    function inside(x)
        all(>(0), x) || return false
        all(<(1), x) || return false
        λ = 1 - sum(x)
        return 0 < λ < 1
    end
    @test all(inside(X[i, :]) for i in 1:np)
    # Ensure weights are strictly positive
    @test all(>(0), W)

    # Expect polynomials of orders P<= to be integrated exactly
    P = N
    pmin = CartesianIndex(ntuple(d -> 0, D))
    pmax = CartesianIndex(ntuple(d -> P, D))
    for p in pmin:pmax
        sum(p.I) <= P || continue
        f(x) = prod(x[i]^p[i] for i in 1:D)
        computed = integrate(f, X, W)
        computed::T
        # Mathematica suggests this neat result
        expected = prod(T.(factorial.(p.I))) / factorial(D + sum(p.I))
        expected::T
        @test abs(computed - expected) <= 10 * eps(T)
    end
end


text:
# This is part of Quiz.jl. License is MIT.

abstract type QUIZ end                  # Module's top-level
abstract type   qEngine <: QUIZ end     # Quiz engine types
abstract type   qSing   <: QUIZ end     # Single problem types
abstract type   qElem   <: QUIZ end     # Elements for single problems
abstract type     qBare <: qElem end    # Bare items to compose the problem
abstract type     qColl <: qElem end    # Collection types of qBare eltypes

@doc QUIZ """
`abstract type QUIZ end`\n
`Quiz.jl`'s top-level type.
"""

@doc qEngine """
`abstract type qEngine <: QUIZ end`\n
`Quiz.jl`'s engine abstract type.
"""

@doc qSing """
`abstract type qSing <: QUIZ end`\n
`Quiz.jl`'s single problem abstract type.
"""

@doc qElem """
`abstract type qElem <: QUIZ end`\n
`Quiz.jl`'s elements for single problems.
"""

export QUIZ
export qEngine, qSing, qElem
export qBare, qColl

"""
qSing
    qProb,

qProb: Concrete container of qElem's that fully specify a single problem.

qElem
    qBare
        qT, qC, qV, qO, qF, qE
    qColl
        q𝕋

qColl: Concrete collections of qBare eltypes:
    - q𝕋: table

qSing: Concrete types of single items to compose a single problem:
    - qT: text
    - qC: constant
    - qV: variable (input parameter range allowing for permutations)
    - qO: output (values pertaining to questions)
    - qF: figure
    - qE: equation
"""


text:
module ArbTools

using Nemo
using Printf
using RecipesBase
using Colors

ArbReal = arb

import Base: length, show, inv, sqrt, log, log1p, atan, asin, acos,
    exp, sin, cos, tan, sinpi, cospi, sinh, cosh, sinc, +, -, *, /, ^,
    div, sincos, isnan

import Nemo: rsqrt, tanpi, gamma, lgamma, rgamma, digamma, compose,
    derivative, integral, zeta

export setinterval, getinterval

export arb_series

export enclosemaximum
export maximumpoly

export isolateroots

include("arb.jl")
include("acb.jl")
include("arb_series.jl")

include("trace.jl")
include("enclosemaximum.jl")
include("maximumtaylor.jl")
include("maximumpoly.jl")

include("isolateroots.jl")

include("integrate.jl")

end # module


text:
# https://juliagraphics.github.io/Gtk.jl/latest/manual/layout/#Box-1

using Gtk

# win = GtkWindow("New title")
# hbox = GtkBox(:h)  # :h makes a horizontal layout, :v a vertical layout
# push!(win, hbox)
# cancel = GtkButton("Cancel")
# ok = GtkButton("OK")
# push!(hbox, cancel)
# push!(hbox, ok)
# showall(win)

win = GtkWindow("New Title")
hbox = GtkBox(:h)    # :h makes a horizontal layout, :v a vertical layout
push!(win, hbox)
cancel = GtkButton("Cancel")
ok = GtkButton("OK")
push!(hbox, cancel)
push!(hbox, ok)

showall(win)

# julia> length(hbox)
# 2
#
# julia> get_gtk_property(hbox[1], :label, String)
# "Cancel"
#
# julia> get_gtk_property(hbox[2], :label, String)
# "OK"

label = get_gtk_property(hbox[1], :label, String)
label = get_gtk_property(hbox[1], "label", String)
label = get_gtk_property(hbox[2], :label, String)

# set_gtk_property!(hbox,:expand,ok,true)
# set_gtk_property!(hbox,:spacing,10)

set_gtk_property!(hbox, :expand, ok, true)
set_gtk_property!(hbox, :spaceing, 10)

# destroy(hbox)
# ok = GtkButton("OK")
# cancel = GtkButton("Cancel")
# hbox = GtkButtonBox(:h)
# push!(win, hbox)
# push!(hbox, cancel)
# push!(hbox, ok)

destroy(hbox)
ok = GtkButton("OK")
cancel = GtkButton("Cancel")
hbox = GtkButtonBox(:h)
push!(win, hbox)
push!(hbox, cancel)
push!(hbox, ok)
showall(win)


text:
# Copyright (c) 2013: Joey Huchette and contributors
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

mutable struct CallbackContext
    model::Optimizer
    ptr::Ptr{Cvoid}
    context_id::Clong
end
Base.cconvert(::Type{Ptr{Cvoid}}, x::CallbackContext) = x
Base.unsafe_convert(::Type{Ptr{Cvoid}}, x::CallbackContext) = x.ptr::Ptr{Cvoid}
Base.broadcastable(x::CallbackContext) = Ref(x)

mutable struct _CallbackUserData
    model::Optimizer
    callback::Function
end
Base.cconvert(::Type{Ptr{Cvoid}}, x::_CallbackUserData) = x
function Base.unsafe_convert(::Type{Ptr{Cvoid}}, x::_CallbackUserData)
    return pointer_from_objref(x)::Ptr{Cvoid}
end

function _cplex_callback_wrapper(
    p_context::Ptr{Cvoid},
    context_id::Clong,
    p_user_data::Ptr{Cvoid},
)
    user_data = unsafe_pointer_to_objref(p_user_data)::_CallbackUserData
    try
        user_data.callback(
            CallbackContext(user_data.model, p_context, context_id),
            context_id,
        )
    catch ex
        CPXcallbackabort(p_context)
        if !(ex isa InterruptException)
            rethrow(ex)
        end
    end
    return Cint(0)
end

"""
    column(cb_data::CallbackContext, x::MOI.VariableIndex)

Return the 1-indexed column associated with `x` in a callback.

The C API requires 0-indexed columns.
"""
function column(cb_data::CallbackContext, x::MOI.VariableIndex)
    return _info(cb_data.model, x).column
end

"""
    CallbackFunction(
        context_mask::UInt16 =
            CPX_CALLBACKCONTEXT_THREAD_UP |
            CPX_CALLBACKCONTEXT_THREAD_DOWN |
            CPX_CALLBACKCONTEXT_LOCAL_PROGRESS |
            CPX_CALLBACKCONTEXT_GLOBAL_PROGRESS |
            CPX_CALLBACKCONTEXT_CANDIDATE |
            CPX_CALLBACKCONTEXT_RELAXATION |
            CPX_CALLBACKCONTEXT_BRANCHING
        )
    )

Set a generic CPLEX callback function. Use `context_mask` to control where the
callback is called from.

Callback must be a function with signature:

    callback(cb_data::CallbackContext, context_id::Clong)

Before accessing `MOI.CallbackVariablePrimal`, you must call
`CPLEX.load_callback_variable_primal(cb_data, context_id)`.
"""
struct CallbackFunction <: MOI.AbstractCallback
    context_mask::UInt16

    function CallbackFunction(
        context_mask::UInt16 = CPX_CALLBACKCONTEXT_THREAD_UP |
                               CPX_CALLBACKCONTEXT_THREAD_DOWN |
                               CPX_CALLBACKCONTEXT_LOCAL_PROGRESS |
                               CPX_CALLBACKCONTEXT_GLOBAL_PROGRESS |
                               CPX_CALLBACKCONTEXT_CANDIDATE |
                               CPX_CALLBACKCONTEXT_RELAXATION |
                               CPX_CALLBACKCONTEXT_BRANCHING,
    )
        return new(context_mask)
    end
end

function MOI.set(model::Optimizer, cb::CallbackFunction, f::Function)
    @static if Sys.iswindows() && Sys.WORD_SIZE == 32
        error(
            "Callbacks are not supported on Win32. Use 64-bit Julia with " *
            "64-bit CPLEX.",
        )
    end
    if MOI.get(model, MOI.NumberOfThreads()) != 1
        @warn(
            "When using callbacks, make sure to set `NumberOfThreads` to `1` " *
            "using `MOI.set(model, MOI.NumberOfThreads(), 1)`. Bad things can" *
            " happen if you don't! Only ignore this message if you know what " *
            "you are doing."
        )
    end
    cpx_callback = @cfunction(
        _cplex_callback_wrapper,
        Cint,
        (Ptr{Cvoid}, Clong, Ptr{Cvoid})
    )
    user_data = _CallbackUserData(
        model,
        (context, context_id) -> begin
            model.callback_state = _CB_GENERIC
            f(context::CallbackContext, context_id::Clong)
            model.callback_state = _CB_NONE
        end,
    )
    ret = CPXcallbacksetfunc(
        model.env,
        model.lp,
        cb.context_mask,
        cpx_callback,
        user_data,
    )
    _check_ret(model, ret)
    model.generic_callback = user_data
    model.has_generic_callback = true
    return
end
MOI.supports(::Optimizer, ::CallbackFunction) = true

"""
    load_callback_variable_primal(cb_data, context_id)

Load the solution during a callback so that it can be accessed using
`MOI.CallbackVariablePrimal`.
"""
function load_callback_variable_primal(
    cb_data::CallbackContext,
    context_id::Clong,
)
    model = cb_data.model::Optimizer
    N = length(model.variable_info)
    resize!(model.callback_variable_primal, N)
    if context_id == CPX_CALLBACKCONTEXT_CANDIDATE
        ret = CPXcallbackgetcandidatepoint(
            cb_data,
            model.callback_variable_primal,
            Cint(0),
            Cint(N - 1),
            C_NULL,
        )
    elseif context_id == CPX_CALLBACKCONTEXT_RELAXATION
        ret = CPXcallbackgetrelaxationpoint(
            cb_data,
            model.callback_variable_primal,
            Cint(0),
            Cint(N - 1),
            C_NULL,
        )
    else
        error(
            "`load_callback_variable_primal` can only be called at " *
            "CPX_CALLBACKCONTEXT_CANDIDATE or CPX_CALLBACKCONTEXT_RELAXATION.",
        )
    end
    return
end

# ==============================================================================
#    MOI callbacks
# ==============================================================================

function _default_moi_callback(model::Optimizer)
    if MOI.get(model, MOI.NumberOfThreads()) != 1
        # The current callback system isn't thread-safe. As a work-around, set
        # the number of threads to 1, regardless of what the user intended.
        # We only do this for the MOI callbacks, since we assume the user knows
        # what they are doing if they use a solver-dependent callback.
        MOI.set(model, MOI.NumberOfThreads(), 1)
    end
    return (cb_data, cb_context) -> begin
        if cb_context == CPX_CALLBACKCONTEXT_CANDIDATE
            ispoint_p = Ref{Cint}()
            ret = CPXcallbackcandidateispoint(cb_data, ispoint_p)
            _check_ret(cb_data.model, ret)
            if ispoint_p[] == 0
                return  # No candidate point available
            end
            load_callback_variable_primal(cb_data, cb_context)
            if model.lazy_callback !== nothing
                model.callback_state = _CB_LAZY
                model.lazy_callback(cb_data)
            end
        elseif cb_context == CPX_CALLBACKCONTEXT_RELAXATION
            load_callback_variable_primal(cb_data, cb_context)
            if model.user_cut_callback !== nothing
                model.callback_state = _CB_USER_CUT
                model.user_cut_callback(cb_data)
            end
            if model.heuristic_callback !== nothing
                model.callback_state = _CB_HEURISTIC
                model.heuristic_callback(cb_data)
            end
        end
        model.callback_state = _CB_NONE
    end
end

function MOI.get(
    model::Optimizer,
    ::MOI.CallbackVariablePrimal{CallbackContext},
    x::MOI.VariableIndex,
)
    return model.callback_variable_primal[_info(model, x).column]
end

function MOI.get(::Optimizer, attr::MOI.CallbackNodeStatus{CallbackContext})
    if attr.callback_data.context_id == CPX_CALLBACKCONTEXT_CANDIDATE
        return MOI.CALLBACK_NODE_STATUS_INTEGER
    elseif attr.callback_data.context_id == CPX_CALLBACKCONTEXT_RELAXATION
        return MOI.CALLBACK_NODE_STATUS_FRACTIONAL
    end
    return MOI.CALLBACK_NODE_STATUS_UNKNOWN
end

# ==============================================================================
#    MOI.LazyConstraint
# ==============================================================================

function MOI.set(model::Optimizer, ::MOI.LazyConstraintCallback, cb::Function)
    model.lazy_callback = cb
    return
end
MOI.supports(::Optimizer, ::MOI.LazyConstraintCallback) = true

function MOI.submit(
    model::Optimizer,
    cb::MOI.LazyConstraint{CallbackContext},
    f::MOI.ScalarAffineFunction{Float64},
    s::Union{
        MOI.LessThan{Float64},
        MOI.GreaterThan{Float64},
        MOI.EqualTo{Float64},
    },
)
    if model.callback_state == _CB_USER_CUT
        throw(MOI.InvalidCallbackUsage(MOI.UserCutCallback(), cb))
    elseif model.callback_state == _CB_HEURISTIC
        throw(MOI.InvalidCallbackUsage(MOI.HeuristicCallback(), cb))
    elseif !iszero(f.constant)
        throw(
            MOI.ScalarFunctionConstantNotZero{Float64,typeof(f),typeof(s)}(
                f.constant,
            ),
        )
    end
    indices, coefficients = _indices_and_coefficients(model, f)
    sense, rhs = _sense_and_rhs(s)
    ret = CPXcallbackrejectcandidate(
        cb.callback_data,
        Cint(1),
        Cint(length(coefficients)),
        Ref(rhs),
        Ref{Cchar}(sense),
        Ref{Cint}(0),
        indices,
        coefficients,
    )
    return
end
MOI.supports(::Optimizer, ::MOI.LazyConstraint{CallbackContext}) = true

# ==============================================================================
#    MOI.UserCutCallback
# ==============================================================================

function MOI.set(model::Optimizer, ::MOI.UserCutCallback, cb::Function)
    model.user_cut_callback = cb
    return
end
MOI.supports(::Optimizer, ::MOI.UserCutCallback) = true

function MOI.submit(
    model::Optimizer,
    cb::MOI.UserCut{CallbackContext},
    f::MOI.ScalarAffineFunction{Float64},
    s::Union{
        MOI.LessThan{Float64},
        MOI.GreaterThan{Float64},
        MOI.EqualTo{Float64},
    },
)
    if model.callback_state == _CB_LAZY
        throw(MOI.InvalidCallbackUsage(MOI.LazyConstraintCallback(), cb))
    elseif model.callback_state == _CB_HEURISTIC
        throw(MOI.InvalidCallbackUsage(MOI.HeuristicCallback(), cb))
    elseif !iszero(f.constant)
        throw(
            MOI.ScalarFunctionConstantNotZero{Float64,typeof(f),typeof(s)}(
                f.constant,
            ),
        )
    end
    rmatind, rmatval = _indices_and_coefficients(model, f)
    sense, rhs = _sense_and_rhs(s)
    ret = CPXcallbackaddusercuts(
        cb.callback_data,
        Cint(1),
        Cint(length(rmatval)),
        Ref(rhs),
        Ref(sense),
        Ref{Cint}(0),
        rmatind,
        rmatval,
        Ref{Cint}(CPX_USECUT_FILTER),
        Ref{Cint}(0),
    )
    _check_ret(model, ret)
    return
end
MOI.supports(::Optimizer, ::MOI.UserCut{CallbackContext}) = true

# ==============================================================================
#    MOI.HeuristicCallback
# ==============================================================================

function MOI.set(model::Optimizer, ::MOI.HeuristicCallback, cb::Function)
    model.heuristic_callback = cb
    return
end
MOI.supports(::Optimizer, ::MOI.HeuristicCallback) = true

function MOI.submit(
    model::Optimizer,
    cb::MOI.HeuristicSolution{CallbackContext},
    variables::Vector{MOI.VariableIndex},
    values::MOI.Vector{Float64},
)
    if model.callback_state == _CB_LAZY
        throw(MOI.InvalidCallbackUsage(MOI.LazyConstraintCallback(), cb))
    elseif model.callback_state == _CB_USER_CUT
        throw(MOI.InvalidCallbackUsage(MOI.UserCutCallback(), cb))
    end
    ret = CPXcallbackpostheursoln(
        cb.callback_data,
        Cint(length(variables)),
        Cint[_info(model, var).column - 1 for var in variables],
        values,
        NaN,
        CPXCALLBACKSOLUTION_SOLVE,
    )
    _check_ret(model, ret)
    return MOI.HEURISTIC_SOLUTION_UNKNOWN
end
MOI.supports(::Optimizer, ::MOI.HeuristicSolution{CallbackContext}) = true


text:
#using ODE
#using Gillespie
using Sundials
using PyCall
#using RCall

@pyimport numpy.random as py_rand
#@rimport CircStats

"""
deterministic simulations (ODE solver)
currently the function only uses the solver ode45 from the package ODE

input:
		- func: function to integrate
		- init_vals: an array of floats for the initial values of each species
		- time_points: a range of floats of the time points for the integration
		- params: an array of floats of the parameters for the integration

output:
		- an array of arrays, with the solutions for each species at each given time point
"""
function abc_ode(func::Function,init_vals::Array{Float64,1},time_points::FloatRange{Float64},params::Array{Float64,1})

	#solve using ODE package
	#convert to function with no params argument
	#g = (t,v) -> func(t,v,params)
	#(t,res) = ode45(g,init_vals,time_points,points=:specified,reltol = 1e-5, abstol = 1e-8)

	#solve using Sundails package
	g = (t,v,v_dot) -> func(t,v,v_dot,params)
	d = Sundials.cvode(g,init_vals,collect(time_points))

	#convert Array{Float64,2} to Array{Array{Float64,1},1}
	res = Array{Array{Float64,1},1}(size(d)[1])
	for i in 1:size(d)[1]
		index=i
		arr = Array{Float64,1}(size(d)[2])
		for j in 1:size(d)[2]
			arr[j] = d[index]
	    	index+=size(d)[1]
	   	end
		res[i] = arr
    end

	return res

end

"""
stochastic simulations (SDE solver)
"""
function abc_sde()
	###TODO
end

"""
stochastic simulations (Gillespie Algorithm)

input:
		- func: function to integrate
		- stoch_matrix : ???
		- init_vals: an array of floats for the initial values of each species
		- time_end: last time point
		- params: an array of floats of the parameters for the integration

output:
		- an array of arrays, with the solutions for each species at each given time point

"""
function abc_gillespie(func::Function,stoch_matrix::Array{Int64,2},init_vals::Array{Int64,1},time_end::Float64,params::Array{Float64,1})

	result = ssa(init_vals, func, stoch_matrix, params, time_end)
	data = ssa_data(result)

	#TODO: convert data from  DataFrames.DataFrame to array of arrays so that the it includes only the species data points and is the same type as for the ODE return value

	return data

end

"""
1D diffusion gradient reaction function to integrate

Intput:
	# x = cell position (as distance from wound)

	# t = absolute time after wounding (in seconds)

	# diffusion parameters:
	# D = diffusion coefficient
	# A = source strength
	# tau0 = time the wound is signalling (in seconds)
"""
function diffusion_gradient(tau::Float64,D::Float64,A::Float64,t::Float64,x::Float64)

	#x0 is the position of the wount - can be set differently if needed
	x0=0
	res = A*exp(-(x-x0)^2/(4*D*(t-tau)))/sqrt(4.0*D*(t-tau)*1pi)
	return res
end

"""
1D diffusion gradient reaction - solve the reaction's function
Intput:
	# x = cell position (as distance from wound)
	# r = cell radios

	# t = absolute time after wounding (in seconds)

	# diffusion parameters:
	# D = diffusion coefficient
	# A = source strength
	# tau0 = time the wound is signalling (in seconds)
"""
function diffusion_gradient_solver(tau0::Float64,D::Float64,A::Float64,t::Float64,x::Float64,r::Float64)

	x1 = x-r
	x2 = x+r
	integrand1 = (tau) -> diffusion_gradient(tau,D,A,t,x1)
	integrand2 = (tau) -> diffusion_gradient(tau,D,A,t,x2)
	conc1 = quadgk(integrand1,0,min(tau0,t))
	conc2 = quadgk(integrand2,0,min(tau0,t))

	return (conc1[1],conc2[1])

end

"""
simulate a random walk for one cell - basic model of (w,b,p)
"""
function simulate_random_walk(n_steps::Int64,start_point::Tuple{Float64,Float64},w::Float64,p::Float64,b::Float64,bias_angle::Float64)

	#save two result vectors: step lengths and step angles
	steps_length = Array{Float64,1}(n_steps)
	steps_angles = Array{Float64,1}(n_steps)
	steps_points = Array{Tuple{Float64,Float64},1}(n_steps+1)
	steps_points[1]=start_point

	#define distributions

	#step length
	dt = 0.05
	truncated_normal = Truncated(Normal(0,1), 0, Inf)

	steps_length[1] , steps_angles[1] , steps_points[2] = simulate_random_walk_one_step(truncated_normal,dt,w,p,Inf,b,bias_angle,start_point)

	for i in 2:n_steps

		steps_length[i] , steps_angles[i] , steps_points[i+1] = simulate_random_walk_one_step(truncated_normal,dt,w,p,steps_angles[i-1],b,bias_angle,steps_points[i-1])

	end

	return Random_Walk(steps_length,steps_angles,steps_points)

end

"""
simulate a random walk model where the bias is changed every time step
Input:
		- n_steps and b and bias_angle need to be of same length
"""
function simulate_random_walk_changing_bias(n_steps::Int64,start_point::Tuple{Float64,Float64},w::Float64,p::Float64,b::Array{Float64,1},bias_angle::Array{Float64,1})

	#save two result vectors: step lengths and step angles
	steps_length = Array{Float64,1}(n_steps)
	steps_angles = Array{Float64,1}(n_steps)
	steps_points = Array{Tuple{Float64,Float64},1}(n_steps+1)
	steps_points[1]=start_point

	#define distributions

	#step length
	dt = 0.001
	truncated_normal = Truncated(Normal(0,1), 0, Inf)

	steps_length[1] , steps_angles[1] , steps_points[2] = simulate_random_walk_one_step(truncated_normal,dt,w,p,Inf,b[1],bias_angle[1],start_point)

	for i in 2:n_steps
		steps_length[i] , steps_angles[i] , steps_points[i+1] = simulate_random_walk_one_step(truncated_normal,dt,w,p,steps_angles[i-1],b[i],bias_angle[i],steps_points[i-1])
	end

	return Random_Walk(steps_length,steps_angles,steps_points)

end

"""
simulate one step of a random walk
"""
function simulate_random_walk_one_step(truncated_normal::Distributions.Truncated{Distributions.Normal{Int64},Distributions.Continuous},
										dt::Float64,
										w::Float64,
										p::Float64,
										persistant_angle::Float64,
										b::Float64,
										bias_angle::Float64,
										prev_point::Tuple{Float64,Float64})

	#randomly choose the next step length, according to the step distribution
	next_length = rand(truncated_normal) * sqrt(dt)

	#randomly choose if the next step is going to be biased or persistent
	rand_w = rand()

	if rand_w<=w
		#biased walk
		if b==0.0
			#sample angle from wrapped uniform distribution
			next_angle = pi*rand(Uniform(-1,1))
		elseif b==1.0
			next_angle = bias_angle
		else
			#sample angle from wrapped normal distribution with mean=bias direction and sigma=-2log(b)
			bias_kappa = 1/(-2*log(b))

			#generate next random angle by calling python's Von Mises distribution
			next_angle = py_rand.vonmises(bias_angle,bias_kappa)

			# """
			# #generate next random angle using Julia's distributions package
			# bias_d = VonMises(bias_angle,bias_kappa)
			# next_angle = rand(sampler(bias_d))
			# """

			# """
			# #generate the next random angle by calling R's wrapped normal distribution
			# next_angle = R"rwrpnorm(1, $bias_angle, $b)"[1]
			# #convert the angle (R returns values between 0 and 2pi)
			# next_angle = (next_angle + pi)%(2*pi)-pi
			# """
		end
	else
		#persistant walk
		if p==0.0
			#sample angle from wrapped uniform distribution
			next_angle = pi*rand(Uniform(-1,1))
		elseif p==1.0
			next_angle = persistant_angle
		else
			#sample angle from wrapped normal distribution with mean=prev_angle direction and sigma=-2log(b)
			if persistant_angle==Inf
				#there is no previous angle to base the distribution on
				next_angle = pi*rand(Uniform(-1,1))
			else
				persistant_kappa = 1/(-2*log(p))

				#generate next random angle by calling python's Von Mises distribution
				next_angle = py_rand.vonmises(persistant_angle,persistant_kappa)

				# """
				# #generate next random angle using Julia's distributions package
				# persistant_d = VonMises(steps_angles[i-1],persistant_kappa)
				# next_angle = rand(sampler(persistant_d))
				# """

				# """
				# #generate the next random angle by calling R's wrapped normal distribution
				# prev_angle = persistant_angle
				# next_angle = R"rwrpnorm(1, $prev_angle, $p)"[1]
				# #convert the angle (R returns values between 0 and 2pi)
				# next_angle = (next_angle + pi)%(2*pi)-pi
				# """
			end
		end
	end

	#TODO: check if needs to be - accodring to orientaiton along the x and y axis
	next_x = prev_point[1] + next_length*cos(next_angle)
	next_y = prev_point[2] + next_length*sin(next_angle)

	next_point = (next_x,next_y)

	return (next_length,next_angle,next_point)
end

"""
simulate a random walk model with a 1D diffusion gradient
Input:
	n_steps : number of steps
	start_point : center of the cell
	cell_r : cell radios
	bias_angle : the angle of the wound
	params: model parameters:
							- w - weight of bias (w-1 is weight of persistance)
							- mb - scaling parameter for bias
							- b0 - baseline observed bias
							- mp - scaling parameter for persistance
							- p0 - baseline observed persistance

							- A - concentration strength
							- D - diffusion coefficient
							- tau0 - time the wound is signalling

							- Kd - receptor binding parameter (dissociation constant)
							- R0 - number of receptors
"""
function simulate_random_walk_attractant_diffusion_gradient(n_steps::Int64,start_point::Tuple{Float64,Float64},cell_r::Float64,bias_angle::Float64,params::Array{Float64,1};start_time::Int64=0.0)

	#save two result vectors: step lengths and step angles
	steps_length = Array{Float64,1}(n_steps)
	steps_angles = Array{Float64,1}(n_steps)
	steps_points = Array{Tuple{Float64,Float64},1}(n_steps+1)
	steps_points[1]=start_point

	#save the ob and op values
	bias=Array{Float64,1}()
    persistance=Array{Float64,1}()

	#save parameters
	(w,mb,b0,mp,p0,A,D,tau0,Kd,R0) = params

	#step length distribution
	#TODO: fit a better one to the data
	dt = 0.001
	truncated_normal = Truncated(Normal(2,1), 0, Inf)
	#gamma_dist = Gamma(0.75,2.5) #control gamma fit
	#gamma_dist = Gamma(1.0,2.5) #DMSO gamma fit
	#gamma_dist = Gamma(1.0,3.0) #p38 gamma fit
	#gamma_dist = Gamma(0.7,3.2) #JNK gamma fit
	#gamma_dist = Gamma(0.9,3.5) #green data
	#gamma_dist = Gamma(0.65,3.5) #red data

	#each step is 30 seconds of time
	current_time = start_time

	for i in 1:n_steps

		current_time += 20.0 #each step is 20 seconds

		# if current_time>50400
		# 	println("time exceeded : ", current_time)
		# end

		#get next step size
		#next_length = rand(truncated_normal) * sqrt(dt)
		next_length = rand(truncated_normal)
		#next_length = rand(gamma_dist)

		steps_length[i] = next_length

		#randomly choose if the next step is going to be biased or persistent
		rand_w = rand()

		#the distance from the wound is the distance along the y position
		if i==1
			y=start_point[2]
		else
			y=steps_points[i-1][2]
		end

		#get concentration at front and back
		(af,ar) = diffusion_gradient_solver(tau0,D,A,current_time,y,cell_r)

		#get internal concentration of the cell
		Raf = receptor_attractant_conc(Kd,R0,af)
		Rar = receptor_attractant_conc(Kd,R0,ar)

		if (Raf-Rar)<0
		 	#DEBUG
		 	#println("NEGATIVE")
		end

		if rand_w<=w

			#observed bias
			ob = mb*(Raf-Rar)+b0
			#println(ob)
			push!(bias,ob)

			 if ob<0.0 || ob>1.0
			 	#the parameters should be rejected, so no walk will be generated
			 	#println("ob is ", ob, " -> discard params")
			 	return false
			 end

			#bias_angle = convert_y_angle_to_x(bias_angle)

			if ob==0.0
				#sample angle from wrapped uniform distribution
				next_angle = pi*rand(Uniform(-1,1))
			elseif ob==1.0
				#next angle is the bias angle
				next_angle = bias_angle
			else
				#get the next angle according to the observed bias and the von mises distribution
				bias_kappa = 1/(-2*log(ob))
				next_angle = py_rand.vonmises(bias_angle,bias_kappa)
			end

		else
			#observed persistant
			op = mp*(Raf-Rar)+p0
			push!(persistance,op)

			 if op<0.0 || op>1.0
			 	#the parameters should be rejected, so no walk will be generated
			 	#println("op is ", op, " -> discard params")
			 	return false
			 end

			if i==1
				prev_angle = pi*rand(Uniform(-1,1))
			else
				prev_angle = steps_angles[i-1]
			end

			#prev_angle = convert_y_angle_to_x(prev_angle)

			if op==0.0
				#sample angle from wrapped uniform distribution
				next_angle = pi*rand(Uniform(-1,1))
			elseif op==1.0
				#next angle is the bias angle
				next_angle = prev_angle
			else
				#get the next angle according to the observed bias and the von mises distribution
				persistant_kappa = 1/(-2*log(op))
				next_angle = py_rand.vonmises(prev_angle,persistant_kappa)
			end

		end

		#convert angle to the negative y-axis
		#TODO: check correctness!!!

		#steps_angles[i] = convert_x_angle_to_y(next_angle)
		steps_angles[i] = next_angle

		#get the next center position of the cell ( according to the angle with the x-axis)
		#TODO: check correctness!!!

		if i==1
			prev_point = start_point
		else
			prev_point = steps_points[i-1]
		end

		next_x = prev_point[1] + (next_length)*cos(next_angle)
		next_y = prev_point[2] + (next_length)*sin(next_angle)

		next_point = (next_x,next_y)
		steps_points[i+1] = next_point

		if next_y<0
			#simulation has reached the y-position of the wound, stop simulations
			return Random_Walk(steps_length[1:i],steps_angles[1:i],steps_points[1:i],bias,persistance)
		end

	end

	return Random_Walk(steps_length,steps_angles,steps_points,bias,persistance)

end

"""
return the internal concentration of the cell
Input:
		- Kd - receptor binding parameter (dissociation constant)
		- R0 - number of receptors
		- ac - attractant concentration
"""
function receptor_attractant_conc(Kd::Float64,R0::Float64,ac::Float64)

	return 0.5*(Kd+R0+ac) - sqrt(0.25*((Kd+R0+ac)^2)-(R0*ac))
end

"""
convert the angle with the x-axis to an angle with the negative y-axis
"""
function convert_x_angle_to_y(x_angle::Float64)

	#TODO: FIX!

	y_angle = x_angle
	if y_angle>=0
		if y_angle<=pi/2
			y_angle += pi/2
		else
			y_angle = -x_angle
		end
	else
		if y_angle>=-pi/2
			y_angle = -x_angle
		else
			y_angle += pi/2
		end
	end

	return y_angle

end

"""
convert an angle with the negative y-axis to the x-axis
"""
function convert_y_angle_to_x(y_angle::Float64)

	#TODO: FIX!

	x_angle = y_angle
	if x_angle>=0

		if x_angle<=pi/2
			x_angle = -y_angle
		else
			x_angle -= pi/2
		end

	else

		if x_angle>=-pi/2
			x_angle -= pi/2
		else
			x_angle = -y_angle
		end

	end

	return x_angle

end


text:
#############################################################################
# SCS.jl
# Wrapper around the SCS solver https://github.com/cvxgrp/scs
# See http://github.com/JuliaOpt/SCS.jl
#############################################################################
# test/options.jl
# Tests the ability to pass options
#############################################################################

using MathProgBase

# The normal test
A = [1.0 1.0 0.0 0.0 0.0;
     0.0 1.0 0.0 0.0 1.0;
     0.0 0.0 1.0 1.0 1.0]
collb = [0.0, 0.0, 0.0, 0.0, 0.0]
obj   = [3.0, 4.0, 4.0, 9.0, 5.0]
rowub = [ 5.0,  3.0,  9.0]
s = SCSSolver()
m = MathProgBase.ConicModel(s)
MathProgBase.loadproblem!(m, -obj, A, rowub, [(:NonNeg,1:3)],[(:NonNeg,1:5)])
MathProgBase.optimize!(m)

@test isapprox(MathProgBase.getobjval(m), -99.0, atol=1e-6, rtol=0.0)
@test !isapprox(MathProgBase.getobjval(m), -99.0, atol=1e-7, rtol=0.0)

# With eps = 1e-12, solution should be far more accurate
s = SCSSolver(eps=1e-12)
m = MathProgBase.ConicModel(s)
MathProgBase.loadproblem!(m, -obj, A, rowub, [(:NonNeg,1:3)],[(:NonNeg,1:5)])
MathProgBase.optimize!(m)
@test isapprox(MathProgBase.getobjval(m), -99.0, atol=1e-10, rtol=0.0)

# With a warmstart from the eps = 1e-8 solution, solution should be extremely accurate even after 1 iteration
SCS.addoption!(m, :warm_start, true)
SCS.addoption!(m, :max_iters, 1)
MathProgBase.optimize!(m)
@test isapprox(MathProgBase.getobjval(m), -99.0, atol=1e-10, rtol=0.0)

# Now let's do the same warmstart, but on a new instance of the same problem
primal_sol = m.primal_sol
dual_sol = m.dual_sol
slack = m.slack
s = SCSSolver(max_iters=1)
m = MathProgBase.ConicModel(s)
MathProgBase.loadproblem!(m, -obj, A, rowub, [(:NonNeg,1:3)],[(:NonNeg,1:5)])
MathProgBase.setwarmstart!(m, primal_sol; dual_sol = dual_sol, slack = slack)
MathProgBase.optimize!(m)
@test isapprox(MathProgBase.getobjval(m), -99.0, atol=1e-10, rtol=0.0)

# tests for incorrect options
s = SCSSolver(eps=1e-12, epps=1.0)
m = MathProgBase.ConicModel(s)
MathProgBase.loadproblem!(m, -obj, A, rowub, [(:NonNeg,1:3)],[(:NonNeg,1:5)])

@test_throws ArgumentError MathProgBase.optimize!(m)

err = try
    MathProgBase.optimize!(m)
catch ex
    ex
end
@test err.msg == "Unrecognized option passed to SCS: epps;
Recognized options are: linear_solver, normalize, scale, rho_x, max_iters, eps, alpha, cg_rate, verbose, warm_start and acceleration_lookback."

# tests for incorrect options
s = SCSSolver(linear_solver="AAA", eps=1e-12, epps=1.0)
m = MathProgBase.ConicModel(s)
MathProgBase.loadproblem!(m, -obj, A, rowub, [(:NonNeg,1:3)],[(:NonNeg,1:5)])

@test_throws ArgumentError MathProgBase.optimize!(m)

err = try
    MathProgBase.optimize!(m)
catch ex
    ex
end

@test err.msg == "Unrecognized linear_solver passed to SCS: AAA;\nRecognized options are: SCS.Direct, SCS.Indirect."


text:
struct DistanceSphereToBox{M,R1,S,m} <: TaskMap{M,R1,S}
    radius::S
    box::Rect{m,S}
    rotation::RotMatrix{m,S}
end
DistanceSphereToBox{M,R1}(r, b::Rect{m,S}, R) where {M,R1,S,m} =
    DistanceSphereToBox{M,R1,S,dim(M)}(r, b, R)

struct DistanceSphereToBoxT{M,R1,S,m} <: TaskMapT{M,R1,S}
    base_map::DistanceSphereToBox{M,R1,S,m}
end
DistanceSphereToBoxT{M,R1}(r, b::Rect{m,S}, R) where {M,R1,S,m} =
    DistanceSphereToBoxT{M,R1,S,m}(DistanceSphereToBox{M,R1}(r, b, R))

function min_dist_dim(rect::Rect{N,T1}, p::Vec{N,T2}, dim::Int) where {N,T1,T2}
    T = promote_type(T1,T2)
    return max(zero(T), max(minimum(rect)[dim] - p[dim], p[dim] - maximum(rect)[dim]))
end

function min_euclideansq(rect::Rect{N,T1}, p::Union{Rect{N,T2},Vec{N,T2}}) where {N,T1,T2}
    T = promote_type(T1,T2)
    minimum_dist = T(0.0)
    for dim in 1:length(p)
        d = min_dist_dim(rect, p, dim)
        minimum_dist += d*d
    end
    return minimum_dist
end

function min_euclidean(rect::Rect{N,T1}, p::Union{Rect{N,T2},Vec{N,T2}}) where {N,T1,T2}
    return sqrt(min_euclideansq(rect, p))
end

function task_map_emb(::EmbRep, ::EmbRep, xme, task_map::DistanceSphereToBox{ℝ{m},R1}) where m
    box = task_map.box
    center = @. box.origin + box.widths/2
    xme_rectframe  = task_map.rotation' * (xme - center) + center
    SA[min_euclidean(box, Vec(xme_rectframe)) - task_map.radius]
end

# Updateable version
struct DistanceSphereToBoxDynamic{M,R1,S,m} <: TaskMap{M,R1,S}
    radius::S
    box::Vector{Rect{m,S}}
    rotation::Vector{RotMatrix{m,S}}
end
DistanceSphereToBoxDynamic{M,R1}(r, b::Vector{Rect{m,S}}, R) where {M,R1,S,m} =
    DistanceSphereToBoxDynamic{M,R1,S,dim(M)}(r, b, R)

struct DistanceSphereToBoxDynamicT{M,R1,S,m} <: TaskMapT{M,R1,S}
    base_map::DistanceSphereToBoxDynamic{M,R1,S,m}
end
DistanceSphereToBoxDynamicT{M,R1}(r, b::Vector{Rect{m,S}}, R) where {M,R1,S,m} =
    DistanceSphereToBoxDynamicT{M,R1,S,m}(DistanceSphereToBoxDynamic{M,R1}(r, b, R))

function task_map_emb(::EmbRep, ::EmbRep, xme,
        task_map::DistanceSphereToBoxDynamic{ℝ{m},R1}) where m
    box = task_map.box[1]
    center = @. box.origin + box.widths/2
    xme_rectframe  = task_map.rotation[1]' * (xme - center) + center
    SA[min_euclidean(box, Vec(xme_rectframe)) - task_map.radius]
end