text:
using GPUCompiler

if !@isdefined(TestRuntime)
    include("../util.jl")
end

# create a native test compiler, and generate reflection methods for it

function native_job(@nospecialize(func), @nospecialize(types); kernel::Bool=false, kwargs...)
    source = FunctionSpec(func, Base.to_tuple_type(types), kernel)
    target = NativeCompilerTarget(always_inline=true)
    params = TestCompilerParams()
    CompilerJob(target, source, params), kwargs
end

function native_code_typed(@nospecialize(func), @nospecialize(types); kwargs...)
    job, kwargs = native_job(func, types; kwargs...)
    GPUCompiler.code_typed(job; kwargs...)
end

function native_code_warntype(io::IO, @nospecialize(func), @nospecialize(types); kwargs...)
    job, kwargs = native_job(func, types; kwargs...)
    GPUCompiler.code_warntype(io, job; kwargs...)
end

function native_code_llvm(io::IO, @nospecialize(func), @nospecialize(types); kwargs...)
    job, kwargs = native_job(func, types; kwargs...)
    GPUCompiler.code_llvm(io, job; kwargs...)
end

function native_code_native(io::IO, @nospecialize(func), @nospecialize(types); kwargs...)
    job, kwargs = native_job(func, types; kwargs...)
    GPUCompiler.code_native(io, job; kwargs...)
end

# aliases without ::IO argument
for method in (:code_warntype, :code_llvm, :code_native)
    native_method = Symbol("native_$(method)")
    @eval begin
        $native_method(@nospecialize(func), @nospecialize(types); kwargs...) =
            $native_method(stdout, func, types; kwargs...)
    end
end

# simulates codegen for a kernel function: validates by default
function native_code_execution(@nospecialize(func), @nospecialize(types); kwargs...)
    job, kwargs = native_job(func, types; kernel=true, kwargs...)
    GPUCompiler.compile(:asm, job; kwargs...)
end

module LazyCodegen
    using LLVM
    using GPUCompiler

    import ..native_job

    # We have one global JIT and TM
    const jit = Ref{LLVM.OrcJIT}()
    const tm  = Ref{LLVM.TargetMachine}()

    function __init__()
        optlevel = LLVM.API.LLVMCodeGenLevelDefault

        tm[] = GPUCompiler.JITTargetMachine(optlevel=optlevel)
        LLVM.asm_verbosity!(tm[], true)

        jit[] = LLVM.OrcJIT(tm[]) # takes ownership of tm
        atexit() do
            dispose(jit[])
        end
    end

    import GPUCompiler: deferred_codegen_jobs, CompilerJob
    mutable struct CallbackContext
        job::CompilerJob
        stub::String
        compiled::Bool
    end

    const outstanding = IdDict{CallbackContext, Nothing}()

    # Setup the lazy callback for creating a module
    function callback(orc_ref::LLVM.API.LLVMOrcJITStackRef, callback_ctx::Ptr{Cvoid})
        orc = LLVM.OrcJIT(orc_ref)
        cc = Base.unsafe_pointer_to_objref(callback_ctx)::CallbackContext

        @assert !cc.compiled
        job = cc.job

        ir, entry = GPUCompiler.codegen(:llvm, job; validate=false)
        entry_name = name(entry)

        jitted_mod = compile!(orc, ir)

        addr = addressin(orc, jitted_mod, entry_name)
        ptr  = pointer(addr)

        cc.compiled = true
        delete!(outstanding, cc)

        # 4. Update the stub pointer to point to the recently compiled module
        set_stub!(orc, cc.stub, ptr)

        # 5. Return the address of the implementation, since we are going to call it now
        return UInt64(reinterpret(UInt, ptr))
    end

    @inline assume(cond::Bool) = Base.llvmcall(("""
    declare void @llvm.assume(i1)

    define void @entry(i8) #0 {
        %cond = icmp eq i8 %0, 1
        call void @llvm.assume(i1 %cond)
        ret void
    }

    attributes #0 = { alwaysinline }""", "entry"),
    Nothing, Tuple{Bool}, cond)

    @generated function deferred_codegen(::Val{f}, ::Val{tt}) where {f,tt}
        job, _ = native_job(f, tt)

        cc = CallbackContext(job, String(gensym(:trampoline)), false)
        outstanding[cc] = nothing

        c_callback = @cfunction(callback, UInt64, (LLVM.API.LLVMOrcJITStackRef, Ptr{Cvoid}))

        orc = jit[]
        initial_addr = callback!(orc, c_callback, pointer_from_objref(cc))
        create_stub!(orc, cc.stub, initial_addr)
        addr = address(orc, cc.stub)
        trampoline = pointer(addr)
        id = Base.reinterpret(Int, trampoline)

        deferred_codegen_jobs[id] = job

        quote
            ptr = ccall("extern deferred_codegen", llvmcall, Ptr{Cvoid}, (Ptr{Cvoid},), $trampoline)
            assume(ptr != C_NULL)
            return ptr
        end
    end

    @generated function abi_call(f::Ptr{Cvoid}, rt::Type{RT}, tt::Type{T}, args::Vararg{Any, N}) where {T, RT, N}
        argtt    = tt.parameters[1]
        rettype  = rt.parameters[1]
        argtypes = DataType[argtt.parameters...]

        argexprs = Union{Expr, Symbol}[]
        ccall_types = DataType[]

        before = :()
        after = :(ret)

        # Note this follows: emit_call_specfun_other
        LLVM.Interop.JuliaContext() do ctx
            T_jlvalue = LLVM.StructType(LLVMType[], ctx)
            T_prjlvalue = LLVM.PointerType(T_jlvalue, #= AddressSpace::Tracked =# 10)

            for (source_i, source_typ) in enumerate(argtypes)
                if GPUCompiler.isghosttype(source_typ) || Core.Compiler.isconstType(source_typ)
                    continue
                end

                argexpr = :(args[$source_i])

                isboxed = GPUCompiler.deserves_argbox(source_typ)
                et = isboxed ? T_prjlvalue : convert(LLVMType, source_typ, ctx)

                if isboxed
                    push!(ccall_types, Any)
                elseif isa(et, LLVM.SequentialType) # et->isAggregateType
                    push!(ccall_types, Ptr{source_typ})
                    argexpr = Expr(:call, GlobalRef(Base, :Ref), argexpr)
                else
                    push!(ccall_types, source_typ)
                end
                push!(argexprs, argexpr)
            end

            if GPUCompiler.isghosttype(rettype) || Core.Compiler.isconstType(rettype)
                # Do nothing...
                # In theory we could set `rettype` to `T_void`, but ccall will do that for us
            # elseif jl_is_uniontype?
            elseif !GPUCompiler.deserves_retbox(rettype)
                rt = convert(LLVMType, rettype, ctx)
                if !isa(rt, LLVM.VoidType) && GPUCompiler.deserves_sret(rettype, rt)
                    before = :(sret = Ref{$rettype}())
                    pushfirst!(argexprs, :(sret))
                    pushfirst!(ccall_types, Ptr{rettype})
                    rettype = Nothing
                    after = :(sret[])
                end
            else
                # rt = T_prjlvalue
            end
        end

        quote
            $before
            ret = ccall(f, $rettype, ($(ccall_types...),), $(argexprs...))
            $after
        end
    end

    @inline function call_delayed(f::F, args...) where F
        tt = Tuple{map(Core.Typeof, args)...}
        rt = Core.Compiler.return_type(f, tt)
        ptr = deferred_codegen(Val(f), Val(tt))
        abi_call(ptr, rt, tt, args...)
    end
end


text:
"""
$(SIGNATURES)

Maxwellian in discrete form

* @args: particle velocity quadrature points
* @args: density, velocity and inverse of temperature
* @return: Maxwellian distribution function

1V
"""
maxwellian(u, ρ, U, λ) = @. ρ * sqrt(λ / π) * exp(-λ * (u - U)^2)

"""
$(SIGNATURES)
"""
maxwellian(u, prim::AV) = maxwellian(u, prim[1], prim[2], prim[end]) # in case of input with length 4/5

"""
$(SIGNATURES)

2V
"""
maxwellian(u, v, ρ, U, V, λ) = @. ρ * (λ / π) * exp(-λ * ((u - U)^2 + (v - V)^2))

"""
$(SIGNATURES)
"""
maxwellian(u, v, prim::AV) = maxwellian(u, v, prim[1], prim[2], prim[3], prim[end]) # in case of input with length 5

"""
$(SIGNATURES)

3V
"""
maxwellian(u, v, w, ρ, U, V, W, λ) =
    @. ρ * sqrt((λ / π)^3) * exp(-λ * ((u - U)^2 + (v - V)^2 + (w - W)^2))

"""
$(SIGNATURES)
"""
maxwellian(u, v, w, prim::AV) =
    maxwellian(u, v, w, prim[1], prim[2], prim[3], prim[4], prim[5])

"""
$(SIGNATURES)

Reduced Maxwellian distribution related to energy
"""
energy_maxwellian(h, λ, K) = @. h * K / (2.0 * λ)

"""
$(SIGNATURES)
"""
energy_maxwellian(h, prim::AV, K) = energy_maxwellian(h, prim[end], K)

"""
$(SIGNATURES)

In-place Maxwellian

1V
"""
function maxwellian!(M::AV, u::AV, ρ, U, λ)
    @. M = ρ * sqrt(λ / π) * exp(-λ * (u - U)^2)
    return nothing
end

"""
$(SIGNATURES)
"""
maxwellian!(M::AV, u::AV, prim::AV) = maxwellian!(M, u, prim[1], prim[2], prim[end])

"""
$(SIGNATURES)

Rykov model
"""
function maxwellian!(
    Ht::T,
    Bt::T,
    Rt::T,
    Hr::T,
    Br::T,
    Rr::T,
    u::AV,
    prim::AV,
    K,
    Kr,
) where {T<:AV}

    @. Ht = prim[1] * sqrt(prim[4] / π) * exp(-prim[4] * (u - prim[2])^2)
    @. Bt = Ht * K / (2.0 * prim[4])
    @. Rt = Ht * Kr / (2.0 * prim[5])

    @. Hr = prim[1] * sqrt(prim[3] / π) * exp(-prim[3] * (u - prim[2])^2)
    @. Br = Hr * K / (2.0 * prim[3])
    @. Rr = Hr * Kr / (2.0 * prim[3])

    return nothing

end

"""
$(SIGNATURES)

2V
"""
function maxwellian!(M::AA, u::T, v::T, ρ, U, V, λ) where {T<:AA}
    @. M = ρ * (λ / π) * exp(-λ * ((u - U)^2 + (v - V)^2))
    return nothing
end

"""
$(SIGNATURES)
"""
maxwellian!(M::AA, u::T, v::T, prim::AV) where {T<:AA} =
    maxwellian!(M, u, v, prim[1], prim[2], prim[3], prim[end])

"""
$(SIGNATURES)

Rykov model
"""
function maxwellian!(
    Ht::T1,
    Bt::T1,
    Rt::T1,
    Hr::T1,
    Br::T1,
    Rr::T1,
    u::T2,
    v::T2,
    prim::AV,
    K,
    Kr,
) where {T1<:AA,T2<:AA}

    @. Ht = prim[1] * (prim[5] / π) * exp(-prim[5] * ((u - prim[2])^2 + (v - prim[3])^2))
    @. Bt = Ht * K / (2.0 * prim[5])
    @. Rt = Ht * Kr / (2.0 * prim[6])

    @. Hr = prim[1] * (prim[4] / π) * exp(-prim[4] * ((u - prim[2])^2 + (v - prim[3])^2))
    @. Br = Hr * K / (2.0 * prim[4])
    @. Rr = Hr * Kr / (2.0 * prim[4])

    return nothing

end

"""
$(SIGNATURES)

3V
"""
function maxwellian!(M::AA, u::T, v::T, w::T, ρ, U, V, W, λ) where {T<:AA}
    @. M = ρ * sqrt((λ / π)^3) * exp(-λ * ((u - U)^2 + (v - V)^2 + (w - W)^2))
    return nothing
end

"""
$(SIGNATURES)
"""
maxwellian!(M::AA, u::T, v::T, w::T, prim::AV) where {T<:AA} =
    maxwellian!(M, u, v, w, prim[1], prim[2], prim[3], prim[4], prim[5])

"""
$(SIGNATURES)

Multi-component Maxwellian in discrete form

1V
"""
function mixture_maxwellian(u::AM, prim::AM)
    mixM = similar(u)
    mixture_maxwellian!(mixM, u, prim)

    return mixM
end

"""
$(SIGNATURES)

2V
"""
function mixture_maxwellian(u::X, v::X, prim::AM) where {X<:AA}
    mixM = similar(u)
    mixture_maxwellian!(mixM, u, v, prim)

    return mixM
end

"""
$(SIGNATURES)

3V
"""
function mixture_maxwellian(u::X, v::X, w::X, prim::AM) where {X<:AA}
    mixM = similar(u)
    mixture_maxwellian!(mixM, u, v, w, prim)

    return mixM
end

"""
$(SIGNATURES)

In-place multi-component Maxwellian

1V
"""
function mixture_maxwellian!(M::AM, u::AA{T,2}, prim::AM) where {T}
    @inbounds for j in axes(M, 2)
        _M = @view M[:, j]
        maxwellian!(_M, u[:, j], prim[:, j])
    end

    return nothing
end

"""
$(SIGNATURES)

2V
"""
function mixture_maxwellian!(
    M::AA{T1,3},
    u::T2,
    v::T2,
    prim::AM,
) where {T1,T2<:AA{T3,3}} where {T3}

    @inbounds for k in axes(M, 3)
        _M = @view M[:, :, k]
        maxwellian!(_M, u[:, :, k], v[:, :, k], prim[:, k])
    end

    return nothing

end

"""
$(SIGNATURES)
"""
function mixture_maxwellian!(M::AM, u::T, v::T, prim::AM) where {T<:AM}

    @inbounds for k in axes(M, 2)
        _M = @view M[:, k]
        maxwellian!(_M, u[:, k], v[:, k], prim[:, k])
    end

    return nothing

end

"""
$(SIGNATURES)

3V
"""
function mixture_maxwellian!(
    M::AA{T1,4},
    u::T2,
    v::T2,
    w::T2,
    prim::AM,
) where {T1,T2<:AA{T3,4}} where {T3}

    @inbounds for l in axes(M, 4)
        _M = @view M[:, :, :, l]
        maxwellian!(_M, u[:, :, :, l], v[:, :, :, l], w[:, :, :, l], prim[:, l])
    end

    return nothing

end

function mixture_maxwellian!(M::AM, u::T, v::T, w::T, prim::AM) where {T<:AM}

    @inbounds for l in axes(M, 2)
        _M = @view M[:, l]
        maxwellian!(_M, u[:, l], v[:, l], w[:, l], prim[:, l])
    end

    return nothing

end

ext: jl


text:
@inline context(ctx::Context) = ctx
@inline context(conn::Connection) = conn.context
@inline context(stmt::Stmt) = context(stmt.connection)
@inline context(pool::Pool) = pool.context
@inline context(variable::Variable) = context(variable.connection)
@inline context(lob::Lob) = context(lob.parent)
@inline context(v::JuliaOracleValue) = v.context
@inline context(v::ExternOracleValue) = context(v.parent)
@inline context(queue::Queue) = context(queue.connection)
@inline context(obj_type::ObjectType) = context(obj_type.connection)
@inline context(msg_props::Message) = context(msg_props.parent)

function Context()
    error_info_ref = Ref{OraErrorInfo}()
    context_handle_ref = Ref{Ptr{Cvoid}}()
    result = dpiContext_createWithParams(ORA_MAJOR_VERSION, ORA_MINOR_VERSION, context_handle_ref, error_info_ref)

    if result == ORA_FAILURE
        error_info = error_info_ref[]

        if error_info.code == 0
            error_message = unsafe_string(error_info.message, error_info.message_length)

            if startswith(error_message, "DPI-1020")
                error("$error_message. Please, run `Pkg.build(\"Oracle\")` and restart Julia.")
            else
                throw(error_info)
            end
        end

        throw(error_info)
    end
    @assert result == ORA_SUCCESS

    return Context(context_handle_ref[])
end

function client_version(ctx::Context) :: OraVersionInfo
    version_info_ref = Ref{OraVersionInfo}()
    result = dpiContext_getClientVersion(ctx.handle, version_info_ref)
    error_check(ctx, result)
    return version_info_ref[]
end


text:
@testset "ClassifiedCollections" begin

    # classified lists
    @testset "ClassifiedLists" begin
        c = classified_lists(String, Int)

        push!(c, "low", 1)
        push!(c, "low", 2)
        push!(c, "low", 3)
        push!(c, "high", 4)
        push!(c, "high", 5)

        @test haskey(c, "low")
        @test haskey(c, "high")
        @test !haskey(c, "mid")

        @test c["low"] == [1, 2, 3]
        @test c["high"] == [4, 5]

        @test length(c) == 2
        @test sort(collect(keys(c))) == String["high","low"]

        pop!(c,"low")
        @test !haskey(c,"low")

    end # @testset ClassifiedLists

    # classified sets
    @testset "ClassifiedSets" begin
        c = classified_sets(String, Int)

        push!(c, "low", 1)
        push!(c, "low", 2)
        push!(c, "low", 3)
        push!(c, "low", 1)
        push!(c, "low", 2)

        push!(c, "high", 4)
        push!(c, "high", 5)
        push!(c, "high", 5)

        @test haskey(c, "low")
        @test haskey(c, "high")
        @test !haskey(c, "mid")

        @test isa(c["low"], Set{Int})
        @test isa(c["high"], Set{Int})

        @test sort(collect(c["low"])) == [1, 2, 3]
        @test sort(collect(c["high"])) == [4, 5]
    end # @testset ClassifiedSets

    # classified counters
    @testset "ClassifiedCounters" begin

        c = classified_counters(String, Float64)

        push!(c, "low", 1.)
        push!(c, "low", 2.)
        push!(c, "low", 3.)
        push!(c, "low", 1.)
        push!(c, "low", 2.)
        push!(c, "low", 2.)

        push!(c, "high", 4.)
        push!(c, "high", 5.)
        push!(c, "high", 5.)

        @test haskey(c, "low")
        @test haskey(c, "high")
        @test !haskey(c, "mid")

        cl = c["low"]
        ch = c["high"]

        @test isa(cl, Accumulator{Float64, Int})
        @test isa(ch, Accumulator{Float64, Int})

        @test cl[1.] == 2
        @test cl[2.] == 3
        @test cl[3.] == 1
        @test ch[4.] == 1
        @test ch[5.] == 2
    end # @testset ClassifiedCounters

end # @testset ClassifiedCollections


text:
## Code for our Types and conversion methods

## Symbol class for controlling dispatch
abstract type SymbolicObject <: Number end

## Basic types defined here
"""

* `Sym("x")`
* `Sym(:x)`
* `Sym("x, y, z")`

The `Sym` type is an immutable type holding a reference to an
underlying python-based SymPy object. Many methods are extended to the
`Sym` type. Instances can be constructed in many ways. The one caveat
is the variables can not be function names in base.

"""
struct Sym <: SymbolicObject
    x::PyCall.PyObject
end
Sym(s::SymbolicObject) = s

## sets
struct SymSet <: SymbolicObject
    x::PyCall.PyObject
end

## complex float
## this cause issue with printing on non-complex objects
#mpctype = sympy.mpmath["ctx_mp_python"]
#pytype_mapping(mpctype, Sym)

## some typealiases
const SymOrReal =  Union{Sym,Real}
const SymOrNumber =  Union{Sym,Number}
const SymOrString =  Union{Sym,AbstractString}
const SymbolicTypes = Union{AbstractString,Symbol,SymbolicObject}

## in #83, @stevengj suggests using
PyCall.PyObject(x::SymbolicObject) = x.x

## Promotion
## promote up to symbolic so that math ops work
promote_rule(::Type{T}, ::Type{S})  where {T<:SymbolicObject, S<:Number}= T

## Conversion
convert(::Type{T}, o::PyCall.PyObject) where {T <: SymbolicObject} = T(o)
convert(::Type{PyObject}, s::Sym) = s.x

function convert(::Type{Tuple}, o::PyCall.PyObject)
    ## check that o is a tuple?
    ## PyCall.pytypeof(o)
    n = o[:__len__]()
    ntuple(i -> o[:__getitem__](i-1), n)
end

## rational
convert(::Type{T}, x::Rational) where {T<:SymbolicObject} = sympy_meth(:Rational, x.num, x.den)::T

## big. Need mpmath installed separately -- not as a SymPy module as that is how it is called in PyCall
convert(::Type{T}, x::BigFloat) where {T<:SymbolicObject} = Sym(PyCall.PyObject(x))::T
convert(::Type{Sym}, x::Complex{BigFloat}) = Sym(PyCall.PyObject(x))::Sym

## real
convert(::Type{S}, x::T) where {S<:SymbolicObject, T <: Real}= sympy_meth(:sympify, x)::S
convert(::Type{T}, x::Sym) where {T <: Real} = convert(T, PyObject(x))

## complex
## IM is SymPy's "i" (sympy["I"], not Python's
## Sym(PyCall.PyObject(im)) which gives 1j.
function convert(::Type{Sym}, x::Complex)
    y = ifelse(isa(x, Complex{Bool}), real(x) + imag(x) * im, x)
    real(y) + imag(y) * IM
end
convert(::Type{Complex}, x::Sym) = complex(map(x -> convert(Float64, x), x[:as_real_imag]())...)::Sym
complex(x::Sym) = convert(Complex, x)
complex(xs::AbstractArray{Sym}) = map(complex, xs)

## string
convert(::Type{Sym}, o::AbstractString) = sympy_meth(:sympify, o)
convert(::Type{Sym}, o::Symbol) = sympy_meth(:sympify, string(o))

## function
convert(::Type{Function}, ex::Sym) = lambdify(ex)

## we usually promote to Sym objects, but here we want to promote to functions
## so [x, sin] -> will be plottable as two functions
Base.promote_rule(::Type{T}, ::Type{S}) where {T<:SymbolicObject, S<:Function} = S


text:
using DataFrames, LinearAlgebra, MixedModels, Random, RData, SparseArrays, Test

if !@isdefined(dat) || !isa(dat, Dict{Symbol, DataFrame})
    dat = Dict(Symbol(k) => v for (k, v) in load(joinpath(dirname(@__FILE__), "dat.rda")))
end

@testset "scalarRe" begin
    dyestuff = dat[:Dyestuff]
    pastes = dat[:Pastes]
    batch = dyestuff[:G]
    sf = ScalarFactorReTerm(batch, :G)
    b2 = categorical(reverse(string.(batch)))
    reverse!(b2.refs)
    bb = ScalarFactorReTerm(b2, :G)
    sf1 = ScalarFactorReTerm(pastes[:G], :G)
    sf2 = ScalarFactorReTerm(pastes[:H], :H)
    Yield = dyestuff[:Y]
    ScalarFactorReTerm(repeat(1:3, inner=2), :G) # just to test a constructor

    @testset "size" begin
        @test size(sf) == (30, 6)
        @test size(sf,1) == 30
        @test size(sf,2) == 6
        @test size(sf,3) == 1
        @test size(sf1) == (60, 30)
        @test size(sf2) == (60, 10)
    end

    @testset "order" begin
        @test bb.refs == sf.refs
    end

    @testset "utilities" begin
        @test MixedModels.levs(sf) == string.('A':'F')
        @test MixedModels.nlevs(sf) == 6
        @test MixedModels.vsize(sf) == 1
        @test MixedModels.nrandomeff(sf) == 6
        @test eltype(sf) == Float64
        @test sparse(sf) == sparse(Int32[1:30;], convert(Vector{Int32},sf.refs), ones(30))
        fsf = Matrix(sf)
        @test size(fsf) == (30, 6)
        @test count(!iszero, fsf) == 30
        @test sort!(unique(fsf)) == [0.0, 1.0]
        @test cond(sf) == 1.0
        @test MixedModels.nθ(sf) == 1
        @test getθ(sf) == ones(1)
        @test MixedModels.getθ!(Vector{Float64}(undef, 1), sf) == ones(1)
        @test lowerbd(sf) == zeros(1)
        @test eltype(sf) == Float64
        @test getθ(setθ!(sf, 0.5)) == [0.5]
        @test_throws DimensionMismatch MixedModels.getθ!(Float64[], sf)
        @test_throws MethodError setθ!(sf, ones(2))
    end

    @testset "products" begin
        @test MatrixTerm(ones(30))'sf == fill(5.0, (1, 6))
        @test mul!(Array{Float64}(undef, (size(sf1, 2), size(sf2, 2))), sf1', sf2) == Array(sf1'sf2)

        crp = sf'sf
        @test isa(crp, Diagonal{Float64})
        crp1 = copy(crp)
        @test crp1 == crp
        @test crp[2,6] == 0
        @test crp[6,6] == 5
        @test size(crp) == (6,6)
        @test crp.diag == fill(5.,6)
        rhs = MatrixTerm(Yield)'sf
        @test rhs == reshape([7525.0,7640.0,7820.0,7490.0,8000.0,7350.0], (1, 6))
        @test ldiv!(crp, copy(rhs)') == reshape([1505.,1528.,1564.,1498.,1600.,1470.], (6, 1))

        @test isa(sf1'sf1, Diagonal{Float64})
        @test isa(sf2'sf2, Diagonal{Float64})
        @test isa(sf2'sf1,SparseMatrixCSC{Float64})

        @test lmul!(Λ(sf)', ones(6)) == fill(0.5, 6)
        @test rmul!(ones(6, 6), Λ(sf)) == fill(0.5, (6, 6))

        @test mul!(Matrix{Float64}(undef, 1, 6), Λ(sf), ones(1, 6)) == fill(0.5, (1,6))
    end

    @testset "reweight!" begin
        wts = rand(MersenneTwister(1234321), size(sf, 1))
        @test vec(MixedModels.reweight!(sf, wts).wtz) == wts
    end
end

@testset "vectorRe" begin
    slp = dat[:sleepstudy]
    corr = VectorFactorReTerm(slp[:G], vcat(ones(1, size(slp, 1)), slp[:U]'),
        :G, ["(Intercept)", "U"], [2])
    nocorr = VectorFactorReTerm(slp[:G], vcat(ones(1, size(slp, 1)), slp[:U]'),
            :G, ["(Intercept)", "U"], [1, 1])
    Reaction = slp[:Y]

    @testset "sizes" begin
        @test size(corr) == (180,36)
        @test size(nocorr) == (180,36)
    end

    @testset "utilities" begin
        @test MixedModels.levs(corr) == levels(slp[:G])
        @test MixedModels.nlevs(corr) == 18
        @test MixedModels.vsize(corr) == 2
        @test MixedModels.nrandomeff(corr) == 36
        @test eltype(corr) == Float64
        @test nnz(sparse(corr)) == 360
        @test cond(corr) == 1.0
        @test MixedModels.nθ(corr) == 3
        @test MixedModels.nθ(nocorr) == 2
        @test getθ(corr) == [1.0, 0.0, 1.0]
        @test getθ(nocorr) == ones(2)
        @test MixedModels.getθ!(Vector{Float64}(undef, 2), nocorr) == ones(2)
        @test lowerbd(nocorr) == zeros(2)
        @test lowerbd(corr) == [0.0, -Inf, 0.0]
        @test getθ(setθ!(corr, fill(0.5, 3))) == [0.5, 0.5, 0.5]
        @test_throws DimensionMismatch MixedModels.getθ!(Vector{Float64}(undef, 2), corr)
        @test_throws DimensionMismatch setθ!(corr, ones(2))
    end

    vrp = corr'corr

    @test isa(vrp, UniformBlockDiagonal{Float64})
    @test size(vrp) == (36, 36)

    @test mul!(Array{Float64}(undef, 36, 36), corr', corr) == Matrix(vrp)
    scl = ScalarFactorReTerm(slp[:G].refs, levels(slp[:G]), Array(slp[:U]), Array(slp[:U]), :G, ["U"], 1.0)

    @test sparse(corr)'sparse(scl) == corr'scl

    b = mul!(Matrix{Float64}(undef, 2,18), Λ(corr).data, ones(2,18))
    @test b == reshape(repeat([0.5, 1.0], outer=18), (2,18))

    @testset "reweight!" begin
        wts = rand(MersenneTwister(1234321), size(corr, 1))
        @test MixedModels.reweight!(corr, wts).wtz[1, :] == wts
        @test corr.z[1, :] == ones(size(corr, 1))
    end

end


text:
## cauchy

function cauchycircleS(cfs::AbstractVector,z::Number,s::Bool)
    ret=zero(Complex{Float64})

    if s
        zm = one(Complex{Float64})

        #odd coefficients are pos
        @simd for k=1:2:length(cfs)
            @inbounds ret += cfs[k]*zm
            zm *= z
        end
    else
        z=1/z
        zm = z

        #even coefficients are neg
        @simd for k=2:2:length(cfs)
            @inbounds ret -= cfs[k]*zm
            zm *= z
        end
    end

    ret
end

function stieltjes(sp::Laurent{DD,RR},f::AbstractVector,z::Number) where {DD<:Circle,RR}
    d=domain(sp)
    if !d.orientation
        return -stieltjes(reverseorientation(Fun(sp,f)),z)
    end

    z=mappoint(d,Circle(),z)
    -2π*im*cauchycircleS(f,z,abs(z) < 1)
end

function stieltjes(sp::Laurent{DD,RR},f::AbstractVector,z::Directed{s}) where {DD<:Circle,RR,s}
    d=domain(sp)
    if !d.orientation
        return -stieltjes(reverseorientation(Fun(sp,f)),reverseorientation(z))
    end

    z=mappoint(d,Circle(),z)
    -2π*im*cauchycircleS(f,undirected(z),orientation(z))
end

stieltjes(sp::Fourier{DD,RR}, f::AbstractVector, z::Number) where {DD<:Circle,RR} =
    stieltjes(Laurent(domain(sp)),coefficients(f,sp,Laurent(domain(sp))),z)

# we implement cauchy ±1 as canonical
# TODO: reimplement directly
hilbert(sp::Laurent{DD,RR}, f::AbstractVector, z::Number) where {DD<:Circle,RR} =
    (stieltjes(sp,f,(z)⁺)+stieltjes(sp,f,(z)⁻))/(-2π)

## stieltjesintegral and logkernel

function stieltjesintegral(sp::Laurent{DD,RR},f::AbstractVector,z::Number) where {DD<:Circle,RR}
    d=domain(sp)
    @assert d==Circle()  #TODO: radius
    ζ=Fun(d)
    r=stieltjes(integrate(f-f[2]/ζ),z)
    abs(z)<1 ? r : r+2π*im*f[2]*log(z)
end

stieltjesintegral(sp::Fourier{DD,RR},f::AbstractVector,z::Number) where {DD<:Circle,RR}=stieltjesintegral(Fun(Fun(sp,f),Laurent),z)

function logkernel(sp::Fourier{DD,RR},g::AbstractVector,z::Number) where {DD<:Circle,RR}
    d=domain(sp)
    c,r=d.center,d.radius
    z=z-c
    if abs(z) ≤r
        ret=2r*log(r)*g[1]
        for j=2:2:length(g)
            k=div(j,2)
            ret+=-g[j]*sin(k*angle(z))*abs(z)^k/(k*r^(k-1))
        end
        for j=3:2:length(g)
            k=div(j,2)
            ret+=-g[j]*cos(k*angle(z))*abs(z)^k/(k*r^(k-1))
        end
        ret
    else
        ret=2r*logabs(z)*g[1]
        for j=2:2:length(g)
            k=div(j,2)
            ret+=-g[j]*sin(k*angle(z))*r^(k+1)/(k*abs(z)^k)
        end
        for j=3:2:length(g)
            k=div(j,2)
            ret+=-g[j]*cos(k*angle(z))*r^(k+1)/(k*abs(z)^k)
        end
        ret
    end
end

logkernel(sp::Laurent{DD,RR},g::AbstractVector,z) where {DD<:Circle,RR} = logkernel(Fun(Fun(sp,g),Fourier),z)


text:
using Rubin
using Tests
using Elliptic
using HypergeometricFunctions
using Polylogarithms
using SpecialFunctions

@test integrate((a+a*cos(x))^-1*tan(x)^4, x) == :((1/2)*a^-1*arctanh(sin(x))+1//3*a^-1*tan(x)^3+-1//2*a^-1*sec(x)*tan(x))
@test integrate((a+a*cos(x))^-1*tan(x)^3, x) == :((1/2)*a^-1*sec(x)^2+-1*a^-1*sec(x))
@test integrate((a+a*cos(x))^-1*tan(x)^2, x) == :(a^-1*tan(x)+-1*a^-1*arctanh(sin(x)))
@test integrate((a+a*cos(x))^-1*tan(x), x) == :(a^-1*log(1+cos(x))+-1*a^-1*log(cos(x)))
@test integrate((a+a*cos(x))^-1*cot(x), x) == :(-1//2*a^-1*csc(x)^2+-1//2*a^-1*arctanh(cos(x))+(1/2)*a^-1*cot(x)*csc(x))
@test integrate((a+a*cos(x))^-1*cot(x)^2, x) == :(-1*a^-1*csc(x)+-1//3*a^-1*cot(x)^3+1//3*a^-1*csc(x)^3)
@test integrate((a+a*cos(x))^-1*cot(x)^3, x) == :(-1//4*a^-1*cot(x)^4+3//8*a^-1*arctanh(cos(x))+-3//8*a^-1*cot(x)*csc(x)+1//4*a^-1*cot(x)^3*csc(x))
@test integrate((a+a*cos(x))^-1*cot(x)^4, x) == :(a^-1*csc(x)+-2//3*a^-1*csc(x)^3+-1//5*a^-1*cot(x)^5+1//5*a^-1*csc(x)^5)
@test integrate((1+cos(3x))^-2*tan(3x), x) == :(-1*(3+3*cos(3x))^-1+-1//3*log(cos(3x))+1//3*log(1+cos(3x)))
@test integrate((a+b*cos(x))^-1*tan(x)^4, x) == :(-1//3*a^-3*(-3*b^2+4*a^2)*tan(x)+1//3*a^-1*sec(x)^2*tan(x)+(1/2)*b*a^-4*(-2*b^2+3*a^2)*arctanh(sin(x))+2*a^-4*(a+b)^3//2*(a+-1b)^3//2*arctan((a+b)^-1//2*(a+-1b)^(1/2)*tan((1/2)*x))+-1//2*b*a^-2*sec(x)*tan(x))
@test integrate((a+b*cos(x))^-1*tan(x)^3, x) == :((1/2)*a^-1*sec(x)^2+a^-3*(a^2+-1*b^2)*log(cos(x))+-1*b*a^-2*sec(x)+-1*a^-3*(a^2+-1*b^2)*log(a+b*cos(x)))
@test integrate((a+b*cos(x))^-1*tan(x)^2, x) == :(a^-1*tan(x)+-1*b*a^-2*arctanh(sin(x))+-2*a^-2*(a+b)^(1/2)*(a+-1b)^(1/2)*arctan((a+b)^-1//2*(a+-1b)^(1/2)*tan((1/2)*x)))
@test integrate((a+b*cos(x))^-1*tan(x), x) == :(a^-1*log(a+b*cos(x))+-1*a^-1*log(cos(x)))
@test integrate((a+b*cos(x))^-1*cot(x), x) == :((-2b+2a)^-1*log(1+cos(x))+(2a+2b)^-1*log(1+-1*cos(x))+-1*a*(a^2+-1*b^2)^-1*log(a+b*cos(x)))
@test integrate((a+b*cos(x))^-1*cot(x)^2, x) == :(b*(a^2+-1*b^2)^-1*csc(x)+-1*a*(a^2+-1*b^2)^-1*cot(x)+-2*a^2*(a+b)^-3//2*(a+-1b)^-3//2*arctan((a+b)^-1//2*(a+-1b)^(1/2)*tan((1/2)*x)))
@test integrate((a+b*cos(x))^-1*cot(x)^3, x) == :(a^3*(a^2+-1*b^2)^-2*log(a+b*cos(x))+-1*(-2*b^2+2*a^2)^-1*csc(x)^2*(a+-1*b*cos(x))+-1//4*(a+b)^-2*(b+2a)*log(1+-1*cos(x))+-1//4*(a+-1b)^-2*(-1b+2a)*log(1+cos(x)))
@test integrate((a+b*cos(x))^-1*cot(x)^4, x) == :(b*(-3*b^2+3*a^2)^-1*csc(x)^3+a^3*(a^2+-1*b^2)^-2*cot(x)+-1*a*(-3*b^2+3*a^2)^-1*cot(x)^3+-1*b*(a^2+-1*b^2)^-1*csc(x)+-1*b*a^2*(a^2+-1*b^2)^-2*csc(x)+2*a^4*(a+b)^-5//2*(a+-1b)^-5//2*arctan((a+b)^-1//2*(a+-1b)^(1/2)*tan((1/2)*x)))
@test integrate((3+-1*cos(x))^-1//2*cot(x), x) == :(-1//2*arctanh((1/2)*(3+-1*cos(x))^(1/2))+-1//2*2^(1/2)*arctanh((1/2)*2^(1/2)*(3+-1*cos(x))^(1/2)))
@test integrate((a+b*cos(x))^(1/2)*tan(x), x) == :(-2*(a+b*cos(x))^(1/2)+2*a^(1/2)*arctanh(a^-1//2*(a+b*cos(x))^(1/2)))
@test integrate((a+b*cos(x))^-1//2*tan(x), x) == :(2*a^-1//2*arctanh(a^-1//2*(a+b*cos(x))^(1/2)))

max_stars_repo_path: test/4 Trig functions/4.2 Cosine/4.2.1.3 (g tan)^p (a+b cos)^m.jl


text:
"""
    KRatioOptimizer

Defines an optimizeks(kro::KRatioOptimizer, krs::Vector{KRatio})::Vector{KRatio} method which takes a vector
of k-ratios which may have redundant data (more than one KRatio per element) and trims it down to a vector
of k-ratios with one KRatio per element.
"""
abstract type KRatioOptimizer end

"""
    SimpleKRatioOptimizer(overvoltage, favor::Vector{CharXRay} = CharXRay[], minOvervoltage=1.2)

Implements a simple optimizer based on shell first, overvoltage next and brightness last.  Once it picks
an optimum set of lines for an element, it will not change.  You can force the selection of specific
lines by including the brightest in a family in `favor`.  Lines with overvoltages less than `minOvervoltage`
are not considered as available.
"""
struct SimpleKRatioOptimizer <: KRatioOptimizer
    overvoltage::Float64
    scores::Dict{Vector{CharXRay},AbstractFloat}
    favor::Vector{CharXRay}
    minOver::Float64

    function SimpleKRatioOptimizer(overvoltage, favor::Vector{CharXRay} = CharXRay[], minOvervoltage=1.2)
        elms = map(element, favor)
        @assert length(unique(elms)) == length(elms) "SimpleKRatioOptimizer: Please specify only one characteristic X-ray per element to favor."
        new(overvoltage, Dict{Vector{CharXRay}, AbstractFloat}(), favor, minOvervoltage)
    end
end

Base.show(io::IO, skro::SimpleKRatioOptimizer) = print(io,"SimpleKRatioOptimizer[U>$(skro.overvoltage), favor=[$(skro.favor)]]")

function optimizeks(skro::SimpleKRatioOptimizer, krs::AbstractVector{T})::Vector{T} where  T <: Union{KRatio, KRatios}
    function score(kr) # Larger is better....
        br = brightest(kr.xrays)
        sc = any(f->f in kr.xrays, skro.favor) && (kr.unkProps[:BeamEnergy] / edgeenergy(br) > skro.minOver) ? #
            1.0e100 : get(skro.scores, kr.xrays, -1.0)
        if sc == -1.0
            ov = min(kr.stdProps[:BeamEnergy], kr.unkProps[:BeamEnergy]) / edgeenergy(br)
            sc = if ov > skro.minOver
                5.0 - n(shell(br)) - # Line K->8, L->6, M->4, N->2
                    skro.overvoltage / ov + # Overvoltage (<1 if ov > over)
                    0.1*sum(weight.(NormalizeToUnity, kr.xrays)) # line weight (favor brighter)
            else
                0.0 # Avoid it...
            end
            skro.scores[kr.xrays]=sc
        end
        return sc
    end
    return map(collect(elms(krs))) do elm
        elmkrs=filter(k->k.element==elm, krs)
        maxi = findmax(score.(elmkrs))
        maxi[1] == 0.0 && @warn "There is no k-ratio with an overvoltage greater than $(skro.minOver) for $elm."
        elmkrs[maxi[2]]
    end
end

max_stars_repo_name: NicholasWMRitchie/NeXLMatrixCorrection

max_issues_repo_name: NicholasWMRitchie/NeXLMatrixCorrection

max_forks_repo_name: NicholasWMRitchie/NeXLMatrixCorrection


text:
"""
    βs(X::T, Y::T) where {T<:BinaryNetwork}

Components of β-diversity as measured on species.
"""
function βs(X::T, Y::T) where {T<:BinaryNetwork}
    a = richness(intersect(X,Y))
    b = richness(Y)-a
    c = richness(X)-a
    return (a=a, b=b, c=c)
end

"""
Overlapping species (bipartite)
"""
function βos(X::T, Y::T) where {T<:BipartiteNetwork}
    core = intersect(X,Y)
    Xi = X[species(core; dims=1),species(core; dims=2)]
    Yi = Y[species(core; dims=1),species(core; dims=2)]
    a = links(intersect(Xi,Yi))
    b = links(Xi)-a
    c = links(Yi)-a
    return (a=a, b=b, c=c)
end

"""
Overlapping species (unipartite)
"""
function βos(X::T, Y::T) where {T<:UnipartiteNetwork}
    core = intersect(X,Y)
    Xi = X[species(core)]
    Yi = Y[species(core)]
    a = links(intersect(Xi,Yi))
    b = links(Xi)-a
    c = links(Yi)-a
    return (a=a, b=b, c=c)
end

"""
Whole network
"""
function βwn(X::T, Y::T) where {T<:BinaryNetwork}
    a = links(intersect(X,Y))
    c = links(X)-a
    b = links(Y)-a
    return (a=a, b=b, c=c)
end


text:
"""Julia program to implement Interpolation Search algorithm.
It is simliar to a person searching for a name in a telephone directory
The algorithm calculates where in the remaining array space the required item would be,
based on the values at the bounds of the search space and the required element, usually via a linear interpolation
"""

function interpolation_search(arr, n, ele)
    low = 1
    high = n
    while low <= high && ele >= arr[low] && ele <= arr[high]
        if(low == high)
            # If element found
            if(arr[low] == ele)
                return true
            end
            return false
        end
        # Formula to return the value of 'pos'
        pos = low + (((high - low) ÷ (arr[high] - arr[low]))* (ele - arr[low]))
        # If the pos element is the required element, return that index
        if(arr[pos] == ele)
            return true
        end
        # Else find the section of the array in which the required element belongs to
        if(arr[pos] < ele)
            low = pos + 1
        else
            high = pos - 1
        end
    end
    return false
end

print("How many numbers are present in the array? ")
n = readline()
n = parse(Int, n)
if (n <= 0)
    println("Array is Empty!!!")
    exit()
end
arr = Int[]
print("Enter the numbers: ")
arr = [parse(Int, num) for num in split(readline())]
print("Which number do you want to search in the array? ")
ele = readline()
ele = parse(Int, ele)
# Sort the array in ascending order
arr = sort(arr)
res = interpolation_search(arr, n, ele)
if (res == 0)
    print("The number $ele is not present in the array")
else
    print("The number $ele is present in the array.")
end

"""
Time Complexity - O(log(n)), where 'n' is the size of the array
Space Complexity - O(n)

SAMPLE INPUT AND OUTPUT

SAMPLE I

How many numbers are present in the array? 5
Enter the numbers: 1 2 3 4 5
Which number do you want to search in the array? 6
The number 6 is not present in the array

SAMPLE II

How many numbers are present in the array? 3
Enter the numbers: 3 1 2
Which number do you want to search in the array? 2
The number 2 is present in the array.

"""

max_stars_repo_name: TechSpiritSS/NeoAlgo


text:
module antfarm
using Revise
using Debugger

## internal

include("multiagent/simulation.jl")
include("multiagent/simulation_init.jl")
include("multiagent/simulation_utils.jl")
include("fmp/fmp.jl")
include("analysis/gridlock_analysis.jl")
include("analysis/dispersal_simulation.jl")
include("analysis/simulation_wrapper.jl")

end # module

max_stars_repo_name: rmcsqrd/antfarm

max_issues_repo_name: rmcsqrd/antfarm

max_forks_repo_name: rmcsqrd/antfarm


text:
"""
Generic function to estimate decay time on a single vector

Args:
	PP (scipy matrix): approximation of the discretized matrix to use
	v (numpy vector):
	my_norm (function):
	M, K (reals with RNDU): constants such that ||P-PP|| ≤ K` and ||P^n|| ≤ M`,
	where P is the exact discretized matrix.

Returns:
	n such that ||PP^n v|| ≤ target`.

Raises:
	ValueError if insufficient precision is detected
"""
function vector_decay_time(B::Basis, PP, v, M=0, K=0, target = 0.5)
	current_norm = weak_norm(B, v)
	error_on_computed_norm = zero(M)
	n = 0
	MK = M * K
	error_propagation_constant = K #this constant is K at the first step and MK afterwards; see notes
	while current_norm + error_on_computed_norm >= target
		n += 1
		v = PP * v
		current_norm = weak_norm(B, v)
		error_on_computed_norm += error_propagation_constant * current_norm
		error_propagation_constant = MK
		if error_on_computed_norm > target
			raise() # ValueError, 'Insufficient precision'
		end
	end
	return n
end

"""
Number of iterations needed to contract all vectors in `basis.contracting_pairs()` to a given target alpha
"""
function decay_time(D::Dynamic, B::Basis, P::AbstractMatrix{Interval{T}}, alpha = 0.5, n_jobs = 1) where {T}

	PP = mid.(P)

	M = bound_on_norms_of_powers(B, D, project_left=True, project_right=True)
	K = numerical_error(B, D, P, PP)
	my_norm(v) = weak_norm(B, v)

	#alpha = Rdown(alpha) # we need a lower bound to alpha*s to make sure we contract below it
	#decay_times = Parallel(n_jobs=n_jobs, verbose=1)(delayed(vector_decay_time)(PP, v, my_norm, M, K, alpha*s) for v, s in basis.contracting_pairs())

	return max(decay_times)
end


text:
# ---
# layout: post
# title: "Random noise"
# date: 2018-10-16 00:00:00 +0000
# categories: blog
# mathjax: true
# ---

# >This is another post in the ongoing series in which I try to learn 2D vector graphics using Julia. It doesn't contain any revelations or new material, and you should visit the following sites if you're looking for a good introduction to the subject of noise in graphics:

# - [red blob games](https://www.redblobgames.com/articles/noise/introduction.html){:target="_blank"}
# - [adrianbiagioli](http://flafla2.github.io/2014/08/09/perlinnoise.html){:target="_blank"}
# - [khan academy](https://www.khanacademy.org/partner-content/pixar/pattern/perlin-noise/v/patterns-9){:target="_blank"}

# I'm using Julia version 1.0 if you want to play along; you can find the source files and notebooks on the github. If you do, you'll need the packages [Luxor](https://github.com/JuliaGraphics/Luxor.jl){:target="_blank"}, [Colors](https://github.com/JuliaGraphics/Colors.jl){:target="_blank"}, and [ColorSchemes](https://github.com/JuliaGraphics/ColorSchemes.jl){:target="_blank"}. I used  [Literate.jl](https://github.com/fredrikekre/Literate.jl){:target="_blank"} to produce the Markdown and Jupyter notebook versions.

# ### Random versus Noise

# Luxor provides a function called `noise()`. This can accept a single floating-point number as input, and it returns a value between 0.0 and 1.0.

using Luxor

noise(0.0)
0.8261106995884773

noise(1.0)
0.5

noise(2.0)
0.4460609053497945

# It will be easier to draw some graphs. Here's a quick throwaway function to draw a simple graph.

function graph(a, width = 800;
        startnumber     = 0,
        endnumber       = 1,
        style           = :line,
        margin          = 30)
    setline(1)
    bars(a, labels      =false,
            xwidth      = (width - 2margin)/length(a),
            yheight     =40,
            barfunction = (bottom::Point, top::Point, value;
        extremes=extrema(a), barnumber=0, bartotal=0) ->
            begin
                if style == :line
                    line(bottom, top, :stroke)
                else
                    circle(top, 1, :fill)
                end
            end)
    sethue("black")
    label(string(startnumber), :S,
        Point(0, 0), offset=10)
    label(string(endnumber), :S,
        Point(width - 2margin, 0), offset=10)
end

function drawgraph(startvalue, endvalue, filename)
    Drawing(800, 150, filename)
    background("white")
    origin()
    ## move to top left corner
    margin=30
    translate(BoundingBox()[1] + (margin, boxheight(BoundingBox()/2)))
    sethue("black")
    graph(noise.(range(startvalue, endvalue, length=200)),
        startnumber=startvalue,
        endnumber=endvalue)
    finish(); preview()
end

# To test this out, graph 200 random integers:

Drawing(800, 150, "images/noise/graph-random.png")
background("white")
origin() ## move to top left corner
margin=30
translate(BoundingBox()[1] + (margin, boxheight(BoundingBox()/2)))
sethue("black")
graph(rand(200))
finish(); preview()

# ![image label](IMAGEFOLDER/graph-random.png)

# To start with, let's graph the output of the `noise()` function for the first 200 integers:

drawgraph(0, 200, "images/noise/graph-0-200.png")

# ![image label](IMAGEFOLDER/graph-0-200.png)

# It looks very random. But let's look at 200 values between 0 and 10:

drawgraph(0, 10, "images/noise/graph-0-10.png")

#- Literate needs this to trigger changes of context

# ![image label](IMAGEFOLDER/graph-0-10.png){: .center-image}

# There's some randomness, but it's smoother, and looks more natural.

# Zoom and enhance, between 0 and 5:

drawgraph(0, 5, "images/noise/graph-5.png")

#-

# ![image label](IMAGEFOLDER/graph-5.png){: .center-image}

# You can see that the left half of the 0 to 10 graph has been stretched.

# Between 0 and 1:

drawgraph(0, 1, "images/noise/graph-0-1.png")

#-

# ![image label](IMAGEFOLDER/graph-0-1.png){: .center-image}

# One more for luck:

drawgraph(0, 0.5, "images/noise/graph-0-05.png")

#-

# ![image label](IMAGEFOLDER/graph-0-05.png){: .center-image}

# The more often you sample the noise space (ie the shorter the gaps between the set of values passed to `noise()`, the closer together the output values will be.

# So the `noise()` function provides gently changing undulations rather than the unpredictable jumps of randomness. Here's a slowly changing color pattern in the LCHab color space, using a set of noisy values to choose the hue. We'll change the length of each line as well just for fun:

using Colors
@svg begin  ## switch to SVG for better graphic quality
    rate = .009
    setline(2)
    for x in -250:2:250
        yn = noise(x * rate)
        hue = rescale(yn, 0, 1, 0, 359)
        sethue(LCHab(50, 100, hue))
        ht = rescale(yn, 0, 1, 10, 100)
        line(Point(x, -ht/2), Point(x, ht/2), :stroke)
    end
end 600 200 "images/noise/colorbars.svg"

#-

# ![image label](IMAGEFOLDER/colorbars.svg){: .center-image}

# You can use noisy values to specify other changing parameters. For example, let's place some pebbles at random, and control their size using a noisy distribution, to give the illusion of a naturally changing distribution.

function drawpebble(pt, radius)
    sethue("grey60")
    @layer begin
        transform([rand(0.5:0.1:1) 0 0 rand(0.5:0.1:1) 0 0])
        circle(pt, radius, :fill)
        for i in 1:-0.02:0.2
            sethue(rescale(i, 1, 0,   0.5 + rand(0:0.1:0.3), 1.0),
                   rescale(i, 1, 0.1, 0.4 + rand(0:0.1:0.5), 1.0),
                   rescale(i, 1, 0.3, 0.5, 1.0))
            setopacity(1 - i)
            circle(pt + (-2i, -2i), i * radius, :fill)
        end
    end
end

@png begin
    ## switch to PNG, SVG can't handle this
    background("palegoldenrod")
    pebblesize = 12
    for i in 1:6000
        pt = Point(rand(-400:400), rand(-200:200))
        n = noise(pt.x * 0.002)
        drawpebble(pt, pebblesize * n)
    end
end 800 400 "images/noise/pebbles.png"

#-

# ![image label](IMAGEFOLDER/pebbles.png){: .center-image}

# ### Detail and persistence

# The `noise()` function has two optional keyword arguments that let you tweak the knobs of the noise generator.

# The first is `detail`, an integer. Increasing it from the default value of 1 upwards will add finer detail to the basic noise. The second is `persistence`, a floating-point value between 0 and 1 (or more).

# 'detail` is graphed here with values from 1 to 12. As the level increases, you can see that the same overall noise contours are gradually modulated with finer variations.

function detailgraph()
    @svg begin
        margin=30
        translate(BoundingBox()[1] + (margin, 0))
        setline(.5)
        sethue("black")
        stopat = 2
        r = range(0, length=400, stop=stopat)
        for detail in 1:2:12
            translate(0, 100)
            sethue("red")
            graph(noise.(r), style=:circle, endnumber=stopat)
            sethue("black")
            text("detail = $detail, persistence = 0.9", Point(200, 15))
            graph(noise.(r, detail=detail, persistence=0.9), endnumber=stopat)
        end
    end 800 650 "images/noise/detail-graph.svg"
 end
detailgraph()

#-

# ![image label](IMAGEFOLDER/detail-graph.svg){: .center-image}

# You can see the original noisy curve (in red) behind each more detailed graph. The noise generator is doubling the frequency but halving the amplitude every time you go one level higher. Noise, like music, can have octaves of higher frequencies mixed with lower fundamental frequency. The `detail` keyword is adding one or more octaves of noise.

# The `persistence` argument defaults to zero. The value controls the amplitude of each successive octave of noise, with higher values of persistence producing higher levels of finer detail, as the values persist for longer.

function persistencegraph()
    @svg begin
        setline(.5)
        sethue("black")
        margin=30
        translate(BoundingBox()[1] + (margin, 0))
        stopat = 10
        r = range(0, length=400, stop=stopat)
        for p in 0:0.25:2
            translate(0, 70)
            sethue("red")
            graph(noise.(r, detail=4, persistence=0),
                endnumber=stopat, style=:circle)
            sethue("black")
            text("detail = 4, persistence = $p", Point(200, 15))
            graph(noise.(r, detail=4, persistence=p),
                endnumber=stopat)
        end
    end 800 675 "images/noise/persistence-graph.svg"
end

persistencegraph()

#-

# ![image label](IMAGEFOLDER/persistence-graph.svg){: .center-image}

# Here, the detail is kept at 4, and the persistence varies from 0 upwards. As the persistence increases, the effects accumulate, until the original curve is barely visible.

# There are many uses for noisy input, such as generating varying shapes that don't have that undesirable 'too random' quality.

using ColorSchemes

function treerings()
    @svg begin
        nrate = 0.01
        npoints= 500
        nrings = 400
        rad = 20
        setline(0.5)
        for ring in nrings:-5:1
            pts = Point[]
            for i in 1:npoints
                push!(pts, polar(rad + (ring * noise(i * nrate)),
                    rescale(i, 1, npoints, 0, 2pi)))
            end
            sethue(get(ColorSchemes.sienna, noise(ring * 5nrate)))
            poly(pts, :fill, close=false)
            sethue("black")
            poly(pts, :stroke, close=false)
        end
    end 800 800 "images/noise/treerings.svg"
end

treerings()

#-

# ![image label](IMAGEFOLDER/treerings.svg){: .center-image}

# Here's a more questionable idea, using noise to control the setting of a line of text.

function drawtextline(t, point, fsize; rate=0.1)
    for (n, c) in enumerate(split(t, ""))
        f = fsize * noise(n * rate, persistence=0, detail=4)
        fontsize(f)
        te = textextents(c)
        text(c, point)
        point = Point(point.x + te[5] * 0.98, 0) ## tightness is tight
        move(point)
    end
end

@png begin
    fontface("Bodoni")
    drawtextline("variablefontsizetextsettingiscool,orisit?",
        O - (380, 0), 50, rate=.11)
end 800 120 "images/noise/text-setting.png"

#-

# ![image label](IMAGEFOLDER/text-setting.png){: .center-image}

# I then used the `readpng()` and `placeimage()` functions to add a background image (the original Tenniel illustration), with the following result:

# ![image label](IMAGEFOLDER/jabberwocky.png){: .center-image}

# ### Why noise?

# The first use of computer graphics in movies is generally considered to be Tron (1981).

# Tron lies at the very beginning of the history of CGI in the movies, and the technology available to the artists, mathematicians, and programmers making Tron was amazingly underpowered compared with the computing power that we have today on our wrists, let alone on our phones.

# Ken Perlin was a mathematician and programmer who worked on Tron, and (I think after the film was released) he realised that there was room for using mathematical techniques for realistic-looking surfaces and textures, such as terrain.

# ![image label](IMAGEFOLDER/tron-800.png){: .center-image}

# Ken's noise, or Perlin Noise as it became known, was quickly adodpted as the best way to generate naturalistic surfaces.

# I think the reason why natural scenes appear to us as variable but not completely random is due to the (possibly hidden) larger scale processes that make smaller and more visible details clump together, and appear to work together and change gradually. For example, clouds, mountains, and pebble beaches have large scale structure controlled by unseen forces like heat, pressure, and gravity. We mostly see the objects that are subject to these forces, rather than the forces themselves.

# ### Moving into 2D

# So far the noise we've been producing has been one-dimensional, although we've been using 2D graphics to draw it.

# The `noise()` function can accept two floating-point numbers as input. These effectively define a rectangular grid of varying noise values: the x and y inputs produce a third value which requires representation.

# A simple way of doing this is to draw a table and vary the color of each square, giving a type of heat map.

using ColorSchemes

@svg begin
    nrows = 40
    ncols = 40
    cellwidth = 15
    cellheight = 15
    table = Table(nrows, ncols, cellwidth, cellheight)
    rate = 0.1
    fontsize(5)
    for row in 1:nrows
        for col in 1:ncols
            zvalue = noise(row * rate, col * rate)
            sethue(get(ColorSchemes.temperaturemap, zvalue))
            box(table[row, col], table.colwidths[1], table.rowheights[1], :fill)
            sethue("black")
            text(string(round(zvalue, digits=1)), table[row, col], halign=:center, valign=:middle)
        end
    end
end 800 700 "images/noise/table.svg"

#-

# ![image label](IMAGEFOLDER/table.svg){: .center-image}

# Alternatively, we can create a 3D surface and use the noise values for the height at each point. Normally this would require a visit to Julia's colorful and generally awesome Swedish-nightclub-themed Package manager, Pkg:

#-

# ![image label](IMAGEFOLDER/imstillfree.gif){: .center-image}

# to download some of the cool plotting packages available, not least Simon Danisch's impressive [Makie.jl](https://github.com/JuliaPlots/Makie.jl){:target="_blank"}.

# But, just to be contrary, I decided to whip up a simple isometric projection:

function project(x, y, z;
        scalingfactor = 3, heightmultiplier = -1)
    ## negative because y is positive downwards!
    u = (x - y)/sqrt(2)
    v = (x + 2(heightmultiplier * z) + y)/sqrt(6)
    return Point(scalingfactor * u, scalingfactor * v)
end

project(t; kwargs...) = project(t[1], t[2], t[3]; kwargs...)

function generatenoisearray(sx=100, sy=100;
    rate=0.5,
    detail=1,
    persistence=0)
    a = Array{Float64}(undef, sx, sy)
    for x in 1:sx
        for y in 1:sy
            a[x, y] = noise(x * rate, y * rate,
                detail=detail, persistence=persistence)
        end
    end
    return a
end

function isograph(a)
    @svg begin
        background("grey30")
        translate(0, -300)
        setline(0.5)
        sx, sy = size(a)
        scalingfactor = 5
        heightmultiplier = -6
        for x in 1:sx-1
            newpath()
            move(project(x, sy, -10,
                    scalingfactor = scalingfactor,
                    heightmultiplier = heightmultiplier))
            for y in sy-1:-1:1
                toppolygon = project.([
                    (x,     y,     a[x, y]),
                    (x + 1, y,     a[x + 1, y]),
                    (x + 1, y + 1, a[x + 1, y + 1]),
                    (x,     y + 1, a[x, y + 1])],
                        scalingfactor = scalingfactor,
                        heightmultiplier = heightmultiplier)
                centroid = polycentroid(toppolygon)
                line(centroid)
            end
            line(project(x, 1, -10,
                scalingfactor = scalingfactor,
                heightmultiplier = heightmultiplier))
            sethue("grey20")
            fillpreserve()
            sethue("grey85")
            strokepath()
        end
    end 800 700 "images/noise/isograph.svg"
end

isograph(generatenoisearray(80, 80, rate=0.08))

#-

# ![image label](IMAGEFOLDER/isograph.svg){: .center-image}

# (This image reminds me of the famous Joy Division LP cover and T-shirt image, which features plots of the first ever pulsar discovered by Jocelyn Bell Burnell and Antony Hewish in 1967. This then became the basis of many entertaining blog posts, such as [this one](https://adamcap.com/2011/05/19/history-of-joy-division-unknown-pleasures-album-art/){:target="_blank"}.)

# ![image label](IMAGEFOLDER/joydivision.jpg){: .center-image}

# A more conventional surface rendering is also possible:

function isograph(a)
    @png begin
        background("grey20")
        translate(0, -200)
        setline(0.5)
        sx, sy = size(a)
        for x in 1:sx-1
            for y in sy-1:-1:1
                toppolygon = project.([
                    (x,     y,     a[x, y]),
                    (x + 1, y,     a[x + 1, y]),
                    (x + 1, y + 1, a[x + 1, y + 1]),
                    (x,     y + 1, a[x, y + 1])],
                        heightmultiplier=-10,
                        scalingfactor=5)
                sethue("black")
                poly(toppolygon, close=true, :stroke)
                sethue(get(ColorSchemes.inferno, a[x, y]))
                poly(toppolygon, close=true, :fill)
            end
        end
    end 800 500 "images/noise/isosurface-2.png"
end

isograph(generatenoisearray(100, 100, rate=0.08))

#-

# ![image label](IMAGEFOLDER/isosurface-2.png){: .center-image}

# ### What, more dimensions?

# So far we've been generating 2D noise. The `noise()` function can also accept three floating-point numbers as input. This produces noise values in 3D space, where each 3D point can have a noise value between 0 and 1. Rendering these point clouds is definitely a job for something other than a simple 2D graphics system. But, while we're here, let's have a go:

function buildarray(a::AbstractArray; rate=20)
    sx, sy, sz = size(a)
    for x in 1:sx
        for y in 1:sy
            for z in 1:sz
                a[x, y, z] = noise(x * rate, y * rate, z * rate)
            end
        end
    end
    return a
end

function iso3d(a)
    background("grey20")
    sethue("gray80")
    setline(0.15)
    rule.([Point(0, y) for y in -400:10:400])
    sx, sy, sz = size(a)
    for x in 1:sx
        for y in 1:sy
            for z in 1:sz
                noisevalue = a[x, y, z]
                sethue(get(ColorSchemes.plasma, noisevalue))
                pt = project(x, y, z, scalingfactor=8)
                setopacity(noisevalue)
                circle(pt, rescale(noisevalue, 0, 1, 0.05, 6), :fill)
            end
        end
    end
end

const A = Array{Float64, 3}(undef, 50, 50, 50)

@png begin
    iso3d(buildarray(A, rate=0.05))
end 800 800 "images/noise/isosolid.png"

#-

# ![image label](IMAGEFOLDER/isosolid.png){: .center-image}

# The noise values nearer 1 look like hot plasma, whereas values nearer 0 are almost translucent. It suggests what you might expect to see from a real volume visualization tool.

# ## Journey to Algorithmia

# The final images in this post combine 2D noise and 1D noise; 2D noise for the sky, and 1D noise to create the contours.

# There's a `initnoise()` function. This is broadly the equivalent of the `Random.seed!()` function in Julia. This is useful when you want the noise to vary from image to image.

function layer(leftminheight, rightminheight, noiserate;
        detail=1, persistence=0)
    c1, c2, c3, c4 = box(BoundingBox(), vertices=true)
    ip1 = between(c1, c2, leftminheight)
    ip2 = between(c4, c3, rightminheight)
    topedge = Point[]
    initnoise(rand(1:12))
    for x in ip1.x:2:ip2.x
        ypos = between(ip1, ip2, rescale(x, ip1.x, ip2.x, 0, 1)).y
        ypos *= noise(x/noiserate, detail=detail, persistence=persistence)
        push!(topedge, Point(x, ypos))
    end
    p = [c1, topedge..., c4]
    poly(p, :fill, close=true)
end

function clouds()
    tiles = Tiler(boxwidth(BoundingBox()),
                  boxheight(BoundingBox()),
                  800, 800, margin=0)
    @layer begin
        transform([3 0 0 1 0 0])
        setopacity(0.3)
        noiserate = 0.03
        for (pos, n) in tiles
            nv = noise(pos.x * noiserate,
                       pos.y * noiserate,
                    detail=4, persistence=.4)
            setgray(nv)
            box(pos, tiles.tilewidth, tiles.tileheight, :fill)
        end
    end
end

function colorblend(fromcolor, tocolor, n=0.5)
    f = clamp(n, 0, 1)
    nc1 = convert(RGBA, fromcolor)
    nc2 = convert(RGBA, tocolor)
    from📕, from📗, from📘, from💡 =
        convert.(Float64, (nc1.r, nc1.g, nc1.b, nc1.alpha))
    to📕, to📗, to📘, to💡 =
        convert.(Float64, (nc2.r, nc2.g, nc2.b, nc1.alpha))
    new📕 = (f * (to📕 - from📕)) + from📕
    new📗 = (f * (to📗 - from📗)) + from📗
    new📘 = (f * (to📘 - from📘)) + from📘
    new💡 = (f * (to💡 - from💡)) + from💡
    return RGBA(new📕, new📗, new📘, new💡)
end

function landscape(scheme, filename)
    Drawing(800, 300, "$(filename).png")
    origin()
    ## sky is gradient mesh
    bb = BoundingBox()
    mesh1 = mesh(box(bb, vertices=true), [
      get(scheme, rand()),
      get(scheme, rand()),
      get(scheme, rand()),
      get(scheme, rand())
      ])
    setmesh(mesh1)
    box(bb, :fill)
    ## clouds are 2D noise
    clouds()
    ## the sun is a disk placed at random
    @layer begin
        setopacity(0.25)
        sethue(get(scheme, .95))
        sunposition = boxtop(bb) +
            (rand(-boxwidth(bb)/3:boxwidth(bb)/3), boxheight(bb)/10)
        circle(sunposition, boxdiagonal(bb)/30, :fill)
    end
    setopacity(0.8)
    ## how many layers
    len = 6
    noiselevels =  range(100, length=len, stop=10)
    detaillevels = 1:len
    persistencelevels = range(0.5, length=len, stop=0.95 )
    for (n, i) in enumerate(range(1, length=len, stop=0))
        ## avoid extremes of range
        sethue(colorblend(get(scheme, .05), get(scheme, .95), i))
        layer(i - rand()/2, i - rand()/2,
            noiselevels[n], detail=detaillevels[n],
            persistence=persistencelevels[n])
    end
    finish()
    preview()
end

landscape(ColorSchemes.leonardo, "images/noise/landscape-leonardo")

#-

# ![image label](IMAGEFOLDER/landscape-leonardo.png){: .center-image}

landscape(ColorSchemes.starrynight, "images/noise/landscapes-starrynight")

#-

# ![image label](IMAGEFOLDER/landscapes-starrynight.png){: .center-image}

# I generated a few hundred of these (there are over 300 colorschemes that can be selected at random) and, scrolling through them quickly, I found that sometimes the results were good, sometimes they weren't. Randomness—and noise—can be hard to predict.

# [2018-10-16]

# ![cormullion signing off](http://steampiano.net/cormullionknot.gif?noise){: .center-image}

#=                                                                             #src
using Literate                                                                 #src
# preprocess for notebooks                                                     #src
function setimagefolder(content)                                               #src
    content = replace(content, "IMAGEFOLDER" => "$IMAGEFOLDER")                #src
    return content                                                             #src
end                                                                            #src

# for Jupyter notebook, put images in subfolder                                #src
IMAGEFOLDER = "images/noise"                                                   #src

Literate.notebook("source/noise.jl", "notebooks", preprocess = setimagefolder) #src

# for Markdown/Jekyll notebook, put images in "/images"                        #src

IMAGEFOLDER = "/images/noise"                                                  #src

Literate.markdown("source/noise.jl", ".", name="_posts/2018-10-11-noise",      #src
 preprocess = setimagefolder,                                                  #src
 codefence = "{% highlight julia %}" => "{% endhighlight julia %}",            #src
 documenter=false)                                                             #src

=#                                                                             #src


text:
# add `Ref`s
lazyable(x, ::Val) = Ref(x)
lazyable(x::Base.RefValue, ::Val) = x
lazyable(x::AbstractLayeredArray{layer}, ::Val{layer}) where {layer} = x
@generated function lazyable(x::Broadcasted{<: Any, <: Any, <: Any, Tup}, ::Val) where {Tup <: Tuple}
    any(T -> T <: AbstractLayeredArray, Tup.parameters) && return :(@unreachable)
    quote
        Ref(copy(x))
    end
end
@generated function lazyables(f, args...)
    layer = maximum(layerof, args)
    exps = [:(lazyable(args[$i], Val($layer))) for i in 1:length(args)]
    quote
        @_inline_meta
        $layer, tuple($(exps...))
    end
end

function return_eltype(f, args...)
    T = Base._return_type(_propagate_lazy, eltypes((f,args...)))
    if T == Union{}
        f(map(first, args)...)
        @unreachable
    end
    T
end
_eltype(x::AbstractLayeredArray) = eltype(x)
_eltype(x::Base.RefValue) = eltype(x)
_eltype(x) = typeof(x)
eltypes(::Tuple{}) = Tuple{}
eltypes(t::Tuple{Any}) = Tuple{_eltype(t[1])}
eltypes(t::Tuple{Any, Any}) = Tuple{_eltype(t[1]), _eltype(t[2])}
eltypes(t::Tuple) = Tuple{_eltype(t[1]), eltypes(Base.tail(t)).types...}

struct LazyLayeredArray{layer, T, N, BC <: Broadcasted{<: AbstractArrayStyle{N}}} <: AbstractLayeredArray{layer, T, N}
    bc::BC
    function LazyLayeredArray{layer, T, N, BC}(bc::BC) where {layer, T, N, BC}
        new{layer::Int, T, N, BC}(bc)
    end
end

const LazyLayeredVector{layer, T} = LazyLayeredArray{layer, T, 1, <: Broadcasted{LayeredArrayStyle{1}}}
const LazyLayeredMatrix{layer, T} = LazyLayeredArray{layer, T, 2, <: Broadcasted{LayeredArrayStyle{2}}}

@inline function LazyLayeredArray{layer, T}(bc::BC) where {layer, T, N, BC <: Broadcasted{<: AbstractArrayStyle{N}}}
    LazyLayeredArray{layer, T, N, BC}(bc)
end

function LazyLayeredArray(f, args...)
    layer, args′ = lazyables(f, args...)
    T = return_eltype(f, args′...)
    LazyLayeredArray{layer, T}(Broadcast.broadcasted(f, args′...))
end

Base.size(x::LazyLayeredArray) = size(x.bc)
Base.axes(x::LazyLayeredArray) = axes(x.bc)

# this propagates lazy operation when any AbstractLayeredArray is found
# otherwise just normally call function `f`.
@generated function _propagate_lazy(f, args...)
    any([t <: AbstractLayeredArray for t in args]) ?
    :(LazyLayeredArray(f, args...)) : :(f(args...))
end
_propagate_lazy(f, arg) = f(arg) # this prevents too much propagation

@inline @propagate_inbounds _getindex(x::AbstractLayeredArray, i) = x[Broadcast.newindex(x,i)]
@inline @propagate_inbounds _getindex(x::AbstractLayeredVector, i::CartesianIndex{1}) = x[i[1]]
@inline @propagate_inbounds _getindex(x::Base.RefValue, i) = x[]
@inline @propagate_inbounds _getindex_broadcast(x::Tuple{Any}, i) = (_getindex(x[1], i),)
@inline @propagate_inbounds _getindex_broadcast(x::Tuple{Any, Any}, i) = (_getindex(x[1], i), _getindex(x[2], i))
@inline @propagate_inbounds _getindex_broadcast(x::Tuple, i) = (_getindex(x[1], i), _getindex_broadcast(Base.tail(x), i)...)
@inline function Base.getindex(x::LazyLayeredArray, I::Int...)
    @boundscheck checkbounds(x, I...)
    bc = x.bc
    @inbounds _propagate_lazy(bc.f, _getindex_broadcast(bc.args, CartesianIndex(I))...)
end

for op in (:+, :-)
    @eval function Base.$op(x::AbstractLayeredArray, y::AbstractLayeredArray)
        layerof(x) == layerof(y) || throw(ArgumentError("layers must match in $($op) operation"))
        LayeredArray{layerof(x)}(map($op, x, y))
    end
end


text: pkg"add Documenter"


text:
##  Copyright 2021 Myrl Marmarelis
##
##  Licensed under the Apache License, Version 2.0 (the "License");
##  you may not use this file except in compliance with the License.
##  You may obtain a copy of the License at
##
##    http://www.apache.org/licenses/LICENSE-2.0
##
##  Unless required by applicable law or agreed to in writing, software
##  distributed under the License is distributed on an "AS IS" BASIS,
##  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
##  See the License for the specific language governing permissions and
##  limitations under the License.

# The prior tends to constrain the range and scale of the inferred distributions,
# but they can still float around during the Gibbs sampling.
# I think that tighter decisions can be made by looking for metrics that are
# translation- and scale- invariant across all the means of the mixture, at least.
# Metrics that are functions of the covariances can be examined later. It is also
# not entirely necessary to do this.

function sample_possible_rankings(mixtures::AbstractVector{V},
    ::Val{R}, sample_size::Int) where {R,N,M,T,V<:VoterMixture{N,M,T}}
  utility = zeros(T, N)
  ranking = zeros(Int, N)
  memberships = zeros(Int, sample_size)
  choices = map(1:sample_size) do i
    mixture = sample(mixtures)
    membership = Categorical(mixture.shares) |> rand
    memberships[i] = membership
    cohort = mixture.cohorts[membership]
    rand!(cohort, utility)
    sortperm!(ranking, utility, order=Base.Reverse)
    convert(RankedChoice{R}, @view ranking[1:R])
  end
  IssueVote{R}(choices, N), memberships
end

function bootstrap_statistic(statistic::Function, mixtures::AbstractVector{V};
    sample_size::Int = length(mixtures)) where {V<:VoterMixture}
  mean(1:sample_size) do _
    sample(mixtures) |> statistic
  end
end

using Statistics
using Distributions
using StatsFuns

Statistics.mean(mixture::VoterMixture{N,M,T}) where {N,M,T} =
  sum(zip(mixture.cohorts, mixture.shares)) do params
    cohort, share = params
    means = convert(SVector{N,T}, mean(cohort))
    means .* share
  end

# https://math.stackexchange.com/questions/195911/calculation-of-the-covariance-of-gaussian-mixtures
Statistics.cov(mixture::VoterMixture{N,M,T}) where {N,M,T} =
  let mixture_means = mean(mixture)
    sum(zip(mixture.cohorts, mixture.shares)) do params
      cohort, share = params
      # these are effectively no-ops with the current data structures in place
      means = convert(SVector{N,T}, mean(cohort))
      covariance = convert(SMatrix{N,N,T}, cov(cohort))
      displacement = means - mixture_means # how much off-center
      cohort_cov = covariance + displacement * displacement'
      cohort_cov .* share
    end
  end

# normalize PER MCMC draw to remove drift!
# fix the normalization function further below for a more general scale/shift correction?
Statistics.cor(mixture::VoterMixture{N,M,T}) where {N,M,T} =
  let covariance = cov(mixture)
    covariance ./ sqrt.(diag(covariance) * diag(covariance)')
  end

# probability of pairwise event (utlity[i] > utility[j]) --- equality has infinitesimal
function contest_probability(mixture::VoterMixture, i::Int, j::Int; log::Bool=false)
  means = mean(mixture)
  covariances = cov(mixture)
  diff_mean = means[i] - means[j]
  first_var = covariances[i, i]
  second_var = covariances[j, j] # slightly inclined to spell as "secnd" for alignment
  covariance = covariances[i, j]
  deviation = sqrt( first_var + second_var - 2covariance )
  (log ? normlogcdf : normcdf)( diff_mean / deviation ) # locally consistent styling..
end

# I know that this works with two blocks of variables, at least
# we can simply average this over monte carlo posteriors but NOT mixture components!
function mutual_information(cohort::VoterCohort{T},
    blocks::AbstractVector{Int}...)::T where T
  @assert length(blocks) > 1
  covariance = cov(cohort)
  numerator = sum(blocks) do block
    @view(covariance[block, block]) |> logdet
  end
  coverage = vcat(blocks...)
  @views denominator = covariance[coverage, coverage] |> logdet #logdetcov(cohort)
  0.5(numerator - denominator)
end

# set a low quantile to find the "safest" candidate, i.e. most moderately agreeable, least polarizing
function estimate_quantiles(mixtures::AbstractVector{V}, portion::Union{Vector{T},T},
    n_draws::Y)::AbstractVector{Utility{N,T}} where {N,T,Y<:Real,V<:VoterMixture}
  n_quantiles = length(portion)
  quantiles = zeros(T, N, n_quantiles) # align for reinterpretation, even though quantile! will have to skip around
  for candidate in 1:N
    sample = T[]
    remainder = zero(Y)
    for mixture in mixtures
      remainder += n_draws
      n_round_draws = floor(Int, remainder)
      remainder -= convert(Y, n_round_draws)
      for draw in 1:n_round_draws
        utility = draw_marginal(mixture, candidate)
        push!(sample, utility) # does it even help to pre-allocate?
      end
    end
    candidate_quantiles = @view quantiles[candidate, :]
    quantile!(candidate_quantiles, sample, portion)
  end
  reinterpret(Utility{N,T}, vec(quantiles))
end

# these functions are designed to be mapped across a mixture sample as a pipeline

function normalize_utility(utilities::Utility{N,T}, norm_power::T=T(2))::Utility{N,T} where {N,T}
  centered = utilities .- mean(utilities)
  norm = mean(u -> abs(u) ^ norm_power, utilities) ^ (1/norm_power)
  centered ./ norm
end

function estimate_quantiles_old(mixture::VoterMixture{N,M,T}, quantile::T)::Utility{N,T} where {N,M,T}
  # obtain the marginals. covariance matrices are magical in that the diagonal tells you all you need
  # to know about individual variables (following the positive-definite construction)
  overall_quantiles = @SVector zeros(T, N)
  for (share, cohort) in zip(mixture.shares, mixture.cohorts)
    quantiles = marginalize_quantiles(Val(N), cohort, quantile)
    overall_quantiles += quantiles * share
  end
  overall_quantiles
end

function marginalize_quantiles(::Val{N}, cohort::VoterCohort{T}, quantile::T)::Utility{N,T} where {N,T}
  means = mean(cohort) # hope these don't allocate new arrays...
  covariance = cov(cohort)
  SVector{N}( begin # @SVector [...] is more precarious, syntactically
    stdev = covariance[i, i] |> sqrt
    normal_quantile(means[i], stdev, quantile)
  end for i in 1:N )
end

using SpecialFunctions: erfinv

function normal_quantile(mu::T, sigma::T, quantile::T)::T where T <: Real
  mu + sigma*sqrt(2)*erfinv(2quantile - 1)
end


text:
@testset "Collectors" begin

    # every 500 ms from 1 sec to 5 sec
    dts = map(x -> DateTime(2000, 1, 1) + Millisecond(x), 1000:500:5000)
    data = [100.0, 110.0, 99.0, 102.0, 105.0, 105.0, 105.0, 120.0, 110.0]

    begin
        # periodic_collector
        f, collected = periodic_collector(Float64, Second(1))

        for i in 1:length(dts)
            f(dts[i], data[i])
        end

        @test length(collected.values) == 5
        @test all(map(x -> x[1], collected.values) .== map(x -> DateTime(2000, 1, 1) + Second(x), 1:5))
        @test collected.last_dt == dts[end]
    end

    begin
        # predicate_collector
        predicate = (collected, dt, value) -> (dt - collected.last_dt) >= Second(1)
        f, collected = predicate_collector(Float64, predicate, 0.0)

        for i in 1:length(dts)
            f(dts[i], data[i])
        end

        @test length(collected.values) == 5
        @test all(map(x -> x[1], collected.values) .== map(x -> DateTime(2000, 1, 1) + Second(x), 1:5))
        @test collected.last_dt == dts[end]
    end

    begin
        # min_value_collector
        f, collected = min_value_collector(Float64)

        for i in 1:length(dts)
            f(dts[i], data[i])
        end

        @test collected.min_value == minimum(data)
        @test collected.dt == dts[indexin(minimum(data), data)][1]
    end

    begin
        # max_value_collector
        f, collected = max_value_collector(Float64)

        for i in 1:length(dts)
            f(dts[i], data[i])
        end

        @test collected.max_value == maximum(data)
        @test collected.dt == dts[indexin(maximum(data), data)][1]
    end

    begin
        # drawdown_collector (P&L)
        p = (v, dt, equity) -> dt - v.last_dt >= Second(1)
        f, collected = drawdown_collector(PnL::DrawdownMode, p)

        for i in 1:length(dts)
            f(dts[i], data[i])
        end

        @test length(collected.values) == 5
        @test all(map(x -> x[1], collected.values) .== map(x -> DateTime(2000, 1, 1) + Second(x), 1:5))
        @test collected.last_dt == dts[end]
        @test all(map(x -> x[2], collected.values) .== [0.0, -11, -5, -5, -10])
    end

    begin
        # drawdown_collector (%)
        p = (dv, dt, equity) -> dt - dv.last_dt >= Second(1)
        f, collected = drawdown_collector(Percentage::DrawdownMode, p)

        for i in 1:length(dts)
            f(dts[i], data[i])
        end

        @test length(collected.values) == 5
        @test all(map(x -> x[1], collected.values) .== map(x -> DateTime(2000, 1, 1) + Second(x), 1:5))
        @test collected.last_dt == dts[end]
        @test all(map(x -> x[2], collected.values) .≈ [0.0, -11/110, -5/110, -5/110, -10/120])
    end

end


text:
function URLify(str :: ASCIIString, real_length :: Int)
  char_array :: Array{Char} = map(single_string -> single_string[1], split(str, ""))

  idx :: Int = length(char_array)

  for i = real_length:-1:1
    if char_array[i] != ' '             # non-space
      char_array[idx] = char_array[i]   # just copy
      idx -= 1
    else                                # space
      char_array[idx] = '0'             # fill "%20"
      char_array[idx - 1] = '2'
      char_array[idx - 2] = '%'
      idx -= 3
    end
  end

  return ASCIIString(char_array)
end

# Time  : O(n)
# Space : O(1)
# Note  : no extra space if assume the input is Array{Char}

using Base: Test

@test URLify("Mr Jhon Smith    ", 13) == "Mr%20Jhon%20Smith"

max_stars_repo_name: Kiandr/crackingcodinginterview

max_issues_repo_name: Kiandr/crackingcodinginterview

max_forks_repo_name: Kiandr/crackingcodinginterview


text:
#=
Wavefunction related routines

The goal is to recover the full wavefunction grid from the planewave coefficients
=#
import FFTW
using FFTW
using LinearAlgebra

export WaveFunction, ChargeDensity, chargedensity, sliceband, slicekpoint, slicespin, slicespinor

"""
    coeff_to_recip(coeff_array, nwaves_at_kp, grid_coords, ngx, ngy, ngz)

Convert planewave coefficients in the reduced representation to grids
"""
function coeff_to_recip(coeff_array, nwaves_at_kp, grid_coords, ngx, ngy, ngz)
    npw, nspinor, band_max, nkpts, nspins = size(coeff_array)
    grid_size = Int[ngx, ngy, ngz]
    indices = coord_to_indices(grid_coords, grid_size)
    # Grid data index with ngx, ngy, ngz, nspinor, band_max, nkpts, nspins
    # Note that the kpoints order may be different from those read from other sections
    # so be careful with that.....
    grid = zeros(ComplexF64, ngx, ngy, ngz, nspinor, band_max, nkpts, nspins)

    for is in 1:nspins
        for ik in 1:nkpts
            for ib in 1:band_max, ispinor in 1:nspinor
                for ipw in 1:nwaves_at_kp[ik]
                    grid[indices[1, ipw, ik], indices[2, ipw, ik],indices[3, ipw, ik],
                         ispinor, ib, ik, is] = coeff_array[ipw, ispinor, ib, ik, is]
                end
            end
        end
    end
    grid
end

"""
    coeff_to_recip(coeff_array, data_dict::Dict)

Convert plane wave coefficients to a reciprocal space grid.
"""
function coeff_to_recip(coeff_array, data_dict::Dict)
    wave = data_dict[:WAVEFUNCTION_DATA]
    coeff_to_recip(coeff_array, wave.nwaves_at_kp, wave.pw_grid_coord, wave.ngx, wave.ngy, wave.ngz)
end

"""
Convert full grid representation back to reciprocal coefficients
The original specifications of the wave function needs to be reused
"""
function recip_to_coeff(grid, coeff_array_orig, nwaves_at_kp, grid_coords, ngx, ngy, ngz)
    npw, nspinor, band_max, nkpts, nspins = size(coeff_array_orig)
    grid_size = Int[ngx, ngy, ngz]
    indices = coord_to_indices(grid_coords, grid_size)
    # Grid data index with ngx, ngy, ngz, nspinor, band_max, nkpts, nspins
    # Note that the kpoints order may be different from those read from other sections
    # so be careful with that.....
    # Allocate the coefficients
    coeff_array = similar(coeff_array_orig)
    fill!(coeff_array, 0.)

    for is in 1:nspins
        for ik in 1:nkpts
            for ib in 1:band_max, ispinor in 1:nspinor
                for ipw in 1:nwaves_at_kp[ik]
                    coeff_array[ipw, ispinor, ib, ik, is] = grid[indices[1, ipw, ik], indices[2, ipw, ik],indices[3, ipw, ik],
                         ispinor, ib, ik, is]
                end
            end
        end
    end
    coeff_array
end

function recip_to_coeff(grid, data_dict)
    wave = data_dict[:WAVEFUNCTION_DATA]
    recip_to_coeff(grid, wave.coeffs, wave.nwaves_at_kp, wave.pw_grid_coord, wave.ngx, wave.ngy, wave.ngz)
end

"""
Convert G-vector coordinates (unit of reciprocal lattice basis) to 1-based grid indices
"""
function coord_to_indices(grid_coords, grid_size)
    indices = copy(grid_coords)
    _, npw, nkpts = size(grid_coords)
    # The grid_coords have sizes in -n/2 to n/2 - we move it to 1 to n
    for k in 1:nkpts
        for j in 1:npw
            for i in 1:3
                if indices[i, j, k] < 0
                    indices[i, j, k] += grid_size[i]
                end
            end
        end
    end
    # Switch to 1 based indexing
    indices .+= 1
end

"""
Load wavefunction using data read from the checkpoint file.

Returns the reciprocal representation on the grids.
Needs to be FFTed to get to the real space representation.
"""
function wavefunction(data_dict::Dict)
    wave = data_dict[:WAVEFUNCTION_DATA]
    coeff_to_recip(wave.coeffs, wave.nwaves_at_kp, wave.pw_grid_coord, wave.ngx, wave.ngy, wave.ngz)
end

"""
Type for storing a wavefunction
"""
mutable struct WaveFunction
    wave::Array{ComplexF64, 7}
    real_space::Bool
    ngx::Int
    ngy::Int
    ngz::Int
    nspinor::Int
    nbands::Int
    nkpts::Int
    nspins::Int
    basis::Matrix{Float64}
    "Occupations in (ib, ik, is)"
    occupations::Array{Float64, 3}
    "Eigenvalues in (ib, ik, is)"
    eigenvalues::Array{Float64, 3}
    "Weight of the kpoints"
    kpoint_weights::Array{Float64, 1}
    "Coordinates of the kpoints"
    kpoints::Array{Float64, 2}
end

struct ChargeDensity{T}
    density::Array{T, 4}
    ngx::Int
    ngy::Int
    ngz::Int
    nspins::Int
    basis::Matrix{Float64}
end

"""
    ChargeDensity(density)

Construction from a density array.
"""
function ChargeDensity(density, basis=diagm([10., 10., 10.]))
    ChargeDensity(density, size(density)..., basis)
end

"""
    ChargeDensity(ngx, ngy, ngz)

Construction an empty ChargeDensity.
"""
function ChargeDensity{T}(ngx::Int, ngy::Int, ngz::Int) where {T}
    ChargeDensity(zeros(T, ngx, ngy, ngz))
end

"""
    WaveFunction(data_dict::Dict)

Construct the WaveFunction from a Dict of read data from castep_bin/check file
"""
function WaveFunction(data_dict::Dict)
    wave = wavefunction(data_dict)
    ngx, ngy, ngz, nspinor, nbands, nkpts, nspins = size(wave)
    basis = data_dict[:REAL_LATTICE] .* Bohr2Ang
    eigenvalues = data_dict[:EIGENVALUES] .* Ha2eV
    occupations = data_dict[:OCCUPATIONS]

    kpoints_weights = data_dict[:KPOINT_WEIGHTS]
    kpoints = data_dict[:KPOINTS]
    WaveFunction(wave, false, ngx, ngy, ngz, nspinor, nbands, nkpts, nspins,
                 basis, occupations, eigenvalues, kpoints_weights, kpoints)
end

slicespinor(wave::WaveFunction, ispinor) = selectdim(wave.wave, 4, ispinor)

sliceband(wave::WaveFunction, ib) = selectdim(wave.wave, 5, ib)

slicekpoint(wave::WaveFunction, ik) = selectdim(wave.wave, 6, ik)

slicespin(wave::WaveFunction, is) = selectdim(wave.wave, 7, is)

"""
    FFTW.fft!(wavef::WaveFunction)
In place FFT to transfer the wavefunction to the frequency space

The in-place transformation is apply slice by slice
"""
function FFTW.fft!(wavef::WaveFunction)
    @assert wavef.real_space "WaveFunction is already in the frequency space"
    xyz_slice = zeros(ComplexF64, wavef.ngx, wavef.ngy, wavef.ngz)
    for i=1:wavef.nspins,j=1:wavef.nkpts,k=1:wavef.nbands,l=1:wavef.nspinor
        xyz_slice .= wavef.wave[:, :, :, l, k, j, i]
        wavef.wave[:, :, :, l, k, j, i] = FFTW.fft(xyz_slice)
    end
    wavef.real_space = false
    wavef
end

"""
    fft(wavef::WaveFunction)
FFT to transfer the wavefunction to the frequency space
"""
function FFTW.fft(wavef::WaveFunction)
    @assert wavef.real_space "WaveFunction is already in the frequency space"
    new_wave = deepcopy(wavef)
    FFTW.fft!(new_wave)
end

"""
In place FFT to transfer the wavefunction to the real space

The in-place transformation is apply slice by slice
"""
function FFTW.ifft!(wavef::WaveFunction)
    @assert ~wavef.real_space "WaveFunction is already in the real space"
    xyz_slice = zeros(ComplexF64, wavef.ngx, wavef.ngy, wavef.ngz)
    for i=1:wavef.nspins,j=1:wavef.nkpts,k=1:wavef.nbands,l=1:wavef.nspinor
        xyz_slice .= wavef.wave[:, :, :, l, k, j, i]
        wavef.wave[:, :, :, l, k, j, i] = FFTW.ifft(xyz_slice)
    end
    wavef.real_space = true
    wavef
end

"""
    fft(wavef::WaveFunction)
FFT to transfer the wavefunction to the frequency space
"""
function FFTW.ifft(wavef::WaveFunction)
    @assert ~wavef.real_space "WaveFunction is already in the real space"
    new_wave = deepcopy(wavef)
    FFTW.ifft!(new_wave)
end

ensure_recip!(wavef::WaveFunction) = wavef.real_space && fft!(wavef)
ensure_real!(wavef::WaveFunction) = wavef.real_space || ifft!(wavef)

"""
    chargedensity(wavef::WaveFunction, kidx, bidx, sidx)

Construct charge density for given kpoints, band, and spin indices
"""
function chargedensity(wavef::WaveFunction, bidx, kidx, sidx; weighted=true, ivbm=-1)
    orig_state = wavef.real_space
    ensure_real!(wavef)
    density = zeros(wavef.ngx, wavef.ngy, wavef.ngz, wavef.nspins)
    # Force occupations
    if ivbm > 1
        occ = zeros(size(wavef.occupations))
        for ib =1:ivbm
            occ[ib, : , :] .= 1.0
        end
    else
        occ = wavef.occupations
    end
    for is in sidx,ik in kidx,ib in bidx,ispinor=1:wavef.nspinor
        for z=1:wavef.ngz, y=1:wavef.ngy, x=1:wavef.ngx
            val = norm(wavef.wave[x, y, z, ispinor, ib, ik, is])
            # Do we weight by kpoint weights and occupations?
            if weighted
                val = val * occ[ib, ik, is] * wavef.kpoint_weights[ik]
            end
            density[x, y, z, is] += val
        end
    end
    orig_state == false && ensure_recip!(wavef)
    ChargeDensity(density, wavef.basis)
end

"""
    chargedensity(wavef::WaveFunction)

Compute charge density from the wavefunction

Note that the density does not include the augmentation charges involved in ultrasoft pseudopotentials
"""
function chargedensity(wavef::WaveFunction; ivbm=-1)
    @assert ivbm > 1 || sum(wavef.occupations[:, 1, 1]) > 0. "VBM index must be passed for unweighted WaveFunction"
    chargedensity(wavef, 1:wavef.nbands, 1:wavef.nkpts, 1:wavef.nspins;ivbm=ivbm)
end

"""
    chargedensity(wavef::WaveFunction, ib, ik, is)

Get the charge density for a single slice of band, kpoint, spin channel
"""
function chargedensity(wavef::WaveFunction, ib::Int, ik::Int, is::Int)
    chargedensity(wavef, ib:ib, ik:ik, is:is;weighted=false)
end

"""
Interpolate charge density from `c1` into `c2` with a finner grid
"""
function fine_grid_interpolate(c1::ChargeDensity{T}, c2::ChargeDensity{T}) where {T}
    # Zero the grid to be interpolated into
    fill!(0., c2.density)
    ngx, ngy, ngz = c1.ngx, c1.ngy, c1.ngz
    ngxf, ngyf, ngzf = c2.ngx, c2.ngy, c2.ngz
    @assert (ngxf > ngx) && (ngyf > ngy) && (ngzf > ngz)

    workspace = zeros(ComplexF64, ngx, ngy, ngz)
    workspace_fine = zeros(ComplexF64, ngxf, ngyf, ngzf)
    for is in 1:c1.nspins
        workspace .= c1.density[:, :, :, is]
        # To reciprocal space
        fft!(workspace)
        for z = 1:ngz, y=1:ngy, x=1:ngx
            # Compute the index in the interpolated grid in the reciprocal space
            x < ngx/2+1 ? fx = x : fx = ngxf - ngx + x
            y < ngy/2+1 ? fy = y : fy = ngyf - ngy + y
            z < ngz/2+1 ? fz = z : fz = ngzf - ngz + z
            work_space_fine[fx, fy, fz] = workspace[x, y, z]
        end
        # Back to real space
        ifft!(workspace_fine)
        # Function barrier for copying back the density
        copy_density(c2, workspace_fine, is)
    end
end

"""
Copy density from a temporary array in to the spin channel
"""
function copy_density(density::ChargeDensity{ComplexF64}, tmp, spin)
    density.density[:, :, :, spin] .= tmp
end

"""
Copy density from a temporary array in to the spin channel
"""
function copy_density(density::ChargeDensity{Float64}, tmp, spin)
    density.density[:, :, :, spin] .= real.(tmp)
end

"""
Computed weighted density from some given weights by bands
"""
function weighted_density_by_bands(fname, weights)
    data = read_castep_check(fname)
    nbands = data[:NBANDS]
    nspins = data[:NSPINS]
    nkpts = data[:NKPTS]
    @assert size(weights, 1) nbands
    wavef = WaveFunction(data)
    ifft!(wavef)
    chargedensity(wavef, ensure_three_dims(weights, nkpts, nspins), data[:KPOINTS_WEIGHTS])
end

ensure_three_dims(weights::AbstractArray{T, 1}, d1, d2) where {T} = repeat(weights[:, :, :], d1, d2)
ensure_three_dims(weights::AbstractArray{T, 2}, d1, d2) where {T} = repeat(weights[:, :, :], d2)
ensure_three_dims(weights::AbstractArray{T, 3}, d1, d2) where {T} = weights

"""
Update wavefunction of an existing file
"""
function CastepBin.update_wavefunction(src, dst, wavef::WaveFunction, data::Dict=read_castep_check(src))
    @assert wavef.real_space == false "Wave must be in the reciprocal space"
    new_coeff = recip_to_coeff(wavef.wave, data)
    CastepBin.update_wavefunction(src, dst, new_coeff)
end

"""
Update wavefunction of an existing file
"""
function CastepBin.update_wavefunction!(src, wavef::WaveFunction, data::Dict=read_castep_check(src))
    @assert wavef.real_space == false "Wave must be in the reciprocal space"
    new_coeff = recip_to_coeff(wavef.wave, data)
    CastepBin.update_wavefunction!(src, new_coeff)
end


text:
push!(LOAD_PATH, pwd()*"/../../src/")
using Distributions
using PyPlot
using BVB

x = linspace(-15, 15, 500)
init_approx_d = Normal(0, 20)

lower = -10.
upper = 10.
max_iter = 200

# example 1
@show true_d = TriangularDist(-10, 10, 0)
ll(x) = logpdf(true_d, x)

# include("run_laplacian_gradient_boosting.jl")
include("run_laplacian_gradient_boosting_exclusive.jl")
show()

println("Done.")
readline(STDIN)

max_stars_repo_name: richardkwo/boosting-VB

max_issues_repo_name: richardkwo/boosting-VB

max_forks_repo_name: richardkwo/boosting-VB


text:
module GGPlot2

__precompile__(false)

using FileIO
using ImageShow
using Random
using Reexport
@reexport using RCall
__temp__ = rimport("ggplot2")
@reexport using .__temp__

include("io.jl")

end # module


text:
#= ===================
   LOGISTIC CLASSIFIER
   =================== =#

@with_kw_noshow mutable struct LogisticClassifier <: MMI.Probabilistic
    lambda::Real             = 1.0
    gamma::Real              = 0.0
    penalty::SymStr          = :l2
    fit_intercept::Bool      = true
    penalize_intercept::Bool = false
    solver::Option{Solver}   = nothing
    multi_class::Bool        = false
end

glr(m::LogisticClassifier) =
    LogisticRegression(m.lambda, m.gamma;
                       penalty=Symbol(m.penalty),
                       multi_class=m.multi_class,
                       fit_intercept=m.fit_intercept,
                       penalize_intercept=m.penalize_intercept)

descr(::Type{LogisticClassifier}) = "Classifier corresponding to the loss function ``L(y, Xθ) + λ|θ|₂²/2 + γ|θ|₁`` where `L` is the logistic loss."

#= ======================
   MULTINOMIAL CLASSIFIER
   ====================== =#

@with_kw_noshow mutable struct MultinomialClassifier <: MMI.Probabilistic
    lambda::Real             = 1.0
    gamma::Real              = 0.0
    penalty::SymStr          = :l2
    fit_intercept::Bool      = true
    penalize_intercept::Bool = false
    solver::Option{Solver}   = nothing
end

glr(m::MultinomialClassifier) =
    MultinomialRegression(m.lambda, m.gamma;
                          penalty=Symbol(m.penalty),
                          fit_intercept=m.fit_intercept,
                          penalize_intercept=m.penalize_intercept)

descr(::Type{MultinomialClassifier}) =
    "Classifier corresponding to the loss function " *
    "``L(y, Xθ) + λ|θ|₂²/2 + γ|θ|₁`` where `L` is the multinomial loss."