text:
using Test
using QuantumOpticsBase
using Random, SparseArrays, LinearAlgebra

# Custom operator type for testing error msg
mutable struct TestOperator{BL<:Basis,BR<:Basis} <: AbstractOperator{BL,BR}; end

@testset "operators-sparse" begin

Random.seed!(0)

D(op1::AbstractOperator, op2::AbstractOperator) = abs(tracedistance_nh(dense(op1), dense(op2)))
D(x1::StateVector, x2::StateVector) = norm(x2-x1)
sprandop(b1, b2) = sparse(randoperator(b1, b2))
sprandop(b) = sprandop(b, b)

b1a = GenericBasis(2)
b1b = GenericBasis(3)
b2a = GenericBasis(1)
b2b = GenericBasis(4)
b3a = GenericBasis(1)
b3b = GenericBasis(5)

b_l = b1a⊗b2a⊗b3a
b_r = b1b⊗b2b⊗b3b

# Test creation
@test_throws DimensionMismatch DenseOperator(b1a, spzeros(ComplexF64, 3, 2))
@test_throws DimensionMismatch DenseOperator(b1a, b1b, spzeros(ComplexF64, 3, 2))
op1 = SparseOperator(b1a, b1b, sparse([1 1 1; 1 1 1]))
op2 = sparse(DenseOperator(b1b, b1a, [1 1; 1 1; 1 1]))
@test op1 == dagger(op2)

# Test transpose
@test transpose(op2) == conj(op2')

# Test copy
op1 = sparse(randoperator(b1a))
op2 = copy(op1)
@test !(op1.data === op2.data)
op2.data[1,1] = complex(10.)
@test op1.data[1,1] != op2.data[1,1]

# Arithmetic operations
# =====================
op_zero = SparseOperator(b_l, b_r)
op1 = sprandop(b_l, b_r)
op2 = sprandop(b_l, b_r)
op3 = sprandop(b_l, b_r)
op1_ = dense(op1)
op2_ = dense(op2)
op3_ = dense(op3)

x1 = randstate(b_r)
x2 = randstate(b_r)

xbra1 = dagger(randstate(b_l))
xbra2 = dagger(randstate(b_l))

# Addition
@test_throws DimensionMismatch op1 + dagger(op2)
@test 1e-14 > D(op1+op2, op1_+op2_)
@test 1e-14 > D(op1+op2, op1+op2_)
@test 1e-14 > D(op1+op2, op1_+op2)

# Subtraction
@test_throws DimensionMismatch op1 - dagger(op2)
@test 1e-14 > D(op1-op2, op1_-op2_)
@test 1e-14 > D(op1-op2, op1-op2_)
@test 1e-14 > D(op1-op2, op1_-op2)
@test 1e-14 > D(op1+(-op2), op1_ - op2_)
@test 1e-14 > D(op1+(-1*op2), op1_ - op2_)

# Test multiplication
@test_throws DimensionMismatch op1*op2
@test 1e-11 > D(op1*(x1 + 0.3*x2), op1_*(x1 + 0.3*x2))
@test 1e-11 > D(op1*x1 + 0.3*op1*x2, op1_*x1 + 0.3*op1_*x2)
@test 1e-11 > D((op1+op2)*(x1+0.3*x2), (op1_+op2_)*(x1+0.3*x2))

@test 1e-11 > D((xbra1 + 0.3*xbra2)*op1, (xbra1 + 0.3*xbra2)*op1_)
@test 1e-11 > D(xbra1*op1 + 0.3*xbra2*op1, xbra1*op1_ + 0.3*xbra2*op1_)
@test 1e-11 > D((xbra1+0.3*xbra2)*(op1+op2), (xbra1+0.3*xbra2)*(op1_+op2_))

@test 1e-12 > D(op1*dagger(0.3*op2), op1_*dagger(0.3*op2_))
@test 1e-12 > D(0.3*dagger(op2*dagger(op1)), 0.3*dagger(op2_*dagger(op1_)))
@test 1e-12 > D((op1 + op2)*dagger(0.3*op3), (op1_ + op2_)*dagger(0.3*op3_))
@test 1e-12 > D(0.3*op1*dagger(op3) + 0.3*op2*dagger(op3), 0.3*op1_*dagger(op3_) + 0.3*op2_*dagger(op3_))

# Test division
@test 1e-14 > D(op1/7, op1_/7)

# Conjugation
tmp = copy(op1)
conj!(tmp)
@test tmp == conj(op1) && conj(tmp.data) == op1.data

# Test identityoperator
Idense = identityoperator(DenseOpType, b_r)
I = identityoperator(SparseOpType, b_r)
@test isa(I, SparseOpType)
@test dense(I) == Idense
@test 1e-11 > D(I*x1, x1)
@test I == identityoperator(SparseOpType, b1b) ⊗ identityoperator(SparseOpType, b2b) ⊗ identityoperator(SparseOpType, b3b)

Idense = identityoperator(DenseOpType, b_l)
I = identityoperator(SparseOpType, b_l)
@test isa(I, SparseOpType)
@test dense(I) == Idense
@test 1e-11 > D(xbra1*I, xbra1)
@test I == identityoperator(SparseOpType, b1a) ⊗ identityoperator(SparseOpType, b2a) ⊗ identityoperator(SparseOpType, b3a)

# Test tr and normalize
op = sparse(DenseOperator(GenericBasis(3), [1 3 2;5 2 2;-1 2 5]))
@test 8 == tr(op)
op_normalized = normalize(op)
@test 8 == tr(op)
@test 1 == tr(op_normalized)
# op_ = normalize!(op)
# @test op_ === op
# @test 1 == tr(op)

# Test partial tr
b1 = GenericBasis(3)
b2 = GenericBasis(5)
b3 = GenericBasis(7)
b_l = b1 ⊗ b2 ⊗ b3
op1 = sprandop(b1)
op2 = sprandop(b2)
op3 = sprandop(b3)
op123 = op1 ⊗ op2 ⊗ op3
op123_ = dense(op123)

@test 1e-14 > D(ptrace(op123_, 3), ptrace(op123, 3))
@test 1e-14 > D(ptrace(op123_, 2), ptrace(op123, 2))
@test 1e-14 > D(ptrace(op123_, 1), ptrace(op123, 1))

@test 1e-14 > D(ptrace(op123_, [2,3]), ptrace(op123, [2,3]))
@test 1e-14 > D(ptrace(op123_, [1,3]), ptrace(op123, [1,3]))
@test 1e-14 > D(ptrace(op123_, [1,2]), ptrace(op123, [1,2]))

@test_throws ArgumentError ptrace(op123, [1,2,3])

# Test expect
state = randstate(b_l)
@test expect(op123, state) ≈ expect(op123_, state)

state = randoperator(b_l)
@test expect(op123, state) ≈ expect(op123_, state)

@test_throws QuantumOpticsBase.IncompatibleBases expect(op1, op2)

# Tensor product
# ==============
b1a = GenericBasis(2)
b1b = GenericBasis(3)
b2a = GenericBasis(1)
b2b = GenericBasis(4)
b3a = GenericBasis(1)
b3b = GenericBasis(5)
b_l = b1a ⊗ b2a ⊗ b3a
b_r = b1b ⊗ b2b ⊗ b3b
op1a = sprandop(b1a, b1b)
op1b = sprandop(b1a, b1b)
op2a = sprandop(b2a, b2b)
op2b = sprandop(b2a, b2b)
op3a = sprandop(b3a, b3b)
op1a_ = dense(op1a)
op1b_ = dense(op1b)
op2a_ = dense(op2a)
op2b_ = dense(op2b)
op3a_ = dense(op3a)
op123 = op1a ⊗ op2a ⊗ op3a
op123_ = op1a_ ⊗ op2a_ ⊗ op3a_
@test op123.basis_l == b_l
@test op123.basis_r == b_r

# Associativity
@test 1e-13 > D((op1a ⊗ op2a) ⊗ op3a, (op1a_ ⊗ op2a_) ⊗ op3a_)
@test 1e-13 > D(op1a ⊗ (op2a ⊗ op3a), op1a_ ⊗ (op2a_ ⊗ op3a_))
@test 1e-13 > D(op1a ⊗ (op2a ⊗ op3a), op1a_ ⊗ (op2a_ ⊗ op3a))

# Linearity
@test 1e-13 > D(op1a ⊗ (0.3*op2a), op1a_ ⊗ (0.3*op2a_))
@test 1e-13 > D(0.3*(op1a ⊗ op2a), 0.3*(op1a_ ⊗ op2a_))
@test 1e-13 > D((0.3*op1a) ⊗ op2a, (0.3*op1a_) ⊗ op2a_)
@test 1e-13 > D(0.3*(op1a ⊗ op2a), 0.3*(op1a_ ⊗ op2a_))
@test 1e-13 > D(0.3*(op1a ⊗ op2a), 0.3*(op1a ⊗ op2a_))

# Distributivity
@test 1e-13 > D(op1a ⊗ (op2a + op2b), op1a_ ⊗ (op2a_ + op2b_))
@test 1e-13 > D(op1a ⊗ op2a + op1a ⊗ op2b, op1a_ ⊗ op2a_ + op1a_ ⊗ op2b_)
@test 1e-13 > D((op2a + op2b) ⊗ op3a, (op2a_ + op2b_) ⊗ op3a_)
@test 1e-13 > D(op2a ⊗ op3a + op2b ⊗ op3a, op2a_ ⊗ op3a_ + op2b_ ⊗ op3a_)
@test 1e-13 > D(op2a ⊗ op3a + op2b ⊗ op3a, op2a ⊗ op3a_ + op2b_ ⊗ op3a_)

# Mixed-product property
@test 1e-13 > D((op1a ⊗ op2a) * dagger(op1b ⊗ op2b), (op1a_ ⊗ op2a_) * dagger(op1b_ ⊗ op2b_))
@test 1e-13 > D((op1a*dagger(op1b)) ⊗ (op2a*dagger(op2b)), (op1a_*dagger(op1b_)) ⊗ (op2a_*dagger(op2b_)))
@test 1e-13 > D((op1a*dagger(op1b)) ⊗ (op2a*dagger(op2b)), (op1a_*dagger(op1b)) ⊗ (op2a_*dagger(op2b_)))

# Transpose
@test 1e-13 > D(dagger(op1a ⊗ op2a), dagger(op1a_ ⊗ op2a_))
@test 1e-13 > D(dagger(op1a ⊗ op2a), dagger(op1a ⊗ op2a_))
@test 1e-13 > D(dagger(op1a) ⊗ dagger(op2a), dagger(op1a_) ⊗ dagger(op2a_))

# Permute systems
op1 = sprandop(b1a, b1b)
op2 = sprandop(b2a, b2b)
op3 = sprandop(b3a, b3b)
op123 = op1⊗op2⊗op3

op132 = op1⊗op3⊗op2
@test 1e-14 > D(permutesystems(op123, [1, 3, 2]), op132)

op213 = op2⊗op1⊗op3
@test 1e-14 > D(permutesystems(op123, [2, 1, 3]), op213)

op231 = op2⊗op3⊗op1
@test 1e-14 > D(permutesystems(op123, [2, 3, 1]), op231)

op312 = op3⊗op1⊗op2
@test 1e-14 > D(permutesystems(op123, [3, 1, 2]), op312)

op321 = op3⊗op2⊗op1
@test 1e-14 > D(permutesystems(op123, [3, 2, 1]), op321)

# Test diagonaloperator
b = GenericBasis(4)
I = identityoperator(b)

@test diagonaloperator(b, [1, 1, 1, 1]) == I
@test diagonaloperator(b, [1., 1., 1., 1.]) == I
@test diagonaloperator(b, [1im, 1im, 1im, 1im]) == 1im*I
@test diagonaloperator(b, [0:3;]) == sparse(DenseOperator(b, Diagonal([0:3;])))

# Test gemv
op = sprandop(b_l, b_r)
op_ = dense(op)
xket = randstate(b_l)
xbra = dagger(xket)

state = randstate(b_r)
result_ = randstate(b_l)
result = deepcopy(result_)
QuantumOpticsBase.mul!(result,op,state,complex(1.0),complex(0.))
@test 1e-13 > D(result, op_*state)

result = deepcopy(result_)
alpha = complex(1.5)
beta = complex(2.1)
QuantumOpticsBase.mul!(result,op,state,alpha,beta)
@test 1e-13 > D(result, alpha*op_*state + beta*result_)

state = dagger(randstate(b_l))
result_ = dagger(randstate(b_r))
result = deepcopy(result_)
QuantumOpticsBase.mul!(result,state,op,complex(1.0),complex(0.))
@test 1e-13 > D(result, state*op_)

result = deepcopy(result_)
alpha = complex(1.5)
beta = complex(2.1)
QuantumOpticsBase.mul!(result,state,op,alpha,beta)
@test 1e-13 > D(result, alpha*state*op_ + beta*result_)

# Test gemm small version
b1 = GenericBasis(3)
b2 = GenericBasis(5)
b3 = GenericBasis(7)

op = sprandop(b1, b2)
op_ = dense(op)

state = randoperator(b2, b3)
result_ = randoperator(b1, b3)
result = deepcopy(result_)
QuantumOpticsBase.mul!(result,op,state,complex(1.),complex(0.))
@test 1e-12 > D(result, op_*state)

result = deepcopy(result_)
alpha = complex(1.5)
beta = complex(2.1)
QuantumOpticsBase.mul!(result,op,state,alpha,beta)
@test 1e-12 > D(result, alpha*op_*state + beta*result_)

state = randoperator(b3, b1)
result_ = randoperator(b3, b2)
result = deepcopy(result_)
QuantumOpticsBase.mul!(result,state,op,complex(1.),complex(0.))
@test 1e-12 > D(result, state*op_)

result = deepcopy(result_)
alpha = complex(1.5)
beta = complex(2.1)
QuantumOpticsBase.mul!(result,state,op,alpha,beta)
@test 1e-12 > D(result, alpha*state*op_ + beta*result_)

state = randoperator(b1, b3)
result_ = randoperator(b2, b3)
result = deepcopy(result_)
result = deepcopy(result_)
alpha = complex(1.5)
beta = complex(2.1)
QuantumOpticsBase.mul!(result,op',state,alpha,beta) # gemm! with lazy adjoint sparse
@test 1e-12 > D(result, alpha*op_'*state + beta*result_)

state = randoperator(b1, b2)
result_ = randoperator(b1, b1)
result = deepcopy(result_)
result = deepcopy(result_)
alpha = complex(1.5)
beta = complex(2.1)
QuantumOpticsBase.mul!(result,state,op',alpha,beta) # gemm! with lazy adjoint sparse
@test 1e-12 > D(result, alpha*state*op_' + beta*result_)

# Test gemm big version
b1 = GenericBasis(50)
b2 = GenericBasis(60)
b3 = GenericBasis(55)

op = sprandop(b1, b2)
op_ = dense(op)

state = randoperator(b2, b3)
result_ = randoperator(b1, b3)
result = deepcopy(result_)
QuantumOpticsBase.mul!(result,op,state,complex(1.),complex(0.))
@test 1e-11 > D(result, op_*state)

result = deepcopy(result_)
alpha = complex(1.5)
beta = complex(2.1)
QuantumOpticsBase.mul!(result,op,state,alpha,beta)
@test 1e-11 > D(result, alpha*op_*state + beta*result_)

state = randoperator(b3, b1)
result_ = randoperator(b3, b2)
result = deepcopy(result_)
QuantumOpticsBase.mul!(result,state,op,complex(1.),complex(0.))
@test 1e-11 > D(result, state*op_)

result = deepcopy(result_)
alpha = complex(1.5)
beta = complex(2.1)
QuantumOpticsBase.mul!(result,state,op,alpha,beta)
@test 1e-11 > D(result, alpha*state*op_ + beta*result_)

state = randoperator(b1, b2)
result_ = randoperator(b1, b1)
result = deepcopy(result_)
result = deepcopy(result_)
alpha = complex(1.5)
beta = complex(2.1)
QuantumOpticsBase.mul!(result,state,op',alpha,beta) # gemm! with lazy adjoint sparse
@test 1e-11 > D(result, alpha*state*op_' + beta*result_)

# Test remaining uncovered code
@test_throws DimensionMismatch SparseOperator(b1, b2, zeros(10, 10))
dat = sprandop(b1, b1).data
@test SparseOperator(b1, dat) == SparseOperator(b1, Matrix{ComplexF64}(dat))

@test_throws ArgumentError sparse(TestOperator{Basis,Basis}())

@test 2*SparseOperator(b1, dat) == SparseOperator(b1, dat)*2
@test copy(op1) == deepcopy(op1)

# Test Hermitian
bspin = SpinBasis(1//2)
bnlevel = NLevelBasis(2)
@test ishermitian(SparseOperator(bspin, bspin, sparse([1.0 im; -im 2.0]))) == true
@test ishermitian(SparseOperator(bspin, bnlevel, sparse([1.0 im; -im 2.0]))) == false

# Test broadcasting
@test_throws DimensionMismatch op1 .+ op2
@test op1 .+ op1 == op1 + op1
op1 .= DenseOperator(op1)
@test isa(op1, SparseOpType)
# @test isa(op1 .+ DenseOperator(op1), DenseOpType) # Broadcasting of sparse .+ dense matrix results in sparse
op3 = sprandop(FockBasis(1),FockBasis(2))
@test_throws QuantumOpticsBase.IncompatibleBases op1 .+ op3
@test_throws QuantumOpticsBase.IncompatibleBases op1 .= op3
op_ = copy(op1)
op_ .+= op1
@test op_ == 2*op1
@test_throws ErrorException cos.(op_)

end # testset


text:
using StatsBase: var
function solver_default(
    X,
    data::AbstractArray{T,2};
    stderror = false,
) where {T<:Union{Missing,<:Number}}
    minfo = []
    sizehint!(minfo, size(data,1))
    beta = zeros(size(data,1),size(X,2)) # had issues with undef
    for ch = 1:size(data, 1)
	 	dd = view(data, ch, :)
        ix = @. !ismissing(dd)
	    # use the previous channel as a starting point
        ch == 1 || copyto!(view(beta, ch, :), view(beta, ch-1, :))

		beta[ch,:],h = lsmr!(@view(beta[ch, :]), (X[ix,:]), @view(data[ch, ix]),log=true)

        push!(minfo, h)
    end

    if stderror
        stderror = calculate_stderror(X, data, beta)
        modelfit = Unfold.LinearModelFit(beta, ["lsmr",minfo], stderror)
    else
        modelfit = Unfold.LinearModelFit(beta, ["lsmr",minfo])
    end
    return modelfit
end

function calculate_stderror(Xdc, data::Matrix{T}, beta) where {T<:Union{Missing,<:Number}}

    # remove missings
    ix = any(.!ismissing.(data), dims = 1)[1, :]
    if length(ix) != size(data, 2)
        @warn(
            "Limitation: Missing data are calculated over all channels for standard error"
        )
    end
    data = data[:, ix]
    Xdc = Xdc[ix, :]

    # Hat matrix only once
    hat_prime = inv(Matrix(Xdc' * Xdc))
    # Calculate residual variance
    @warn(
        "Autocorrelation was NOT taken into account. Therefore SE are UNRELIABLE. Use at your own discretion"
    )

    se = Array{Float64}(undef, size(data, 1), size(Xdc, 2))
    for ch = 1:size(data, 1)
        residualVar = var(data[ch, :] .- Xdc * beta[ch, :])
        @assert(!isnan(residualVar), "residual Variance was NaN")
        hat = hat_prime .* residualVar
        #se = sqrt(diag(cfg.contrast(:,:)*hat*cfg.contrast(:,:)'));
        se[ch, :] = sqrt.(diag(hat))
    end
    return se
end
function calculate_stderror(
    X,
    data::AbstractArray{T,3},
    beta,
) where {T<:Union{Missing,<:Number}}
    #function calculate_stderror(Xdc,data::AbstractArray{T,2},beta) where {T<:Union{Missing, <:Number}}

    # Hat matrix
    hat_prime = inv(Matrix(X' * X))
    # Calculate residual variance
    @warn(
        "Autocorrelation was NOT taken into account. Therefore SE are UNRELIABLE. Use at your own discretion"
    )

    se = Array{Float64}(undef, size(data, 1), size(data, 2), size(X, 2))
    for ch = 1:size(data, 1)
        for t = 1:size(data, 2)
            ix = .!ismissing.(data[ch, t, :])
            residualVar = var(data[ch, t, ix] .- X[ix, :] * beta[ch, t, :])
            @assert(!isnan(residualVar), "residual Variance was NaN")
            hat = hat_prime .* residualVar
            #se = sqrt(diag(cfg.contrast(:,:)*hat*cfg.contrast(:,:)'));
            se[ch, t, :] = sqrt.(diag(hat))
        end
    end
    return se
end
function solver_default(
    X,
    data::AbstractArray{T,3};
    stderror = false,
) where {T<:Union{Missing,<:Number}}
    #beta = Array{Union{Missing,Number}}(undef, size(data, 1), size(data, 2), size(X, 2))
    beta = zeros(Union{Missing,Number},size(data, 1), size(data, 2), size(X, 2))
    @showprogress 0.1 for ch = 1:size(data, 1)
        for t = 1:size(data, 2)
            @debug("$(ndims(data,)),$t,$ch")

            dd = view(data, ch, t,:)
            ix = @. !ismissing(dd)

            beta[ch,t,:] = @view(X[ix,:]) \ @view(data[ch,t,ix])
            # qr(X) was slower on Februar 2022
        end
    end

    if stderror
        stderror = calculate_stderror(X, data, beta)
        modelfit = LinearModelFit(beta, ["solver_default"], stderror)
    else
        modelfit = LinearModelFit(beta, ["solver_default"])
    end
    return modelfit
end

solver_b2b(X, data, cross_val_reps) = solver_b2b(X, data, cross_val_reps = cross_val_reps)
function solver_b2b(
    X,
    data::AbstractArray{T,3};
    cross_val_reps = 10,
) where {T<:Union{Missing,<:Number}}

    X, data = dropMissingEpochs(X, data)

    E = zeros(size(data, 2), size(X, 2), size(X, 2))
    W = Array{Float64}(undef, size(data, 2), size(X, 2), size(data, 1))

    prog = Progress(size(data, 2) * cross_val_reps, 0.1)
    for t = 1:size(data, 2)

        for m = 1:cross_val_reps
            k_ix = collect(Kfold(size(data, 3), 2))
            Y1 = @view data[:, t, k_ix[1]]
            Y2 = @view data[:, t, k_ix[2]]
            X1 = @view X[k_ix[1], :]
            X2 = @view X[k_ix[2], :]

            G = (Y1' \ X1)
            H = X2 \ (Y2' * G)

            E[t, :, :] +=  Diagonal(H[diagind(H)])
            ProgressMeter.next!(prog; showvalues = [(:time, t), (:cross_val_rep, m)])
        end
        E[t, :, :] = E[t, :, :] ./ cross_val_reps
        W[t, :, :] = (X * E[t, :, :])' / data[:, t, :]

    end

    # extract diagonal
    beta = mapslices(diag, E, dims = [2, 3])
    # reshape to conform to ch x time x pred
    beta = permutedims(beta, [3 1 2])
    modelinfo = Dict("W" => W, "E" => E, "cross_val_reps" => cross_val_reps) # no history implemented (yet?)
    return LinearModelFit(beta, modelinfo)
end

function solver_robust(
    X,
    data::AbstractArray{T,3};
	estimator = MEstimator{TukeyLoss}(),
    rlmOptions = (initial_scale=:mad,)
) where {T<:Union{Missing,<:Number}}
    #beta = zeros(Union{Missing,Number},size(data, 1), size(data, 2), size(X, 2))
	beta = zeros(T,size(data, 1), size(data, 2), size(X, 2))
    @showprogress 0.1 for ch = 1:size(data, 1)
        for t = 1:size(data, 2)
            #@debug("$(ndims(data,)),$t,$ch")

            dd = view(data, ch, t,:)
            ix = @. !ismissing(dd)
			# init beta
			#if ch == 1 && t==1
			#elseif ch > 1 && t == 1
			#	copyto!(view(beta, ch, 1,:), view(beta, ch-1, 1,:))
			#else
			# 	copyto!(view(beta, ch,t, :), view(beta, ch,t-1, :))
			#end

			X_local = disallowmissing((X[ix,:])) # view crashes robust model here. XXX follow up once
            # https://github.com/JuliaStats/GLM.jl/issues/470 received a satisfying result
			y_local = disallowmissing(@view(data[ch,t,ix]))

            # if not specified otherwise
            # this results in strange regularizing effects
			#if :initial_coef ∉ keys(rlmOptions)
            #    rlmOptions = merge(rlmOptions,(initial_coef=@view(beta[ch,t,:]),))
            #end
			m = rlm(X_local,y_local,estimator;rlmOptions...)
			beta[ch,t,:] .= coef(m)

        end
    end

	modelfit = Unfold.LinearModelFit(beta, ["solver_robust"])

    return modelfit
end

ext: jl

max_stars_repo_name: unfoldtoolbox/jlunfold

max_issues_repo_name: unfoldtoolbox/jlunfold

max_forks_repo_name: unfoldtoolbox/jlunfold


text:
macro load_examples()
  #Add some String Dict lookup
  Dict_example=quote
    function custom_completions{K<:AbstractString,V}(::typeof(getindex),a::Dict{K,V},s::AbstractString)
      return UTF8String[utf8(x) for x in keys(a)]
    end
  end

  #Add lookup of names in NetCDF variables
  NetCDF_example=quote
    using NetCDF
    function custom_completions(f::typeof(ncread),filename::AbstractString,varname::AbstractString)
      !isfile(filename) && return UTF8STring[]
      try
        nc=NetCDF.open(filename)
        return collect(keys(nc.vars))
      end
      return UTF8String[]
    end

    function custom_completions(f::typeof(ncgetatt),filename::AbstractString,varname::AbstractString)
      !isfile(filename) && return UTF8STring[]
      try
        nc=NetCDF.open(filename)
        return collect(keys(nc.vars))
      end
      return UTF8String[]
    end

    function custom_completions(f::typeof(ncgetatt),filename::AbstractString,varname::AbstractString,attname::AbstractString)
      !isfile(filename) && return UTF8STring[]
      try
        nc=NetCDF.open(filename)
        return collect(keys(nc.vars[varname].atts))
      end
      return UTF8String[]
    end
  end
  HDF5_example=quote
    #Add group and variable name lookup in HDF5
    using HDF5
    function custom_completions(f::typeof(h5read),filename::AbstractString,var::AbstractString)
      !isfile(filename) && return UTF8STring[]
      try
        #First try to split expression in var
        m=match(r"/",reverse(var))
        if m==nothing
          h=h5open(filename,"r")
          nlist=names(h)
          return UTF8String[isa(h[nn],HDF5.HDF5Group) ? nn*"/" : nn for nn in nlist]
        else
          h=h5open(filename,"r")
          sleft=var[1:(end-m.offset)]
          d=h[sleft]
          nlist=names(d)
          return UTF8String[isa(d[nn],HDF5.HDF5Group) ? nn*"/" : nn for nn in nlist]
        end
      end
      return UTF8String[]
    end
  end
  retex=Dict_example
  if Pkg.installed("NetCDF")!=nothing
    retex=quote
      $retex
      $NetCDF_example
    end
  end
  if Pkg.installed("HDF5")!=nothing
    retex=quote
      $retex
      $HDF5_example
    end
  end
  esc(retex)
end


text:
function gen_symmetric_matrix(n)
  a = randn(n,n)
  (a + a')/2
end

for n in [10, 100, 500, 1000, 2000]
  a = gen_symmetric_matrix(n)

  tic()
    de, ve = eig(a)
  teig = toq()

  k = 5
  while k <= n/2
    tic()
      des, ves = eigs(a; nev = k)
    teigs = toq()

    #dnorm = norm( de-des[1:k] )
    #vnorm = norm( ve-ves )

    println("$(n): k = $(k), eigd norm = , eigv norm = , eigs time factor = $(teig/teigs), eig time = $(teig), eigs time = $(teigs)")

    k *= 10
  end
end


text:
merge_ext!(S::SeisData, Ω::Int64, rest::Array{Int64, 1}) = nothing

function merge_ext!(C::T1, D::T2) where {T1<:GphysChannel, T2<:GphysChannel}
  if T1 == T2
    ff = setdiff(fieldnames(T1), SeisIO.datafields)
    for f in ff
      setfield!(C, f, deepcopy(getfield(D, f)))
    end
  end
  return nothing
end

function dup_check!(subgrp::Array{Int64, 1}, to_delete::Array{Int64, 1}, T::Array{Array{Int64, 2}, 1}, X::Array{FloatArray, 1})
  N = length(subgrp)
  if N > 1
    # Check for duplicates
    sort!(subgrp)
    u = falses(N)
    while N > 1
      t1 = getindex(T, getindex(subgrp, N))
      x1 = getindex(X, getindex(subgrp, N))

      # if a channel is already known to be a duplicate, skip it
      if getindex(u, N) == false
        j = N
        while j > 1
          j = j-1
          t2 = getindex(T, getindex(subgrp, j))
          x2 = getindex(X, getindex(subgrp, j))
          if t1 == t2 && x1 == x2
            setindex!(u, true, j)
          end
        end
      end

      N = N-1
    end

    # flag duplicates for deletion
    for i in 1:length(u)
      if u[i]
        push!(to_delete, subgrp[i])
      end
    end

    # remove duplicates from merge targets
    deleteat!(subgrp, u)
  end
  return length(subgrp)
end

function get_δt(t::Int64, Δ::Int64)
  δts = rem(t, Δ)
  if δts > div(Δ,2)
    δts -= Δ
  end
  return δts
end

function check_alignment(Ti::Array{Int64,1}, Tj::Array{Int64,1}, Xi::Array{T,1}, Xj::Array{T,1}, Δ::Int64) where {T<:AbstractFloat}
  z = zero(Int64)

  # Rare case, but possible: they're the same data
  if first(Ti) == first(Tj) && last(Ti) == last(Tj) && Xi == Xj
    return Ti, Xi, false, z
  else

    # Find a time correction δt that gets applied to Ti, ts[i], te[i], etc.
    # Find overlapping region, if one exists
    L = length(Ti)
    for n = 1:4
      if n > L-1
        break
      end

      xi_f = view(Xi, n+1:L)
      xj_f = view(Xj, 1:L-n)
      if isapprox(xi_f, xj_f)
        return Tj[1:L-n], xj_f, false, 1*n
      end

      xi_b = view(Xi, 1:L-n)
      xj_b = view(Xj, n+1:L)
      if isapprox(xi_b, xj_b)
        return Ti[1:L-n], xi_f, false, -1*n
      end
    end
  end
  return vcat(Ti,Tj), vcat(Xi,Xj), true, z
end

function xtmerge!(t::Array{Int64,1}, x::Array{T,1}, d::Int64) where {T<:AbstractFloat}
  # Sanity check
  (length(t) == length(x)) || error(string("Badly set times (Nt=", length(t), ",Nx=", length(x), "); can't merge!"))

  # Sort
  i = sortperm(t)
  sort!(t)
  x[:] = x[i]

  # Check for duplicates
  J0 = findall((diff(t).==0).*(diff(x).==0))
  while !isempty(J0)
    deleteat!(x, J0)
    deleteat!(t, J0)
    J0 = findall(diff(t) .== 0)
  end

  J0 = findall(diff(t) .< d)
  while !isempty(J0)
    J1 = J0.+1
    K = [isnan.(x[J0]) isnan.(x[J1])]

    # Average nearly-overlapping x that are either both NaN or neither Nan
    ii = findall(K[:,1].==K[:,2])
    i0 = J0[ii]
    i1 = J1[ii]
    t[i0] = div.(t[i0].+t[i1], 2)
    x[i0] = 0.5.*(x[i0].+x[i1])

    # Delete nearly-overlapping x with only one NaN (and delete all x ∈ i1)
    i3 = findall(K[:,1].*(K[:,2].==false))
    i4 = findall((K[:,1].==false).*K[:,2])
    II = sort([J0[i4]; J1[i3]; i1])
    deleteat!(t, II)
    deleteat!(x, II)

    J0 = findall(diff(t) .< d)
  end
  return nothing
end

# get hash of each non-empty loc, resp, units, fs ( fs should never be empty )
function get_subgroups( LOC   ::Array{InstrumentPosition,1},
                        FS    ::Array{Float64,1},
                        RESP  ::Array{InstrumentResponse,1},
                        UNITS ::Array{String,1},
                        group ::Array{Int64,1} )
  zh = zero(UInt64)
  N_grp = length(group)
  N_grp == 1 && return([group])

  H = Array{UInt,2}(undef, N_grp, 4)
  for (n,g) in enumerate(group)
    H[n,1] = hash(getindex(FS, g))
    H[n,2] = isempty(getindex(LOC, g)) ? zh : hash(getindex(LOC, g))
    H[n,3] = isempty(getindex(RESP, g)) ? zh : hash(getindex(RESP, g))
    H[n,4] = isempty(getindex(UNITS, g)) ? zh : hash(getindex(UNITS, g))
  end

  # If an entire column is unset, we don't care
  H = H[:, findall([sum(H[:,i])>0 for i in 1:4])]
  (Nh,Nc) = size(H)

  # Find unique rows of H; sort
  H_filled = sum(H .!= zh, dims=2)[:]
  H_inds = sortperm(H_filled, rev=true)
  N_subgrp = length(H_inds)
  N_subgrp == 1 && return([group])

  subgrp_inds = Array{Array{Int64,1},1}(undef, N_subgrp)

  H = H[H_inds, :]
  group = group[H_inds]
  H_sub = deepcopy(H)

  for i = 1:N_subgrp
    subgrp_inds[i] = Array{Int64,1}(undef, 0)
    subgrp_hash = H_sub[i, :]
    for j = N_grp:-1:1
      m = prod([H[j,k] in (zh, subgrp_hash[k]) for k=1:Nc])
      if m
        push!(subgrp_inds[i], group[j])
        deleteat!(group, j)
        H = H[setdiff(1:end, j), :]
      end
      N_grp = length(group)
    end
    if N_grp == 0
      subgrp_inds = subgrp_inds[1:i]
      N_subgrp = length(subgrp_inds)
      break
    end
  end
  deleteat!(subgrp_inds, [isempty(subgrp_inds[j]) for j=1:length(subgrp_inds)])
  return subgrp_inds
end

function get_next_pair(W::Array{Int64,2})
  L = size(W,1)
  i = 1

  # dest loop
  while i < L
    si = getindex(W, i)
    ei = getindex(W, i + L)
    j = i + 1

    # src loop
    while j ≤ L
      sj = getindex(W, j)
      ej = getindex(W, j+L)
      if min(si ≤ ej, ei ≥ sj) == true
        # src    # dest
        return vcat(W[j,:], j), vcat(W[i,:], i)
      end
      j = j + 1
    end
    i = i + 1
  end
  return zeros(Int64, 7), zeros(Int64, 7)
end

function get_merge_w(Δ::Int64, subgrp::Array{Int64,1}, T::Array{Array{Int64, 2}, 1}, X::Array{FloatArray, 1})

  N = length(subgrp)
  w_tmp = Array{Array{Int64,2}, 1}(undef, N)
  te    = Array{Int64, 1}(undef, N)
  for i = 1:N
    m     = getindex(subgrp, i)
    w_m   = t_win(getindex(T, m), Δ)
    n_w   = size(w_m, 1)

    # Store: w_start, w_end, channel_number, window_number, x_start, x_end
    w = hcat(w_m, Array{Int64,2}(undef, n_w, 4))
    j  = 0
    ws = 0
    we = 0
    while j < n_w
      j       = j + 1
      we      = max(we, getindex(w_m, j, 2))
      setindex!(w, m, j, 3)                     # channel number
      setindex!(w, j, j, 4)                     # window number
      setindex!(w, ws+1, j, 5)                  # x_start
      ws      = ws + div(w[j,2]-w[j,1], Δ)+1    # x_end
      setindex!(w, ws, j, 6)
    end
    setindex!(te, we, i)
    setindex!(w_tmp, w, i)
  end
  W   = vcat(w_tmp...)
  ii  = sortperm(W[:,2], rev=true)
  W   = W[ii,:]

  #= added 2019-11-19:
  fixes a rare off-by-one bug with slightly-offset windows (issue #29)
  issue creator only sees bug in one file from 20 years of data

  clumsy fix based on sound principles:
  * force X to start an integer # of samples from the epoch
    - prevents a discrepancy between length(X) and length(T)
  * add the offset back to the start time of the merged channel data
  =#

  # Let Ω be the channel number in subgrp with the last end time
  Ω = subgrp[argmax(te)]

  return W, Ω
end

function segment_merge(Δ::Int64, Ω::Int64, W::Array{Int64, 2}, X::Array{FloatArray, 1})
  nW = size(W,1)
  i = argmin(W[:,1])
  δts = get_δt(W[i,1], Δ)
  if δts != 0
    for i in 1:size(W,1)
      W[i,1] -= δts
      W[i,2] -= δts
    end
  end
  (src, dest) = get_next_pair(W)

  while (src, dest) != (zeros(Int64, 7), zeros(Int64, 7))
    ts_i = src[1];  te_i = src[2];  p = src[3];  p_i = src[4];  os_p = src[5];  W_p = src[7]
    ts_j = dest[1]; te_j = dest[2]; q = dest[3]; q_i = dest[4]; os_q = dest[5]; W_q = dest[7]
    ts_max = max(ts_i, ts_j); ts_max -= get_δt(ts_max, Δ)
    te_min = min(te_i, te_j); te_min -= get_δt(te_min, Δ)
    nov = 1 + div(te_min - ts_max, Δ)
    Xq = getindex(X, q)

    # (1) determine the times and indices of overlap within each pair
    # a. determine sample times of overlap
    Ti = collect(ts_max:Δ:te_min)
    Tj = deepcopy(Ti)

    # b. get sample indices within each overlap window
    # i
    xsi_i = round(Int64, (ts_max - ts_i)/Δ) + os_p
    xei_i = xsi_i + nov - 1
    # j
    xsi_j = round(Int64, (ts_max - ts_j)/Δ) + os_q
    xei_j = xsi_j + nov - 1

    # (2) Extract sample windows
    Xi = getindex(getindex(X, p), xsi_i:xei_i)
    Xj = getindex(getindex(X, q), xsi_j:xei_j)
    lxp = length(getindex(X, p))
    lxq = length(getindex(X, q))

    # ================================================================
    # check for duplicate windows
    if (ts_i == ts_j) && (te_i == te_j) && (Xi == Xj)
      # delete time window
      W = W[setdiff(1:end, W_p), :]
    else
      # Check for misalignment:
      τ, χ, do_xtmerge, δj = check_alignment(Ti, Tj, Xi, Xj, Δ)
      if do_xtmerge
        xtmerge!(τ, χ, div(Δ,2))
      end
      if δj != 0
        xsi_i += δj
        xei_j -= δj
      end

      # (3) Merge X,T into S[q]
      deleteat!(Xq, xsi_j:xei_j)
      if xsi_j == 1
        prepend!(Xq, χ)
      else
        splice!(Xq, xsi_j:xsi_j-1, χ)
      end

      # (4) Adjust start, end indices of windows ≥ q_i in q
      # structure: w_start, w_end, channel_number, window_number, x_start, x_end
      nxq = length(Xq) - lxq
      i = 0
      while i < nW
        i += 1
        if W[i, 3] == q && W[i, 4] ≥ q_i
          W[i, 1] += nxq*Δ
          W[i, 2] += nxq*Δ
          W[i, 5] += nxq
          W[i, 6] += nxq
        end
      end

      #= if xsi_i ≤ os_p (which is always true, at this point in
      the control flow), we decrease W[W_p, 1:2] =#
      nxp = xei_i-xsi_i+1
      W[W_p, 1] -= δj*Δ
      W[W_p, 2] -= (nxp + δj)*Δ

      #= Control for when window P is emptied; the above two statements
      make this possible =#
      if (W[W_p, 2] < W[W_p, 1])
        W = W[setdiff(1:end, W_p), :]
      else
        W[W_p, 6] -= nxp
      end
    end
    # Sort by end time, to ensure we pick the window with latest end next
    k = sortperm(W[:, 2], rev=true)
    W = W[k, :]
    nW = size(W, 1)

    # Repeat until no further merges are possible
    (src, dest) = get_next_pair(W)
  end

  kk = sortperm(W[:, 1])
  W = W[kk, :]
  #= At this point, we have nothing left that can be merged. So we're going
  to arrange T[subgrp] and X[subgrp] in windows using t_win =#

  n = size(W, 1)
  nx = broadcast(+, getindex(W, :, 6).-getindex(W, :,5), 1)
  X_Ω = Array{eltype(X[Ω]),1}(undef, sum(nx))
  xi = 1
  i = 0
  while i < n
    i   = i + 1
    p   = getindex(W, i, 3)
    lx  = getindex(nx, i)
    copyto!(X_Ω, xi, getindex(X, p), getindex(W, i, 5), lx)
    xi  = xi + lx
  end

  # Shrink W and eliminate windows with no actual gap between them
  m = trues(n)
  while n > 1
    if W[n-1, 2] + Δ == W[n,1]
      W[n-1, 2] = W[n,2]
      m[n] = false
    end
    n = n - 1
  end
  W = W[m,[1,2]]
  broadcast!(+, W, W, δts)
  T_Ω = w_time(W, Δ)
  return T_Ω, X_Ω
end

# merges into the channel with the most recent data
function merge_non_ts!(S::GphysData, subgrp::Array{Int64,1})
  te = [maximum(t[:,2]) for t in S.t[subgrp]]
  Ω = subgrp[argmax(te)]
  T = vcat(S.t[subgrp]...)[:,2]
  X = vcat(S.x[subgrp]...)
  Z = unique(collect(zip(T,X)))
  T = first.(Z)
  X = last.(Z)
  ii = sortperm(T)
  S.t[Ω] = hcat(collect(1:1:length(T)), T[ii])
  S.x[Ω] = X[ii]
  return Ω
end

function merge_non_ts!(C::GphysChannel, D::GphysChannel)
  T = vcat(C.t, D.t)[:, 2]
  X = vcat(C.x, D.x)
  Z = unique(collect(zip(T,X)))
  T = first.(Z)
  X = last.(Z)
  ii = sortperm(T)
  C.t = hcat(collect(1:length(T)), T[ii])
  C.x = X[ii]
  return nothing
end


text:
# TODO: this should be more robust, ideally.
using Logging
import WebIO, HTTP, WebSockets

const pages = Dict{String,Any}()
const server = Ref{Any}()
const port = Ref{Int}(9000)

function isrunning(server)
    isassigned(server) && !istaskdone(server[].serve_task)
end

# Be ready for the deprecation of tohtml :P
# once the no unsafe-script lands
if isdefined(WebIO, :tohtml)
    const tohtml = WebIO.tohtml
else
    tohtml(io, app) = show(io, MIME"text/html"(), app)
end

function routepages(req)
    target = req.target[2:end]

    haskey(pages, target) || return missing

    io = IOBuffer()
    webio_script = WebIO.wio_asseturl("/webio/dist/bundle.js")
    ws_script = WebIO.wio_asseturl("/providers/websocket_connection.js")
    print(io, """
        <!doctype html>
        <html>
        <head>
        <meta charset="UTF-8">
        <script>var websocket_url = 'ws://localhost:$(port[])/webio_websocket/'</script>
        <script src=$(repr(webio_script))></script>
        <script src=$(repr(ws_script))></script>
        </head>
        <body>
    """)
    tohtml(io, pages[target])
    print(io, """
        </body>
        </html>
    """)

    return HTTP.Response(
        200,
        ["Content-Type" => "text/html"],
        body = String(take!(io))
    )
end

function Base.show(io::IO, ::MIME"application/juno+plotpane", n::Union{WebIO.Node, WebIO.Scope, WebIO.AbstractWidget})
    global pages, server
    id = rand(UInt128)
    pages[string(id)] = WebIO.render(n)

    if !isrunning(server)
        server[] = setup_server()
    end

    print(io, "<meta http-equiv=\"refresh\" content=\"0; url=http://localhost:$(port[])/$(id)\"/>")
end

function setup_server()
    port[] = rand(8000:9000)

    server = nothing
    # STFU
    with_logger(Logging.NullLogger()) do
        server = WebIO.WebIOServer(
            routepages,
            http_port = port[],
            singleton = false
        )
    end

    return server
end


text:
export ColoringParams, colorize, blend, reorderColor, existing_cmaps, cmap,
       overlay, linearDoge, complexColoring

struct ColoringParams
  cmin::Float64
  cmax::Float64
  cmap::String
end

function ColoringParams(cmin::Real, cmax::Real, i::Int)
  all_maps = existing_cmaps()
  colormap = all_maps[i+1]
  return ColoringParams(cmin, cmax, colormap)
end

#####################
# Weighting function
#####################

"""
    weighted_color_mean(w1, c1, c2)
Returns the color `w1*c1 + (1-w1)*c2` that is the weighted mean of `c1` and
`c2`, where `c1` has a weight 0 ≤ `w1` ≤ 1.
"""
function weighted_color_mean(w1::Real, c1::Colorant, c2::Colorant)
    weight1 = convert(promote_type(eltype(c1), eltype(c2)),w1)
    weight2 = weight1 >= 0 && weight1 <= 1 ? oftype(weight1,1-weight1) : throw(DomainError())
    mapc((x,y)->weight1*x+weight2*y, c1, c2)
end

#####################
# Colorizing an image
#####################

function getBounds( data, coloring, minval, maxval)
  rangeval = maxval - minval
  lowerBound = minval + coloring.cmin * rangeval
  higherBound = minval + coloring.cmax * rangeval
  return lowerBound, higherBound
end

function colorize(data::AbstractArray{T}, coloring::Vector{ColoringParams},
                  minval::Vector, maxval::Vector, blendChannels, activeChannel, complexBlending) where {T<:Real}
  if blendChannels
    if complexBlending
      if size(data,1) != 2
        error("Data needs to be of length 2")
      end
      return complexColoring(sliceColorDim(data,1), sliceColorDim(data,2))
    else
      return colorize(data, coloring, minval, maxval)
    end
  else
    return colorize(data, coloring, minval, maxval, activeChannel)
  end
end

function colorize(data::AbstractArray{T}, coloring::Vector{ColoringParams},
                 minval::Vector, maxval::Vector) where {T<:Real}
  I = [colorize(sliceColorDim(data,i), coloring[i], minval[i], maxval[i]) for i = 1:size(data,1)]
  return linearDodge(I)
end

function colorize(data::AbstractArray{T}, coloring::Vector{ColoringParams},
                 minval, maxval, chan::Int) where {T<:Real}
  return colorize(sliceColorDim(data,chan), coloring[chan], minval[chan], maxval[chan])
end

function colorize(data::AbstractArray{T}, c::ColoringParams,
                  minval::Real, maxval::Real) where {T<:Real}

  lowerBound, higherBound = getBounds(data, c, minval, maxval)

  cdata = colorize( data, lowerBound, higherBound, cmap(c.cmap),
                    normalize=false)

  return cdata
end

"""
  coloredimage = colorize(inputimage,wmin,wmax,cmap)

This function uses the windowing parameters `w1` and `w2` together with the
color map `cmap` to colorize the `inputimage`. Using the keywordargument
`normalize` the data will be scaled such that all data lies in [0,1].
"""
function colorize(inputimage::AbstractArray, wmin::T, wmax::T, cmap;
                  normalize=true) where {T<:Number}
  if normalize
    image = copy(inputimage)
    image ./= maximum(image)
  else
    image = inputimage
  end

  return _colorize(image, wmin, wmax, cmap)
end

function _colorize(image::AbstractArray,wmin,wmax,cmap::Vector{T}) where {T<:Colorant}
  C = similar(image, T)
  _colorize!(image,C,wmin,wmax,cmap)
  return C
end

function _colorize!(image,C,wmin,wmax,cmap)
  for l=1:length(image)
    C[l] = _colorize(image[l],wmin,wmax,cmap)
  end
  C
end

function _colorize(x::Number,wmin,wmax,cmap)
  L = length(cmap)
  y = (x-wmin) / (wmax-wmin)*(L-1) + 1

  if y <= 1 || wmax == wmin
    c = first(cmap)
  elseif y >= L
    c = last(cmap)
  else
    yidx = floor(Int, y)
    c = weighted_color_mean(y-yidx,cmap[yidx+1],cmap[yidx])
  end
  c
end

##################
#  Linear dodge  #
##################
"""
	imageC = linearDodge(imageA,imageB)

Blends imageA and imageB using the linear dodge mode, which sums pixelwise with
subsequent clamping.

This function may be used for images of the same size: size(imageA)==size(imageB)
or for a static imageA and a time series of static images imageB.
"""
function linearDodge(images::Vector)
	out = zero(images[1])
	for imag in images
		out .+= imag
	end
	return clamp01!(out)
end

linearDodge(images...) = linearDodge([images...])

complexColoring(C::Array{T,2}; colormap=ColorSchemes.phase, normalizeG=true, g=0.7, amax=maximum(abs.(C))) where T<:Complex = complexColoring(abs.(C),angle.(C),colormap=colormap, normalizeG=normalizeG, g=g, amax=amax)

function complexColoring(amp, phase; colormap=ColorSchemes.phase, normalizeG::Bool=true, g=0.7, amax=maximum(amp))
    if normalizeG
        colormap = ColorScheme(normalizeGray.(colormap,g))
    end
    # normalize amplitude
    I = amp./amax
    clamp01!(I) # using desaturation colormap: make sure to map intensity between 0 and 1
    rawImage = similar(I,RGBA{N0f8})
    for n=1:length(I)
      rawImage[n] = I[n]*get(colormap,phase[n],(-π,Float64(π)))
    end
    return convert.(RGBA{N0f8},rawImage)
end

function overlay(imageBG::AbstractArray{T,D}, imageFG::AbstractArray{U,D},
                  translucentColor) where {T,U,D}
  C = similar(imageFG, RGBA{N0f8})
  overlay!(imageBG, imageFG, C, 0, translucentColor)
  return C
end

# This version operates on a 3D bg and a 4D fg
function overlay(imageBG::AbstractArray{T,3}, imageFG::AbstractArray{U,4},
                  translucentColor) where {T,U}
  C = similar(imageFG, RGBA{N0f8})
  offset = 0
  for l=1:size(C,4)
    overlay!(imageBG, imageFG, C, offset, translucentColor)
    offset += length(imageBG)
  end
  return C
end

function overlay!(imageBG, imageFG, C, offset, translucentColor)
  for l=1:length(imageBG)
    IFG = imageFG[l+offset]
    IBG = imageBG[l]
    C[l+offset] = ( red(IFG) == red(translucentColor) &&
                    green(IFG) == green(translucentColor) &&
		                blue(IFG) == blue(translucentColor)) ? IBG : IFG
  end
  return C
end

##################
# Alpha blending
##################

# This version works on the RGB level and alpha
# is associated to the foreground color.
function blend(bg::C, fg::C, alpha) where {C<:Color}
  #alpha = oftype(eltype(C),alpha)
  output = alpha*fg + (1-alpha)*bg
end

# Type stable?
function blend(bg::C, fg::C) where {C<:TransparentColor}
  α = alpha(bg)
  β  = alpha(fg) * (one(α)-α)
  return mapc((x,y)->α*x + β*y, bg, fg)
end

# this should be replaced by weighted_color_mean(w1, c1, c2) as soon as
# weighted_color_mean supports all colors since that is what it is
function blend(bg::Color, fg::TransparentColor)
  output = blend(bg, color(fg), alpha(fg))
end

# This version operates on two equally dimensioned images
function blend(imageBG::AbstractArray{T,D}, imageFG::AbstractArray{U,D}) where {T,U,D}
  C = similar(imageFG, RGBA{N0f8})
  blend!(imageBG, imageFG, C)
  return C
end

# This version operates on a 3D bg and a 4D fg
function blend(imageBG::AbstractArray{T,3}, imageFG::AbstractArray{U,4}) where {T,U}
  C = similar(imageFG, RGBA{N0f8})
  offset = 0
  for l=1:size(C,4)
    blend!(imageBG, imageFG, C, offset)
    offset += length(imageBG)
  end
  return C
end

function blend!(imageBG, imageFG, C, offset=0)
  for l=1:length(imageBG)
    C[l+offset] = blend(imageBG[l],imageFG[l+offset])
  end
  return C
end

##################
# reorder color
##################

"""
	reorderColor(image)

separates the color channels of `imgage`.
"""
function reorderColor(image::AbstractArray{T}) where {T<:Colorant}
  imageOut = similar(image, UInt8, size(image)..., 3)
  for l=1:length(image)
    imageOut[l] = reinterpret(UInt8,red(image[l]))
    imageOut[l+length(image)] = reinterpret(UInt8,green(image[l]))
    imageOut[l+2*length(image)] = reinterpret(UInt8,blue(image[l]))
  end
  imageOut
end

function reorderColor(image::AbstractArray{UInt8})
  imageOut = similar(image, RGB{N0f8}, size(image)[1:end-1]...)
  for l=1:length(imageOut)
    imageOut[l] = RGB{N0f8}(reinterpret(N0f8,image[l]),reinterpret(Ufixed8,image[l+length(imageOut)]),reinterpret(N0f8,image[l+2*length(imageOut)]))
  end
  imageOut
end

##################
# Color maps
##################

"""
    `normalizeGray(colormap,g=0.45)`

Transform a color lightness value in Lab space with the aim to output a color
of similar color but with the specified gray value.
"""
function normalizeGray(c::T,g=0.7) where {T<:Colorant}
    cluv = Lab(c)
    # test range
    minG = Gray.(Lab(0,cluv.a,cluv.b)).val
    maxG = Gray.(Lab(100,cluv.a,cluv.b)).val
    if g < minG
      g = minG
      @warn "Normalization of gray value failed, g ≈ $(ceil(minG,digits=3)) is used. Use g >= $(ceil(minG,digits=3)) for a consistent normalization."
    elseif g > maxG
      g = maxG
      @warn "Normalization of gray value failed, g ≈ $(floor(maxG,digits=3)) is used. Use g <= $(floor(maxG,digits=3)) for a consistent normalization."
    end

    l = find_zero(l->Gray.(Lab(l,cluv.a,cluv.b)).val-g, (0, 100))

    return T(Lab(l,cluv.a,cluv.b))
end

"""
    cmap(color)

Creates a colormap from a single `color`, with a color gradient from
transparent black to the opaque `color`.
"""
function cmap(color::T) where {T<:Colorant}
	c1 = RGBA{N0f8}(0,0,0,0)
	c2 = convert(RGBA{N0f8},color)
	return RGBA{N0f8}[c1,c2]
end

"""
    cmap(colors)

Maps a vector of tuples (length 3 or 4) to `RGBA{N0f8}`.
"""
cmap(colors::Vector{T}) where {T<:Tuple} = map(x->RGBA{N0f8}(x...),colors)

# New matplotlib colormaps by Nathaniel J. Smith, Stefan van der Walt,
# and (in the case of viridis) Eric Firing.
#
# This file and the colormaps in it are released under the CC0 license /
# public domain dedication. We would appreciate credit if you use or
# redistribute these colormaps, but do not impose any legal restrictions.
#
# To the extent possible under law, the persons who associated CC0 with
# mpl-colormaps have waived all copyright and related or neighboring rights
# to mpl-colormaps.
#
# You should have received a copy of the CC0 legalcode along with this
# work.  If not, see <http://creativecommons.org/publicdomain/zero/1.0/>.
# see https://github.com/BIDS/colormap/blob/master/colormaps.py
function cmap(str::AbstractString)::Vector{RGBA{Normed{UInt8,8}}}
  if str == "gray"
    return cmap(colorant"white")
  elseif str == "blue"
    return cmap(colorant"blue")
  elseif str == "green"
    return cmap(colorant"green")
  elseif str == "red"
    return cmap(colorant"red")
  elseif str == "UKE"
    return cmap(RGBA(0,73/255,146/255))
  elseif str == "redgreen"
    return cmap([(1.0,0,0,1),(0,1.0,0,1)])
  elseif str == "redgreenalpha"
    return cmap([(1.0,0,0,0), (1.0,0,0,1), (0.0,0.5,0,1), (0.0,1,0,1)])
  elseif str == "greenorangealpha"
    return cmap([(0,1.0,0,0.0),(1.0,0.498,0,1)])
  elseif str == "greenblueredalpha"
    return cmap([(0,1.0,0,0),(0,0,1.0,0.5),(1.0,0,0,1)])
  elseif str == "bluegreenredalpha"
    return cmap([(0,0.2,1.0,0),(0,1.0,0,0.5),(1.0,0.498,0,1)])
  elseif str == "yellowredalpha"
    return cmap([(1.0,1.0,0,0),(1.0,0,0,1)])
  elseif str == "tobi"
    return cmap([(1.0,1.0,0,0), (1.0,0.498,0,1),(1.0,0.0,0,1),(0.0,1,0,1), (0.0,0,1,1)])
  elseif str == "tobitr"
    return cmap([(1.0,1.0,0,0), (1.0,0.498,0,0.6),(1.0,0.0,0,0.6),(0.0,1,0,0.6), (0.0,0,1,0.6)])
  elseif str == "redwhite"
    return cmap([(1.0,0,0,0),(1.0,1.0,1.0,1.0)])
  elseif str == "timemap"
    return cmap([(0,0,1.0,1.0),(0,1.0,1.0,1.0),(0,1.0,0,1.0),(1.0,1.0,0,1.0),(1.0,0.5,0,1.0),(1.0,0,0,1.0),(139/255,35/255,35/255,1.0)])
  elseif str == "UPthreshold"
    return cmap([(0,0,0,0),(0,0,1.0,1.0),(0,1.0,0,1.0),(1.0,1.0,0,1.0),(1.0,0,0,1.0),(139/255,35/255,35/255,1.0),(0,0,0,0)])
  elseif str == "perfusionmap"
    return cmap([(0,0,0,0),(0,0,1.0,1.0),(0,1.0,0,1.0),(1.0,1.0,0,1.0),(1.0,0,0,1.0),(139/255,35/255,35/255,1.0)])
  elseif str == "Liver"
    return cmap([(0,0.2,1.0,0.6), (0.0,1.0,0.0,0.6),(1.0,1.0,0.0,0.6),(1.0,0.0,0.0,0.6)])
  elseif str == "viridis"
    return cmap([(0.267004,0.004874,0.329415,0.0),(0.26851,0.009605,0.335427,0.003922),(0.269944,0.014625,0.341379,0.007843),(0.271305,0.019942,0.347269,0.011765),(0.272594,0.025563,0.353093,0.015686),(0.273809,0.031497,0.358853,0.019608),(0.274952,0.037752,0.364543,0.023529),(0.276022,0.044167,0.370164,0.027451),(0.277018,0.050344,0.375715,0.031373),(0.277941,0.056324,0.381191,0.035294),(0.278791,0.062145,0.386592,0.039216),(0.279566,0.067836,0.391917,0.043137),(0.280267,0.073417,0.397163,0.047059),(0.280894,0.078907,0.402329,0.05098),(0.281446,0.08432,0.407414,0.054902),(0.281924,0.089666,0.412415,0.058824),(0.282327,0.094955,0.417331,0.062745),(0.282656,0.100196,0.42216,0.066667),(0.28291,0.105393,0.426902,0.070588),(0.283091,0.110553,0.431554,0.07451),(0.283197,0.11568,0.436115,0.078431),(0.283229,0.120777,0.440584,0.082353),(0.283187,0.125848,0.44496,0.086275),(0.283072,0.130895,0.449241,0.090196),(0.282884,0.13592,0.453427,0.094118),(0.282623,0.140926,0.457517,0.098039),(0.28229,0.145912,0.46151,0.101961),(0.281887,0.150881,0.465405,0.105882),(0.281412,0.155834,0.469201,0.109804),(0.280868,0.160771,0.472899,0.113725),(0.280255,0.165693,0.476498,0.117647),(0.279574,0.170599,0.479997,0.121569),(0.278826,0.17549,0.483397,0.12549),(0.278012,0.180367,0.486697,0.129412),(0.277134,0.185228,0.489898,0.133333),(0.276194,0.190074,0.493001,0.137255),(0.275191,0.194905,0.496005,0.141176),(0.274128,0.199721,0.498911,0.145098),(0.273006,0.20452,0.501721,0.14902),(0.271828,0.209303,0.504434,0.152941),(0.270595,0.214069,0.507052,0.156863),(0.269308,0.218818,0.509577,0.160784),(0.267968,0.223549,0.512008,0.164706),(0.26658,0.228262,0.514349,0.168627),(0.265145,0.232956,0.516599,0.172549),(0.263663,0.237631,0.518762,0.176471),(0.262138,0.242286,0.520837,0.180392),(0.260571,0.246922,0.522828,0.184314),(0.258965,0.251537,0.524736,0.188235),(0.257322,0.25613,0.526563,0.192157),(0.255645,0.260703,0.528312,0.196078),(0.253935,0.265254,0.529983,0.2),(0.252194,0.269783,0.531579,0.203922),(0.250425,0.27429,0.533103,0.207843),(0.248629,0.278775,0.534556,0.211765),(0.246811,0.283237,0.535941,0.215686),(0.244972,0.287675,0.53726,0.219608),(0.243113,0.292092,0.538516,0.223529),(0.241237,0.296485,0.539709,0.227451),(0.239346,0.300855,0.540844,0.231373),(0.237441,0.305202,0.541921,0.235294),(0.235526,0.309527,0.542944,0.239216),(0.233603,0.313828,0.543914,0.243137),(0.231674,0.318106,0.544834,0.247059),(0.229739,0.322361,0.545706,0.25098),(0.227802,0.326594,0.546532,0.254902),(0.225863,0.330805,0.547314,0.258824),(0.223925,0.334994,0.548053,0.262745),(0.221989,0.339161,0.548752,0.266667),(0.220057,0.343307,0.549413,0.270588),(0.21813,0.347432,0.550038,0.27451),(0.21621,0.351535,0.550627,0.278431),(0.214298,0.355619,0.551184,0.282353),(0.212395,0.359683,0.55171,0.286275),(0.210503,0.363727,0.552206,0.290196),(0.208623,0.367752,0.552675,0.294118),(0.206756,0.371758,0.553117,0.298039),(0.204903,0.375746,0.553533,0.301961),(0.203063,0.379716,0.553925,0.305882),(0.201239,0.38367,0.554294,0.309804),(0.19943,0.387607,0.554642,0.313725),(0.197636,0.391528,0.554969,0.317647),(0.19586,0.395433,0.555276,0.321569),(0.1941,0.399323,0.555565,0.32549),(0.192357,0.403199,0.555836,0.329412),(0.190631,0.407061,0.556089,0.333333),(0.188923,0.41091,0.556326,0.337255),(0.187231,0.414746,0.556547,0.341176),(0.185556,0.41857,0.556753,0.345098),(0.183898,0.422383,0.556944,0.34902),(0.182256,0.426184,0.55712,0.352941),(0.180629,0.429975,0.557282,0.356863),(0.179019,0.433756,0.55743,0.360784),(0.177423,0.437527,0.557565,0.364706),(0.175841,0.44129,0.557685,0.368627),(0.174274,0.445044,0.557792,0.372549),(0.172719,0.448791,0.557885,0.376471),(0.171176,0.45253,0.557965,0.380392),(0.169646,0.456262,0.55803,0.384314),(0.168126,0.459988,0.558082,0.388235),(0.166617,0.463708,0.558119,0.392157),(0.165117,0.467423,0.558141,0.396078),(0.163625,0.471133,0.558148,0.4),(0.162142,0.474838,0.55814,0.403922),(0.160665,0.47854,0.558115,0.407843),(0.159194,0.482237,0.558073,0.411765),(0.157729,0.485932,0.558013,0.415686),(0.15627,0.489624,0.557936,0.419608),(0.154815,0.493313,0.55784,0.423529),(0.153364,0.497,0.557724,0.427451),(0.151918,0.500685,0.557587,0.431373),(0.150476,0.504369,0.55743,0.435294),(0.149039,0.508051,0.55725,0.439216),(0.147607,0.511733,0.557049,0.443137),(0.14618,0.515413,0.556823,0.447059),(0.144759,0.519093,0.556572,0.45098),(0.143343,0.522773,0.556295,0.454902),(0.141935,0.526453,0.555991,0.458824),(0.140536,0.530132,0.555659,0.462745),(0.139147,0.533812,0.555298,0.466667),(0.13777,0.537492,0.554906,0.470588),(0.136408,0.541173,0.554483,0.47451),(0.135066,0.544853,0.554029,0.478431),(0.133743,0.548535,0.553541,0.482353),(0.132444,0.552216,0.553018,0.486275),(0.131172,0.555899,0.552459,0.490196),(0.129933,0.559582,0.551864,0.494118),(0.128729,0.563265,0.551229,0.498039),(0.127568,0.566949,0.550556,0.501961),(0.126453,0.570633,0.549841,0.505882),(0.125394,0.574318,0.549086,0.509804),(0.124395,0.578002,0.548287,0.513725),(0.123463,0.581687,0.547445,0.517647),(0.122606,0.585371,0.546557,0.521569),(0.121831,0.589055,0.545623,0.52549),(0.121148,0.592739,0.544641,0.529412),(0.120565,0.596422,0.543611,0.533333),(0.120092,0.600104,0.54253,0.537255),(0.119738,0.603785,0.5414,0.541176),(0.119512,0.607464,0.540218,0.545098),(0.119423,0.611141,0.538982,0.54902),(0.119483,0.614817,0.537692,0.552941),(0.119699,0.61849,0.536347,0.556863),(0.120081,0.622161,0.534946,0.560784),(0.120638,0.625828,0.533488,0.564706),(0.12138,0.629492,0.531973,0.568627),(0.122312,0.633153,0.530398,0.572549),(0.123444,0.636809,0.528763,0.576471),(0.12478,0.640461,0.527068,0.580392),(0.126326,0.644107,0.525311,0.584314),(0.128087,0.647749,0.523491,0.588235),(0.130067,0.651384,0.521608,0.592157),(0.132268,0.655014,0.519661,0.596078),(0.134692,0.658636,0.517649,0.6),(0.137339,0.662252,0.515571,0.603922),(0.14021,0.665859,0.513427,0.607843),(0.143303,0.669459,0.511215,0.611765),(0.146616,0.67305,0.508936,0.615686),(0.150148,0.676631,0.506589,0.619608),(0.153894,0.680203,0.504172,0.623529),(0.157851,0.683765,0.501686,0.627451),(0.162016,0.687316,0.499129,0.631373),(0.166383,0.690856,0.496502,0.635294),(0.170948,0.694384,0.493803,0.639216),(0.175707,0.6979,0.491033,0.643137),(0.180653,0.701402,0.488189,0.647059),(0.185783,0.704891,0.485273,0.65098),(0.19109,0.708366,0.482284,0.654902),(0.196571,0.711827,0.479221,0.658824),(0.202219,0.715272,0.476084,0.662745),(0.20803,0.718701,0.472873,0.666667),(0.214,0.722114,0.469588,0.670588),(0.220124,0.725509,0.466226,0.67451),(0.226397,0.728888,0.462789,0.678431),(0.232815,0.732247,0.459277,0.682353),(0.239374,0.735588,0.455688,0.686275),(0.24607,0.73891,0.452024,0.690196),(0.252899,0.742211,0.448284,0.694118),(0.259857,0.745492,0.444467,0.698039),(0.266941,0.748751,0.440573,0.701961),(0.274149,0.751988,0.436601,0.705882),(0.281477,0.755203,0.432552,0.709804),(0.288921,0.758394,0.428426,0.713725),(0.296479,0.761561,0.424223,0.717647),(0.304148,0.764704,0.419943,0.721569),(0.311925,0.767822,0.415586,0.72549),(0.319809,0.770914,0.411152,0.729412),(0.327796,0.77398,0.40664,0.733333),(0.335885,0.777018,0.402049,0.737255),(0.344074,0.780029,0.397381,0.741176),(0.35236,0.783011,0.392636,0.745098),(0.360741,0.785964,0.387814,0.74902),(0.369214,0.788888,0.382914,0.752941),(0.377779,0.791781,0.377939,0.756863),(0.386433,0.794644,0.372886,0.760784),(0.395174,0.797475,0.367757,0.764706),(0.404001,0.800275,0.362552,0.768627),(0.412913,0.803041,0.357269,0.772549),(0.421908,0.805774,0.35191,0.776471),(0.430983,0.808473,0.346476,0.780392),(0.440137,0.811138,0.340967,0.784314),(0.449368,0.813768,0.335384,0.788235),(0.458674,0.816363,0.329727,0.792157),(0.468053,0.818921,0.323998,0.796078),(0.477504,0.821444,0.318195,0.8),(0.487026,0.823929,0.312321,0.803922),(0.496615,0.826376,0.306377,0.807843),(0.506271,0.828786,0.300362,0.811765),(0.515992,0.831158,0.294279,0.815686),(0.525776,0.833491,0.288127,0.819608),(0.535621,0.835785,0.281908,0.823529),(0.545524,0.838039,0.275626,0.827451),(0.555484,0.840254,0.269281,0.831373),(0.565498,0.84243,0.262877,0.835294),(0.575563,0.844566,0.256415,0.839216),(0.585678,0.846661,0.249897,0.843137),(0.595839,0.848717,0.243329,0.847059),(0.606045,0.850733,0.236712,0.85098),(0.616293,0.852709,0.230052,0.854902),(0.626579,0.854645,0.223353,0.858824),(0.636902,0.856542,0.21662,0.862745),(0.647257,0.8584,0.209861,0.866667),(0.657642,0.860219,0.203082,0.870588),(0.668054,0.861999,0.196293,0.87451),(0.678489,0.863742,0.189503,0.878431),(0.688944,0.865448,0.182725,0.882353),(0.699415,0.867117,0.175971,0.886275),(0.709898,0.868751,0.169257,0.890196),(0.720391,0.87035,0.162603,0.894118),(0.730889,0.871916,0.156029,0.898039),(0.741388,0.873449,0.149561,0.901961),(0.751884,0.874951,0.143228,0.905882),(0.762373,0.876424,0.137064,0.909804),(0.772852,0.877868,0.131109,0.913725),(0.783315,0.879285,0.125405,0.917647),(0.79376,0.880678,0.120005,0.921569),(0.804182,0.882046,0.114965,0.92549),(0.814576,0.883393,0.110347,0.929412),(0.82494,0.88472,0.106217,0.933333),(0.83527,0.886029,0.102646,0.937255),(0.845561,0.887322,0.099702,0.941176),(0.85581,0.888601,0.097452,0.945098),(0.866013,0.889868,0.095953,0.94902),(0.876168,0.891125,0.09525,0.952941),(0.886271,0.892374,0.095374,0.956863),(0.89632,0.893616,0.096335,0.960784),(0.906311,0.894855,0.098125,0.964706),(0.916242,0.896091,0.100717,0.968627),(0.926106,0.89733,0.104071,0.972549),(0.935904,0.89857,0.108131,0.976471),(0.945636,0.899815,0.112838,0.980392),(0.9553,0.901065,0.118128,0.984314),(0.964894,0.902323,0.123941,0.988235),(0.974417,0.90359,0.130215,0.992157),(0.983868,0.904867,0.136897,0.996078),(0.993248,0.906157,0.143936,1.0)])
  elseif str == "plasma"
    return cmap([(0.050383,0.029803,0.527975,0.0),(0.063536,0.028426,0.533124,0.003922),(0.075353,0.027206,0.538007,0.007843),(0.086222,0.026125,0.542658,0.011765),(0.096379,0.025165,0.547103,0.015686),(0.10598,0.024309,0.551368,0.019608),(0.115124,0.023556,0.555468,0.023529),(0.123903,0.022878,0.559423,0.027451),(0.132381,0.022258,0.56325,0.031373),(0.140603,0.021687,0.566959,0.035294),(0.148607,0.021154,0.570562,0.039216),(0.156421,0.020651,0.574065,0.043137),(0.16407,0.020171,0.577478,0.047059),(0.171574,0.019706,0.580806,0.05098),(0.17895,0.019252,0.584054,0.054902),(0.186213,0.018803,0.587228,0.058824),(0.193374,0.018354,0.59033,0.062745),(0.200445,0.017902,0.593364,0.066667),(0.207435,0.017442,0.596333,0.070588),(0.21435,0.016973,0.599239,0.07451),(0.221197,0.016497,0.602083,0.078431),(0.227983,0.016007,0.604867,0.082353),(0.234715,0.015502,0.607592,0.086275),(0.241396,0.014979,0.610259,0.090196),(0.248032,0.014439,0.612868,0.094118),(0.254627,0.013882,0.615419,0.098039),(0.261183,0.013308,0.617911,0.101961),(0.267703,0.012716,0.620346,0.105882),(0.274191,0.012109,0.622722,0.109804),(0.280648,0.011488,0.625038,0.113725),(0.287076,0.010855,0.627295,0.117647),(0.293478,0.010213,0.62949,0.121569),(0.299855,0.009561,0.631624,0.12549),(0.30621,0.008902,0.633694,0.129412),(0.312543,0.008239,0.6357,0.133333),(0.318856,0.007576,0.63764,0.137255),(0.32515,0.006915,0.639512,0.141176),(0.331426,0.006261,0.641316,0.145098),(0.337683,0.005618,0.643049,0.14902),(0.343925,0.004991,0.64471,0.152941),(0.35015,0.004382,0.646298,0.156863),(0.356359,0.003798,0.64781,0.160784),(0.362553,0.003243,0.649245,0.164706),(0.368733,0.002724,0.650601,0.168627),(0.374897,0.002245,0.651876,0.172549),(0.381047,0.001814,0.653068,0.176471),(0.387183,0.001434,0.654177,0.180392),(0.393304,0.001114,0.655199,0.184314),(0.399411,0.000859,0.656133,0.188235),(0.405503,0.000678,0.656977,0.192157),(0.41158,0.000577,0.65773,0.196078),(0.417642,0.000564,0.65839,0.2),(0.423689,0.000646,0.658956,0.203922),(0.429719,0.000831,0.659425,0.207843),(0.435734,0.001127,0.659797,0.211765),(0.441732,0.00154,0.660069,0.215686),(0.447714,0.00208,0.66024,0.219608),(0.453677,0.002755,0.66031,0.223529),(0.459623,0.003574,0.660277,0.227451),(0.46555,0.004545,0.660139,0.231373),(0.471457,0.005678,0.659897,0.235294),(0.477344,0.00698,0.659549,0.239216),(0.48321,0.00846,0.659095,0.243137),(0.489055,0.010127,0.658534,0.247059),(0.494877,0.01199,0.657865,0.25098),(0.500678,0.014055,0.657088,0.254902),(0.506454,0.016333,0.656202,0.258824),(0.512206,0.018833,0.655209,0.262745),(0.517933,0.021563,0.654109,0.266667),(0.523633,0.024532,0.652901,0.270588),(0.529306,0.027747,0.651586,0.27451),(0.534952,0.031217,0.650165,0.278431),(0.54057,0.03495,0.64864,0.282353),(0.546157,0.038954,0.64701,0.286275),(0.551715,0.043136,0.645277,0.290196),(0.557243,0.047331,0.643443,0.294118),(0.562738,0.051545,0.641509,0.298039),(0.568201,0.055778,0.639477,0.301961),(0.573632,0.060028,0.637349,0.305882),(0.579029,0.064296,0.635126,0.309804),(0.584391,0.068579,0.632812,0.313725),(0.589719,0.072878,0.630408,0.317647),(0.595011,0.07719,0.627917,0.321569),(0.600266,0.081516,0.625342,0.32549),(0.605485,0.085854,0.622686,0.329412),(0.610667,0.090204,0.619951,0.333333),(0.615812,0.094564,0.61714,0.337255),(0.620919,0.098934,0.614257,0.341176),(0.625987,0.103312,0.611305,0.345098),(0.631017,0.107699,0.608287,0.34902),(0.636008,0.112092,0.605205,0.352941),(0.640959,0.116492,0.602065,0.356863),(0.645872,0.120898,0.598867,0.360784),(0.650746,0.125309,0.595617,0.364706),(0.65558,0.129725,0.592317,0.368627),(0.660374,0.134144,0.588971,0.372549),(0.665129,0.138566,0.585582,0.376471),(0.669845,0.142992,0.582154,0.380392),(0.674522,0.147419,0.578688,0.384314),(0.67916,0.151848,0.575189,0.388235),(0.683758,0.156278,0.57166,0.392157),(0.688318,0.160709,0.568103,0.396078),(0.69284,0.165141,0.564522,0.4),(0.697324,0.169573,0.560919,0.403922),(0.701769,0.174005,0.557296,0.407843),(0.706178,0.178437,0.553657,0.411765),(0.710549,0.182868,0.550004,0.415686),(0.714883,0.187299,0.546338,0.419608),(0.719181,0.191729,0.542663,0.423529),(0.723444,0.196158,0.538981,0.427451),(0.72767,0.200586,0.535293,0.431373),(0.731862,0.205013,0.531601,0.435294),(0.736019,0.209439,0.527908,0.439216),(0.740143,0.213864,0.524216,0.443137),(0.744232,0.218288,0.520524,0.447059),(0.748289,0.222711,0.516834,0.45098),(0.752312,0.227133,0.513149,0.454902),(0.756304,0.231555,0.509468,0.458824),(0.760264,0.235976,0.505794,0.462745),(0.764193,0.240396,0.502126,0.466667),(0.76809,0.244817,0.498465,0.470588),(0.771958,0.249237,0.494813,0.47451),(0.775796,0.253658,0.491171,0.478431),(0.779604,0.258078,0.487539,0.482353),(0.783383,0.2625,0.483918,0.486275),(0.787133,0.266922,0.480307,0.490196),(0.790855,0.271345,0.476706,0.494118),(0.794549,0.27577,0.473117,0.498039),(0.798216,0.280197,0.469538,0.501961),(0.801855,0.284626,0.465971,0.505882),(0.805467,0.289057,0.462415,0.509804),(0.809052,0.293491,0.45887,0.513725),(0.812612,0.297928,0.455338,0.517647),(0.816144,0.302368,0.451816,0.521569),(0.819651,0.306812,0.448306,0.52549),(0.823132,0.311261,0.444806,0.529412),(0.826588,0.315714,0.441316,0.533333),(0.830018,0.320172,0.437836,0.537255),(0.833422,0.324635,0.434366,0.541176),(0.836801,0.329105,0.430905,0.545098),(0.840155,0.33358,0.427455,0.54902),(0.843484,0.338062,0.424013,0.552941),(0.846788,0.342551,0.420579,0.556863),(0.850066,0.347048,0.417153,0.560784),(0.853319,0.351553,0.413734,0.564706),(0.856547,0.356066,0.410322,0.568627),(0.85975,0.360588,0.406917,0.572549),(0.862927,0.365119,0.403519,0.576471),(0.866078,0.36966,0.400126,0.580392),(0.869203,0.374212,0.396738,0.584314),(0.872303,0.378774,0.393355,0.588235),(0.875376,0.383347,0.389976,0.592157),(0.878423,0.387932,0.3866,0.596078),(0.881443,0.392529,0.383229,0.6),(0.884436,0.397139,0.37986,0.603922),(0.887402,0.401762,0.376494,0.607843),(0.89034,0.406398,0.37313,0.611765),(0.89325,0.411048,0.369768,0.615686),(0.896131,0.415712,0.366407,0.619608),(0.898984,0.420392,0.363047,0.623529),(0.901807,0.425087,0.359688,0.627451),(0.904601,0.429797,0.356329,0.631373),(0.907365,0.434524,0.35297,0.635294),(0.910098,0.439268,0.34961,0.639216),(0.9128,0.444029,0.346251,0.643137),(0.915471,0.448807,0.34289,0.647059),(0.918109,0.453603,0.339529,0.65098),(0.920714,0.458417,0.336166,0.654902),(0.923287,0.463251,0.332801,0.658824),(0.925825,0.468103,0.329435,0.662745),(0.928329,0.472975,0.326067,0.666667),(0.930798,0.477867,0.322697,0.670588),(0.933232,0.48278,0.319325,0.67451),(0.93563,0.487712,0.315952,0.678431),(0.93799,0.492667,0.312575,0.682353),(0.940313,0.497642,0.309197,0.686275),(0.942598,0.502639,0.305816,0.690196),(0.944844,0.507658,0.302433,0.694118),(0.947051,0.512699,0.299049,0.698039),(0.949217,0.517763,0.295662,0.701961),(0.951344,0.52285,0.292275,0.705882),(0.953428,0.52796,0.288883,0.709804),(0.95547,0.533093,0.28549,0.713725),(0.957469,0.53825,0.282096,0.717647),(0.959424,0.543431,0.278701,0.721569),(0.961336,0.548636,0.275305,0.72549),(0.963203,0.553865,0.271909,0.729412),(0.965024,0.559118,0.268513,0.733333),(0.966798,0.564396,0.265118,0.737255),(0.968526,0.5697,0.261721,0.741176),(0.970205,0.575028,0.258325,0.745098),(0.971835,0.580382,0.254931,0.74902),(0.973416,0.585761,0.25154,0.752941),(0.974947,0.591165,0.248151,0.756863),(0.976428,0.596595,0.244767,0.760784),(0.977856,0.602051,0.241387,0.764706),(0.979233,0.607532,0.238013,0.768627),(0.980556,0.613039,0.234646,0.772549),(0.981826,0.618572,0.231287,0.776471),(0.983041,0.624131,0.227937,0.780392),(0.984199,0.629718,0.224595,0.784314),(0.985301,0.63533,0.221265,0.788235),(0.986345,0.640969,0.217948,0.792157),(0.987332,0.646633,0.214648,0.796078),(0.98826,0.652325,0.211364,0.8),(0.989128,0.658043,0.2081,0.803922),(0.989935,0.663787,0.204859,0.807843),(0.990681,0.669558,0.201642,0.811765),(0.991365,0.675355,0.198453,0.815686),(0.991985,0.681179,0.195295,0.819608),(0.992541,0.68703,0.19217,0.823529),(0.993032,0.692907,0.189084,0.827451),(0.993456,0.69881,0.186041,0.831373),(0.993814,0.704741,0.183043,0.835294),(0.994103,0.710698,0.180097,0.839216),(0.994324,0.716681,0.177208,0.843137),(0.994474,0.722691,0.174381,0.847059),(0.994553,0.728728,0.171622,0.85098),(0.994561,0.734791,0.168938,0.854902),(0.994495,0.74088,0.166335,0.858824),(0.994355,0.746995,0.163821,0.862745),(0.994141,0.753137,0.161404,0.866667),(0.993851,0.759304,0.159092,0.870588),(0.993482,0.765499,0.156891,0.87451),(0.993033,0.77172,0.154808,0.878431),(0.992505,0.777967,0.152855,0.882353),(0.991897,0.784239,0.151042,0.886275),(0.991209,0.790537,0.149377,0.890196),(0.990439,0.796859,0.14787,0.894118),(0.989587,0.803205,0.146529,0.898039),(0.988648,0.809579,0.145357,0.901961),(0.987621,0.815978,0.144363,0.905882),(0.986509,0.822401,0.143557,0.909804),(0.985314,0.828846,0.142945,0.913725),(0.984031,0.835315,0.142528,0.917647),(0.982653,0.841812,0.142303,0.921569),(0.98119,0.848329,0.142279,0.92549),(0.979644,0.854866,0.142453,0.929412),(0.977995,0.861432,0.142808,0.933333),(0.976265,0.868016,0.143351,0.937255),(0.974443,0.874622,0.144061,0.941176),(0.97253,0.88125,0.144923,0.945098),(0.970533,0.887896,0.145919,0.94902),(0.968443,0.894564,0.147014,0.952941),(0.966271,0.901249,0.14818,0.956863),(0.964021,0.90795,0.14937,0.960784),(0.961681,0.914672,0.15052,0.964706),(0.959276,0.921407,0.151566,0.968627),(0.956808,0.928152,0.152409,0.972549),(0.954287,0.934908,0.152921,0.976471),(0.951726,0.941671,0.152925,0.980392),(0.949151,0.948435,0.152178,0.984314),(0.946602,0.95519,0.150328,0.988235),(0.944152,0.961916,0.146861,0.992157),(0.941896,0.96859,0.140956,0.996078),(0.940015,0.975158,0.131326,1.0)])
  elseif str == "inferno"
    return cmap([(0.001462,0.000466,0.013866,0.0),(0.002267,0.00127,0.01857,0.003922),(0.003299,0.002249,0.024239,0.007843),(0.004547,0.003392,0.030909,0.011765),(0.006006,0.004692,0.038558,0.015686),(0.007676,0.006136,0.046836,0.019608),(0.009561,0.007713,0.055143,0.023529),(0.011663,0.009417,0.06346,0.027451),(0.013995,0.011225,0.071862,0.031373),(0.016561,0.013136,0.080282,0.035294),(0.019373,0.015133,0.088767,0.039216),(0.022447,0.017199,0.097327,0.043137),(0.025793,0.019331,0.10593,0.047059),(0.029432,0.021503,0.114621,0.05098),(0.033385,0.023702,0.123397,0.054902),(0.037668,0.025921,0.132232,0.058824),(0.042253,0.028139,0.141141,0.062745),(0.046915,0.030324,0.150164,0.066667),(0.051644,0.032474,0.159254,0.070588),(0.056449,0.034569,0.168414,0.07451),(0.06134,0.03659,0.177642,0.078431),(0.066331,0.038504,0.186962,0.082353),(0.071429,0.040294,0.196354,0.086275),(0.076637,0.041905,0.205799,0.090196),(0.081962,0.043328,0.215289,0.094118),(0.087411,0.044556,0.224813,0.098039),(0.09299,0.045583,0.234358,0.101961),(0.098702,0.046402,0.243904,0.105882),(0.104551,0.047008,0.25343,0.109804),(0.110536,0.047399,0.262912,0.113725),(0.116656,0.047574,0.272321,0.117647),(0.122908,0.047536,0.281624,0.121569),(0.129285,0.047293,0.290788,0.12549),(0.135778,0.046856,0.299776,0.129412),(0.142378,0.046242,0.308553,0.133333),(0.149073,0.045468,0.317085,0.137255),(0.15585,0.044559,0.325338,0.141176),(0.162689,0.043554,0.333277,0.145098),(0.169575,0.042489,0.340874,0.14902),(0.176493,0.041402,0.348111,0.152941),(0.183429,0.040329,0.354971,0.156863),(0.190367,0.039309,0.361447,0.160784),(0.197297,0.0384,0.367535,0.164706),(0.204209,0.037632,0.373238,0.168627),(0.211095,0.03703,0.378563,0.172549),(0.217949,0.036615,0.383522,0.176471),(0.224763,0.036405,0.388129,0.180392),(0.231538,0.036405,0.3924,0.184314),(0.238273,0.036621,0.396353,0.188235),(0.244967,0.037055,0.400007,0.192157),(0.25162,0.037705,0.403378,0.196078),(0.258234,0.038571,0.406485,0.2),(0.26481,0.039647,0.409345,0.203922),(0.271347,0.040922,0.411976,0.207843),(0.27785,0.042353,0.414392,0.211765),(0.284321,0.043933,0.416608,0.215686),(0.290763,0.045644,0.418637,0.219608),(0.297178,0.04747,0.420491,0.223529),(0.303568,0.049396,0.422182,0.227451),(0.309935,0.051407,0.423721,0.231373),(0.316282,0.05349,0.425116,0.235294),(0.32261,0.055634,0.426377,0.239216),(0.328921,0.057827,0.427511,0.243137),(0.335217,0.06006,0.428524,0.247059),(0.3415,0.062325,0.429425,0.25098),(0.347771,0.064616,0.430217,0.254902),(0.354032,0.066925,0.430906,0.258824),(0.360284,0.069247,0.431497,0.262745),(0.366529,0.071579,0.431994,0.266667),(0.372768,0.073915,0.4324,0.270588),(0.379001,0.076253,0.432719,0.27451),(0.385228,0.078591,0.432955,0.278431),(0.391453,0.080927,0.433109,0.282353),(0.397674,0.083257,0.433183,0.286275),(0.403894,0.08558,0.433179,0.290196),(0.410113,0.087896,0.433098,0.294118),(0.416331,0.090203,0.432943,0.298039),(0.422549,0.092501,0.432714,0.301961),(0.428768,0.09479,0.432412,0.305882),(0.434987,0.097069,0.432039,0.309804),(0.441207,0.099338,0.431594,0.313725),(0.447428,0.101597,0.43108,0.317647),(0.453651,0.103848,0.430498,0.321569),(0.459875,0.106089,0.429846,0.32549),(0.4661,0.108322,0.429125,0.329412),(0.472328,0.110547,0.428334,0.333333),(0.478558,0.112764,0.427475,0.337255),(0.484789,0.114974,0.426548,0.341176),(0.491022,0.117179,0.425552,0.345098),(0.497257,0.119379,0.424488,0.34902),(0.503493,0.121575,0.423356,0.352941),(0.50973,0.123769,0.422156,0.356863),(0.515967,0.12596,0.420887,0.360784),(0.522206,0.12815,0.419549,0.364706),(0.528444,0.130341,0.418142,0.368627),(0.534683,0.132534,0.416667,0.372549),(0.54092,0.134729,0.415123,0.376471),(0.547157,0.136929,0.413511,0.380392),(0.553392,0.139134,0.411829,0.384314),(0.559624,0.141346,0.410078,0.388235),(0.565854,0.143567,0.408258,0.392157),(0.572081,0.145797,0.406369,0.396078),(0.578304,0.148039,0.404411,0.4),(0.584521,0.150294,0.402385,0.403922),(0.590734,0.152563,0.40029,0.407843),(0.59694,0.154848,0.398125,0.411765),(0.603139,0.157151,0.395891,0.415686),(0.60933,0.159474,0.393589,0.419608),(0.615513,0.161817,0.391219,0.423529),(0.621685,0.164184,0.388781,0.427451),(0.627847,0.166575,0.386276,0.431373),(0.633998,0.168992,0.383704,0.435294),(0.640135,0.171438,0.381065,0.439216),(0.64626,0.173914,0.378359,0.443137),(0.652369,0.176421,0.375586,0.447059),(0.658463,0.178962,0.372748,0.45098),(0.66454,0.181539,0.369846,0.454902),(0.670599,0.184153,0.366879,0.458824),(0.676638,0.186807,0.363849,0.462745),(0.682656,0.189501,0.360757,0.466667),(0.688653,0.192239,0.357603,0.470588),(0.694627,0.195021,0.354388,0.47451),(0.700576,0.197851,0.351113,0.478431),(0.7065,0.200728,0.347777,0.482353),(0.712396,0.203656,0.344383,0.486275),(0.718264,0.206636,0.340931,0.490196),(0.724103,0.20967,0.337424,0.494118),(0.729909,0.212759,0.333861,0.498039),(0.735683,0.215906,0.330245,0.501961),(0.741423,0.219112,0.326576,0.505882),(0.747127,0.222378,0.322856,0.509804),(0.752794,0.225706,0.319085,0.513725),(0.758422,0.229097,0.315266,0.517647),(0.76401,0.232554,0.311399,0.521569),(0.769556,0.236077,0.307485,0.52549),(0.775059,0.239667,0.303526,0.529412),(0.780517,0.243327,0.299523,0.533333),(0.785929,0.247056,0.295477,0.537255),(0.791293,0.250856,0.29139,0.541176),(0.796607,0.254728,0.287264,0.545098),(0.801871,0.258674,0.283099,0.54902),(0.807082,0.262692,0.278898,0.552941),(0.812239,0.266786,0.274661,0.556863),(0.817341,0.270954,0.27039,0.560784),(0.822386,0.275197,0.266085,0.564706),(0.827372,0.279517,0.26175,0.568627),(0.832299,0.283913,0.257383,0.572549),(0.837165,0.288385,0.252988,0.576471),(0.841969,0.292933,0.248564,0.580392),(0.846709,0.297559,0.244113,0.584314),(0.851384,0.30226,0.239636,0.588235),(0.855992,0.307038,0.235133,0.592157),(0.860533,0.311892,0.230606,0.596078),(0.865006,0.316822,0.226055,0.6),(0.869409,0.321827,0.221482,0.603922),(0.873741,0.326906,0.216886,0.607843),(0.878001,0.33206,0.212268,0.611765),(0.882188,0.337287,0.207628,0.615686),(0.886302,0.342586,0.202968,0.619608),(0.890341,0.347957,0.198286,0.623529),(0.894305,0.353399,0.193584,0.627451),(0.898192,0.358911,0.18886,0.631373),(0.902003,0.364492,0.184116,0.635294),(0.905735,0.37014,0.17935,0.639216),(0.90939,0.375856,0.174563,0.643137),(0.912966,0.381636,0.169755,0.647059),(0.916462,0.387481,0.164924,0.65098),(0.919879,0.393389,0.16007,0.654902),(0.923215,0.399359,0.155193,0.658824),(0.92647,0.405389,0.150292,0.662745),(0.929644,0.411479,0.145367,0.666667),(0.932737,0.417627,0.140417,0.670588),(0.935747,0.423831,0.13544,0.67451),(0.938675,0.430091,0.130438,0.678431),(0.941521,0.436405,0.125409,0.682353),(0.944285,0.442772,0.120354,0.686275),(0.946965,0.449191,0.115272,0.690196),(0.949562,0.45566,0.110164,0.694118),(0.952075,0.462178,0.105031,0.698039),(0.954506,0.468744,0.099874,0.701961),(0.956852,0.475356,0.094695,0.705882),(0.959114,0.482014,0.089499,0.709804),(0.961293,0.488716,0.084289,0.713725),(0.963387,0.495462,0.079073,0.717647),(0.965397,0.502249,0.073859,0.721569),(0.967322,0.509078,0.068659,0.72549),(0.969163,0.515946,0.063488,0.729412),(0.970919,0.522853,0.058367,0.733333),(0.97259,0.529798,0.053324,0.737255),(0.974176,0.53678,0.048392,0.741176),(0.975677,0.543798,0.043618,0.745098),(0.977092,0.55085,0.03905,0.74902),(0.978422,0.557937,0.034931,0.752941),(0.979666,0.565057,0.031409,0.756863),(0.980824,0.572209,0.028508,0.760784),(0.981895,0.579392,0.02625,0.764706),(0.982881,0.586606,0.024661,0.768627),(0.983779,0.593849,0.02377,0.772549),(0.984591,0.601122,0.023606,0.776471),(0.985315,0.608422,0.024202,0.780392),(0.985952,0.61575,0.025592,0.784314),(0.986502,0.623105,0.027814,0.788235),(0.986964,0.630485,0.030908,0.792157),(0.987337,0.63789,0.034916,0.796078),(0.987622,0.64532,0.039886,0.8),(0.987819,0.652773,0.045581,0.803922),(0.987926,0.66025,0.05175,0.807843),(0.987945,0.667748,0.058329,0.811765),(0.987874,0.675267,0.065257,0.815686),(0.987714,0.682807,0.072489,0.819608),(0.987464,0.690366,0.07999,0.823529),(0.987124,0.697944,0.087731,0.827451),(0.986694,0.70554,0.095694,0.831373),(0.986175,0.713153,0.103863,0.835294),(0.985566,0.720782,0.112229,0.839216),(0.984865,0.728427,0.120785,0.843137),(0.984075,0.736087,0.129527,0.847059),(0.983196,0.743758,0.138453,0.85098),(0.982228,0.751442,0.147565,0.854902),(0.981173,0.759135,0.156863,0.858824),(0.980032,0.766837,0.166353,0.862745),(0.978806,0.774545,0.176037,0.866667),(0.977497,0.782258,0.185923,0.870588),(0.976108,0.789974,0.196018,0.87451),(0.974638,0.797692,0.206332,0.878431),(0.973088,0.805409,0.216877,0.882353),(0.971468,0.813122,0.227658,0.886275),(0.969783,0.820825,0.238686,0.890196),(0.968041,0.828515,0.249972,0.894118),(0.966243,0.836191,0.261534,0.898039),(0.964394,0.843848,0.273391,0.901961),(0.962517,0.851476,0.285546,0.905882),(0.960626,0.859069,0.29801,0.909804),(0.95872,0.866624,0.31082,0.913725),(0.956834,0.874129,0.323974,0.917647),(0.954997,0.881569,0.337475,0.921569),(0.953215,0.888942,0.351369,0.92549),(0.951546,0.896226,0.365627,0.929412),(0.950018,0.903409,0.380271,0.933333),(0.948683,0.910473,0.395289,0.937255),(0.947594,0.917399,0.410665,0.941176),(0.946809,0.924168,0.426373,0.945098),(0.946392,0.930761,0.442367,0.94902),(0.946403,0.937159,0.458592,0.952941),(0.946903,0.943348,0.47497,0.956863),(0.947937,0.949318,0.491426,0.960784),(0.949545,0.955063,0.50786,0.964706),(0.95174,0.960587,0.524203,0.968627),(0.954529,0.965896,0.540361,0.972549),(0.957896,0.971003,0.556275,0.976471),(0.961812,0.975924,0.571925,0.980392),(0.966249,0.980678,0.587206,0.984314),(0.971162,0.985282,0.602154,0.988235),(0.976511,0.989753,0.61676,0.992157),(0.982257,0.994109,0.631017,0.996078),(0.988362,0.998364,0.644924,1.0)])
  elseif str == "magma"
    return cmap([(0.001462,0.000466,0.013866,0.0),(0.002258,0.001295,0.018331,0.003922),(0.003279,0.002305,0.023708,0.007843),(0.004512,0.00349,0.029965,0.011765),(0.00595,0.004843,0.03713,0.015686),(0.007588,0.006356,0.044973,0.019608),(0.009426,0.008022,0.052844,0.023529),(0.011465,0.009828,0.06075,0.027451),(0.013708,0.011771,0.068667,0.031373),(0.016156,0.01384,0.076603,0.035294),(0.018815,0.016026,0.084584,0.039216),(0.021692,0.01832,0.09261,0.043137),(0.024792,0.020715,0.100676,0.047059),(0.028123,0.023201,0.108787,0.05098),(0.031696,0.025765,0.116965,0.054902),(0.03552,0.028397,0.125209,0.058824),(0.039608,0.03109,0.133515,0.062745),(0.04383,0.03383,0.141886,0.066667),(0.048062,0.036607,0.150327,0.070588),(0.05232,0.039407,0.158841,0.07451),(0.056615,0.04216,0.167446,0.078431),(0.060949,0.044794,0.176129,0.082353),(0.06533,0.047318,0.184892,0.086275),(0.069764,0.049726,0.193735,0.090196),(0.074257,0.052017,0.20266,0.094118),(0.078815,0.054184,0.211667,0.098039),(0.083446,0.056225,0.220755,0.101961),(0.088155,0.058133,0.229922,0.105882),(0.092949,0.059904,0.239164,0.109804),(0.097833,0.061531,0.248477,0.113725),(0.102815,0.06301,0.257854,0.117647),(0.107899,0.064335,0.267289,0.121569),(0.113094,0.065492,0.276784,0.12549),(0.118405,0.066479,0.286321,0.129412),(0.123833,0.067295,0.295879,0.133333),(0.12938,0.067935,0.305443,0.137255),(0.135053,0.068391,0.315,0.141176),(0.140858,0.068654,0.324538,0.145098),(0.146785,0.068738,0.334011,0.14902),(0.152839,0.068637,0.343404,0.152941),(0.159018,0.068354,0.352688,0.156863),(0.165308,0.067911,0.361816,0.160784),(0.171713,0.067305,0.370771,0.164706),(0.178212,0.066576,0.379497,0.168627),(0.184801,0.065732,0.387973,0.172549),(0.19146,0.064818,0.396152,0.176471),(0.198177,0.063862,0.404009,0.180392),(0.204935,0.062907,0.411514,0.184314),(0.211718,0.061992,0.418647,0.188235),(0.218512,0.061158,0.425392,0.192157),(0.225302,0.060445,0.431742,0.196078),(0.232077,0.059889,0.437695,0.2),(0.238826,0.059517,0.443256,0.203922),(0.245543,0.059352,0.448436,0.207843),(0.25222,0.059415,0.453248,0.211765),(0.258857,0.059706,0.45771,0.215686),(0.265447,0.060237,0.46184,0.219608),(0.271994,0.060994,0.46566,0.223529),(0.278493,0.061978,0.46919,0.227451),(0.284951,0.063168,0.472451,0.231373),(0.291366,0.064553,0.475462,0.235294),(0.29774,0.066117,0.478243,0.239216),(0.304081,0.067835,0.480812,0.243137),(0.310382,0.069702,0.483186,0.247059),(0.316654,0.07169,0.48538,0.25098),(0.322899,0.073782,0.487408,0.254902),(0.329114,0.075972,0.489287,0.258824),(0.335308,0.078236,0.491024,0.262745),(0.341482,0.080564,0.492631,0.266667),(0.347636,0.082946,0.494121,0.270588),(0.353773,0.085373,0.495501,0.27451),(0.359898,0.087831,0.496778,0.278431),(0.366012,0.090314,0.49796,0.282353),(0.372116,0.092816,0.499053,0.286275),(0.378211,0.095332,0.500067,0.290196),(0.384299,0.097855,0.501002,0.294118),(0.390384,0.100379,0.501864,0.298039),(0.396467,0.102902,0.502658,0.301961),(0.402548,0.10542,0.503386,0.305882),(0.408629,0.10793,0.504052,0.309804),(0.414709,0.110431,0.504662,0.313725),(0.420791,0.11292,0.505215,0.317647),(0.426877,0.115395,0.505714,0.321569),(0.432967,0.117855,0.50616,0.32549),(0.439062,0.120298,0.506555,0.329412),(0.445163,0.122724,0.506901,0.333333),(0.451271,0.125132,0.507198,0.337255),(0.457386,0.127522,0.507448,0.341176),(0.463508,0.129893,0.507652,0.345098),(0.46964,0.132245,0.507809,0.34902),(0.47578,0.134577,0.507921,0.352941),(0.481929,0.136891,0.507989,0.356863),(0.488088,0.139186,0.508011,0.360784),(0.494258,0.141462,0.507988,0.364706),(0.500438,0.143719,0.50792,0.368627),(0.506629,0.145958,0.507806,0.372549),(0.512831,0.148179,0.507648,0.376471),(0.519045,0.150383,0.507443,0.380392),(0.52527,0.152569,0.507192,0.384314),(0.531507,0.154739,0.506895,0.388235),(0.537755,0.156894,0.506551,0.392157),(0.544015,0.159033,0.506159,0.396078),(0.550287,0.161158,0.505719,0.4),(0.556571,0.163269,0.50523,0.403922),(0.562866,0.165368,0.504692,0.407843),(0.569172,0.167454,0.504105,0.411765),(0.57549,0.16953,0.503466,0.415686),(0.581819,0.171596,0.502777,0.419608),(0.588158,0.173652,0.502035,0.423529),(0.594508,0.175701,0.501241,0.427451),(0.600868,0.177743,0.500394,0.431373),(0.607238,0.179779,0.499492,0.435294),(0.613617,0.181811,0.498536,0.439216),(0.620005,0.18384,0.497524,0.443137),(0.626401,0.185867,0.496456,0.447059),(0.632805,0.187893,0.495332,0.45098),(0.639216,0.189921,0.49415,0.454902),(0.645633,0.191952,0.49291,0.458824),(0.652056,0.193986,0.491611,0.462745),(0.658483,0.196027,0.490253,0.466667),(0.664915,0.198075,0.488836,0.470588),(0.671349,0.200133,0.487358,0.47451),(0.677786,0.202203,0.485819,0.478431),(0.684224,0.204286,0.484219,0.482353),(0.690661,0.206384,0.482558,0.486275),(0.697098,0.208501,0.480835,0.490196),(0.703532,0.210638,0.479049,0.494118),(0.709962,0.212797,0.477201,0.498039),(0.716387,0.214982,0.47529,0.501961),(0.722805,0.217194,0.473316,0.505882),(0.729216,0.219437,0.471279,0.509804),(0.735616,0.221713,0.46918,0.513725),(0.742004,0.224025,0.467018,0.517647),(0.748378,0.226377,0.464794,0.521569),(0.754737,0.228772,0.462509,0.52549),(0.761077,0.231214,0.460162,0.529412),(0.767398,0.233705,0.457755,0.533333),(0.773695,0.236249,0.455289,0.537255),(0.779968,0.238851,0.452765,0.541176),(0.786212,0.241514,0.450184,0.545098),(0.792427,0.244242,0.447543,0.54902),(0.798608,0.24704,0.444848,0.552941),(0.804752,0.249911,0.442102,0.556863),(0.810855,0.252861,0.439305,0.560784),(0.816914,0.255895,0.436461,0.564706),(0.822926,0.259016,0.433573,0.568627),(0.828886,0.262229,0.430644,0.572549),(0.834791,0.26554,0.427671,0.576471),(0.840636,0.268953,0.424666,0.580392),(0.846416,0.272473,0.421631,0.584314),(0.852126,0.276106,0.418573,0.588235),(0.857763,0.279857,0.415496,0.592157),(0.86332,0.283729,0.412403,0.596078),(0.868793,0.287728,0.409303,0.6),(0.874176,0.291859,0.406205,0.603922),(0.879464,0.296125,0.403118,0.607843),(0.884651,0.30053,0.400047,0.611765),(0.889731,0.305079,0.397002,0.615686),(0.8947,0.309773,0.393995,0.619608),(0.899552,0.314616,0.391037,0.623529),(0.904281,0.31961,0.388137,0.627451),(0.908884,0.324755,0.385308,0.631373),(0.913354,0.330052,0.382563,0.635294),(0.917689,0.3355,0.379915,0.639216),(0.921884,0.341098,0.377376,0.643137),(0.925937,0.346844,0.374959,0.647059),(0.929845,0.352734,0.372677,0.65098),(0.933606,0.358764,0.370541,0.654902),(0.937221,0.364929,0.368567,0.658824),(0.940687,0.371224,0.366762,0.662745),(0.944006,0.377643,0.365136,0.666667),(0.94718,0.384178,0.363701,0.670588),(0.95021,0.39082,0.362468,0.67451),(0.953099,0.397563,0.361438,0.678431),(0.955849,0.4044,0.360619,0.682353),(0.958464,0.411324,0.360014,0.686275),(0.960949,0.418323,0.35963,0.690196),(0.96331,0.42539,0.359469,0.694118),(0.965549,0.432519,0.359529,0.698039),(0.967671,0.439703,0.35981,0.701961),(0.96968,0.446936,0.360311,0.705882),(0.971582,0.45421,0.36103,0.709804),(0.973381,0.46152,0.361965,0.713725),(0.975082,0.468861,0.363111,0.717647),(0.97669,0.476226,0.364466,0.721569),(0.97821,0.483612,0.366025,0.72549),(0.979645,0.491014,0.367783,0.729412),(0.981,0.498428,0.369734,0.733333),(0.982279,0.505851,0.371874,0.737255),(0.983485,0.51328,0.374198,0.741176),(0.984622,0.520713,0.376698,0.745098),(0.985693,0.528148,0.379371,0.74902),(0.9867,0.535582,0.38221,0.752941),(0.987646,0.543015,0.38521,0.756863),(0.988533,0.550446,0.388365,0.760784),(0.989363,0.557873,0.391671,0.764706),(0.990138,0.565296,0.395122,0.768627),(0.990871,0.572706,0.398714,0.772549),(0.991558,0.580107,0.402441,0.776471),(0.992196,0.587502,0.406299,0.780392),(0.992785,0.594891,0.410283,0.784314),(0.993326,0.602275,0.41439,0.788235),(0.993834,0.609644,0.418613,0.792157),(0.994309,0.616999,0.42295,0.796078),(0.994738,0.62435,0.427397,0.8),(0.995122,0.631696,0.431951,0.803922),(0.99548,0.639027,0.436607,0.807843),(0.99581,0.646344,0.441361,0.811765),(0.996096,0.653659,0.446213,0.815686),(0.996341,0.660969,0.45116,0.819608),(0.99658,0.668256,0.456192,0.823529),(0.996775,0.675541,0.461314,0.827451),(0.996925,0.682828,0.466526,0.831373),(0.997077,0.690088,0.471811,0.835294),(0.997186,0.697349,0.477182,0.839216),(0.997254,0.704611,0.482635,0.843137),(0.997325,0.711848,0.488154,0.847059),(0.997351,0.719089,0.493755,0.85098),(0.997351,0.726324,0.499428,0.854902),(0.997341,0.733545,0.505167,0.858824),(0.997285,0.740772,0.510983,0.862745),(0.997228,0.747981,0.516859,0.866667),(0.997138,0.75519,0.522806,0.870588),(0.997019,0.762398,0.528821,0.87451),(0.996898,0.769591,0.534892,0.878431),(0.996727,0.776795,0.541039,0.882353),(0.996571,0.783977,0.547233,0.886275),(0.996369,0.791167,0.553499,0.890196),(0.996162,0.798348,0.55982,0.894118),(0.995932,0.805527,0.566202,0.898039),(0.99568,0.812706,0.572645,0.901961),(0.995424,0.819875,0.57914,0.905882),(0.995131,0.827052,0.585701,0.909804),(0.994851,0.834213,0.592307,0.913725),(0.994524,0.841387,0.598983,0.917647),(0.994222,0.84854,0.605696,0.921569),(0.993866,0.855711,0.612482,0.92549),(0.993545,0.862859,0.619299,0.929412),(0.99317,0.870024,0.626189,0.933333),(0.992831,0.877168,0.633109,0.937255),(0.99244,0.88433,0.640099,0.941176),(0.992089,0.89147,0.647116,0.945098),(0.991688,0.898627,0.654202,0.94902),(0.991332,0.905763,0.661309,0.952941),(0.99093,0.912915,0.668481,0.956863),(0.99057,0.920049,0.675675,0.960784),(0.990175,0.927196,0.682926,0.964706),(0.989815,0.934329,0.690198,0.968627),(0.989434,0.94147,0.697519,0.972549),(0.989077,0.948604,0.704863,0.976471),(0.988717,0.955742,0.712242,0.980392),(0.988367,0.962878,0.719649,0.984314),(0.988033,0.970012,0.727077,0.988235),(0.987691,0.977154,0.734536,0.992157),(0.987387,0.984288,0.742002,0.996078),(0.987053,0.991438,0.749504,1.0)])
  else
    error("Colormap ", str, " not implemented!")
  end
end

function existing_cmaps()
  ["gray", "blue", "green", "red", "UKE", "redgreen", "redgreenalpha",
   "greenorangealpha", "greenblueredalpha", "bluegreenredalpha",
   "yellowredalpha", "tobi", "tobitr", "redwhite","timemap","UPthreshold",
   "perfusionmap","magma","inferno","plasma","Liver","viridis"]
end


text:
#-------------------------------------------------------------------------------
# Test suite for 'utils.jl' file
#-------------------------------------------------------------------------------

@testset "utils" begin

    @testset "verify_initial_room" begin

        # Test parameters
        n = 2
        Δ = 2.0
        a = [1.0, 3.0]
        b = [5.0, 10.0]
        # Test
        @test( LOWDER.verify_initial_room(n, Δ, a, b) == true )

        # Test parameters
        n = 3
        Δ = 1.0
        a = [0.0, 0.0, 0.0]
        b = [1.0, 5.0, 12.0]
        # Test
        @test( LOWDER.verify_initial_room(n, Δ, a, b) == false )

        # Test parameters
        n = 5
        Δ = 3.0
        a = [- 2.0, - 4.0, 0.0, 1.0, 1.5]
        b = [1.0, - 2.0, 1.0, 2.5, 6.0]
        # Test
        @test( LOWDER.verify_initial_room(n, Δ, a, b) == false )

    end

    @testset "correct_guess_bounds!" begin

        # Test parameters
        n = 2
        δ = 2.0
        a = [1.0, 3.0]
        b = [5.0, 10.0]
        x = [3.0, 5.0]
        ao = [0.0, 0.0]
        bo = [0.0, 0.0]
        # Expected return
        x_return = [3.0, 5.0]
        ao_return = [- 2.0, - 2.0]
        bo_return = [2.0, 5.0]
        # Test
        LOWDER.correct_guess_bounds!(n, δ, a, b, x, ao, bo)
        @test( x == x_return )
        @test( ao == ao_return)
        @test( bo == bo_return)

        # Test parameters
        n = 2
        δ = 2.0
        a = [1.0, 3.0]
        b = [5.0, 10.0]
        x = [0.0, 0.0]
        ao = [0.0, 0.0]
        bo = [0.0, 0.0]
        # Expected return
        x_return = [1.0, 3.0]
        ao_return = [0.0, 0.0]
        bo_return = [4.0, 7.0]
        # Test
        LOWDER.correct_guess_bounds!(n, δ, a, b, x, ao, bo)
        @test( x == x_return )
        @test( ao == ao_return)
        @test( bo == bo_return)

        # Test parameters
        n = 2
        δ = 2.0
        a = [1.0, 3.0]
        b = [5.0, 10.0]
        x = [10.0, 10.0]
        ao = [0.0, 0.0]
        bo = [0.0, 0.0]
        # Expected return
        x_return = [5.0, 10.0]
        ao_return = [- 4.0, - 7.0]
        bo_return = [0.0, 0.0]
        # Test
        LOWDER.correct_guess_bounds!(n, δ, a, b, x, ao, bo)
        @test( x == x_return )
        @test( ao == ao_return)
        @test( bo == bo_return)

        # Test parameters
        n = 2
        δ = 2.0
        a = [1.0, 3.0]
        b = [5.0, 10.0]
        x = [0.0, 20.0]
        ao = [0.0, 0.0]
        bo = [0.0, 0.0]
        # Expected return
        x_return = [1.0, 10.0]
        ao_return = [0.0, - 7.0]
        bo_return = [4.0, 0.0]
        # Test
        LOWDER.correct_guess_bounds!(n, δ, a, b, x, ao, bo)
        @test( x == x_return )
        @test( ao == ao_return)
        @test( bo == bo_return)

        # Test parameters
        n = 2
        δ = 2.0
        a = [1.0, 3.0]
        b = [5.0, 10.0]
        x = [2.0, 5.0]
        ao = [0.0, 0.0]
        bo = [0.0, 0.0]
        # Expected return
        x_return = [3.0, 5.0]
        ao_return = [- 2.0, - 2.0]
        bo_return = [2.0, 5.0]
        # Test
        LOWDER.correct_guess_bounds!(n, δ, a, b, x, ao, bo)
        @test( x == x_return )
        @test( ao == ao_return)
        @test( bo == bo_return)

        # Test parameters
        n = 2
        δ = 2.0
        a = [1.0, 3.0]
        b = [5.0, 10.0]
        x = [3.0, 4.0]
        ao = [0.0, 0.0]
        bo = [0.0, 0.0]
        # Expected return
        x_return = [3.0, 5.0]
        ao_return = [- 2.0, - 2.0]
        bo_return = [2.0, 5.0]
        # Test
        LOWDER.correct_guess_bounds!(n, δ, a, b, x, ao, bo)
        @test( x == x_return )
        @test( ao == ao_return)
        @test( bo == bo_return)

        # Test parameters
        n = 2
        δ = 2.0
        a = [1.0, 3.0]
        b = [5.0, 10.0]
        x = [2.0, 4.0]
        ao = [0.0, 0.0]
        bo = [0.0, 0.0]
        # Expected return
        x_return = [3.0, 5.0]
        ao_return = [- 2.0, - 2.0]
        bo_return = [2.0, 5.0]
        # Test
        LOWDER.correct_guess_bounds!(n, δ, a, b, x, ao, bo)
        @test( x == x_return )
        @test( ao == ao_return)
        @test( bo == bo_return)

        # Test parameters
        n = 2
        δ = 2.0
        a = [1.0, 3.0]
        b = [5.0, 10.0]
        x = [4.0, 5.0]
        ao = [0.0, 0.0]
        bo = [0.0, 0.0]
        # Expected return
        x_return = [3.0, 5.0]
        ao_return = [- 2.0, - 2.0]
        bo_return = [2.0, 5.0]
        # Test
        LOWDER.correct_guess_bounds!(n, δ, a, b, x, ao, bo)
        @test( x == x_return )
        @test( ao == ao_return)
        @test( bo == bo_return)

        # Test parameters
        n = 2
        δ = 2.0
        a = [1.0, 3.0]
        b = [5.0, 10.0]
        x = [3.0, 9.0]
        ao = [0.0, 0.0]
        bo = [0.0, 0.0]
        # Expected return
        x_return = [3.0, 8.0]
        ao_return = [- 2.0, - 5.0]
        bo_return = [2.0, 2.0]
        # Test
        LOWDER.correct_guess_bounds!(n, δ, a, b, x, ao, bo)
        @test( x == x_return )
        @test( ao == ao_return)
        @test( bo == bo_return)

        # Test parameters
        n = 2
        δ = 2.0
        a = [1.0, 3.0]
        b = [5.0, 10.0]
        x = [4.0, 9.0]
        ao = [0.0, 0.0]
        bo = [0.0, 0.0]
        # Expected return
        x_return = [3.0, 8.0]
        ao_return = [- 2.0, - 5.0]
        bo_return = [2.0, 2.0]
        # Test
        LOWDER.correct_guess_bounds!(n, δ, a, b, x, ao, bo)
        @test( x == x_return )
        @test( ao == ao_return)
        @test( bo == bo_return)

    end

end

max_stars_repo_name: aschwertner/LOWDER

max_issues_repo_name: aschwertner/LOWDER

max_forks_repo_name: aschwertner/LOWDER


text:
# ---
# title: 670. Maximum Swap
# id: problem670
# author: Indigo
# date: 2021-06-28
# difficulty: Medium
# categories: Array, Math
# link: <https://leetcode.com/problems/maximum-swap/description/>
# hidden: true
# ---
#
# Given a non-negative integer, you could swap two digits **at most** once to
# get the maximum valued number. Return the maximum valued number you could get.
#
# **Example 1:**
#
#
#
#     Input: 2736
#     Output: 7236
#     Explanation: Swap the number 2 and the number 7.
#
#
# **Example 2:**
#
#
#
#     Input: 9973
#     Output: 9973
#     Explanation: No swap.
#
#
# **Note:**
#
#   1. The given number is in the range [0, 108]
#
#
## @lc code=start
using LeetCode

function maximum_swap670(nums::Int)
    dgs = digits(nums)
    for i in 1:length(dgs)
        cur = @view(dgs[i:end])
        max_idx = argmax(cur)
        max_idx == lastindex(cur) && continue
        dgs[end], dgs[max_idx] = dgs[max_idx], dgs[end]
        break
    end
    return [10 .^ (0:length(dgs)-1);]' * dgs
end
## @lc code=end


text:
# Autogenerated wrapper script for boostpython_jll for x86_64-linux-musl-cxx03
export libboost_python

using boost_jll
using Python_jll
JLLWrappers.@generate_wrapper_header("boostpython")
JLLWrappers.@declare_library_product(libboost_python, "libboost_python38.so")
function __init__()
    JLLWrappers.@generate_init_header(boost_jll, Python_jll)
    JLLWrappers.@init_library_product(
        libboost_python,
        "lib/libboost_python38.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Test, DeviationsLH

dlh = DeviationsLH;

function change_table_test()
    @testset "Change table basics" begin
        n = 3;
        dv = make_deviation_vector(n);

        paramNameV = [:aa, :bbb, :cc, :d];
        nParam = length(paramNameV);

        ct = dlh.ChangeTable(dv, paramNameV);
        @test dlh.param_names(ct) == paramNameV
        @test dlh.n_params(ct) == nParam
        @test dlh.n_devs(ct) == length(dv)

        for j = 1 : nParam
            # Offset ensures that the first deviation is the same as the base
            dv2 = make_deviation_vector(n; offset = 0.5 * (j-1));
            dlh.set_param_values!(ct, j, dv2; scalarDev = nothing);
        end

        dlh.show_table(ct)
        for transposed = [false, true]
            dlh.show_largest_change_table(ct, 3; transposed = transposed)
        end

        pNameV = dlh.find_unchanged_devs(ct; rtol = 0.02);
        println("Unchanged deviations: $pNameV")
        @test isa(pNameV, typeof(paramNameV))
        @test length(pNameV) >= 1
	end
end

@testset "Change Table" begin
    change_table_test()
end

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